{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"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":31192,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"id":"c2e60daa280e2cf0","cell_type":"code","source":"!pip install Biopython","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-14T10:34:50.287641Z","iopub.execute_input":"2025-11-14T10:34:50.287986Z","iopub.status.idle":"2025-11-14T10:34:57.144120Z","shell.execute_reply.started":"2025-11-14T10:34:50.287961Z","shell.execute_reply":"2025-11-14T10:34:57.143183Z"}},"outputs":[],"execution_count":null},{"id":"a3dd16b7-e374-410a-aede-4290040c12ec","cell_type":"code","source":"import os\nimport gc\nfrom collections import defaultdict, Counter\nimport pandas as pd\nimport numpy as np\nfrom tqdm.auto import tqdm\nfrom Bio import SeqIO\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-14T10:35:11.496478Z","iopub.execute_input":"2025-11-14T10:35:11.497440Z","iopub.status.idle":"2025-11-14T10:35:12.463962Z","shell.execute_reply.started":"2025-11-14T10:35:11.497404Z","shell.execute_reply":"2025-11-14T10:35:12.463009Z"}},"outputs":[],"execution_count":null},{"id":"dd8d5914d399e553","cell_type":"code","source":"# =============================================================================\n# CONFIGURATION\n# =============================================================================\nBASE_PATH = '/kaggle/input/cafa-6-protein-function-prediction'\nTRAIN_DIR = f'{BASE_PATH}/Train'\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-14T10:35:14.670176Z","iopub.execute_input":"2025-11-14T10:35:14.670695Z","iopub.status.idle":"2025-11-14T10:35:14.676224Z","shell.execute_reply.started":"2025-11-14T10:35:14.670667Z","shell.execute_reply":"2025-11-14T10:35:14.674892Z"}},"outputs":[],"execution_count":null},{"id":"794e7718cf2a85a6","cell_type":"code","source":"# =============================================================================\n# 1. PARSE GO ONTOLOGY (Enhanced with namespace extraction)\n# =============================================================================\ndef parse_obo_detailed(go_obo_path):\n    \"\"\"Parse OBO file and extract detailed graph structure and metadata\"\"\"\n    parents = defaultdict(set)\n    children = defaultdict(set)\n    term_to_namespace = {}  # MF, BP, or CC\n    term_to_name = {}\n    obsolete_terms = set()\n\n    if not os.path.exists(go_obo_path):\n        return parents, children, term_to_namespace, term_to_name, obsolete_terms\n\n    with open(go_obo_path, \"r\") as f:\n        cur_id = None\n        cur_namespace = None\n        cur_name = None\n        is_obsolete = False\n\n        for line in f:\n            line = line.strip()\n\n            if line == \"[Term]\":\n                # Save previous term if it exists\n                if cur_id and not is_obsolete:\n                    term_to_namespace[cur_id] = cur_namespace\n                    term_to_name[cur_id] = cur_name\n                elif cur_id and is_obsolete:\n                    obsolete_terms.add(cur_id)\n\n                # Reset for new term\n                cur_id = None\n                cur_namespace = None\n                cur_name = None\n                is_obsolete = False\n\n            elif line.startswith(\"id: \"):\n                cur_id = line.split(\"id: \")[1].strip()\n\n            elif line.startswith(\"name: \"):\n                cur_name = line.split(\"name: \")[1].strip()\n\n            elif line.startswith(\"namespace: \"):\n                ns = line.split(\"namespace: \")[1].strip()\n                if ns == \"molecular_function\":\n                    cur_namespace = \"MF\"\n                elif ns == \"biological_process\":\n                    cur_namespace = \"BP\"\n                elif ns == \"cellular_component\":\n                    cur_namespace = \"CC\"\n\n            elif line.startswith(\"is_obsolete: true\"):\n                is_obsolete = True\n\n            elif line.startswith(\"is_a: \"):\n                pid = line.split()[1].strip()\n                if cur_id and not is_obsolete:\n                    parents[cur_id].add(pid)\n                    children[pid].add(cur_id)\n\n            elif line.startswith(\"relationship: part_of \"):\n                parts = line.split()\n                if len(parts) >= 3:\n                    pid = parts[2].strip()\n                    if cur_id and not is_obsolete:\n                        parents[cur_id].add(pid)\n                        children[pid].add(cur_id)\n\n    print(f\"✓ Parsed OBO: {len(parents)} terms with parents\")\n    print(f\"✓ Namespaces: MF={sum(1 for ns in term_to_namespace.values() if ns == 'MF')}, \"\n          f\"BP={sum(1 for ns in term_to_namespace.values() if ns == 'BP')}, \"\n          f\"CC={sum(1 for ns in term_to_namespace.values() if ns == 'CC')}\")\n    print(f\"✓ Obsolete terms excluded: {len(obsolete_terms)}\")\n\n    return parents, children, term_to_namespace, term_to_name, obsolete_terms\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-14T10:35:16.610119Z","iopub.execute_input":"2025-11-14T10:35:16.610438Z","iopub.status.idle":"2025-11-14T10:35:16.621586Z","shell.execute_reply.started":"2025-11-14T10:35:16.610414Z","shell.execute_reply":"2025-11-14T10:35:16.620777Z"}},"outputs":[],"execution_count":null},{"id":"d9f609777b2408f3","cell_type":"code","source":"# =============================================================================\n# 2. DAG ANALYSIS FUNCTIONS\n# =============================================================================\ndef get_roots(parents, term_to_namespace):\n    \"\"\"Find root nodes (nodes with no parents) per ontology\"\"\"\n    roots = defaultdict(list)\n    all_terms = set(parents.keys()) | set(term_to_namespace.keys())\n\n    for term in all_terms:\n        if term not in parents or len(parents[term]) == 0:\n            namespace = term_to_namespace.get(term, 'Unknown')\n            roots[namespace].append(term)\n\n    return roots\n\n\ndef get_leaves(children, term_to_namespace):\n    \"\"\"Find leaf nodes (nodes with no children) per ontology\"\"\"\n    leaves = defaultdict(list)\n    all_terms = set(children.keys()) | set(term_to_namespace.keys())\n\n    for term in all_terms:\n        if term not in children or len(children[term]) == 0:\n            namespace = term_to_namespace.get(term, 'Unknown')\n            leaves[namespace].append(term)\n\n    return leaves\n\n\ndef compute_depth(term, parents, memo=None):\n    \"\"\"Compute maximum depth of term in DAG (distance to root)\"\"\"\n    if memo is None:\n        memo = {}\n\n    if term in memo:\n        return memo[term]\n\n    if term not in parents or len(parents[term]) == 0:\n        memo[term] = 0\n        return 0\n\n    max_depth = max(compute_depth(p, parents, memo) for p in parents[term]) + 1\n    memo[term] = max_depth\n    return max_depth\n\n\ndef analyze_dag_structure(parents, children, term_to_namespace):\n    \"\"\"Comprehensive DAG analysis per ontology\"\"\"\n    stats = {}\n\n    for ont in ['MF', 'BP', 'CC']:\n        ont_terms = [t for t, ns in term_to_namespace.items() if ns == ont]\n\n        if not ont_terms:\n            continue\n\n        # Basic counts\n        n_terms = len(ont_terms)\n\n        # Compute depths\n        depths = {}\n        for term in ont_terms:\n            depths[term] = compute_depth(term, parents)\n\n        depth_values = list(depths.values())\n\n        # Count edges\n        n_edges = sum(len(parents.get(t, set())) for t in ont_terms)\n\n        # Branching factor (children per non-leaf node)\n        non_leaf_branching = []\n        for term in ont_terms:\n            if term in children and len(children[term]) > 0:\n                non_leaf_branching.append(len(children[term]))\n\n        # Roots and leaves\n        roots = [t for t in ont_terms if t not in parents or len(parents[t]) == 0]\n        leaves = [t for t in ont_terms if t not in children or len(children[t]) == 0]\n\n        stats[ont] = {\n            'n_terms': n_terms,\n            'n_edges': n_edges,\n            'n_roots': len(roots),\n            'n_leaves': len(leaves),\n            'max_depth': max(depth_values) if depth_values else 0,\n            'mean_depth': np.mean(depth_values) if depth_values else 0,\n            'median_depth': np.median(depth_values) if depth_values else 0,\n            'mean_branching': np.mean(non_leaf_branching) if non_leaf_branching else 0,\n            'max_branching': max(non_leaf_branching) if non_leaf_branching else 0,\n        }\n\n    return stats\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-14T10:35:19.372055Z","iopub.execute_input":"2025-11-14T10:35:19.372358Z","iopub.status.idle":"2025-11-14T10:35:19.386114Z","shell.execute_reply.started":"2025-11-14T10:35:19.372333Z","shell.execute_reply":"2025-11-14T10:35:19.385155Z"}},"outputs":[],"execution_count":null},{"id":"6048b5207e6ec9b5","cell_type":"code","source":"# =============================================================================\n# 3. LOAD ALL DATA\n# =============================================================================\nprint(\"\\n\" + \"=\" * 80)\nprint(\"LOADING CAFA6 DATASET\")\nprint(\"=\" * 80)\n\nprint(\"\\n[1/5] Loading GO ontology...\")\nparents, children, term_to_namespace, term_to_name, obsolete = parse_obo_detailed(\n    f'{TRAIN_DIR}/go-basic.obo'\n)\n\nprint(\"\\n[2/5] Loading train_terms.tsv...\")\ntrain_terms = pd.read_csv(f'{TRAIN_DIR}/train_terms.tsv', sep='\\t')\nprint(f\"✓ Loaded {len(train_terms):,} annotations\")\n\nprint(\"\\n[3/5] Loading train_taxonomy.tsv...\")\ntrain_taxonomy = pd.read_csv(f'{TRAIN_DIR}/train_taxonomy.tsv', sep='\\t', header=None)\ntrain_taxonomy.columns = ['EntryID', 'TaxonomyID']\nprint(f\"✓ Loaded {len(train_taxonomy):,} protein-taxonomy mappings\")\n\nprint(\"\\n[4/5] Loading IA.tsv...\")\nia_df = pd.read_csv(f'{BASE_PATH}/IA.tsv', sep='\\t', header=None)\nia_df.columns = ['GO_term', 'IA_weight']\nprint(f\"✓ Loaded {len(ia_df):,} IA weights\")\n\nprint(\"\\n[5/5] Loading train_sequences.fasta...\")\nn_seqs = 0\nseq_lengths = []\nfor record in tqdm(SeqIO.parse(f'{TRAIN_DIR}/train_sequences.fasta', 'fasta'), desc=\"Reading sequences\"):\n    n_seqs += 1\n    seq_lengths.append(len(record.seq))\n    if n_seqs >= 50000:  # Sample for speed\n        break\nprint(f\"✓ Sampled {n_seqs:,} sequences (mean length: {np.mean(seq_lengths):.0f}aa)\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-14T10:35:25.725407Z","iopub.execute_input":"2025-11-14T10:35:25.725869Z","iopub.status.idle":"2025-11-14T10:35:27.503661Z","shell.execute_reply.started":"2025-11-14T10:35:25.725840Z","shell.execute_reply":"2025-11-14T10:35:27.502568Z"}},"outputs":[],"execution_count":null},{"id":"93d47664e2dd2da3","cell_type":"code","source":"# =============================================================================\n# 4. COMPUTE DATASET STATISTICS\n# =============================================================================\nprint(\"\\n\" + \"=\" * 80)\nprint(\"COMPUTING STATISTICS\")\nprint(\"=\" * 80)\n\n# Basic counts\nn_proteins = train_terms['EntryID'].nunique()\nn_annotations = len(train_terms)\nn_unique_go_terms = train_terms['term'].nunique()\nn_species = train_taxonomy['TaxonomyID'].nunique()\n\n# Annotations per ontology\naspect_map = {'C': 'CC', 'F': 'MF', 'P': 'BP'}\ntrain_terms['Ontology'] = train_terms['aspect'].map(aspect_map)\nannotations_per_ont = train_terms['Ontology'].value_counts()\n\n# Annotations per protein\nannotations_per_protein = train_terms.groupby('EntryID').size()\n\n# Proteins per GO term\nproteins_per_go = train_terms.groupby('term').size()\n\n# DAG structure analysis\nprint(\"\\n[DAG Analysis] Computing depths and structure...\")\ndag_stats = analyze_dag_structure(parents, children, term_to_namespace)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-14T10:35:28.508077Z","iopub.execute_input":"2025-11-14T10:35:28.508409Z","iopub.status.idle":"2025-11-14T10:35:29.557354Z","shell.execute_reply.started":"2025-11-14T10:35:28.508385Z","shell.execute_reply":"2025-11-14T10:35:29.556303Z"}},"outputs":[],"execution_count":null},{"id":"e75125dcc6b9084c","cell_type":"code","source":"# =============================================================================\n# 5. CREATE SUMMARY TABLE\n# =============================================================================\nprint(\"\\n\" + \"=\" * 80)\nprint(\"DATASET STATISTICS SUMMARY\")\nprint(\"=\" * 80)\n\nsummary = {\n    'Metric': [],\n    'Value': [],\n    'Notes': []\n}\n\n# Overall statistics\nsummary['Metric'].extend([\n    'Total proteins',\n    'Total annotations',\n    'Unique GO terms (used)',\n    'Species (organisms)',\n    'Mean annotations/protein',\n    'Median annotations/protein',\n    'Mean proteins/GO term',\n    'Median proteins/GO term',\n])\n\nsummary['Value'].extend([\n    f\"{n_proteins:,}\",\n    f\"{n_annotations:,}\",\n    f\"{n_unique_go_terms:,}\",\n    f\"{n_species:,}\",\n    f\"{annotations_per_protein.mean():.1f}\",\n    f\"{annotations_per_protein.median():.0f}\",\n    f\"{proteins_per_go.mean():.1f}\",\n    f\"{proteins_per_go.median():.0f}\",\n])\n\nsummary['Notes'].extend([\n    'Unique protein identifiers',\n    'Total protein-GO associations',\n    'GO terms present in training',\n    'Unique taxonomy IDs',\n    'Sparsity indicator',\n    'Typical annotation density',\n    'Label frequency distribution',\n    'Most terms are rare',\n])\n\n# Per-ontology statistics\nfor ont in ['MF', 'BP', 'CC']:\n    ont_data = train_terms[train_terms['Ontology'] == ont]\n    ont_proteins = ont_data['EntryID'].nunique()\n    ont_annotations = len(ont_data)\n    ont_terms = ont_data['term'].nunique()\n\n    summary['Metric'].extend([\n        f'{ont} proteins',\n        f'{ont} annotations',\n        f'{ont} unique terms',\n    ])\n\n    summary['Value'].extend([\n        f\"{ont_proteins:,}\",\n        f\"{ont_annotations:,}\",\n        f\"{ont_terms:,}\",\n    ])\n\n    summary['Notes'].extend([\n        f'{ont_annotations / n_annotations * 100:.1f}% of total',\n        f'Mean {ont_annotations / ont_proteins:.1f} per protein',\n        f'{ont_terms / n_unique_go_terms * 100:.1f}% of used terms',\n    ])\n\n# DAG structure statistics\nsummary['Metric'].append('--- DAG STRUCTURE ---')\nsummary['Value'].append('')\nsummary['Notes'].append('')\n\nfor ont in ['MF', 'BP', 'CC']:\n    if ont not in dag_stats:\n        continue\n\n    stats = dag_stats[ont]\n\n    summary['Metric'].extend([\n        f'{ont} total terms (ontology)',\n        f'{ont} max depth',\n        f'{ont} mean depth',\n        f'{ont} edges (is_a + part_of)',\n        f'{ont} root nodes',\n        f'{ont} leaf nodes',\n        f'{ont} mean branching factor',\n    ])\n\n    summary['Value'].extend([\n        f\"{stats['n_terms']:,}\",\n        f\"{stats['max_depth']}\",\n        f\"{stats['mean_depth']:.2f}\",\n        f\"{stats['n_edges']:,}\",\n        f\"{stats['n_roots']}\",\n        f\"{stats['n_leaves']:,}\",\n        f\"{stats['mean_branching']:.2f}\",\n    ])\n\n    summary['Notes'].extend([\n        'Full ontology size',\n        'Longest path to root',\n        'Average hierarchy level',\n        'Parent-child relationships',\n        'Most general terms',\n        'Most specific terms',\n        'Avg children per non-leaf',\n    ])\n\nsummary_df = pd.DataFrame(summary)\nprint(summary_df.to_string(index=False))\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-14T10:37:05.660440Z","iopub.execute_input":"2025-11-14T10:37:05.660822Z","iopub.status.idle":"2025-11-14T10:37:05.938969Z","shell.execute_reply.started":"2025-11-14T10:37:05.660798Z","shell.execute_reply":"2025-11-14T10:37:05.937956Z"}},"outputs":[],"execution_count":null},{"id":"ab2f62eb8390583d","cell_type":"code","source":"# =============================================================================\n# 6. PRELIMINARY SUMMARY (TRAIN ONLY)\n# =============================================================================\n# This summary will be used later for the comprehensive LaTeX table\nprint(\"\\n[Preliminary] Train set statistics computed\")\n","metadata":{},"outputs":[],"execution_count":null},{"id":"db3961c89588ea13","cell_type":"code","source":"# =============================================================================\n# 7. VISUALIZATIONS\n# =============================================================================\nfig, axes = plt.subplots(2, 3, figsize=(18, 12))\n\n# 1. Annotations per ontology\naxes[0, 0].bar(['MF', 'BP', 'CC'],\n               [annotations_per_ont.get('MF', 0),\n                annotations_per_ont.get('BP', 0),\n                annotations_per_ont.get('CC', 0)],\n               color=['#FF6B6B', '#4ECDC4', '#45B7D1'], edgecolor='black')\naxes[0, 0].set_title('Annotations by Ontology', fontweight='bold')\naxes[0, 0].set_ylabel('Count')\n\n# 2. DAG depth distribution\nfor i, ont in enumerate(['MF', 'BP', 'CC']):\n    if ont in dag_stats:\n        axes[0, 1].bar(i, dag_stats[ont]['mean_depth'],\n                       color=['#FF6B6B', '#4ECDC4', '#45B7D1'][i],\n                       edgecolor='black', alpha=0.7)\n        axes[0, 1].text(i, dag_stats[ont]['mean_depth'],\n                        f\"{dag_stats[ont]['mean_depth']:.1f}\",\n                        ha='center', va='bottom', fontweight='bold')\naxes[0, 1].set_xticks([0, 1, 2])\naxes[0, 1].set_xticklabels(['MF', 'BP', 'CC'])\naxes[0, 1].set_title('Mean DAG Depth by Ontology', fontweight='bold')\naxes[0, 1].set_ylabel('Mean Depth')\n\n# 3. Annotations per protein (histogram)\naxes[0, 2].hist(annotations_per_protein, bins=50, color='#95E1D3', edgecolor='black', alpha=0.7)\naxes[0, 2].axvline(annotations_per_protein.mean(), color='red', linestyle='--', linewidth=2)\naxes[0, 2].set_title('Annotations per Protein', fontweight='bold')\naxes[0, 2].set_xlabel('# Annotations')\naxes[0, 2].set_ylabel('# Proteins')\naxes[0, 2].set_yscale('log')\n\n# 4. Proteins per GO term (log scale)\naxes[1, 0].hist(proteins_per_go, bins=50, color='#FFEAA7', edgecolor='black', alpha=0.7)\naxes[1, 0].set_title('Proteins per GO Term', fontweight='bold')\naxes[1, 0].set_xlabel('# Proteins')\naxes[1, 0].set_ylabel('# GO Terms')\naxes[1, 0].set_yscale('log')\n\n# 5. Top 10 species\ntop_species = train_taxonomy['TaxonomyID'].value_counts().head(10)\naxes[1, 1].barh(range(len(top_species)), top_species.values, color='#74B9FF', edgecolor='black')\naxes[1, 1].set_yticks(range(len(top_species)))\naxes[1, 1].set_yticklabels([f\"TaxID {tid}\" for tid in top_species.index], fontsize=8)\naxes[1, 1].set_title('Top 10 Species', fontweight='bold')\naxes[1, 1].set_xlabel('# Proteins')\naxes[1, 1].invert_yaxis()\n\n# 6. DAG complexity comparison\nmetrics = ['n_terms', 'n_edges', 'n_leaves']\nx = np.arange(len(metrics))\nwidth = 0.25\n\nfor i, ont in enumerate(['MF', 'BP', 'CC']):\n    if ont in dag_stats:\n        values = [dag_stats[ont][m] for m in metrics]\n        # Normalize for visualization\n        values_norm = [v / max(dag_stats[o][m] for o in ['MF', 'BP', 'CC'] if o in dag_stats)\n                       for m, v in zip(metrics, values)]\n        axes[1, 2].bar(x + i * width, values_norm, width,\n                       label=ont, color=['#FF6B6B', '#4ECDC4', '#45B7D1'][i],\n                       edgecolor='black', alpha=0.7)\n\naxes[1, 2].set_xticks(x + width)\naxes[1, 2].set_xticklabels(['Terms', 'Edges', 'Leaves'])\naxes[1, 2].set_title('DAG Complexity (Normalized)', fontweight='bold')\naxes[1, 2].set_ylabel('Normalized Value')\naxes[1, 2].legend()\n\nplt.tight_layout()\nplt.savefig('cafa6_dataset_analysis.png', dpi=300, bbox_inches='tight')\nplt.show()\n\nprint(\"\\n✅ Analysis complete! Visualizations saved to cafa6_dataset_analysis.png\")\n","metadata":{},"outputs":[],"execution_count":null},{"id":"59d392b231375599","cell_type":"code","source":"# =============================================================================\n# 8. TEST SET STATISTICS\n# =============================================================================\nprint(\"\\n\" + \"=\" * 80)\nprint(\"LOADING TEST SET DATA\")\nprint(\"=\" * 80)\n\nTEST_DIR = f'{BASE_PATH}/Test'\n\nprint(\"\\n[1/2] Loading testsuperset.fasta...\")\ntest_seq_lengths = []\ntest_protein_ids = []\nfor record in tqdm(SeqIO.parse(f'{TEST_DIR}/testsuperset.fasta', 'fasta'), desc=\"Reading test sequences\"):\n    test_protein_ids.append(record.id)\n    test_seq_lengths.append(len(record.seq))\nn_test_proteins = len(test_protein_ids)\nprint(f\"✓ Loaded {n_test_proteins:,} test sequences\")\nprint(f\"  Mean length: {np.mean(test_seq_lengths):.0f} aa\")\nprint(f\"  Median length: {np.median(test_seq_lengths):.0f} aa\")\nprint(f\"  Min length: {min(test_seq_lengths):.0f} aa\")\nprint(f\"  Max length: {max(test_seq_lengths):.0f} aa\")\n\nprint(\"\\n[2/2] Loading testsuperset-taxon-list.tsv...\")\ntest_taxonomy = pd.read_csv(f'{TEST_DIR}/testsuperset-taxon-list.tsv', sep='\\t', header=None)\ntest_taxonomy.columns = ['EntryID', 'TaxonomyID']\nprint(f\"✓ Loaded {len(test_taxonomy):,} protein-taxonomy mappings\")\nn_test_species = test_taxonomy['TaxonomyID'].nunique()\nprint(f\"  Unique species: {n_test_species:,}\")\n","metadata":{},"outputs":[],"execution_count":null},{"id":"589014e9bdd05e31","cell_type":"code","source":"# =============================================================================\n# 9. TEST SET STATISTICS SUMMARY\n# =============================================================================\nprint(\"\\n\" + \"=\" * 80)\nprint(\"TEST SET STATISTICS SUMMARY\")\nprint(\"=\" * 80)\n\ntest_summary = {\n    'Metric': [],\n    'Value': [],\n    'Notes': []\n}\n\ntest_summary['Metric'].extend([\n    'Total test proteins',\n    'Unique species',\n    'Mean sequence length',\n    'Median sequence length',\n    'Min sequence length',\n    'Max sequence length',\n    'Species overlap with train',\n])\n\n# Calculate species overlap\ntrain_species = set(train_taxonomy['TaxonomyID'].unique())\ntest_species = set(test_taxonomy['TaxonomyID'].unique())\nspecies_overlap = len(train_species & test_species)\nspecies_overlap_pct = species_overlap / n_test_species * 100 if n_test_species > 0 else 0\n\ntest_summary['Value'].extend([\n    f\"{n_test_proteins:,}\",\n    f\"{n_test_species:,}\",\n    f\"{np.mean(test_seq_lengths):.0f} aa\",\n    f\"{np.median(test_seq_lengths):.0f} aa\",\n    f\"{min(test_seq_lengths):.0f} aa\",\n    f\"{max(test_seq_lengths):.0f} aa\",\n    f\"{species_overlap:,} ({species_overlap_pct:.1f}%)\",\n])\n\ntest_summary['Notes'].extend([\n    'Proteins to predict',\n    'Unique taxonomy IDs',\n    'Average protein length',\n    'Typical protein length',\n    'Shortest protein',\n    'Longest protein',\n    'Species also in training set',\n])\n\ntest_summary_df = pd.DataFrame(test_summary)\nprint(test_summary_df.to_string(index=False))\n\n# Top species in test set\nprint(\"\\n\" + \"=\" * 80)\nprint(\"TOP 10 SPECIES IN TEST SET\")\nprint(\"=\" * 80)\ntop_test_species = test_taxonomy['TaxonomyID'].value_counts().head(10)\nfor i, (taxid, count) in enumerate(top_test_species.items(), 1):\n    in_train = \"✓\" if taxid in train_species else \"✗\"\n    print(f\"{i:2d}. TaxID {taxid:>10} : {count:>6,} proteins [{in_train} in train]\")\n","metadata":{},"outputs":[],"execution_count":null},{"id":"845d940ff2aa354d","cell_type":"code","source":"# =============================================================================\n# 9b. COMPUTE COMBINED DATASET STATISTICS (TRAIN + TEST)\n# =============================================================================\nprint(\"\\n\" + \"=\" * 80)\nprint(\"COMPUTING COMBINED DATASET STATISTICS\")\nprint(\"=\" * 80)\n\n# Combine protein IDs from train and test (removing duplicates)\ntrain_protein_ids = set(train_terms['EntryID'].unique())\ntest_protein_ids_set = set(test_protein_ids)\nall_protein_ids = train_protein_ids | test_protein_ids_set\nn_combined_proteins = len(all_protein_ids)\nn_overlap_proteins = len(train_protein_ids & test_protein_ids_set)\n\n# Combine species (removing duplicates)\nall_species = train_species | test_species\nn_combined_species = len(all_species)\n\nprint(f\"✓ Combined unique proteins: {n_combined_proteins:,}\")\nprint(f\"  Train proteins: {len(train_protein_ids):,}\")\nprint(f\"  Test proteins: {len(test_protein_ids_set):,}\")\nprint(f\"  Overlapping proteins: {n_overlap_proteins:,}\")\nprint(f\"✓ Combined unique species: {n_combined_species:,}\")\n","metadata":{},"outputs":[],"execution_count":null},{"id":"25541445545caff1","cell_type":"code","source":"# =============================================================================\n# 9c. COMPREHENSIVE LATEX TABLE (TRAIN + TEST + COMBINED)\n# =============================================================================\nprint(\"\\n\" + \"=\" * 80)\nprint(\"GENERATING COMPREHENSIVE LATEX TABLE\")\nprint(\"=\" * 80)\n\nlatex_table = r\"\"\"\n\\begin{table}[h]\n\\centering\n\\caption{CAFA6 Dataset Statistics: Train, Test, and Combined}\n\\label{tab:cafa6_stats}\n\\begin{tabular}{lrrr}\n\\hline\n\\textbf{Metric} & \\textbf{Train} & \\textbf{Test} & \\textbf{Combined} \\\\\n\\hline\n\"\"\"\n\n# Dataset Overview\nlatex_table += r\"\\multicolumn{4}{c}{\\textit{Dataset Overview}} \\\\\" + \"\\n\"\nlatex_table += r\"\\hline\" + \"\\n\"\nlatex_table += f\"Total proteins & {n_proteins:,} & {n_test_proteins:,} & {n_combined_proteins:,} \\\\\\\\\\n\"\nlatex_table += f\"Unique species & {n_species:,} & {n_test_species:,} & {n_combined_species:,} \\\\\\\\\\n\"\nlatex_table += f\"Protein overlap & — & — & {n_overlap_proteins:,} \\\\\\\\\\n\"\nlatex_table += f\"Total annotations & {n_annotations:,} & — & — \\\\\\\\\\n\"\nlatex_table += f\"Unique GO terms & {n_unique_go_terms:,} & — & — \\\\\\\\\\n\"\nlatex_table += r\"\\hline\" + \"\\n\"\n\n# Sequence Statistics\nlatex_table += r\"\\multicolumn{4}{c}{\\textit{Sequence Statistics}} \\\\\" + \"\\n\"\nlatex_table += r\"\\hline\" + \"\\n\"\nlatex_table += f\"Mean sequence length (aa) & {np.mean(seq_lengths):.0f} & {np.mean(test_seq_lengths):.0f} & {np.mean(seq_lengths + test_seq_lengths):.0f} \\\\\\\\\\n\"\nlatex_table += f\"Median sequence length (aa) & {np.median(seq_lengths):.0f} & {np.median(test_seq_lengths):.0f} & {np.median(seq_lengths + test_seq_lengths):.0f} \\\\\\\\\\n\"\nlatex_table += f\"Min sequence length (aa) & {min(seq_lengths):.0f} & {min(test_seq_lengths):.0f} & {min(seq_lengths + test_seq_lengths):.0f} \\\\\\\\\\n\"\nlatex_table += f\"Max sequence length (aa) & {max(seq_lengths):.0f} & {max(test_seq_lengths):.0f} & {max(seq_lengths + test_seq_lengths):.0f} \\\\\\\\\\n\"\nlatex_table += r\"\\hline\" + \"\\n\"\n\n# Annotation Statistics (Train only)\nlatex_table += r\"\\multicolumn{4}{c}{\\textit{Annotation Statistics (Train Only)}} \\\\\" + \"\\n\"\nlatex_table += r\"\\hline\" + \"\\n\"\nlatex_table += f\"Mean annotations/protein & {annotations_per_protein.mean():.1f} & — & — \\\\\\\\\\n\"\nlatex_table += f\"Median annotations/protein & {annotations_per_protein.median():.0f} & — & — \\\\\\\\\\n\"\nlatex_table += f\"Mean proteins/GO term & {proteins_per_go.mean():.1f} & — & — \\\\\\\\\\n\"\nlatex_table += f\"Median proteins/GO term & {proteins_per_go.median():.0f} & — & — \\\\\\\\\\n\"\nlatex_table += r\"\\hline\" + \"\\n\"\n\n# Per-Ontology Statistics\nlatex_table += r\"\\multicolumn{4}{c}{\\textit{Per-Ontology Statistics (Train Only)}} \\\\\" + \"\\n\"\nlatex_table += r\"\\hline\" + \"\\n\"\nfor ont in ['MF', 'BP', 'CC']:\n    ont_data = train_terms[train_terms['Ontology'] == ont]\n    ont_proteins = ont_data['EntryID'].nunique()\n    ont_annotations = len(ont_data)\n    ont_terms = ont_data['term'].nunique()\n\n    ont_name = {'MF': 'Molecular Function', 'BP': 'Biological Process', 'CC': 'Cellular Component'}[ont]\n    latex_table += f\"{ont_name} proteins & {ont_proteins:,} & — & — \\\\\\\\\\n\"\n    latex_table += f\"{ont_name} annotations & {ont_annotations:,} & — & — \\\\\\\\\\n\"\n    latex_table += f\"{ont_name} unique terms & {ont_terms:,} & — & — \\\\\\\\\\n\"\n\nlatex_table += r\"\\hline\" + \"\\n\"\n\n# DAG Structure\nlatex_table += r\"\\multicolumn{4}{c}{\\textit{GO Ontology DAG Structure}} \\\\\" + \"\\n\"\nlatex_table += r\"\\hline\" + \"\\n\"\nfor ont in ['MF', 'BP', 'CC']:\n    if ont not in dag_stats:\n        continue\n\n    stats = dag_stats[ont]\n    ont_name = {'MF': 'Molecular Function', 'BP': 'Biological Process', 'CC': 'Cellular Component'}[ont]\n\n    latex_table += f\"{ont_name} total terms & {stats['n_terms']:,} & — & — \\\\\\\\\\n\"\n    latex_table += f\"{ont_name} max depth & {stats['max_depth']} & — & — \\\\\\\\\\n\"\n    latex_table += f\"{ont_name} mean depth & {stats['mean_depth']:.2f} & — & — \\\\\\\\\\n\"\n    latex_table += f\"{ont_name} edges & {stats['n_edges']:,} & — & — \\\\\\\\\\n\"\n    latex_table += f\"{ont_name} leaf nodes & {stats['n_leaves']:,} & — & — \\\\\\\\\\n\"\n\nlatex_table += r\"\"\"\\hline\n\\end{tabular}\n\\end{table}\n\"\"\"\n\nprint(latex_table)\n\n# Also save to file\nwith open('cafa6_statistics_table.tex', 'w') as f:\n    f.write(latex_table)\nprint(\"\\n✅ LaTeX table saved to cafa6_statistics_table.tex\")\n","metadata":{},"outputs":[],"execution_count":null},{"id":"e6c959c91a8b713e","cell_type":"code","source":"# =============================================================================\n# 10. TRAIN vs TEST COMPARISON VISUALIZATIONS\n# =============================================================================\nfig, axes = plt.subplots(2, 2, figsize=(14, 10))\n\n# 1. Dataset size comparison\naxes[0, 0].bar(['Train', 'Test'], [n_proteins, n_test_proteins],\n               color=['#4ECDC4', '#FF6B6B'], edgecolor='black', alpha=0.7)\naxes[0, 0].set_title('Dataset Size Comparison', fontweight='bold')\naxes[0, 0].set_ylabel('# Proteins')\nfor i, v in enumerate([n_proteins, n_test_proteins]):\n    axes[0, 0].text(i, v, f'{v:,}', ha='center', va='bottom', fontweight='bold')\n\n# 2. Species count comparison\naxes[0, 1].bar(['Train', 'Test'], [n_species, n_test_species],\n               color=['#4ECDC4', '#FF6B6B'], edgecolor='black', alpha=0.7)\naxes[0, 1].set_title('Species Count Comparison', fontweight='bold')\naxes[0, 1].set_ylabel('# Unique Species')\nfor i, v in enumerate([n_species, n_test_species]):\n    axes[0, 1].text(i, v, f'{v:,}', ha='center', va='bottom', fontweight='bold')\n\n# 3. Sequence length distribution comparison\naxes[1, 0].hist(seq_lengths, bins=50, alpha=0.6, label='Train (sampled)',\n                color='#4ECDC4', edgecolor='black')\naxes[1, 0].hist(test_seq_lengths, bins=50, alpha=0.6, label='Test',\n                color='#FF6B6B', edgecolor='black')\naxes[1, 0].axvline(np.mean(seq_lengths), color='#4ECDC4', linestyle='--', linewidth=2)\naxes[1, 0].axvline(np.mean(test_seq_lengths), color='#FF6B6B', linestyle='--', linewidth=2)\naxes[1, 0].set_title('Sequence Length Distribution', fontweight='bold')\naxes[1, 0].set_xlabel('Sequence Length (aa)')\naxes[1, 0].set_ylabel('Frequency')\naxes[1, 0].legend()\naxes[1, 0].set_xlim(0, 2000)\n\n# 4. Species overlap Venn diagram representation\noverlap_data = {\n    'Train only': len(train_species - test_species),\n    'Both': species_overlap,\n    'Test only': len(test_species - train_species)\n}\ncolors_venn = ['#4ECDC4', '#95E1D3', '#FF6B6B']\naxes[1, 1].bar(range(3), overlap_data.values(), color=colors_venn, edgecolor='black', alpha=0.7)\naxes[1, 1].set_xticks(range(3))\naxes[1, 1].set_xticklabels(overlap_data.keys())\naxes[1, 1].set_title('Species Distribution', fontweight='bold')\naxes[1, 1].set_ylabel('# Species')\nfor i, (k, v) in enumerate(overlap_data.items()):\n    axes[1, 1].text(i, v, f'{v:,}', ha='center', va='bottom', fontweight='bold')\n\nplt.tight_layout()\nplt.savefig('cafa6_train_test_comparison.png', dpi=300, bbox_inches='tight')\nplt.show()\n\nprint(\"\\n✅ Test set analysis complete! Comparison saved to cafa6_train_test_comparison.png\")\n","metadata":{},"outputs":[],"execution_count":null},{"id":"797291aeac1cb2ca","cell_type":"code","source":"","metadata":{},"outputs":[],"execution_count":null}]}