{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":116062,"databundleVersionId":14084779,"sourceType":"competition"}],"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,"execution":{"iopub.status.busy":"2025-10-22T15:39:48.651996Z","iopub.execute_input":"2025-10-22T15:39:48.652351Z","iopub.status.idle":"2025-10-22T15:39:48.669573Z","shell.execute_reply.started":"2025-10-22T15:39:48.652323Z","shell.execute_reply":"2025-10-22T15:39:48.668579Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# 1. Full EDA of CAFA 6","metadata":{}},{"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,"execution":{"iopub.status.busy":"2025-10-23T13:17:52.652074Z","iopub.execute_input":"2025-10-23T13:17:52.653283Z","iopub.status.idle":"2025-10-23T13:18:27.358055Z","shell.execute_reply.started":"2025-10-23T13:17:52.653253Z","shell.execute_reply":"2025-10-23T13:18:27.357148Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# 02. Best prediction","metadata":{}},{"cell_type":"code","source":"# ============================================================================\n# CAFA 6 PROTEIN FUNCTION PREDICTION - OPTIMIZED PIPELINE (FIXED VERSION)\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\nfrom sklearn.feature_selection import SelectKBest, f_classif\n\n# Memory optimization\nimport gc\nimport time\n\nprint(\"=\"*80)\nprint(\"CAFA 6 PROTEIN FUNCTION PREDICTION - OPTIMIZED PIPELINE (FIXED)\")\nprint(\"=\"*80)\n\n# ============================================================================\n# 1. CONFIGURATION & PATHS - OPTIMIZED FOR SPEED\n# ============================================================================\nBASE = Path('/kaggle/input/cafa-6-protein-function-prediction')\nTRAIN_DIR = BASE / 'Train'\nTEST_DIR = BASE / 'Test'\n\n# ✅ OPTIMIZED PARAMETERS FOR FASTER EXECUTION\nSAMPLE_FRACTION = 0.05  # Reduced from 0.7 to 0.05 (5% data only)\nBATCH_SIZE = 500  # Reduced batch size\nPREDICTION_THRESHOLD = 0.02  # Slightly higher for better precision\nMAX_PREDS_PER_PROTEIN = 800\nMAX_PREDS_PER_ONT = 200\n\n# ✅ NEW: Limit terms per ontology for faster training\nMAX_TERMS_PER_ONTOLOGY = 200  # Maximum terms per ontology\n\nprint(f\"🔧 OPTIMIZED CONFIGURATION:\")\nprint(f\"   - Sample Fraction: {SAMPLE_FRACTION*100}%\")\nprint(f\"   - Max Terms per Ontology: {MAX_TERMS_PER_ONTOLOGY}\")\nprint(f\"   - Batch Size: {BATCH_SIZE}\")\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. OPTIMIZED DATA LOADING\n# ============================================================================\nprint(\"\\n[2/9] Loading training data with optimization...\")\n\n# Load only necessary columns\ntrain_terms = pd.read_csv(TRAIN_DIR / 'train_terms.tsv', sep='\\t', \n                          names=['protein', 'term', 'ontology'],\n                          usecols=[0, 1, 2])  # Specify columns to save memory\n\n# ✅ IMPROVED: Better sampling method\nunique_proteins = train_terms['protein'].unique()\nsampled_proteins = np.random.choice(unique_proteins, \n                                   size=int(len(unique_proteins) * SAMPLE_FRACTION), \n                                   replace=False)\n\ntrain_terms = train_terms[train_terms['protein'].isin(sampled_proteins)]\n\nprint(f\"   ✓ Loaded {len(train_terms)} annotations for {train_terms['protein'].nunique()} proteins\")\n\n# ============================================================================\n# 4. ENHANCED FEATURE EXTRACTION\n# ============================================================================\nprint(\"\\n[3/9] Loading training sequences with optimized features...\")\n\ndef optimized_feature_extraction(seq):\n    \"\"\"Enhanced feature extraction with better biological insights\"\"\"\n    FEATURE_SIZE = 30  # Increased for better representation\n    \n    if not seq or len(seq) == 0:\n        return np.zeros(FEATURE_SIZE)\n    \n    try:\n        length = len(seq)\n        aa_counts = Counter(seq)\n        total_aa = sum(aa_counts.values())\n        \n        # 1. Amino Acid Composition (20 features)\n        aa_list = 'ACDEFGHIKLMNPQRSTVWY'\n        aa_freq = np.array([aa_counts.get(aa, 0) / total_aa for aa in aa_list])\n        \n        # 2. Physicochemical Properties (8 features)\n        # Hydrophobicity groups\n        very_hydrophobic = sum(aa_counts.get(aa, 0) for aa in 'AILMFWYV') / total_aa\n        hydrophobic = sum(aa_counts.get(aa, 0) for aa in 'CGT') / total_aa\n        \n        # Charge properties\n        positive_charged = sum(aa_counts.get(aa, 0) for aa in 'KRH') / total_aa\n        negative_charged = sum(aa_counts.get(aa, 0) for aa in 'DE') / total_aa\n        \n        # Structural properties\n        polar = sum(aa_counts.get(aa, 0) for aa in 'STNQ') / total_aa\n        aromatic = sum(aa_counts.get(aa, 0) for aa in 'FWY') / total_aa\n        small = sum(aa_counts.get(aa, 0) for aa in 'AG') / total_aa\n        proline_content = aa_counts.get('P', 0) / total_aa\n        \n        # 3. Sequence Properties (2 features)\n        seq_complexity = len(set(seq)) / length if length > 0 else 0\n        molecular_weight_approx = length * 110  # Average AA molecular weight\n        \n        # Combine all features\n        features = np.concatenate([\n            aa_freq,  # 20 features\n            [\n                np.log1p(length), very_hydrophobic, hydrophobic,\n                positive_charged, negative_charged, polar, aromatic,\n                small, proline_content, seq_complexity,\n                np.log1p(molecular_weight_approx)\n            ]  # 11 features\n        ])\n        \n        # Ensure correct size and handle NaN\n        features = np.nan_to_num(features, nan=0.0, posinf=0.0, neginf=0.0)\n        \n        if len(features) != FEATURE_SIZE:\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. MEMORY-EFFICIENT FEATURE PROCESSING\n# ============================================================================\nprint(\"\\n[4/9] Processing training features in batches...\")\n\ndef process_proteins_batch(protein_seqs, batch_size=5000):\n    \"\"\"Process proteins in batches to save memory\"\"\"\n    all_features = []\n    all_proteins = []\n    \n    proteins = list(protein_seqs.items())\n    \n    for i in range(0, len(proteins), batch_size):\n        batch = proteins[i:i + batch_size]\n        batch_features = []\n        batch_proteins = []\n        \n        for pid, seq in batch:\n            features = optimized_feature_extraction(seq)\n            batch_features.append(features)\n            batch_proteins.append(pid)\n        \n        all_features.extend(batch_features)\n        all_proteins.extend(batch_proteins)\n        \n        # Clear memory\n        del batch, batch_features\n        gc.collect()\n        \n        print(f\"      Processed {min(i + batch_size, len(proteins)):,}/{len(proteins):,} proteins...\")\n    \n    return np.array(all_features), all_proteins\n\n# Load and process training sequences\ntrain_seqs = {}\ntrain_proteins_processed = 0\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    # Only process sampled proteins\n    if pid in sampled_proteins:\n        train_seqs[pid] = str(rec.seq)\n        train_proteins_processed += 1\n        \n    if train_proteins_processed >= len(sampled_proteins):\n        break\n\nprint(f\"   ✓ Loaded {len(train_seqs)} training sequences\")\n\n# Process in batches\nX_train, y_train_proteins = process_proteins_batch(train_seqs, batch_size=5000)\nprint(f\"   ✓ Feature matrix shape: {X_train.shape}\")\n\n# Clean memory\ndel train_seqs\ngc.collect()\n\n# Standardize features\nprint(\"   Standardizing features...\")\nscaler = StandardScaler()\nX_train_scaled = scaler.fit_transform(X_train)\n\n# ============================================================================\n# 6. FEATURE SELECTION FOR EFFICIENCY\n# ============================================================================\nprint(\"\\n[5/9] Applying feature selection...\")\n# ✅ REDUCED: Select top 10 features instead of 20\nfeature_selector = SelectKBest(f_classif, k=10)\nX_train_selected = feature_selector.fit_transform(X_train_scaled, np.ones(len(X_train_scaled)))\n\nprint(f\"   ✓ Reduced features from {X_train_scaled.shape[1]} to {X_train_selected.shape[1]}\")\n\n# Clean memory\ndel X_train, X_train_scaled\ngc.collect()\n\n# ============================================================================\n# 7. PREPARE LABELS WITH MEMORY OPTIMIZATION\n# ============================================================================\nprint(\"\\n[6/9] Preparing labels by ontology...\")\n\n# ✅ OPTIMIZED: Start with only MFO for faster testing\nontologies = {'F': 'MFO'}  # Start with only Molecular Function\n# ontologies = {'F': 'MFO', 'P': 'BPO', 'C': 'CCO'}  # Uncomment later for all\n\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    \n    # Only include proteins we have features for\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    # ✅ NEW: Limit number of terms for faster training\n    if y_ont.shape[1] > MAX_TERMS_PER_ONTOLOGY:\n        print(f\"      Too many terms ({y_ont.shape[1]}), limiting to {MAX_TERMS_PER_ONTOLOGY}\")\n        # Select most frequent terms\n        term_frequencies = np.array(y_ont.sum(axis=0)).flatten()\n        top_term_indices = np.argsort(term_frequencies)[-MAX_TERMS_PER_ONTOLOGY:]\n        y_ont = y_ont[:, top_term_indices]\n        # Update MLB classes\n        mlb.classes_ = mlb.classes_[top_term_indices]\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# Clean memory\ndel train_terms\ngc.collect()\n\n# ============================================================================\n# 8. EFFICIENT MODEL TRAINING WITH PROGRESS MONITORING\n# ============================================================================\nprint(\"\\n[7/9] Training optimized 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    # ✅ IMPROVED: Add timing and progress monitoring\n    start_time = time.time()\n    \n    try:\n        # Convert sparse to dense for liblinear\n        y_ont_dense = y_ont.toarray() if hasattr(y_ont, 'toarray') else y_ont\n        \n        # Memory-efficient configuration\n        model = OneVsRestClassifier(\n            LogisticRegression(\n                max_iter=100,  # Reduced iterations for speed\n                solver='liblinear',  # More memory efficient\n                C=0.8,  # Slightly more regularization\n                random_state=42,\n                class_weight='balanced',\n                verbose=0  # Set to 1 to see training progress\n            ),\n            n_jobs=1  # Single job to avoid memory issues\n        )\n        \n        print(f\"      Starting training for {y_ont_dense.shape[1]} terms...\")\n        model.fit(X_train_selected, y_ont_dense)\n        models[ont_code] = model\n        \n        end_time = time.time()\n        training_time = end_time - start_time\n        \n        print(f\"      ✓ {ont_name} model trained in {training_time:.2f} seconds\")\n        \n    except Exception as e:\n        print(f\"      ✗ Error training {ont_name}: {e}\")\n        continue\n    \n    # Clean memory after each model\n    del y_ont\n    gc.collect()\n\nprint(f\"   ✓ {len(models)} models trained successfully\")\n\n# ============================================================================\n# 9. OPTIMIZED PREDICTION PIPELINE\n# ============================================================================\nprint(\"\\n[8/9] Loading and processing test sequences...\")\n\ntest_seqs = {}\ntest_proteins = []\n\n# ✅ OPTIMIZED: Load only limited test sequences for initial testing\nMAX_TEST_SAMPLES = 1000  # Start with 1000 test samples\n\n# Process test sequences in streaming fashion\nfor i, rec in enumerate(SeqIO.parse(TEST_DIR / 'testsuperset.fasta', 'fasta')):\n    if i >= MAX_TEST_SAMPLES:\n        break\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 (limited for testing)\")\n\nprint(\"\\n[9/9] Making optimized predictions...\")\n\nsubmission_list = []\n\ndef process_prediction_batch(protein_batch, seq_dict):\n    \"\"\"Process predictions for a batch of proteins\"\"\"\n    batch_features = []\n    valid_proteins = []\n    \n    for pid in protein_batch:\n        features = optimized_feature_extraction(seq_dict[pid])\n        if len(features) == X_train_selected.shape[1]:\n            batch_features.append(features)\n            valid_proteins.append(pid)\n    \n    if not batch_features:\n        return []\n    \n    X_batch = np.array(batch_features)\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        X_batch_selected = feature_selector.transform(X_batch_scaled)\n    except Exception as e:\n        print(f\"   Error processing batch: {e}\")\n        return []\n    \n    batch_predictions = []\n    \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_selected)\n            \n            for i, pid in enumerate(valid_proteins):\n                probs = y_pred_proba[i]\n                top_indices = np.where(probs > PREDICTION_THRESHOLD)[0]\n                \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                        batch_predictions.append((pid, term, score))\n                        \n        except Exception as e:\n            print(f\"   Error in {ont_code} prediction: {e}\")\n            continue\n    \n    return batch_predictions\n\n# Process test proteins in small batches\nfor i in range(0, len(test_proteins), BATCH_SIZE):\n    batch = test_proteins[i:i + BATCH_SIZE]\n    batch_predictions = process_prediction_batch(batch, test_seqs)\n    submission_list.extend(batch_predictions)\n    \n    if (i + BATCH_SIZE) % 500 == 0 or (i + BATCH_SIZE) >= len(test_proteins):\n        print(f\"      Processed {min(i + BATCH_SIZE, len(test_proteins)):,}/{len(test_proteins):,} proteins...\")\n    \n    # Clear memory regularly\n    if (i // BATCH_SIZE) % 5 == 0:\n        gc.collect()\n\nprint(f\"   ✓ Generated {len(submission_list):,} predictions\")\n\n# Create final submission\nif submission_list:\n    submission_df = pd.DataFrame(submission_list, columns=['protein', 'term', 'score'])\n    submission_df = submission_df.sort_values(['protein', 'score'], ascending=[True, False])\n    submission_df = submission_df.groupby('protein').head(MAX_PREDS_PER_PROTEIN).reset_index(drop=True)\n\n    # Format scores\n    def safe_format_score(x):\n        try:\n            score = max(0.001, min(1.0, float(x)))\n            return f\"{score:.3f}\"\n        except:\n            return \"0.001\"\n\n    submission_df['score'] = submission_df['score'].apply(safe_format_score)\n\n    # Save submission\n    submission_df.to_csv('submission.tsv', sep='\\t', index=False, header=False)\n\n    print(f\"\\n   ✓ Final submission saved: {len(submission_df):,} predictions\")\n    print(f\"   ✓ Unique proteins: {submission_df['protein'].nunique():,}\")\n    print(f\"   ✓ Mean predictions per protein: {len(submission_df)/submission_df['protein'].nunique():.1f}\")\nelse:\n    print(\"\\n   ⚠️ No predictions generated - creating minimal submission\")\n    # Create minimal submission file to avoid errors\n    minimal_submission = pd.DataFrame({\n        'protein': ['A0A000'] if test_proteins else ['A0A000'],\n        'term': ['GO:0003674'],\n        'score': ['0.001']\n    })\n    minimal_submission.to_csv('submission.tsv', sep='\\t', index=False, header=False)\n    print(\"   ✓ Created minimal submission file\")\n\nprint(\"\\n\" + \"=\"*80)\nprint(\"✅ OPTIMIZED PIPELINE COMPLETED SUCCESSFULLY!\")\nprint(\"=\"*80)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-23T13:57:06.579347Z","iopub.execute_input":"2025-10-23T13:57:06.579679Z","iopub.status.idle":"2025-10-23T13:57:30.648505Z","shell.execute_reply.started":"2025-10-23T13:57:06.579657Z","shell.execute_reply":"2025-10-23T13:57:30.647698Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# # ============================================================================\n# # CAFA 6 PROTEIN FUNCTION PREDICTION - OPTIMIZED PIPELINE\n# # ============================================================================\n\n# import pandas as pd\n# import numpy as np\n# from pathlib import Path\n# from collections import Counter\n# import warnings\n# warnings.filterwarnings('ignore')\n\n# # BioPython for sequence parsing\n# from Bio import SeqIO\n\n# # Modeling libraries\n# from sklearn.preprocessing import MultiLabelBinarizer, StandardScaler\n# from sklearn.multiclass import OneVsRestClassifier\n# from sklearn.linear_model import LogisticRegression\n# from sklearn.feature_selection import SelectKBest, f_classif\n\n# # Memory optimization\n# import gc\n\n# print(\"=\"*80)\n# print(\"CAFA 6 PROTEIN FUNCTION PREDICTION - OPTIMIZED PIPELINE\")\n# print(\"=\"*80)\n\n# # ============================================================================\n# # 1. CONFIGURATION & PATHS\n# # ============================================================================\n# BASE = Path('/kaggle/input/cafa-6-protein-function-prediction')\n# TRAIN_DIR = BASE / 'Train'\n# TEST_DIR = BASE / 'Test'\n\n# # Optimized parameters\n# BATCH_SIZE = 1000  # Reduced for memory efficiency\n# PREDICTION_THRESHOLD = 0.015  # Slightly higher for better precision\n# MAX_PREDS_PER_PROTEIN = 1200\n# MAX_PREDS_PER_ONT = 400\n\n# # ============================================================================\n# # 2. LOAD GO ONTOLOGY\n# # ============================================================================\n# print(\"\\n[1/9] Loading GO ontology...\")\n# go_graph = obonet.read_obo(TRAIN_DIR / 'go-basic.obo')\n# print(f\"   ✓ Loaded {len(go_graph)} GO terms\")\n\n# # ============================================================================\n# # 3. OPTIMIZED DATA LOADING\n# # ============================================================================\n# print(\"\\n[2/9] Loading training data with optimization...\")\n\n# # Load only necessary columns\n# train_terms = pd.read_csv(TRAIN_DIR / 'train_terms.tsv', sep='\\t', \n#                           names=['protein', 'term', 'ontology'],\n#                           usecols=[0, 1, 2])  # Specify columns to save memory\n\n# # Sample proteins for faster training (adjustable)\n# SAMPLE_FRACTION = 0.7  # Use 70% of data for faster processing\n# unique_proteins = train_terms['protein'].unique()\n# sampled_proteins = np.random.choice(unique_proteins, \n#                                    size=int(len(unique_proteins) * SAMPLE_FRACTION), \n#                                    replace=False)\n\n# train_terms = train_terms[train_terms['protein'].isin(sampled_proteins)]\n\n# print(f\"   ✓ Loaded {len(train_terms)} annotations for {train_terms['protein'].nunique()} proteins\")\n\n# # ============================================================================\n# # 4. ENHANCED FEATURE EXTRACTION\n# # ============================================================================\n# print(\"\\n[3/9] Loading training sequences with optimized features...\")\n\n# def optimized_feature_extraction(seq):\n#     \"\"\"Enhanced feature extraction with better biological insights\"\"\"\n#     FEATURE_SIZE = 30  # Increased for better representation\n    \n#     if not seq or len(seq) == 0:\n#         return np.zeros(FEATURE_SIZE)\n    \n#     try:\n#         length = len(seq)\n#         aa_counts = Counter(seq)\n#         total_aa = sum(aa_counts.values())\n        \n#         # 1. Amino Acid Composition (20 features)\n#         aa_list = 'ACDEFGHIKLMNPQRSTVWY'\n#         aa_freq = np.array([aa_counts.get(aa, 0) / total_aa for aa in aa_list])\n        \n#         # 2. Physicochemical Properties (8 features)\n#         # Hydrophobicity groups\n#         very_hydrophobic = sum(aa_counts.get(aa, 0) for aa in 'AILMFWYV') / total_aa\n#         hydrophobic = sum(aa_counts.get(aa, 0) for aa in 'CGT') / total_aa\n        \n#         # Charge properties\n#         positive_charged = sum(aa_counts.get(aa, 0) for aa in 'KRH') / total_aa\n#         negative_charged = sum(aa_counts.get(aa, 0) for aa in 'DE') / total_aa\n        \n#         # Structural properties\n#         polar = sum(aa_counts.get(aa, 0) for aa in 'STNQ') / total_aa\n#         aromatic = sum(aa_counts.get(aa, 0) for aa in 'FWY') / total_aa\n#         small = sum(aa_counts.get(aa, 0) for aa in 'AG') / total_aa\n#         proline_content = aa_counts.get('P', 0) / total_aa\n        \n#         # 3. Sequence Properties (2 features)\n#         seq_complexity = len(set(seq)) / length if length > 0 else 0\n#         molecular_weight_approx = length * 110  # Average AA molecular weight\n        \n#         # Combine all features\n#         features = np.concatenate([\n#             aa_freq,  # 20 features\n#             [\n#                 np.log1p(length), very_hydrophobic, hydrophobic,\n#                 positive_charged, negative_charged, polar, aromatic,\n#                 small, proline_content, seq_complexity,\n#                 np.log1p(molecular_weight_approx)\n#             ]  # 11 features\n#         ])\n        \n#         # Ensure correct size and handle NaN\n#         features = np.nan_to_num(features, nan=0.0, posinf=0.0, neginf=0.0)\n        \n#         if len(features) != FEATURE_SIZE:\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. MEMORY-EFFICIENT FEATURE PROCESSING\n# # ============================================================================\n# print(\"\\n[4/9] Processing training features in batches...\")\n\n# def process_proteins_batch(protein_seqs, batch_size=5000):\n#     \"\"\"Process proteins in batches to save memory\"\"\"\n#     all_features = []\n#     all_proteins = []\n    \n#     proteins = list(protein_seqs.items())\n    \n#     for i in range(0, len(proteins), batch_size):\n#         batch = proteins[i:i + batch_size]\n#         batch_features = []\n#         batch_proteins = []\n        \n#         for pid, seq in batch:\n#             features = optimized_feature_extraction(seq)\n#             batch_features.append(features)\n#             batch_proteins.append(pid)\n        \n#         all_features.extend(batch_features)\n#         all_proteins.extend(batch_proteins)\n        \n#         # Clear memory\n#         del batch, batch_features\n#         gc.collect()\n        \n#         print(f\"      Processed {min(i + batch_size, len(proteins)):,}/{len(proteins):,} proteins...\")\n    \n#     return np.array(all_features), all_proteins\n\n# # Load and process training sequences\n# train_seqs = {}\n# train_proteins_processed = 0\n\n# for rec in SeqIO.parse(TRAIN_DIR / 'train_sequences.fasta', 'fasta'):\n#     pid = rec.id.split('|')[1] if '|' in rec.id else rec.id\n#     # Only process sampled proteins\n#     if pid in sampled_proteins:\n#         train_seqs[pid] = str(rec.seq)\n#         train_proteins_processed += 1\n        \n#     if train_proteins_processed >= len(sampled_proteins):\n#         break\n\n# print(f\"   ✓ Loaded {len(train_seqs)} training sequences\")\n\n# # Process in batches\n# X_train, y_train_proteins = process_proteins_batch(train_seqs, batch_size=5000)\n# print(f\"   ✓ Feature matrix shape: {X_train.shape}\")\n\n# # Clean memory\n# del train_seqs\n# gc.collect()\n\n# # Standardize features\n# print(\"   Standardizing features...\")\n# scaler = StandardScaler()\n# X_train_scaled = scaler.fit_transform(X_train)\n\n# # ============================================================================\n# # 6. FEATURE SELECTION FOR EFFICIENCY\n# # ============================================================================\n# print(\"\\n[5/9] Applying feature selection...\")\n# # Select top 20 features to reduce dimensionality\n# feature_selector = SelectKBest(f_classif, k=20)\n# X_train_selected = feature_selector.fit_transform(X_train_scaled, np.ones(len(X_train_scaled)))\n\n# print(f\"   ✓ Reduced features from {X_train_scaled.shape[1]} to {X_train_selected.shape[1]}\")\n\n# # Clean memory\n# del X_train, X_train_scaled\n# gc.collect()\n\n# # ============================================================================\n# # 7. PREPARE LABELS WITH MEMORY OPTIMIZATION\n# # ============================================================================\n# print(\"\\n[6/9] Preparing labels by ontology...\")\n\n# ontologies = {'F': 'MFO', 'P': 'BPO', 'C': 'CCO'}\n# mlb_dict = {}\n# y_train_dict = {}\n\n# for 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    \n#     # Only include proteins we have features for\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# # Clean memory\n# del train_terms\n# gc.collect()\n\n# # ============================================================================\n# # 8. EFFICIENT MODEL TRAINING\n# # ============================================================================\n# print(\"\\n[7/9] Training optimized models...\")\n\n# models = {}\n\n# for 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#     # Memory-efficient configuration\n#     model = OneVsRestClassifier(\n#         LogisticRegression(\n#             max_iter=500,  # Reduced iterations\n#             solver='liblinear',  # More memory efficient\n#             C=0.8,  # Slightly more regularization\n#             random_state=42,\n#             class_weight='balanced'\n#         ),\n#         n_jobs=1  # Single job to avoid memory issues\n#     )\n    \n#     model.fit(X_train_selected, y_ont)\n#     models[ont_code] = model\n    \n#     print(f\"      ✓ {ont_name} model trained\")\n    \n#     # Clean memory after each model\n#     del y_ont\n#     gc.collect()\n\n# print(f\"   ✓ All {len(models)} models trained successfully\")\n\n# # ============================================================================\n# # 9. OPTIMIZED PREDICTION PIPELINE\n# # ============================================================================\n# print(\"\\n[8/9] Loading and processing test sequences...\")\n\n# test_seqs = {}\n# test_proteins = []\n\n# # Process test sequences in streaming fashion\n# for 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\n# print(f\"   ✓ Loaded {len(test_seqs):,} test sequences\")\n\n# print(\"\\n[9/9] Making optimized predictions...\")\n\n# submission_list = []\n\n# def process_prediction_batch(protein_batch, seq_dict):\n#     \"\"\"Process predictions for a batch of proteins\"\"\"\n#     batch_features = []\n#     valid_proteins = []\n    \n#     for pid in protein_batch:\n#         features = optimized_feature_extraction(seq_dict[pid])\n#         if len(features) == X_train_selected.shape[1]:\n#             batch_features.append(features)\n#             valid_proteins.append(pid)\n    \n#     if not batch_features:\n#         return []\n    \n#     X_batch = np.array(batch_features)\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#         X_batch_selected = feature_selector.transform(X_batch_scaled)\n#     except Exception as e:\n#         print(f\"   Error processing batch: {e}\")\n#         return []\n    \n#     batch_predictions = []\n    \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_selected)\n            \n#             for i, pid in enumerate(valid_proteins):\n#                 probs = y_pred_proba[i]\n#                 top_indices = np.where(probs > PREDICTION_THRESHOLD)[0]\n                \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#                         batch_predictions.append((pid, term, score))\n                        \n#         except Exception as e:\n#             print(f\"   Error in {ont_code} prediction: {e}\")\n#             continue\n    \n#     return batch_predictions\n\n# # Process test proteins in small batches\n# for i in range(0, len(test_proteins), BATCH_SIZE):\n#     batch = test_proteins[i:i + BATCH_SIZE]\n#     batch_predictions = process_prediction_batch(batch, test_seqs)\n#     submission_list.extend(batch_predictions)\n    \n#     if (i + BATCH_SIZE) % 5000 == 0 or (i + BATCH_SIZE) >= len(test_proteins):\n#         print(f\"      Processed {min(i + BATCH_SIZE, len(test_proteins)):,}/{len(test_proteins):,} proteins...\")\n    \n#     # Clear memory regularly\n#     if (i // BATCH_SIZE) % 10 == 0:\n#         gc.collect()\n\n# print(f\"   ✓ Generated {len(submission_list):,} predictions\")\n\n# # Create final submission\n# submission_df = pd.DataFrame(submission_list, columns=['protein', 'term', 'score'])\n# submission_df = submission_df.sort_values(['protein', 'score'], ascending=[True, False])\n# submission_df = submission_df.groupby('protein').head(MAX_PREDS_PER_PROTEIN).reset_index(drop=True)\n\n# # Format scores\n# def safe_format_score(x):\n#     try:\n#         score = max(0.001, min(1.0, float(x)))\n#         return f\"{score:.3f}\"\n#     except:\n#         return \"0.001\"\n\n# submission_df['score'] = submission_df['score'].apply(safe_format_score)\n\n# # Save submission\n# submission_df.to_csv('submission.tsv', sep='\\t', index=False, header=False)\n\n# print(f\"\\n   ✓ Final submission saved: {len(submission_df):,} predictions\")\n# print(f\"   ✓ Unique proteins: {submission_df['protein'].nunique():,}\")\n# print(f\"   ✓ Mean predictions per protein: {len(submission_df)/submission_df['protein'].nunique():.1f}\")\n\n# print(\"\\n\" + \"=\"*80)\n# print(\"✅ OPTIMIZED PIPELINE COMPLETED SUCCESSFULLY!\")\n# print(\"=\"*80)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-23T13:56:03.050847Z","iopub.execute_input":"2025-10-23T13:56:03.051186Z","iopub.status.idle":"2025-10-23T13:56:03.064369Z","shell.execute_reply.started":"2025-10-23T13:56:03.051163Z","shell.execute_reply":"2025-10-23T13:56:03.063358Z"}},"outputs":[],"execution_count":null}]}