{"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":"none","dataSources":[{"sourceId":116062,"databundleVersionId":14084779,"sourceType":"competition"},{"sourceId":5499219,"sourceType":"datasetVersion","datasetId":3167603}],"dockerImageVersionId":31153,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"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":"# ============================================================================\n# CAFA 6 PROTEIN FUNCTION PREDICTION - COMPLETE FIXED PIPELINE\n# ============================================================================\n\nimport pandas as pd\nimport numpy as np\nfrom pathlib import Path\nfrom collections import Counter\nimport warnings\nwarnings.filterwarnings('ignore')\n\n# BioPython for sequence parsing\nfrom Bio import SeqIO\n\n# Modeling libraries\nfrom sklearn.preprocessing import MultiLabelBinarizer, StandardScaler\nfrom sklearn.multiclass import OneVsRestClassifier\nfrom sklearn.linear_model import LogisticRegression\n\n# Network analysis for GO\nimport networkx as nx\nimport obonet\n\nprint(\"=\"*80)\nprint(\"CAFA 6 PROTEIN FUNCTION PREDICTION - COMPLETE FIXED PIPELINE\")\nprint(\"=\"*80)\n\n# ============================================================================\n# 1. CONFIGURATION & PATHS\n# ============================================================================\nBASE = Path('/kaggle/input/cafa-6-protein-function-prediction')\nTRAIN_DIR = BASE / 'Train'\nTEST_DIR = BASE / 'Test'\n\n# Model parameters\nBATCH_SIZE = 2500\nPREDICTION_THRESHOLD = 0.012\nMAX_PREDS_PER_PROTEIN = 1500\nMAX_PREDS_PER_ONT = 500\n\n# ============================================================================\n# 2. LOAD GO ONTOLOGY\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# ============================================================================\n# 3. LOAD TRAINING DATA\n# ============================================================================\nprint(\"\\n[2/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\"   ✓ Loaded {len(train_terms)} annotations for {train_terms['protein'].nunique()} proteins\")\n\n# Load training sequences\nprint(\"\\n[3/9] Loading training sequences...\")\ntrain_seqs = {}\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    train_seqs[pid] = str(rec.seq)\n    if len(train_seqs) % 20000 == 0:\n        print(f\"      Loaded {len(train_seqs):,} sequences...\")\n\nprint(f\"   ✓ Loaded {len(train_seqs)} training sequences\")\n\n# ============================================================================\n# 4. FEATURE EXTRACTION WITH COMPLETE NAN PROTECTION\n# ============================================================================\nprint(\"\\n[4/9] Defining unified feature extraction...\")\n\ndef extract_sequence_features(seq):\n    \"\"\"Extract features from amino acid sequence with full NaN protection\"\"\"\n    # Define feature size - must be consistent!\n    FEATURE_SIZE = 25\n    \n    if not seq or len(seq) == 0:\n        return np.zeros(FEATURE_SIZE)\n    \n    try:\n        # Amino acid composition (20 features)\n        aa_counts = Counter(seq)\n        length = len(seq)\n        \n        # Standard amino acids\n        aa_freq = np.array([aa_counts.get(aa, 0) / length \n                           for aa in 'ACDEFGHIKLMNPQRSTVWY'])\n        \n        # Basic properties (5 features)\n        hydrophobic = sum(aa_counts.get(aa, 0) for aa in 'AILMFWYV') / length\n        charged = sum(aa_counts.get(aa, 0) for aa in 'DEKR') / length\n        polar = sum(aa_counts.get(aa, 0) for aa in 'STNQ') / length\n        aromatic = sum(aa_counts.get(aa, 0) for aa in 'FWY') / length\n        \n        # Combine all features\n        features = np.concatenate([\n            aa_freq,  # 20 features\n            [np.log1p(length), hydrophobic, charged, polar, aromatic]  # 5 features\n        ])\n        \n        # Replace any NaN or inf values with 0\n        features = np.nan_to_num(features, nan=0.0, posinf=0.0, neginf=0.0)\n        \n        # Ensure exactly FEATURE_SIZE features\n        if len(features) != FEATURE_SIZE:\n            print(f\"Warning: Feature size mismatch. Expected {FEATURE_SIZE}, got {len(features)}\")\n            features = np.resize(features, FEATURE_SIZE)\n        \n        return features\n        \n    except Exception as e:\n        print(f\"Error in feature extraction: {e}\")\n        return np.zeros(FEATURE_SIZE)\n\n# ============================================================================\n# 5. EXTRACT FEATURES FOR TRAINING\n# ============================================================================\nprint(\"\\n[5/9] Extracting features for training proteins...\")\nX_train_list = []\ny_train_proteins = []\n\nfor i, (pid, seq) in enumerate(train_seqs.items()):\n    features = extract_sequence_features(seq)\n    X_train_list.append(features)\n    y_train_proteins.append(pid)\n    \n    if (i + 1) % 20000 == 0:\n        print(f\"      Processed {i+1:,} proteins...\")\n\nX_train = np.array(X_train_list)\nprint(f\"   ✓ Feature matrix shape: {X_train.shape}\")\n\n# Final NaN check on training data\nif np.any(np.isnan(X_train)) or np.any(np.isinf(X_train)):\n    print(\"   Warning: NaN or inf values found in training features, replacing...\")\n    X_train = np.nan_to_num(X_train, nan=0.0, posinf=0.0, neginf=0.0)\n\n# Standardize features\nprint(\"   Standardizing features...\")\nscaler = StandardScaler()\nX_train_scaled = scaler.fit_transform(X_train)\n\n# ============================================================================\n# 6. PREPARE LABELS BY ONTOLOGY\n# ============================================================================\nprint(\"\\n[6/9] Preparing labels by ontology...\")\n\nontologies = {'F': 'MFO', 'P': 'BPO', 'C': 'CCO'}\nmlb_dict = {}\ny_train_dict = {}\n\nfor ont_code, ont_name in ontologies.items():\n    print(f\"   Processing {ont_name}...\")\n    \n    ont_terms = train_terms[train_terms['ontology'] == ont_code]\n    protein_terms = ont_terms.groupby('protein')['term'].apply(list).to_dict()\n    labels_list = [protein_terms.get(pid, []) for pid in y_train_proteins]\n    \n    mlb = MultiLabelBinarizer(sparse_output=True)\n    y_ont = mlb.fit_transform(labels_list)\n    \n    mlb_dict[ont_code] = mlb\n    y_train_dict[ont_code] = y_ont\n    \n    print(f\"      {ont_name}: {y_ont.shape[1]} unique terms\")\n\n# ============================================================================\n# 7. TRAIN MODELS\n# ============================================================================\nprint(\"\\n[7/9] Training models...\")\n\nmodels = {}\n\nfor ont_code, ont_name in ontologies.items():\n    print(f\"   Training {ont_name} model...\")\n    \n    y_ont = y_train_dict[ont_code]\n    \n    if y_ont.shape[1] == 0:\n        print(f\"      Skipping {ont_name} (no terms)\")\n        continue\n    \n    # Simple, robust configuration\n    model = OneVsRestClassifier(\n        LogisticRegression(\n            max_iter=1000,\n            solver='lbfgs',\n            C=1.0,\n            random_state=42,\n            n_jobs=1\n        ),\n        n_jobs=-1\n    )\n    \n    model.fit(X_train_scaled, y_ont)\n    models[ont_code] = model\n    \n    print(f\"      ✓ {ont_name} model trained\")\n\nprint(f\"   ✓ All {len(models)} models trained successfully\")\n\n# ============================================================================\n# 8. LOAD TEST DATA AND MAKE PREDICTIONS\n# ============================================================================\nprint(\"\\n[8/9] Loading test sequences...\")\n\ntest_seqs = {}\ntest_proteins = []\n\nfor rec in SeqIO.parse(TEST_DIR / 'testsuperset.fasta', 'fasta'):\n    pid = rec.id.split('|')[1] if '|' in rec.id else rec.id\n    test_seqs[pid] = str(rec.seq)\n    test_proteins.append(pid)\n\nprint(f\"   ✓ Loaded {len(test_seqs):,} test sequences\")\n\nprint(\"\\n[9/9] Making predictions...\")\n\n# Process predictions with robust error handling\nsubmission_list = []\nskipped_proteins = 0\n\nfor batch_start in range(0, len(test_proteins), BATCH_SIZE):\n    batch_end = min(batch_start + BATCH_SIZE, len(test_proteins))\n    batch_pids = test_proteins[batch_start:batch_end]\n    \n    # Extract features for batch with validation\n    X_batch = []\n    valid_pids = []\n    \n    for pid in batch_pids:\n        features = extract_sequence_features(test_seqs[pid])\n        \n        # Validate features\n        if features.shape[0] != X_train.shape[1]:\n            print(f\"   Warning: Protein {pid} has wrong feature dimension, skipping...\")\n            skipped_proteins += 1\n            continue\n            \n        # Check for NaN or inf\n        if np.any(np.isnan(features)) or np.any(np.isinf(features)):\n            features = np.nan_to_num(features, nan=0.0, posinf=0.0, neginf=0.0)\n        \n        X_batch.append(features)\n        valid_pids.append(pid)\n    \n    if not X_batch:\n        print(f\"   Warning: No valid proteins in batch {batch_start}-{batch_end}\")\n        continue\n    \n    X_batch = np.array(X_batch)\n    \n    # Final NaN check before scaling\n    if np.any(np.isnan(X_batch)) or np.any(np.isinf(X_batch)):\n        print(f\"   Warning: NaN/inf found in batch, cleaning...\")\n        X_batch = np.nan_to_num(X_batch, nan=0.0, posinf=0.0, neginf=0.0)\n    \n    try:\n        X_batch_scaled = scaler.transform(X_batch)\n    except Exception as e:\n        print(f\"   Error scaling batch {batch_start}-{batch_end}: {e}\")\n        continue\n    \n    # Make predictions for each ontology\n    for ont_code in models:\n        model = models[ont_code]\n        mlb = mlb_dict[ont_code]\n        \n        try:\n            y_pred_proba = model.predict_proba(X_batch_scaled)\n            \n            for i, pid in enumerate(valid_pids):\n                probs = y_pred_proba[i]\n                \n                # Get predictions above threshold\n                top_indices = np.where(probs > PREDICTION_THRESHOLD)[0]\n                \n                # Sort by probability and limit\n                if len(top_indices) > 0:\n                    sorted_indices = top_indices[np.argsort(probs[top_indices])[::-1]]\n                    sorted_indices = sorted_indices[:MAX_PREDS_PER_ONT]\n                    \n                    for idx in sorted_indices:\n                        term = mlb.classes_[idx]\n                        score = probs[idx]\n                        submission_list.append((pid, term, score))\n                        \n        except Exception as e:\n            print(f\"   Error predicting {ont_code} for batch {batch_start}-{batch_end}: {e}\")\n            continue\n    \n    if (batch_end % 10000 == 0) or (batch_end == len(test_proteins)):\n        print(f\"      Processed {batch_end:,}/{len(test_proteins):,} proteins...\")\n\nprint(f\"   ✓ Generated {len(submission_list):,} predictions\")\nif skipped_proteins > 0:\n    print(f\"   ⚠ Skipped {skipped_proteins} proteins due to feature extraction issues\")\n\n# Create submission dataframe\nsubmission_df = pd.DataFrame(submission_list, columns=['protein', 'term', 'score'])\n\n# Sort and limit to 1500 per protein\nsubmission_df = submission_df.sort_values(['protein', 'score'], ascending=[True, False])\nsubmission_df = submission_df.groupby('protein').head(MAX_PREDS_PER_PROTEIN).reset_index(drop=True)\n\n# Format scores with validation\ndef format_score(x):\n    try:\n        score = float(x)\n        if score <= 0:\n            return \"0.001\"\n        elif score > 1:\n            return \"1.0\"\n        else:\n            return f\"{score:.3g}\"\n    except:\n        return \"0.001\"\n\nsubmission_df['score'] = submission_df['score'].apply(format_score)\n\n# Save submission\nsubmission_df.to_csv('submission.tsv', sep='\\t', index=False, header=False)\n\nprint(f\"\\n   ✓ Final submission saved: {len(submission_df):,} predictions\")\nprint(\"\\n   Sample predictions:\")\nprint(submission_df.head(10).to_string(index=False))\n\n# Calculate statistics\nprotein_counts = submission_df['protein'].value_counts()\nprint(f\"\\n   Prediction statistics:\")\nprint(f\"      Mean predictions per protein: {protein_counts.mean():.1f}\")\nprint(f\"      Median predictions per protein: {protein_counts.median():.0f}\")\nprint(f\"      Max predictions per protein: {protein_counts.max()}\")\nprint(f\"      Min predictions per protein: {protein_counts.min()}\")\nprint(f\"      Total unique proteins: {len(protein_counts):,}\")\n\nprint(\"\\n\" + \"=\"*80)\nprint(\"✅ COMPLETE FIXED PIPELINE FINISHED SUCCESSFULLY!\")\nprint(\"=\"*80)\nprint(f\"\\nFeatures: 25 robust features with NaN protection\")\nprint(f\"Models: Ontology-specific LogisticRegression\")\nprint(f\"Output: submission.tsv with {len(submission_df):,} predictions\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}