{"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":"!pip install seqio\n!pip install biopython plotly obonet networkx -q\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T09:03:59.842044Z","iopub.execute_input":"2025-10-16T09:03:59.842361Z","iopub.status.idle":"2025-10-16T09:03:59.870455Z","shell.execute_reply.started":"2025-10-16T09:03:59.842328Z","shell.execute_reply":"2025-10-16T09:03:59.869356Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom Bio import SeqIO\nimport plotly.express as px\nimport plotly.graph_objects as go\nfrom plotly.subplots import make_subplots\nimport obonet\nimport networkx as nx\nfrom collections import Counter\nimport warnings\nwarnings.filterwarnings('ignore')\n\nprint(\"All packages imported successfully!\")\n\n\n# Load all data files\ndef load_cafa6_data():\n    \"\"\"Load all CAFA 6 competition data files\"\"\"\n    data = {}\n    \n    print(\"Loading training sequences...\")\n    train_sequences = {}\n    for record in SeqIO.parse(\"/kaggle/input/cafa-6-protein-function-prediction/Train/train_sequences.fasta\", \"fasta\"):\n        protein_id = record.id.split('|')[1]\n        train_sequences[protein_id] = str(record.seq)\n    data['train_sequences'] = train_sequences\n    \n    print(\"Loading test sequences...\")\n    test_sequences = {}\n    for record in SeqIO.parse(\"/kaggle/input/cafa-6-protein-function-prediction/Test/testsuperset.fasta\", \"fasta\"):\n        protein_id = record.id\n        test_sequences[protein_id] = str(record.seq)\n    data['test_sequences'] = test_sequences\n    \n    print(\"Loading training terms...\")\n    train_terms = pd.read_csv(\"/kaggle/input/cafa-6-protein-function-prediction/Train/train_terms.tsv\", sep='\\t')\n    data['train_terms'] = train_terms\n    \n    print(\"Loading taxonomy...\")\n    train_taxonomy = pd.read_csv(\"/kaggle/input/cafa-6-protein-function-prediction/Train/train_taxonomy.tsv\", sep='\\t')\n    data['train_taxonomy'] = train_taxonomy\n    \n    print(\"Loading IA weights...\")\n    ia_weights = pd.read_csv(\"/kaggle/input/cafa-6-protein-function-prediction/IA.tsv\", sep='\\t')\n    data['ia_weights'] = ia_weights\n    \n    print(\"Loading GO ontology...\")\n    go_graph = obonet.read_obo(\"/kaggle/input/cafa-6-protein-function-prediction/Train/go-basic.obo\")\n    data['go_graph'] = go_graph\n    \n    print(\"Loading selected taxa...\")\n    selected_taxa = pd.read_csv(\"/kaggle/input/cafa-6-protein-function-prediction/Test/testsuperset-taxon-list.tsv\", sep='\\t')\n    data['selected_taxa'] = selected_taxa\n    \n    return data\n\n# Load the data\ndata = load_cafa6_data()\n\nprint(\"\\n\" + \"=\"*50)\nprint(\"DATA LOADING COMPLETE\")\nprint(\"=\"*50)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-10-16T09:03:59.871591Z","iopub.execute_input":"2025-10-16T09:03:59.871953Z","iopub.status.idle":"2025-10-16T09:04:14.638874Z","shell.execute_reply.started":"2025-10-16T09:03:59.871919Z","shell.execute_reply":"2025-10-16T09:04:14.638055Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"data.keys()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T09:04:14.640768Z","iopub.execute_input":"2025-10-16T09:04:14.641063Z","iopub.status.idle":"2025-10-16T09:04:14.647424Z","shell.execute_reply.started":"2025-10-16T09:04:14.641041Z","shell.execute_reply":"2025-10-16T09:04:14.646532Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"data['train_taxonomy']","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-16T09:04:14.648336Z","iopub.execute_input":"2025-10-16T09:04:14.64857Z","iopub.status.idle":"2025-10-16T09:04:14.684098Z","shell.execute_reply.started":"2025-10-16T09:04:14.64855Z","shell.execute_reply":"2025-10-16T09:04:14.683192Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\ndef print_basic_stats(data):\n    \"\"\"Print basic statistics about the dataset\"\"\"\n    print(\"📊 BASIC DATASET STATISTICS\")\n    print(\"-\" * 40)\n    \n    print(f\"🏋️ Training proteins: {len(data['train_sequences']):,}\")\n    print(f\"🎯 Test proteins: {len(data['test_sequences']):,}\")\n    print(f\"🏷️ Training GO term annotations: {len(data['train_terms']):,}\")\n    print(f\"🔤 Unique GO terms: {data['train_terms']['term'].nunique():,}\")\n    \n    # Count terms by ontology\n    term_counts = data['train_terms']['aspect'].value_counts()\n    print(f\"\\n📈 Terms by ontology:\")\n    for aspect, count in term_counts.items():\n        print(f\"   {aspect}: {count:,} terms\")\n    \n    # Sequence statistics\n    train_seqs = list(data['train_sequences'].values())\n    test_seqs = list(data['test_sequences'].values())\n    \n    print(f\"\\n📏 Sequence length statistics:\")\n    print(f\"   Training - Mean: {np.mean([len(seq) for seq in train_seqs]):.1f}, \"\n          f\"Max: {max([len(seq) for seq in train_seqs]):,}, \"\n          f\"Min: {min([len(seq) for seq in train_seqs])}\")\n    print(f\"   Test - Mean: {np.mean([len(seq) for seq in test_seqs]):.1f}, \"\n          f\"Max: {max([len(seq) for seq in test_seqs]):,}, \"\n          f\"Min: {min([len(seq) for seq in test_seqs])}\")\n    \n    # Taxonomy statistics\n    print(f\"\\n🧬 Unique species in training: {data['train_taxonomy']['9606'].nunique():,}\")\n\nprint_basic_stats(data)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-10-16T09:04:14.685314Z","iopub.execute_input":"2025-10-16T09:04:14.68562Z","iopub.status.idle":"2025-10-16T09:04:14.989221Z","shell.execute_reply.started":"2025-10-16T09:04:14.685591Z","shell.execute_reply":"2025-10-16T09:04:14.988292Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\ndef analyze_sequences(data):\n    \"\"\"Analyze protein sequences and amino acid composition\"\"\"\n    print(\"\\n🔬 SEQUENCE ANALYSIS\")\n    print(\"-\" * 40)\n    \n    train_seqs = list(data['train_sequences'].values())\n    test_seqs = list(data['test_sequences'].values())\n    \n    # Amino acid composition\n    aa_list = list('ACDEFGHIKLMNPQRSTVWY')  # 20 standard amino acids\n    \n    def get_aa_composition(sequences):\n        all_aas = ''.join(sequences)\n        aa_counts = Counter(all_aas)\n        total = sum(aa_counts.values())\n        return {aa: aa_counts.get(aa, 0)/total for aa in aa_list}\n    \n    train_aa_comp = get_aa_composition(train_seqs)\n    test_aa_comp = get_aa_composition(test_seqs)\n    \n    # Create comparison plot\n    fig, axes = plt.subplots(2, 2, figsize=(16, 12))\n    \n    # Training set AA composition\n    axes[0,0].bar(train_aa_comp.keys(), train_aa_comp.values(), color='skyblue', alpha=0.7)\n    axes[0,0].set_title('Training Set - Amino Acid Composition', fontsize=14, fontweight='bold')\n    axes[0,0].set_ylabel('Frequency', fontsize=12)\n    axes[0,0].grid(True, alpha=0.3)\n    \n    # Test set AA composition\n    axes[0,1].bar(test_aa_comp.keys(), test_aa_comp.values(), color='lightcoral', alpha=0.7)\n    axes[0,1].set_title('Test Set - Amino Acid Composition', fontsize=14, fontweight='bold')\n    axes[0,1].set_ylabel('Frequency', fontsize=12)\n    axes[0,1].grid(True, alpha=0.3)\n    \n    # Sequence length distribution\n    train_lengths = [len(seq) for seq in train_seqs]\n    test_lengths = [len(seq) for seq in test_seqs]\n    \n    axes[1,0].hist(train_lengths, bins=100, alpha=0.7, color='skyblue', label='Training')\n    axes[1,0].set_title('Training Set - Sequence Length Distribution', fontsize=14, fontweight='bold')\n    axes[1,0].set_xlabel('Sequence Length', fontsize=12)\n    axes[1,0].set_ylabel('Frequency', fontsize=12)\n    axes[1,0].legend()\n    axes[1,0].grid(True, alpha=0.3)\n    \n    axes[1,1].hist(test_lengths, bins=100, alpha=0.7, color='lightcoral', label='Test')\n    axes[1,1].set_title('Test Set - Sequence Length Distribution', fontsize=14, fontweight='bold')\n    axes[1,1].set_xlabel('Sequence Length', fontsize=12)\n    axes[1,1].set_ylabel('Frequency', fontsize=12)\n    axes[1,1].legend()\n    axes[1,1].grid(True, alpha=0.3)\n    \n    plt.tight_layout()\n    plt.show()\n    \n    # Print statistics\n    print(f\"📊 Sequence Length Statistics:\")\n    print(f\"   Training - Mean: {np.mean(train_lengths):.1f}, Std: {np.std(train_lengths):.1f}\")\n    print(f\"   Test - Mean: {np.mean(test_lengths):.1f}, Std: {np.std(test_lengths):.1f}\")\n    \n    # Most common amino acids\n    print(f\"\\n🏆 Top 5 Amino Acids:\")\n    sorted_train = sorted(train_aa_comp.items(), key=lambda x: x[1], reverse=True)[:5]\n    for aa, freq in sorted_train:\n        print(f\"   {aa}: {freq:.3f}\")\n\nanalyze_sequences(data)\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-10-16T09:04:14.99022Z","iopub.execute_input":"2025-10-16T09:04:14.990545Z","iopub.status.idle":"2025-10-16T09:04:25.523747Z","shell.execute_reply.started":"2025-10-16T09:04:14.990517Z","shell.execute_reply":"2025-10-16T09:04:25.522744Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\\def analyze_go_terms(data):\n    \"\"\"Analyze GO terms and their distribution\"\"\"\n    print(\"\\n🏷️ GO TERM ANALYSIS\")\n    print(\"-\" * 40)\n    \n    train_terms = data['train_terms']\n    \n    # Terms per protein\n    terms_per_protein = train_terms.groupby('EntryID')['term'].count()\n    \n    # Create comprehensive visualization\n    fig, axes = plt.subplots(2, 3, figsize=(18, 12))\n    \n    # Terms per protein distribution\n    axes[0,0].hist(terms_per_protein, bins=50, alpha=0.7, color='teal')\n    axes[0,0].set_title('GO Terms per Protein', fontsize=14, fontweight='bold')\n    axes[0,0].set_xlabel('Number of Terms')\n    axes[0,0].set_ylabel('Number of Proteins')\n    axes[0,0].grid(True, alpha=0.3)\n    \n    # Terms by ontology (pie chart)\n    ontology_counts = train_terms['aspect'].value_counts()\n    colors = ['#ff9999', '#66b3ff', '#99ff99']\n    axes[0,1].pie(ontology_counts.values, labels=ontology_counts.index, autopct='%1.1f%%', \n                  colors=colors, startangle=90)\n    axes[0,1].set_title('GO Terms by Ontology', fontsize=14, fontweight='bold')\n    \n    # Terms by ontology (bar chart)\n    axes[0,2].bar(ontology_counts.index, ontology_counts.values, color=colors, alpha=0.7)\n    axes[0,2].set_title('GO Terms Count by Ontology', fontsize=14, fontweight='bold')\n    axes[0,2].set_ylabel('Number of Terms')\n    axes[0,2].grid(True, alpha=0.3)\n    \n    # Top terms in each ontology\n    ontologies = ['BPO', 'MFO', 'CCO']\n    colors_ont = ['lightblue', 'lightgreen', 'lightcoral']\n    \n    for i, ontology in enumerate(ontologies):\n        ontology_terms = train_terms[train_terms['aspect'] == ontology]\n        top_terms = ontology_terms['term'].value_counts().head(10)\n        \n        axes[1,i].barh(range(len(top_terms)), top_terms.values, color=colors_ont[i], alpha=0.7)\n        axes[1,i].set_title(f'Top 10 {ontology} Terms', fontsize=12, fontweight='bold')\n        axes[1,i].set_yticks(range(len(top_terms)))\n        axes[1,i].set_yticklabels([f'GO:{term}' for term in top_terms.index], fontsize=8)\n        axes[1,i].set_xlabel('Frequency')\n        axes[1,i].grid(True, alpha=0.3)\n    \n    plt.tight_layout()\n    plt.show()\n    \n    # Print detailed statistics\n    print(f\"📊 GO Term Statistics:\")\n    print(f\"   Average terms per protein: {terms_per_protein.mean():.2f}\")\n    print(f\"   Max terms per protein: {terms_per_protein.max()}\")\n    print(f\"   Min terms per protein: {terms_per_protein.min()}\")\n    print(f\"   Std of terms per protein: {terms_per_protein.std():.2f}\")\n    \n    # Label sparsity calculation\n    total_possible_annotations = len(data['train_sequences']) * train_terms['term'].nunique()\n    actual_annotations = len(train_terms)\n    sparsity = (1 - actual_annotations / total_possible_annotations) * 100\n    print(f\"   Label matrix sparsity: {sparsity:.4f}%\")\n    \n    # Most annotated proteins\n    print(f\"\\n🏆 Top 5 Most Annotated Proteins:\")\n    top_proteins = terms_per_protein.sort_values(ascending=False).head()\n    for protein, count in top_proteins.items():\n        print(f\"   {protein}: {count} terms\")\n\nanalyze_go_terms(data)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-10-16T09:04:25.524776Z","iopub.execute_input":"2025-10-16T09:04:25.5252Z","iopub.status.idle":"2025-10-16T09:04:26.717849Z","shell.execute_reply.started":"2025-10-16T09:04:25.525176Z","shell.execute_reply":"2025-10-16T09:04:26.716971Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\ndef analyze_taxonomy(data):\n    \"\"\"Analyze species distribution in the dataset\"\"\"\n    print(\"\\n🧬 TAXONOMY ANALYSIS\")\n    print(\"-\" * 40)\n    \n    taxonomy = data['train_taxonomy']\n    selected_taxa = data['selected_taxa']\n    \n    # Species distribution\n    species_counts = taxonomy['9606'].value_counts()\n    \n    # Create visualization\n    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))\n    \n    # Top 20 species\n    top_species = species_counts.head(20)\n    ax1.bar(range(len(top_species)), top_species.values, color='purple', alpha=0.7)\n    ax1.set_title('Top 20 Species by Protein Count', fontsize=14, fontweight='bold')\n    ax1.set_xlabel('Species Rank')\n    ax1.set_ylabel('Number of Proteins')\n    ax1.set_xticks(range(len(top_species)))\n    ax1.set_xticklabels([f'Taxon {taxon}' for taxon in top_species.index], rotation=45, ha='right')\n    ax1.grid(True, alpha=0.3)\n    \n    # Species distribution (log scale)\n    ax2.hist(species_counts.values, bins=50, alpha=0.7, color='orange', log=True)\n    ax2.set_title('Species Distribution (Log Scale)', fontsize=14, fontweight='bold')\n    ax2.set_xlabel('Proteins per Species')\n    ax2.set_ylabel('Number of Species (log)')\n    ax2.grid(True, alpha=0.3)\n    \n    plt.tight_layout()\n    plt.show()\n    \n    # Print statistics\n    print(f\"📊 Taxonomy Statistics:\")\n    print(f\"   Total number of species: {taxonomy['9606'].nunique():,}\")\n    print(f\"   Total number of selected taxa: {len(selected_taxa):,}\")\n    \n    print(f\"\\n🏆 Top 10 Species:\")\n    for i, (taxon, count) in enumerate(species_counts.head(10).items(), 1):\n        print(f\"   {i:2d}. Taxon {taxon}: {count:,} proteins\")\n    \n    # Coverage statistics\n    coverage_stats = species_counts.describe()\n    print(f\"\\n📈 Species Coverage:\")\n    print(f\"   Mean proteins per species: {coverage_stats['mean']:.1f}\")\n    print(f\"   Std proteins per species: {coverage_stats['std']:.1f}\")\n    print(f\"   Max proteins per species: {coverage_stats['max']:,}\")\n    print(f\"   Min proteins per species: {coverage_stats['min']}\")\n\nanalyze_taxonomy(data)\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-10-16T09:04:26.718946Z","iopub.execute_input":"2025-10-16T09:04:26.719292Z","iopub.status.idle":"2025-10-16T09:04:27.581484Z","shell.execute_reply.started":"2025-10-16T09:04:26.71926Z","shell.execute_reply":"2025-10-16T09:04:27.580569Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"data['ia_weights']","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-10-16T09:04:27.583738Z","iopub.execute_input":"2025-10-16T09:04:27.58404Z","iopub.status.idle":"2025-10-16T09:04:27.59533Z","shell.execute_reply.started":"2025-10-16T09:04:27.584018Z","shell.execute_reply":"2025-10-16T09:04:27.59441Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n# def analyze_ia_weights(data):\n#     \"\"\"Analyze information accretion weights\"\"\"\n#     print(\"\\n⚖️ INFORMATION ACCRETION ANALYSIS\")\n#     print(\"-\" * 40)\n    \n#     ia_weights = data['ia_weights']\n#     train_terms = data['train_terms']\n    \n#     # Merge with train terms to get ontology information\n#     ia_with_ontology = ia_weights.merge(\n#         train_terms[['term', 'aspect']].drop_duplicates(), \n#         on='term', \n#         how='left'\n#     )\n    \n#     # Create visualization\n#     fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(16, 12))\n    \n#     # IA distribution by ontology\n#     ontologies = ['BPO', 'MFO', 'CCO']\n#     colors = ['lightblue', 'lightgreen', 'lightcoral']\n    \n#     for i, ontology in enumerate(ontologies):\n#         ontology_ia = ia_with_ontology[ia_with_ontology['aspect'] == ontology]['ia']\n#         ax1.hist(ontology_ia, bins=50, alpha=0.6, color=colors[i], label=ontology)\n    \n#     ax1.set_title('IA Distribution by Ontology', fontsize=14, fontweight='bold')\n#     ax1.set_xlabel('Information Accretion')\n#     ax1.set_ylabel('Frequency')\n#     ax1.legend()\n#     ax1.grid(True, alpha=0.3)\n    \n#     # Box plot of IA by ontology\n#     ia_data = [ia_with_ontology[ia_with_ontology['aspect'] == ont]['ia'] for ont in ontologies]\n#     ax2.boxplot(ia_data, labels=ontologies, patch_artist=True,\n#                 boxprops=dict(facecolor='lightgray', color='black'),\n#                 medianprops=dict(color='red'))\n#     ax2.set_title('IA Distribution Box Plot', fontsize=14, fontweight='bold')\n#     ax2.set_ylabel('Information Accretion')\n#     ax2.grid(True, alpha=0.3)\n    \n#     # Cumulative distribution of IA\n#     for i, ontology in enumerate(ontologies):\n#         ontology_ia = ia_with_ontology[ia_with_ontology['aspect'] == ontology]['ia']\n#         sorted_ia = np.sort(ontology_ia)\n#         y = np.arange(1, len(sorted_ia) + 1) / len(sorted_ia)\n#         ax3.plot(sorted_ia, y, label=ontology, color=colors[i], linewidth=2)\n    \n#     ax3.set_title('Cumulative Distribution of IA', fontsize=14, fontweight='bold')\n#     ax3.set_xlabel('Information Accretion')\n#     ax3.set_ylabel('Cumulative Probability')\n#     ax3.legend()\n#     ax3.grid(True, alpha=0.3)\n    \n#     # IA vs Term Frequency\n#     term_freq = train_terms['term'].value_counts()\n#     ia_freq = ia_weights.merge(term_freq.rename('frequency'), left_on='term', right_index=True)\n    \n#     ax4.scatter(ia_freq['frequency'], ia_freq['ia'], alpha=0.5, color='purple')\n#     ax4.set_xscale('log')\n#     ax4.set_title('IA vs Term Frequency', fontsize=14, fontweight='bold')\n#     ax4.set_xlabel('Term Frequency (log scale)')\n#     ax4.set_ylabel('Information Accretion')\n#     ax4.grid(True, alpha=0.3)\n    \n#     plt.tight_layout()\n#     plt.show()\n    \n#     # Print IA statistics\n#     print(f\"📊 Information Accretion Statistics:\")\n#     for ontology in ontologies:\n#         ontology_ia = ia_with_ontology[ia_with_ontology['aspect'] == ontology]['ia']\n#         print(f\"   {ontology}:\")\n#         print(f\"      Mean IA: {ontology_ia.mean():.4f}\")\n#         print(f\"      Std IA: {ontology_ia.std():.4f}\")\n#         print(f\"      Max IA: {ontology_ia.max():.4f}\")\n#         print(f\"      Min IA: {ontology_ia.min():.4f}\")\n#         print(f\"      Median IA: {ontology_ia.median():.4f}\")\n\n# analyze_ia_weights(data)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-10-16T09:04:27.596366Z","iopub.execute_input":"2025-10-16T09:04:27.596679Z","iopub.status.idle":"2025-10-16T09:04:27.917605Z","shell.execute_reply.started":"2025-10-16T09:04:27.596655Z","shell.execute_reply":"2025-10-16T09:04:27.91628Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\ndef analyze_go_structure(data):\n    \"\"\"Analyze the GO ontology graph structure\"\"\"\n    print(\"\\n🌳 GO ONTOLOGY STRUCTURE ANALYSIS\")\n    print(\"-\" * 40)\n    \n    go_graph = data['go_graph']\n    \n    # Basic graph statistics\n    print(f\"📊 GO Graph Statistics:\")\n    print(f\"   Number of nodes (terms): {len(go_graph.nodes):,}\")\n    print(f\"   Number of edges (relationships): {len(go_graph.edges):,}\")\n    \n    # Root nodes\n    root_nodes = {\n        'BPO': 'GO:0008150',  # biological_process\n        'MFO': 'GO:0003674',  # molecular_function\n        'CCO': 'GO:0005575'   # cellular_component\n    }\n    \n    # Calculate depth and number of children for each term\n    def get_subtree_size(graph, root):\n        \"\"\"Get the size of subtree under root\"\"\"\n        descendants = nx.descendants(graph, root)\n        return len(descendants) + 1  # +1 for root itself\n    \n    print(f\"\\n🌲 Ontology Tree Sizes:\")\n    for ont_name, root_id in root_nodes.items():\n        size = get_subtree_size(go_graph, root_id)\n        print(f\"   {ont_name} ({root_id}): {size:,} terms\")\n    \n    # Degree distribution\n    in_degrees = [d for n, d in go_graph.in_degree()]\n    out_degrees = [d for n, d in go_graph.out_degree()]\n    \n    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))\n    \n    ax1.hist(in_degrees, bins=50, alpha=0.7, color='green', log=True)\n    ax1.set_title('In-Degree Distribution (Log Scale)', fontsize=14, fontweight='bold')\n    ax1.set_xlabel('In-Degree')\n    ax1.set_ylabel('Frequency (log)')\n    ax1.grid(True, alpha=0.3)\n    \n    ax2.hist(out_degrees, bins=50, alpha=0.7, color='blue', log=True)\n    ax2.set_title('Out-Degree Distribution (Log Scale)', fontsize=14, fontweight='bold')\n    ax2.set_xlabel('Out-Degree')\n    ax2.set_ylabel('Frequency (log)')\n    ax2.grid(True, alpha=0.3)\n    \n    plt.tight_layout()\n    plt.show()\n    \n    # Print degree statistics\n    print(f\"\\n📈 Degree Statistics:\")\n    print(f\"   In-degree - Mean: {np.mean(in_degrees):.2f}, Max: {max(in_degrees)}\")\n    print(f\"   Out-degree - Mean: {np.mean(out_degrees):.2f}, Max: {max(out_degrees)}\")\n    \n    # Find terms with highest degree\n    high_in_degree = sorted(go_graph.in_degree(), key=lambda x: x[1], reverse=True)[:5]\n    high_out_degree = sorted(go_graph.out_degree(), key=lambda x: x[1], reverse=True)[:5]\n    \n    print(f\"\\n🏆 Terms with Highest In-Degree:\")\n    for term, degree in high_in_degree:\n        print(f\"   {term}: {degree}\")\n    \n    print(f\"\\n🏆 Terms with Highest Out-Degree:\")\n    for term, degree in high_out_degree:\n        print(f\"   {term}: {degree}\")\n\nanalyze_go_structure(data)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-10-16T09:05:13.019608Z","iopub.execute_input":"2025-10-16T09:05:13.019973Z","iopub.status.idle":"2025-10-16T09:05:14.507523Z","shell.execute_reply.started":"2025-10-16T09:05:13.01995Z","shell.execute_reply":"2025-10-16T09:05:14.506631Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Ultra memory-efficient multi-label analysis\ndef analyze_multilabel_ultra_efficient(data):\n    \"\"\"Ultra memory-efficient multi-label analysis\"\"\"\n    print(\"\\n🎯 MULTI-LABEL CHARACTERISTICS (ULTRA EFFICIENT)\")\n    print(\"-\" * 55)\n    \n    train_terms = data['train_terms']\n    \n    # Process in chunks to avoid memory issues\n    chunk_size = 10000\n    total_chunks = (len(train_terms) + chunk_size - 1) // chunk_size\n    \n    print(\"Processing data in chunks...\")\n    \n    # Initialize counters\n    protein_term_counts = {}\n    ontology_counts = {'BPO': set(), 'MFO': set(), 'CCO': set()}\n    term_counts = {}\n    \n    # Process data in chunks\n    for i in range(total_chunks):\n        start_idx = i * chunk_size\n        end_idx = min((i + 1) * chunk_size, len(train_terms))\n        chunk = train_terms.iloc[start_idx:end_idx]\n        \n        # Update protein-term counts\n        for _, row in chunk.iterrows():\n            protein = row['EntryID']\n            term = row['term']\n            aspect = row['aspect']\n            \n            # Count terms per protein\n            if protein not in protein_term_counts:\n                protein_term_counts[protein] = 0\n            protein_term_counts[protein] += 1\n            \n            # Track proteins per ontology\n            ontology_counts[aspect].add(protein)\n            \n            # Count term frequencies\n            if term not in term_counts:\n                term_counts[term] = 0\n            term_counts[term] += 1\n        \n        if (i + 1) % 10 == 0:\n            print(f\"Processed {end_idx:,} rows...\")\n    \n    # Convert to series for analysis\n    terms_per_protein = pd.Series(protein_term_counts)\n    \n    # Calculate statistics\n    label_cardinality = terms_per_protein.mean()\n    total_terms = len(term_counts)\n    label_density = label_cardinality / total_terms\n    \n    # Create simplified visualization\n    fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(18, 5))\n    \n    # Histogram of terms per protein (sampled)\n    if len(terms_per_protein) > 50000:\n        sampled = terms_per_protein.sample(n=30000, random_state=42)\n        ax1.hist(sampled, bins=30, alpha=0.7, color='red')\n        ax1.set_title(f'Terms per Protein\\n(30K sample of {len(terms_per_protein):,})', fontweight='bold')\n    else:\n        ax1.hist(terms_per_protein, bins=30, alpha=0.7, color='red')\n        ax1.set_title('Terms per Protein', fontweight='bold')\n    ax1.set_xlabel('Number of Terms')\n    ax1.set_ylabel('Number of Proteins')\n    ax1.grid(True, alpha=0.3)\n    \n    # Ontology distribution\n    ontology_sizes = {ont: len(proteins) for ont, proteins in ontology_counts.items()}\n    ax2.bar(ontology_sizes.keys(), ontology_sizes.values(), \n            color=['lightblue', 'lightgreen', 'lightcoral'], alpha=0.7)\n    ax2.set_title('Proteins per Ontology', fontweight='bold')\n    ax2.set_ylabel('Number of Proteins')\n    ax2.grid(True, alpha=0.3)\n    \n    # Term frequency distribution (log scale)\n    term_freq = pd.Series(term_counts)\n    ax3.hist(term_freq.values, bins=50, alpha=0.7, color='purple', log=True)\n    ax3.set_title('GO Term Frequency Distribution\\n(Log Scale)', fontweight='bold')\n    ax3.set_xlabel('Term Frequency')\n    ax3.set_ylabel('Number of Terms (log)')\n    ax3.grid(True, alpha=0.3)\n    \n    plt.tight_layout()\n    plt.show()\n    \n    # Print statistics\n    print(f\"📊 Multi-label Statistics:\")\n    print(f\"   Total proteins: {len(protein_term_counts):,}\")\n    print(f\"   Total GO terms: {total_terms:,}\")\n    print(f\"   Total annotations: {len(train_terms):,}\")\n    print(f\"   Label cardinality: {label_cardinality:.2f} terms/protein\")\n    print(f\"   Label density: {label_density:.6f}\")\n    \n    print(f\"\\n🏷️ Ontology Coverage:\")\n    for ont, proteins in ontology_counts.items():\n        coverage = len(proteins) / len(protein_term_counts) * 100\n        print(f\"   {ont}: {len(proteins):,} proteins ({coverage:.1f}%)\")\n    \n    # Calculate overlaps\n    bpo_set = ontology_counts['BPO']\n    mfo_set = ontology_counts['MFO']\n    cco_set = ontology_counts['CCO']\n    \n    all_three = len(bpo_set & mfo_set & cco_set)\n    print(f\"   All three ontologies: {all_three:,} proteins ({all_three/len(protein_term_counts)*100:.1f}%)\")\n    \n    print(f\"\\n📈 Term Frequency Statistics:\")\n    print(f\"   Most frequent term: {term_freq.max()} occurrences\")\n    print(f\"   Average term frequency: {term_freq.mean():.1f}\")\n    print(f\"   Terms appearing only once: {(term_freq == 1).sum():,}\")\n    \n    # Show top terms\n    top_terms = term_freq.nlargest(10)\n    print(f\"\\n🏆 Top 10 Most Frequent GO Terms:\")\n    for term, count in top_terms.items():\n        print(f\"   {term}: {count:,} proteins\")\n\n# Use this if the optimized version still causes memory issues\nanalyze_multilabel_ultra_efficient(data)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-10-16T09:05:17.855022Z","iopub.execute_input":"2025-10-16T09:05:17.855342Z","iopub.status.idle":"2025-10-16T09:05:17.892467Z","shell.execute_reply.started":"2025-10-16T09:05:17.855321Z","shell.execute_reply":"2025-10-16T09:05:17.890927Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Create interactive visualizations\ndef create_interactive_plots(data):\n    \"\"\"Create interactive plots using Plotly\"\"\"\n    print(\"\\n📊 INTERACTIVE VISUALIZATIONS\")\n    print(\"-\" * 40)\n    \n    train_terms = data['train_terms']\n    \n    # Interactive terms per protein distribution\n    terms_per_protein = train_terms.groupby('EntryID')['term'].count()\n    \n    fig1 = px.histogram(terms_per_protein, nbins=50, \n                       title='Interactive: GO Terms per Protein Distribution',\n                       labels={'value': 'Number of Terms', 'count': 'Number of Proteins'},\n                       opacity=0.7)\n    fig1.show()\n    \n    # Interactive ontology comparison\n    ontology_stats = train_terms.groupby('aspect').agg({\n        'term': ['count', 'nunique'],\n        'EntryID': 'nunique'\n    }).round(2)\n    ontology_stats.columns = ['Total Annotations', 'Unique Terms', 'Unique Proteins']\n    ontology_stats = ontology_stats.reset_index()\n    \n    fig2 = px.bar(ontology_stats, x='aspect', y='Unique Terms',\n                 title='Interactive: Unique GO Terms by Ontology',\n                 color='aspect',\n                 text='Unique Terms')\n    fig2.show()\n    \n    # Sequence length interactive plot\n    train_lengths = [len(seq) for seq in data['train_sequences'].values()]\n    test_lengths = [len(seq) for seq in data['test_sequences'].values()]\n    \n    fig3 = go.Figure()\n    fig3.add_trace(go.Box(y=train_lengths, name='Training Set', boxpoints='outliers'))\n    fig3.add_trace(go.Box(y=test_lengths, name='Test Set', boxpoints='outliers'))\n    fig3.update_layout(title='Interactive: Sequence Length Distribution Comparison',\n                      yaxis_title='Sequence Length')\n    fig3.show()\n    \n    print(\"✅ Interactive plots created successfully!\")\n\n# Uncomment to run interactive plots (may be heavy for large datasets)\ncreate_interactive_plots(data)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-10-16T08:50:38.413581Z","iopub.execute_input":"2025-10-16T08:50:38.41531Z","iopub.status.idle":"2025-10-16T08:50:43.935675Z","shell.execute_reply.started":"2025-10-16T08:50:38.41527Z","shell.execute_reply":"2025-10-16T08:50:43.934544Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Generate key insights summary\ndef generate_insights_summary(data):\n    \"\"\"Generate a summary of key insights from the EDA\"\"\"\n    print(\"\\n\" + \"=\"*60)\n    print(\"🔍 KEY INSIGHTS SUMMARY\")\n    print(\"=\"*60)\n    \n    train_terms = data['train_terms']\n    train_seqs = list(data['train_sequences'].values())\n    test_seqs = list(data['test_sequences'].values())\n    \n    # Calculate key metrics\n    terms_per_protein = train_terms.groupby('EntryID')['term'].count()\n    avg_terms_per_protein = terms_per_protein.mean()\n    \n    train_avg_len = np.mean([len(seq) for seq in train_seqs])\n    test_avg_len = np.mean([len(seq) for seq in test_seqs])\n    \n    # Unique terms by ontology\n    bpo_terms = train_terms[train_terms['aspect'] == 'BPO']['term'].nunique()\n    mfo_terms = train_terms[train_terms['aspect'] == 'MFO']['term'].nunique()\n    cco_terms = train_terms[train_terms['aspect'] == 'CCO']['term'].nunique()\n    \n    print(\"\\n📈 COMPETITION SCALE:\")\n    print(f\"   • {len(data['train_sequences']):,} training proteins\")\n    print(f\"   • {len(data['test_sequences']):,} test proteins\")\n    print(f\"   • {len(train_terms):,} total GO term annotations\")\n    print(f\"   • {train_terms['term'].nunique():,} unique GO terms\")\n    \n    print(\"\\n🎯 MULTI-LABEL COMPLEXITY:\")\n    print(f\"   • Average {avg_terms_per_protein:.1f} terms per protein\")\n    print(f\"   • Extremely sparse label matrix (>99.9% sparse)\")\n    print(f\"   • Hierarchical relationships between terms\")\n    \n    print(\"\\n🧬 ONTOLOGY DISTRIBUTION:\")\n    print(f\"   • Biological Process (BPO): {bpo_terms:,} terms\")\n    print(f\"   • Molecular Function (MFO): {mfo_terms:,} terms\")\n    print(f\"   • Cellular Component (CCO): {cco_terms:,} terms\")\n    \n    print(\"\\n📏 SEQUENCE CHARACTERISTICS:\")\n    print(f\"   • Training avg length: {train_avg_len:.1f} amino acids\")\n    print(f\"   • Test avg length: {test_avg_len:.1f} amino acids\")\n    print(f\"   • Wide range of sequence lengths (tens to thousands)\")\n    \n    print(\"\\n⚖️ EVALUATION COMPLEXITY:\")\n    print(f\"   • Weighted by Information Accretion (IA)\")\n    print(f\"   • Hierarchical precision/recall\")\n    print(f\"   • Three separate subontology evaluations\")\n    \n    print(\"\\n🚀 RECOMMENDATIONS FOR MODELING:\")\n    print(f\"   • Use protein language models (ESM, ProtBERT)\")\n    print(f\"   • Implement hierarchical multi-label classification\")\n    print(f\"   • Handle extreme class imbalance\")\n    print(f\"   • Consider taxonomic information\")\n    print(f\"   • Use IA weights in loss function\")\n    \n    print(\"\\n⚠️  CHALLENGES:\")\n    print(f\"   • Extreme multi-label classification\")\n    print(f\"   • Hierarchical label relationships\")\n    print(f\"   • Sparse and imbalanced annotations\")\n    print(f\"   • Computational complexity\")\n    print(f\"   • Prospective evaluation (future test set)\")\n\ngenerate_insights_summary(data)\n\n# %% [markdown]\n# ## 10. Data Quality Checks\n\n# %% [code]\n# Data quality validation\ndef perform_quality_checks(data):\n    \"\"\"Perform data quality checks and validation\"\"\"\n    print(\"\\n🔍 DATA QUALITY CHECKS\")\n    print(\"-\" * 40)\n    \n    # Check for missing values\n    print(\"1. Missing Values Check:\")\n    train_terms_missing = data['train_terms'].isnull().sum()\n    taxonomy_missing = data['train_taxonomy'].isnull().sum()\n    \n    print(f\"   Train terms - Missing values: {train_terms_missing.sum()}\")\n    print(f\"   Taxonomy - Missing values: {taxonomy_missing.sum()}\")\n    \n    # Check sequence validity\n    print(\"\\n2. Sequence Validity Check:\")\n    valid_aa = set('ACDEFGHIKLMNPQRSTVWY')\n    \n    def check_sequence_validity(sequences):\n        invalid_chars = []\n        for seq in sequences:\n            invalid_in_seq = set(seq) - valid_aa\n            invalid_chars.extend(invalid_in_seq)\n        return set(invalid_chars)\n    \n    train_invalid = check_sequence_validity(data['train_sequences'].values())\n    test_invalid = check_sequence_validity(data['test_sequences'].values())\n    \n    print(f\"   Training - Invalid characters: {train_invalid if train_invalid else 'None'}\")\n    print(f\"   Test - Invalid characters: {test_invalid if test_invalid else 'None'}\")\n    \n    # Check ID consistency\n    print(\"\\n3. ID Consistency Check:\")\n    train_proteins = set(data['train_sequences'].keys())\n    term_proteins = set(data['train_terms']['EntryID'].unique())\n    taxonomy_proteins = set(data['train_taxonomy']['EntryID'])\n    \n    print(f\"   Proteins in sequences: {len(train_proteins):,}\")\n    print(f\"   Proteins in terms: {len(term_proteins):,}\")\n    print(f\"   Proteins in taxonomy: {len(taxonomy_proteins):,}\")\n    \n    missing_in_terms = train_proteins - term_proteins\n    missing_in_taxonomy = train_proteins - taxonomy_proteins\n    \n    print(f\"   Proteins missing term annotations: {len(missing_in_terms)}\")\n    print(f\"   Proteins missing taxonomy: {len(missing_in_taxonomy)}\")\n    \n    # Check GO term validity\n    print(\"\\n4. GO Term Validity Check:\")\n    valid_terms = set(data['go_graph'].nodes())\n    used_terms = set(data['train_terms']['term'].unique())\n    \n    invalid_terms = used_terms - valid_terms\n    print(f\"   Used GO terms: {len(used_terms):,}\")\n    print(f\"   Valid GO terms in ontology: {len(valid_terms):,}\")\n    print(f\"   Invalid GO terms: {len(invalid_terms)}\")\n    \n    if invalid_terms:\n        print(f\"   Example invalid terms: {list(invalid_terms)[:5]}\")\n    \n    print(\"\\n✅ Data quality checks completed!\")\n\nperform_quality_checks(data)\n\n# This comprehensive EDA reveals that CAFA 6 is an extremely challenging multi-label classification problem with:\n# \n# - **Large scale**: Hundreds of thousands of proteins and GO terms\n# - **High complexity**: Hierarchical, multi-ontology predictions\n# - **Extreme sparsity**: Very few positive labels per protein\n# - **Complex evaluation**: Information accretion weighted metrics\n# - **Prospective test set**: Future annotations as ground truth\n\n# Final summary cell\nprint(\"\\n\" + \"=\"*70)\nprint(\"🎉 CAFA 6 EDA COMPLETED SUCCESSFULLY!\")\nprint(\"=\"*70)\nprint(\"\\nNext steps:\")\nprint(\"1. Implement baseline models (sequence similarity, ESM embeddings)\")\nprint(\"2. Develop hierarchical multi-label classification approaches\")\nprint(\"3. Incorporate taxonomic and evolutionary information\")\nprint(\"4. Handle class imbalance with appropriate loss functions\")\nprint(\"5. Validate using the proper evaluation metrics\")\nprint(\"\\nGood luck with the competition! 🚀\")\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-10-16T09:05:48.993511Z","iopub.execute_input":"2025-10-16T09:05:48.993895Z","iopub.status.idle":"2025-10-16T09:05:51.541458Z","shell.execute_reply.started":"2025-10-16T09:05:48.993866Z","shell.execute_reply":"2025-10-16T09:05:51.540051Z"}},"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},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}