{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.11"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":116062,"databundleVersionId":14084779,"sourceType":"competition"}],"dockerImageVersionId":31154,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"id":"27b29f73-7228-4760-90ab-1d241e9a308d","cell_type":"markdown","source":"# 🧬 CAFA-6 Protein Function Prediction — CPU Baseline\n\n**Runtime:** ~1–10 min depending on `N_TRAIN`  \n**Validation:** Macro & Micro F1, Average Precision  \n**Submission:** `submission.tsv` (protein_id, GO_term, confidence)","metadata":{}},{"id":"ed7700b0-d5f4-49dc-a9a7-8f12422db4e0","cell_type":"code","source":"# ==============================\n# 📦 Imports & Config\n# ==============================\nimport os, re, time, warnings\nimport numpy as np\nimport pandas as pd\nfrom tqdm import tqdm\nfrom Bio import SeqIO\nfrom sklearn.preprocessing import MultiLabelBinarizer\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.multiclass import OneVsRestClassifier\nfrom sklearn.metrics import f1_score, average_precision_score\nfrom xgboost import XGBClassifier\n\nwarnings.filterwarnings(\"ignore\", category=UserWarning, module=\"sklearn.multiclass\")\nwarnings.filterwarnings(\"ignore\", category=UserWarning, module=\"xgboost\")\n\n# Paths\nBASE = \"/kaggle/input/cafa-6-protein-function-prediction\"\nTRAIN_TERMS = f\"{BASE}/Train/train_terms.tsv\"\nTRAIN_SEQ   = f\"{BASE}/Train/train_sequences.fasta\"\nTEST_SEQ    = f\"{BASE}/Test/testsuperset.fasta\"\n\n# Config\nRNG = 42\nN_TRAIN = 40000       # increase gradually (e.g. 20000, 40000, or None for full)\nMIN_POS = 10          # minimum positives per GO term\nTEST_LIMIT = 500      # limit for quick debug (None for full)\nTHRESH_GRID = np.round(np.arange(0.1, 0.91, 0.05), 2)\nAA = list(\"ACDEFGHIKLMNPQRSTVWY\")\n\nnp.random.seed(RNG)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-28T09:07:52.639242Z","iopub.execute_input":"2025-10-28T09:07:52.639608Z","iopub.status.idle":"2025-10-28T09:07:52.64953Z","shell.execute_reply.started":"2025-10-28T09:07:52.639584Z","shell.execute_reply":"2025-10-28T09:07:52.64809Z"}},"outputs":[],"execution_count":null},{"id":"b09ff97d-c99d-4f7d-b953-8badf2703852","cell_type":"markdown","source":"## 1️⃣ Load GO terms","metadata":{}},{"id":"cde409a7-1343-45f5-b4af-50de987ca888","cell_type":"code","source":"train_terms = pd.read_csv(TRAIN_TERMS, sep=\"\\t\")\ntrain_terms.columns = [c.strip().replace(\"\\ufeff\", \"\") for c in train_terms.columns]\ntrain_terms = train_terms.rename(columns={\"EntryID\": \"protein_id\", \"term\": \"go_term\"})\ntrain_terms = train_terms.groupby(\"protein_id\")[\"go_term\"].apply(list).reset_index()\ntrain_terms[\"go_terms\"] = train_terms[\"go_term\"].apply(lambda x: x if isinstance(x, list) else [x])\ntrain_terms = train_terms.drop(columns=[\"go_term\"])\nprint(\"✅ Grouped per-protein:\", train_terms.shape)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-28T09:07:52.651276Z","iopub.execute_input":"2025-10-28T09:07:52.651616Z","iopub.status.idle":"2025-10-28T09:07:55.346344Z","shell.execute_reply.started":"2025-10-28T09:07:52.65158Z","shell.execute_reply":"2025-10-28T09:07:55.344874Z"}},"outputs":[],"execution_count":null},{"id":"c6a33a88-d763-4f9e-afb0-f654cd9ad410","cell_type":"markdown","source":"## 2️⃣ Load FASTA files","metadata":{}},{"id":"8b677bec-496b-49a3-a443-c38189b2b3b5","cell_type":"code","source":"def extract_accession(desc):\n    if \"|\" in desc:\n        parts = desc.split(\"|\")\n        if len(parts) > 1:\n            return parts[1].strip()\n    return desc.split()[0]\n\ndef fasta_to_dict(path):\n    seqs = {}\n    for rec in tqdm(SeqIO.parse(path, \"fasta\"), desc=f\"Reading {os.path.basename(path)}\"):\n        acc = extract_accession(rec.description)\n        seqs[acc] = str(rec.seq)\n    return seqs\n\ntrain_seqs = fasta_to_dict(TRAIN_SEQ)\ntest_seqs  = fasta_to_dict(TEST_SEQ)\n\noverlap = len(set(train_terms[\"protein_id\"]) & set(train_seqs.keys()))\nprint(f\"✅ Overlap FASTA vs TSV: {overlap}\")\nprint(f\"✅ Train sequences: {len(train_seqs):,} | Test sequences: {len(test_seqs):,}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-28T09:07:55.347489Z","iopub.execute_input":"2025-10-28T09:07:55.347807Z","iopub.status.idle":"2025-10-28T09:07:57.424009Z","shell.execute_reply.started":"2025-10-28T09:07:55.347782Z","shell.execute_reply":"2025-10-28T09:07:57.422289Z"}},"outputs":[],"execution_count":null},{"id":"d48bc593-d6e8-44d7-8369-7c4604ca836b","cell_type":"markdown","source":"## 3️⃣ Encode Amino Acid Composition Features","metadata":{}},{"id":"d1ce661e-bcf6-4ee7-a63a-3bc5839fea66","cell_type":"code","source":"def seq_to_features(seq):\n    seq = seq.upper()\n    counts = np.array([seq.count(a) for a in AA], dtype=float)\n    return counts / max(len(seq), 1)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-28T09:07:57.42705Z","iopub.execute_input":"2025-10-28T09:07:57.427402Z","iopub.status.idle":"2025-10-28T09:07:57.433687Z","shell.execute_reply.started":"2025-10-28T09:07:57.427375Z","shell.execute_reply":"2025-10-28T09:07:57.43257Z"}},"outputs":[],"execution_count":null},{"id":"a3d424b1-0135-4220-8cae-4e890bc7be13","cell_type":"markdown","source":"## 4️⃣ Build Training Matrix","metadata":{}},{"id":"9dc531ec-1011-49a7-93bd-7011c3b49063","cell_type":"code","source":"if N_TRAIN is not None:\n    subset = train_terms.sample(min(N_TRAIN, len(train_terms)), random_state=RNG)\nelse:\n    subset = train_terms\n\nX, ylist, keep_ids = [], [], []\nfor pid, goterms in tqdm(zip(subset[\"protein_id\"], subset[\"go_terms\"]), total=len(subset), desc=\"Encoding training sequences\"):\n    if pid in train_seqs:\n        X.append(seq_to_features(train_seqs[pid]))\n        ylist.append(tuple(goterms))\n        keep_ids.append(pid)\n\nX = np.vstack(X)\nmlb = MultiLabelBinarizer()\ny_raw = mlb.fit_transform(ylist)\n\nmask = (y_raw.sum(axis=0) >= MIN_POS) & (y_raw.sum(axis=0) < y_raw.shape[0])\ny = y_raw[:, mask]\nmlb.classes_ = mlb.classes_[mask]\n\nprint(f\"✅ Training proteins: {X.shape[0]:,}\")\nprint(f\"✅ GO labels kept: {y.shape[1]:,}\")\n\nXtr, Xval, ytr, yval = train_test_split(X, y, test_size=0.2, random_state=RNG)\nprint(f\"➡️ Train: {Xtr.shape} | Val: {Xval.shape}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-28T09:07:57.43478Z","iopub.execute_input":"2025-10-28T09:07:57.435327Z","iopub.status.idle":"2025-10-28T09:08:24.296087Z","shell.execute_reply.started":"2025-10-28T09:07:57.435291Z","shell.execute_reply":"2025-10-28T09:08:24.294711Z"}},"outputs":[],"execution_count":null},{"id":"6df6dd81-9dc8-47be-a6f4-2fb1fb96ed20","cell_type":"markdown","source":"## 5️⃣ Train One-vs-Rest XGBoost (CPU)","metadata":{}},{"id":"c822ec23-0efb-4ca7-9828-e07e04c83049","cell_type":"code","source":"clf = OneVsRestClassifier(\n    XGBClassifier(\n        n_estimators=150,\n        learning_rate=0.1,\n        max_depth=6,\n        subsample=0.9,\n        colsample_bytree=0.9,\n        tree_method=\"hist\",\n        n_jobs=-1,\n        verbosity=0\n    ),\n    n_jobs=-1\n)\n\nprint(f\"⚙️ Training on {Xtr.shape[0]} proteins × {ytr.shape[1]} GO terms (CPU)...\")\nstart = time.time()\nclf.fit(Xtr, ytr)\nprint(f\"✅ Training done in {(time.time() - start)/60:.2f} min\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-28T09:08:24.296967Z","iopub.execute_input":"2025-10-28T09:08:24.297297Z"}},"outputs":[],"execution_count":null},{"id":"dac2c5a2-a1d8-4575-a97b-4e5ab315fb30","cell_type":"markdown","source":"## 6️⃣ Validation Metrics (Fixed `predict_proba` + AP)","metadata":{}},{"id":"4d9ab21a-742a-450c-815e-5c5ce55992f3","cell_type":"code","source":"# ✅ Safe multilabel probability stacking\nraw_proba = clf.predict_proba(Xval)\nif isinstance(raw_proba, list):\n    val_proba = np.column_stack([p[:, 1] for p in raw_proba])\nelif isinstance(raw_proba, np.ndarray):\n    val_proba = raw_proba\nelse:\n    raise TypeError(f\"Unsupported type from predict_proba: {type(raw_proba)}\")\n\nassert val_proba.shape == yval.shape, f\"Shape mismatch: {val_proba.shape} vs {yval.shape}\"\n\n# Threshold tuning\nbest_thr, best_f1 = 0.5, -1\nfor thr in THRESH_GRID:\n    preds = (val_proba >= thr).astype(int)\n    f1 = f1_score(yval, preds, average=\"macro\", zero_division=0)\n    if f1 > best_f1:\n        best_f1, best_thr = f1, thr\n\nval_pred_labels = (val_proba >= best_thr).astype(int)\nval_f1_macro = f1_score(yval, val_pred_labels, average=\"macro\", zero_division=0)\nval_f1_micro = f1_score(yval, val_pred_labels, average=\"micro\", zero_division=0)\n\nimport warnings\nwith warnings.catch_warnings():\n    warnings.simplefilter(\"ignore\", category=UserWarning)\n    val_ap_macro = average_precision_score(yval, val_proba, average=\"macro\")\n\nprint(f\"🎯 Best global threshold: {best_thr:.2f}\")\nprint(f\"✅ Validation Macro-F1: {val_f1_macro:.4f}\")\nprint(f\"✅ Validation Micro-F1: {val_f1_micro:.4f}\")\nprint(f\"✅ Validation AP (macro): {val_ap_macro:.4f}\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"3336a375-d453-4b92-a6d0-dcf6b57353ac","cell_type":"markdown","source":"## 7️⃣ Predict Test Proteins & Save Submission","metadata":{}},{"id":"f911dcbd-b646-4579-8e37-813613d72b8c","cell_type":"code","source":"if TEST_LIMIT is not None:\n    test_ids = list(test_seqs.keys())[:TEST_LIMIT]\nelse:\n    test_ids = list(test_seqs.keys())\n\nprint(f\"🧪 Encoding test proteins: {len(test_ids):,}\")\nXtest = np.vstack([seq_to_features(test_seqs[pid]) for pid in test_ids])\n\nraw_pred = clf.predict_proba(Xtest)\nif isinstance(raw_pred, list):\n    y_pred = np.column_stack([p[:, 1] for p in raw_pred])\nelse:\n    y_pred = raw_pred\n\nrows = []\nfor i, pid in enumerate(test_ids):\n    mask = y_pred[i] >= best_thr\n    if not np.any(mask):\n        continue\n    for go, conf in zip(mlb.classes_[mask], y_pred[i, mask]):\n        rows.append([pid, go, round(float(conf), 6)])\n\nsub = pd.DataFrame(rows, columns=[\"protein_id\", \"GO_term\", \"confidence\"])\nsub.to_csv(\"submission.tsv\", sep=\"\\t\", index=False)\nprint(f\"✅ submission.tsv saved: {sub.shape}\")\nsub.head()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"9a8c8ab1-9f15-48e1-884d-63b78014d27a","cell_type":"markdown","source":"---\n\n## 8️⃣ (Optional) CAFA-style Weighted Evaluation (Propagation + IA)\n\n**Requirements** (add as Kaggle datasets or local files):\n- `go-basic.obo` (Gene Ontology DAG)\n- `ia_weights.tsv` with columns: `term` and `ia` (information accretion)\n\nThis block computes a **simplified** CAFA-style weighted F1 on your validation split by:\n1) propagating predicted & ground-truth terms to ancestors\n2) weighting confusion matrix entries by `ia(term)`\n3) macro-averaging across proteins\n\nIf files are missing or `goatools` isn't installed, this cell will skip gracefully.","metadata":{}},{"id":"2ca3988c-5a0b-4898-a470-08069189fa5f","cell_type":"code","source":"try:\n    from goatools.obo_parser import GODag\n\n    # Set your dataset paths here if you uploaded them\n    GO_OBO = \"/kaggle/input/go-basic/go-basic.obo\"         # adjust if different\n    IA_FILE = \"/kaggle/input/ia-weights/ia_weights.tsv\"    # adjust if different\n\n    # Load ontology + IA weights\n    go = GODag(GO_OBO)\n    ia_df = pd.read_csv(IA_FILE, sep=\"\\t\")\n    ia = dict(zip(ia_df[\"term\"], ia_df[\"ia\"]))\n\n    def propagate_terms(terms, go):\n        expanded = set(terms)\n        for t in list(terms):\n            if t in go:\n                expanded |= {a.id for a in go[t].get_all_parents()}\n        return expanded\n\n    def weighted_confusion(pred, truth, ia):\n        tp = fp = fn = 0.0\n        for term in pred & truth:\n            tp += ia.get(term, 1.0)\n        for term in pred - truth:\n            fp += ia.get(term, 1.0)\n        for term in truth - pred:\n            fn += ia.get(term, 1.0)\n        return tp, fp, fn\n\n    def weighted_f1_ia(tp, fp, fn):\n        prec = tp / (tp + fp + 1e-12)\n        rec  = tp / (tp + fn + 1e-12)\n        if prec + rec == 0:\n            return 0.0, prec, rec\n        f1 = 2 * prec * rec / (prec + rec)\n        return f1, prec, rec\n\n    # Build per-protein dicts from validation split\n    # WARNING: here we don't know the original protein IDs for Xval directly;\n    # keep_ids was collected in the same order as X before splitting.\n    # We'll reconstruct approximate IDs for the validation set using indices.\n    # Safer approach: store indices during train_test_split.\n    # For simplicity, use index mapping from train_test_split.\n    from sklearn.model_selection import StratifiedKFold\n\n    # We'll approximate: rebuild indices by splitting indices array\n    idx_all = np.arange(X.shape[0])\n    # To preserve split, redo split with same random_state\n    Xtr_idx, Xval_idx = train_test_split(idx_all, test_size=0.2, random_state=RNG)\n\n    # Now map back to keep_ids\n    pid_val = [keep_ids[i] for i in Xval_idx]\n\n    # Build predicted/true term sets for validation proteins\n    pred_dict, truth_dict = {}, {}\n    thr = best_thr\n    for i, pid in enumerate(pid_val):\n        # predicted terms for this sample\n        preds = [mlb.classes_[j] for j in np.where(val_proba[i] >= thr)[0]]\n        trues = [mlb.classes_[j] for j in np.where(yval[i] == 1)[0]]\n        if preds:\n            pred_dict[pid] = set(preds)\n        if trues:\n            truth_dict[pid] = set(trues)\n\n    results = []\n    for pid in truth_dict.keys():\n        pset = propagate_terms(pred_dict.get(pid, set()), go)\n        tset = propagate_terms(truth_dict[pid], go)\n        tp, fp, fn = weighted_confusion(pset, tset, ia)\n        f1, prec, rec = weighted_f1_ia(tp, fp, fn)\n        results.append((f1, prec, rec))\n\n    if results:\n        f1_mean = float(np.mean([r[0] for r in results]))\n        p_mean  = float(np.mean([r[1] for r in results]))\n        r_mean  = float(np.mean([r[2] for r in results]))\n        print(f\"🎯 CAFA-style Weighted F1(ia): {f1_mean:.4f}\")\n        print(f\"✅ Weighted Precision(ia): {p_mean:.4f}\")\n        print(f\"✅ Weighted Recall(ia): {r_mean:.4f}\")\n    else:\n        print(\"No overlapping validation proteins with truth/prediction sets to evaluate.\")\n\nexcept Exception as e:\n    print(\"⚠️ Skipping CAFA-style evaluation.\", str(e))","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"09c8ea0c-b767-4fd9-9dfa-09100272c0d5","cell_type":"markdown","source":"---\n\n## ✅ Tips & Next Steps\n\n- Increase `N_TRAIN` up to `None` (full ~82k) for best baseline accuracy.\n- Adjust `MIN_POS` to control label frequency threshold (e.g., 5, 10, 20).\n- For faster experimentation, keep `TEST_LIMIT=500`; for submission, set to `None`.\n- For GPU speedups, switch XGBoost to `tree_method=\"gpu_hist\", predictor=\"gpu_predictor\"` (not needed on CPU baseline).\n- Consider adding richer features (k-mers, ESM/ProtBERT embeddings) for big gains.\n\nGood luck & happy modeling! 🚀","metadata":{}}]}