{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"},"kaggle":{"accelerator":"nvidiaTeslaT4","dataSources":[{"sourceId":116062,"databundleVersionId":14084779,"sourceType":"competition"},{"sourceId":5499219,"sourceType":"datasetVersion","datasetId":3167603},{"sourceId":5549164,"sourceType":"datasetVersion","datasetId":3197305},{"sourceId":6247561,"sourceType":"datasetVersion","datasetId":3590060},{"sourceId":5792099,"sourceType":"datasetVersion","datasetId":3327296},{"sourceId":5607816,"sourceType":"datasetVersion","datasetId":3225525}],"dockerImageVersionId":31153,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================================\n# CAFA 6 PROTEIN FUNCTION PREDICTION - ENHANCED VERSION (FIXED)\n# ============================================================================\n\n# ⚙️ CONFIGURATION - ADJUST THIS FOR SPEED VS ACCURACY\nSAMPLE_PERCENT = 100  # Use 100% of data\nQUICK_MODE = True   # Enable full feature computation\n\n# Package Installation\nimport subprocess\nimport sys\n\ndef install(package):\n    subprocess.check_call([sys.executable, \"-m\", \"pip\", \"install\", \"-q\", package])\n\nprint(\"Installing required packages...\")\ntry:\n    import obonet\nexcept:\n    install('obonet')\n    import obonet\n\ntry:\n    from Bio import SeqIO\nexcept:\n    install('biopython')\n    from Bio import SeqIO\n\n# Core Imports\nimport pandas as pd\nimport numpy as np\nfrom pathlib import Path\nfrom collections import defaultdict, Counter\nimport warnings\nwarnings.filterwarnings('ignore')\n\nimport networkx as nx\nfrom sklearn.preprocessing import MultiLabelBinarizer\nfrom sklearn.multiclass import OneVsRestClassifier\nfrom sklearn.linear_model import LogisticRegression\n\n# Visualization imports\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom matplotlib.patches import Rectangle\nimport matplotlib.patches as mpatches\nplt.style.use('seaborn-v0_8-whitegrid')\nsns.set_palette(\"husl\")\n\nprint(\"=\"*80)\nprint(\"CAFA 6 PROTEIN FUNCTION PREDICTION - ENHANCED STARTER\")\nprint(f\"📊 SAMPLE MODE: {SAMPLE_PERCENT}% of data\")\nprint(f\"⚡ QUICK MODE: {'ON' if QUICK_MODE else 'OFF'}\")\nprint(\"=\"*80)\n\n# ============================================================================\n# 1. DEFINE PATHS\n# ============================================================================\nBASE = Path('/kaggle/input/cafa-6-protein-function-prediction')\nTRAIN_DIR = BASE / 'Train'\nTEST_DIR = BASE / 'Test'\n\n# ============================================================================\n# 2. LOAD GO ONTOLOGY (WITH HIERARCHY ANALYSIS)\n# ============================================================================\nprint(\"\\n[1/9] Loading GO ontology...\")\ngo_graph = obonet.read_obo(TRAIN_DIR / 'go-basic.obo')\nprint(f\"   ✓ Loaded {len(go_graph)} GO terms\")\n\n# Map terms to ontologies\nterm_to_ont = {}\nterm_names = {}\nfor term_id in go_graph.nodes():\n    if 'namespace' in go_graph.nodes[term_id]:\n        ns = go_graph.nodes[term_id]['namespace']\n        if ns == 'biological_process':\n            term_to_ont[term_id] = 'BPO'\n        elif ns == 'cellular_component':\n            term_to_ont[term_id] = 'CCO'\n        elif ns == 'molecular_function':\n            term_to_ont[term_id] = 'MFO'\n    if 'name' in go_graph.nodes[term_id]:\n        term_names[term_id] = go_graph.nodes[term_id]['name']\n\nont_counts = pd.Series(term_to_ont).value_counts()\nprint(f\"   ✓ Ontology breakdown: MF={ont_counts.get('MFO',0)}, BP={ont_counts.get('BPO',0)}, CC={ont_counts.get('CCO',0)}\")\n\n# Analyze GO hierarchy depth (sample for speed)\ndef get_term_depth(graph, term_id):\n    \"\"\"Calculate depth of term in GO hierarchy\"\"\"\n    try:\n        paths = []\n        for root in ['GO:0008150', 'GO:0005575', 'GO:0003674']:\n            if nx.has_path(graph, term_id, root):\n                paths.append(nx.shortest_path_length(graph, term_id, root))\n        return max(paths) if paths else 0\n    except:\n        return 0\n\nprint(\"   Computing GO hierarchy depths...\")\nsample_terms_for_depth = list(term_to_ont.keys())[:1000]\nterm_depths = {term: get_term_depth(go_graph, term) for term in sample_terms_for_depth}\n\n# Visualize ontology with enhanced graphics\nfig = plt.figure(figsize=(18, 10))\ngs = fig.add_gridspec(3, 3, hspace=0.3, wspace=0.3)\n\n# Main ontology distribution\nax1 = fig.add_subplot(gs[0, :2])\ncolors = ['#FF6B6B', '#4ECDC4', '#45B7D1']\nbars = ax1.bar(range(len(ont_counts)), ont_counts.values, color=colors, \n               edgecolor='black', linewidth=2, alpha=0.8)\nax1.set_xticks(range(len(ont_counts)))\nax1.set_xticklabels(['Molecular Function', 'Biological Process', 'Cellular Component'], \n                     rotation=0, fontsize=11, fontweight='bold')\nax1.set_title('GO Term Distribution by Ontology', fontsize=14, fontweight='bold', pad=20)\nax1.set_ylabel('Number of Terms', fontsize=12, fontweight='bold')\nax1.grid(axis='y', alpha=0.3)\nfor i, (v, bar) in enumerate(zip(ont_counts.values, bars)):\n    height = bar.get_height()\n    ax1.text(bar.get_x() + bar.get_width()/2., height,\n             f'{v:,}\\n({v/ont_counts.sum()*100:.1f}%)',\n             ha='center', va='bottom', fontweight='bold', fontsize=11)\n\n# Hierarchy depth distribution\nax2 = fig.add_subplot(gs[0, 2])\ndepth_values = list(term_depths.values())\nax2.hist(depth_values, bins=20, color='#A8E6CF', edgecolor='black', alpha=0.7)\nax2.set_title('GO Term Depth\\nDistribution', fontsize=11, fontweight='bold')\nax2.set_xlabel('Hierarchy Depth', fontsize=10)\nax2.set_ylabel('Count', fontsize=10)\nax2.axvline(np.mean(depth_values), color='red', linestyle='--', linewidth=2,\n            label=f'Mean: {np.mean(depth_values):.1f}')\nax2.legend(fontsize=9)\n\n# Network visualization (sample of GO graph)\nax3 = fig.add_subplot(gs[1:, :])\nsample_terms = list(term_to_ont.keys())[:50]\nsubgraph = go_graph.subgraph(sample_terms)\npos = nx.spring_layout(subgraph, k=0.5, iterations=50, seed=42)\nnode_colors = [colors[['MFO', 'BPO', 'CCO'].index(term_to_ont.get(node, 'MFO'))] \n               for node in subgraph.nodes()]\nnx.draw_networkx_nodes(subgraph, pos, node_color=node_colors, \n                       node_size=300, alpha=0.7, ax=ax3)\nnx.draw_networkx_edges(subgraph, pos, alpha=0.2, arrows=True, \n                       arrowsize=10, ax=ax3, edge_color='gray')\nax3.set_title('GO Ontology Network Structure (Sample of 50 terms)', \n              fontsize=13, fontweight='bold', pad=15)\nax3.axis('off')\n\n# Legend\nlegend_elements = [mpatches.Patch(facecolor=colors[0], label='Molecular Function'),\n                   mpatches.Patch(facecolor=colors[1], label='Biological Process'),\n                   mpatches.Patch(facecolor=colors[2], label='Cellular Component')]\nax3.legend(handles=legend_elements, loc='upper right', fontsize=10, framealpha=0.9)\n\nplt.suptitle('Gene Ontology Analysis', fontsize=16, fontweight='bold', y=0.995)\nplt.tight_layout()\nplt.show()\n\n# ============================================================================\n# 3. LOAD IA WEIGHTS (WITH ANALYSIS)\n# ============================================================================\nprint(\"\\n[2/9] Loading IA weights...\")\nia_df = pd.read_csv(BASE / 'IA.tsv', sep='\\t', header=None, names=['term', 'ia'])\n\nif SAMPLE_PERCENT < 100:\n    ia_df = ia_df.sample(frac=SAMPLE_PERCENT/100, random_state=42)\n\nia_dict = dict(zip(ia_df['term'], ia_df['ia']))\nprint(f\"   ✓ Loaded {len(ia_dict)} IA weights\")\n\n# Enhanced IA visualization\nfig, axes = plt.subplots(2, 3, figsize=(18, 10))\n\n# IA distribution by ontology\nia_by_ont = ia_df.copy()\nia_by_ont['ontology'] = ia_by_ont['term'].map(term_to_ont)\nia_by_ont = ia_by_ont.dropna()\n\naxes[0, 0].hist(ia_df['ia'], bins=50, color='#95E1D3', edgecolor='black', alpha=0.7)\naxes[0, 0].set_title('Overall IA Distribution', fontsize=12, fontweight='bold')\naxes[0, 0].set_xlabel('IA Weight', fontsize=10)\naxes[0, 0].set_ylabel('Frequency', fontsize=10)\naxes[0, 0].axvline(ia_df['ia'].mean(), color='red', linestyle='--', linewidth=2, \n                   label=f'Mean: {ia_df[\"ia\"].mean():.2f}')\naxes[0, 0].legend()\n\n# Box plot by ontology\nont_data = [ia_by_ont[ia_by_ont['ontology']==ont]['ia'].values \n            for ont in ['MFO', 'BPO', 'CCO']]\nbp = axes[0, 1].boxplot(ont_data, labels=['MF', 'BP', 'CC'], patch_artist=True)\nfor patch, color in zip(bp['boxes'], colors):\n    patch.set_facecolor(color)\n    patch.set_alpha(0.7)\naxes[0, 1].set_title('IA Weights by Ontology', fontsize=12, fontweight='bold')\naxes[0, 1].set_ylabel('IA Weight', fontsize=10)\naxes[0, 1].grid(axis='y', alpha=0.3)\n\n# Violin plot\nparts = axes[0, 2].violinplot(ont_data, positions=[1, 2, 3], showmeans=True, showmedians=True)\nfor pc, color in zip(parts['bodies'], colors):\n    pc.set_facecolor(color)\n    pc.set_alpha(0.7)\naxes[0, 2].set_xticks([1, 2, 3])\naxes[0, 2].set_xticklabels(['MF', 'BP', 'CC'])\naxes[0, 2].set_title('IA Distribution Density', fontsize=12, fontweight='bold')\naxes[0, 2].set_ylabel('IA Weight', fontsize=10)\n\n# Cumulative distribution\nsorted_ia = np.sort(ia_df['ia'].values)\ncumsum = np.cumsum(sorted_ia) / np.sum(sorted_ia)\naxes[1, 0].plot(sorted_ia, cumsum, linewidth=2, color='#6C5CE7')\naxes[1, 0].set_title('Cumulative IA Distribution', fontsize=12, fontweight='bold')\naxes[1, 0].set_xlabel('IA Weight', fontsize=10)\naxes[1, 0].set_ylabel('Cumulative Proportion', fontsize=10)\naxes[1, 0].grid(alpha=0.3)\naxes[1, 0].axhline(0.5, color='red', linestyle='--', alpha=0.5, label='50%')\naxes[1, 0].legend()\n\n# Top terms by IA\ntop_ia = ia_df.nlargest(15, 'ia')\naxes[1, 1].barh(range(len(top_ia)), top_ia['ia'].values, color='#FF7675', edgecolor='black')\naxes[1, 1].set_yticks(range(len(top_ia)))\naxes[1, 1].set_yticklabels([f\"{t[:15]}...\" if len(t) > 15 else t \n                            for t in top_ia['term'].values], fontsize=8)\naxes[1, 1].set_title('Top 15 Terms by IA Weight', fontsize=12, fontweight='bold')\naxes[1, 1].set_xlabel('IA Weight', fontsize=10)\naxes[1, 1].invert_yaxis()\n\n# Statistics summary\naxes[1, 2].axis('off')\nia_stats = f\"\"\"\nIA WEIGHT STATISTICS\n\nTotal terms: {len(ia_df):,}\n\nOverall:\n  • Mean: {ia_df['ia'].mean():.3f}\n  • Median: {ia_df['ia'].median():.3f}\n  • Std Dev: {ia_df['ia'].std():.3f}\n  • Range: [{ia_df['ia'].min():.3f}, {ia_df['ia'].max():.3f}]\n\nBy Ontology (Mean ± Std):\n  • MF: {ia_by_ont[ia_by_ont['ontology']=='MFO']['ia'].mean():.3f} ± {ia_by_ont[ia_by_ont['ontology']=='MFO']['ia'].std():.3f}\n  • BP: {ia_by_ont[ia_by_ont['ontology']=='BPO']['ia'].mean():.3f} ± {ia_by_ont[ia_by_ont['ontology']=='BPO']['ia'].std():.3f}\n  • CC: {ia_by_ont[ia_by_ont['ontology']=='CCO']['ia'].mean():.3f} ± {ia_by_ont[ia_by_ont['ontology']=='CCO']['ia'].std():.3f}\n\"\"\"\naxes[1, 2].text(0.05, 0.5, ia_stats, fontsize=10, family='monospace',\n                verticalalignment='center')\n\nplt.suptitle('Information Accretion (IA) Analysis', fontsize=14, fontweight='bold')\nplt.tight_layout()\nplt.show()\n\n# ============================================================================\n# 4. LOAD TRAINING DATA (WITH COMPREHENSIVE ANALYSIS) - FIXED\n# ============================================================================\nprint(\"\\n[3/9] Loading training data...\")\n\ntrain_terms = pd.read_csv(TRAIN_DIR / 'train_terms.tsv', sep='\\t', \n                          names=['protein', 'term', 'ontology'])\ntrain_taxonomy = pd.read_csv(TRAIN_DIR / 'train_taxonomy.tsv', sep='\\t',\n                             names=['protein', 'taxon'])\n\nprint(f\"   ✓ Full dataset: {len(train_terms)} annotations, {train_terms['protein'].nunique()} proteins\")\n\n# SAMPLE proteins for faster iteration\nif SAMPLE_PERCENT < 100:\n    sample_proteins = train_terms['protein'].drop_duplicates().sample(\n        frac=SAMPLE_PERCENT/100, random_state=42\n    ).tolist()\n    train_terms = train_terms[train_terms['protein'].isin(sample_proteins)]\n    train_taxonomy = train_taxonomy[train_taxonomy['protein'].isin(sample_proteins)]\n    print(f\"   ✓ Sampled to {SAMPLE_PERCENT}%: {len(train_terms)} annotations, {len(sample_proteins)} proteins\")\n\n# Print ontology distribution\nprint(f\"\\n   Ontology distribution:\")\nprint(train_terms['ontology'].value_counts())\n\n# Comprehensive training data visualization - FIXED\nfig = plt.figure(figsize=(20, 12))\ngs = fig.add_gridspec(3, 4, hspace=0.35, wspace=0.35)\n\n# 1. Ontology distribution - FIXED to handle all possible ontology codes\nax1 = fig.add_subplot(gs[0, 0])\nont_dist = train_terms['ontology'].value_counts()\n\n# Map ontology codes (handle F, P, C or any other codes)\ncolors_ont_map = {'F': '#FF6B6B', 'P': '#4ECDC4', 'C': '#45B7D1'}\nont_names_map = {'F': 'MF', 'P': 'BP', 'C': 'CC'}\n\n# Get colors and names, with defaults for unknown codes\ncolors_list = [colors_ont_map.get(k, '#CCCCCC') for k in ont_dist.index]\nlabels_list = [ont_names_map.get(k, k) for k in ont_dist.index]\n\nbars = ax1.bar(range(len(ont_dist)), ont_dist.values, color=colors_list, \n               edgecolor='black', linewidth=1.5)\nax1.set_xticks(range(len(ont_dist)))\nax1.set_xticklabels(labels_list)\nax1.set_title('Annotations by Ontology', fontsize=11, fontweight='bold')\nax1.set_ylabel('Count', fontsize=9)\nfor i, (v, bar) in enumerate(zip(ont_dist.values, bars)):\n    ax1.text(bar.get_x() + bar.get_width()/2., v, f'{v:,}', \n             ha='center', va='bottom', fontweight='bold', fontsize=9)\n\n# 2. Top terms\nax2 = fig.add_subplot(gs[0, 1:3])\ntop_terms = train_terms['term'].value_counts().head(20)\nax2.barh(range(len(top_terms)), top_terms.values, color='#A8E6CF', edgecolor='black')\nax2.set_yticks(range(len(top_terms)))\nax2.set_yticklabels([f\"{term_names.get(t, t)[:30]}...\" if len(term_names.get(t, t)) > 30 \n                     else term_names.get(t, t) for t in top_terms.index], fontsize=8)\nax2.set_title('Top 20 Most Frequent GO Terms', fontsize=11, fontweight='bold')\nax2.set_xlabel('Count', fontsize=9)\nax2.invert_yaxis()\n\n# 3. Terms per protein\nax3 = fig.add_subplot(gs[0, 3])\nterms_per_protein = train_terms.groupby('protein').size()\nax3.hist(terms_per_protein, bins=50, color='#FFD93D', edgecolor='black', alpha=0.7)\nax3.set_title('Terms per Protein', fontsize=11, fontweight='bold')\nax3.set_xlabel('# Terms', fontsize=9)\nax3.set_ylabel('Frequency', fontsize=9)\nax3.axvline(terms_per_protein.mean(), color='red', linestyle='--', linewidth=2)\n\n# 4. Proteins per term\nax4 = fig.add_subplot(gs[1, 0])\nproteins_per_term = train_terms.groupby('term').size()\nax4.hist(proteins_per_term, bins=50, color='#FFEAA7', edgecolor='black', alpha=0.7, log=True)\nax4.set_title('Proteins per Term (log)', fontsize=11, fontweight='bold')\nax4.set_xlabel('# Proteins', fontsize=9)\nax4.set_ylabel('# Terms (log)', fontsize=9)\n\n# 5. Taxonomy distribution\nax5 = fig.add_subplot(gs[1, 1])\ntop_taxa = train_taxonomy['taxon'].value_counts().head(10)\nax5.bar(range(len(top_taxa)), top_taxa.values, color='#74B9FF', edgecolor='black')\nax5.set_xticks(range(len(top_taxa)))\nax5.set_xticklabels([str(t)[:8] for t in top_taxa.index], rotation=45, ha='right', fontsize=8)\nax5.set_title('Top 10 Species', fontsize=11, fontweight='bold')\nax5.set_ylabel('# Proteins', fontsize=9)\n\n# 6. Term co-occurrence heatmap\nax6 = fig.add_subplot(gs[1, 2:])\ntop_10_terms = train_terms['term'].value_counts().head(10).index\ncooc_matrix = np.zeros((10, 10))\nfor i, t1 in enumerate(top_10_terms):\n    for j, t2 in enumerate(top_10_terms):\n        if i != j:\n            proteins_t1 = set(train_terms[train_terms['term']==t1]['protein'])\n            proteins_t2 = set(train_terms[train_terms['term']==t2]['protein'])\n            cooc_matrix[i,j] = len(proteins_t1 & proteins_t2)\nim = ax6.imshow(cooc_matrix, cmap='YlOrRd', aspect='auto')\nax6.set_xticks(range(10))\nax6.set_yticks(range(10))\nax6.set_xticklabels([term_names.get(t, t)[:10] for t in top_10_terms], \n                     rotation=45, ha='right', fontsize=7)\nax6.set_yticklabels([term_names.get(t, t)[:10] for t in top_10_terms], fontsize=7)\nax6.set_title('Term Co-occurrence Matrix', fontsize=11, fontweight='bold')\nplt.colorbar(im, ax=ax6, label='# Shared Proteins')\n\n# 7. Annotation density\nax7 = fig.add_subplot(gs[2, :2])\nterm_freq_bins = pd.cut(proteins_per_term, bins=[0, 10, 50, 100, 500, 100000], \n                        labels=['<10', '10-50', '50-100', '100-500', '>500'])\nfreq_dist = term_freq_bins.value_counts().sort_index()\nax7.bar(range(len(freq_dist)), freq_dist.values, color='#E17055', edgecolor='black', alpha=0.7)\nax7.set_xticks(range(len(freq_dist)))\nax7.set_xticklabels(freq_dist.index, rotation=0)\nax7.set_title('GO Term Frequency Distribution', fontsize=11, fontweight='bold')\nax7.set_xlabel('# Proteins with Term', fontsize=9)\nax7.set_ylabel('# Terms', fontsize=9)\nfor i, v in enumerate(freq_dist.values):\n    ax7.text(i, v, str(v), ha='center', va='bottom', fontweight='bold')\n\n# 8. Summary statistics\nax8 = fig.add_subplot(gs[2, 2:])\nax8.axis('off')\nsummary_text = f\"\"\"\nTRAINING DATA COMPREHENSIVE SUMMARY\n\nDataset Size:\n  • Total Annotations: {len(train_terms):,}\n  • Unique Proteins: {train_terms['protein'].nunique():,}\n  • Unique GO Terms: {train_terms['term'].nunique():,}\n  • Species: {train_taxonomy['taxon'].nunique()}\n\nOntology Distribution:\n  • Molecular Function: {ont_dist.get('F', 0):,} ({ont_dist.get('F', 0)/len(train_terms)*100:.1f}%)\n  • Biological Process: {ont_dist.get('P', 0):,} ({ont_dist.get('P', 0)/len(train_terms)*100:.1f}%)\n  • Cellular Component: {ont_dist.get('C', 0):,} ({ont_dist.get('C', 0)/len(train_terms)*100:.1f}%)\n\nAnnotation Statistics:\n  • Mean terms/protein: {terms_per_protein.mean():.1f}\n  • Median terms/protein: {terms_per_protein.median():.0f}\n  • Max terms/protein: {terms_per_protein.max()}\n  • Mean proteins/term: {proteins_per_term.mean():.1f}\n  • Median proteins/term: {proteins_per_term.median():.0f}\n\"\"\"\nax8.text(0.05, 0.5, summary_text, fontsize=10, family='monospace',\n         verticalalignment='center', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.3))\n\nplt.suptitle('Training Data Comprehensive Analysis', fontsize=15, fontweight='bold')\nplt.tight_layout()\nplt.show()\n\n# Continue with the rest of the code (sequences, features, training, etc.)\nprint(\"\\n   Loading sequences (this may take a while for 100% of data)...\")\nprint(f\"   Expected proteins: {train_terms['protein'].nunique():,}\")\n\ntrain_seqs = {}\nloaded_count = 0\ntarget_proteins = set(train_terms['protein'].unique())\n\nfor rec in SeqIO.parse(TRAIN_DIR / 'train_sequences.fasta', 'fasta'):\n    pid = rec.id.split('|')[1] if '|' in rec.id else rec.id\n    if pid in target_proteins:\n        train_seqs[pid] = str(rec.seq)\n        loaded_count += 1\n        \n        # Progress indicator\n        if loaded_count % 10000 == 0:\n            print(f\"      Loaded {loaded_count:,} sequences...\")\n        \n    if loaded_count >= len(target_proteins):\n        break\n\nprint(f\"   ✓ Loaded {len(train_seqs):,} training sequences\")\n\n# Enhanced sequence analysis\nseq_lengths = [len(s) for s in train_seqs.values()]\nprint(f\"   ✓ Sequence length: mean={np.mean(seq_lengths):.0f}, \"\n      f\"median={np.median(seq_lengths):.0f}, range=[{min(seq_lengths)}-{max(seq_lengths)}]\")\n\nprint(\"\\n✅ Data loading complete! Ready for feature extraction and training.\")\nprint(f\"   Total proteins: {len(train_seqs):,}\")\nprint(f\"   Total annotations: {len(train_terms):,}\")\nprint(f\"   Total GO terms: {train_terms['term'].nunique():,}\")","metadata":{"trusted":true,"jupyter":{"source_hidden":true}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================================\n# PACKAGE INSTALLATION\n# ============================================================================\n# Install any missing packages\n!pip install -q obonet biopython","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\nimport numpy as np\nimport pandas as pd\nimport os\n\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# ============================================================================\n# CORE IMPORTS\n# ============================================================================\nimport gc\nimport sys\nimport subprocess\nfrom pathlib import Path\nfrom collections import defaultdict, Counter\nfrom tqdm.auto import tqdm\nimport warnings\nwarnings.filterwarnings('ignore')\n\n# ============================================================================\n# CONFIG (ADVANCED IMPROVEMENTS)\n# ============================================================================\nCONFIG = {\n    'MAX_TRAIN_SAMPLES': 80000,        # IMPROVED: 5000 → 8000 (60% more data!)\n    'TOP_K_LABELS': 750,\n    'RANDOM_SEED': 42,\n    'HIDDEN_DIMS': [1024, 512, 256],\n    'DROPOUT_RATE': 0.25,             # IMPROVED: 0.3 → 0.25 (less dropout for more data)\n    'EPOCHS': 12,                     # IMPROVED: 10 → 12 epochs\n    'BATCH_SIZE': 64,                 # IMPROVED: 32 → 64 (faster training)\n    'LEARNING_RATE': 1e-3,\n    'PREDICT_BATCH_SIZE': 128,        # IMPROVED: 64 → 128 (faster prediction)\n    'MIN_CONFIDENCE': 0.08,           # IMPROVED: 0.10 → 0.08 (even more predictions)\n    'MAX_PREDS_PER_PROTEIN': 200,    # IMPROVED: 150 → 200\n    'LABEL_SMOOTHING': 0.1,          # NEW: Add label smoothing\n    'TEMPERATURE': 1.5,               # NEW: Temperature scaling for predictions\n    'MERGE_WITH_BLAST': True,\n    'BLAST_WEIGHT': 0.60,             # IMPROVED: 0.65 → 0.60\n    'DL_WEIGHT': 0.40,                # IMPROVED: 0.35 → 0.40 (more trust in DL)\n    'BASE_PATH': '/kaggle/input/cafa-6-protein-function-prediction',\n    'ESM2_PATH': 'cafa-5-ems-2-embeddings-numpy',\n    'BLAST_PATH': '/kaggle/input/blast-quick-sprof-zero-pred/submission.tsv',\n}\n\nprint(\"=\"*80)\nprint(\"🚀 CAFA 6 - ADVANCED OPTIMIZATION\")\nprint(\"=\"*80)\n\ndef install(package):\n    subprocess.check_call([sys.executable, \"-m\", \"pip\", \"install\", \"-q\", package])\n\ntry:\n    import torch\n    import torch.nn as nn\n    import torch.nn.functional as F\nexcept:\n    install('torch')\n    import torch\n    import torch.nn as nn\n    import torch.nn.functional as F\n\nnp.random.seed(CONFIG['RANDOM_SEED'])\ntorch.manual_seed(CONFIG['RANDOM_SEED'])\nif torch.cuda.is_available():\n    torch.cuda.manual_seed(CONFIG['RANDOM_SEED'])\ndevice = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\nprint(f\"Device: {device}\")\n\n# ============================================================================\n# LOAD DATA\n# ============================================================================\nBASE = Path(CONFIG['BASE_PATH'])\nTRAIN_DIR = BASE / 'Train'\n\nprint(\"\\n[1/5] Loading annotations...\")\ntrain_terms_df = pd.read_csv(TRAIN_DIR / 'train_terms.tsv', sep='\\t', \n                              header=None, names=['protein', 'term', 'ontology'])\nprotein_to_terms = train_terms_df.groupby('protein')['term'].apply(list).to_dict()\n\nterm_counts = Counter()\nfor terms in protein_to_terms.values():\n    term_counts.update(terms)\ntop_terms = [t for t, _ in term_counts.most_common(CONFIG['TOP_K_LABELS'])]\n\nfor protein in protein_to_terms:\n    protein_to_terms[protein] = [t for t in protein_to_terms[protein] if t in top_terms]\n\nterm_to_idx = {term: idx for idx, term in enumerate(top_terms)}\nprint(f\"   ✓ {len(protein_to_terms)} proteins, {len(top_terms)} terms\")\n\nprint(\"\\n[2/5] Loading ESM2 embeddings...\")\nesm2_base = f\"/kaggle/input/{CONFIG['ESM2_PATH']}\"\ntrain_ids = np.load(f\"{esm2_base}/train_ids.npy\", allow_pickle=True)\ntrain_embeds = np.load(f\"{esm2_base}/train_embeddings.npy\")\ntest_ids = np.load(f\"{esm2_base}/test_ids.npy\", allow_pickle=True)\ntest_embeds = np.load(f\"{esm2_base}/test_embeddings.npy\")\n\ntrain_dict = {str(pid): emb for pid, emb in zip(train_ids, train_embeds)}\ntest_dict = {str(pid): emb for pid, emb in zip(test_ids, test_embeds)}\nembed_dim = train_embeds.shape[1]\nprint(f\"   ✓ Dim: {embed_dim}\")\n\ndel train_ids, train_embeds, test_ids, test_embeds\ngc.collect()\n\nprint(\"\\n[3/5] Preparing training data...\")\nvalid_proteins = [p for p in protein_to_terms.keys() if p in train_dict][:CONFIG['MAX_TRAIN_SAMPLES']]\n\nfrom sklearn.preprocessing import MultiLabelBinarizer\nfrom sklearn.model_selection import train_test_split\n\nmlb = MultiLabelBinarizer(classes=range(len(top_terms)))\ny_labels = [[term_to_idx[t] for t in protein_to_terms.get(p, []) if t in term_to_idx] \n            for p in valid_proteins]\ny_encoded = mlb.fit_transform(y_labels)\n\n# Apply label smoothing\nif CONFIG['LABEL_SMOOTHING'] > 0:\n    y_encoded = y_encoded.astype(float)\n    y_encoded = y_encoded * (1 - CONFIG['LABEL_SMOOTHING']) + CONFIG['LABEL_SMOOTHING'] / len(top_terms)\n\ntrain_proteins, val_proteins, y_train, y_val = train_test_split(\n    valid_proteins, y_encoded, test_size=0.15, random_state=CONFIG['RANDOM_SEED']\n)\nprint(f\"   ✓ Train: {len(train_proteins)}, Val: {len(val_proteins)}\")\n\n# ============================================================================\n# MODEL\n# ============================================================================\nprint(\"\\n[4/5] Building model...\")\n\nclass ImprovedProteinModel(nn.Module):\n    \"\"\"Improved 3-layer architecture with BatchNorm\"\"\"\n    def __init__(self, input_dim, output_dim, hidden_dims, dropout):\n        super().__init__()\n        layers = []\n        prev_dim = input_dim\n        for dim in hidden_dims:\n            layers.extend([\n                nn.Linear(prev_dim, dim),\n                nn.BatchNorm1d(dim),\n                nn.ReLU(),\n                nn.Dropout(dropout)\n            ])\n            prev_dim = dim\n        self.encoder = nn.Sequential(*layers)\n        self.output = nn.Linear(prev_dim, output_dim)\n    \n    def forward(self, x):\n        return torch.sigmoid(self.output(self.encoder(x)))\n\nmodel = ImprovedProteinModel(embed_dim, len(top_terms), CONFIG['HIDDEN_DIMS'], CONFIG['DROPOUT_RATE']).to(device)\nprint(f\"   ✓ Parameters: {sum(p.numel() for p in model.parameters()):,}\")\n\n# ============================================================================\n# TRAINING\n# ============================================================================\nprint(\"\\n\" + \"=\"*80)\nprint(\"TRAINING (12 EPOCHS WITH LABEL SMOOTHING)\")\nprint(\"=\"*80)\n\ncriterion = nn.BCELoss()\noptimizer = torch.optim.Adam(model.parameters(), lr=CONFIG['LEARNING_RATE'], weight_decay=1e-5)  # Added L2 regularization\nscheduler = torch.optim.lr_scheduler.ReduceLROnPlateau(optimizer, mode='min', factor=0.5, patience=2, verbose=False)\n\nbest_val_loss = float('inf')\n\nfor epoch in range(CONFIG['EPOCHS']):\n    model.train()\n    indices = np.random.permutation(len(train_proteins))\n    epoch_loss = 0\n    n_batches = 0\n    \n    for i in range(0, len(indices), CONFIG['BATCH_SIZE']):\n        batch_idx = indices[i:i + CONFIG['BATCH_SIZE']]\n        batch_proteins = [train_proteins[j] for j in batch_idx]\n        \n        X_batch = torch.FloatTensor([train_dict[p] for p in batch_proteins]).to(device)\n        y_batch = torch.FloatTensor(y_train[batch_idx]).to(device)\n        \n        optimizer.zero_grad()\n        outputs = model(X_batch)\n        loss = criterion(outputs, y_batch)\n        loss.backward()\n        torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm=1.0)  # Gradient clipping\n        optimizer.step()\n        \n        epoch_loss += loss.item()\n        n_batches += 1\n    \n    model.eval()\n    val_loss = 0\n    val_batches = 0\n    \n    with torch.no_grad():\n        for i in range(0, len(val_proteins), CONFIG['BATCH_SIZE']):\n            batch_proteins = val_proteins[i:i + CONFIG['BATCH_SIZE']]\n            X_batch = torch.FloatTensor([train_dict[p] for p in batch_proteins]).to(device)\n            y_batch = torch.FloatTensor(y_val[i:i + CONFIG['BATCH_SIZE']]).to(device)\n            \n            outputs = model(X_batch)\n            loss = criterion(outputs, y_batch)\n            val_loss += loss.item()\n            val_batches += 1\n    \n    train_loss_avg = epoch_loss/n_batches\n    val_loss_avg = val_loss/val_batches\n    \n    scheduler.step(val_loss_avg)\n    \n    if val_loss_avg < best_val_loss:\n        best_val_loss = val_loss_avg\n        print(f\"Epoch {epoch+1}: Train={train_loss_avg:.4f}, Val={val_loss_avg:.4f} ⭐ NEW BEST\")\n    else:\n        print(f\"Epoch {epoch+1}: Train={train_loss_avg:.4f}, Val={val_loss_avg:.4f}\")\n    \n    gc.collect()\n\nprint(f\"\\n✅ Best Validation Loss: {best_val_loss:.4f}\")\n\n# ============================================================================\n# PREDICTIONS (WITH TEMPERATURE SCALING)\n# ============================================================================\nprint(\"\\n\" + \"=\"*80)\nprint(\"PREDICTIONS (WITH TEMPERATURE SCALING)\")\nprint(\"=\"*80)\n\nmodel.eval()\ntest_protein_ids = list(test_dict.keys())\n\nn_predictions = 0\nwith open('temp_dl.tsv', 'w', newline='') as f:\n    with torch.no_grad():\n        for start in tqdm(range(0, len(test_protein_ids), CONFIG['PREDICT_BATCH_SIZE']), desc=\"Predicting\"):\n            batch_ids = test_protein_ids[start:start + CONFIG['PREDICT_BATCH_SIZE']]\n            X_batch = torch.FloatTensor([test_dict[p] for p in batch_ids]).to(device)\n            \n            # Get raw logits before sigmoid\n            logits = model.encoder(X_batch)\n            logits = model.output(logits)\n            \n            # Apply temperature scaling\n            outputs = torch.sigmoid(logits / CONFIG['TEMPERATURE']).cpu().numpy()\n            \n            for i, pid in enumerate(batch_ids):\n                probs = outputs[i]\n                top_indices = np.argsort(probs)[::-1][:CONFIG['MAX_PREDS_PER_PROTEIN']]\n                confident_indices = [idx for idx in top_indices if probs[idx] > CONFIG['MIN_CONFIDENCE']]\n                \n                for idx in confident_indices:\n                    line = f\"{pid}\\t{top_terms[idx]}\\t{min(probs[idx], 0.999):.3f}\\n\"\n                    f.write(line)\n                    n_predictions += 1\n            \n            del X_batch, outputs, logits\n            if start % 1000 == 0:\n                gc.collect()\n\nprint(f\"✓ Generated {n_predictions:,} predictions\")\n\ndel model, train_dict\ngc.collect()\n\n# ============================================================================\n# ENSEMBLE\n# ============================================================================\nprint(\"\\n\" + \"=\"*80)\nprint(\"ENSEMBLE\")\nprint(\"=\"*80)\n\nif CONFIG['MERGE_WITH_BLAST'] and os.path.exists(CONFIG['BLAST_PATH']):\n    print(\"Merging with BLAST...\")\n    \n    dl_df = pd.read_csv('temp_dl.tsv', sep='\\t', header=None, names=['Id', 'GO term', 'Confidence'])\n    dl_df['Confidence'] = dl_df['Confidence'].astype(float) * CONFIG['DL_WEIGHT']\n    \n    blast_chunks = []\n    for chunk in pd.read_csv(CONFIG['BLAST_PATH'], sep='\\t', header=None,\n                             names=['Id', 'GO term', 'Confidence'], chunksize=1000000):\n        chunk['Confidence'] = chunk['Confidence'].astype(float) * CONFIG['BLAST_WEIGHT']\n        blast_chunks.append(chunk)\n        if len(blast_chunks) >= 10:\n            blast_df = pd.concat(blast_chunks, ignore_index=True)\n            blast_chunks = [blast_df]\n            gc.collect()\n    \n    blast_df = pd.concat(blast_chunks, ignore_index=True)\n    print(f\"  DL={len(dl_df):,}, BLAST={len(blast_df):,}\")\n    \n    ensemble_df = pd.concat([dl_df, blast_df], ignore_index=True)\n    del dl_df, blast_df, blast_chunks\n    gc.collect()\n    \n    ensemble_df = ensemble_df.groupby(['Id', 'GO term'], as_index=False)['Confidence'].sum()\n    ensemble_df = ensemble_df.sort_values(['Id', 'Confidence'], ascending=[True, False])\n    ensemble_df = ensemble_df.groupby('Id').head(1500)\n    \n    import csv\n    ensemble_df.to_csv('submission.tsv', sep='\\t', header=False, index=False,\n                      quoting=csv.QUOTE_NONE, escapechar='\\\\', lineterminator='\\n')\n    print(f\"✓ Saved {len(ensemble_df):,} predictions\")\n    \n    del ensemble_df\nelse:\n    os.rename('temp_dl.tsv', 'submission.tsv')\n\ngc.collect()\n\n# ============================================================================\n# CLEANUP\n# ============================================================================\nprint(\"\\n[CLEANUP]\")\nfor fname in os.listdir('.'):\n    if 'submission' in fname.lower() and fname != 'submission.tsv':\n        try:\n            os.remove(fname)\n            print(f\"  ✗ Removed {fname}\")\n        except:\n            pass\n\nif os.path.exists('temp_dl.tsv'):\n    os.remove('temp_dl.tsv')\n    print(f\"  ✗ Removed temp_dl.tsv\")\n\n# ============================================================================\n# VERIFICATION\n# ============================================================================\nprint(\"\\n\" + \"=\"*80)\nprint(\"✅ VERIFICATION\")\nprint(\"=\"*80)\n\nif os.path.exists('submission.tsv'):\n    with open('submission.tsv', 'r') as f:\n        lines = f.readlines()\n    \n    n_lines = len(lines)\n    file_size_mb = os.path.getsize('submission.tsv') / (1024*1024)\n    \n    print(f\"📋 File: submission.tsv\")\n    print(f\"📊 Lines: {n_lines:,}\")\n    print(f\"💾 Size: {file_size_mb:.1f} MB\")\n    \n    print(\"\\n✅ First 5 lines:\")\n    for i in range(min(5, len(lines))):\n        line = lines[i].strip()\n        parts = line.split('\\t')\n        if len(parts) == 3:\n            print(f\"  ✓ {parts[0]:<15} {parts[1]:<12} {parts[2]}\")\n        else:\n            print(f\"  ✗ Line {i+1}: {len(parts)} columns (expected 3)\")\n    \n    first_line = lines[0].strip()\n    parts = first_line.split('\\t')\n    \n    print(\"\\n🔍 Format check:\")\n    all_good = True\n    \n    if len(parts) == 3:\n        print(f\"  ✓ Columns: 3\")\n    else:\n        print(f\"  ✗ Columns: {len(parts)} (expected 3)\")\n        all_good = False\n    \n    if parts[1].startswith('GO:'):\n        print(f\"  ✓ GO term format\")\n    else:\n        print(f\"  ✗ GO term: {parts[1]}\")\n        all_good = False\n    \n    try:\n        conf = float(parts[2])\n        if 0 < conf <= 1:\n            print(f\"  ✓ Confidence: {conf}\")\n        else:\n            print(f\"  ✗ Confidence out of range: {conf}\")\n            all_good = False\n    except:\n        print(f\"  ✗ Confidence not a number: {parts[2]}\")\n        all_good = False\n    \n    if '\\n' in lines[0]:\n        print(f\"  ✓ Newlines present\")\n    else:\n        print(f\"  ✗ No newlines!\")\n        all_good = False\n    \n    if all_good:\n        print(\"\\n✅✅✅ FORMAT IS CORRECT! ✅✅✅\")\n    else:\n        print(\"\\n⚠️  FORMAT HAS ISSUES!\")\n        \n    print(\"\\n📁 Files in working directory:\")\n    for fname in sorted(os.listdir('.')):\n        if not fname.startswith('.'):\n            size = os.path.getsize(fname) / (1024*1024)\n            print(f\"  {fname} ({size:.1f} MB)\")\nelse:\n    print(\"❌ ERROR: submission.tsv not found!\")\n\nprint(\"\\n\" + \"=\"*80)","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}