{"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":"nvidiaTeslaT4","dataSources":[{"sourceId":116062,"databundleVersionId":14084779,"sourceType":"competition"}],"dockerImageVersionId":31154,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Topics Covered in This Notebook\n\nThis notebook presents an exploratory data analysis (EDA) for the CAFA 6 Protein Function Prediction competition. It validates the dataset schema and integrity, characterizes label imbalance across GO subontologies (MF, BP, CC), examines information accumulation (IA) weights, analyzes protein sequence length distributions, investigates taxonomic skewness, and summarizes the properties of the Gene Ontology (GO) graph relevant to modeling.\n\n## Objectives\n- Establish a reliable understanding of data formats, sizes, and constraints.\n- Quantify class imbalance and its effects on evaluation and calibration.\n- Investigate IA weights to predict metric sensitivity to rare/deep terms.\n- Evaluate sequence length distributions to inform feature extraction and clustering.\n- Examine taxonomic composition, promote taxon awareness prerequisites, and validation splits.\n- Summarize the GO DAG structure (depth, relationships) to guide propagation and thresholding.\n\n## Data Sources\n- Training: sequences (`train_sequences.fasta`), labels (`train_terms.tsv`), taxonomy (`train_taxonomy.tsv`)\n- Ontology: GO graph (`go-basic.obo`)\n- Evaluation weights: knowledge base (`IA.tsv`)\n- Test superset: sequences (`testsuperset.fasta`) and taxa (`testsuperset-taxon-list.tsv`)\n- Reference: `sample_submission.tsv` for submission format\n\n## Key Questions Addressed\n- How imbalanced are GO term labels overall and per subontology?\n- How do IA weights relate to term frequency and ontology depth?\n- What are typical and extreme sequence lengths?\n- How skewed is the taxonomy, and how much train–test taxon overlap exists?\n- What ontology characteristics (depth/co-occurrence) should influence modeling choices?\n\n## Methods (High-Level)\n- Lightweight parsing of TSV/FASTA/OBO to extract schema and counts\n- Frequency and density plots for labels, IA, and sequence lengths\n- Overlap and distribution comparisons for taxonomy and IDs\n- Approximate GO depth via shortest distance from subontology roots\n\n## Outputs and Takeaways\n- Clear summary of label imbalance and IA-weighted evaluation implications\n- Practical guidance for batching/featurization from sequence length analysis\n- Rationale for taxon-aware priors and stratified validation\n- Motivation for ontology-aware prediction propagation and joint thresholding\n\n## Reproducibility\n- No comments in code cells; markdown cells explain outputs\n- Figures and tables generated directly from the provided competition data","metadata":{}},{"cell_type":"code","source":"import os\nimport gc\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom collections import Counter, defaultdict\nfrom itertools import islice","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T10:52:03.698000Z","iopub.execute_input":"2025-10-16T10:52:03.698604Z","iopub.status.idle":"2025-10-16T10:52:03.702388Z","shell.execute_reply.started":"2025-10-16T10:52:03.698580Z","shell.execute_reply":"2025-10-16T10:52:03.701595Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"base_dir = \"/kaggle/input/cafa-6-protein-function-prediction\"\ntrain_dir = os.path.join(base_dir, \"Train\")\ntest_dir = os.path.join(base_dir, \"Test\")\n\npaths = {\n    \"go_obo\": os.path.join(train_dir, \"go-basic.obo\"),\n    \"train_fasta\": os.path.join(train_dir, \"train_sequences.fasta\"),\n    \"train_terms\": os.path.join(train_dir, \"train_terms.tsv\"),\n    \"train_tax\": os.path.join(train_dir, \"train_taxonomy.tsv\"),\n    \"ia\": os.path.join(base_dir, \"IA.tsv\"),\n    \"test_fasta\": os.path.join(test_dir, \"testsuperset.fasta\"),\n    \"selected_taxa\": os.path.join(test_dir, \"testsuperset-taxon-list.tsv\"),\n    \"sample_sub\": os.path.join(base_dir, \"sample_submission.tsv\"),\n}\n\n{key: (os.path.exists(p), p) for key, p in paths.items()}","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T10:52:13.680812Z","iopub.execute_input":"2025-10-16T10:52:13.681610Z","iopub.status.idle":"2025-10-16T10:52:13.723604Z","shell.execute_reply.started":"2025-10-16T10:52:13.681587Z","shell.execute_reply":"2025-10-16T10:52:13.723014Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"sample_sub = pd.read_csv(paths[\"sample_sub\"], sep=\"\\t\", header=None, names=[\"EntryID\",\"term_or_text\",\"score\",\"text_opt\"], nrows=1000)\nsample_sub.head(8)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T10:52:23.049237Z","iopub.execute_input":"2025-10-16T10:52:23.049830Z","iopub.status.idle":"2025-10-16T10:52:23.090265Z","shell.execute_reply.started":"2025-10-16T10:52:23.049808Z","shell.execute_reply":"2025-10-16T10:52:23.089710Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Sample submission has three primary columns: protein ID, GO term ID or the literal `Text`, and a probability score in (0, 1]. Text predictions include a fourth field for the description. Scores should have up to three significant figures, and each protein can have at most 1500 terms across MF/BP/CC.","metadata":{}},{"cell_type":"code","source":"terms_df = pd.read_csv(paths[\"train_terms\"], sep=\"\\t\")\nterms_df.shape, terms_df.head(10)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T10:53:14.189251Z","iopub.execute_input":"2025-10-16T10:53:14.189570Z","iopub.status.idle":"2025-10-16T10:53:14.414010Z","shell.execute_reply.started":"2025-10-16T10:53:14.189549Z","shell.execute_reply":"2025-10-16T10:53:14.413397Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Training label table loaded. Columns are expected to be `EntryID`, `term`, and `aspect` where aspect ∈ {F, P, C} corresponding to MF, BP, and CC subontologies.","metadata":{}},{"cell_type":"code","source":"n_proteins = terms_df[\"EntryID\"].nunique()\nn_terms = terms_df[\"term\"].nunique()\naspect_counts = terms_df[\"aspect\"].value_counts().rename(index={\"F\":\"MF\",\"P\":\"BP\",\"C\":\"CC\"})\nn_rows = len(terms_df)\nsummary_basic = {\n    \"rows\": n_rows,\n    \"unique_proteins\": n_proteins,\n    \"unique_terms\": n_terms,\n    \"labels_per_protein_mean\": terms_df.groupby(\"EntryID\")[\"term\"].nunique().mean(),\n}\nsummary_basic, aspect_counts.to_dict()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T10:53:35.544561Z","iopub.execute_input":"2025-10-16T10:53:35.544825Z","iopub.status.idle":"2025-10-16T10:53:35.818200Z","shell.execute_reply.started":"2025-10-16T10:53:35.544807Z","shell.execute_reply":"2025-10-16T10:53:35.817580Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Basic summary shows the size of the multi-label problem, unique proteins and GO terms, and the count of label rows per subontology. The average number of unique terms per protein indicates label cardinality.","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(5,3))\nsns.countplot(x=terms_df[\"aspect\"].map({\"F\":\"MF\",\"P\":\"BP\",\"C\":\"CC\"}))\nplt.title(\"Label Rows by GO Subontology\")\nplt.xlabel(\"Subontology\")\nplt.ylabel(\"Count\")\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T10:54:08.564271Z","iopub.execute_input":"2025-10-16T10:54:08.564858Z","iopub.status.idle":"2025-10-16T10:54:09.115273Z","shell.execute_reply.started":"2025-10-16T10:54:08.564834Z","shell.execute_reply":"2025-10-16T10:54:09.114481Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Label distribution by subontology provides an initial view of which branch (MF/BP/CC) dominates annotations, informing sampling strategies and evaluation slicing.","metadata":{}},{"cell_type":"code","source":"term_freq = terms_df[\"term\"].value_counts()\nterm_freq_desc = term_freq.describe(percentiles=[0.5,0.9,0.99]).to_dict()\nterm_freq_desc","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T10:54:40.058769Z","iopub.execute_input":"2025-10-16T10:54:40.059033Z","iopub.status.idle":"2025-10-16T10:54:40.109614Z","shell.execute_reply.started":"2025-10-16T10:54:40.059013Z","shell.execute_reply":"2025-10-16T10:54:40.108933Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Term frequency statistics quantify class imbalance and long-tail behavior. High skew suggests that naive baselines (global priors) may overfit head classes while underrepresenting rare terms.","metadata":{}},{"cell_type":"code","source":"top_terms = term_freq.head(20)\nplt.figure(figsize=(7,4))\nsns.barplot(x=top_terms.values, y=top_terms.index, orient=\"h\")\nplt.title(\"Top 20 Most Frequent GO Terms\")\nplt.xlabel(\"Frequency\")\nplt.ylabel(\"GO Term\")\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T10:55:14.377679Z","iopub.execute_input":"2025-10-16T10:55:14.377950Z","iopub.status.idle":"2025-10-16T10:55:14.758000Z","shell.execute_reply.started":"2025-10-16T10:55:14.377931Z","shell.execute_reply":"2025-10-16T10:55:14.757388Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"The top terms contribute a sizable fraction of labels. Downstream models should control for dominance of high-frequency terms to avoid inflated precision on shallow ontology nodes.","metadata":{}},{"cell_type":"code","source":"ia_df = pd.read_csv(paths[\"ia\"], sep=\"\\t\", header=None, names=[\"term\",\"ia\"])\nia_df[\"ia\"].describe(percentiles=[0.5,0.9,0.99]).to_dict()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T10:55:52.179287Z","iopub.execute_input":"2025-10-16T10:55:52.179590Z","iopub.status.idle":"2025-10-16T10:55:52.224686Z","shell.execute_reply.started":"2025-10-16T10:55:52.179569Z","shell.execute_reply":"2025-10-16T10:55:52.223880Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Information accretion weights characterize the relative importance of terms. The long tail at higher IA values often corresponds to deeper, rarer, and harder-to-predict terms, which are emphasized in weighted metrics.","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(6,3.8))\nsns.histplot(ia_df[\"ia\"], bins=80, log_scale=(False, False))\nplt.title(\"Distribution of IA Weights\")\nplt.xlabel(\"IA\")\nplt.ylabel(\"Count\")\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T10:56:07.269036Z","iopub.execute_input":"2025-10-16T10:56:07.269573Z","iopub.status.idle":"2025-10-16T10:56:07.559115Z","shell.execute_reply.started":"2025-10-16T10:56:07.269549Z","shell.execute_reply":"2025-10-16T10:56:07.558410Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"IA weight distribution helps anticipate evaluation sensitivity to rare, specific terms. Calibration and thresholding should account for IA to optimize the weighted F1.","metadata":{}},{"cell_type":"code","source":"# Efficient FASTA header parsing and sequence length distribution\ndef fasta_lengths(path, max_records=None):\n    lengths = []\n    ids = []\n    cnt = 0\n    with open(path, \"r\", encoding=\"utf-8\", errors=\"ignore\") as f:\n        seq_len = 0\n        cur_id = None\n        for line in f:\n            if line.startswith(\">\"):\n                if cur_id is not None:\n                    lengths.append(seq_len)\n                    ids.append(cur_id)\n                cur_id = line.strip().split(\"|\")[1] if \"|\" in line else line.strip()[1:].split()[0]\n                seq_len = 0\n                cnt += 1\n                if max_records and cnt > max_records:\n                    break\n            else:\n                seq_len += len(line.strip())\n        if cur_id is not None:\n            lengths.append(seq_len)\n            ids.append(cur_id)\n    return pd.DataFrame({\"EntryID\": ids, \"seq_len\": lengths})\n\ntrain_len_df = fasta_lengths(paths[\"train_fasta\"])\ntrain_len_df[\"seq_len\"].describe(percentiles=[0.5,0.9,0.99]).to_dict()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T10:56:22.199679Z","iopub.execute_input":"2025-10-16T10:56:22.200331Z","iopub.status.idle":"2025-10-16T10:56:23.170407Z","shell.execute_reply.started":"2025-10-16T10:56:22.200309Z","shell.execute_reply":"2025-10-16T10:56:23.169782Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Sequence length summary offers a sense of the modeling space for sequence encoders and similarity tools. Extreme lengths may require batching strategies or truncated embeddings.","metadata":{}},{"cell_type":"code","source":"tax_df = pd.read_csv(paths[\"train_tax\"], sep=\"\\t\", header=None, names=[\"EntryID\",\"taxon_id\"])\ntax_counts = tax_df[\"taxon_id\"].value_counts()\ntax_counts.head(10).to_dict()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T10:56:35.233713Z","iopub.execute_input":"2025-10-16T10:56:35.234227Z","iopub.status.idle":"2025-10-16T10:56:35.291109Z","shell.execute_reply.started":"2025-10-16T10:56:35.234207Z","shell.execute_reply":"2025-10-16T10:56:35.290506Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Top taxa by protein counts indicate potential benefits of taxon-aware priors or stratified evaluation splits. Species imbalance may bias learned representations.","metadata":{}},{"cell_type":"code","source":"top_taxon = tax_counts.head(20)\nplt.figure(figsize=(7,4))\nsns.barplot(x=top_taxon.values, y=top_taxon.index, orient=\"h\")\nplt.title(\"Top 20 Taxa by Training Proteins\")\nplt.xlabel(\"Count\")\nplt.ylabel(\"Taxon ID\")\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T10:56:50.668341Z","iopub.execute_input":"2025-10-16T10:56:50.668860Z","iopub.status.idle":"2025-10-16T10:56:50.927676Z","shell.execute_reply.started":"2025-10-16T10:56:50.668838Z","shell.execute_reply":"2025-10-16T10:56:50.926989Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Dominant taxa may warrant per-taxon calibration or separate baselines to avoid conflating phylogenetic signal with functional signal.","metadata":{}},{"cell_type":"code","source":"labels_per_protein = terms_df.groupby(\"EntryID\")[\"term\"].nunique()\nlabels_per_protein.describe(percentiles=[0.5,0.9,0.99]).to_dict()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T10:57:01.814603Z","iopub.execute_input":"2025-10-16T10:57:01.815150Z","iopub.status.idle":"2025-10-16T10:57:02.002965Z","shell.execute_reply.started":"2025-10-16T10:57:01.815128Z","shell.execute_reply":"2025-10-16T10:57:02.002272Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Label cardinality per protein summarizes how many GO terms are typically associated with a single protein. This impacts thresholding and the expected sparsity of predictions.","metadata":{}},{"cell_type":"code","source":"sample_entries = set(terms_df[\"EntryID\"].drop_duplicates().sample(5000, random_state=42)) if terms_df[\"EntryID\"].nunique() > 5000 else set(terms_df[\"EntryID\"].unique())\nsub = terms_df[terms_df[\"EntryID\"].isin(sample_entries)]\npair_counts = Counter()\nfor eid, grp in sub.groupby(\"EntryID\"):\n    terms = grp[\"term\"].values\n    if len(terms) <= 1:\n        continue\n    for i in range(len(terms)):\n        for j in range(i+1, len(terms)):\n            a, b = terms[i], terms[j]\n            if a > b:\n                a, b = b, a\n            pair_counts[(a,b)] += 1\nco_top = sorted(pair_counts.items(), key=lambda x: x[1], reverse=True)[:20]\nco_top[:10]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T10:57:26.068254Z","iopub.execute_input":"2025-10-16T10:57:26.068857Z","iopub.status.idle":"2025-10-16T10:57:26.589464Z","shell.execute_reply.started":"2025-10-16T10:57:26.068833Z","shell.execute_reply":"2025-10-16T10:57:26.588809Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Frequent co-occurring terms reveal ontology or biological relationships commonly annotated together. This can motivate structured predictors or joint calibration strategies.","metadata":{}},{"cell_type":"code","source":"freq_df = term_freq.rename(\"freq\").reset_index().rename(columns={\"index\":\"term\"})\nia_join = freq_df.merge(ia_df, on=\"term\", how=\"left\").dropna()\ncorr = ia_join[[\"freq\",\"ia\"]].corr().iloc[0,1]\ncorr","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T10:57:45.011759Z","iopub.execute_input":"2025-10-16T10:57:45.012429Z","iopub.status.idle":"2025-10-16T10:57:45.047150Z","shell.execute_reply.started":"2025-10-16T10:57:45.012405Z","shell.execute_reply":"2025-10-16T10:57:45.046490Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Correlation between term frequency and IA weights typically trends negative: rarer (deeper) terms tend to have higher IA. This highlights the need to avoid shallow-term bias.","metadata":{}},{"cell_type":"code","source":"test_ids = []\nwith open(paths[\"test_fasta\"], \"r\", encoding=\"utf-8\", errors=\"ignore\") as f:\n    for line in f:\n        if line.startswith(\">\"):\n            pid = line.strip().split(\"|\")[1] if \"|\" in line else line.strip()[1:].split()[0]\n            test_ids.append(pid)\ntest_ids = pd.Series(test_ids, name=\"EntryID\")\ntrain_ids = terms_df[\"EntryID\"].drop_duplicates()\noverlap = len(set(train_ids).intersection(set(test_ids)))\ntotal_test = test_ids.nunique()\n{\"test_unique\": total_test, \"overlap_with_train\": overlap, \"overlap_rate\": round(overlap / max(1,total_test), 4)}","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T10:57:56.758045Z","iopub.execute_input":"2025-10-16T10:57:56.758326Z","iopub.status.idle":"2025-10-16T10:57:59.621943Z","shell.execute_reply.started":"2025-10-16T10:57:56.758307Z","shell.execute_reply":"2025-10-16T10:57:59.621295Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Training–test superset overlap indicates how often the same UniProt IDs appear across splits. Although evaluations are prospective, overlap patterns can influence similarity-based baselines.","metadata":{}},{"cell_type":"code","source":"def parse_obo_ids(path, max_nodes=None):\n    parents = defaultdict(list)\n    aspects = {}\n    cur_id = None\n    cur_ns = None\n    cnt = 0\n    with open(path, \"r\", encoding=\"utf-8\", errors=\"ignore\") as f:\n        for line in f:\n            s = line.strip()\n            if s == \"[Term]\":\n                if cur_id is not None:\n                    cnt += 1\n                    if max_nodes and cnt > max_nodes:\n                        break\n                cur_id = None\n                cur_ns = None\n            elif s.startswith(\"id: GO:\"):\n                cur_id = s.split(\"id: \")[1]\n            elif s.startswith(\"namespace:\"):\n                cur_ns = s.split(\"namespace: \")[1]\n                if cur_id:\n                    aspects[cur_id] = cur_ns\n            elif s.startswith(\"is_a: GO:\") and cur_id:\n                p = s.split(\"is_a: \")[1].split()[0]\n                parents[cur_id].append(p)\n            elif s.startswith(\"relationship: part_of GO:\") and cur_id:\n                p = s.split(\"relationship: part_of \")[1].split()[0]\n                parents[cur_id].append(p)\n    return parents, aspects\n\nparents, aspects_ns = parse_obo_ids(paths[\"go_obo\"])\nnum_nodes = len(aspects_ns)\nnum_edges = sum(len(v) for v in parents.values())\n{\"nodes\": num_nodes, \"edges\": num_edges}","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T10:58:12.897004Z","iopub.execute_input":"2025-10-16T10:58:12.897753Z","iopub.status.idle":"2025-10-16T10:58:13.831526Z","shell.execute_reply.started":"2025-10-16T10:58:12.897727Z","shell.execute_reply":"2025-10-16T10:58:13.830848Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"GO graph summary provides scale of the DAG and supports simple diagnostics without pulling full ontology toolchains. Relationships include `is_a` and `part_of`.","metadata":{}},{"cell_type":"code","source":"from collections import deque\n\nroots = {\n    \"biological_process\": \"GO:0008150\",\n    \"cellular_component\": \"GO:0005575\",\n    \"molecular_function\": \"GO:0003674\",\n}\n\ndef shortest_depth_from_roots(parents, roots_set):\n    depth = {}\n    indeg = defaultdict(int)\n    children = defaultdict(list)\n    for child, pars in parents.items():\n        for p in pars:\n            children[p].append(child)\n            indeg[child] += 1\n        if child not in indeg:\n            indeg[child] = indeg.get(child, 0)\n    q = deque()\n    for r in roots_set:\n        depth[r] = 0\n        q.append(r)\n    while q:\n        u = q.popleft()\n        for v in children.get(u, []):\n            dv = depth.get(v, None)\n            nu = depth[u] + 1\n            if dv is None or nu < dv:\n                depth[v] = nu\n            indeg[v] -= 1\n            if indeg[v] <= 0:\n                q.append(v)\n    return depth\n\ndepth_bp = shortest_depth_from_roots(parents, {roots[\"biological_process\"]})\ndepth_cc = shortest_depth_from_roots(parents, {roots[\"cellular_component\"]})\ndepth_mf = shortest_depth_from_roots(parents, {roots[\"molecular_function\"]})\n\ndepth_stats = {\n    \"BP_median_depth\": float(pd.Series([d for k,d in depth_bp.items()]).median()) if depth_bp else None,\n    \"CC_median_depth\": float(pd.Series([d for k,d in depth_cc.items()]).median()) if depth_cc else None,\n    \"MF_median_depth\": float(pd.Series([d for k,d in depth_mf.items()]).median()) if depth_mf else None,\n}\ndepth_stats","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T10:58:26.882072Z","iopub.execute_input":"2025-10-16T10:58:26.882350Z","iopub.status.idle":"2025-10-16T10:58:27.020499Z","shell.execute_reply.started":"2025-10-16T10:58:26.882330Z","shell.execute_reply":"2025-10-16T10:58:27.019591Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Approximate depth distributions per subontology indicate hierarchical granularity. Deeper terms are typically rarer, with higher IA, and drive weighted evaluation sensitivity.","metadata":{}},{"cell_type":"code","source":"del sample_sub\ngc.collect()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T10:58:41.739547Z","iopub.execute_input":"2025-10-16T10:58:41.740140Z","iopub.status.idle":"2025-10-16T10:58:41.829667Z","shell.execute_reply.started":"2025-10-16T10:58:41.740117Z","shell.execute_reply":"2025-10-16T10:58:41.829009Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## EDA Summary and Modeling Implications\n\n- Significant class imbalance with a heavy head and a long tail of rare terms.\n- IA prioritizes deeper, rarer GO terms; calibration should consider IA-weighted objectives.\n- Wide sequence length variation suggests careful batching and featurization choices.\n- Taxonomic skew indicates benefits from taxon-aware priors and stratified validation.\n- GO DAG depth and co-occurrence patterns motivate ancestor propagation and joint thresholding.","metadata":{}},{"cell_type":"code","source":"aspect_map = {\"F\":\"MF\", \"P\":\"BP\", \"C\":\"CC\"}\nterms_df[\"_aspect_name\"] = terms_df[\"aspect\"].map(aspect_map)\n\nfreq_by_aspect = (\n    terms_df.groupby(\"_aspect_name\")[\"term\"]\n    .value_counts()\n    .rename(\"freq\")\n    .reset_index()\n)\n\ng = sns.FacetGrid(freq_by_aspect, col=\"_aspect_name\", sharex=False, sharey=False, height=3.2)\ng.map_dataframe(sns.histplot, x=\"freq\", bins=80)\ng.set_axis_labels(\"Term frequency\", \"Count\")\ng.set_titles(\"{col_name}\")\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T11:05:08.528644Z","iopub.execute_input":"2025-10-16T11:05:08.529187Z","iopub.status.idle":"2025-10-16T11:05:09.588890Z","shell.execute_reply.started":"2025-10-16T11:05:08.529162Z","shell.execute_reply":"2025-10-16T11:05:09.588263Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Aspect-wise term frequency histograms illustrate class imbalance per subontology. Each subontology shows a heavy head with a long tail of rare terms, implying that naive global priors will be insufficient without balancing or calibration.","metadata":{}},{"cell_type":"code","source":"freq_df = terms_df[\"term\"].value_counts().rename(\"freq\").reset_index().rename(columns={\"index\":\"term\"})\nia_join = freq_df.merge(ia_df, on=\"term\", how=\"left\").dropna()\n\ncorr_pearson = ia_join[[\"freq\",\"ia\"]].corr(method=\"pearson\").iloc[0,1]\ncorr_spearman = ia_join[[\"freq\",\"ia\"]].corr(method=\"spearman\").iloc[0,1]\n{\"pearson_freq_IA\": float(corr_pearson), \"spearman_freq_IA\": float(corr_spearman)}","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T11:05:33.540794Z","iopub.execute_input":"2025-10-16T11:05:33.541346Z","iopub.status.idle":"2025-10-16T11:05:33.623511Z","shell.execute_reply.started":"2025-10-16T11:05:33.541325Z","shell.execute_reply":"2025-10-16T11:05:33.622749Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"IA generally increases as terms get rarer (negative correlation with frequency). This reinforces the need to optimize for IA-weighted objectives or apply thresholds that avoid shallow-term bias.","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(6,4))\nsns.scatterplot(data=ia_join.sample(min(50000, len(ia_join)), random_state=42), x=\"freq\", y=\"ia\", alpha=0.25, s=10)\nplt.xscale(\"log\")\nplt.title(\"IA vs Term Frequency\")\nplt.xlabel(\"Term frequency (log scale)\")\nplt.ylabel(\"IA\")\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T11:05:46.271987Z","iopub.execute_input":"2025-10-16T11:05:46.272513Z","iopub.status.idle":"2025-10-16T11:05:46.795075Z","shell.execute_reply.started":"2025-10-16T11:05:46.272489Z","shell.execute_reply":"2025-10-16T11:05:46.794385Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"The IA–frequency scatter (log-scaled x-axis) highlights that many high-IA terms are very rare. Models should be attentive to recall on these rare but important terms.","metadata":{}},{"cell_type":"code","source":"ns_map = {\n    \"biological_process\": \"BP\",\n    \"cellular_component\": \"CC\",\n    \"molecular_function\": \"MF\",\n}\n\ndef term_depth(go_id):\n    if go_id in depth_bp:\n        return depth_bp.get(go_id, np.nan)\n    if go_id in depth_cc:\n        return depth_cc.get(go_id, np.nan)\n    if go_id in depth_mf:\n        return depth_mf.get(go_id, np.nan)\n    return np.nan\n\nunique_terms = terms_df[\"term\"].drop_duplicates()\ndepth_df = pd.DataFrame({\"term\": unique_terms})\ndepth_df[\"depth\"] = depth_df[\"term\"].map(term_depth)\n\naspect_ns_df = pd.DataFrame({\"term\": list(aspects_ns.keys()), \"ns\": [aspects_ns[k] for k in aspects_ns.keys()]})\naspect_ns_df[\"aspect\"] = aspect_ns_df[\"ns\"].map(ns_map)\n\ndepth_ns = depth_df.merge(aspect_ns_df[[\"term\",\"aspect\"]], on=\"term\", how=\"left\")\ndepth_ns = depth_ns.dropna(subset=[\"depth\",\"aspect\"])\ndepth_ns.head(3)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T11:06:15.966345Z","iopub.execute_input":"2025-10-16T11:06:15.967117Z","iopub.status.idle":"2025-10-16T11:06:16.069090Z","shell.execute_reply.started":"2025-10-16T11:06:15.967093Z","shell.execute_reply":"2025-10-16T11:06:16.068279Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"We derive an approximate depth for labeled terms by shortest distance from each subontology root. This provides a coarse measure of hierarchical specificity per term.","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(7,3.5))\nsns.kdeplot(data=depth_ns, x=\"depth\", hue=\"aspect\", common_norm=False, fill=True)\nplt.title(\"Approximate Term Depth by Subontology\")\nplt.xlabel(\"Depth from subontology root (approx.)\")\nplt.ylabel(\"Density\")\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T11:06:39.999912Z","iopub.execute_input":"2025-10-16T11:06:40.000692Z","iopub.status.idle":"2025-10-16T11:06:40.499182Z","shell.execute_reply.started":"2025-10-16T11:06:40.000667Z","shell.execute_reply":"2025-10-16T11:06:40.498492Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Depth distributions differ across MF/BP/CC. Deeper terms are more specific and tend to carry higher IA, which shifts evaluation sensitivity toward accurate predictions in the tail.","metadata":{}},{"cell_type":"code","source":"ia_depth = ia_join.merge(depth_ns[[\"term\",\"depth\"]], on=\"term\", how=\"left\").dropna()\nsample = ia_depth.sample(min(100000, len(ia_depth)), random_state=42)\n\nplt.figure(figsize=(6.2,4.2))\nplt.hexbin(sample[\"depth\"], sample[\"ia\"], gridsize=50, cmap=\"viridis\")\ncb = plt.colorbar()\ncb.set_label(\"Counts\")\nplt.xlabel(\"Depth\")\nplt.ylabel(\"IA\")\nplt.title(\"IA vs Depth (Hexbin)\")\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T11:06:57.472834Z","iopub.execute_input":"2025-10-16T11:06:57.473106Z","iopub.status.idle":"2025-10-16T11:06:57.757649Z","shell.execute_reply.started":"2025-10-16T11:06:57.473085Z","shell.execute_reply":"2025-10-16T11:06:57.756826Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"IA grows with depth, confirming that deeper ontology nodes are weighted more heavily. This motivates ancestor propagation and careful score calibration to avoid over-penalizing shallow-only predictions.","metadata":{}},{"cell_type":"code","source":"entry_aspects = terms_df.groupby(\"EntryID\")[\"_aspect_name\"].unique().reset_index()\nentry_aspects = entry_aspects.explode(\"_aspect_name\").rename(columns={\"_aspect_name\":\"aspect\"})\nlen_by_aspect = train_len_df.merge(entry_aspects, on=\"EntryID\", how=\"left\").dropna(subset=[\"aspect\"])\n\nplt.figure(figsize=(6.5,4))\nsns.boxplot(data=len_by_aspect, x=\"aspect\", y=\"seq_len\", showfliers=False)\nplt.title(\"Sequence Length by Subontology Presence\")\nplt.xlabel(\"Aspect\")\nplt.ylabel(\"Sequence length (aa)\")\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T11:07:13.953950Z","iopub.execute_input":"2025-10-16T11:07:13.954415Z","iopub.status.idle":"2025-10-16T11:07:17.209962Z","shell.execute_reply.started":"2025-10-16T11:07:13.954389Z","shell.execute_reply":"2025-10-16T11:07:17.209320Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Sequence lengths vary across proteins associated with different subontologies. While not causal, this can affect batching and memory usage for sequence encoders.","metadata":{}},{"cell_type":"code","source":"test_taxa = pd.read_csv(paths[\"selected_taxa\"], sep=\"\\t\", header=None, names=[\"taxon_id\"])\ntrain_tax_counts = tax_df[\"taxon_id\"].value_counts()\ntest_tax_counts = test_taxa[\"taxon_id\"].value_counts()\n\noverlap_taxa = set(train_tax_counts.index).intersection(set(test_tax_counts.index))\noverlap_stats = {\n    \"train_taxa_unique\": int(train_tax_counts.size),\n    \"test_taxa_unique\": int(test_tax_counts.size),\n    \"overlap_taxa\": int(len(overlap_taxa)),\n    \"overlap_rate\": round(len(overlap_taxa) / max(1, test_tax_counts.size), 4),\n}\noverlap_stats","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T11:07:35.936902Z","iopub.execute_input":"2025-10-16T11:07:35.937639Z","iopub.status.idle":"2025-10-16T11:07:35.986204Z","shell.execute_reply.started":"2025-10-16T11:07:35.937616Z","shell.execute_reply":"2025-10-16T11:07:35.985672Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Taxonomic overlap between train and test-superset indicates potential transferability of taxon-aware priors and nearest-neighbor baselines. High overlap supports leveraging taxonomy-conditioned statistics.","metadata":{}},{"cell_type":"code","source":"top_train_tax = train_tax_counts.head(20)\ntop_test_tax = test_tax_counts.head(20)\n\nfig, axes = plt.subplots(1, 2, figsize=(12,4), sharey=True)\nsns.barplot(x=top_train_tax.values, y=top_train_tax.index, ax=axes[0], orient=\"h\")\naxes[0].set_title(\"Top 20 Train Taxa\")\naxes[0].set_xlabel(\"Count\")\naxes[0].set_ylabel(\"Taxon ID\")\n\nsns.barplot(x=top_test_tax.values, y=top_test_tax.index, ax=axes[1], orient=\"h\")\naxes[1].set_title(\"Top 20 Test-superset Taxa\")\naxes[1].set_xlabel(\"Count\")\naxes[1].set_ylabel(\"Taxon ID\")\n\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T11:08:39.172690Z","iopub.execute_input":"2025-10-16T11:08:39.173254Z","iopub.status.idle":"2025-10-16T11:08:39.663284Z","shell.execute_reply.started":"2025-10-16T11:08:39.173235Z","shell.execute_reply":"2025-10-16T11:08:39.662571Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Comparing top taxa surfaces distribution shifts. If the test-superset emphasizes different taxa, priors and similarity searches should account for this to maintain generalization.","metadata":{}},{"cell_type":"code","source":"terms_per_protein_aspect = (\n    terms_df.groupby([\"EntryID\",\"_aspect_name\"])[\"term\"].nunique().rename(\"n_terms\").reset_index()\n)\n\nplt.figure(figsize=(6.5,4))\nsns.violinplot(data=terms_per_protein_aspect, x=\"_aspect_name\", y=\"n_terms\", inner=\"quartile\", cut=0)\nplt.title(\"Number of Terms per Protein by Aspect\")\nplt.xlabel(\"Aspect\")\nplt.ylabel(\"Distinct GO terms\")\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T11:09:22.250986Z","iopub.execute_input":"2025-10-16T11:09:22.251581Z","iopub.status.idle":"2025-10-16T11:09:23.336620Z","shell.execute_reply.started":"2025-10-16T11:09:22.251559Z","shell.execute_reply":"2025-10-16T11:09:23.335869Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Proteins often carry multiple terms per subontology. This density guides expectations for the number of predictions per protein and informs per-aspect thresholding.","metadata":{}},{"cell_type":"code","source":"missing_ia = set(terms_df[\"term\"].unique()) - set(ia_df[\"term\"].unique())\nlen(missing_ia)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T11:09:34.941002Z","iopub.execute_input":"2025-10-16T11:09:34.941291Z","iopub.status.idle":"2025-10-16T11:09:34.998911Z","shell.execute_reply.started":"2025-10-16T11:09:34.941271Z","shell.execute_reply":"2025-10-16T11:09:34.998314Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"A small number of labeled terms may not have IA values due to ontology version or filtering. During submission assembly, ensure that parent propagation and scoring handle such cases gracefully.","metadata":{}}]}