{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.11.0"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":116062,"databundleVersionId":14875579,"sourceType":"competition"}],"dockerImageVersionId":31234,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Protein GO Function Prediction: EDA, Visual and LR\n**Task:** Multi-label classification of Gene Ontology (GO) terms from protein sequences\n\n### Table of Contents\n1. **Section 1: Setup and Configuration**\n2. **Section 2: Data Loading and Preprocessing**\n3. **Section 3: Exploratory Data Analysis (EDA)**\n4. **Section 4: Feature Engineering**\n5. **Section 5: Model Training**\n6. **Section 6: Model Evaluation and Visualization**\n7. **Section 7: Inference and Submission**","metadata":{}},{"cell_type":"markdown","source":"---\n# Section 1: Setup and Configuration\n---","metadata":{}},{"cell_type":"code","source":"# 1.1 Install Dependencies\n!pip install biopython -q\n!pip install obonet -q","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 1.2 Import Libraries\nimport pandas as pd\nimport numpy as np\nfrom Bio import SeqIO\nfrom tqdm.auto import tqdm\nfrom collections import Counter\n\n# Visualization\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom matplotlib.gridspec import GridSpec\nfrom matplotlib.colors import LinearSegmentedColormap\n\n# Statistics\nfrom scipy.stats import pearsonr, spearmanr, linregress\n\n# Modeling\nfrom sklearn.preprocessing import MultiLabelBinarizer\nfrom sklearn.multiclass import OneVsRestClassifier\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import precision_recall_curve, average_precision_score\n\nimport warnings\nwarnings.filterwarnings('ignore')\npd.set_option('display.max_columns', None)\n\nprint(\"All libraries imported successfully!\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 1.3 Configuration and CNS Color Palette\n\nclass CFG:\n    \"\"\"Global configuration for file paths and model parameters.\"\"\"\n    # File Paths (Kaggle input directory structure for CAFA 6)\n    DATA_DIR = '/kaggle/input/cafa-6-protein-function-prediction'\n    TRAIN_TERMS_PATH = f'{DATA_DIR}/Train/train_terms.tsv'\n    TRAIN_SEQUENCES_PATH = f'{DATA_DIR}/Train/train_sequences.fasta'\n    TRAIN_TAXONOMY_PATH = f'{DATA_DIR}/Train/train_taxonomy.tsv'\n    TEST_SEQUENCES_PATH = f'{DATA_DIR}/Test/testsuperset.fasta'\n    IA_PATH = f'{DATA_DIR}/IA.tsv'  # Fixed: IA.tsv not IA.txt\n    \n    # Modeling Parameters\n    PROBABILITY_THRESHOLD = 0.02\n    RANDOM_STATE = 42\n    TEST_SIZE = 0.2\n\n# CNS-Level Color Palettes (Nature/Science/Cell Style)\nCNS_COLORS = {\n    'primary': '#2C3E50',\n    'secondary': '#E74C3C',\n    'tertiary': '#3498DB',\n    'quaternary': '#27AE60',\n    'accent': '#9B59B6',\n    'neutral': '#95A5A6',\n    'background': '#FAFBFC',\n    'text': '#2C3E50',\n}\n\nNATURE_PALETTE = ['#E64B35', '#4DBBD5', '#00A087', '#3C5488', '#F39B7F', '#8491B4', '#91D1C2', '#DC0000', '#7E6148', '#B09C85']\nSCIENCE_PALETTE = ['#3B4992', '#EE0000', '#008B45', '#631879', '#008280', '#BB0021', '#5F559B', '#A20056', '#808180', '#1B1919']\nCELL_PALETTE = ['#BC3C29', '#0072B5', '#E18727', '#20854E', '#7876B1', '#6F99AD', '#FFDC91', '#EE4C97']\n\n# Taxonomy ID to Species Name Mapping\nTAXON_NAMES = {\n    9606: 'H. sapiens', 10090: 'M. musculus', 7955: 'D. rerio',\n    3702: 'A. thaliana', 559292: 'S. cerevisiae', 7227: 'D. melanogaster',\n    36329: 'P. falciparum', 83333: 'E. coli', 10116: 'R. norvegicus', 9913: 'B. taurus'\n}\n\nprint(\"Configuration loaded successfully!\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 1.4 Plotting Utility Functions\n\ndef style_axis(ax, title=None, xlabel=None, ylabel=None, hide_spines=['top', 'right']):\n    \"\"\"Apply consistent CNS styling to an axis.\"\"\"\n    ax.set_facecolor(CNS_COLORS['background'])\n    for spine in hide_spines:\n        ax.spines[spine].set_visible(False)\n    for spine in ['left', 'bottom']:\n        if spine not in hide_spines:\n            ax.spines[spine].set_linewidth(1.2)\n            ax.spines[spine].set_color(CNS_COLORS['text'])\n    ax.tick_params(labelsize=10, colors=CNS_COLORS['text'])\n    if title:\n        ax.set_title(title, fontsize=13, weight='bold', fontfamily='serif', color=CNS_COLORS['text'], pad=10)\n    if xlabel:\n        ax.set_xlabel(xlabel, fontsize=11, weight='bold', fontfamily='serif')\n    if ylabel:\n        ax.set_ylabel(ylabel, fontsize=11, weight='bold', fontfamily='serif')\n\ndef add_figure_labels(fig, title, subtitle, signature='Created by yujuan.'):\n    \"\"\"Add title, subtitle, and signature to a figure.\"\"\"\n    fig.text(0.01, 1.0, title, {'font': 'serif', 'size': 24, 'weight': 'bold', 'color': CNS_COLORS['text']}, alpha=0.95)\n    fig.text(0.01, 0.95, subtitle, {'font': 'serif', 'size': 12, 'color': '#5D6D7E'}, alpha=0.9)\n    fig.text(0.99, 0.01, signature, {'font': 'serif', 'size': 10, 'weight': 'bold', 'ha': 'right', 'color': CNS_COLORS['text']}, alpha=0.8)\n\ndef create_figure(rows=1, cols=1, figsize=(20, 10), dpi=100, **kwargs):\n    \"\"\"Create a figure with CNS styling.\"\"\"\n    fig = plt.figure(figsize=figsize, dpi=dpi)\n    fig.patch.set_facecolor(CNS_COLORS['background'])\n    gs = GridSpec(rows, cols, figure=fig, **kwargs)\n    return fig, gs\n\nprint(\"Plotting utilities loaded!\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n# Section 2: Data Loading and Preprocessing\n---","metadata":{}},{"cell_type":"code","source":"# 2.1 Load Tabular Datasets\nprint(\"Loading tabular datasets...\")\ntrain_terms_df = pd.read_csv(CFG.TRAIN_TERMS_PATH, sep='\\t')\ntrain_taxonomy_df = pd.read_csv(CFG.TRAIN_TAXONOMY_PATH, sep='\\t', header=None, names=['EntryID', 'taxonomyID'])\nia_df = pd.read_csv(CFG.IA_PATH, sep='\\t', names=['term', 'ia'])\n\nprint(f\"  - train_terms: {train_terms_df.shape}\")\nprint(f\"  - train_taxonomy: {train_taxonomy_df.shape}\")\nprint(f\"  - ia_df: {ia_df.shape}\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 2.2 Load FASTA Sequences\ndef load_fasta_to_dataframe(file_path, is_train=True):\n    \"\"\"Parse FASTA file and return a DataFrame with EntryID and sequence.\"\"\"\n    records = []\n    for record in tqdm(SeqIO.parse(file_path, \"fasta\"), desc=f\"Parsing {file_path.split('/')[-1]}\"):\n        entry_id = record.id.split('|')[1] if is_train and '|' in record.id else record.id.split()[0]\n        records.append({'EntryID': entry_id, 'sequence': str(record.seq)})\n    return pd.DataFrame(records)\n\nprint(\"\\nLoading sequence datasets...\")\ntrain_sequences_df = load_fasta_to_dataframe(CFG.TRAIN_SEQUENCES_PATH, is_train=True)\nprint(f\"  - train_sequences: {train_sequences_df.shape}\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 2.3 Consolidate Data into Master DataFrame\nprint(\"\\nConsolidating data...\")\nprotein_labels = train_terms_df.groupby('EntryID')['term'].apply(list).reset_index(name='labels')\ntrain_df = pd.merge(train_sequences_df, train_taxonomy_df, on='EntryID', how='left')\ntrain_df = pd.merge(train_df, protein_labels, on='EntryID', how='inner')\ntrain_df['seq_length'] = train_df['sequence'].str.len()\ntrain_df['num_labels'] = train_df['labels'].str.len()\ntrain_df['species_name'] = train_df['taxonomyID'].map(TAXON_NAMES).fillna('Other')\n\nprint(f\"  - Final train_df: {train_df.shape}\")\nprint(\"\\nData loading complete!\")\ntrain_df.head(3)","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n# Section 3: Exploratory Data Analysis (EDA)\n---","metadata":{}},{"cell_type":"code","source":"# 3.1 Dataset Overview Dashboard\nfig, gs = create_figure(2, 3, figsize=(20, 12), hspace=0.35, wspace=0.3)\n\n# Panel A: Dataset Statistics\nax1 = fig.add_subplot(gs[0, 0])\nstats_text = f\"Total Proteins\\n{train_df['EntryID'].nunique():,}\\n\\nTotal Annotations\\n{len(train_terms_df):,}\\n\\nUnique GO Terms\\n{train_terms_df['term'].nunique():,}\"\nax1.text(0.5, 0.5, stats_text, ha='center', va='center', fontsize=18, fontweight='bold', color=CNS_COLORS['text'], fontfamily='serif', linespacing=1.8)\nax1.set_facecolor('#F0F3F5')\nfor spine in ax1.spines.values(): spine.set_visible(False)\nax1.set_xticks([]); ax1.set_yticks([])\nax1.set_title('Dataset Statistics', fontsize=13, weight='bold', fontfamily='serif', pad=15, color=CNS_COLORS['text'])\n\n# Panel B: Ontology Distribution\nax2 = fig.add_subplot(gs[0, 1])\nont_counts = train_terms_df['aspect'].value_counts()\nont_labels = {'P': 'Biological Process', 'C': 'Cellular Component', 'F': 'Molecular Function'}\nwedges, texts, autotexts = ax2.pie(ont_counts.values, labels=[ont_labels.get(x, x) for x in ont_counts.index],\n                                    autopct='%1.1f%%', startangle=90, colors=NATURE_PALETTE[:3],\n                                    wedgeprops=dict(width=0.5, edgecolor='white', linewidth=2.5), pctdistance=0.75)\nfor at in autotexts: at.set(fontsize=10, weight='bold', color='white', fontfamily='serif')\nstyle_axis(ax2, title='Ontology Distribution')\n\n# Panel C: Top Species\nax3 = fig.add_subplot(gs[0, 2])\ntop_taxons = train_df['species_name'].value_counts().nlargest(8)\nax3.barh(range(len(top_taxons)), top_taxons.values, color=[CELL_PALETTE[i % len(CELL_PALETTE)] for i in range(len(top_taxons))], height=0.7, edgecolor='white', linewidth=1.5)\nax3.set_yticks(range(len(top_taxons)))\nax3.set_yticklabels([f'$\\it{{{name}}}$' if name != 'Other' else name for name in top_taxons.index], fontsize=9, fontfamily='serif')\nfor i, val in enumerate(top_taxons.values): ax3.text(val + max(top_taxons.values)*0.02, i, f'{val:,}', va='center', fontsize=9, weight='bold', fontfamily='serif')\nax3.invert_yaxis(); ax3.set_xlim(0, max(top_taxons.values) * 1.15)\nstyle_axis(ax3, title='Top 8 Species', xlabel='Protein Count')\n\n# Panel D: Sequence Length Distribution\nax4 = fig.add_subplot(gs[1, :])\nsns.histplot(train_df['seq_length'], ax=ax4, color=SCIENCE_PALETTE[0], bins=50, kde=True, alpha=0.7, edgecolor='white', linewidth=0.8, line_kws={'linewidth': 2.5, 'color': SCIENCE_PALETTE[1]})\nmedian_len = train_df['seq_length'].median()\nax4.axvline(median_len, color=NATURE_PALETTE[0], linestyle='--', linewidth=2.5, label=f'Median: {median_len:.0f} aa')\nax4.set_xlim(0, 2500)\nax4.legend(fontsize=10, frameon=True, fancybox=True, framealpha=0.9, loc='upper right')\nstyle_axis(ax4, title='Sequence Length Distribution', xlabel='Protein Sequence Length (amino acids)', ylabel='Frequency')\n\nadd_figure_labels(fig, 'CAFA 6 Dataset Overview: Scale, Structure, and Composition',\n                  'Biological Process dominates annotations. Human proteins are most abundant. Most proteins are under 1000 aa.')\nplt.tight_layout(rect=[0, 0.02, 1, 0.92])\nplt.show()","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 3.2 Amino Acid Composition Analysis\nfig, gs = create_figure(2, 3, figsize=(20, 12), hspace=0.4, wspace=0.3)\n\n# Count amino acids\nall_sequences = \"\".join(train_sequences_df['sequence'])\naa_counts = Counter(all_sequences)\naa_df = pd.DataFrame.from_dict(aa_counts, orient='index', columns=['count']).sort_values('count', ascending=False)\naa_df = aa_df[aa_df.index.isin(list('ACDEFGHIKLMNPQRSTVWY'))]\n\n# Panel A: Global Distribution\nax1 = fig.add_subplot(gs[0, :])\nn_aa = len(aa_df)\ngradient_colors = [plt.cm.RdYlBu_r(0.2 + 0.6 * i / n_aa) for i in range(n_aa)]\nax1.bar(aa_df.index, aa_df['count'], color=gradient_colors, width=0.75, edgecolor='white', linewidth=1.2, zorder=3)\nax1.plot(range(len(aa_df)), aa_df['count'], color=CNS_COLORS['text'], linewidth=2.5, marker='o', markersize=6, markerfacecolor='white', markeredgecolor=CNS_COLORS['text'], markeredgewidth=1.5, zorder=4)\nax1.grid(axis='y', linestyle='--', alpha=0.3, zorder=0)\nstyle_axis(ax1, title='Global Amino Acid Frequency Distribution', ylabel='Total Count')\n\n# Panels B-D: Grouped by Chemical Properties\naa_groups = {'Hydrophobic': ('AILMFWYV', CELL_PALETTE[0]), 'Polar / Neutral': ('STNQCGP', CELL_PALETTE[1]), 'Charged (+/-)': ('RHKDE', CELL_PALETTE[3])}\nfor ax, (group_name, (aas, color)) in zip([fig.add_subplot(gs[1, i]) for i in range(3)], aa_groups.items()):\n    group_counts = {aa: aa_counts.get(aa, 0) for aa in aas}\n    group_df = pd.DataFrame.from_dict(group_counts, orient='index', columns=['count']).sort_values('count', ascending=False)\n    ax.bar(group_df.index, group_df.values.flatten(), color=color, width=0.7, edgecolor='white', linewidth=1.5, alpha=0.85)\n    for i, (idx, row) in enumerate(group_df.iterrows()): ax.text(i, row['count'] + max(group_df['count'])*0.02, f'{row[\"count\"]:,}', ha='center', fontsize=9, weight='bold', fontfamily='serif')\n    style_axis(ax, title=group_name, ylabel='Count')\n    ax.set_title(group_name, fontsize=13, weight='bold', fontfamily='serif', pad=10, color=color)\n\nadd_figure_labels(fig, 'The Building Blocks of Life: Amino Acid Composition Analysis',\n                  'Leucine (L) and Alanine (A) are most abundant. Hydrophobic residues dominate, essential for protein folding.')\nplt.tight_layout(rect=[0, 0.02, 1, 0.92])\nplt.show()","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 3.3 Statistical Correlations Analysis\nfig, gs = create_figure(2, 3, figsize=(20, 12), hspace=0.35, wspace=0.3)\n\n# Panel A: Length vs Labels Regression\nax1 = fig.add_subplot(gs[0, 0])\nr_pearson, _ = pearsonr(train_df['seq_length'], train_df['num_labels'])\nr_spearman, _ = spearmanr(train_df['seq_length'], train_df['num_labels'])\nsample_data = train_df.sample(min(1000, len(train_df)), random_state=CFG.RANDOM_STATE)\nsns.regplot(x='seq_length', y='num_labels', data=sample_data, ax=ax1,\n            scatter_kws={'alpha': 0.4, 'color': NATURE_PALETTE[3], 's': 35, 'edgecolor': 'white'},\n            line_kws={'color': NATURE_PALETTE[0], 'linewidth': 2.5})\nax1.text(0.05, 0.95, f'Pearson r = {r_pearson:.3f}\\nSpearman ρ = {r_spearman:.3f}', transform=ax1.transAxes, fontsize=10, weight='bold', va='top', fontfamily='serif',\n         bbox=dict(boxstyle='round,pad=0.5', facecolor='white', edgecolor=NATURE_PALETTE[3], alpha=0.9, linewidth=1.5))\nax1.set_xlim(0, 3000); ax1.set_ylim(0, 100)\nstyle_axis(ax1, title='Length vs. Functional Annotations', xlabel='Sequence Length', ylabel='Number of Labels')\n\n# Panel B: Boxplot by Ontology\nax2 = fig.add_subplot(gs[0, 1])\nont_label_counts = train_terms_df.groupby(['EntryID', 'aspect']).size().reset_index(name='count')\nbp = ax2.boxplot([ont_label_counts[ont_label_counts['aspect'] == asp]['count'].values for asp in ['P', 'C', 'F']],\n                  labels=['BPO', 'CCO', 'MFO'], patch_artist=True, widths=0.6,\n                  medianprops=dict(color='black', linewidth=2), whiskerprops=dict(color=CNS_COLORS['text'], linewidth=1.5),\n                  capprops=dict(color=CNS_COLORS['text'], linewidth=1.5), flierprops=dict(marker='o', markerfacecolor=CNS_COLORS['neutral'], markersize=4, alpha=0.5))\nfor patch, color in zip(bp['boxes'], [CELL_PALETTE[0], CELL_PALETTE[1], CELL_PALETTE[3]]): patch.set(facecolor=color, alpha=0.8, edgecolor='black', linewidth=1.5)\nstyle_axis(ax2, title='Label Distribution by Ontology', xlabel='Ontology Aspect', ylabel='Labels per Protein')\n\n# Panel C: IA Distribution\nax3 = fig.add_subplot(gs[0, 2])\nsns.histplot(ia_df['ia'], ax=ax3, color=SCIENCE_PALETTE[0], bins=40, kde=True, alpha=0.75, edgecolor='white', linewidth=0.8, line_kws={'linewidth': 2.5, 'color': SCIENCE_PALETTE[1]})\nmedian_ia = ia_df['ia'].median()\nax3.axvline(median_ia, color=SCIENCE_PALETTE[1], linestyle='--', linewidth=2.5, label=f\"Median: {median_ia:.2f}\")\nax3.legend(fontsize=10, frameon=True)\nstyle_axis(ax3, title='IA Score Distribution', xlabel='Information Accretion (IA)', ylabel='Frequency')\n\n# Panel D: Correlation Heatmap\nax4 = fig.add_subplot(gs[1, :])\ncorr_data = train_df[['seq_length', 'num_labels']].copy()\ncorr_data['avg_ia'] = train_df['labels'].apply(lambda x: ia_df[ia_df['term'].isin(x)]['ia'].mean() if len(x) > 0 else 0)\ncorr_matrix = corr_data.corr()\ncns_cmap = LinearSegmentedColormap.from_list('cns', ['#FFFFFF', '#4DBBD5', '#3C5488', '#2C3E50'])\nsns.heatmap(corr_matrix, annot=True, cmap=cns_cmap, ax=ax4, fmt='.3f', annot_kws={'size': 14, 'weight': 'bold', 'fontfamily': 'serif'},\n            linewidths=3, linecolor='white', cbar_kws={'shrink': 0.5, 'label': 'Correlation'},\n            xticklabels=['Sequence Length', 'Num Labels', 'Avg IA Score'], yticklabels=['Sequence Length', 'Num Labels', 'Avg IA Score'], vmin=-1, vmax=1)\nstyle_axis(ax4, title='Feature Correlation Matrix')\n\nadd_figure_labels(fig, 'Statistical Deep Dive: Correlations and Distributions',\n                  'Weak positive correlation (r≈0.2) between length and labels. IA scores show right-skewed distribution.')\nplt.tight_layout(rect=[0, 0.02, 1, 0.92])\nplt.show()","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 3.4 Species Analysis\nfig, gs = create_figure(2, 2, figsize=(20, 12), hspace=0.35, wspace=0.25)\n\n# Panel A: Top 10 Species\nax1 = fig.add_subplot(gs[0, 0])\ntop_species = train_df['species_name'].value_counts().nlargest(10)\nax1.barh(range(len(top_species)), top_species.values, color=[NATURE_PALETTE[i % len(NATURE_PALETTE)] for i in range(len(top_species))], height=0.7, edgecolor='white', linewidth=1.5)\nax1.set_yticks(range(len(top_species)))\nax1.set_yticklabels([f'$\\it{{{name}}}$' if name != 'Other' else name for name in top_species.index], fontsize=10, fontfamily='serif')\nfor i, val in enumerate(top_species.values): ax1.text(val + max(top_species.values)*0.02, i, f'{val:,}', va='center', fontsize=10, weight='bold', fontfamily='serif')\nax1.invert_yaxis(); ax1.set_xlim(0, max(top_species.values) * 1.15)\nstyle_axis(ax1, title='Top 10 Species by Protein Count', xlabel='Number of Proteins')\n\n# Panel B & C: Species vs Length & Labels\ntop_5_species = train_df['species_name'].value_counts().nlargest(5).index\ndf_top5 = train_df[train_df['species_name'].isin(top_5_species)]\n\nax2 = fig.add_subplot(gs[0, 1])\nspecies_avg_len = df_top5.groupby('species_name')['seq_length'].mean().reindex(top_5_species)\nax2.bar(range(len(species_avg_len)), species_avg_len.values, color=[CELL_PALETTE[i] for i in range(len(species_avg_len))], width=0.65, edgecolor='black', linewidth=1.2)\nax2.set_xticks(range(len(species_avg_len)))\nax2.set_xticklabels([f'$\\it{{{name[:12]}}}$' for name in species_avg_len.index], fontsize=9, rotation=20, ha='right', fontfamily='serif')\nfor i, val in enumerate(species_avg_len.values): ax2.text(i, val + 15, f'{val:.0f}', ha='center', fontsize=10, weight='bold', fontfamily='serif')\nstyle_axis(ax2, title='Mean Protein Length by Species', ylabel='Average Sequence Length (aa)')\n\nax3 = fig.add_subplot(gs[1, 0])\nspecies_avg_labels = df_top5.groupby('species_name')['num_labels'].mean().reindex(top_5_species)\nax3.bar(range(len(species_avg_labels)), species_avg_labels.values, color=[SCIENCE_PALETTE[i] for i in range(len(species_avg_labels))], width=0.65, edgecolor='black', linewidth=1.2)\nax3.set_xticks(range(len(species_avg_labels)))\nax3.set_xticklabels([f'$\\it{{{name[:12]}}}$' for name in species_avg_labels.index], fontsize=9, rotation=20, ha='right', fontfamily='serif')\nfor i, val in enumerate(species_avg_labels.values): ax3.text(i, val + 0.5, f'{val:.1f}', ha='center', fontsize=10, weight='bold', fontfamily='serif')\nstyle_axis(ax3, title='Mean Functional Annotations by Species', ylabel='Average GO Annotations')\n\n# Panel D: Species Pie Chart\nax4 = fig.add_subplot(gs[1, 1])\ntop_5_counts = train_df[train_df['species_name'].isin(top_5_species)].groupby('species_name').size().reindex(top_5_species)\nother_count = len(train_df) - top_5_counts.sum()\npie_data = list(top_5_counts.values) + [other_count]\npie_labels = list(top_5_species) + ['Other Species']\nwedges, texts, autotexts = ax4.pie(pie_data, labels=None, autopct='%1.1f%%', startangle=90, colors=NATURE_PALETTE[:5] + [CNS_COLORS['neutral']],\n                                    wedgeprops=dict(width=0.55, edgecolor='white', linewidth=2.5), pctdistance=0.75)\nfor at in autotexts: at.set(fontsize=10, weight='bold', color='white', fontfamily='serif')\nax4.legend(wedges, [f'$\\it{{{n[:15]}}}$' if n != 'Other Species' else n for n in pie_labels], title='Species', loc='center left', bbox_to_anchor=(1.0, 0.5), fontsize=9, title_fontsize=11, frameon=True)\nstyle_axis(ax4, title='Species Distribution Overview')\n\nadd_figure_labels(fig, 'Species Spotlight: Taxonomic Distribution Analysis',\n                  'Human proteins dominate. Different species show varying protein lengths and annotation depths.')\nplt.tight_layout(rect=[0, 0.02, 1, 0.92])\nplt.show()","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 3.5 GO Term Frequency Analysis\nfig, gs = create_figure(2, 3, figsize=(20, 13), hspace=0.4, wspace=0.3)\n\n# Panel A: Top 15 GO Terms\nax1 = fig.add_subplot(gs[0, :2])\ntop_terms = train_terms_df['term'].value_counts().nlargest(15)\ngradient_colors = [plt.cm.RdYlBu_r(i/len(top_terms)) for i in range(len(top_terms))]\nax1.barh(range(len(top_terms)), top_terms.values, color=gradient_colors, height=0.75, edgecolor='white', linewidth=1)\nax1.set_yticks(range(len(top_terms)))\nax1.set_yticklabels(top_terms.index, fontsize=9, fontfamily='monospace')\nfor i, val in enumerate(top_terms.values): ax1.text(val + max(top_terms.values)*0.01, i, f'{val:,}', va='center', fontsize=9, weight='bold', fontfamily='serif')\nax1.invert_yaxis(); ax1.set_xlim(0, max(top_terms.values) * 1.12)\nstyle_axis(ax1, title='Top 15 Most Frequent GO Terms', xlabel='Protein Count')\n\n# Panel B: Label Count Distribution\nax2 = fig.add_subplot(gs[0, 2])\nlabel_dist = train_df['num_labels'].value_counts().sort_index()\nax2.scatter(label_dist.index, label_dist.values, color=SCIENCE_PALETTE[0], alpha=0.7, s=60, edgecolor='white', linewidth=0.8, zorder=3)\nax2.plot(label_dist.index, label_dist.values, color=SCIENCE_PALETTE[0], alpha=0.3, linewidth=1.5, zorder=2)\nax2.fill_between(label_dist.index, label_dist.values, alpha=0.15, color=SCIENCE_PALETTE[0], zorder=1)\nax2.set_yscale('log')\nax2.grid(True, alpha=0.3, linestyle='--', zorder=0)\nstyle_axis(ax2, title='Label Count Distribution', xlabel='Number of GO Labels', ylabel='Protein Count (log scale)')\n\n# Panels C-E: Top 5 Terms per Ontology\nontology_info = {'P': ('Biological Process (BPO)', CELL_PALETTE[0]), 'C': ('Cellular Component (CCO)', CELL_PALETTE[1]), 'F': ('Molecular Function (MFO)', CELL_PALETTE[3])}\nfor i, (aspect, (full_name, color)) in enumerate(ontology_info.items()):\n    ax = fig.add_subplot(gs[1, i])\n    ont_top = train_terms_df[train_terms_df['aspect'] == aspect]['term'].value_counts().nlargest(5)\n    ax.barh(range(len(ont_top)), ont_top.values, color=color, height=0.65, edgecolor='white', linewidth=1.5, alpha=0.85)\n    ax.set_yticks(range(len(ont_top)))\n    ax.set_yticklabels([t[:18] + '...' if len(t) > 18 else t for t in ont_top.index], fontsize=8, fontfamily='monospace')\n    for j, val in enumerate(ont_top.values): ax.text(val + max(ont_top.values)*0.03, j, f'{val}', va='center', fontsize=9, weight='bold', fontfamily='serif')\n    ax.invert_yaxis(); ax.set_xlim(0, max(ont_top.values) * 1.18)\n    style_axis(ax, xlabel='Count')\n    ax.set_title(full_name, fontsize=11, weight='bold', fontfamily='serif', pad=8, color=color)\n\nadd_figure_labels(fig, 'GO Term Analysis: Functional Annotation Landscape',\n                  \"A few GO terms dominate — generic functions like 'protein binding' appear in thousands of proteins.\")\nplt.tight_layout(rect=[0, 0.02, 1, 0.92])\nplt.show()","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 3.6 Information Content Analysis\nfig, gs = create_figure(1, 2, figsize=(20, 10), wspace=0.25)\n\n# Prepare data\nterm_counts = train_terms_df['term'].value_counts().reset_index()\nterm_counts.columns = ['term', 'protein_count']\nmerged_df = pd.merge(term_counts, ia_df, on='term')\nsample_df = merged_df.sample(n=min(3000, len(merged_df)), random_state=CFG.RANDOM_STATE)\n\n# Panel A: IA vs Frequency\nax1 = fig.add_subplot(gs[0, 0])\nscatter = ax1.scatter(sample_df['protein_count'], sample_df['ia'], c=sample_df['ia'], cmap='RdYlBu_r', alpha=0.6, s=50, edgecolor='white', linewidth=0.5)\nax1.set_xscale('log')\nlog_counts = np.log10(sample_df['protein_count'])\nslope, intercept, r_value, _, _ = linregress(log_counts, sample_df['ia'])\nx_line = np.logspace(0, np.log10(sample_df['protein_count'].max()), 100)\nax1.plot(x_line, slope * np.log10(x_line) + intercept, color=NATURE_PALETTE[0], linewidth=3, linestyle='--', label=f'Trend (r² = {r_value**2:.3f})')\nax1.legend(fontsize=10, frameon=True, loc='upper right')\nax1.grid(True, alpha=0.3, linestyle='--')\ncbar = plt.colorbar(scatter, ax=ax1, shrink=0.8, pad=0.02)\ncbar.set_label('IA Score', fontsize=10, fontfamily='serif')\nstyle_axis(ax1, title='Frequency vs. Information Content', xlabel='Number of Proteins with Term (log scale)', ylabel='Information Accretion (IA) Score')\n\n# Panel B: IA Distribution with Quantiles\nax2 = fig.add_subplot(gs[0, 1])\nsns.histplot(ia_df['ia'], ax=ax2, color=SCIENCE_PALETTE[0], bins=50, kde=True, alpha=0.7, edgecolor='white', linewidth=0.8, line_kws={'linewidth': 2.5, 'color': SCIENCE_PALETTE[1]})\nfor q, color in zip([0.25, 0.5, 0.75], [NATURE_PALETTE[1], NATURE_PALETTE[0], NATURE_PALETTE[2]]):\n    q_val = ia_df['ia'].quantile(q)\n    ax2.axvline(q_val, color=color, linestyle='--', linewidth=2.5, label=f'Q{int(q*100)}: {q_val:.2f}')\nax2.legend(fontsize=10, frameon=True, loc='upper right')\nstyle_axis(ax2, title='IA Score Distribution with Quantiles', xlabel='Information Accretion (IA) Score', ylabel='Frequency')\n\nadd_figure_labels(fig, 'Rarity is Valuable: Information Content Analysis of GO Terms',\n                  \"Strong inverse relationship between frequency and IA. General terms have low IA; rare functions have high value.\")\nplt.tight_layout(rect=[0, 0.02, 1, 0.92])\nplt.show()","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n# Section 4: Feature Engineering\n---","metadata":{}},{"cell_type":"code","source":"# 4.1 Define Feature Extraction Function\ndef extract_sequence_features(seq):\n    \"\"\"Extract 71-dimensional physicochemical features from a protein sequence.\"\"\"\n    if not seq or len(seq) == 0:\n        return np.zeros(71)\n    \n    aa_counts = Counter(seq)\n    length = len(seq)\n    \n    # Amino acid frequencies (20 features)\n    aa_freq = np.array([aa_counts.get(aa, 0) / length for aa in 'ACDEFGHIKLMNPQRSTVWY'])\n    \n    # Physicochemical properties (4 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    aa_weights = {'A': 89, 'C': 121, 'D': 133, 'E': 147, 'F': 165, 'G': 75, 'H': 155, 'I': 131, 'K': 146, 'L': 131, 'M': 149, 'N': 132, 'P': 115, 'Q': 146, 'R': 174, 'S': 105, 'T': 119, 'V': 117, 'W': 204, 'Y': 181}\n    mol_weight = sum(aa_counts.get(aa, 0) * aa_weights.get(aa, 0) for aa in aa_counts)\n    \n    # CTD features (7 features)\n    groups = {'hydrophobic': set('AILMFWYV'), 'polar': set('STNQ'), 'positive': set('RK'), 'negative': set('DE'), 'aromatic': set('FWY'), 'aliphatic': set('ILV'), 'small': set('ACDGNPSTV')}\n    ctd_features = [sum(1 for aa in seq if aa in group_aas) / length for group_aas in groups.values()]\n    \n    # Dipeptide frequencies (20 features)\n    dipeptides = Counter([seq[i:i+2] for i in range(length - 1)])\n    top_dipeptides = ['AL', 'LA', 'LE', 'EA', 'AA', 'AS', 'SA', 'EL', 'LL', 'AE', 'SE', 'ES', 'GA', 'AG', 'VA', 'AV', 'LV', 'VL', 'LS', 'SL']\n    dipeptide_freq = np.array([dipeptides.get(dp, 0) / (length-1 if length > 1 else 1) for dp in top_dipeptides])\n    \n    # Tripeptide frequencies (20 features)\n    tripeptides = Counter([seq[i:i+3] for i in range(length - 2)])\n    top_tripeptides = ['ALA', 'LEA', 'EAL', 'LAL', 'AAA', 'LLE', 'ELE', 'ALE', 'GAL', 'ASA', 'VLA', 'LAV', 'SLS', 'LSL', 'GLA', 'LAG', 'AVL', 'VLA', 'SLE', 'LES']\n    tripeptide_freq = np.array([tripeptides.get(tp, 0) / (length-2 if length > 2 else 1) for tp in top_tripeptides])\n    \n    return np.concatenate([aa_freq, np.array([np.log1p(length), hydrophobic, charged, np.log1p(mol_weight)]), np.array(ctd_features), dipeptide_freq, tripeptide_freq])\n\nprint(\"Feature extraction function defined!\")\nprint(f\"Feature vector dimension: 71\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 4.2 Extract Features for Training Data\nprint(\"Extracting features for all training proteins...\")\n\nX_train_list, y_train_proteins = [], []\nprotein_ids_in_scope = set(train_terms_df['EntryID'].unique())\ntrain_sequences_dict = dict(zip(train_sequences_df['EntryID'], train_sequences_df['sequence']))\n\nfor pid, seq in tqdm(train_sequences_dict.items(), desc=\"Processing Train Sequences\"):\n    if pid in protein_ids_in_scope:\n        X_train_list.append(extract_sequence_features(seq))\n        y_train_proteins.append(pid)\n\nX_train_full = np.array(X_train_list)\nprint(f\"\\nTraining feature matrix shape: {X_train_full.shape}\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n# Section 5: Model Training\n---","metadata":{}},{"cell_type":"code","source":"# 5.1 Train-Validation Split\nprint(\"Splitting data into train and validation sets...\")\n\ntrain_indices, val_indices = train_test_split(\n    range(len(y_train_proteins)), test_size=CFG.TEST_SIZE, random_state=CFG.RANDOM_STATE\n)\n\nX_train = X_train_full[train_indices]\nX_val = X_train_full[val_indices]\ntrain_protein_ids = [y_train_proteins[i] for i in train_indices]\nval_protein_ids = [y_train_proteins[i] for i in val_indices]\n\nprint(f\"  - Training set: {len(train_protein_ids)} proteins\")\nprint(f\"  - Validation set: {len(val_protein_ids)} proteins\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 5.2 Train Ontology-Specific Models\nprint(\"Training models for each ontology...\")\n\nmodels, mlb_dict = {}, {}\nval_predictions = {}  # Store validation predictions for evaluation\nval_labels = {}  # Store validation labels\nontology_map = {'BPO': 'P', 'CCO': 'C', 'MFO': 'F'}\n\nfor aspect_name, aspect_code in ontology_map.items():\n    print(f\"\\n--- Processing {aspect_name} ({aspect_code}) ---\")\n    \n    # Filter terms for current ontology\n    ont_terms_df = train_terms_df[train_terms_df['aspect'] == aspect_code]\n    if ont_terms_df.empty:\n        print(f\"No terms found for {aspect_name}. Skipping.\")\n        continue\n    \n    # Build labels\n    protein_terms = ont_terms_df.groupby('EntryID')['term'].apply(list).to_dict()\n    train_labels_list = [protein_terms.get(pid, []) for pid in train_protein_ids]\n    val_labels_list = [protein_terms.get(pid, []) for pid in val_protein_ids]\n    \n    # Binarize labels\n    mlb = MultiLabelBinarizer(sparse_output=True)\n    y_train_ont = mlb.fit_transform(train_labels_list)\n    y_val_ont = mlb.transform(val_labels_list)\n    \n    print(f\"Number of unique {aspect_name} terms: {y_train_ont.shape[1]}\")\n    \n    if y_train_ont.shape[1] > 0:\n        mlb_dict[aspect_code] = mlb\n        \n        # Train model\n        model = OneVsRestClassifier(\n            LogisticRegression(max_iter=500, solver='lbfgs', C=0.5, random_state=CFG.RANDOM_STATE),\n            n_jobs=-1\n        )\n        model.fit(X_train, y_train_ont)\n        models[aspect_code] = model\n        \n        # Store validation predictions\n        val_predictions[aspect_code] = model.predict_proba(X_val)\n        val_labels[aspect_code] = y_val_ont.toarray()\n        \n        print(f\"Model for {aspect_name} trained successfully.\")\n\nprint(\"\\nAll models trained!\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n# Section 6: Model Evaluation and Visualization\n---","metadata":{}},{"cell_type":"code","source":"# 6.1 Model Performance Dashboard\nfig, gs = create_figure(2, 3, figsize=(20, 12), hspace=0.4, wspace=0.3)\n\nontology_colors = {'P': CELL_PALETTE[0], 'C': CELL_PALETTE[1], 'F': CELL_PALETTE[3]}\nontology_names = {'P': 'Biological Process (BPO)', 'C': 'Cellular Component (CCO)', 'F': 'Molecular Function (MFO)'}\n\n# Panel A-C: Precision-Recall Curves for each ontology\nap_scores = {}\nfor i, (aspect_code, aspect_name) in enumerate(ontology_names.items()):\n    ax = fig.add_subplot(gs[0, i])\n    \n    if aspect_code in val_predictions:\n        y_true = val_labels[aspect_code]\n        y_pred = val_predictions[aspect_code]\n        \n        # Calculate micro-average PR curve\n        precision, recall, _ = precision_recall_curve(y_true.ravel(), y_pred.ravel())\n        ap = average_precision_score(y_true, y_pred, average='micro')\n        ap_scores[aspect_code] = ap\n        \n        ax.plot(recall, precision, color=ontology_colors[aspect_code], linewidth=2.5, label=f'AP = {ap:.3f}')\n        ax.fill_between(recall, precision, alpha=0.2, color=ontology_colors[aspect_code])\n        ax.axhline(y=y_true.sum() / y_true.size, color=CNS_COLORS['neutral'], linestyle='--', linewidth=1.5, label='Baseline')\n        ax.legend(fontsize=10, loc='upper right', frameon=True)\n    \n    style_axis(ax, title=aspect_name, xlabel='Recall', ylabel='Precision')\n    ax.set_xlim(0, 1); ax.set_ylim(0, 1)\n\n# Panel D: Average Precision Comparison\nax4 = fig.add_subplot(gs[1, 0])\nap_values = [ap_scores.get(code, 0) for code in ['P', 'C', 'F']]\nbars = ax4.bar(['BPO', 'CCO', 'MFO'], ap_values, color=[ontology_colors[c] for c in ['P', 'C', 'F']], width=0.6, edgecolor='black', linewidth=1.5)\nfor bar, val in zip(bars, ap_values): ax4.text(bar.get_x() + bar.get_width()/2, val + 0.01, f'{val:.3f}', ha='center', fontsize=12, weight='bold', fontfamily='serif')\nstyle_axis(ax4, title='Average Precision by Ontology', xlabel='Ontology', ylabel='Average Precision')\nax4.set_ylim(0, max(ap_values) * 1.2)\n\n# Panel E: Number of Terms per Ontology\nax5 = fig.add_subplot(gs[1, 1])\nterm_counts = [mlb_dict[code].classes_.shape[0] if code in mlb_dict else 0 for code in ['P', 'C', 'F']]\nbars = ax5.bar(['BPO', 'CCO', 'MFO'], term_counts, color=[ontology_colors[c] for c in ['P', 'C', 'F']], width=0.6, edgecolor='black', linewidth=1.5)\nfor bar, val in zip(bars, term_counts): ax5.text(bar.get_x() + bar.get_width()/2, val + 5, f'{val}', ha='center', fontsize=12, weight='bold', fontfamily='serif')\nstyle_axis(ax5, title='Number of GO Terms per Ontology', xlabel='Ontology', ylabel='Number of Terms')\n\n# Panel F: Prediction Score Distribution\nax6 = fig.add_subplot(gs[1, 2])\nfor aspect_code, color in ontology_colors.items():\n    if aspect_code in val_predictions:\n        scores = val_predictions[aspect_code].ravel()\n        scores_nonzero = scores[scores > 0.001]\n        if len(scores_nonzero) > 0:\n            sns.kdeplot(scores_nonzero, ax=ax6, color=color, linewidth=2.5, label=ontology_names[aspect_code][:3], fill=True, alpha=0.2)\nax6.legend(fontsize=10, frameon=True)\nax6.set_xlim(0, 1)\nstyle_axis(ax6, title='Prediction Score Distribution', xlabel='Predicted Probability', ylabel='Density')\n\nadd_figure_labels(fig, 'Model Evaluation: Performance Metrics Dashboard',\n                  'Precision-Recall curves show model performance. CCO typically achieves highest AP due to fewer, more distinct terms.')\nplt.tight_layout(rect=[0, 0.02, 1, 0.92])\nplt.show()","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 6.2 Top Predicted Terms Analysis\nfig, gs = create_figure(1, 3, figsize=(20, 8), wspace=0.3)\n\nfor i, (aspect_code, aspect_name) in enumerate(ontology_names.items()):\n    ax = fig.add_subplot(gs[0, i])\n    \n    if aspect_code in val_predictions and aspect_code in mlb_dict:\n        # Get average prediction scores for each term\n        mean_scores = val_predictions[aspect_code].mean(axis=0)\n        term_scores = pd.DataFrame({'term': mlb_dict[aspect_code].classes_, 'score': mean_scores})\n        top_10 = term_scores.nlargest(10, 'score')\n        \n        colors = [plt.cm.RdYlBu_r(0.2 + 0.6 * j / 10) for j in range(10)]\n        ax.barh(range(len(top_10)), top_10['score'].values, color=colors, height=0.7, edgecolor='white', linewidth=1)\n        ax.set_yticks(range(len(top_10)))\n        ax.set_yticklabels([t[:20] + '...' if len(t) > 20 else t for t in top_10['term']], fontsize=8, fontfamily='monospace')\n        ax.invert_yaxis()\n        \n        for j, val in enumerate(top_10['score'].values):\n            ax.text(val + 0.005, j, f'{val:.3f}', va='center', fontsize=9, weight='bold', fontfamily='serif')\n    \n    style_axis(ax, xlabel='Mean Prediction Score')\n    ax.set_title(aspect_name, fontsize=12, weight='bold', fontfamily='serif', pad=10, color=ontology_colors[aspect_code])\n\nadd_figure_labels(fig, 'Top Predicted GO Terms by Ontology',\n                  'Most confidently predicted GO terms for each ontology. General terms typically rank highest.')\nplt.tight_layout(rect=[0, 0.02, 1, 0.90])\nplt.show()","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n# Section 7: Inference and Submission\n---","metadata":{}},{"cell_type":"code","source":"# 7.1 Retrain Models on Full Data\nprint(\"Retraining models on full dataset...\")\n\nfinal_models, final_mlb_dict = {}, {}\n\nfor aspect_name, aspect_code in ontology_map.items():\n    print(f\"\\n--- Processing {aspect_name} ({aspect_code}) ---\")\n    \n    ont_terms_df = train_terms_df[train_terms_df['aspect'] == aspect_code]\n    if ont_terms_df.empty:\n        continue\n    \n    protein_terms = ont_terms_df.groupby('EntryID')['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_train_ont = mlb.fit_transform(labels_list)\n    \n    if y_train_ont.shape[1] > 0:\n        final_mlb_dict[aspect_code] = mlb\n        model = OneVsRestClassifier(\n            LogisticRegression(max_iter=500, solver='lbfgs', C=0.5, random_state=CFG.RANDOM_STATE),\n            n_jobs=-1\n        )\n        model.fit(X_train_full, y_train_ont)\n        final_models[aspect_code] = model\n        print(f\"Final model for {aspect_name} trained.\")\n\nprint(\"\\nAll final models trained!\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 7.2 Generate Predictions on Test Set\nprint(\"\\nLoading and processing test sequences...\")\ntest_sequences_df = load_fasta_to_dataframe(CFG.TEST_SEQUENCES_PATH, is_train=False)\ntest_protein_ids = test_sequences_df['EntryID'].tolist()\ntest_sequences_dict = dict(zip(test_sequences_df['EntryID'], test_sequences_df['sequence']))\n\nprint(f\"Test set size: {len(test_protein_ids)} proteins\")\n\n# Batch prediction\nBATCH_SIZE = 5000\nsubmission_list = []\n\nfor i in tqdm(range(0, len(test_protein_ids), BATCH_SIZE), desc=\"Predicting on Test Set\"):\n    batch_pids = test_protein_ids[i:i+BATCH_SIZE]\n    batch_seqs = [test_sequences_dict[pid] for pid in batch_pids]\n    X_batch = np.array([extract_sequence_features(seq) for seq in batch_seqs])\n    \n    for aspect_code, model in final_models.items():\n        mlb = final_mlb_dict[aspect_code]\n        y_pred_proba = model.predict_proba(X_batch)\n        \n        for j, pid in enumerate(batch_pids):\n            probs = y_pred_proba[j]\n            candidate_indices = np.where(probs > CFG.PROBABILITY_THRESHOLD)[0]\n            for idx in candidate_indices:\n                submission_list.append((pid, mlb.classes_[idx], probs[idx]))\n\nprint(f\"\\nGenerated {len(submission_list):,} total predictions.\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 7.3 Create Submission File\nprint(\"Creating submission file...\")\n\nsubmission_df = pd.DataFrame(submission_list, columns=['Protein Id', 'GO Term Id', 'Prediction'])\n\n# Apply 1500 prediction limit per protein\nprint(\"Applying 1500 prediction limit per protein...\")\nsubmission_df = submission_df.sort_values(by=['Protein Id', 'Prediction'], ascending=[True, False])\nfinal_submission_df = submission_df.groupby('Protein Id').head(1500).reset_index(drop=True)\n\n# Save submission\nfinal_submission_df.to_csv('submission.tsv', sep='\\t', index=False, header=False)\n\nprint(f\"\\nSubmission file 'submission.tsv' created successfully!\")\nprint(f\"Total predictions in final submission: {len(final_submission_df):,}\")\nprint(f\"Unique proteins with predictions: {final_submission_df['Protein Id'].nunique():,}\")\n\nprint(\"\\nSubmission Preview:\")\ndisplay(final_submission_df.head(10))","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 7.4 Submission Statistics Visualization\nfig, gs = create_figure(1, 3, figsize=(20, 6), wspace=0.3)\n\n# Panel A: Predictions per Protein Distribution\nax1 = fig.add_subplot(gs[0, 0])\npreds_per_protein = final_submission_df.groupby('Protein Id').size()\nsns.histplot(preds_per_protein, ax=ax1, color=SCIENCE_PALETTE[0], bins=30, kde=True, alpha=0.7, edgecolor='white', linewidth=0.8)\nax1.axvline(preds_per_protein.median(), color=NATURE_PALETTE[0], linestyle='--', linewidth=2.5, label=f'Median: {preds_per_protein.median():.0f}')\nax1.legend(fontsize=10, frameon=True)\nstyle_axis(ax1, title='Predictions per Protein', xlabel='Number of Predictions', ylabel='Frequency')\n\n# Panel B: Prediction Score Distribution\nax2 = fig.add_subplot(gs[0, 1])\nsns.histplot(final_submission_df['Prediction'], ax=ax2, color=CELL_PALETTE[1], bins=50, kde=True, alpha=0.7, edgecolor='white', linewidth=0.8)\nax2.axvline(final_submission_df['Prediction'].median(), color=NATURE_PALETTE[0], linestyle='--', linewidth=2.5, label=f'Median: {final_submission_df[\"Prediction\"].median():.3f}')\nax2.legend(fontsize=10, frameon=True)\nstyle_axis(ax2, title='Prediction Score Distribution', xlabel='Prediction Score', ylabel='Frequency')\n\n# Panel C: Summary Statistics\nax3 = fig.add_subplot(gs[0, 2])\nstats_text = f\"Total Predictions\\n{len(final_submission_df):,}\\n\\nUnique Proteins\\n{final_submission_df['Protein Id'].nunique():,}\\n\\nUnique GO Terms\\n{final_submission_df['GO Term Id'].nunique():,}\\n\\nMean Score\\n{final_submission_df['Prediction'].mean():.4f}\"\nax3.text(0.5, 0.5, stats_text, ha='center', va='center', fontsize=16, fontweight='bold', color=CNS_COLORS['text'], fontfamily='serif', linespacing=1.8)\nax3.set_facecolor('#F0F3F5')\nfor spine in ax3.spines.values(): spine.set_visible(False)\nax3.set_xticks([]); ax3.set_yticks([])\nax3.set_title('Submission Summary', fontsize=13, weight='bold', fontfamily='serif', pad=15, color=CNS_COLORS['text'])\n\nadd_figure_labels(fig, 'Submission Statistics: Final Output Analysis',\n                  'Overview of the generated submission file including prediction counts and score distributions.')\nplt.tight_layout(rect=[0, 0.02, 1, 0.88])\nplt.show()","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## Conclusion\n\nThis notebook presents a complete pipeline for the CAFA 6 Protein Function Prediction challenge:\n\n1. **Data Loading**: Parsed FASTA sequences and tabular annotations\n2. **EDA**: Comprehensive visualization of dataset characteristics\n3. **Feature Engineering**: Extracted 71-dimensional physicochemical features\n4. **Modeling**: Trained ontology-specific LogisticRegression models\n5. **Evaluation**: Assessed performance using Precision-Recall metrics\n6. **Submission**: Generated predictions for the test set\n\n**Future Improvements:**\n- Explore protein language model embeddings (ESM2, ProtBERT)\n- Implement ensemble methods\n- Fine-tune probability thresholds per ontology\n","metadata":{}}]}