{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":116062,"databundleVersionId":14084779,"sourceType":"competition"}],"dockerImageVersionId":31153,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# Standard libs\nimport os, json, re\nfrom collections import defaultdict, Counter\n\n# Data & math\nimport numpy as np\nimport pandas as pd\n\n# Modeling\nfrom sklearn.preprocessing import MultiLabelBinarizer\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.neural_network import MLPClassifier\nfrom sklearn.multiclass import OneVsRestClassifier\nfrom sklearn.metrics import precision_score, recall_score, f1_score\n\n# Paths\nBASE = \"/kaggle/input/cafa-6-protein-function-prediction/Train\"\nPATH_OBO   = os.path.join(BASE, \"go-basic.obo\")\nPATH_FASTA = os.path.join(BASE, \"train_sequences.fasta\")\nPATH_TAX   = os.path.join(BASE, \"train_taxonomy.tsv\")\nPATH_TERMS = os.path.join(BASE, \"train_terms.tsv\")\n\nprint(\"Exists:\")\nfor p in [PATH_OBO, PATH_FASTA, PATH_TAX, PATH_TERMS]:\n    print(p, os.path.exists(p))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-26T15:41:21.370174Z","iopub.execute_input":"2025-10-26T15:41:21.370438Z","iopub.status.idle":"2025-10-26T15:41:24.774128Z","shell.execute_reply.started":"2025-10-26T15:41:21.370413Z","shell.execute_reply":"2025-10-26T15:41:24.773290Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def peek_file(path, n=5):\n    print(f\"== Peek {path} ==\")\n    with open(path, 'r', encoding='utf-8', errors='ignore') as f:\n        for i, line in zip(range(n), f):\n            print(line.rstrip())\n    print()\n\npeek_file(PATH_FASTA, 3)\npeek_file(PATH_TERMS, 10)\npeek_file(PATH_TAX, 5)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-26T15:41:24.775023Z","iopub.execute_input":"2025-10-26T15:41:24.775444Z","iopub.status.idle":"2025-10-26T15:41:24.790661Z","shell.execute_reply.started":"2025-10-26T15:41:24.775411Z","shell.execute_reply":"2025-10-26T15:41:24.789927Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## FASTA Reader + ID Normalization","metadata":{}},{"cell_type":"code","source":"# Simple FASTA reader\n\ndef read_fasta(path):\n    seqs = {}\n    with open(path, 'r') as fh:\n        sid, buf = None, []\n        for line in fh:\n            line = line.strip()\n            if not line:\n                continue\n            if line.startswith('>'):\n                if sid is not None:\n                    seqs[sid] = ''.join(buf).upper()\n                    buf = []\n                sid = line[1:].split()[0]\n            else:\n                buf.append(line)\n        if sid is not None:\n            seqs[sid] = ''.join(buf).upper()\n    return seqs\n\n# Robust ID normalization to align across files\n_UNIPROT_LIKE = re.compile(r\"^[A-Z0-9]{5,12}$\")\n\ndef normalize_seq_id(raw: str) -> str:\n    if raw is None:\n        return \"\"\n    s = str(raw).strip()\n    if not s:\n        return \"\"\n    s = s.split()[0]              # take token before whitespace\n    parts = s.split('|')          # common in FASTA (sp|Q9..., tr|A0A...)\n    def score(tok):\n        t = tok.split('.')[0]     # drop version suffixes\n        if _UNIPROT_LIKE.match(t):\n            return (2, len(t))\n        if t.isalnum():\n            return (1, len(t))\n        return (0, len(t))\n    best = max(parts, key=score)\n    best = best.split('.')[0]\n    return best.upper()\n\n# Normalize FASTA ids\nseqs_raw = read_fasta(PATH_FASTA)\nseqs_norm = {normalize_seq_id(k): v for k,v in seqs_raw.items()}\nprint(f\"Loaded {len(seqs_raw)} FASTA entries; {len(seqs_norm)} after normalization\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-26T15:41:24.792621Z","iopub.execute_input":"2025-10-26T15:41:24.792845Z","iopub.status.idle":"2025-10-26T15:41:25.780182Z","shell.execute_reply.started":"2025-10-26T15:41:24.792829Z","shell.execute_reply":"2025-10-26T15:41:25.779456Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Cell 4 — Load train_terms.tsv (Flexible) + Normalize IDs","metadata":{}},{"cell_type":"code","source":"# Handle both: (a) headered TSV or (b) 2 columns without header\n\ndef load_terms_flexible(path):\n    try:\n        df = pd.read_csv(path, sep='\\t')\n        cols = {c.lower(): c for c in df.columns}\n        id_col = None\n        for key in ['seq_id','protein_id','prot','id','entry','uniprot','protein']:\n            if key in cols:\n                id_col = cols[key]; break\n        term_col = None\n        for key in ['go_id','go','term','go_term']:\n            if key in cols:\n                term_col = cols[key]; break\n        if id_col is None or term_col is None:\n            raise ValueError('Missing headers, fallback to no‑header mode')\n        out = df[[id_col, term_col]].rename(columns={id_col:'seq_id', term_col:'term'})\n    except Exception:\n        df = pd.read_csv(path, sep='\\t', header=None)\n        assert df.shape[1] >= 2, \"train_terms.tsv must have ≥ 2 columns\"\n        out = df.iloc[:, :2]\n        out.columns = ['seq_id','term']\n    out = out.dropna(subset=['seq_id','term'])\n    out['seq_id'] = out['seq_id'].astype(str)\n    out['term'] = out['term'].astype(str)\n    return out\n\nterms_long = load_terms_flexible(PATH_TERMS)\nterms_long['seq_id_norm'] = terms_long['seq_id'].apply(normalize_seq_id)\n\n# Aggregate to list of GO terms per protein\nterms_df = terms_long.groupby('seq_id_norm')['term'].apply(list).reset_index()\nterms_df = terms_df.rename(columns={'seq_id_norm':'seq_id'})\nprint(\"Unique proteins in terms:\", len(terms_df))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-26T15:41:25.781034Z","iopub.execute_input":"2025-10-26T15:41:25.781353Z","iopub.status.idle":"2025-10-26T15:41:28.491784Z","shell.execute_reply.started":"2025-10-26T15:41:25.781326Z","shell.execute_reply":"2025-10-26T15:41:28.491075Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Cell 5 — Load Taxonomy (Optional) + Normalize IDs","metadata":{}},{"cell_type":"code","source":"try:\n    tax_df = pd.read_csv(PATH_TAX, sep='\\t')\nexcept Exception:\n    tax_df = pd.read_csv(PATH_TAX)\n\n# Try to detect columns\ntc = {c.lower(): c for c in tax_df.columns}\nseq_col = tc.get('seq_id') or tc.get('protein_id') or list(tax_df.columns)[0]\nif 'tax_id' in tc:\n    tax_col = tc['tax_id']\nelif 'taxonomy_id' in tc:\n    tax_col = tc['taxonomy_id']\nelse:\n    # fallback: second column\n    tax_col = list(tax_df.columns)[1]\n\n# Normalize IDs\ntax_df = tax_df[[seq_col, tax_col]].rename(columns={seq_col:'seq_id', tax_col:'tax_id'})\ntax_df['seq_id'] = tax_df['seq_id'].astype(str)\ntax_df['seq_id_norm'] = tax_df['seq_id'].apply(normalize_seq_id)\n\n# Map: normalized id -> tax_id\ntax_map = dict(zip(tax_df['seq_id_norm'], tax_df['tax_id']))\nprint(\"Taxonomy entries:\", len(tax_map))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-26T15:41:28.492519Z","iopub.execute_input":"2025-10-26T15:41:28.492795Z","iopub.status.idle":"2025-10-26T15:41:28.678267Z","shell.execute_reply.started":"2025-10-26T15:41:28.492776Z","shell.execute_reply":"2025-10-26T15:41:28.677680Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Cell 6 — Align Terms with FASTA (Normalized IDs)","metadata":{}},{"cell_type":"code","source":"before = len(terms_df)\nterms_df = terms_df[terms_df['seq_id'].isin(seqs_norm.keys())].copy()\nprint(\"After aligning with FASTA:\", terms_df.shape, f\"(from {before})\")\n\nif len(terms_df) == 0:\n    print(\"Alignment is empty. Re‑inspect headers and normalization (see Cell 2 peeks).\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-26T15:41:28.679030Z","iopub.execute_input":"2025-10-26T15:41:28.679560Z","iopub.status.idle":"2025-10-26T15:41:28.724322Z","shell.execute_reply.started":"2025-10-26T15:41:28.679509Z","shell.execute_reply":"2025-10-26T15:41:28.723575Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Cell 7 — Parse GO OBO & Build Parent Map","metadata":{}},{"cell_type":"code","source":"def parse_go_obo(path):\n    term_to_namespace = {}\n    parents = defaultdict(set)\n    with open(path, 'r', encoding='utf-8') as f:\n        cur_id, cur_ns = None, None\n        for raw in f:\n            line = raw.strip()\n            if line == '[Term]':\n                cur_id, cur_ns = None, None\n            elif line.startswith('id: GO:'):\n                cur_id = line.split('id: ')[1].strip()\n            elif line.startswith('namespace:'):\n                cur_ns = line.split('namespace:')[1].strip()\n            elif line.startswith('is_a:') and cur_id is not None:\n                parent = line.split('is_a:')[1].split('!')[0].strip()\n                parents[cur_id].add(parent)\n            elif line == '' and cur_id is not None:\n                if cur_ns:\n                    code = {'biological_process':'BP', 'molecular_function':'MF', 'cellular_component':'CC'}.get(cur_ns, cur_ns)\n                    term_to_namespace[cur_id] = code\n                cur_id, cur_ns = None, None\n    return term_to_namespace, parents\n\nterm_to_ns, go_parents = parse_go_obo(PATH_OBO)\nprint(\"GO terms parsed:\", len(term_to_ns))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-26T15:41:28.725166Z","iopub.execute_input":"2025-10-26T15:41:28.725492Z","iopub.status.idle":"2025-10-26T15:41:29.269485Z","shell.execute_reply.started":"2025-10-26T15:41:28.725465Z","shell.execute_reply":"2025-10-26T15:41:29.268852Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Cell 8 — True‑Path Propagation (Ancestors)","metadata":{}},{"cell_type":"code","source":"from functools import lru_cache\n\n@lru_cache(maxsize=None)\ndef ancestors(go_id):\n    ancs = {go_id}\n    for p in go_parents.get(go_id, []):\n        ancs |= ancestors(p)\n    return ancs\n\ndef propagate_terms(go_list):\n    out = set()\n    for t in go_list:\n        out |= ancestors(t)\n    return list(out)\n\nterms_df['propagated'] = terms_df['term'].apply(propagate_terms)\nprint(terms_df.head())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-26T15:41:29.270144Z","iopub.execute_input":"2025-10-26T15:41:29.270345Z","iopub.status.idle":"2025-10-26T15:41:29.958186Z","shell.execute_reply.started":"2025-10-26T15:41:29.270328Z","shell.execute_reply":"2025-10-26T15:41:29.957386Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Cell 9 — Namespace Filter + Frequency Pruning","metadata":{}},{"cell_type":"code","source":"# Choose namespaces to keep\nKEEP_NS = {\"MF\", \"BP\", \"CC\"}   # e.g., set to {\"MF\"} for just Molecular Function\n\n# Apply namespace filter\n\ndef filter_by_namespace(lst, keep=KEEP_NS):\n    return [t for t in lst if term_to_ns.get(t) in keep]\n\nterms_df['terms_ns'] = terms_df['propagated'].apply(filter_by_namespace)\n\n# Frequency threshold (tune). Start with small value to avoid zero‑labels.\nMIN_FREQ = 5\n\ncnt = Counter()\nfor ts in terms_df['terms_ns']:\n    cnt.update(set(ts))          # presence per protein\n\nkept_terms = {t for t,c in cnt.items() if c >= MIN_FREQ}\nprint(f\"Kept {len(kept_terms)} GO terms with freq >= {MIN_FREQ}\")\n\nterms_df['terms_kept'] = terms_df['terms_ns'].apply(lambda ts: [t for t in ts if t in kept_terms])\nterms_df = terms_df[terms_df['terms_kept'].str.len() > 0].copy()\nprint(\"Remaining proteins with labels:\", len(terms_df))\n\nif len(terms_df) == 0:\n    print(\"No proteins left. Try LOWER MIN_FREQ (e.g., 2) or restrict KEEP_NS (e.g., {'MF'}).\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-26T15:41:29.960373Z","iopub.execute_input":"2025-10-26T15:41:29.960613Z","iopub.status.idle":"2025-10-26T15:41:31.275964Z","shell.execute_reply.started":"2025-10-26T15:41:29.960595Z","shell.execute_reply":"2025-10-26T15:41:31.274920Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Cell 10 — Feature Engineering (AA composition + 3‑mer hash + taxonomy hash)","metadata":{}},{"cell_type":"code","source":"AMINO_ACIDS = list(\"ACDEFGHIKLMNPQRSTVWY\")\nAA2IDX = {aa:i for i,aa in enumerate(AMINO_ACIDS)}\n\ndef aa_composition(seq):\n    vec = np.zeros(len(AMINO_ACIDS), dtype=np.float16)\n    L = len(seq)\n    if L == 0:\n        return vec\n    for ch in seq:\n        i = AA2IDX.get(ch)\n        if i is not None:\n            vec[i] += 1.0\n    vec /= L\n    return vec\n\ndef kmer_hash(seq, k=3, buckets=512):\n    vec = np.zeros(buckets, dtype=np.float16)\n    L = len(seq)\n    if L < k:\n        return vec\n    for i in range(L-k+1):\n        km = seq[i:i+k]\n        h = 0\n        for ch in km:\n            h = (h * 27 + (ord(ch) % 27)) % (2**32)\n        vec[h % buckets] += 1.0\n    vec /= (L - k + 1)\n    return vec\n\ndef tax_hash_feature(tax_id, buckets=64):\n    try:\n        x = int(tax_id)\n    except Exception:\n        return np.zeros(buckets, dtype=np.float16)\n    rng = np.random.RandomState((x * 2654435761) % (2**32))\n    v = rng.randn(buckets).astype(np.float32)\n    v /= (np.linalg.norm(v) + 1e-8)\n    return v\n\n# Build taxonomy map already computed in Cell 5: tax_map\n\nids = []\nX_list = []\nfor sid in terms_df['seq_id']:\n    seq = seqs_norm.get(sid, \"\")\n    f1 = aa_composition(seq)\n    f2 = kmer_hash(seq, k=3, buckets=512)\n    f3 = tax_hash_feature(tax_map.get(sid, None), buckets=64)\n    X_list.append(np.concatenate([f1, f2, f3], axis=0))\n    ids.append(sid)\n\nX = np.stack(X_list, axis=0)\nprint(\"Feature matrix shape:\", X.shape)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-26T15:41:31.276793Z","iopub.execute_input":"2025-10-26T15:41:31.277348Z","iopub.status.idle":"2025-10-26T15:45:08.038095Z","shell.execute_reply.started":"2025-10-26T15:41:31.277324Z","shell.execute_reply":"2025-10-26T15:45:08.037370Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Cell 11 — Multi‑Label Binarization","metadata":{}},{"cell_type":"code","source":"id2labels = dict(zip(terms_df['seq_id'], terms_df['terms_kept']))\nlabel_lists = [id2labels[i] for i in ids]\n\nmlb = MultiLabelBinarizer(sparse_output=True)  # sparse to reduce memory\nY = mlb.fit_transform(label_lists)\nprint(\"Label matrix shape:\", Y.shape, \"sparse:\", hasattr(Y, 'tocsr'))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-26T15:47:00.364853Z","iopub.execute_input":"2025-10-26T15:47:00.365332Z","iopub.status.idle":"2025-10-26T15:47:00.975528Z","shell.execute_reply.started":"2025-10-26T15:47:00.365303Z","shell.execute_reply":"2025-10-26T15:47:00.974711Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Cell 12 — Train/Validation Split","metadata":{}},{"cell_type":"code","source":"import numpy as np\nfrom sklearn.model_selection import train_test_split\n\nidx = np.arange(X.shape[0])\n# Compute cardinality per sample (works for dense or sparse Y)\nif hasattr(Y, 'sum') and not isinstance(Y, np.ndarray):\n    cardinality = np.asarray(Y.sum(axis=1)).ravel()\nelse:\n    cardinality = Y.sum(axis=1).ravel()\n\n# Helper: check all bins have at least 2 members\ndef valid_strata(z):\n    _, counts = np.unique(z, return_counts=True)\n    return (counts >= 2).all()\n\n# Attempt 1: exact cardinality\nstrat = cardinality.copy()\nif not valid_strata(strat):\n    # Attempt 2: quantile bins\n    # Start with up to 10 bins and decrease until valid\n    nbins = min(10, max(2, int(np.unique(cardinality).size)))\n    ok = False\n    while nbins >= 2:\n        try:\n            # Use pandas qcut for quantile bins; fallback to cut if duplicates\n            import pandas as pd\n            # qcut may fail when there are too many identical values; handle with duplicates='drop'\n            try:\n                strat = pd.qcut(cardinality, q=nbins, labels=False, duplicates='drop')\n            except Exception:\n                strat = pd.cut(cardinality, bins=nbins, labels=False, include_lowest=True)\n            strat = np.asarray(strat)\n            if valid_strata(strat):\n                ok = True\n                break\n        except Exception:\n            pass\n        nbins -= 1\n    if not ok:\n        strat = None  # fall back to random split\n\nidx_tr, idx_va = train_test_split(\n    idx, test_size=0.15, random_state=42, stratify=strat\n)\nprint(\"Index splits:\", idx_tr.shape, idx_va.shape, \"| stratified:\", strat is not None)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-26T15:49:10.083789Z","iopub.execute_input":"2025-10-26T15:49:10.084375Z","iopub.status.idle":"2025-10-26T15:49:10.143248Z","shell.execute_reply.started":"2025-10-26T15:49:10.084353Z","shell.execute_reply":"2025-10-26T15:49:10.142583Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Cell 13 — Train Classifier (LogReg or MLP)","metadata":{}},{"cell_type":"code","source":"MODEL = 'mlp'  # 'logreg' or 'mlp'\n\nif MODEL == 'logreg':\n    clf = OneVsRestClassifier(\n        LogisticRegression(max_iter=200, solver='saga', random_state=42)\n    )\nelif MODEL == 'mlp':\n    clf = OneVsRestClassifier(\n        MLPClassifier(hidden_layer_sizes=(512,256), activation='relu', batch_size=128,\n                      learning_rate_init=1e-3, max_iter=50, random_state=42)\n    )\nelse:\n    raise ValueError(\"MODEL must be 'logreg' or 'mlp'\")\n\n# Slice *only when fitting* to keep memory low\nclf.fit(X[idx_tr], (Y[idx_tr] if hasattr(Y, 'tocsr') else Y[idx_tr]))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-26T15:49:14.408740Z","iopub.execute_input":"2025-10-26T15:49:14.409009Z","execution_failed":"2025-10-26T16:44:57.130Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Cell 14 — Evaluation (micro/macro + threshold sweep ≈ Fmax)","metadata":{}},{"cell_type":"code","source":"def basic_metrics(Y_true, Y_pred):\n    return {\n        'micro': {\n            'precision': float(precision_score(Y_true, Y_pred, average='micro', zero_division=0)),\n            'recall': float(recall_score(Y_true, Y_pred, average='micro', zero_division=0)),\n            'f1': float(f1_score(Y_true, Y_pred, average='micro', zero_division=0)),\n        },\n        'macro': {\n            'precision': float(precision_score(Y_true, Y_pred, average='macro', zero_division=0)),\n            'recall': float(recall_score(Y_true, Y_pred, average='macro', zero_division=0)),\n            'f1': float(f1_score(Y_true, Y_pred, average='macro', zero_division=0)),\n        }\n    }\n\nX_val = X[idx_va]\nY_val = (Y[idx_va] if hasattr(Y, 'tocsr') else Y[idx_va])\n\n# Get probabilities or decision scores\nif hasattr(clf, 'predict_proba'):\n    Yp = clf.predict_proba(X_val)\n    if isinstance(Yp, list):\n        Yp = np.vstack([p[:,1] for p in Yp]).T\nelse:\n    def sigmoid(z): return 1/(1+np.exp(-z))\n    scores = clf.decision_function(X_val)\n    if scores.ndim == 1: scores = scores.reshape(-1,1)\n    Yp = sigmoid(scores)\n\n# Sweep thresholds to approximate Fmax\ngrid = np.linspace(0.1, 0.9, 9)\nbest = {\"f1_micro\": -1, \"thr\": None}\nfor thr in grid:\n    Yhat = (Yp >= thr).astype(int)\n    m = basic_metrics(Y_val, Yhat)\n    if m['micro']['f1'] > best['f1_micro']:\n        best = {\"f1_micro\": m['micro']['f1'], \"thr\": float(thr), \"metrics\": m}\n\nprint(json.dumps(best, indent=2))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-26T15:45:08.804634Z","iopub.status.idle":"2025-10-26T15:45:08.805167Z","shell.execute_reply.started":"2025-10-26T15:45:08.805039Z","shell.execute_reply":"2025-10-26T15:45:08.805056Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Cell 15 — Save Artifacts to /kaggle/working\n","metadata":{}},{"cell_type":"code","source":"import joblib, os\nos.makedirs('/kaggle/working', exist_ok=True)\n\njoblib.dump({\n    'clf': clf,\n    'mlb': mlb,\n    'best_threshold': best.get('thr', 0.5),\n    'classes': list(mlb.classes_),\n    'config': {\n        'MODEL': MODEL,\n        'KEEP_NS': list(KEEP_NS),\n        'MIN_FREQ': MIN_FREQ\n    }\n}, '/kaggle/working/cafa6_model.joblib')\n\nwith open('/kaggle/working/val_report.json', 'w') as f:\n    json.dump(best, f, indent=2)\nwith open('/kaggle/working/label_classes.json', 'w') as f:\n    json.dump(list(mlb.classes_), f)\n\nprint(\"Saved: /kaggle/working/cafa6_model.joblib and reports\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-26T15:45:08.805806Z","iopub.status.idle":"2025-10-26T15:45:08.806044Z","shell.execute_reply.started":"2025-10-26T15:45:08.805916Z","shell.execute_reply":"2025-10-26T15:45:08.805925Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Cell 16 — Predict on Any FASTA (Using This Pipeline)","metadata":{}},{"cell_type":"code","source":"import os, io\nimport numpy as np\nimport pandas as pd\n\n# --- Config ---\nTEST_FASTA = \"/kaggle/input/cafa-6-protein-function-prediction/Test/testsuperset.fasta\"\nSUB_PATH   = \"/kaggle/working/submission.txt\"\nTHRESH     = 0.10   # do not output any (seq, term) with prob < THRESH\nBATCH_SIZE = 2048   # tune for RAM; smaller -> safer\n\n# --- Safety checks ---\nassert 'clf' in globals() and 'mlb' in globals(), \"Train the model first (Cells 13–15).\"\n\n# --- FASTA streaming reader that yields batches of (norm_id, seq) ---\n\ndef fasta_stream_batches(path, batch_size=1024):\n    buf_ids, buf_seqs = [], []\n    sid, seq_parts = None, []\n    with open(path, 'r') as fh:\n        for line in fh:\n            line = line.strip()\n            if not line:\n                continue\n            if line.startswith('>'):\n                if sid is not None:\n                    buf_ids.append(normalize_seq_id(sid))\n                    buf_seqs.append(''.join(seq_parts).upper())\n                    seq_parts = []\n                    if len(buf_ids) >= batch_size:\n                        yield buf_ids, buf_seqs\n                        buf_ids, buf_seqs = [], []\n                sid = line[1:]\n            else:\n                seq_parts.append(line)\n        # flush last record\n        if sid is not None:\n            buf_ids.append(normalize_seq_id(sid))\n            buf_seqs.append(''.join(seq_parts).upper())\n    if buf_ids:\n        yield buf_ids, buf_seqs\n\n# --- Helper: build features for a batch ---\n\ndef featurize_batch(batch_ids, batch_seqs):\n    X_list = []\n    for i, sid in enumerate(batch_ids):\n        seq = batch_seqs[i]\n        f1 = aa_composition(seq)\n        f2 = kmer_hash(seq, k=3, buckets=512)\n        # taxonomy feature if available\n        f3 = tax_hash_feature(tax_map.get(sid, None), buckets=64) if 'tax_map' in globals() else np.zeros(64, dtype=np.float32)\n        X_list.append(np.concatenate([f1, f2, f3], axis=0))\n    Xb = np.stack(X_list, axis=0)\n    return Xb\n\n# --- Prediction writer ---\n\ndef write_predictions(ids, Yprob, classes, fh, thresh=0.10):\n    # For each protein, emit only terms with prob >= thresh\n    for i, sid in enumerate(ids):\n        probs = Yprob[i]\n        idx = np.where(probs >= thresh)[0]\n        if idx.size == 0:\n            continue  # nothing above threshold -> skip\n        for j in idx:\n            term = classes[j]\n            score = float(probs[j])\n            # Write exactly 3 lines per record\n            fh.write(f\"{sid}\n\")\n            fh.write(f\"{term}\n\")\n            fh.write(f\"{score:.6f}\n\")\n\n# --- Main streaming inference loop ---\n\ntotal_seqs = 0\nemitted = 0\n\nwith open(SUB_PATH, 'w') as fout:\n    for ids_batch, seqs_batch in fasta_stream_batches(TEST_FASTA, batch_size=BATCH_SIZE):\n        total_seqs += len(ids_batch)\n        Xb = featurize_batch(ids_batch, seqs_batch)\n        # Predict probabilities (handle OneVsRest returning list)\n        if hasattr(clf, 'predict_proba'):\n            Ypb = clf.predict_proba(Xb)\n            if isinstance(Ypb, list):\n                Ypb = np.vstack([p[:, 1] for p in Ypb]).T\n        else:\n            def sigmoid(z): return 1/(1+np.exp(-z))\n            scores = clf.decision_function(Xb)\n            if scores.ndim == 1:\n                scores = scores.reshape(-1, 1)\n            Ypb = sigmoid(scores)\n        # Write thresholded predictions\n        before = fout.tell()\n        write_predictions(ids_batch, Ypb, mlb.classes_, fout, thresh=THRESH)\n        after = fout.tell()\n        if after > before:\n            # crude estimate of how many lines written: (after-before) bytes isn't reliable\n            # but we can approximate emitted pairs by counting above-threshold indices\n            emitted += int((Ypb >= THRESH).sum())\n\nprint(f\"Done. Processed {total_seqs} sequences. Emitted >= {THRESH} predictions: {emitted}\")\nprint(f\"Submission path: {SUB_PATH}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-26T15:45:08.807305Z","iopub.status.idle":"2025-10-26T15:45:08.807513Z","shell.execute_reply.started":"2025-10-26T15:45:08.807415Z","shell.execute_reply":"2025-10-26T15:45:08.807424Z"}},"outputs":[],"execution_count":null}]}