{"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":14875579,"sourceType":"competition"}],"dockerImageVersionId":31192,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# CAFA 6 Protein Function Prediction: FULL EDA \n\n\n## Reference: \n\n- https://www.kaggle.com/code/taylorsamarel/cafa-6-protein-function-starter-eda-model\n- https://www.kaggle.com/code/leonidkulyk/eda-cafa5-pfp-interactive-dags-plotly\n- https://www.kaggle.com/code/thedrcat/cafa-5-protein-function-prediction-eda/notebook\n\n\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"import subprocess\nimport sys\nimport warnings\nfrom pathlib import Path\nfrom collections import defaultdict, Counter\nfrom typing import Dict, List, Tuple\nimport pickle\nimport gc\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\nimport pandas as pd\nimport numpy as np\nimport networkx as nx\nimport obonet\nfrom Bio import SeqIO\nfrom Bio.SeqUtils.ProtParam import ProteinAnalysis\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport plotly.express as px\nimport plotly.graph_objects as go\nfrom plotly.subplots import make_subplots\n\nfrom sklearn.preprocessing import MultiLabelBinarizer\nfrom sklearn.decomposition import PCA\nfrom sklearn.manifold import TSNE\nfrom scipy.stats import chi2_contingency, spearmanr\nfrom scipy.cluster.hierarchy import dendrogram, linkage\n\nwarnings.filterwarnings('ignore')\nplt.style.use('fivethirtyeight')\nsns.set_palette(\"Reds\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-16T14:51:35.002347Z","iopub.execute_input":"2025-12-16T14:51:35.003122Z","iopub.status.idle":"2025-12-16T14:51:48.835299Z","shell.execute_reply.started":"2025-12-16T14:51:35.003074Z","shell.execute_reply":"2025-12-16T14:51:48.834508Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\nclass Config:\n    \"\"\"Centralized configuration\"\"\"\n    # Paths\n    BASE_PATH = Path('/kaggle/input/cafa-6-protein-function-prediction')\n    TRAIN_DIR = BASE_PATH / 'Train'\n    TEST_DIR = BASE_PATH / 'Test'\n    \n    # Sampling parameters\n    SAMPLE_PERCENT = 10  # Start with 10% for faster iteration\n    QUICK_MODE = True\n    RANDOM_SEED = 42\n    \n    # Analysis parameters\n    TOP_N_TERMS = 20\n    TOP_N_TAXA = 15\n    MIN_TERM_FREQUENCY = 5\n    \n    # Visualization\n    FIGURE_DPI = 100\n    COLOR_PALETTE = {\n        'MFO': '#FF6B6B',\n        'BPO': '#4ECDC4', \n        'CCO': '#45B7D1',\n        'F': '#FF6B6B',\n        'P': '#4ECDC4',\n        'C': '#45B7D1'\n    }\n    \n    # Cache\n    CACHE_DIR = Path('./cache')\n    USE_CACHE = True\n\n# Create cache directory\nConfig.CACHE_DIR.mkdir(exist_ok=True)\n\nprint(\"=\"*80)\nprint(\"CAFA 6 PROTEIN FUNCTION PREDICTION - OPTIMIZED EDA\")\nprint(f\"Sample: {Config.SAMPLE_PERCENT}% | Quick Mode: {Config.QUICK_MODE}\")\nprint(\"=\"*80)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-16T14:51:48.836496Z","iopub.execute_input":"2025-12-16T14:51:48.837281Z","iopub.status.idle":"2025-12-16T14:51:48.844830Z","shell.execute_reply.started":"2025-12-16T14:51:48.837254Z","shell.execute_reply":"2025-12-16T14:51:48.843887Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## UTILITY FUNCTIONS","metadata":{}},{"cell_type":"code","source":"\ndef cache_data(func):\n    \"\"\"Decorator for caching expensive operations\"\"\"\n    def wrapper(*args, **kwargs):\n        if not Config.USE_CACHE:\n            return func(*args, **kwargs)\n        \n        cache_file = Config.CACHE_DIR / f\"{func.__name__}_{Config.SAMPLE_PERCENT}.pkl\"\n        \n        if cache_file.exists():\n            print(f\"  Loading cached data for {func.__name__}\")\n            with open(cache_file, 'rb') as f:\n                return pickle.load(f)\n        \n        result = func(*args, **kwargs)\n        \n        with open(cache_file, 'wb') as f:\n            pickle.dump(result, f)\n        \n        return result\n    return wrapper\n\ndef memory_usage():\n    \"\"\"Print current memory usage\"\"\"\n    import psutil\n    process = psutil.Process()\n    mem = process.memory_info().rss / 1024 ** 3\n    print(f\"   Memory usage: {mem:.2f} GB\")\n\ndef create_subplot_figure(rows, cols, titles, figsize=(20, 12)):\n    \"\"\"Create standardized subplot figure\"\"\"\n    fig = make_subplots(\n        rows=rows, cols=cols,\n        subplot_titles=titles,\n        vertical_spacing=0.1,\n        horizontal_spacing=0.1\n    )\n    return fig","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-16T14:51:48.845769Z","iopub.execute_input":"2025-12-16T14:51:48.846014Z","iopub.status.idle":"2025-12-16T14:51:48.868421Z","shell.execute_reply.started":"2025-12-16T14:51:48.845995Z","shell.execute_reply":"2025-12-16T14:51:48.867576Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## DATA LOADING WITH OPTIMIZATION","metadata":{}},{"cell_type":"code","source":"\n@cache_data\ndef load_go_ontology():\n    \"\"\"Load and analyze GO ontology with caching\"\"\"\n    print(\"\\n[1/8]  Loading GO Ontology...\")\n    \n    graph = obonet.read_obo(Config.TRAIN_DIR / 'go-basic.obo')\n    print(f\"   ✓ Loaded {len(graph):,} GO terms\")\n    \n    # Extract metadata\n    term_to_ont = {}\n    term_names = {}\n    term_definitions = {}\n    \n    for term_id in graph.nodes():\n        node_data = graph.nodes[term_id]\n        \n        # Map namespace to ontology\n        if 'namespace' in node_data:\n            ns = node_data['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        \n        # Store names and definitions\n        if 'name' in node_data:\n            term_names[term_id] = node_data['name']\n        if 'def' in node_data:\n            term_definitions[term_id] = node_data['def']\n    \n    return {\n        'graph': graph,\n        'term_to_ont': term_to_ont,\n        'term_names': term_names,\n        'term_definitions': term_definitions\n    }\n\n@cache_data\ndef load_ia_weights():\n    \"\"\"Load IA weights with caching\"\"\"\n    print(\"\\n[2/8] ⚖️ Loading IA Weights...\")\n    \n    ia_df = pd.read_csv(Config.BASE_PATH / 'IA.tsv', sep='\\t', \n                        header=None, names=['term', 'ia'])\n    \n    print(f\"    Loaded {len(ia_df):,} IA weights\")\n    print(f\"    Range: [{ia_df['ia'].min():.3f}, {ia_df['ia'].max():.3f}]\")\n    \n    return ia_df\n\n@cache_data\ndef load_training_data():\n    \"\"\"Load training annotations and taxonomy\"\"\"\n    print(\"\\n[3/8]  Loading Training Data...\")\n    \n    # Load annotations\n    train_terms = pd.read_csv(\n        Config.TRAIN_DIR / 'train_terms.tsv', \n        sep='\\t',\n        names=['protein', 'term', 'aspect']\n    )\n    \n    # Load taxonomy\n    train_taxonomy = pd.read_csv(\n        Config.TRAIN_DIR / 'train_taxonomy.tsv',\n        sep='\\t',\n        names=['protein', 'taxon']\n    )\n    \n    print(f\"    Full dataset: {len(train_terms):,} annotations\")\n    print(f\"    Proteins: {train_terms['protein'].nunique():,}\")\n    print(f\"    GO terms: {train_terms['term'].nunique():,}\")\n    \n    # Sampling\n    if Config.SAMPLE_PERCENT < 100:\n        sample_proteins = train_terms['protein'].drop_duplicates().sample(\n            frac=Config.SAMPLE_PERCENT/100, \n            random_state=Config.RANDOM_SEED\n        ).tolist()\n        \n        train_terms = train_terms[train_terms['protein'].isin(sample_proteins)]\n        train_taxonomy = train_taxonomy[train_taxonomy['protein'].isin(sample_proteins)]\n        \n        print(f\"    Sampled to {Config.SAMPLE_PERCENT}%: {len(train_terms):,} annotations\")\n    \n    return train_terms, train_taxonomy\n\ndef load_sequences_efficiently(protein_ids: set, max_proteins: int = None):\n    \"\"\"\n    Load sequences efficiently with progress tracking\n    \n    Args:\n        protein_ids: Set of protein IDs to load\n        max_proteins: Maximum number of proteins to load (None for all)\n    \"\"\"\n    print(\"\\n[4/8]  Loading Protein Sequences...\")\n    \n    sequences = {}\n    loaded_count = 0\n    target_count = min(len(protein_ids), max_proteins) if max_proteins else len(protein_ids)\n    \n    with open(Config.TRAIN_DIR / 'train_sequences.fasta', 'r') as fasta_file:\n        for record in SeqIO.parse(fasta_file, 'fasta'):\n            # Extract protein ID\n            pid = record.id.split('|')[1] if '|' in record.id else record.id\n            if pid in protein_ids:\n                sequences[pid] = str(record.seq)\n                loaded_count += 1\n                \n                # Progress indicator\n                if loaded_count % 5000 == 0:\n                    progress = (loaded_count / target_count) * 100\n                    print(f\"    Progress: {loaded_count:,}/{target_count:,} ({progress:.1f}%)\")\n                \n                if max_proteins and loaded_count >= max_proteins:\n                    break\n    \n    print(f\"    Loaded {len(sequences):,} sequences\")\n    return sequences\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-16T14:51:48.869433Z","iopub.execute_input":"2025-12-16T14:51:48.869863Z","iopub.status.idle":"2025-12-16T14:51:48.899539Z","shell.execute_reply.started":"2025-12-16T14:51:48.869830Z","shell.execute_reply":"2025-12-16T14:51:48.898739Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## FEATURE EXTRACTION","metadata":{}},{"cell_type":"code","source":"\n\ndef extract_sequence_features(sequence: str) -> Dict:\n    \"\"\"\n    Extract comprehensive sequence features\n    \n    Features:\n    - Length\n    - Molecular weight\n    - Aromaticity\n    - Instability index\n    - Isoelectric point\n    - Gravy (hydrophobicity)\n    - Secondary structure fractions\n    - Amino acid composition\n    \"\"\"\n    try:\n        analyzer = ProteinAnalysis(sequence)\n        \n        # Handle potential errors in secondary structure calculation\n        try:\n            helix, turn, sheet = analyzer.secondary_structure_fraction()\n        except:\n            helix, turn, sheet = 0, 0, 0\n        \n        features = {\n            'length': len(sequence),\n            'molecular_weight': analyzer.molecular_weight(),\n            'aromaticity': analyzer.aromaticity(),\n            'instability_index': analyzer.instability_index(),\n            'isoelectric_point': analyzer.isoelectric_point(),\n            'gravy': analyzer.gravy(),\n            'helix_fraction': helix,\n            'turn_fraction': turn,\n            'sheet_fraction': sheet,\n        }\n        \n        # Add amino acid composition\n        aa_comp = analyzer.get_amino_acids_percent()\n        for aa, percent in aa_comp.items():\n            features[f'aa_{aa}'] = percent\n        \n        return features\n    \n    except Exception as e:\n        # Return default features on error\n        return {\n            'length': len(sequence),\n            'molecular_weight': np.nan,\n            'aromaticity': np.nan,\n            'instability_index': np.nan,\n            'isoelectric_point': np.nan,\n            'gravy': np.nan,\n            'helix_fraction': np.nan,\n            'turn_fraction': np.nan,\n            'sheet_fraction': np.nan,\n        }\n\ndef compute_sequence_features_batch(sequences: Dict[str, str], sample_size: int = None):\n    \"\"\"Compute features for multiple sequences with sampling\"\"\"\n    print(\"\\n[5/8]  Computing Sequence Features...\")\n    \n    if sample_size:\n        sample_proteins = list(sequences.keys())[:sample_size]\n        sequences = {k: sequences[k] for k in sample_proteins}\n    \n    features_list = []\n    \n    for i, (protein_id, seq) in enumerate(sequences.items()):\n        if i % 1000 == 0 and i > 0:\n            print(f\"    Processed {i:,}/{len(sequences):,} proteins\")\n        \n        features = extract_sequence_features(seq)\n        features['protein'] = protein_id\n        features_list.append(features)\n    \n    features_df = pd.DataFrame(features_list)\n    print(f\"    Computed features for {len(features_df):,} proteins\")\n    \n    return features_df","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-16T14:51:48.901729Z","iopub.execute_input":"2025-12-16T14:51:48.901999Z","iopub.status.idle":"2025-12-16T14:51:48.921733Z","shell.execute_reply.started":"2025-12-16T14:51:48.901973Z","shell.execute_reply":"2025-12-16T14:51:48.920742Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## VISUALIZATIONS","metadata":{}},{"cell_type":"code","source":"\ndef visualize_go_hierarchy(go_data, train_terms, ia_df):\n    \"\"\"Enhanced GO hierarchy visualization\"\"\"\n    print(\"\\n[6/8]  Visualizing GO Hierarchy...\")\n    \n    graph = go_data['graph']\n    term_to_ont = go_data['term_to_ont']\n    term_names = go_data['term_names']\n    \n    # Compute hierarchy depths (sample for efficiency)\n    sample_terms = list(term_to_ont.keys())[:500]\n    term_depths = {}\n    \n    for term in sample_terms:\n        try:\n            paths = []\n            for root in ['GO:0008150', 'GO:0005575', 'GO:0003674']:\n                if nx.has_path(graph, term, root):\n                    paths.append(nx.shortest_path_length(graph, term, root))\n            term_depths[term] = max(paths) if paths else 0\n        except:\n            term_depths[term] = 0\n    \n    # Create comprehensive figure\n    fig = plt.figure(figsize=(20, 14))\n    gs = fig.add_gridspec(3, 3, hspace=0.3, wspace=0.3)\n    \n    # 1. Ontology distribution\n    ax1 = fig.add_subplot(gs[0, 0])\n    ont_counts = pd.Series(term_to_ont).value_counts()\n    colors = [Config.COLOR_PALETTE.get(k, '#CCCCCC') for k in ont_counts.index]\n    \n    bars = ax1.bar(range(len(ont_counts)), ont_counts.values, \n                   color=colors, edgecolor='black', linewidth=2, alpha=0.8)\n    ax1.set_xticks(range(len(ont_counts)))\n    ax1.set_xticklabels(ont_counts.index, fontsize=11, fontweight='bold')\n    ax1.set_title('GO Terms by Ontology', fontsize=13, fontweight='bold', pad=15)\n    ax1.set_ylabel('Count', fontsize=11)\n    ax1.grid(axis='y', alpha=0.3)\n    \n    for bar, val in zip(bars, ont_counts.values):\n        height = bar.get_height()\n        ax1.text(bar.get_x() + bar.get_width()/2., height,\n                f'{val:,}\\n({val/ont_counts.sum()*100:.1f}%)',\n                ha='center', va='bottom', fontweight='bold', fontsize=10)\n    \n    # 2. Hierarchy depth distribution\n    ax2 = fig.add_subplot(gs[0, 1])\n    depth_values = list(term_depths.values())\n    ax2.hist(depth_values, bins=20, color='#A8E6CF', edgecolor='black', alpha=0.7)\n    ax2.set_title('GO Term Hierarchy Depth', fontsize=12, fontweight='bold')\n    ax2.set_xlabel('Depth Level', fontsize=10)\n    ax2.set_ylabel('Frequency', fontsize=10)\n    ax2.axvline(np.mean(depth_values), color='red', linestyle='--', \n                linewidth=2, label=f'Mean: {np.mean(depth_values):.1f}')\n    ax2.legend()\n    \n    # 3. IA distribution by ontology\n    ax3 = fig.add_subplot(gs[0, 2])\n    ia_by_ont = ia_df.copy()\n    ia_by_ont['ontology'] = ia_by_ont['term'].map(term_to_ont)\n    ia_by_ont = ia_by_ont.dropna()\n    \n    ont_order = ['MFO', 'BPO', 'CCO']\n    ia_data = [ia_by_ont[ia_by_ont['ontology']==ont]['ia'].values \n               for ont in ont_order if ont in ia_by_ont['ontology'].unique()]\n    \n    bp = ax3.boxplot(ia_data, labels=[o for o in ont_order if o in ia_by_ont['ontology'].unique()],\n                     patch_artist=True, showmeans=True)\n    \n    for patch, ont in zip(bp['boxes'], ont_order):\n        if ont in ia_by_ont['ontology'].unique():\n            patch.set_facecolor(Config.COLOR_PALETTE[ont])\n            patch.set_alpha(0.7)\n    \n    ax3.set_title('IA Weights by Ontology', fontsize=12, fontweight='bold')\n    ax3.set_ylabel('IA Weight', fontsize=10)\n    ax3.grid(axis='y', alpha=0.3)\n    \n    # 4. Term frequency in training data\n    ax4 = fig.add_subplot(gs[1, :2])\n    term_freq = train_terms['term'].value_counts().head(Config.TOP_N_TERMS)\n    \n    # Map to names and ontologies\n    term_labels = [f\"{term_names.get(t, t)[:40]}\" for t in term_freq.index]\n    term_colors = [Config.COLOR_PALETTE.get(term_to_ont.get(t, 'MFO'), '#CCCCCC') \n                   for t in term_freq.index]\n    \n    ax4.barh(range(len(term_freq)), term_freq.values, color=term_colors, \n             edgecolor='black', alpha=0.7)\n    ax4.set_yticks(range(len(term_freq)))\n    ax4.set_yticklabels(term_labels, fontsize=9)\n    ax4.set_title(f'Top {Config.TOP_N_TERMS} Most Frequent GO Terms', \n                  fontsize=12, fontweight='bold')\n    ax4.set_xlabel('Frequency in Training Set', fontsize=10)\n    ax4.invert_yaxis()\n    \n    # 5. Network visualization\n    ax5 = fig.add_subplot(gs[1, 2])\n    sample_for_network = list(term_to_ont.keys())[:30]\n    subgraph = graph.subgraph(sample_for_network)\n    \n    pos = nx.spring_layout(subgraph, k=0.8, iterations=50, seed=Config.RANDOM_SEED)\n    node_colors = [Config.COLOR_PALETTE.get(term_to_ont.get(node, 'MFO'), '#CCCCCC') \n                   for node in subgraph.nodes()]\n    \n    nx.draw_networkx_nodes(subgraph, pos, node_color=node_colors,\n                          node_size=200, alpha=0.7, ax=ax5)\n    nx.draw_networkx_edges(subgraph, pos, alpha=0.2, arrows=True,\n                          arrowsize=8, ax=ax5, edge_color='gray')\n    \n    ax5.set_title('GO Network Sample', fontsize=12, fontweight='bold')\n    ax5.axis('off')\n    \n    # 6. Statistics panel\n    ax6 = fig.add_subplot(gs[2, :])\n    ax6.axis('off')\n    \n    stats_text = f\"\"\"\n    GO ONTOLOGY STATISTICS\n    \n    Total Terms: {len(term_to_ont):,}\n    • Molecular Function: {ont_counts.get('MFO', 0):,}\n    • Biological Process: {ont_counts.get('BPO', 0):,}\n    • Cellular Component: {ont_counts.get('CCO', 0):,}\n    \n    Hierarchy:\n    • Mean depth: {np.mean(depth_values):.2f}\n    • Max depth: {max(depth_values)}\n    • Min depth: {min(depth_values)}\n    \n    IA Weights:\n    • Overall range: [{ia_df['ia'].min():.3f}, {ia_df['ia'].max():.3f}]\n    • Mean ± Std: {ia_df['ia'].mean():.3f} ± {ia_df['ia'].std():.3f}\n    \"\"\"\n    \n    ax6.text(0.1, 0.5, stats_text, fontsize=11, family='monospace',\n            verticalalignment='center',\n            bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.3))\n    \n    plt.suptitle('Gene Ontology Comprehensive Analysis', \n                 fontsize=16, fontweight='bold', y=0.995)\n    plt.tight_layout()\n    plt.savefig('go_analysis.png', dpi=Config.FIGURE_DPI, bbox_inches='tight')\n    plt.show()\n    \n    print(\"    GO hierarchy visualization complete\")\n\ndef visualize_training_data(train_terms, train_taxonomy, go_data):\n    \"\"\"Enhanced training data visualization\"\"\"\n    print(\"\\n[7/8]  Visualizing Training Data...\")\n    \n    term_names = go_data['term_names']\n    term_to_ont = go_data['term_to_ont']\n    \n    fig = plt.figure(figsize=(22, 14))\n    gs = fig.add_gridspec(4, 4, hspace=0.35, wspace=0.35)\n    \n    # 1. Ontology distribution in training\n    ax1 = fig.add_subplot(gs[0, 0])\n    aspect_counts = train_terms['aspect'].value_counts()\n    colors_list = [Config.COLOR_PALETTE.get(k, '#CCCCCC') for k in aspect_counts.index]\n    \n    bars = ax1.bar(range(len(aspect_counts)), aspect_counts.values,\n                   color=colors_list, edgecolor='black', linewidth=1.5, alpha=0.8)\n    ax1.set_xticks(range(len(aspect_counts)))\n    ax1.set_xticklabels(aspect_counts.index, fontsize=10, fontweight='bold')\n    ax1.set_title('Training Annotations by Ontology', fontsize=11, fontweight='bold')\n    ax1.set_ylabel('Count', fontsize=10)\n    \n    for bar, val in zip(bars, aspect_counts.values):\n        ax1.text(bar.get_x() + bar.get_width()/2., val,\n                f'{val:,}', ha='center', va='bottom', fontweight='bold', fontsize=9)\n    \n    # 2. Terms per protein distribution\n    ax2 = fig.add_subplot(gs[0, 1])\n    terms_per_protein = train_terms.groupby('protein').size()\n    ax2.hist(terms_per_protein, bins=50, color='#FFD93D', edgecolor='black', alpha=0.7)\n    ax2.set_title('Terms per Protein', fontsize=11, fontweight='bold')\n    ax2.set_xlabel('Number of Terms', fontsize=9)\n    ax2.set_ylabel('Frequency', fontsize=9)\n    ax2.axvline(terms_per_protein.mean(), color='red', linestyle='--', linewidth=2,\n                label=f'Mean: {terms_per_protein.mean():.1f}')\n    ax2.axvline(terms_per_protein.median(), color='blue', linestyle='--', linewidth=2,\n                label=f'Median: {terms_per_protein.median():.0f}')\n    ax2.legend(fontsize=8)\n    \n    # 3. Proteins per term distribution (log scale)\n    ax3 = fig.add_subplot(gs[0, 2])\n    proteins_per_term = train_terms.groupby('term').size()\n    ax3.hist(proteins_per_term, bins=50, color='#FFEAA7', \n             edgecolor='black', alpha=0.7, log=True)\n    ax3.set_title('Proteins per Term (log scale)', fontsize=11, fontweight='bold')\n    ax3.set_xlabel('Number of Proteins', fontsize=9)\n    ax3.set_ylabel('Number of Terms (log)', fontsize=9)\n    \n    # 4. Taxonomy distribution\n    ax4 = fig.add_subplot(gs[0, 3])\n    top_taxa = train_taxonomy['taxon'].value_counts().head(10)\n    ax4.bar(range(len(top_taxa)), top_taxa.values, color='#74B9FF', edgecolor='black')\n    ax4.set_xticks(range(len(top_taxa)))\n    ax4.set_xticklabels([str(t)[:8] for t in top_taxa.index], \n                        rotation=45, ha='right', fontsize=8)\n    ax4.set_title('Top 10 Species', fontsize=11, fontweight='bold')\n    ax4.set_ylabel('Number of Proteins', fontsize=9)\n    \n    # 5. Term frequency bins\n    ax5 = fig.add_subplot(gs[1, :2])\n    freq_bins = pd.cut(proteins_per_term, \n                       bins=[0, 10, 50, 100, 500, 10000],\n                       labels=['1-10', '11-50', '51-100', '101-500', '>500'])\n    freq_dist = freq_bins.value_counts().sort_index()\n    \n    ax5.bar(range(len(freq_dist)), freq_dist.values, \n            color='#E17055', edgecolor='black', alpha=0.7)\n    ax5.set_xticks(range(len(freq_dist)))\n    ax5.set_xticklabels(freq_dist.index, fontsize=10)\n    ax5.set_title('GO Term Frequency Distribution', fontsize=11, fontweight='bold')\n    ax5.set_xlabel('Proteins Annotated with Term', fontsize=10)\n    ax5.set_ylabel('Number of Terms', fontsize=10)\n    \n    for i, val in enumerate(freq_dist.values):\n        ax5.text(i, val, f'{val:,}', ha='center', va='bottom', fontweight='bold')\n    \n    # 6. Cross-ontology annotation patterns\n    ax6 = fig.add_subplot(gs[1, 2:])\n    protein_ont_counts = train_terms.groupby('protein')['aspect'].value_counts().unstack(fill_value=0)\n    \n    if protein_ont_counts.shape[0] >= 2:\n        # Create scatter plot\n        aspects = protein_ont_counts.columns.tolist()\n        if len(aspects) >= 2:\n            ax6.scatter(protein_ont_counts[aspects[0]], \n                       protein_ont_counts[aspects[1]],\n                       alpha=0.5, s=20, color='#6C5CE7')\n            ax6.set_xlabel(f'{aspects[0]} annotations', fontsize=10)\n            ax6.set_ylabel(f'{aspects[1]} annotations', fontsize=10)\n            ax6.set_title('Cross-Ontology Annotation Pattern', fontsize=11, fontweight='bold')\n            ax6.grid(alpha=0.3)\n    \n    # 7. Top terms by ontology\n    ax7 = fig.add_subplot(gs[2, :])\n    top_terms_by_ont = train_terms.groupby(['aspect', 'term']).size().reset_index(name='count')\n    \n    y_pos = 0\n    colors_by_aspect = []\n    labels = []\n    values = []\n    \n    for aspect in ['F', 'P', 'C']:\n        if aspect in top_terms_by_ont['aspect'].values:\n            top_5 = top_terms_by_ont[top_terms_by_ont['aspect']==aspect].nlargest(5, 'count')\n            for _, row in top_5.iterrows():\n                term_name = term_names.get(row['term'], row['term'])[:35]\n                labels.append(f\"{aspect}: {term_name}\")\n                values.append(row['count'])\n                colors_by_aspect.append(Config.COLOR_PALETTE.get(aspect, '#CCCCCC'))\n    \n    ax7.barh(range(len(labels)), values, color=colors_by_aspect, \n             edgecolor='black', alpha=0.7)\n    ax7.set_yticks(range(len(labels)))\n    ax7.set_yticklabels(labels, fontsize=8)\n    ax7.set_title('Top 5 Terms per Ontology', fontsize=12, fontweight='bold')\n    ax7.set_xlabel('Frequency', fontsize=10)\n    ax7.invert_yaxis()\n\n    \n    ","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-16T14:51:48.922796Z","iopub.execute_input":"2025-12-16T14:51:48.923062Z","iopub.status.idle":"2025-12-16T14:51:48.961202Z","shell.execute_reply.started":"2025-12-16T14:51:48.923040Z","shell.execute_reply":"2025-12-16T14:51:48.960411Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\ndef visualize_sequence_features(features_df, train_terms):\n    \"\"\"Visualize protein sequence features\"\"\"\n    print(\"\\n[8/8]  Visualizing Sequence Features...\")\n    \n    # Merge with annotations\n    merged_df = features_df.merge(\n        train_terms.groupby('protein')['aspect'].apply(lambda x: ','.join(x.unique())),\n        left_on='protein', right_index=True, how='left'\n    )\n    \n    fig = plt.figure(figsize=(22, 14))\n    gs = fig.add_gridspec(3, 4, hspace=0.35, wspace=0.35)\n    \n    # 1. Sequence length distribution\n    ax1 = fig.add_subplot(gs[0, 0])\n    ax1.hist(features_df['length'], bins=50, color='#95E1D3', \n             edgecolor='black', alpha=0.7)\n    ax1.set_title('Sequence Length Distribution', fontsize=11, fontweight='bold')\n    ax1.set_xlabel('Length (amino acids)', fontsize=9)\n    ax1.set_ylabel('Frequency', fontsize=9)\n    ax1.axvline(features_df['length'].mean(), color='red', linestyle='--', linewidth=2)\n    \n    # 2. Molecular weight\n    ax2 = fig.add_subplot(gs[0, 1])\n    ax2.hist(features_df['molecular_weight'].dropna(), bins=50, \n             color='#FAD390', edgecolor='black', alpha=0.7)\n    ax2.set_title('Molecular Weight Distribution', fontsize=11, fontweight='bold')\n    ax2.set_xlabel('Molecular Weight (Da)', fontsize=9)\n    ax2.set_ylabel('Frequency', fontsize=9)\n    \n    # 3. Isoelectric point\n    ax3 = fig.add_subplot(gs[0, 2])\n    ax3.hist(features_df['isoelectric_point'].dropna(), bins=50,\n             color='#A8E6CF', edgecolor='black', alpha=0.7)\n    ax3.set_title('Isoelectric Point Distribution', fontsize=11, fontweight='bold')\n    ax3.set_xlabel('pI', fontsize=9)\n    ax3.set_ylabel('Frequency', fontsize=9)\n    \n    # 4. GRAVY (hydrophobicity)\n    ax4 = fig.add_subplot(gs[0, 3])\n    ax4.hist(features_df['gravy'].dropna(), bins=50,\n             color='#FF6B6B', edgecolor='black', alpha=0.7)\n    ax4.set_title('GRAVY (Hydrophobicity)', fontsize=11, fontweight='bold')\n    ax4.set_xlabel('GRAVY Score', fontsize=9)\n    ax4.set_ylabel('Frequency', fontsize=9)\n    ax4.axvline(0, color='black', linestyle='--', linewidth=1, alpha=0.5)\n    \n    # 5. Secondary structure composition\n    ax5 = fig.add_subplot(gs[1, 0])\n    structure_data = features_df[['helix_fraction', 'turn_fraction', 'sheet_fraction']].dropna()\n    if len(structure_data) > 0:\n        structure_means = structure_data.mean()\n        ax5.bar(['Helix', 'Turn', 'Sheet'], structure_means.values,\n               color=['#FF6B6B', '#4ECDC4', '#45B7D1'], edgecolor='black', alpha=0.7)\n        ax5.set_title('Average Secondary Structure', fontsize=11, fontweight='bold')\n        ax5.set_ylabel('Fraction', fontsize=9)\n        ax5.set_ylim([0, 1])\n    \n    # 6. Aromaticity vs Instability\n    ax6 = fig.add_subplot(gs[1, 1])\n    valid_data = features_df[['aromaticity', 'instability_index']].dropna()\n    if len(valid_data) > 0:\n        ax6.scatter(valid_data['aromaticity'], valid_data['instability_index'],\n                   alpha=0.5, s=10, color='#6C5CE7')\n        ax6.set_xlabel('Aromaticity', fontsize=9)\n        ax6.set_ylabel('Instability Index', fontsize=9)\n        ax6.set_title('Aromaticity vs Instability', fontsize=11, fontweight='bold')\n        ax6.axhline(40, color='red', linestyle='--', alpha=0.5, label='Stability threshold')\n        ax6.legend(fontsize=8)\n    \n    # 7. Length vs Molecular Weight\n    ax7 = fig.add_subplot(gs[1, 2])\n    valid_data = features_df[['length', 'molecular_weight']].dropna()\n    if len(valid_data) > 0:\n        ax7.scatter(valid_data['length'], valid_data['molecular_weight'],\n                   alpha=0.5, s=10, color='#00B894')\n        ax7.set_xlabel('Sequence Length', fontsize=9)\n        ax7.set_ylabel('Molecular Weight (Da)', fontsize=9)\n        ax7.set_title('Length vs Molecular Weight', fontsize=11, fontweight='bold')\n        \n        # Add correlation\n        corr = valid_data['length'].corr(valid_data['molecular_weight'])\n        ax7.text(0.05, 0.95, f'Corr: {corr:.3f}', transform=ax7.transAxes,\n                fontsize=9, verticalalignment='top',\n                bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))\n    \n    # 8. Top amino acid composition\n    ax8 = fig.add_subplot(gs[1, 3])\n    aa_cols = [col for col in features_df.columns if col.startswith('aa_')]\n    if aa_cols:\n        aa_means = features_df[aa_cols].mean().sort_values(ascending=False).head(10)\n        aa_names = [col.replace('aa_', '') for col in aa_means.index]\n        \n        ax8.barh(range(len(aa_means)), aa_means.values,\n                color='#FDCB6E', edgecolor='black', alpha=0.7)\n        ax8.set_yticks(range(len(aa_means)))\n        ax8.set_yticklabels(aa_names, fontsize=9)\n        ax8.set_title('Top 10 Amino Acids', fontsize=11, fontweight='bold')\n        ax8.set_xlabel('Average Composition (%)', fontsize=9)\n        ax8.invert_yaxis()\n    \n    # 9. Feature correlation heatmap\n    ax9 = fig.add_subplot(gs[2, :2])\n    numeric_features = ['length', 'molecular_weight', 'aromaticity', \n                       'instability_index', 'isoelectric_point', 'gravy']\n    corr_data = features_df[numeric_features].dropna()\n    \n    if len(corr_data) > 0:\n        corr_matrix = corr_data.corr()\n        im = ax9.imshow(corr_matrix, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto')\n        \n        ax9.set_xticks(range(len(numeric_features)))\n        ax9.set_yticks(range(len(numeric_features)))\n        ax9.set_xticklabels([f.replace('_', ' ').title()[:10] for f in numeric_features],\n                           rotation=45, ha='right', fontsize=8)\n        ax9.set_yticklabels([f.replace('_', ' ').title()[:10] for f in numeric_features],\n                           fontsize=8)\n        ax9.set_title('Feature Correlation Matrix', fontsize=11, fontweight='bold')\n        \n        # Add correlation values\n        for i in range(len(numeric_features)):\n            for j in range(len(numeric_features)):\n                text = ax9.text(j, i, f'{corr_matrix.iloc[i, j]:.2f}',\n                              ha=\"center\", va=\"center\", color=\"black\", fontsize=7)\n        \n        plt.colorbar(im, ax=ax9, label='Correlation')\n    \n\n    \n    return features_df","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-16T14:51:48.962296Z","iopub.execute_input":"2025-12-16T14:51:48.962869Z","iopub.status.idle":"2025-12-16T14:51:48.988827Z","shell.execute_reply.started":"2025-12-16T14:51:48.962836Z","shell.execute_reply":"2025-12-16T14:51:48.987849Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n\n# Load data\ngo_data = load_go_ontology()\nia_df = load_ia_weights()\ntrain_terms, train_taxonomy = load_training_data()\n\n# Load sequences\nprotein_ids = set(train_terms['protein'].unique())\nsequences = load_sequences_efficiently(protein_ids)\n\n# Visualize GO and training data\nvisualize_go_hierarchy(go_data, train_terms, ia_df)\nvisualize_training_data(train_terms, train_taxonomy, go_data)\n\n# Compute and visualize sequence features\nif not Config.QUICK_MODE:\n    features_df = compute_sequence_features_batch(sequences, sample_size=5000)\n    visualize_sequence_features(features_df, train_terms)\nelse:\n    print(\"\\n Quick mode: Skipping intensive feature computation\")\n    print(\"   Set QUICK_MODE=False to enable full feature analysis\")\n\n# Memory cleanup\ngc.collect()\nmemory_usage()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-16T14:51:48.989831Z","iopub.execute_input":"2025-12-16T14:51:48.990162Z","iopub.status.idle":"2025-12-16T14:52:01.835400Z","shell.execute_reply.started":"2025-12-16T14:51:48.990134Z","shell.execute_reply":"2025-12-16T14:52:01.834462Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}