{"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":"nvidiaTeslaT4","dataSources":[{"sourceId":116062,"databundleVersionId":14084779,"sourceType":"competition"}],"dockerImageVersionId":31154,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# šŸ§¬ CAFA 6: Complete Exploratory Data Analysis\n- Author: Dr.Hany Ghazal (PhD)\n- Competition: CAFA 6 - Protein Function Prediction\n- Last Updated: October 2025","metadata":{}},{"cell_type":"code","source":"!pip install biopython obonet","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T19:58:58.568580Z","iopub.execute_input":"2025-10-19T19:58:58.569100Z","iopub.status.idle":"2025-10-19T19:59:04.013898Z","shell.execute_reply.started":"2025-10-19T19:58:58.569071Z","shell.execute_reply":"2025-10-19T19:59:04.013042Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nfrom collections import Counter, defaultdict\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom Bio import SeqIO\nimport networkx as nx\nfrom typing import Dict, List, Tuple, Set\nimport warnings\nwarnings.filterwarnings('ignore')\n\n# Set plotting style\nsns.set_style(\"whitegrid\")\nplt.rcParams['figure.figsize'] = (12, 6)\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n for filename in filenames:\n print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-10-19T19:59:04.015496Z","iopub.execute_input":"2025-10-19T19:59:04.015773Z","iopub.status.idle":"2025-10-19T19:59:05.296041Z","shell.execute_reply.started":"2025-10-19T19:59:04.015748Z","shell.execute_reply":"2025-10-19T19:59:05.295426Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# First, we need to write functions for data loading","metadata":{}},{"cell_type":"code","source":"data_dir = \"/kaggle/input/cafa-6-protein-function-prediction/\"","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T19:59:05.296744Z","iopub.execute_input":"2025-10-19T19:59:05.297153Z","iopub.status.idle":"2025-10-19T19:59:05.300798Z","shell.execute_reply.started":"2025-10-19T19:59:05.297097Z","shell.execute_reply":"2025-10-19T19:59:05.300017Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# File paths\ntrain_seq_file = f\"{data_dir}Train/train_sequences.fasta\"\ntest_seq_file = f\"{data_dir}Test/testsuperset.fasta\"\nterms_file = f\"{data_dir}Train/train_terms.tsv\"\ntaxonomy_file = f\"{data_dir}Train/train_taxonomy.tsv\"\nia_file = f\"{data_dir}IA.tsv\"\nobo_file = f\"{data_dir}Train/go-basic.obo\"","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T19:59:05.302612Z","iopub.execute_input":"2025-10-19T19:59:05.302815Z","iopub.status.idle":"2025-10-19T19:59:05.320216Z","shell.execute_reply.started":"2025-10-19T19:59:05.302799Z","shell.execute_reply":"2025-10-19T19:59:05.319541Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================\n# 1. DATA LOADING FUNCTIONS\n# ============================================\n\ndef load_fasta_sequences(fasta_file: str) -> Dict[str, str]:\n \"\"\"\n Load protein sequences from FASTA file.\n \n Args:\n fasta_file: Path to FASTA file\n \n Returns:\n Dictionary mapping protein ID to sequence\n \"\"\"\n sequences = {}\n \n try:\n for record in SeqIO.parse(fasta_file, \"fasta\"):\n # Extract UniProt ID from header (format: sp|ID|NAME)\n header_parts = record.id.split('|')\n if len(header_parts) >= 2:\n protein_id = header_parts[1]\n else:\n protein_id = record.id\n \n sequences[protein_id] = str(record.seq)\n \n print(f\"Loaded {len(sequences):,} sequences from {fasta_file}\")\n \n except FileNotFoundError:\n print(f\"Error: File {fasta_file} not found\")\n # Return sample data for demo\n sequences = {\n 'P9WHI7': 'MTKPTQVLVRLEQVMAEGDLRTLVPLLKAGQYRDDVTFGGPTVHPGVRVEGHTL',\n 'P04637': 'MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFT'\n }\n print(\"Using sample data for demonstration\")\n \n return sequences\n\n\ndef load_go_annotations(terms_file: str) -> pd.DataFrame:\n \"\"\"\n Load GO term annotations.\n \n Args:\n terms_file: Path to train_terms.tsv\n \n Returns:\n DataFrame with columns: protein_id, go_term, ontology\n \"\"\"\n try:\n df = pd.read_csv(terms_file, sep='\\t', \n names=['protein_id', 'go_term', 'ontology'])\n print(f\"Loaded {len(df):,} GO annotations for {df['protein_id'].nunique():,} proteins\")\n \n except FileNotFoundError:\n print(f\"Error: File {terms_file} not found\")\n # Sample data for demo\n df = pd.DataFrame({\n 'protein_id': ['P9WHI7', 'P9WHI7', 'P04637', 'P04637'],\n 'go_term': ['GO:0009274', 'GO:0071944', 'GO:1990837', 'GO:0031625'],\n 'ontology': ['BPO', 'CCO', 'MFO', 'MFO']\n })\n print(\"Using sample data for demonstration\")\n \n return df\n\n\ndef load_taxonomy(taxonomy_file: str) -> pd.DataFrame:\n \"\"\"\n Load taxonomy information.\n \n Args:\n taxonomy_file: Path to train_taxonomy.tsv\n \n Returns:\n DataFrame with protein_id and taxon_id\n \"\"\"\n try:\n df = pd.read_csv(taxonomy_file, sep='\\t', \n names=['protein_id', 'taxon_id'])\n print(f\"Loaded taxonomy for {len(df):,} proteins across {df['taxon_id'].nunique()} species\")\n \n except FileNotFoundError:\n print(f\"Error: File {taxonomy_file} not found\")\n # Sample data\n df = pd.DataFrame({\n 'protein_id': ['P9WHI7', 'P04637'],\n 'taxon_id': [1773, 9606]\n })\n print(\"Using sample data for demonstration\")\n \n return df\n\n\ndef load_information_accretion(ia_file: str) -> pd.DataFrame:\n \"\"\"\n Load information accretion weights.\n \n Args:\n ia_file: Path to IA.tsv\n \n Returns:\n DataFrame with go_term and ia_weight\n \"\"\"\n try:\n df = pd.read_csv(ia_file, sep='\\t', names=['go_term', 'ia_weight'])\n print(f\"Loaded IA weights for {len(df):,} GO terms\")\n \n except FileNotFoundError:\n print(f\"Error: File {ia_file} not found\")\n # Sample data\n df = pd.DataFrame({\n 'go_term': ['GO:0009274', 'GO:0071944', 'GO:1990837', 'GO:0031625'],\n 'ia_weight': [0.5, 0.3, 0.8, 0.6]\n })\n print(\"Using sample data for demonstration\")\n \n return df\n\n\ndef parse_obo_file(obo_file: str) -> nx.DiGraph:\n \"\"\"\n Parse GO ontology OBO file and create directed graph.\n \n Args:\n obo_file: Path to go-basic.obo\n \n Returns:\n NetworkX directed graph of GO ontology\n \"\"\"\n try:\n import obonet\n graph = obonet.read_obo(obo_file)\n print(f\"Loaded GO ontology: {graph.number_of_nodes():,} terms, {graph.number_of_edges():,} relationships\")\n return graph\n \n except (FileNotFoundError, ImportError):\n print(f\"Error loading {obo_file} or obonet not installed\")\n print(\"Creating sample graph for demonstration\")\n \n # Create sample GO graph\n G = nx.DiGraph()\n sample_terms = [\n ('GO:0008150', 'biological_process'),\n ('GO:0009987', 'cellular process'),\n ('GO:0009274', 'peptidoglycan-based cell wall'),\n ('GO:0003674', 'molecular_function'),\n ('GO:0005575', 'cellular_component'),\n ]\n \n for term_id, name in sample_terms:\n G.add_node(term_id, name=name)\n \n # Add some edges\n G.add_edge('GO:0008150', 'GO:0009987')\n G.add_edge('GO:0009987', 'GO:0009274')\n \n return G","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T19:59:05.320902Z","iopub.execute_input":"2025-10-19T19:59:05.321139Z","iopub.status.idle":"2025-10-19T19:59:05.339085Z","shell.execute_reply.started":"2025-10-19T19:59:05.321103Z","shell.execute_reply":"2025-10-19T19:59:05.338456Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 1. Load all data\nprint(\"\\n\" + \"=\"*60)\nprint(\"1. LOADING DATA\")\nprint(\"=\"*60)\n\ntrain_sequences = load_fasta_sequences(train_seq_file)\ntest_sequences = load_fasta_sequences(test_seq_file)\nannotations = load_go_annotations(terms_file)\ntaxonomy = load_taxonomy(taxonomy_file)\nia_weights = load_information_accretion(ia_file)\ngo_graph = parse_obo_file(obo_file)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T19:59:05.339855Z","iopub.execute_input":"2025-10-19T19:59:05.340341Z","iopub.status.idle":"2025-10-19T19:59:13.024229Z","shell.execute_reply.started":"2025-10-19T19:59:05.340323Z","shell.execute_reply":"2025-10-19T19:59:13.023588Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Sequence Analysis Functions","metadata":{}},{"cell_type":"code","source":"def analyze_sequence_statistics(sequences: Dict[str, str]) -> pd.DataFrame:\n \"\"\"\n Calculate basic statistics for protein sequences.\n \n Args:\n sequences: Dictionary mapping protein_id to sequence\n \n Returns:\n DataFrame with sequence statistics\n \"\"\"\n stats = []\n \n for protein_id, seq in sequences.items():\n stats.append({\n 'protein_id': protein_id,\n 'length': len(seq),\n 'molecular_weight': calculate_molecular_weight(seq),\n 'num_A': seq.count('A'),\n 'num_C': seq.count('C'),\n 'num_D': seq.count('D'),\n 'num_E': seq.count('E'),\n 'num_F': seq.count('F'),\n 'num_G': seq.count('G'),\n 'num_H': seq.count('H'),\n 'num_I': seq.count('I'),\n 'num_K': seq.count('K'),\n 'num_L': seq.count('L'),\n 'num_M': seq.count('M'),\n 'num_N': seq.count('N'),\n 'num_P': seq.count('P'),\n 'num_Q': seq.count('Q'),\n 'num_R': seq.count('R'),\n 'num_S': seq.count('S'),\n 'num_T': seq.count('T'),\n 'num_V': seq.count('V'),\n 'num_W': seq.count('W'),\n 'num_Y': seq.count('Y'),\n })\n \n df = pd.DataFrame(stats)\n \n print(\"\\nSequence Statistics Summary:\")\n print(f\"Total proteins: {len(df):,}\")\n print(f\"Length - Mean: {df['length'].mean():.1f}, Median: {df['length'].median():.1f}\")\n print(f\"Length - Min: {df['length'].min()}, Max: {df['length'].max()}\")\n print(f\"Molecular Weight - Mean: {df['molecular_weight'].mean():.1f} Da\")\n \n return df\n\ndef calculate_molecular_weight(sequence: str) -> float:\n \"\"\"Calculate approximate molecular weight of protein sequence.\"\"\"\n aa_weights = {\n 'A': 89.1, 'C': 121.2, 'D': 133.1, 'E': 147.1, 'F': 165.2,\n 'G': 75.1, 'H': 155.2, 'I': 131.2, 'K': 146.2, 'L': 131.2,\n 'M': 149.2, 'N': 132.1, 'P': 115.1, 'Q': 146.2, 'R': 174.2,\n 'S': 105.1, 'T': 119.1, 'V': 117.1, 'W': 204.2, 'Y': 181.2\n }\n \n weight = sum(aa_weights.get(aa, 0) for aa in sequence)\n return weight - (len(sequence) - 1) * 18.015 # Subtract water molecules\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T19:59:13.024859Z","iopub.execute_input":"2025-10-19T19:59:13.025036Z","iopub.status.idle":"2025-10-19T19:59:13.033466Z","shell.execute_reply.started":"2025-10-19T19:59:13.025021Z","shell.execute_reply":"2025-10-19T19:59:13.032805Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport seaborn as sns\nimport pandas as pd\nimport numpy as np\nfrom collections import Counter, defaultdict\nfrom typing import Dict, List, Tuple\nimport networkx as nx\nimport warnings\nwarnings.filterwarnings('ignore')\n\n# Set global plotting style\nsns.set_style(\"whitegrid\")\nsns.set_palette(\"husl\")\nplt.rcParams['figure.figsize'] = (14, 6)\nplt.rcParams['font.size'] = 10\nplt.rcParams['axes.titlesize'] = 14\nplt.rcParams['axes.labelsize'] = 12\n\n\n# ============================================================================\n# 1. SEQUENCE VISUALIZATIONS\n# ============================================================================\n\ndef plot_sequence_length_distribution(sequences: Dict[str, str], \n title: str = \"Sequence Length Distribution\",\n save_path: str = None):\n \"\"\"\n Plot distribution of protein sequence lengths.\n \n Args:\n sequences: Dictionary mapping protein_id to sequence\n title: Plot title\n save_path: Optional path to save figure\n \"\"\"\n lengths = [len(seq) for seq in sequences.values()]\n \n fig, axes = plt.subplots(1, 2, figsize=(16, 6))\n \n # Histogram\n axes[0].hist(lengths, bins=50, color='steelblue', alpha=0.7, edgecolor='black')\n axes[0].axvline(np.mean(lengths), color='red', linestyle='--', \n linewidth=2, label=f'Mean: {np.mean(lengths):.0f}')\n axes[0].axvline(np.median(lengths), color='orange', linestyle='--', \n linewidth=2, label=f'Median: {np.median(lengths):.0f}')\n axes[0].set_xlabel('Sequence Length (amino acids)', fontsize=12)\n axes[0].set_ylabel('Frequency', fontsize=12)\n axes[0].set_title(title, fontsize=14, fontweight='bold')\n axes[0].legend()\n axes[0].grid(alpha=0.3)\n \n # Box plot\n axes[1].boxplot(lengths, vert=True, patch_artist=True,\n boxprops=dict(facecolor='lightblue', alpha=0.7),\n medianprops=dict(color='red', linewidth=2))\n axes[1].set_ylabel('Sequence Length (amino acids)', fontsize=12)\n axes[1].set_title('Box Plot', fontsize=14, fontweight='bold')\n axes[1].grid(alpha=0.3)\n \n plt.tight_layout()\n \n if save_path:\n plt.savefig(save_path, dpi=300, bbox_inches='tight')\n \n plt.show()\n \n # Print statistics\n print(f\"\\n{'='*60}\")\n print(f\"SEQUENCE LENGTH STATISTICS\")\n print(f\"{'='*60}\")\n print(f\"Total sequences: {len(lengths):,}\")\n print(f\"Mean length: {np.mean(lengths):.2f} AA\")\n print(f\"Median length: {np.median(lengths):.2f} AA\")\n print(f\"Std deviation: {np.std(lengths):.2f} AA\")\n print(f\"Min length: {min(lengths)} AA\")\n print(f\"Max length: {max(lengths):,} AA\")\n print(f\"25th percentile: {np.percentile(lengths, 25):.0f} AA\")\n print(f\"75th percentile: {np.percentile(lengths, 75):.0f} AA\")\n print(f\"{'='*60}\\n\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T19:59:13.034219Z","iopub.execute_input":"2025-10-19T19:59:13.034459Z","iopub.status.idle":"2025-10-19T19:59:13.053160Z","shell.execute_reply.started":"2025-10-19T19:59:13.034435Z","shell.execute_reply":"2025-10-19T19:59:13.052537Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"analyze_sequence_statistics(train_sequences)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T19:59:13.053966Z","iopub.execute_input":"2025-10-19T19:59:13.054189Z","iopub.status.idle":"2025-10-19T19:59:17.791044Z","shell.execute_reply.started":"2025-10-19T19:59:13.054170Z","shell.execute_reply":"2025-10-19T19:59:17.790505Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plot_sequence_length_distribution(train_sequences, save_path=\"/kaggle/working/train_sequence_length_distribution.png\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T19:59:17.794742Z","iopub.execute_input":"2025-10-19T19:59:17.795098Z","iopub.status.idle":"2025-10-19T19:59:19.729215Z","shell.execute_reply.started":"2025-10-19T19:59:17.795071Z","shell.execute_reply":"2025-10-19T19:59:19.728373Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"analyze_sequence_statistics(test_sequences)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T19:59:19.730061Z","iopub.execute_input":"2025-10-19T19:59:19.730440Z","iopub.status.idle":"2025-10-19T19:59:31.035703Z","shell.execute_reply.started":"2025-10-19T19:59:19.730416Z","shell.execute_reply":"2025-10-19T19:59:31.035072Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plot_sequence_length_distribution(test_sequences, save_path=\"/kaggle/working/test_sequence_length_distribution.png\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T19:59:31.036463Z","iopub.execute_input":"2025-10-19T19:59:31.036718Z","iopub.status.idle":"2025-10-19T19:59:33.565484Z","shell.execute_reply.started":"2025-10-19T19:59:31.036690Z","shell.execute_reply":"2025-10-19T19:59:33.564881Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def analyze_amino_acid_composition(sequences: Dict[str, str]) -> pd.DataFrame:\n \"\"\"\n Calculate amino acid composition statistics.\n \n Returns:\n DataFrame with AA frequencies\n \"\"\"\n total_counts = Counter()\n \n for seq in sequences.values():\n total_counts.update(seq)\n \n total_aa = sum(total_counts.values())\n \n composition = pd.DataFrame([\n {'amino_acid': aa, 'count': count, 'frequency': count/total_aa}\n for aa, count in sorted(total_counts.items())\n ])\n \n print(\"\\nAmino Acid Composition (Top 10):\")\n print(composition.nlargest(10, 'frequency')[['amino_acid', 'frequency']])\n \n return composition\n\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T19:59:33.566158Z","iopub.execute_input":"2025-10-19T19:59:33.566420Z","iopub.status.idle":"2025-10-19T19:59:33.571675Z","shell.execute_reply.started":"2025-10-19T19:59:33.566399Z","shell.execute_reply":"2025-10-19T19:59:33.571103Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"analyze_amino_acid_composition(train_sequences)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T19:59:33.572330Z","iopub.execute_input":"2025-10-19T19:59:33.572549Z","iopub.status.idle":"2025-10-19T19:59:35.783571Z","shell.execute_reply.started":"2025-10-19T19:59:33.572525Z","shell.execute_reply":"2025-10-19T19:59:35.782922Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"analyze_amino_acid_composition(test_sequences)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T19:59:35.784319Z","iopub.execute_input":"2025-10-19T19:59:35.784585Z","iopub.status.idle":"2025-10-19T19:59:40.689618Z","shell.execute_reply.started":"2025-10-19T19:59:35.784568Z","shell.execute_reply":"2025-10-19T19:59:40.689019Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def plot_amino_acid_composition(sequences: Dict[str, str], \n top_n: int = 20,\n save_path: str = None):\n \"\"\"\n Plot amino acid composition across all sequences.\n \n Args:\n sequences: Dictionary mapping protein_id to sequence\n top_n: Number of amino acids to display (default 20 for all)\n save_path: Optional path to save figure\n \"\"\"\n # Count amino acids\n total_counts = Counter()\n for seq in sequences.values():\n total_counts.update(seq)\n \n total_aa = sum(total_counts.values())\n \n # Create DataFrame\n aa_data = pd.DataFrame([\n {'AA': aa, 'Count': count, 'Frequency': count/total_aa*100}\n for aa, count in sorted(total_counts.items())\n ]).sort_values('Frequency', ascending=False).head(top_n)\n \n # Create figure\n fig, axes = plt.subplots(1, 2, figsize=(16, 6))\n \n # Bar chart\n colors = plt.cm.viridis(np.linspace(0, 1, len(aa_data)))\n axes[0].bar(aa_data['AA'], aa_data['Frequency'], \n color=colors, alpha=0.8, edgecolor='black')\n axes[0].set_xlabel('Amino Acid', fontsize=12)\n axes[0].set_ylabel('Frequency (%)', fontsize=12)\n axes[0].set_title('Amino Acid Composition', fontsize=14, fontweight='bold')\n axes[0].grid(axis='y', alpha=0.3)\n \n # Add percentage labels on bars\n for i, (aa, freq) in enumerate(zip(aa_data['AA'], aa_data['Frequency'])):\n axes[0].text(i, freq + 0.2, f'{freq:.1f}%', \n ha='center', va='bottom', fontsize=8)\n \n # Pie chart (top 10)\n top10 = aa_data.head(10)\n axes[1].pie(top10['Frequency'], labels=top10['AA'], autopct='%1.1f%%',\n colors=colors[:10], startangle=90)\n axes[1].set_title('Top 10 Amino Acids', fontsize=14, fontweight='bold')\n \n plt.tight_layout()\n \n if save_path:\n plt.savefig(save_path, dpi=300, bbox_inches='tight')\n \n plt.show()\n \n # Print statistics\n print(f\"\\n{'='*60}\")\n print(f\"AMINO ACID COMPOSITION\")\n print(f\"{'='*60}\")\n print(f\"Total amino acids: {total_aa:,}\")\n print(f\"\\nTop 10 Most Frequent:\")\n for i, row in aa_data.head(10).iterrows():\n print(f\" {row['AA']}: {row['Count']:>12,} ({row['Frequency']:>5.2f}%)\")\n print(f\"{'='*60}\\n\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T19:59:40.690276Z","iopub.execute_input":"2025-10-19T19:59:40.690482Z","iopub.status.idle":"2025-10-19T19:59:40.699597Z","shell.execute_reply.started":"2025-10-19T19:59:40.690462Z","shell.execute_reply":"2025-10-19T19:59:40.699070Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plot_amino_acid_composition(train_sequences, save_path=\"/kaggle/working/train_amino_acid_composition.png\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T19:59:40.700332Z","iopub.execute_input":"2025-10-19T19:59:40.700664Z","iopub.status.idle":"2025-10-19T19:59:44.535557Z","shell.execute_reply.started":"2025-10-19T19:59:40.700639Z","shell.execute_reply":"2025-10-19T19:59:44.534818Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plot_amino_acid_composition(test_sequences, save_path=\"/kaggle/working/test_amino_acid_composition.png\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T19:59:44.536320Z","iopub.execute_input":"2025-10-19T19:59:44.536518Z","iopub.status.idle":"2025-10-19T19:59:51.438905Z","shell.execute_reply.started":"2025-10-19T19:59:44.536503Z","shell.execute_reply":"2025-10-19T19:59:51.438184Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def plot_sequence_comparison(train_seqs: Dict[str, str], \n test_seqs: Dict[str, str],\n save_path: str = None):\n \"\"\"\n Compare train and test sequence distributions.\n \n Args:\n train_seqs: Training sequences\n test_seqs: Test sequences\n save_path: Optional path to save figure\n \"\"\"\n train_lengths = [len(seq) for seq in train_seqs.values()]\n test_lengths = [len(seq) for seq in test_seqs.values()]\n \n fig, axes = plt.subplots(2, 2, figsize=(16, 12))\n \n # Overlapping histograms\n axes[0, 0].hist(train_lengths, bins=50, color='blue', alpha=0.5, \n label=f'Train (n={len(train_lengths):,})', edgecolor='black')\n axes[0, 0].hist(test_lengths, bins=50, color='red', alpha=0.5, \n label=f'Test (n={len(test_lengths):,})', edgecolor='black')\n axes[0, 0].set_xlabel('Sequence Length', fontsize=12)\n axes[0, 0].set_ylabel('Frequency', fontsize=12)\n axes[0, 0].set_title('Length Distribution Comparison', fontsize=14, fontweight='bold')\n axes[0, 0].legend()\n axes[0, 0].grid(alpha=0.3)\n \n # Box plots\n axes[0, 1].boxplot([train_lengths, test_lengths], \n labels=['Train', 'Test'],\n patch_artist=True,\n boxprops=dict(alpha=0.7))\n axes[0, 1].set_ylabel('Sequence Length', fontsize=12)\n axes[0, 1].set_title('Box Plot Comparison', fontsize=14, fontweight='bold')\n axes[0, 1].grid(alpha=0.3)\n \n # Violin plots\n data_to_plot = [train_lengths, test_lengths]\n parts = axes[1, 0].violinplot(data_to_plot, positions=[1, 2], \n showmeans=True, showmedians=True)\n axes[1, 0].set_xticks([1, 2])\n axes[1, 0].set_xticklabels(['Train', 'Test'])\n axes[1, 0].set_ylabel('Sequence Length', fontsize=12)\n axes[1, 0].set_title('Violin Plot Comparison', fontsize=14, fontweight='bold')\n axes[1, 0].grid(alpha=0.3)\n \n # Statistical comparison table\n axes[1, 1].axis('off')\n stats_data = [\n ['Metric', 'Train', 'Test', 'Difference'],\n ['Count', f'{len(train_lengths):,}', f'{len(test_lengths):,}', '-'],\n ['Mean', f'{np.mean(train_lengths):.1f}', f'{np.mean(test_lengths):.1f}', \n f'{np.mean(train_lengths) - np.mean(test_lengths):+.1f}'],\n ['Median', f'{np.median(train_lengths):.1f}', f'{np.median(test_lengths):.1f}',\n f'{np.median(train_lengths) - np.median(test_lengths):+.1f}'],\n ['Std Dev', f'{np.std(train_lengths):.1f}', f'{np.std(test_lengths):.1f}',\n f'{np.std(train_lengths) - np.std(test_lengths):+.1f}'],\n ['Min', f'{min(train_lengths)}', f'{min(test_lengths)}', '-'],\n ['Max', f'{max(train_lengths):,}', f'{max(test_lengths):,}', '-'],\n ]\n \n table = axes[1, 1].table(cellText=stats_data, loc='center', cellLoc='center')\n table.auto_set_font_size(False)\n table.set_fontsize(10)\n table.scale(1, 2)\n \n # Color header row\n for i in range(4):\n table[(0, i)].set_facecolor('#4ECDC4')\n table[(0, i)].set_text_props(weight='bold', color='white')\n \n axes[1, 1].set_title('Statistical Comparison', fontsize=14, fontweight='bold', pad=20)\n \n plt.tight_layout()\n \n if save_path:\n plt.savefig(save_path, dpi=300, bbox_inches='tight')\n \n plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T19:59:51.439822Z","iopub.execute_input":"2025-10-19T19:59:51.440085Z","iopub.status.idle":"2025-10-19T19:59:51.647983Z","shell.execute_reply.started":"2025-10-19T19:59:51.440061Z","shell.execute_reply":"2025-10-19T19:59:51.647195Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plot_sequence_comparison(train_sequences, test_sequences, save_path=\"/kaggle/working/sequence_comparison.png\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T19:59:51.648792Z","iopub.execute_input":"2025-10-19T19:59:51.649045Z","iopub.status.idle":"2025-10-19T19:59:59.148558Z","shell.execute_reply.started":"2025-10-19T19:59:51.649020Z","shell.execute_reply":"2025-10-19T19:59:59.147794Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Target Label Analysis Functions (For Annotations)","metadata":{}},{"cell_type":"code","source":"def plot_labels_per_protein(annotations: pd.DataFrame,\n max_labels: int = 150,\n save_path: str = None):\n \"\"\"\n Plot distribution of labels per protein.\n \n Args:\n annotations: DataFrame with protein_id, go_term, ontology columns\n max_labels: Maximum number of labels to display on x-axis\n save_path: Optional path to save figure\n \"\"\"\n labels_per_protein = annotations.groupby('protein_id').size()\n \n fig, axes = plt.subplots(2, 2, figsize=(16, 12))\n \n # Histogram\n axes[0, 0].hist(labels_per_protein, bins=50, color='coral', \n alpha=0.7, edgecolor='black')\n axes[0, 0].axvline(labels_per_protein.mean(), color='red', linestyle='--',\n linewidth=2, label=f'Mean: {labels_per_protein.mean():.1f}')\n axes[0, 0].axvline(labels_per_protein.median(), color='orange', linestyle='--',\n linewidth=2, label=f'Median: {labels_per_protein.median():.1f}')\n axes[0, 0].set_xlabel('Number of Labels', fontsize=12)\n axes[0, 0].set_ylabel('Number of Proteins', fontsize=12)\n axes[0, 0].set_title('Labels per Protein Distribution', fontsize=14, fontweight='bold')\n axes[0, 0].legend()\n axes[0, 0].grid(alpha=0.3)\n \n # Cumulative distribution\n sorted_labels = np.sort(labels_per_protein)\n cumulative = np.arange(1, len(sorted_labels) + 1) / len(sorted_labels) * 100\n axes[0, 1].plot(sorted_labels, cumulative, color='steelblue', linewidth=2)\n axes[0, 1].set_xlabel('Number of Labels', fontsize=12)\n axes[0, 1].set_ylabel('Cumulative Percentage (%)', fontsize=12)\n axes[0, 1].set_title('Cumulative Distribution', fontsize=14, fontweight='bold')\n axes[0, 1].grid(alpha=0.3)\n axes[0, 1].axhline(50, color='red', linestyle='--', alpha=0.5, label='50%')\n axes[0, 1].axhline(90, color='orange', linestyle='--', alpha=0.5, label='90%')\n axes[0, 1].legend()\n \n # Bar chart by ranges\n ranges = [(1, 5), (6, 10), (11, 20), (21, 30), (31, 50), \n (51, 75), (76, 100), (101, max_labels)]\n range_counts = []\n range_labels = []\n \n for start, end in ranges:\n count = ((labels_per_protein >= start) & (labels_per_protein <= end)).sum()\n range_counts.append(count)\n if end == max_labels:\n range_labels.append(f'{start}+')\n else:\n range_labels.append(f'{start}-{end}')\n \n colors_gradient = plt.cm.RdYlGn_r(np.linspace(0.2, 0.8, len(range_counts)))\n axes[1, 0].bar(range_labels, range_counts, color=colors_gradient, \n alpha=0.8, edgecolor='black')\n axes[1, 0].set_xlabel('Label Range', fontsize=12)\n axes[1, 0].set_ylabel('Number of Proteins', fontsize=12)\n axes[1, 0].set_title('Proteins by Label Range', fontsize=14, fontweight='bold')\n axes[1, 0].tick_params(axis='x', rotation=45)\n axes[1, 0].grid(axis='y', alpha=0.3)\n \n # Statistics table\n axes[1, 1].axis('off')\n stats_data = [\n ['Metric', 'Value'],\n ['Total Proteins', f'{len(labels_per_protein):,}'],\n ['Mean Labels', f'{labels_per_protein.mean():.2f}'],\n ['Median Labels', f'{labels_per_protein.median():.0f}'],\n ['Std Dev', f'{labels_per_protein.std():.2f}'],\n ['Min Labels', f'{labels_per_protein.min()}'],\n ['Max Labels', f'{labels_per_protein.max()}'],\n ['25th Percentile', f'{labels_per_protein.quantile(0.25):.0f}'],\n ['75th Percentile', f'{labels_per_protein.quantile(0.75):.0f}'],\n ['Proteins with >50 labels', f'{(labels_per_protein > 50).sum():,}'],\n ]\n \n table = axes[1, 1].table(cellText=stats_data, loc='center', cellLoc='left',\n colWidths=[0.6, 0.4])\n table.auto_set_font_size(False)\n table.set_fontsize(11)\n table.scale(1, 2.5)\n \n for i in range(2):\n table[(0, i)].set_facecolor('#FF6B6B')\n table[(0, i)].set_text_props(weight='bold', color='white')\n \n axes[1, 1].set_title('Statistics Summary', fontsize=14, fontweight='bold', pad=20)\n \n plt.tight_layout()\n \n if save_path:\n plt.savefig(save_path, dpi=300, bbox_inches='tight')\n \n plt.show()\n \n # Print detailed statistics\n print(f\"\\n{'='*60}\")\n print(f\"LABELS PER PROTEIN STATISTICS\")\n print(f\"{'='*60}\")\n print(f\"Total proteins: {len(labels_per_protein):,}\")\n print(f\"Mean labels: {labels_per_protein.mean():.2f}\")\n print(f\"Median labels: {labels_per_protein.median():.0f}\")\n print(f\"Std deviation: {labels_per_protein.std():.2f}\")\n print(f\"Min labels: {labels_per_protein.min()}\")\n print(f\"Max labels: {labels_per_protein.max()}\")\n print(f\"\\nPercentiles:\")\n for p in [10, 25, 50, 75, 90, 95, 99]:\n print(f\" {p}th: {labels_per_protein.quantile(p/100):.0f}\")\n print(f\"{'='*60}\\n\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T19:59:59.149422Z","iopub.execute_input":"2025-10-19T19:59:59.149649Z","iopub.status.idle":"2025-10-19T19:59:59.166925Z","shell.execute_reply.started":"2025-10-19T19:59:59.149632Z","shell.execute_reply":"2025-10-19T19:59:59.166370Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plot_labels_per_protein(annotations, save_path=\"/kaggle/working/labels_per_protein.png\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T19:59:59.167669Z","iopub.execute_input":"2025-10-19T19:59:59.167887Z","iopub.status.idle":"2025-10-19T20:00:02.456823Z","shell.execute_reply.started":"2025-10-19T19:59:59.167871Z","shell.execute_reply":"2025-10-19T20:00:02.456187Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def plot_ontology_distribution(annotations: pd.DataFrame,\n save_path: str = None):\n \"\"\"\n Plot distribution of annotations across ontologies.\n \n Args:\n annotations: DataFrame with protein_id, go_term, ontology columns\n save_path: Optional path to save figure\n \"\"\"\n # Create a clean copy of the dataframe, skipping the header row\n if annotations.iloc[0, 0] == 'EntryID' and annotations.iloc[0, 1] == 'term' and annotations.iloc[0, 2] == 'aspect':\n clean_annotations = annotations.iloc[1:].copy()\n clean_annotations.columns = ['protein_id', 'go_term', 'ontology']\n else:\n clean_annotations = annotations.copy()\n \n # Map aspect codes to ontology codes\n aspect_to_ontology = {\n 'P': 'BPO', # Biological Process\n 'F': 'MFO', # Molecular Function \n 'C': 'CCO' # Cellular Component\n }\n \n # Convert aspect codes to ontology codes\n clean_annotations['ontology'] = clean_annotations['ontology'].map(aspect_to_ontology)\n \n # Check if annotations DataFrame is empty\n if clean_annotations.empty:\n print(\"Warning: annotations DataFrame is empty. No data to plot.\")\n return\n \n # Check for required columns\n required_columns = ['protein_id', 'go_term', 'ontology']\n missing_columns = [col for col in required_columns if col not in clean_annotations.columns]\n if missing_columns:\n print(f\"Warning: Missing required columns: {missing_columns}\")\n return\n \n # Filter out any rows with missing values in required columns\n clean_annotations = clean_annotations.dropna(subset=required_columns)\n \n # Check if we have any data left after filtering\n if clean_annotations.empty:\n print(\"Warning: No valid data after filtering missing values.\")\n return\n \n fig, axes = plt.subplots(2, 2, figsize=(16, 12))\n \n # Count by ontology\n ontology_counts = clean_annotations['ontology'].value_counts()\n ontology_order = ['BPO', 'MFO', 'CCO']\n \n # Reindex with the desired order, but only include existing ontologies\n existing_ontologies = [ont for ont in ontology_order if ont in ontology_counts.index]\n ontology_counts = ontology_counts.reindex(ontology_order, fill_value=0)\n \n ontology_names = {\n 'BPO': 'Biological Process',\n 'MFO': 'Molecular Function',\n 'CCO': 'Cellular Component'\n }\n \n colors = {'BPO': '#EF4444', 'MFO': '#3B82F6', 'CCO': '#10B981'}\n \n # Bar chart - annotations by ontology\n axes[0, 0].bar(ontology_counts.index, ontology_counts.values,\n color=[colors[ont] for ont in ontology_counts.index],\n alpha=0.8, edgecolor='black', linewidth=2)\n axes[0, 0].set_xlabel('Ontology', fontsize=12)\n axes[0, 0].set_ylabel('Number of Annotations', fontsize=12)\n axes[0, 0].set_title('Annotations by Ontology', fontsize=14, fontweight='bold')\n axes[0, 0].grid(axis='y', alpha=0.3)\n \n # Add count labels on bars\n for i, (ont, count) in enumerate(ontology_counts.items()):\n axes[0, 0].text(i, count, f'{count:,}\\n({count/ontology_counts.sum()*100:.1f}%)',\n ha='center', va='bottom', fontsize=10, fontweight='bold')\n \n # Pie chart\n axes[0, 1].pie(ontology_counts.values, \n labels=[f\"{ont}\\n{ontology_names[ont]}\" for ont in ontology_counts.index],\n autopct='%1.1f%%', colors=[colors[ont] for ont in ontology_counts.index],\n startangle=90, textprops={'fontsize': 11, 'weight': 'bold'})\n axes[0, 1].set_title('Ontology Distribution', fontsize=14, fontweight='bold')\n \n # Unique terms per ontology\n unique_terms = clean_annotations.groupby('ontology')['go_term'].nunique()\n unique_terms = unique_terms.reindex(ontology_order, fill_value=0)\n \n axes[1, 0].bar(unique_terms.index, unique_terms.values,\n color=[colors[ont] for ont in unique_terms.index],\n alpha=0.8, edgecolor='black', linewidth=2)\n axes[1, 0].set_xlabel('Ontology', fontsize=12)\n axes[1, 0].set_ylabel('Number of Unique GO Terms', fontsize=12)\n axes[1, 0].set_title('Unique Terms by Ontology', fontsize=14, fontweight='bold')\n axes[1, 0].grid(axis='y', alpha=0.3)\n \n # Add count labels\n for i, (ont, count) in enumerate(unique_terms.items()):\n axes[1, 0].text(i, count, f'{count:,}', ha='center', va='bottom',\n fontsize=10, fontweight='bold')\n \n # Unique proteins per ontology\n unique_proteins = clean_annotations.groupby('ontology')['protein_id'].nunique()\n unique_proteins = unique_proteins.reindex(ontology_order, fill_value=0)\n \n axes[1, 1].bar(unique_proteins.index, unique_proteins.values,\n color=[colors[ont] for ont in unique_proteins.index],\n alpha=0.8, edgecolor='black', linewidth=2)\n axes[1, 1].set_xlabel('Ontology', fontsize=12)\n axes[1, 1].set_ylabel('Number of Unique Proteins', fontsize=12)\n axes[1, 1].set_title('Unique Proteins by Ontology', fontsize=14, fontweight='bold')\n axes[1, 1].grid(axis='y', alpha=0.3)\n \n # Add count labels\n for i, (ont, count) in enumerate(unique_proteins.items()):\n axes[1, 1].text(i, count, f'{count:,}', ha='center', va='bottom',\n fontsize=10, fontweight='bold')\n \n plt.tight_layout()\n \n if save_path:\n plt.savefig(save_path, dpi=300, bbox_inches='tight')\n \n plt.show()\n \n # Print statistics\n print(f\"\\n{'='*60}\")\n print(f\"ONTOLOGY DISTRIBUTION\")\n print(f\"{'='*60}\")\n total = ontology_counts.sum()\n for ont in ontology_order:\n print(f\"\\n{ont} - {ontology_names[ont]}:\")\n print(f\" Annotations: {ontology_counts[ont]:,} ({ontology_counts[ont]/total*100:.1f}%)\")\n print(f\" Unique terms: {unique_terms[ont]:,}\")\n print(f\" Unique proteins: {unique_proteins[ont]:,}\")\n print(f\"\\nTotal annotations: {total:,}\")\n print(f\"{'='*60}\\n\")\n \n return clean_annotations","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T20:00:02.457646Z","iopub.execute_input":"2025-10-19T20:00:02.457899Z","iopub.status.idle":"2025-10-19T20:00:02.475622Z","shell.execute_reply.started":"2025-10-19T20:00:02.457881Z","shell.execute_reply":"2025-10-19T20:00:02.474860Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plot_ontology_distribution(annotations, save_path=\"/kaggle/working/ontology_distribution.png\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T20:00:02.476258Z","iopub.execute_input":"2025-10-19T20:00:02.476420Z","iopub.status.idle":"2025-10-19T20:00:04.950334Z","shell.execute_reply.started":"2025-10-19T20:00:02.476407Z","shell.execute_reply":"2025-10-19T20:00:04.949499Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def plot_term_frequency(annotations: pd.DataFrame,\n top_n: int = 20,\n min_freq: int = 1,\n save_path: str = None):\n \"\"\"\n Plot GO term frequency distribution.\n \n Args:\n annotations: DataFrame with protein_id, go_term, ontology columns\n top_n: Number of top terms to display\n min_freq: Minimum frequency to include\n save_path: Optional path to save figure\n \"\"\"\n term_freq = annotations['go_term'].value_counts()\n \n fig, axes = plt.subplots(2, 2, figsize=(16, 12))\n \n # Top N terms (horizontal bar chart)\n top_terms = term_freq.head(top_n)\n y_pos = np.arange(len(top_terms))\n \n axes[0, 0].barh(y_pos, top_terms.values, color='steelblue', \n alpha=0.8, edgecolor='black')\n axes[0, 0].set_yticks(y_pos)\n axes[0, 0].set_yticklabels(top_terms.index, fontsize=9)\n axes[0, 0].set_xlabel('Number of Proteins', fontsize=12)\n axes[0, 0].set_title(f'Top {top_n} Most Frequent GO Terms', \n fontsize=14, fontweight='bold')\n axes[0, 0].invert_yaxis()\n axes[0, 0].grid(axis='x', alpha=0.3)\n \n # Add count labels\n for i, v in enumerate(top_terms.values):\n axes[0, 0].text(v, i, f' {v:,}', va='center', fontsize=9)\n \n # Frequency distribution (log scale)\n axes[0, 1].hist(np.log10(term_freq.values), bins=50, \n color='coral', alpha=0.7, edgecolor='black')\n axes[0, 1].set_xlabel('log10(Term Frequency)', fontsize=12)\n axes[0, 1].set_ylabel('Number of Terms', fontsize=12)\n axes[0, 1].set_title('Term Frequency Distribution (Log Scale)', \n fontsize=14, fontweight='bold')\n axes[0, 1].grid(alpha=0.3)\n \n # Cumulative frequency\n sorted_freq = np.sort(term_freq.values)[::-1]\n cumulative = np.cumsum(sorted_freq) / sorted_freq.sum() * 100\n \n axes[1, 0].plot(range(len(cumulative)), cumulative, \n color='green', linewidth=2)\n axes[1, 0].set_xlabel('Number of Terms (Ranked)', fontsize=12)\n axes[1, 0].set_ylabel('Cumulative Coverage (%)', fontsize=12)\n axes[1, 0].set_title('Cumulative Term Coverage', fontsize=14, fontweight='bold')\n axes[1, 0].axhline(80, color='red', linestyle='--', \n alpha=0.5, label='80% coverage')\n axes[1, 0].axhline(90, color='orange', linestyle='--', \n alpha=0.5, label='90% coverage')\n axes[1, 0].legend()\n axes[1, 0].grid(alpha=0.3)\n \n # Class imbalance statistics\n axes[1, 1].axis('off')\n \n imbalance_ratio = term_freq.iloc[0] / term_freq.iloc[-1]\n rare_terms = (term_freq < 10).sum()\n singleton_terms = (term_freq == 1).sum()\n coverage_80 = np.where(cumulative >= 80)[0][0] if np.any(cumulative >= 80) else len(cumulative)\n \n stats_data = [\n ['Metric', 'Value'],\n ['Total Unique Terms', f'{len(term_freq):,}'],\n ['Most Common Term', f'{term_freq.iloc[0]:,}'],\n ['Least Common Term', f'{term_freq.iloc[-1]}'],\n ['Imbalance Ratio', f'{imbalance_ratio:.1f}x'],\n ['Terms with <10 proteins', f'{rare_terms:,} ({rare_terms/len(term_freq)*100:.1f}%)'],\n ['Singleton Terms', f'{singleton_terms:,} ({singleton_terms/len(term_freq)*100:.1f}%)'],\n ['Terms for 80% coverage', f'{coverage_80:,}'],\n ['Mean Frequency', f'{term_freq.mean():.1f}'],\n ['Median Frequency', f'{term_freq.median():.0f}'],\n ]\n \n table = axes[1, 1].table(cellText=stats_data, loc='center', cellLoc='left',\n colWidths=[0.65, 0.35])\n table.auto_set_font_size(False)\n table.set_fontsize(11)\n table.scale(1, 2.5)\n \n for i in range(2):\n table[(0, i)].set_facecolor('#FF6B6B')\n table[(0, i)].set_text_props(weight='bold', color='white')\n \n axes[1, 1].set_title('Class Imbalance Statistics', \n fontsize=14, fontweight='bold', pad=20)\n \n plt.tight_layout()\n \n if save_path:\n plt.savefig(save_path, dpi=300, bbox_inches='tight')\n \n plt.show()\n \n # Print detailed statistics\n print(f\"\\n{'='*60}\")\n print(f\"GO TERM FREQUENCY ANALYSIS\")\n print(f\"{'='*60}\")\n print(f\"Total unique terms: {len(term_freq):,}\")\n print(f\"Total annotations: {term_freq.sum():,}\")\n print(f\"\\nFrequency Statistics:\")\n print(f\" Most common: {term_freq.iloc[0]:,} ({term_freq.index[0]})\")\n print(f\" Least common: {term_freq.iloc[-1]} ({term_freq.index[-1]})\")\n print(f\" Imbalance ratio: {imbalance_ratio:.1f}x\")\n print(f\" Mean: {term_freq.mean():.1f}\")\n print(f\" Median: {term_freq.median():.0f}\")\n print(f\"\\nClass Imbalance:\")\n print(f\" Terms with <10 proteins: {rare_terms:,} ({rare_terms/len(term_freq)*100:.1f}%)\")\n print(f\" Singleton terms: {singleton_terms:,} ({singleton_terms/len(term_freq)*100:.1f}%)\")\n print(f\" Terms for 80% coverage: {coverage_80:,}\")\n print(f\"{'='*60}\\n\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T20:00:04.951170Z","iopub.execute_input":"2025-10-19T20:00:04.951391Z","iopub.status.idle":"2025-10-19T20:00:04.967325Z","shell.execute_reply.started":"2025-10-19T20:00:04.951365Z","shell.execute_reply":"2025-10-19T20:00:04.966472Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plot_term_frequency(annotations, save_path=\"/kaggle/working/term_frequency.png\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T20:00:04.968083Z","iopub.execute_input":"2025-10-19T20:00:04.968411Z","iopub.status.idle":"2025-10-19T20:00:08.250967Z","shell.execute_reply.started":"2025-10-19T20:00:04.968386Z","shell.execute_reply":"2025-10-19T20:00:08.250168Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport networkx as nx\nfrom collections import Counter\nfrom typing import Dict, Tuple, Optional\nimport warnings\nwarnings.filterwarnings('ignore')\n\n\ndef analyze_label_cooccurrence(annotations: pd.DataFrame, \n min_support: int = None,\n auto_threshold: bool = True) -> pd.DataFrame:\n \"\"\"\n Analyze which GO terms frequently co-occur with adaptive thresholding.\n \n Args:\n annotations: DataFrame with protein_id, go_term, ontology columns\n min_support: Minimum number of proteins for pair to be included\n auto_threshold: If True, automatically determine min_support\n \n Returns:\n DataFrame with term pairs and co-occurrence counts\n \"\"\"\n if len(annotations) == 0:\n print(\"Warning: Empty annotations DataFrame\")\n return pd.DataFrame(columns=['term1', 'term2', 'count', 'ont1', 'ont2'])\n \n # Group terms by protein\n protein_terms = annotations.groupby('protein_id')['go_term'].apply(list).to_dict()\n \n # Create ontology mapping for later use\n term_to_ont = annotations.drop_duplicates('go_term').set_index('go_term')['ontology'].to_dict()\n \n # Count pairs\n pair_counts = Counter()\n for terms in protein_terms.values():\n if len(terms) > 1:\n for i, term1 in enumerate(terms):\n for term2 in terms[i+1:]:\n pair = tuple(sorted([term1, term2]))\n pair_counts[pair] += 1\n \n if len(pair_counts) == 0:\n print(\"Warning: No co-occurring term pairs found\")\n return pd.DataFrame(columns=['term1', 'term2', 'count', 'ont1', 'ont2'])\n \n # Auto-determine min_support if requested\n if auto_threshold and min_support is None:\n counts = list(pair_counts.values())\n min_support = max(1, int(np.percentile(counts, 50))) # Use median as threshold\n print(f\"Auto-determined min_support: {min_support} (median of co-occurrence counts)\")\n \n elif min_support is None:\n min_support = 10\n \n # Filter and create DataFrame\n cooccurrence = pd.DataFrame([\n {\n 'term1': pair[0], \n 'term2': pair[1], \n 'count': count,\n 'ont1': term_to_ont.get(pair[0], 'Unknown'),\n 'ont2': term_to_ont.get(pair[1], 'Unknown')\n }\n for pair, count in pair_counts.items()\n if count >= min_support\n ])\n \n if len(cooccurrence) == 0:\n # Try with lower threshold\n new_min = max(1, min_support // 2)\n print(f\"No pairs found with min_support={min_support}. Retrying with {new_min}...\")\n cooccurrence = pd.DataFrame([\n {\n 'term1': pair[0], \n 'term2': pair[1], \n 'count': count,\n 'ont1': term_to_ont.get(pair[0], 'Unknown'),\n 'ont2': term_to_ont.get(pair[1], 'Unknown')\n }\n for pair, count in pair_counts.items()\n if count >= new_min\n ])\n min_support = new_min\n \n if len(cooccurrence) > 0:\n cooccurrence = cooccurrence.sort_values('count', ascending=False)\n \n # Add pair type column\n cooccurrence['pair_type'] = cooccurrence.apply(\n lambda row: 'within_ontology' if row['ont1'] == row['ont2'] else 'cross_ontology',\n axis=1\n )\n \n print(f\"\\n{'='*60}\")\n print(f\"LABEL CO-OCCURRENCE ANALYSIS\")\n print(f\"{'='*60}\")\n print(f\"Total unique term pairs: {len(pair_counts):,}\")\n print(f\"Pairs with ā‰„{min_support} co-occurrences: {len(cooccurrence):,}\")\n print(f\" Within-ontology pairs: {(cooccurrence['pair_type'] == 'within_ontology').sum():,}\")\n print(f\" Cross-ontology pairs: {(cooccurrence['pair_type'] == 'cross_ontology').sum():,}\")\n \n # Statistics by ontology pair type\n print(f\"\\nPairs by Ontology Combination:\")\n for ont1 in ['BPO', 'MFO', 'CCO']:\n for ont2 in ['BPO', 'MFO', 'CCO']:\n if ont1 <= ont2: # Avoid duplicates\n mask = ((cooccurrence['ont1'] == ont1) & (cooccurrence['ont2'] == ont2)) | \\\n ((cooccurrence['ont1'] == ont2) & (cooccurrence['ont2'] == ont1))\n count = mask.sum()\n if count > 0:\n label = f\"{ont1}-{ont2}\" if ont1 != ont2 else ont1\n print(f\" {label}: {count:,} pairs\")\n \n print(f\"\\nTop 10 Co-occurring Term Pairs:\")\n for idx, row in cooccurrence.head(10).iterrows():\n print(f\" {row['term1']} ({row['ont1']}) ā†” {row['term2']} ({row['ont2']}): {row['count']:,}\")\n print(f\"{'='*60}\\n\")\n else:\n print(f\"\\nNo term pairs found with ā‰„{min_support} co-occurrences\")\n \n return cooccurrence\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T20:00:08.251883Z","iopub.execute_input":"2025-10-19T20:00:08.252102Z","iopub.status.idle":"2025-10-19T20:00:08.266435Z","shell.execute_reply.started":"2025-10-19T20:00:08.252086Z","shell.execute_reply":"2025-10-19T20:00:08.265680Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import seaborn as sns\nimport networkx as nx\nfrom matplotlib.patches import FancyBboxPatch\nimport numpy as np\n\ndef visualize_cooccurrence(cooccurrence_df: pd.DataFrame, \n annotations: pd.DataFrame,\n top_n: int = 20,\n save_path: str = None):\n \"\"\"\n Visualize GO term co-occurrence analysis results.\n \n Args:\n cooccurrence_df: DataFrame from analyze_label_cooccurrence\n annotations: Original annotations DataFrame for context\n top_n: Number of top pairs to visualize\n save_path: Optional path to save figure\n \"\"\"\n if cooccurrence_df.empty:\n print(\"No co-occurrence data to visualize.\")\n return\n \n # Create figure with multiple subplots\n fig = plt.figure(figsize=(20, 15))\n \n # Define grid layout\n gs = fig.add_gridspec(3, 2, height_ratios=[1, 1, 0.8])\n \n # 1. Heatmap of top co-occurring pairs\n ax1 = fig.add_subplot(gs[0, 0])\n plot_cooccurrence_heatmap(cooccurrence_df, top_n, ax1)\n \n # 2. Network graph\n ax2 = fig.add_subplot(gs[0, 1])\n plot_cooccurrence_network(cooccurrence_df, top_n, ax2)\n \n # 3. Distribution of co-occurrence counts\n ax3 = fig.add_subplot(gs[1, 0])\n plot_cooccurrence_distribution(cooccurrence_df, ax3)\n \n # 4. Top terms by overall frequency\n ax4 = fig.add_subplot(gs[1, 1])\n plot_term_frequency(annotations, cooccurrence_df, top_n, ax4)\n \n # 5. Statistics table\n ax5 = fig.add_subplot(gs[2, :])\n plot_cooccurrence_stats(cooccurrence_df, annotations, ax5)\n \n plt.tight_layout()\n \n if save_path:\n plt.savefig(save_path, dpi=300, bbox_inches='tight')\n \n plt.show()\n\ndef plot_cooccurrence_heatmap(cooccurrence_df: pd.DataFrame, top_n: int, ax):\n \"\"\"Plot heatmap of top co-occurring term pairs.\"\"\"\n top_pairs = cooccurrence_df.head(top_n).copy()\n \n # Get unique terms from top pairs\n all_terms = set(top_pairs['term1']).union(set(top_pairs['term2']))\n term_list = sorted(list(all_terms))\n \n # Create co-occurrence matrix\n cooccurrence_matrix = pd.DataFrame(0, index=term_list, columns=term_list)\n \n for _, row in top_pairs.iterrows():\n cooccurrence_matrix.loc[row['term1'], row['term2']] = row['count']\n cooccurrence_matrix.loc[row['term2'], row['term1']] = row['count']\n \n # Plot heatmap\n mask = np.triu(np.ones_like(cooccurrence_matrix, dtype=bool), k=1)\n sns.heatmap(cooccurrence_matrix, mask=mask, cmap='YlOrRd', \n annot=True, fmt='d', ax=ax, cbar_kws={'label': 'Co-occurrence Count'})\n \n ax.set_title(f'Top {top_n} GO Term Co-occurrence Heatmap', fontsize=14, fontweight='bold')\n ax.tick_params(axis='x', rotation=45)\n ax.tick_params(axis='y', rotation=0)\n\ndef plot_cooccurrence_network(cooccurrence_df: pd.DataFrame, top_n: int, ax):\n \"\"\"Plot network graph of co-occurring terms.\"\"\"\n top_pairs = cooccurrence_df.head(top_n).copy()\n \n # Create graph\n G = nx.Graph()\n \n # Add edges with weights\n for _, row in top_pairs.iterrows():\n G.add_edge(row['term1'], row['term2'], weight=row['count'])\n \n if len(G.edges()) == 0:\n ax.text(0.5, 0.5, 'No co-occurrence pairs to display', \n ha='center', va='center', transform=ax.transAxes)\n ax.set_title('Co-occurrence Network', fontsize=14, fontweight='bold')\n return\n \n # Calculate node positions\n pos = nx.spring_layout(G, k=1, iterations=50)\n \n # Calculate node sizes based on degree\n degrees = dict(G.degree())\n node_sizes = [300 + degrees[node] * 100 for node in G.nodes()]\n \n # Calculate edge widths based on co-occurrence count\n edge_weights = [G[u][v]['weight'] / 10 for u, v in G.edges()]\n \n # Draw network\n nx.draw_networkx_nodes(G, pos, node_size=node_sizes, \n node_color='lightblue', alpha=0.9, ax=ax)\n nx.draw_networkx_edges(G, pos, width=edge_weights, \n alpha=0.6, edge_color='gray', ax=ax)\n nx.draw_networkx_labels(G, pos, font_size=8, ax=ax)\n \n # Add edge labels for top 5 strongest connections\n edge_labels = {(u, v): f\"{G[u][v]['weight']}\" \n for u, v in list(G.edges())[:5]}\n nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels, \n font_size=7, ax=ax)\n \n ax.set_title(f'Co-occurrence Network (Top {top_n} Pairs)', \n fontsize=14, fontweight='bold')\n ax.axis('off')\n\ndef plot_cooccurrence_distribution(cooccurrence_df: pd.DataFrame, ax):\n \"\"\"Plot distribution of co-occurrence counts.\"\"\"\n counts = cooccurrence_df['count']\n \n ax.hist(counts, bins=30, alpha=0.7, color='skyblue', \n edgecolor='black', linewidth=0.5)\n ax.axvline(counts.mean(), color='red', linestyle='--', \n linewidth=2, label=f'Mean: {counts.mean():.1f}')\n ax.axvline(counts.median(), color='orange', linestyle='--', \n linewidth=2, label=f'Median: {counts.median():.1f}')\n \n ax.set_xlabel('Co-occurrence Count', fontsize=12)\n ax.set_ylabel('Frequency', fontsize=12)\n ax.set_title('Distribution of Co-occurrence Counts', \n fontsize=14, fontweight='bold')\n ax.legend()\n ax.grid(alpha=0.3)\n\ndef plot_term_frequency(annotations: pd.DataFrame, cooccurrence_df: pd.DataFrame, top_n: int, ax):\n \"\"\"Plot most frequent terms in co-occurring pairs.\"\"\"\n # Get all unique terms from co-occurrence pairs\n all_cooccur_terms = set(cooccurrence_df['term1']).union(set(cooccurrence_df['term2']))\n \n # Calculate overall frequency of these terms\n term_freq = annotations[annotations['go_term'].isin(all_cooccur_terms)]['go_term'].value_counts()\n top_terms = term_freq.head(top_n)\n \n # Get colors based on ontology\n ontology_map = annotations.set_index('go_term')['ontology'].to_dict()\n colors = []\n for term in top_terms.index:\n ont = ontology_map.get(term, 'Unknown')\n if ont == 'BPO': colors.append('#EF4444')\n elif ont == 'MFO': colors.append('#3B82F6')\n elif ont == 'CCO': colors.append('#10B981')\n else: colors.append('gray')\n \n bars = ax.barh(range(len(top_terms)), top_terms.values, color=colors, alpha=0.8)\n ax.set_yticks(range(len(top_terms)))\n ax.set_yticklabels([f\"{term}\" for term in top_terms.index], fontsize=9)\n ax.invert_yaxis()\n ax.set_xlabel('Frequency in Dataset', fontsize=12)\n ax.set_title(f'Top {top_n} Terms in Co-occurrence Pairs', \n fontsize=14, fontweight='bold')\n ax.grid(axis='x', alpha=0.3)\n \n # Add value labels\n for i, (bar, count) in enumerate(zip(bars, top_terms.values)):\n ax.text(bar.get_width() + bar.get_width() * 0.01, bar.get_y() + bar.get_height()/2,\n f'{count:,}', ha='left', va='center', fontsize=8)\n\ndef plot_cooccurrence_stats(cooccurrence_df: pd.DataFrame, annotations: pd.DataFrame, ax):\n \"\"\"Plot statistics table.\"\"\"\n ax.axis('off')\n \n # Calculate statistics\n total_pairs = len(cooccurrence_df)\n total_cooccurrences = cooccurrence_df['count'].sum()\n avg_cooccurrence = cooccurrence_df['count'].mean()\n max_cooccurrence = cooccurrence_df['count'].max()\n \n # Get unique terms involved in co-occurrence\n cooccur_terms = set(cooccurrence_df['term1']).union(set(cooccurrence_df['term2']))\n \n # Calculate statistics by ontology\n ontologies = ['BPO', 'MFO', 'CCO']\n stats_by_ontology = []\n \n for ont in ontologies:\n ont_terms = set(annotations[annotations['ontology'] == ont]['go_term'])\n ont_cooccur_terms = cooccur_terms.intersection(ont_terms)\n ont_pairs = cooccurrence_df[\n (cooccurrence_df['term1'].isin(ont_terms)) | \n (cooccurrence_df['term2'].isin(ont_terms))\n ]\n stats_by_ontology.append({\n 'ontology': ont,\n 'terms_in_cooccurrence': len(ont_cooccur_terms),\n 'pairs_involving_ontology': len(ont_pairs)\n })\n \n # Create table data\n table_data = [\n ['Total Pairs', f'{total_pairs:,}'],\n ['Total Co-occurrences', f'{total_cooccurrences:,}'],\n ['Average Co-occurrence', f'{avg_cooccurrence:.1f}'],\n ['Maximum Co-occurrence', f'{max_cooccurrence:,}'],\n ['Unique Terms in Pairs', f'{len(cooccur_terms):,}'],\n ['', ''],\n ['Ontology', 'Terms in Pairs / Pairs Involving'],\n ]\n \n for stats in stats_by_ontology:\n table_data.append([\n stats['ontology'],\n f\"{stats['terms_in_cooccurrence']:,} / {stats['pairs_involving_ontology']:,}\"\n ])\n \n # Create table\n table = ax.table(cellText=table_data, \n cellLoc='left',\n loc='center',\n bbox=[0.1, 0.1, 0.8, 0.8])\n \n table.auto_set_font_size(False)\n table.set_fontsize(10)\n table.scale(1, 2)\n \n # Style header row\n for i in range(len(table_data)):\n for j in range(len(table_data[0])):\n if i == 0 or i == 6: # Header rows\n table[(i, j)].set_facecolor('#4B5563')\n table[(i, j)].set_text_props(weight='bold', color='white')\n elif i > 6: # Ontology rows\n if j == 0:\n table[(i, j)].set_facecolor('#E5E7EB')\n \n ax.set_title('Co-occurrence Analysis Statistics', \n fontsize=16, fontweight='bold', pad=20)\n\n# Updated wrapper function that combines analysis and visualization\ndef analyze_and_visualize_cooccurrence(annotations: pd.DataFrame, \n min_support: int = 10,\n auto_threshold: bool = True,\n top_n: int = 20,\n save_path: str = None) -> pd.DataFrame:\n \"\"\"\n Analyze GO term co-occurrence and create visualizations.\n \n Args:\n annotations: DataFrame with annotations\n min_support: Minimum number of proteins for pair to be included\n top_n: Number of top pairs to visualize\n save_path: Optional path to save figure\n \n Returns:\n DataFrame with term pairs and co-occurrence counts\n \"\"\"\n # Perform analysis\n cooccurrence_df = analyze_label_cooccurrence(annotations, min_support, auto_threshold)\n \n # Create visualizations if we have data\n if not cooccurrence_df.empty:\n visualize_cooccurrence(cooccurrence_df, annotations, top_n, save_path)\n \n return cooccurrence_df","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T20:00:08.270569Z","iopub.execute_input":"2025-10-19T20:00:08.270969Z","iopub.status.idle":"2025-10-19T20:00:08.296775Z","shell.execute_reply.started":"2025-10-19T20:00:08.270954Z","shell.execute_reply":"2025-10-19T20:00:08.296115Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Pre-process your annotations DataFrame first\ndef prepare_annotations_for_analysis(annotations):\n \"\"\"Prepare annotations DataFrame for co-occurrence analysis.\"\"\"\n clean_annotations = annotations.copy()\n \n # Skip header row if present\n if clean_annotations.iloc[0, 0] == 'EntryID':\n clean_annotations = clean_annotations.iloc[1:].copy()\n clean_annotations.columns = ['protein_id', 'go_term', 'ontology']\n \n # Map aspect codes to ontology codes\n aspect_to_ontology = {\n 'P': 'BPO', # Biological Process\n 'F': 'MFO', # Molecular Function \n 'C': 'CCO' # Cellular Component\n }\n \n clean_annotations['ontology'] = clean_annotations['ontology'].map(aspect_to_ontology)\n return clean_annotations\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T20:00:08.297621Z","iopub.execute_input":"2025-10-19T20:00:08.297838Z","iopub.status.idle":"2025-10-19T20:00:08.315740Z","shell.execute_reply.started":"2025-10-19T20:00:08.297821Z","shell.execute_reply":"2025-10-19T20:00:08.314889Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Usage:\nannotations = prepare_annotations_for_analysis(annotations)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T20:00:08.316650Z","iopub.execute_input":"2025-10-19T20:00:08.316999Z","iopub.status.idle":"2025-10-19T20:00:08.382550Z","shell.execute_reply.started":"2025-10-19T20:00:08.316970Z","shell.execute_reply":"2025-10-19T20:00:08.381930Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"analyze_and_visualize_cooccurrence(\n annotations, \n auto_threshold=True,\n top_n=20,\n save_path='/kaggle/working/co-occurence.png')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T20:00:08.383259Z","iopub.execute_input":"2025-10-19T20:00:08.383436Z","iopub.status.idle":"2025-10-19T20:00:19.121017Z","shell.execute_reply.started":"2025-10-19T20:00:08.383422Z","shell.execute_reply":"2025-10-19T20:00:19.120276Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Graph Analysis Functions (For Ontologies)","metadata":{}},{"cell_type":"code","source":"import networkx as nx\nimport numpy as np\nfrom typing import Dict\nimport matplotlib.pyplot as plt\n\ndef analyze_go_hierarchy(graph: nx.DiGraph, annotations: pd.DataFrame) -> Dict:\n \"\"\"\n Analyze GO ontology graph structure.\n \n Args:\n graph: NetworkX graph instance of GO ontology\n annotations: DataFrame with annotations (to find annotated terms)\n \n Returns:\n Dictionary with graph statistics\n \"\"\"\n stats = {}\n \n # Basic graph statistics\n stats['num_nodes'] = graph.number_of_nodes()\n stats['num_edges'] = graph.number_of_edges()\n \n # Find annotated terms\n annotated_terms = set(annotations['go_term'].unique())\n annotated_in_graph = annotated_terms & set(graph.nodes())\n \n stats['annotated_terms_in_graph'] = len(annotated_in_graph)\n stats['annotated_terms_not_in_graph'] = len(annotated_terms - annotated_in_graph)\n \n # Calculate term depths (distance from root) - FIXED VERSION\n root_terms = {\n 'GO:0008150', # biological_process\n 'GO:0003674', # molecular_function\n 'GO:0005575', # cellular_component\n }\n \n # Find actual roots (nodes with no incoming edges)\n actual_roots = [node for node in graph.nodes() if graph.in_degree(node) == 0]\n stats['actual_roots'] = actual_roots\n \n depths = {}\n valid_roots_found = 0\n \n for root in root_terms:\n if root in graph:\n try:\n # Use reverse graph for depth calculation (distance TO root)\n reverse_graph = graph.reverse()\n depth_dict = nx.single_source_shortest_path_length(reverse_graph, root)\n if depth_dict: # Only include if we found reachable nodes\n depths[root] = depth_dict\n valid_roots_found += 1\n print(f\"Found root {root} with {len(depth_dict)} reachable terms\")\n except nx.NetworkXError as e:\n print(f\"Warning: Root term {root} error: {e}\")\n continue\n \n # Alternative: Calculate depth from any root\n if not depths:\n print(\"Trying alternative depth calculation from any root...\")\n # Use nodes with in_degree 0 as roots\n for root in actual_roots[:3]: # Use first 3 actual roots\n try:\n reverse_graph = graph.reverse()\n depth_dict = nx.single_source_shortest_path_length(reverse_graph, root)\n if len(depth_dict) > 10: # Only include meaningful roots\n depths[root] = depth_dict\n print(f\"Using root {root} with {len(depth_dict)} reachable terms\")\n break\n except:\n continue\n \n if depths:\n all_depths = []\n for depth_dict in depths.values():\n all_depths.extend(depth_dict.values())\n \n if all_depths:\n stats['term_depth'] = {\n 'mean': np.mean(all_depths),\n 'max': max(all_depths),\n 'min': min(all_depths),\n 'std': np.std(all_depths),\n 'num_measured_terms': len(all_depths)\n }\n print(f\"Depth calculation: measured {len(all_depths)} terms\")\n else:\n stats['term_depth'] = {'error': 'No depth values calculated'}\n else:\n stats['term_depth'] = {'error': 'No valid roots found for depth calculation'}\n \n # Additional graph analysis\n if graph.number_of_nodes() > 0:\n # Check if graph is a DAG (should be for ontology)\n stats['is_dag'] = nx.is_directed_acyclic_graph(graph)\n \n # Calculate connectivity on a sample for large graphs\n if graph.number_of_nodes() > 10000:\n print(\"Large graph detected, using sampling for connectivity analysis...\")\n # Sample nodes for connectivity analysis\n sample_nodes = list(graph.nodes())[:1000]\n subgraph = graph.subgraph(sample_nodes)\n stats['is_weakly_connected'] = nx.is_weakly_connected(subgraph)\n stats['is_strongly_connected'] = nx.is_strongly_connected(subgraph)\n stats['connectivity_note'] = 'Based on 1000 node sample'\n else:\n stats['is_weakly_connected'] = nx.is_weakly_connected(graph)\n stats['is_strongly_connected'] = nx.is_strongly_connected(graph)\n \n # Node degree analysis (optimized for large graphs)\n if graph.number_of_nodes() > 50000:\n print(\"Large graph detected, sampling for degree analysis...\")\n sample_size = min(10000, graph.number_of_nodes())\n sample_nodes = np.random.choice(list(graph.nodes()), sample_size, replace=False)\n in_degrees = [graph.in_degree(node) for node in sample_nodes]\n out_degrees = [graph.out_degree(node) for node in sample_nodes]\n stats['degree_sample_size'] = sample_size\n else:\n in_degrees = [d for n, d in graph.in_degree()]\n out_degrees = [d for n, d in graph.out_degree()]\n \n stats['degree_analysis'] = {\n 'in_degree_mean': np.mean(in_degrees) if in_degrees else 0,\n 'in_degree_max': max(in_degrees) if in_degrees else 0,\n 'out_degree_mean': np.mean(out_degrees) if out_degrees else 0,\n 'out_degree_max': max(out_degrees) if out_degrees else 0,\n 'in_degree_median': np.median(in_degrees) if in_degrees else 0,\n 'out_degree_median': np.median(out_degrees) if out_degrees else 0,\n }\n \n # Component analysis\n if graph.number_of_nodes() <= 50000:\n weakly_components = list(nx.weakly_connected_components(graph))\n stats['weakly_connected_components'] = len(weakly_components)\n if weakly_components:\n stats['largest_component_size'] = len(max(weakly_components, key=len))\n else:\n stats['component_note'] = 'Skipped for large graph'\n \n print(\"\\nGO Ontology Graph Analysis:\")\n print(f\"Total terms in ontology: {stats['num_nodes']:,}\")\n print(f\"Total relationships: {stats['num_edges']:,}\")\n print(f\"Annotated terms in training: {stats['annotated_terms_in_graph']:,}\")\n print(f\"Annotated terms not in graph: {stats['annotated_terms_not_in_graph']:,}\")\n \n if 'term_depth' in stats and 'error' not in stats['term_depth']:\n depth_stats = stats['term_depth']\n print(f\"\\nTerm depth statistics:\")\n print(f\" Mean: {depth_stats['mean']:.2f}\")\n print(f\" Max: {depth_stats['max']}\")\n print(f\" Min: {depth_stats['min']}\")\n print(f\" Std: {depth_stats['std']:.2f}\")\n print(f\" Terms measured: {depth_stats['num_measured_terms']:,}\")\n elif 'term_depth' in stats:\n print(f\"\\nDepth analysis: {stats['term_depth']['error']}\")\n \n if 'degree_analysis' in stats:\n deg = stats['degree_analysis']\n print(f\"\\nDegree analysis:\")\n print(f\" Mean in-degree: {deg['in_degree_mean']:.2f}\")\n print(f\" Median in-degree: {deg['in_degree_median']:.1f}\")\n print(f\" Max in-degree: {deg['in_degree_max']}\")\n print(f\" Mean out-degree: {deg['out_degree_mean']:.2f}\")\n print(f\" Median out-degree: {deg['out_degree_median']:.1f}\")\n print(f\" Max out-degree: {deg['out_degree_max']}\")\n if 'degree_sample_size' in stats:\n print(f\" (Based on {stats['degree_sample_size']:,} node sample)\")\n \n print(f\"\\nGraph properties:\")\n print(f\" Is DAG: {stats.get('is_dag', 'N/A')}\")\n print(f\" Weakly connected: {stats.get('is_weakly_connected', 'N/A')}\")\n print(f\" Strongly connected: {stats.get('is_strongly_connected', 'N/A')}\")\n if 'weakly_connected_components' in stats:\n print(f\" Weakly connected components: {stats['weakly_connected_components']}\")\n print(f\" Largest component size: {stats.get('largest_component_size', 'N/A'):,}\")\n \n return stats\n\n# OPTIMIZED Visualization function for GO hierarchy\ndef visualize_go_hierarchy(graph: nx.DiGraph, annotations: pd.DataFrame, \n save_path: str = None, max_plot_terms: int = 1000):\n \"\"\"\n Visualize GO ontology graph structure - OPTIMIZED VERSION.\n \n Args:\n graph: NetworkX graph instance of GO ontology\n annotations: DataFrame with annotations\n save_path: Optional path to save figure\n max_plot_terms: Maximum number of terms to include in plots\n \"\"\"\n print(\"Generating optimized visualizations...\")\n \n fig, axes = plt.subplots(2, 2, figsize=(16, 12))\n \n # 1. Optimized depth distribution\n ax1 = axes[0, 0]\n plot_depth_distribution_optimized(graph, ax1)\n \n # 2. Optimized degree distribution\n ax2 = axes[0, 1]\n plot_degree_distribution_optimized(graph, ax2)\n \n # 3. Statistics summary\n ax3 = axes[1, 0]\n plot_hierarchy_stats_optimized(graph, annotations, ax3)\n \n # 4. Annotated terms overview (fast alternative to network plot)\n ax4 = axes[1, 1]\n plot_annotated_overview(graph, annotations, ax4)\n \n plt.tight_layout()\n \n if save_path:\n plt.savefig(save_path, dpi=300, bbox_inches='tight')\n print(f\"Saved visualization to {save_path}\")\n \n plt.show()\n\ndef plot_depth_distribution_optimized(graph: nx.DiGraph, ax):\n \"\"\"Optimized depth distribution plot.\"\"\"\n print(\"Calculating depth distribution...\")\n \n # Use reverse graph for depth calculation\n reverse_graph = graph.reverse()\n \n # Try standard roots first\n root_terms = ['GO:0008150', 'GO:0003674', 'GO:0005575']\n depths = []\n \n for root in root_terms:\n if root in reverse_graph:\n try:\n depth_dict = nx.single_source_shortest_path_length(reverse_graph, root)\n if depth_dict:\n depths.extend(list(depth_dict.values()))\n print(f\" Root {root}: {len(depth_dict)} terms\")\n break # Use first successful root\n except:\n continue\n \n # Fallback: use any node with in_degree 0 as root\n if not depths:\n roots = [node for node in graph.nodes() if graph.in_degree(node) == 0]\n if roots:\n root = roots[0]\n try:\n depth_dict = nx.single_source_shortest_path_length(reverse_graph, root)\n depths.extend(list(depth_dict.values()))\n print(f\" Using root {root}: {len(depth_dict)} terms\")\n except:\n pass\n \n if depths:\n # Sample if too many points for performance\n if len(depths) > 10000:\n depths = np.random.choice(depths, 10000, replace=False)\n sampling_note = f\" (sampled from {len(depths):,})\"\n else:\n sampling_note = \"\"\n \n ax.hist(depths, bins=50, alpha=0.7, color='skyblue', \n edgecolor='black', linewidth=0.5)\n ax.axvline(np.mean(depths), color='red', linestyle='--', \n label=f'Mean: {np.mean(depths):.2f}')\n ax.set_xlabel('Depth from Root', fontsize=12)\n ax.set_ylabel('Frequency', fontsize=12)\n ax.set_title(f'Term Depth Distribution{sampling_note}', fontsize=14, fontweight='bold')\n ax.legend()\n ax.grid(alpha=0.3)\n else:\n ax.text(0.5, 0.5, 'No depth data available', \n ha='center', va='center', transform=ax.transAxes)\n ax.set_title('Term Depth Distribution', fontsize=14, fontweight='bold')\n\ndef plot_degree_distribution_optimized(graph: nx.DiGraph, ax):\n \"\"\"Optimized degree distribution plot.\"\"\"\n print(\"Calculating degree distribution...\")\n \n # Sample nodes for large graphs\n if graph.number_of_nodes() > 10000:\n sample_size = min(5000, graph.number_of_nodes())\n sample_nodes = np.random.choice(list(graph.nodes()), sample_size, replace=False)\n in_degrees = [graph.in_degree(node) for node in sample_nodes]\n out_degrees = [graph.out_degree(node) for node in sample_nodes]\n sampling_note = f\" (sampled {sample_size:,} nodes)\"\n else:\n in_degrees = [d for n, d in graph.in_degree()]\n out_degrees = [d for n, d in graph.out_degree()]\n sampling_note = \"\"\n \n if in_degrees and out_degrees:\n # Use logarithmic bins for better visualization\n max_degree = max(max(in_degrees), max(out_degrees))\n bins = np.logspace(0, np.log10(max_degree + 1), 30)\n \n ax.hist(in_degrees, bins=bins, alpha=0.7, label='In-degree', \n color='lightcoral', edgecolor='black', linewidth=0.5)\n ax.hist(out_degrees, bins=bins, alpha=0.7, label='Out-degree',\n color='lightgreen', edgecolor='black', linewidth=0.5)\n ax.set_xlabel('Degree', fontsize=12)\n ax.set_ylabel('Frequency', fontsize=12)\n ax.set_title(f'Degree Distribution{sampling_note}', fontsize=14, fontweight='bold')\n ax.legend()\n ax.grid(alpha=0.3)\n ax.set_xscale('log')\n ax.set_yscale('log')\n else:\n ax.text(0.5, 0.5, 'No degree data available', \n ha='center', va='center', transform=ax.transAxes)\n ax.set_title('Degree Distribution', fontsize=14, fontweight='bold')\n\ndef plot_hierarchy_stats_optimized(graph: nx.DiGraph, annotations: pd.DataFrame, ax):\n \"\"\"Optimized hierarchy statistics table.\"\"\"\n ax.axis('off')\n \n # Use cached stats or calculate quickly\n stats = analyze_go_hierarchy(graph, annotations)\n \n table_data = [\n ['Metric', 'Value'],\n ['Total Terms', f\"{stats['num_nodes']:,}\"],\n ['Total Relationships', f\"{stats['num_edges']:,}\"],\n ['Annotated Terms', f\"{stats['annotated_terms_in_graph']:,}\"],\n ['Missing Terms', f\"{stats['annotated_terms_not_in_graph']:,}\"],\n ]\n \n if 'term_depth' in stats and 'error' not in stats['term_depth']:\n depth_stats = stats['term_depth']\n table_data.extend([\n ['Mean Depth', f\"{depth_stats['mean']:.2f}\"],\n ['Max Depth', f\"{depth_stats['max']}\"],\n ['Depth Std', f\"{depth_stats['std']:.2f}\"],\n ['Terms Measured', f\"{depth_stats['num_measured_terms']:,}\"],\n ])\n \n if 'degree_analysis' in stats:\n deg = stats['degree_analysis']\n table_data.extend([\n ['Mean In-degree', f\"{deg['in_degree_mean']:.2f}\"],\n ['Median In-degree', f\"{deg['in_degree_median']:.1f}\"],\n ['Max In-degree', f\"{deg['in_degree_max']}\"],\n ['Mean Out-degree', f\"{deg['out_degree_mean']:.2f}\"],\n ])\n \n table = ax.table(cellText=table_data, \n cellLoc='center',\n loc='center',\n bbox=[0.1, 0.1, 0.8, 0.8])\n \n table.auto_set_font_size(False)\n table.set_fontsize(9) # Slightly smaller font\n table.scale(1, 1.8) # Slightly smaller scale\n \n # Style header row\n for i in range(len(table_data)):\n for j in range(len(table_data[0])):\n if i == 0: # Header row\n table[(i, j)].set_facecolor('#4B5563')\n table[(i, j)].set_text_props(weight='bold', color='white')\n \n ax.set_title('GO Hierarchy Statistics', fontsize=16, fontweight='bold', pad=20)\n\ndef plot_annotated_overview(graph: nx.DiGraph, annotations: pd.DataFrame, ax):\n \"\"\"Fast overview of annotated terms (replaces slow network plot).\"\"\"\n annotated_terms = set(annotations['go_term'].unique())\n annotated_in_graph = annotated_terms & set(graph.nodes())\n \n # Get ontology distribution of annotated terms\n if 'ontology' in annotations.columns:\n ont_data = annotations[annotations['go_term'].isin(annotated_in_graph)]\n ont_counts = ont_data['ontology'].value_counts()\n \n colors = {'BPO': '#EF4444', 'MFO': '#3B82F6', 'CCO': '#10B981'}\n color_list = [colors.get(ont, 'gray') for ont in ont_counts.index]\n \n bars = ax.bar(ont_counts.index, ont_counts.values, color=color_list, alpha=0.8)\n ax.set_xlabel('Ontology', fontsize=12)\n ax.set_ylabel('Number of Annotated Terms', fontsize=12)\n ax.set_title(f'Annotated Terms by Ontology\\n({len(annotated_in_graph):,} terms in graph)', \n fontsize=14, fontweight='bold')\n ax.grid(axis='y', alpha=0.3)\n \n # Add value labels on bars\n for bar, count in zip(bars, ont_counts.values):\n ax.text(bar.get_x() + bar.get_width()/2, bar.get_height(),\n f'{count:,}', ha='center', va='bottom', fontsize=10, fontweight='bold')\n else:\n ax.text(0.5, 0.5, f'{len(annotated_in_graph):,} annotated terms in graph', \n ha='center', va='center', transform=ax.transAxes, fontsize=12)\n ax.set_title('Annotated Terms Overview', fontsize=14, fontweight='bold')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T20:00:19.121953Z","iopub.execute_input":"2025-10-19T20:00:19.122224Z","iopub.status.idle":"2025-10-19T20:00:19.160784Z","shell.execute_reply.started":"2025-10-19T20:00:19.122203Z","shell.execute_reply":"2025-10-19T20:00:19.160045Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(\"Starting optimized analysis...\")\nstats = analyze_go_hierarchy(go_graph, annotations)\nprint(\"\\nGenerating optimized visualizations...\")\nvisualize_go_hierarchy(go_graph, annotations, \"/kaggle/working/go_hierarchy_analysis.png\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T20:00:19.161517Z","iopub.execute_input":"2025-10-19T20:00:19.161790Z","iopub.status.idle":"2025-10-19T20:00:32.499959Z","shell.execute_reply.started":"2025-10-19T20:00:19.161775Z","shell.execute_reply":"2025-10-19T20:00:32.499149Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def find_term_ancestors(graph: nx.DiGraph, term: str) -> Set[str]:\n \"\"\"Find all ancestor terms of a given GO term.\"\"\"\n try:\n return nx.ancestors(graph, term)\n except nx.NetworkXError:\n return set()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T20:00:32.500719Z","iopub.execute_input":"2025-10-19T20:00:32.500921Z","iopub.status.idle":"2025-10-19T20:00:32.505166Z","shell.execute_reply.started":"2025-10-19T20:00:32.500904Z","shell.execute_reply":"2025-10-19T20:00:32.504471Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport networkx as nx\nfrom typing import Set, List, Dict\nimport numpy as np\nfrom collections import deque\n\ndef visualize_term_ancestors(graph: nx.DiGraph, term: str, \n annotations: pd.DataFrame = None,\n max_depth: int = 5,\n save_path: str = None):\n \"\"\"\n Visualize the ancestor hierarchy for a specific GO term - FIXED VERSION.\n \n Args:\n graph: GO ontology graph\n term: GO term to analyze\n annotations: Optional annotations for highlighting\n max_depth: Maximum depth to display\n save_path: Optional path to save figure\n \"\"\"\n print(f\"Finding ancestors for {term}...\")\n ancestors = find_term_ancestors(graph, term)\n \n if not ancestors:\n print(f\"No ancestors found for term {term}\")\n # Still create a minimal visualization\n create_minimal_visualization(term, \"No ancestors found\", save_path)\n return\n \n print(f\"Found {len(ancestors)} ancestors. Creating visualization...\")\n \n # Create figure with multiple views\n fig, axes = plt.subplots(2, 2, figsize=(20, 16))\n \n # 1. Optimized ancestor tree\n ax1 = axes[0, 0]\n plot_ancestor_tree_fixed(graph, term, ancestors, ax1, max_depth)\n \n # 2. Ancestor statistics\n ax2 = axes[0, 1]\n plot_ancestor_stats_fixed(graph, term, ancestors, annotations, ax2)\n \n # 3. Depth distribution\n ax3 = axes[1, 0]\n plot_ancestor_depth_distribution_fixed(graph, term, ancestors, ax3)\n \n # 4. Ontology breakdown\n ax4 = axes[1, 1]\n plot_ancestor_ontology_breakdown_fixed(ancestors, annotations, ax4)\n \n plt.suptitle(f\"Ancestor Analysis for {term}\", fontsize=16, fontweight='bold')\n plt.tight_layout()\n \n if save_path:\n plt.savefig(save_path, dpi=300, bbox_inches='tight')\n print(f\"Saved ancestor visualization to {save_path}\")\n \n plt.show()\n\ndef create_minimal_visualization(term: str, message: str, save_path: str = None):\n \"\"\"Create a minimal visualization when no ancestors are found.\"\"\"\n fig, ax = plt.subplots(figsize=(10, 6))\n ax.text(0.5, 0.5, f\"Term: {term}\\n\\n{message}\", \n ha='center', va='center', transform=ax.transAxes, fontsize=14)\n ax.set_title(f\"Ancestor Analysis for {term}\", fontsize=16, fontweight='bold')\n ax.axis('off')\n \n if save_path:\n plt.savefig(save_path, dpi=300, bbox_inches='tight')\n plt.show()\n\ndef plot_ancestor_tree_fixed(graph: nx.DiGraph, term: str, ancestors: Set[str], \n ax, max_depth: int = 5):\n \"\"\"Fixed ancestor tree plot.\"\"\"\n # Create subgraph with term and its ancestors\n all_nodes = ancestors.union({term})\n \n if len(all_nodes) > 100:\n print(f\"Large ancestor set ({len(all_nodes)} nodes), sampling for visualization...\")\n # Get closest ancestors using BFS\n sampled_nodes = get_closest_ancestors_bfs(graph, term, max_nodes=100)\n subgraph = graph.subgraph(sampled_nodes)\n else:\n subgraph = graph.subgraph(all_nodes)\n \n if subgraph.number_of_nodes() <= 1:\n ax.text(0.5, 0.5, 'No connected ancestors to display', \n ha='center', va='center', transform=ax.transAxes)\n ax.set_title('Ancestor Tree', fontsize=14, fontweight='bold')\n return\n \n # Use spring layout (most reliable)\n try:\n pos = nx.spring_layout(subgraph, k=2, iterations=50, seed=42)\n except:\n # Fallback to random layout\n pos = nx.random_layout(subgraph)\n \n # Calculate node properties\n node_colors = []\n node_sizes = []\n \n for node in subgraph.nodes():\n if node == term:\n node_colors.append('red')\n node_sizes.append(500)\n else:\n node_colors.append('lightblue')\n node_sizes.append(200)\n \n # Draw the graph\n nx.draw_networkx_nodes(subgraph, pos, ax=ax, \n node_color=node_colors, node_size=node_sizes,\n alpha=0.8, edgecolors='black')\n \n # Draw edges if they exist\n if subgraph.number_of_edges() > 0:\n nx.draw_networkx_edges(subgraph, pos, ax=ax,\n edge_color='gray', arrows=True, arrowsize=15,\n alpha=0.6, arrowstyle='->')\n \n # Add labels for important nodes only\n labels = {}\n for node in subgraph.nodes():\n if node == term or subgraph.degree(node) <= 2: # Label term and leaf nodes\n labels[node] = node.split(':')[-1] # Show only the number part\n \n if labels:\n nx.draw_networkx_labels(subgraph, pos, labels, ax=ax,\n font_size=8, font_weight='bold')\n \n ax.set_title(f'Ancestor Tree\\n{len(subgraph.nodes())} of {len(ancestors) + 1} terms shown', \n fontsize=14, fontweight='bold')\n ax.axis('off')\n\ndef get_closest_ancestors_bfs(graph: nx.DiGraph, term: str, max_nodes: int = 100) -> Set[str]:\n \"\"\"Get closest ancestors using BFS.\"\"\"\n visited = {term}\n queue = deque([(term, 0)]) # (node, depth)\n closest_ancestors = set()\n \n while queue and len(visited) < max_nodes:\n current, depth = queue.popleft()\n \n for predecessor in graph.predecessors(current):\n if predecessor not in visited:\n visited.add(predecessor)\n closest_ancestors.add(predecessor)\n if len(visited) < max_nodes:\n queue.append((predecessor, depth + 1))\n \n return visited\n\ndef plot_ancestor_stats_fixed(graph: nx.DiGraph, term: str, ancestors: Set[str],\n annotations: pd.DataFrame, ax):\n \"\"\"Fixed ancestor statistics plot.\"\"\"\n ax.axis('off')\n \n # Basic statistics\n total_ancestors = len(ancestors)\n \n # Check if term exists in graph\n term_in_graph = term in graph\n \n # Find roots in ancestry\n root_terms = {'GO:0008150', 'GO:0003674', 'GO:0005575'}\n roots_in_ancestors = root_terms.intersection(ancestors.union({term}))\n \n # Calculate distances to roots\n distances_to_roots = {}\n if term_in_graph:\n reverse_graph = graph.reverse()\n for root in roots_in_ancestors:\n if root in reverse_graph:\n try:\n dist = nx.shortest_path_length(reverse_graph, root, term)\n distances_to_roots[root] = dist\n except nx.NetworkXNoPath:\n distances_to_roots[root] = 'No path'\n except:\n distances_to_roots[root] = 'Error'\n \n # Get ontology information\n ontology_info = \"Not available\"\n if annotations is not None and 'ontology' in annotations.columns:\n # Get ontology of the target term\n term_ontology = annotations[annotations['go_term'] == term]['ontology']\n if not term_ontology.empty:\n target_ontology = term_ontology.iloc[0]\n else:\n target_ontology = 'Unknown'\n \n # Get ontology distribution of ancestors\n ancestor_ontologies = annotations[annotations['go_term'].isin(ancestors)]\n if not ancestor_ontologies.empty:\n ontology_counts = ancestor_ontologies['ontology'].value_counts()\n ontology_info = f\"Target: {target_ontology}\\n\"\n for ont, count in ontology_counts.items():\n ontology_info += f\" {ont}: {count:,}\\n\"\n else:\n ontology_info = f\"Target: {target_ontology}\\nNo ancestor ontology data\"\n \n # Create statistics table\n stats_data = [\n ['Statistic', 'Value'],\n ['Target Term', term],\n ['Total Ancestors', f'{total_ancestors:,}'],\n ['Term in Graph', 'Yes' if term_in_graph else 'No'],\n ['Roots in Ancestry', f'{len(roots_in_ancestors)}'],\n ]\n \n # Add distance information\n for root, dist in distances_to_roots.items():\n root_name = {\n 'GO:0008150': 'Biological Process',\n 'GO:0003674': 'Molecular Function', \n 'GO:0005575': 'Cellular Component'\n }.get(root, root)\n stats_data.append([f'Distance to {root_name}', f'{dist}'])\n \n stats_data.extend([\n ['', ''],\n ['Ontology Information', ontology_info]\n ])\n \n # Create table\n table = ax.table(cellText=stats_data, \n cellLoc='left',\n loc='center',\n bbox=[0.1, 0.1, 0.8, 0.8])\n \n table.auto_set_font_size(False)\n table.set_fontsize(8) # Smaller font to fit more content\n table.scale(1, 1.2)\n \n # Style table\n for i in range(len(stats_data)):\n for j in range(len(stats_data[0])):\n if i == 0: # Header\n table[(i, j)].set_facecolor('#4B5563')\n table[(i, j)].set_text_props(weight='bold', color='white')\n elif stats_data[i][0] == 'Target Term':\n table[(i, j)].set_facecolor('#FEF3C7')\n \n ax.set_title('Ancestor Statistics', fontsize=14, fontweight='bold')\n\ndef plot_ancestor_depth_distribution_fixed(graph: nx.DiGraph, term: str, \n ancestors: Set[str], ax):\n \"\"\"Fixed depth distribution plot.\"\"\"\n if term not in graph:\n ax.text(0.5, 0.5, 'Term not in graph\\nCannot calculate depths', \n ha='center', va='center', transform=ax.transAxes)\n ax.set_title('Ancestor Depth Distribution', fontsize=14, fontweight='bold')\n return\n \n # Calculate depths using reverse graph (distance TO term)\n try:\n reverse_graph = graph.reverse()\n depths = {}\n \n for ancestor in ancestors:\n if ancestor in reverse_graph:\n try:\n depth = nx.shortest_path_length(reverse_graph, ancestor, term)\n depths[ancestor] = depth\n except nx.NetworkXNoPath:\n continue\n \n depth_values = list(depths.values())\n \n if depth_values:\n # Create histogram\n n, bins, patches = ax.hist(depth_values, bins=min(20, len(set(depth_values))), \n alpha=0.7, color='lightgreen', \n edgecolor='black', linewidth=0.5)\n \n # Add statistics lines\n ax.axvline(np.mean(depth_values), color='red', linestyle='--',\n label=f'Mean: {np.mean(depth_values):.1f}')\n ax.axvline(np.median(depth_values), color='orange', linestyle=':',\n label=f'Median: {np.median(depth_values):.1f}')\n \n ax.set_xlabel('Distance to Target Term', fontsize=12)\n ax.set_ylabel('Number of Ancestors', fontsize=12)\n ax.set_title(f'Ancestor Depth Distribution\\n({len(depth_values):,} measurable ancestors)', \n fontsize=14, fontweight='bold')\n ax.legend()\n ax.grid(alpha=0.3)\n \n # Add some statistics text\n stats_text = f'Min: {min(depth_values)}\\nMax: {max(depth_values)}\\nStd: {np.std(depth_values):.1f}'\n ax.text(0.95, 0.95, stats_text, transform=ax.transAxes, \n ha='right', va='top', bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))\n else:\n ax.text(0.5, 0.5, 'No measurable depths\\n(ancestors not connected to term)', \n ha='center', va='center', transform=ax.transAxes)\n ax.set_title('Ancestor Depth Distribution', fontsize=14, fontweight='bold')\n \n except Exception as e:\n ax.text(0.5, 0.5, f'Error calculating depths:\\n{str(e)}', \n ha='center', va='center', transform=ax.transAxes)\n ax.set_title('Ancestor Depth Distribution', fontsize=14, fontweight='bold')\n\ndef plot_ancestor_ontology_breakdown_fixed(ancestors: Set[str], \n annotations: pd.DataFrame, ax):\n \"\"\"Fixed ontology breakdown plot.\"\"\"\n if annotations is None or 'ontology' not in annotations.columns:\n ax.text(0.5, 0.5, 'No ontology data available', \n ha='center', va='center', transform=ax.transAxes, fontsize=12)\n ax.set_title('Ontology Breakdown', fontsize=14, fontweight='bold')\n return\n \n # Get ontology distribution of ancestors\n ancestor_ontologies = annotations[annotations['go_term'].isin(ancestors)]\n \n if ancestor_ontologies.empty:\n ax.text(0.5, 0.5, 'No ontology data\\nfor ancestors', \n ha='center', va='center', transform=ax.transAxes, fontsize=12)\n ax.set_title('Ontology Breakdown', fontsize=14, fontweight='bold')\n return\n \n ontology_counts = ancestor_ontologies['ontology'].value_counts()\n \n # Colors for ontologies\n colors = {'BPO': '#EF4444', 'MFO': '#3B82F6', 'CCO': '#10B981'}\n color_list = [colors.get(ont, 'gray') for ont in ontology_counts.index]\n \n # Create bar chart\n bars = ax.bar(range(len(ontology_counts)), ontology_counts.values,\n color=color_list, alpha=0.8, edgecolor='black')\n \n ax.set_xticks(range(len(ontology_counts)))\n ax.set_xticklabels([f\"{ont}\\n({count:,})\" for ont, count in ontology_counts.items()])\n ax.set_ylabel('Number of Ancestors', fontsize=12)\n ax.set_title(f'Ancestors by Ontology\\n({len(ancestor_ontologies):,} of {len(ancestors):,} ancestors)', \n fontsize=14, fontweight='bold')\n ax.grid(axis='y', alpha=0.3)\n \n # Add value labels on bars\n for bar, count in zip(bars, ontology_counts.values):\n height = bar.get_height()\n ax.text(bar.get_x() + bar.get_width()/2, height,\n f'{count:,}', ha='center', va='bottom', fontweight='bold')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T20:00:32.506008Z","iopub.execute_input":"2025-10-19T20:00:32.506264Z","iopub.status.idle":"2025-10-19T20:00:32.539253Z","shell.execute_reply.started":"2025-10-19T20:00:32.506243Z","shell.execute_reply":"2025-10-19T20:00:32.538655Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(\"Ancestor visualizations...\")\n\n# Single term analysis (fast)\nvisualize_term_ancestors(\n graph=go_graph,\n term=\"GO:0006915\",\n annotations=annotations,\n save_path=\"/kaggle/working/term_ancestors.png\"\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T20:00:32.539879Z","iopub.execute_input":"2025-10-19T20:00:32.540390Z","iopub.status.idle":"2025-10-19T20:00:39.911160Z","shell.execute_reply.started":"2025-10-19T20:00:32.540368Z","shell.execute_reply":"2025-10-19T20:00:39.910375Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def quick_ancestor_analysis(graph: nx.DiGraph, terms: List[str], \n annotations: pd.DataFrame = None,\n save_path: str = None):\n \"\"\"\n Quick analysis of ancestors for multiple terms - FIXED VERSION.\n \n Args:\n graph: GO ontology graph\n terms: List of GO terms to analyze\n annotations: Optional annotations DataFrame\n \"\"\"\n print(\"Quick Ancestor Analysis\")\n print(\"=\" * 50)\n \n results = []\n for term in terms:\n ancestors = find_term_ancestors(graph, term)\n \n # Get ontology if available\n if annotations is not None and 'ontology' in annotations.columns:\n term_ont = annotations[annotations['go_term'] == term]['ontology']\n ontology = term_ont.iloc[0] if not term_ont.empty else 'Unknown'\n else:\n ontology = 'Unknown'\n \n # FIXED: Check if term's hierarchy connects to roots\n roots_reached = check_roots_connected(graph, term)\n \n results.append({\n 'term': term,\n 'ontology': ontology,\n 'ancestor_count': len(ancestors),\n 'roots': roots_reached\n })\n \n # Create summary visualization\n fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(18, 6))\n \n # 1. Ancestor counts\n terms_list = [f\"{r['term']}\\n({r['ontology']})\" for r in results]\n counts = [r['ancestor_count'] for r in results]\n \n bars = ax1.bar(terms_list, counts, \n color=['#EF4444', '#3B82F6', '#10B981'][:len(terms)],\n alpha=0.8, edgecolor='black')\n ax1.set_ylabel('Number of Ancestors')\n ax1.set_title('Ancestor Counts by Term', fontsize=14, fontweight='bold')\n ax1.tick_params(axis='x', rotation=45)\n ax1.grid(axis='y', alpha=0.3)\n \n # Add value labels\n for bar, count in zip(bars, counts):\n ax1.text(bar.get_x() + bar.get_width()/2, bar.get_height(),\n f'{count:,}', ha='center', va='bottom', fontweight='bold')\n \n # 2. Roots reached - FIXED\n roots_data = [r['roots'] for r in results]\n root_names = ['BP', 'MF', 'CC'] # Simplified root names\n \n # Create stacked bar chart showing which roots are reached\n bottom = np.zeros(len(results))\n colors = ['#EF4444', '#3B82F6', '#10B981'] # BPO, MFO, CCO colors\n \n for i in range(3): # For each possible root\n root_counts = [1 if r['roots'] > i else 0 for r in results]\n ax2.bar(terms_list, root_counts, bottom=bottom, \n color=colors[i], alpha=0.8, edgecolor='black', \n label=root_names[i])\n bottom += root_counts\n \n ax2.set_ylabel('Roots Reached')\n ax2.set_title('Ontology Roots Reached', fontsize=14, fontweight='bold')\n ax2.tick_params(axis='x', rotation=45)\n ax2.grid(axis='y', alpha=0.3)\n ax2.legend(title='Roots')\n ax2.set_ylim(0, 3.5)\n \n # 3. Depth analysis\n depth_data = []\n for term in terms:\n try:\n # Calculate average depth from roots\n avg_depth = calculate_average_depth(graph, term)\n depth_data.append(avg_depth if avg_depth is not None else 0)\n except:\n depth_data.append(0)\n \n bars = ax3.bar(terms_list, depth_data, \n color='lightgreen', alpha=0.8, edgecolor='black')\n ax3.set_ylabel('Average Depth from Roots')\n ax3.set_title('Hierarchical Depth', fontsize=14, fontweight='bold')\n ax3.tick_params(axis='x', rotation=45)\n ax3.grid(axis='y', alpha=0.3)\n \n # Add value labels\n for bar, depth in zip(bars, depth_data):\n ax3.text(bar.get_x() + bar.get_width()/2, bar.get_height(),\n f'{depth:.1f}', ha='center', va='bottom', fontweight='bold')\n \n plt.tight_layout()\n if save_path:\n plt.savefig(save_path, dpi=300, bbox_inches='tight')\n print(f\"Saved comparison visualization to {save_path}\")\n plt.show()\n \n # Print detailed summary\n print(\"\\nDetailed Summary:\")\n print(\"-\" * 60)\n for result in results:\n print(f\"Term: {result['term']} ({result['ontology']})\")\n print(f\" Ancestors: {result['ancestor_count']:,}\")\n print(f\" Roots reached: {result['roots']}/3\")\n print(f\" Root details: {get_root_details(graph, result['term'])}\")\n print()\n\ndef check_roots_connected(graph: nx.DiGraph, term: str) -> int:\n \"\"\"\n Check how many of the three main roots are connected to the term.\n \n Returns:\n Number of roots (0-3) that are connected to the term\n \"\"\"\n root_terms = {\n 'GO:0008150', # biological_process\n 'GO:0003674', # molecular_function \n 'GO:0005575', # cellular_component\n }\n \n roots_connected = 0\n reverse_graph = graph.reverse()\n \n for root in root_terms:\n if root in reverse_graph:\n try:\n # Check if there's a path from root to term\n if nx.has_path(reverse_graph, root, term):\n roots_connected += 1\n except:\n continue\n \n return roots_connected\n\ndef get_root_details(graph: nx.DiGraph, term: str) -> str:\n \"\"\"Get detailed information about which roots are connected.\"\"\"\n root_mapping = {\n 'GO:0008150': 'BP',\n 'GO:0003674': 'MF',\n 'GO:0005575': 'CC'\n }\n \n connected_roots = []\n reverse_graph = graph.reverse()\n \n for root, short_name in root_mapping.items():\n if root in reverse_graph:\n try:\n if nx.has_path(reverse_graph, root, term):\n connected_roots.append(short_name)\n except:\n pass\n \n return \", \".join(connected_roots) if connected_roots else \"None\"\n\ndef calculate_average_depth(graph: nx.DiGraph, term: str) -> float:\n \"\"\"Calculate average depth of term from connected roots.\"\"\"\n root_terms = ['GO:0008150', 'GO:0003674', 'GO:0005575']\n depths = []\n reverse_graph = graph.reverse()\n \n for root in root_terms:\n if root in reverse_graph:\n try:\n depth = nx.shortest_path_length(reverse_graph, root, term)\n depths.append(depth)\n except:\n continue\n \n return np.mean(depths) if depths else None\n\n# Additional diagnostic function\ndef diagnose_term_roots(graph: nx.DiGraph, term: str):\n \"\"\"Diagnose root connectivity for a single term.\"\"\"\n print(f\"\\nDiagnosing root connectivity for {term}:\")\n print(\"-\" * 40)\n \n root_terms = {\n 'GO:0008150': 'Biological Process',\n 'GO:0003674': 'Molecular Function',\n 'GO:0005575': 'Cellular Component'\n }\n \n reverse_graph = graph.reverse()\n \n for root, description in root_terms.items():\n print(f\"\\n{description} ({root}):\")\n \n if root not in reverse_graph:\n print(\" Root not in graph\")\n continue\n \n if term not in reverse_graph:\n print(\" Term not in graph\")\n continue\n \n try:\n has_path = nx.has_path(reverse_graph, root, term)\n print(f\" Path exists: {has_path}\")\n \n if has_path:\n distance = nx.shortest_path_length(reverse_graph, root, term)\n print(f\" Distance: {distance}\")\n \n # Get the path\n path = nx.shortest_path(reverse_graph, root, term)\n print(f\" Path length: {len(path)} nodes\")\n if len(path) <= 8: # Show short paths\n print(f\" Path: {' -> '.join(path[:4])} ... -> {term}\")\n except nx.NetworkXNoPath:\n print(\" No path exists\")\n except Exception as e:\n print(f\" Error: {e}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T20:00:39.912007Z","iopub.execute_input":"2025-10-19T20:00:39.912258Z","iopub.status.idle":"2025-10-19T20:00:39.939405Z","shell.execute_reply.started":"2025-10-19T20:00:39.912239Z","shell.execute_reply":"2025-10-19T20:00:39.938774Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Test the fixed version with diagnostics\nprint(\"Testing root connectivity...\")\n\n# First, diagnose each term\ntest_terms = [\"GO:0005515\", \"GO:0009507\", \"GO:0000122\"]\nfor term in test_terms:\n diagnose_term_roots(go_graph, term)\n\n# Then run the fixed analysis\nquick_ancestor_analysis(\n go_graph, \n test_terms, \n annotations=annotations,\n save_path=\"/kaggle/working/ancestors_comparison.png\"\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T20:00:39.940383Z","iopub.execute_input":"2025-10-19T20:00:39.940651Z","iopub.status.idle":"2025-10-19T20:00:55.853368Z","shell.execute_reply.started":"2025-10-19T20:00:39.940634Z","shell.execute_reply":"2025-10-19T20:00:55.852689Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def find_term_descendants(graph: nx.DiGraph, term: str) -> Set[str]:\n \"\"\"Find all descendant terms of a given GO term.\"\"\"\n try:\n return nx.descendants(graph, term)\n except nx.NetworkXError:\n return set()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T20:00:55.854145Z","iopub.execute_input":"2025-10-19T20:00:55.854329Z","iopub.status.idle":"2025-10-19T20:00:55.858355Z","shell.execute_reply.started":"2025-10-19T20:00:55.854314Z","shell.execute_reply":"2025-10-19T20:00:55.857675Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport networkx as nx\nfrom typing import Set, List, Dict\nimport numpy as np\nfrom collections import deque, Counter\n\n# Missing helper functions that need to be defined\ndef find_term_descendants(graph: nx.DiGraph, term: str) -> Set[str]:\n \"\"\"Find all descendant terms of a given GO term.\"\"\"\n try:\n return nx.descendants(graph, term)\n except nx.NetworkXError:\n return set()\n\ndef calculate_descendant_depth_stats(graph: nx.DiGraph, term: str, descendants: Set[str]) -> Dict:\n \"\"\"Calculate depth statistics for descendants.\"\"\"\n try:\n depths = nx.single_source_shortest_path_length(graph, term)\n descendant_depths = [depth for node, depth in depths.items() if node in descendants]\n \n if descendant_depths:\n return {\n 'mean': np.mean(descendant_depths),\n 'max': max(descendant_depths),\n 'min': min(descendant_depths),\n 'std': np.std(descendant_depths),\n 'count': len(descendant_depths)\n }\n except:\n pass\n return {}\n\ndef calculate_branching_stats(graph: nx.DiGraph, term: str, descendants: Set[str]) -> Dict:\n \"\"\"Calculate branching statistics for the descendant tree.\"\"\"\n all_nodes = descendants.union({term})\n subgraph = graph.subgraph(all_nodes)\n \n out_degrees = [subgraph.out_degree(node) for node in all_nodes]\n leaves = sum(1 for degree in out_degrees if degree == 0)\n \n return {\n 'mean': np.mean(out_degrees) if out_degrees else 0,\n 'max': max(out_degrees) if out_degrees else 0,\n 'leaves': leaves,\n 'internal_nodes': len(all_nodes) - leaves\n }\n\ndef create_minimal_descendant_visualization(term: str, message: str, save_path: str = None):\n \"\"\"Create a minimal visualization when no descendants are found.\"\"\"\n fig, ax = plt.subplots(figsize=(10, 6))\n ax.text(0.5, 0.5, f\"Term: {term}\\n\\n{message}\", \n ha='center', va='center', transform=ax.transAxes, fontsize=14)\n ax.set_title(f\"Descendant Analysis for {term}\", fontsize=16, fontweight='bold')\n ax.axis('off')\n \n if save_path:\n plt.savefig(save_path, dpi=300, bbox_inches='tight')\n plt.show()\n\ndef visualize_term_descendants_optimized(graph: nx.DiGraph, term: str, \n annotations: pd.DataFrame = None,\n save_path: str = None):\n \"\"\"\n Optimized visualization for descendant analysis - only reliable, valuable plots.\n \n Args:\n graph: GO ontology graph\n term: GO term to analyze\n annotations: Optional annotations for highlighting\n save_path: Optional path to save figure\n \"\"\"\n print(f\"Finding descendants for {term}...\")\n descendants = find_term_descendants(graph, term)\n \n if not descendants:\n print(f\"No descendants found for term {term}\")\n create_minimal_descendant_visualization(term, \"No descendants found\", save_path)\n return\n \n print(f\"Found {len(descendants)} descendants. Creating optimized visualization...\")\n \n # Create figure with only reliable plots\n fig, axes = plt.subplots(2, 2, figsize=(18, 12))\n \n # 1. Statistics overview (always reliable)\n ax1 = axes[0, 0]\n plot_descendant_stats_optimized(graph, term, descendants, annotations, ax1)\n \n # 2. Depth distribution (reliable)\n ax2 = axes[0, 1]\n plot_descendant_depth_distribution_optimized(graph, term, descendants, ax2)\n \n # 3. Ontology and annotation analysis (reliable)\n ax3 = axes[1, 0]\n plot_descendant_ontology_analysis(descendants, annotations, ax3)\n \n # 4. Branching analysis (reliable)\n ax4 = axes[1, 1]\n plot_branching_analysis(graph, term, descendants, ax4)\n \n plt.suptitle(f\"Descendant Analysis for {term}\", fontsize=16, fontweight='bold')\n plt.tight_layout()\n \n if save_path:\n plt.savefig(save_path, dpi=300, bbox_inches='tight')\n print(f\"Saved optimized descendant visualization to {save_path}\")\n \n plt.show()\n\ndef plot_descendant_stats_optimized(graph: nx.DiGraph, term: str, descendants: Set[str],\n annotations: pd.DataFrame, ax):\n \"\"\"Optimized statistics plot with only reliable metrics.\"\"\"\n ax.axis('off')\n \n # Basic statistics (always reliable)\n total_descendants = len(descendants)\n \n # Depth statistics (reliable)\n depth_stats = calculate_descendant_depth_stats(graph, term, descendants)\n \n # Branching statistics (reliable)\n branching_stats = calculate_branching_stats(graph, term, descendants)\n \n # Get reliable ontology information\n ontology_info = get_reliable_ontology_info(descendants, annotations, term)\n \n # Annotation statistics (reliable)\n annotation_stats = get_reliable_annotation_stats(descendants, annotations)\n \n # Create statistics table with only reliable metrics\n stats_data = [\n ['Metric', 'Value'],\n ['Target Term', term],\n ['Total Descendants', f'{total_descendants:,}'],\n ['Max Depth', f\"{depth_stats.get('max', 'N/A')}\"],\n ['Mean Depth', f\"{depth_stats.get('mean', 'N/A'):.1f}\"],\n ['Branching Factor', f\"{branching_stats.get('mean', 'N/A'):.1f}\"],\n ['Leaf Nodes', f\"{branching_stats.get('leaves', 'N/A'):,}\"],\n ['Internal Nodes', f\"{branching_stats.get('internal_nodes', 'N/A'):,}\"],\n ]\n \n # Add ontology info if available\n if ontology_info != \"Not available\":\n stats_data.extend([['', ''], ['Ontology', ontology_info]])\n \n # Add annotation stats if available\n if annotation_stats != \"No annotation data\":\n stats_data.extend([['', ''], ['Annotations', annotation_stats]])\n \n # Create table\n table = ax.table(cellText=stats_data, \n cellLoc='left',\n loc='center',\n bbox=[0.1, 0.1, 0.8, 0.8])\n \n table.auto_set_font_size(False)\n table.set_fontsize(9)\n table.scale(1, 1.1)\n \n # Style table\n for i in range(len(stats_data)):\n for j in range(len(stats_data[0])):\n if i == 0: # Header\n table[(i, j)].set_facecolor('#4B5563')\n table[(i, j)].set_text_props(weight='bold', color='white')\n elif stats_data[i][0] == 'Target Term':\n table[(i, j)].set_facecolor('#FEF3C7')\n \n ax.set_title('Descendant Statistics', fontsize=14, fontweight='bold')\n\ndef plot_descendant_depth_distribution_optimized(graph: nx.DiGraph, term: str, descendants: Set[str], ax):\n \"\"\"Optimized depth distribution plot.\"\"\"\n try:\n depths = nx.single_source_shortest_path_length(graph, term)\n descendant_depths = [depth for node, depth in depths.items() if node in descendants and depth > 0]\n \n if descendant_depths and len(set(descendant_depths)) > 1:\n # Create histogram with optimal bin calculation\n unique_depths = len(set(descendant_depths))\n bins = min(15, unique_depths)\n \n n, bins, patches = ax.hist(descendant_depths, bins=bins, \n alpha=0.7, color='#3B82F6', \n edgecolor='black', linewidth=0.5)\n \n # Add statistics\n mean_depth = np.mean(descendant_depths)\n median_depth = np.median(descendant_depths)\n \n ax.axvline(mean_depth, color='red', linestyle='--', linewidth=2,\n label=f'Mean: {mean_depth:.1f}')\n ax.axvline(median_depth, color='orange', linestyle=':', linewidth=2,\n label=f'Median: {median_depth:.1f}')\n \n ax.set_xlabel('Distance from Target Term', fontsize=12)\n ax.set_ylabel('Number of Descendants', fontsize=12)\n ax.set_title(f'Depth Distribution\\n({len(descendant_depths):,} descendants)', \n fontsize=14, fontweight='bold')\n ax.legend()\n ax.grid(alpha=0.3)\n \n # Add statistics box\n stats_text = f'Min: {min(descendant_depths)}\\nMax: {max(descendant_depths)}\\nStd: {np.std(descendant_depths):.1f}'\n ax.text(0.95, 0.95, stats_text, transform=ax.transAxes, \n ha='right', va='top', fontsize=10,\n bbox=dict(boxstyle='round', facecolor='white', alpha=0.9))\n else:\n # Handle cases with insufficient depth variation\n if descendant_depths:\n ax.text(0.5, 0.7, f'All descendants at depth: {descendant_depths[0]}', \n ha='center', va='center', transform=ax.transAxes, fontsize=12)\n else:\n ax.text(0.5, 0.5, 'No depth data available', \n ha='center', va='center', transform=ax.transAxes, fontsize=12)\n ax.set_title('Depth Distribution', fontsize=14, fontweight='bold')\n ax.grid(alpha=0.3)\n \n except Exception as e:\n ax.text(0.5, 0.5, f'Error calculating depths', \n ha='center', va='center', transform=ax.transAxes, fontsize=12)\n ax.set_title('Depth Distribution', fontsize=14, fontweight='bold')\n ax.grid(alpha=0.3)\n\ndef plot_descendant_ontology_analysis(descendants: Set[str], annotations: pd.DataFrame, ax):\n \"\"\"Reliable ontology analysis plot.\"\"\"\n if annotations is None or 'ontology' not in annotations.columns:\n ax.text(0.5, 0.5, 'No ontology data available', \n ha='center', va='center', transform=ax.transAxes, fontsize=12)\n ax.set_title('Ontology Analysis', fontsize=14, fontweight='bold')\n ax.grid(alpha=0.3)\n return\n \n # Get descendant annotations\n descendant_data = annotations[annotations['go_term'].isin(descendants)]\n \n if descendant_data.empty:\n ax.text(0.5, 0.5, 'No ontology data\\nfor descendants', \n ha='center', va='center', transform=ax.transAxes, fontsize=12)\n ax.set_title('Ontology Analysis', fontsize=14, fontweight='bold')\n ax.grid(alpha=0.3)\n return\n \n # Ontology distribution\n ontology_counts = descendant_data['ontology'].value_counts()\n \n # Colors for ontologies\n colors = {'BPO': '#EF4444', 'MFO': '#3B82F6', 'CCO': '#10B981'}\n color_list = [colors.get(ont, 'gray') for ont in ontology_counts.index]\n \n # Create bar chart\n bars = ax.bar(ontology_counts.index, ontology_counts.values,\n color=color_list, alpha=0.8, edgecolor='black')\n \n ax.set_ylabel('Number of Terms', fontsize=12)\n ax.set_title(f'Ontology Distribution\\n({len(descendant_data):,} annotated descendants)', \n fontsize=14, fontweight='bold')\n ax.grid(axis='y', alpha=0.3)\n \n # Add value labels and percentages\n total = ontology_counts.sum()\n for i, (bar, count) in enumerate(zip(bars, ontology_counts.values)):\n percentage = (count / total) * 100\n ax.text(bar.get_x() + bar.get_width()/2, bar.get_height(),\n f'{count:,}\\n({percentage:.1f}%)', ha='center', va='bottom', \n fontsize=10, fontweight='bold')\n\ndef plot_branching_analysis(graph: nx.DiGraph, term: str, descendants: Set[str], ax):\n \"\"\"Reliable branching analysis plot.\"\"\"\n all_nodes = descendants.union({term})\n subgraph = graph.subgraph(all_nodes)\n \n # Calculate out-degrees (branching factors)\n out_degrees = [subgraph.out_degree(node) for node in all_nodes]\n \n if out_degrees:\n # Create histogram of branching factors\n unique_degrees = len(set(out_degrees))\n bins = min(10, unique_degrees)\n \n n, bins, patches = ax.hist(out_degrees, bins=bins, \n alpha=0.7, color='#10B981',\n edgecolor='black', linewidth=0.5)\n \n # Calculate statistics\n mean_degree = np.mean(out_degrees)\n max_degree = max(out_degrees)\n leaves = sum(1 for degree in out_degrees if degree == 0)\n leaf_percentage = (leaves / len(out_degrees)) * 100\n \n ax.axvline(mean_degree, color='red', linestyle='--', linewidth=2,\n label=f'Mean: {mean_degree:.1f}')\n \n ax.set_xlabel('Branching Factor (Out-degree)', fontsize=12)\n ax.set_ylabel('Number of Terms', fontsize=12)\n ax.set_title(f'Branching Analysis\\n{leaves:,} leaves ({leaf_percentage:.1f}%)', \n fontsize=14, fontweight='bold')\n ax.legend()\n ax.grid(alpha=0.3)\n \n # Add statistics\n stats_text = f'Max: {max_degree}\\nLeaves: {leaves:,}\\nTotal: {len(out_degrees):,}'\n ax.text(0.95, 0.95, stats_text, transform=ax.transAxes, \n ha='right', va='top', fontsize=10,\n bbox=dict(boxstyle='round', facecolor='white', alpha=0.9))\n else:\n ax.text(0.5, 0.5, 'No branching data', \n ha='center', va='center', transform=ax.transAxes, fontsize=12)\n ax.set_title('Branching Analysis', fontsize=14, fontweight='bold')\n ax.grid(alpha=0.3)\n\ndef get_reliable_ontology_info(descendants: Set[str], annotations: pd.DataFrame, term: str) -> str:\n \"\"\"Get reliable ontology information.\"\"\"\n if annotations is None or 'ontology' not in annotations.columns:\n return \"Not available\"\n \n # Get target term ontology\n term_ontology = \"Unknown\"\n term_ont_data = annotations[annotations['go_term'] == term]['ontology']\n if not term_ont_data.empty:\n term_ontology = term_ont_data.iloc[0]\n \n # Get descendant ontology distribution\n descendant_ontologies = annotations[annotations['go_term'].isin(descendants)]\n if not descendant_ontologies.empty:\n ontology_counts = descendant_ontologies['ontology'].value_counts()\n total_annotated = len(descendant_ontologies)\n \n info = f\"Target: {term_ontology}\\n\"\n for ont, count in ontology_counts.items():\n percentage = (count / total_annotated) * 100\n info += f\"{ont}: {count:,} ({percentage:.1f}%)\\n\"\n return info.strip()\n \n return f\"Target: {term_ontology}\\nNo descendant data\"\n\ndef get_reliable_annotation_stats(descendants: Set[str], annotations: pd.DataFrame) -> str:\n \"\"\"Get reliable annotation statistics.\"\"\"\n if annotations is None:\n return \"Not available\"\n \n descendant_annotations = annotations[annotations['go_term'].isin(descendants)]\n \n if descendant_annotations.empty:\n return \"No annotation data\"\n \n stats = []\n annotated_terms = descendant_annotations['go_term'].nunique()\n total_annotations = len(descendant_annotations)\n unique_proteins = descendant_annotations['protein_id'].nunique()\n \n stats.append(f\"Annotated terms: {annotated_terms:,}\")\n stats.append(f\"Total annotations: {total_annotations:,}\")\n stats.append(f\"Unique proteins: {unique_proteins:,}\")\n \n # Calculate coverage\n if len(descendants) > 0:\n coverage = (annotated_terms / len(descendants)) * 100\n stats.append(f\"Coverage: {coverage:.1f}%\")\n \n return \"\\n\".join(stats)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T20:00:55.859114Z","iopub.execute_input":"2025-10-19T20:00:55.859431Z","iopub.status.idle":"2025-10-19T20:00:55.894613Z","shell.execute_reply.started":"2025-10-19T20:00:55.859409Z","shell.execute_reply":"2025-10-19T20:00:55.894016Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"visualize_term_descendants_optimized(\n graph=go_graph,\n term=\"GO:0006915\", \n annotations=annotations,\n save_path=\"/kaggle/working/term_descendants.png\"\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T20:00:55.895307Z","iopub.execute_input":"2025-10-19T20:00:55.895587Z","iopub.status.idle":"2025-10-19T20:00:59.003985Z","shell.execute_reply.started":"2025-10-19T20:00:55.895572Z","shell.execute_reply":"2025-10-19T20:00:59.003194Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def compare_term_descendants_optimized(graph: nx.DiGraph, \n terms: List[str], \n annotations: pd.DataFrame = None,\n max_terms_display: int = 8,\n save_path: str = None) -> pd.DataFrame:\n \"\"\"\n Optimized comparison of GO term descendants with enhanced visualizations.\n \n Args:\n graph: GO ontology graph\n terms: List of GO terms to compare\n annotations: Optional annotations DataFrame for additional insights\n max_terms_display: Maximum number of terms to display in detail\n save_path: Optional path to save figure\n \n Returns:\n DataFrame with comparison statistics\n \"\"\"\n print(\"šŸ” OPTIMIZED TERM DESCENDANT COMPARISON\")\n print(\"=\" * 60)\n \n # Validate input terms\n valid_terms = [term for term in terms if term in graph]\n if len(valid_terms) != len(terms):\n invalid_terms = set(terms) - set(valid_terms)\n print(f\"āš ļø Warning: {len(invalid_terms)} terms not found in graph: {invalid_terms}\")\n \n if not valid_terms:\n print(\"āŒ No valid terms to analyze\")\n return pd.DataFrame()\n \n # Limit display terms for readability\n display_terms = valid_terms[:max_terms_display]\n if len(valid_terms) > max_terms_display:\n print(f\"šŸ“Š Displaying first {max_terms_display} of {len(valid_terms)} terms for clarity\")\n \n # Collect comprehensive statistics\n results = []\n for term in valid_terms:\n descendants = find_term_descendants(graph, term)\n \n # Get ontology information\n ontology = get_term_ontology(term, annotations)\n \n # Calculate comprehensive statistics\n depth_stats = calculate_descendant_depth_stats(graph, term, descendants)\n branching_stats = calculate_branching_stats(graph, term, descendants)\n connectivity_stats = calculate_connectivity_stats(graph, term, descendants)\n \n # Annotation statistics\n annotation_stats = calculate_annotation_stats(descendants, annotations)\n \n results.append({\n 'term': term,\n 'ontology': ontology,\n 'descendant_count': len(descendants),\n 'mean_depth': depth_stats.get('mean', 0),\n 'max_depth': depth_stats.get('max', 0),\n 'depth_std': depth_stats.get('std', 0),\n 'branching_factor': branching_stats.get('mean', 0),\n 'max_branching': branching_stats.get('max', 0),\n 'leaf_count': branching_stats.get('leaves', 0),\n 'internal_nodes': branching_stats.get('internal_nodes', 0),\n 'graph_density': connectivity_stats.get('density', 0),\n 'annotated_descendants': annotation_stats.get('annotated_count', 0),\n 'annotation_coverage': annotation_stats.get('coverage', 0),\n 'avg_annotations_per_term': annotation_stats.get('avg_annotations', 0)\n })\n \n # Create results DataFrame\n results_df = pd.DataFrame(results)\n \n # Create enhanced visualization\n create_comparison_visualization(results_df, display_terms, save_path)\n \n # Print comprehensive summary\n print_comparison_summary(results_df, valid_terms)\n \n return results_df\n\ndef get_term_ontology(term: str, annotations: pd.DataFrame) -> str:\n \"\"\"Get ontology for a term from annotations.\"\"\"\n if annotations is None or 'ontology' not in annotations.columns:\n return \"Unknown\"\n \n term_data = annotations[annotations['go_term'] == term]\n if not term_data.empty:\n return term_data['ontology'].iloc[0]\n return \"Unknown\"\n\ndef calculate_connectivity_stats(graph: nx.DiGraph, term: str, descendants: Set[str]) -> Dict:\n \"\"\"Calculate connectivity statistics for descendant subgraph.\"\"\"\n all_nodes = descendants.union({term})\n subgraph = graph.subgraph(all_nodes)\n \n stats = {}\n \n if subgraph.number_of_nodes() > 1:\n stats['density'] = nx.density(subgraph)\n \n # Calculate connectivity measures\n try:\n if nx.is_weakly_connected(subgraph):\n stats['diameter'] = nx.diameter(subgraph)\n else:\n # For disconnected graphs, calculate for largest component\n largest_cc = max(nx.weakly_connected_components(subgraph), key=len)\n if len(largest_cc) > 1:\n largest_subgraph = subgraph.subgraph(largest_cc)\n stats['diameter_largest_cc'] = nx.diameter(largest_subgraph)\n except:\n stats['diameter'] = 'N/A'\n \n return stats\n\ndef calculate_annotation_stats(descendants: Set[str], annotations: pd.DataFrame) -> Dict:\n \"\"\"Calculate annotation statistics for descendants.\"\"\"\n stats = {\n 'annotated_count': 0,\n 'coverage': 0,\n 'avg_annotations': 0\n }\n \n if annotations is None or len(descendants) == 0:\n return stats\n \n descendant_annotations = annotations[annotations['go_term'].isin(descendants)]\n \n if not descendant_annotations.empty:\n annotated_terms = descendant_annotations['go_term'].nunique()\n total_annotations = len(descendant_annotations)\n \n stats['annotated_count'] = annotated_terms\n stats['coverage'] = (annotated_terms / len(descendants)) * 100\n stats['avg_annotations'] = total_annotations / annotated_terms if annotated_terms > 0 else 0\n \n return stats\n\ndef create_comparison_visualization(results_df: pd.DataFrame, \n display_terms: List[str], \n save_path: str = None):\n \"\"\"Create enhanced comparison visualization.\"\"\"\n if results_df.empty:\n print(\"āŒ No data to visualize\")\n return\n \n # Filter to display terms only\n display_df = results_df[results_df['term'].isin(display_terms)].copy()\n \n # Create enhanced visualization layout\n fig = plt.figure(figsize=(20, 16))\n \n # Define grid layout\n gs = fig.add_gridspec(3, 3, height_ratios=[1, 1, 0.8])\n \n # 1. Main comparison metrics\n ax1 = fig.add_subplot(gs[0, 0])\n plot_descendant_comparison(display_df, ax1)\n \n # 2. Depth analysis\n ax2 = fig.add_subplot(gs[0, 1])\n plot_depth_analysis(display_df, ax2)\n \n # 3. Branching complexity\n ax3 = fig.add_subplot(gs[0, 2])\n plot_branching_analysis_comparison(display_df, ax3)\n \n # 4. Structural metrics\n ax4 = fig.add_subplot(gs[1, 0])\n plot_structural_metrics(display_df, ax4)\n \n # 5. Annotation coverage\n ax5 = fig.add_subplot(gs[1, 1])\n plot_annotation_coverage(display_df, ax5)\n \n # 6. Term relationships\n ax6 = fig.add_subplot(gs[1, 2])\n plot_term_relationships(display_df, ax6)\n \n # 7. Comprehensive statistics table\n ax7 = fig.add_subplot(gs[2, :])\n plot_comprehensive_statistics(display_df, ax7)\n \n plt.suptitle('GO Term Descendant Comparison Analysis', fontsize=18, fontweight='bold', y=0.98)\n plt.tight_layout()\n \n if save_path:\n plt.savefig(save_path, dpi=300, bbox_inches='tight', facecolor='white')\n print(f\"šŸ’¾ Comparison visualization saved to: {save_path}\")\n \n plt.show()\n\ndef plot_descendant_comparison(df: pd.DataFrame, ax):\n \"\"\"Plot main descendant comparison with enhanced styling.\"\"\"\n terms_display = [f\"{row['term']}\\n({row['ontology']})\" for _, row in df.iterrows()]\n counts = df['descendant_count'].values\n \n # Create gradient color based on count\n colors = plt.cm.viridis(np.linspace(0.2, 0.8, len(counts)))\n \n bars = ax.bar(terms_display, counts, color=colors, alpha=0.8, \n edgecolor='black', linewidth=1.5)\n \n ax.set_ylabel('Number of Descendants', fontsize=12, fontweight='bold')\n ax.set_title('Descendant Count Comparison', fontsize=14, fontweight='bold')\n ax.tick_params(axis='x', rotation=45, labelsize=10)\n ax.grid(axis='y', alpha=0.3, linestyle='--')\n \n # Enhanced value labels\n for bar, count in zip(bars, counts):\n height = bar.get_height()\n ax.text(bar.get_x() + bar.get_width()/2, height * 1.02,\n f'{count:,}', ha='center', va='bottom', \n fontsize=11, fontweight='bold',\n bbox=dict(boxstyle='round,pad=0.2', facecolor='white', alpha=0.8))\n \n # Add ranking\n for i, (bar, count) in enumerate(zip(bars, counts)):\n rank = f\"#{i+1}\"\n ax.text(bar.get_x() + bar.get_width()/2, -height * 0.1,\n rank, ha='center', va='top', fontsize=9, \n fontweight='bold', color='red')\n\ndef plot_depth_analysis(df: pd.DataFrame, ax):\n \"\"\"Plot comprehensive depth analysis.\"\"\"\n terms_display = [f\"{row['term']}\\n({row['ontology']})\" for _, row in df.iterrows()]\n \n x = np.arange(len(terms_display))\n width = 0.35\n \n # Mean and max depth bars\n bars1 = ax.bar(x - width/2, df['mean_depth'], width, \n label='Mean Depth', color='#FF6B6B', alpha=0.8)\n bars2 = ax.bar(x + width/2, df['max_depth'], width, \n label='Max Depth', color='#4ECDC4', alpha=0.8)\n \n ax.set_ylabel('Depth', fontsize=12, fontweight='bold')\n ax.set_title('Hierarchy Depth Analysis', fontsize=14, fontweight='bold')\n ax.set_xticks(x)\n ax.set_xticklabels(terms_display, rotation=45, ha='right', fontsize=9)\n ax.legend(fontsize=10)\n ax.grid(axis='y', alpha=0.3, linestyle='--')\n \n # Add value labels\n for bars, values in [(bars1, df['mean_depth']), (bars2, df['max_depth'])]:\n for bar, value in zip(bars, values):\n ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() * 1.02,\n f'{value:.1f}', ha='center', va='bottom', \n fontsize=9, fontweight='bold')\n\ndef plot_branching_analysis_comparison(df: pd.DataFrame, ax):\n \"\"\"Plot branching factor comparison with leaf information.\"\"\"\n terms_display = [f\"{row['term']}\\n({row['ontology']})\" for _, row in df.iterrows()]\n \n x = np.arange(len(terms_display))\n width = 0.6\n \n # Branching factor bars\n bars = ax.bar(x, df['branching_factor'], width, \n color='#45B7D1', alpha=0.8, label='Branching Factor')\n \n ax.set_ylabel('Average Branching Factor', fontsize=12, fontweight='bold')\n ax.set_title('Branching Complexity', fontsize=14, fontweight='bold')\n ax.set_xticks(x)\n ax.set_xticklabels(terms_display, rotation=45, ha='right', fontsize=9)\n ax.grid(axis='y', alpha=0.3, linestyle='--')\n \n # Add branching factor values\n for bar, value in zip(bars, df['branching_factor']):\n ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() * 1.02,\n f'{value:.2f}', ha='center', va='bottom', \n fontsize=10, fontweight='bold')\n \n # Add leaf percentage as text above bars\n for i, (_, row) in enumerate(df.iterrows()):\n leaf_pct = (row['leaf_count'] / row['descendant_count']) * 100\n ax.text(i, row['branching_factor'] * 1.15,\n f'Leaves: {leaf_pct:.1f}%', ha='center', va='bottom',\n fontsize=8, fontweight='bold', color='red',\n bbox=dict(boxstyle='round,pad=0.2', facecolor='yellow', alpha=0.7))\n\ndef plot_structural_metrics(df: pd.DataFrame, ax):\n \"\"\"Plot structural metrics comparison.\"\"\"\n metrics = ['Internal Nodes', 'Leaf Nodes']\n x = np.arange(len(df))\n width = 0.35\n \n bars1 = ax.bar(x - width/2, df['internal_nodes'], width, \n label='Internal Nodes', color='#96CEB4', alpha=0.8)\n bars2 = ax.bar(x + width/2, df['leaf_count'], width, \n label='Leaf Nodes', color='#FFEAA7', alpha=0.8)\n \n ax.set_ylabel('Number of Nodes', fontsize=12, fontweight='bold')\n ax.set_title('Node Type Distribution', fontsize=14, fontweight='bold')\n ax.set_xticks(x)\n ax.set_xticklabels([row['ontology'] for _, row in df.iterrows()], fontsize=10)\n ax.legend(fontsize=10)\n ax.grid(axis='y', alpha=0.3, linestyle='--')\n \n # Add ratio annotations\n for i, (_, row) in enumerate(df.iterrows()):\n total = row['internal_nodes'] + row['leaf_count']\n if total > 0:\n leaf_ratio = (row['leaf_count'] / total) * 100\n ax.text(i, max(row['internal_nodes'], row['leaf_count']) * 1.05,\n f'Leaf Ratio: {leaf_ratio:.1f}%', ha='center', va='bottom',\n fontsize=8, fontweight='bold')\n\ndef plot_annotation_coverage(df: pd.DataFrame, ax):\n \"\"\"Plot annotation coverage metrics.\"\"\"\n if 'annotation_coverage' not in df.columns:\n ax.text(0.5, 0.5, 'No annotation data available', \n ha='center', va='center', fontsize=12, transform=ax.transAxes)\n ax.set_title('Annotation Coverage', fontsize=14, fontweight='bold')\n return\n \n terms_display = [f\"{row['term']}\\n({row['ontology']})\" for _, row in df.iterrows()]\n \n # Create stacked bar chart for annotation metrics\n coverage = df['annotation_coverage'].fillna(0)\n missing_coverage = 100 - coverage\n \n bars1 = ax.bar(terms_display, coverage, label='Annotated', \n color='#6A0572', alpha=0.8)\n bars2 = ax.bar(terms_display, missing_coverage, bottom=coverage, \n label='Not Annotated', color='#AB83A1', alpha=0.5)\n \n ax.set_ylabel('Coverage (%)', fontsize=12, fontweight='bold')\n ax.set_title('Annotation Coverage', fontsize=14, fontweight='bold')\n ax.set_xticklabels(terms_display, rotation=45, ha='right', fontsize=9)\n ax.legend(fontsize=10)\n ax.set_ylim(0, 100)\n ax.grid(axis='y', alpha=0.3, linestyle='--')\n \n # Add coverage percentages\n for bar, cov in zip(bars1, coverage):\n if cov > 0:\n ax.text(bar.get_x() + bar.get_width()/2, cov/2,\n f'{cov:.1f}%', ha='center', va='center',\n fontsize=9, fontweight='bold', color='white')\n\ndef plot_term_relationships(df: pd.DataFrame, ax):\n \"\"\"Plot relationships between different metrics.\"\"\"\n if len(df) < 2:\n ax.text(0.5, 0.5, 'Insufficient data for relationship analysis', \n ha='center', va='center', fontsize=12, transform=ax.transAxes)\n ax.set_title('Term Relationships', fontsize=14, fontweight='bold')\n return\n \n # Create scatter plot of descendant count vs branching factor\n scatter = ax.scatter(df['descendant_count'], df['branching_factor'],\n c=df['mean_depth'], s=df['max_depth']*20, \n cmap='viridis', alpha=0.7, edgecolors='black')\n \n # Add term labels\n for _, row in df.iterrows():\n ax.annotate(row['term'], \n (row['descendant_count'], row['branching_factor']),\n xytext=(5, 5), textcoords='offset points', fontsize=8,\n bbox=dict(boxstyle='round,pad=0.2', facecolor='white', alpha=0.7))\n \n ax.set_xlabel('Number of Descendants', fontsize=12, fontweight='bold')\n ax.set_ylabel('Branching Factor', fontsize=12, fontweight='bold')\n ax.set_title('Structural Relationships', fontsize=14, fontweight='bold')\n ax.grid(alpha=0.3, linestyle='--')\n \n # Add colorbar\n cbar = plt.colorbar(scatter, ax=ax)\n cbar.set_label('Mean Depth', fontsize=10)\n\ndef plot_comprehensive_statistics(df: pd.DataFrame, ax):\n \"\"\"Plot comprehensive statistics table.\"\"\"\n ax.axis('off')\n \n # Prepare table data\n table_data = [['Term', 'Ontology', 'Descendants', 'Mean Depth', \n 'Max Depth', 'Branching', 'Leaves', 'Coverage%']]\n \n for _, row in df.iterrows():\n coverage = row.get('annotation_coverage', 0)\n table_data.append([\n row['term'],\n row['ontology'],\n f\"{row['descendant_count']:,}\",\n f\"{row['mean_depth']:.1f}\",\n f\"{row['max_depth']}\",\n f\"{row['branching_factor']:.2f}\",\n f\"{row['leaf_count']:,}\",\n f\"{coverage:.1f}%\" if pd.notna(coverage) else \"N/A\"\n ])\n \n # Create table\n table = ax.table(cellText=table_data, \n cellLoc='center',\n loc='center',\n bbox=[0, 0, 1, 1])\n \n table.auto_set_font_size(False)\n table.set_fontsize(9)\n table.scale(1, 1.8)\n \n # Style table\n for i in range(len(table_data)):\n for j in range(len(table_data[0])):\n if i == 0: # Header\n table[(i, j)].set_facecolor('#2E86AB')\n table[(i, j)].set_text_props(weight='bold', color='white', size=10)\n else:\n # Color code by ontology\n ontology = table_data[i][1]\n colors = {'BPO': '#FF6B6B', 'MFO': '#4ECDC4', 'CCO': '#45B7D1', 'Unknown': '#F7F7F7'}\n table[(i, j)].set_facecolor(colors.get(ontology, '#F7F7F7'))\n\ndef print_comparison_summary(df: pd.DataFrame, all_terms: List[str]):\n \"\"\"Print comprehensive comparison summary.\"\"\"\n print(\"\\nšŸ“Š COMPREHENSIVE COMPARISON SUMMARY\")\n print(\"=\" * 80)\n \n # Overall statistics\n print(f\"\\nšŸ“ˆ Overall Statistics:\")\n print(f\" ā€¢ Terms analyzed: {len(df)}\")\n print(f\" ā€¢ Total descendants: {df['descendant_count'].sum():,}\")\n print(f\" ā€¢ Average descendants per term: {df['descendant_count'].mean():.0f}\")\n print(f\" ā€¢ Depth range: {df['mean_depth'].min():.1f} - {df['mean_depth'].max():.1f}\")\n \n # Top performers\n max_descendants = df.loc[df['descendant_count'].idxmax()]\n max_depth = df.loc[df['max_depth'].idxmax()]\n max_branching = df.loc[df['branching_factor'].idxmax()]\n \n print(f\"\\nšŸ† Top Performers:\")\n print(f\" ā€¢ Most descendants: {max_descendants['term']} ({max_descendants['descendant_count']:,})\")\n print(f\" ā€¢ Deepest hierarchy: {max_depth['term']} (max depth: {max_depth['max_depth']})\")\n print(f\" ā€¢ Highest branching: {max_branching['term']} ({max_branching['branching_factor']:.2f})\")\n \n # Ontology breakdown\n if 'ontology' in df.columns:\n ontology_counts = df['ontology'].value_counts()\n print(f\"\\nšŸ”¬ Ontology Distribution:\")\n for ont, count in ontology_counts.items():\n ont_df = df[df['ontology'] == ont]\n avg_descendants = ont_df['descendant_count'].mean()\n print(f\" ā€¢ {ont}: {count} terms, avg {avg_descendants:.0f} descendants\")\n \n # Annotation statistics\n if 'annotation_coverage' in df.columns:\n coverage_stats = df['annotation_coverage'].describe()\n print(f\"\\nšŸ“ Annotation Coverage:\")\n print(f\" ā€¢ Mean coverage: {coverage_stats['mean']:.1f}%\")\n print(f\" ā€¢ Range: {coverage_stats['min']:.1f}% - {coverage_stats['max']:.1f}%\")\n \n print(f\"\\n{'=' * 80}\")\n\n# Quick usage examples:\ndef demonstrate_comparison_functions(graph: nx.DiGraph, annotations: pd.DataFrame):\n \"\"\"Demonstrate the optimized comparison function.\"\"\"\n \n # Example terms from different ontologies\n test_terms = [\n \"GO:0008150\", # biological_process\n \"GO:0003674\", # molecular_function\n \"GO:0005575\", # cellular_component\n \"GO:0006915\", # apoptosis\n \"GO:0005524\", # ATP binding\n \"GO:0005737\", # cytoplasm\n ]\n \n print(\"šŸ§Ŗ DEMONSTRATING OPTIMIZED TERM COMPARISON\")\n print(\"=\" * 60)\n \n # Run comparison\n results = compare_term_descendants_optimized(\n graph=graph,\n terms=test_terms,\n annotations=annotations,\n max_terms_display=6,\n save_path=\"/kaggle/working/term_comparison_optimized.png\"\n )\n \n return results\n\n# Main usage:\n# results = compare_term_descendants_optimized(go_graph, terms_list, annotations)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T20:13:34.697588Z","iopub.execute_input":"2025-10-19T20:13:34.697866Z","iopub.status.idle":"2025-10-19T20:13:34.743457Z","shell.execute_reply.started":"2025-10-19T20:13:34.697845Z","shell.execute_reply":"2025-10-19T20:13:34.742692Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Optimized comparison\ntest_terms = [\"GO:0005515\", \"GO:0006915\", \"GO:0000122\"]\ncompare_term_descendants_optimized(\n graph=go_graph,\n terms=test_terms,\n annotations=annotations,\n save_path=\"/kaggle/working/descendant_comparison_optimized.png\"\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T20:13:35.025802Z","iopub.execute_input":"2025-10-19T20:13:35.026441Z","iopub.status.idle":"2025-10-19T20:13:40.057531Z","shell.execute_reply.started":"2025-10-19T20:13:35.026390Z","shell.execute_reply":"2025-10-19T20:13:40.056982Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Taxonomy Analysis Functions","metadata":{}},{"cell_type":"code","source":"def analyze_taxonomy_distribution(taxonomy: pd.DataFrame) -> Dict:\n \"\"\"\n Analyze distribution of proteins across species.\n \n Args:\n taxonomy: DataFrame with protein_id and taxon_id\n \n Returns:\n Dictionary with taxonomy statistics\n \"\"\"\n stats = {}\n \n taxon_counts = taxonomy['taxon_id'].value_counts()\n \n stats['num_species'] = len(taxon_counts)\n stats['num_proteins'] = len(taxonomy)\n stats['proteins_per_species'] = {\n 'mean': taxon_counts.mean(),\n 'median': taxon_counts.median(),\n 'min': taxon_counts.min(),\n 'max': taxon_counts.max()\n }\n \n # Top species\n stats['top_species'] = taxon_counts.head(10).to_dict()\n \n print(\"\\nTaxonomy Distribution:\")\n print(f\"Total species: {stats['num_species']}\")\n print(f\"Total proteins: {stats['num_proteins']:,}\")\n print(f\"Proteins per species - Mean: {stats['proteins_per_species']['mean']:.1f}, \"\n f\"Median: {stats['proteins_per_species']['median']:.1f}\")\n \n print(\"\\nTop 10 Species by Protein Count:\")\n for taxon_id, count in list(stats['top_species'].items())[:10]:\n print(f\" Taxon {taxon_id}: {count:,} proteins\")\n \n return stats","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T20:14:13.539860Z","iopub.execute_input":"2025-10-19T20:14:13.540139Z","iopub.status.idle":"2025-10-19T20:14:13.545983Z","shell.execute_reply.started":"2025-10-19T20:14:13.540101Z","shell.execute_reply":"2025-10-19T20:14:13.545202Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport seaborn as sns\nimport numpy as np\nfrom typing import Dict\n\ndef visualize_taxonomy_distribution(taxonomy: pd.DataFrame, \n top_n: int = 20,\n save_path: str = None) -> Dict:\n \"\"\"\n Visualize distribution of proteins across species with comprehensive plots.\n \n Args:\n taxonomy: DataFrame with protein_id and taxon_id\n top_n: Number of top species to display in detail\n save_path: Optional path to save figure\n \n Returns:\n Dictionary with taxonomy statistics\n \"\"\"\n # Get statistics from the analysis function\n stats = analyze_taxonomy_distribution(taxonomy)\n \n # Create comprehensive visualization\n fig, axes = plt.subplots(2, 3, figsize=(20, 12))\n \n # 1. Top species bar chart\n ax1 = axes[0, 0]\n plot_top_species(stats, top_n, ax1)\n \n # 2. Species distribution histogram\n ax2 = axes[0, 1]\n plot_species_distribution(taxonomy, ax2)\n \n # 3. Cumulative distribution\n ax3 = axes[0, 2]\n plot_cumulative_distribution(taxonomy, ax3)\n \n # 4. Statistics overview\n ax4 = axes[1, 0]\n plot_taxonomy_stats(stats, ax4)\n \n # 5. Protein coverage by species\n ax5 = axes[1, 1]\n plot_protein_coverage(taxonomy, ax5)\n \n # 6. Species rank plot\n ax6 = axes[1, 2]\n plot_species_rank(taxonomy, ax6)\n \n plt.suptitle('Taxonomic Distribution Analysis', fontsize=16, fontweight='bold')\n plt.tight_layout()\n \n if save_path:\n plt.savefig(save_path, dpi=300, bbox_inches='tight')\n print(f\"Saved taxonomy visualization to {save_path}\")\n \n plt.show()\n \n return stats\n\ndef plot_top_species(stats: Dict, top_n: int, ax):\n \"\"\"Plot top N species by protein count.\"\"\"\n top_species = dict(list(stats['top_species'].items())[:top_n])\n \n if not top_species:\n ax.text(0.5, 0.5, 'No species data available', \n ha='center', va='center', transform=ax.transAxes)\n ax.set_title(f'Top {top_n} Species', fontsize=14, fontweight='bold')\n return\n \n species_ids = list(top_species.keys())\n counts = list(top_species.values())\n \n # Create bar chart\n bars = ax.bar(range(len(species_ids)), counts, \n color=plt.cm.viridis(np.linspace(0, 1, len(species_ids))),\n alpha=0.8, edgecolor='black')\n \n ax.set_xlabel('Species (Taxon ID)', fontsize=12)\n ax.set_ylabel('Number of Proteins', fontsize=12)\n ax.set_title(f'Top {top_n} Species by Protein Count', fontsize=14, fontweight='bold')\n ax.set_xticks(range(len(species_ids)))\n ax.set_xticklabels(species_ids, rotation=45, ha='right')\n ax.grid(axis='y', alpha=0.3)\n \n # Add value labels on bars\n for bar, count in zip(bars, counts):\n height = bar.get_height()\n ax.text(bar.get_x() + bar.get_width()/2, height,\n f'{count:,}', ha='center', va='bottom', \n fontsize=9, fontweight='bold')\n \n # Add percentage labels\n total_proteins = stats['num_proteins']\n for i, (bar, count) in enumerate(zip(bars, counts)):\n percentage = (count / total_proteins) * 100\n ax.text(bar.get_x() + bar.get_width()/2, height/2,\n f'{percentage:.1f}%', ha='center', va='center',\n fontsize=8, color='white', fontweight='bold')\n\ndef plot_species_distribution(taxonomy: pd.DataFrame, ax):\n \"\"\"Plot distribution of proteins per species.\"\"\"\n taxon_counts = taxonomy['taxon_id'].value_counts()\n \n if len(taxon_counts) == 0:\n ax.text(0.5, 0.5, 'No species data available', \n ha='center', va='center', transform=ax.transAxes)\n ax.set_title('Species Distribution', fontsize=14, fontweight='bold')\n return\n \n # Create histogram\n counts = taxon_counts.values\n \n # Use log scale for better visualization of long-tailed distributions\n if counts.max() / counts.min() > 100: # Large dynamic range\n bins = np.logspace(np.log10(counts.min()), np.log10(counts.max()), 30)\n ax.set_xscale('log')\n xlabel = 'Proteins per Species (log scale)'\n else:\n bins = 30\n xlabel = 'Proteins per Species'\n \n n, bins, patches = ax.hist(counts, bins=bins, alpha=0.7, \n color='skyblue', edgecolor='black')\n \n # Add statistics lines\n mean_count = counts.mean()\n median_count = np.median(counts)\n \n ax.axvline(mean_count, color='red', linestyle='--', linewidth=2,\n label=f'Mean: {mean_count:.1f}')\n ax.axvline(median_count, color='orange', linestyle=':', linewidth=2,\n label=f'Median: {median_count:.1f}')\n \n ax.set_xlabel(xlabel, fontsize=12)\n ax.set_ylabel('Number of Species', fontsize=12)\n ax.set_title('Distribution of Proteins per Species', fontsize=14, fontweight='bold')\n ax.legend()\n ax.grid(alpha=0.3)\n \n # Add statistics box\n stats_text = f'Total species: {len(counts):,}\\nMin: {counts.min():,}\\nMax: {counts.max():,}'\n ax.text(0.95, 0.95, stats_text, transform=ax.transAxes,\n ha='right', va='top', fontsize=10,\n bbox=dict(boxstyle='round', facecolor='white', alpha=0.9))\n\ndef plot_cumulative_distribution(taxonomy: pd.DataFrame, ax):\n \"\"\"Plot cumulative distribution of proteins across species.\"\"\"\n taxon_counts = taxonomy['taxon_id'].value_counts()\n \n if len(taxon_counts) == 0:\n ax.text(0.5, 0.5, 'No species data available', \n ha='center', va='center', transform=ax.transAxes)\n ax.set_title('Cumulative Distribution', fontsize=14, fontweight='bold')\n return\n \n # Calculate cumulative distribution\n sorted_counts = np.sort(taxon_counts.values)[::-1] # Descending\n cumulative = np.cumsum(sorted_counts)\n cumulative_percentage = (cumulative / cumulative[-1]) * 100\n species_percentage = (np.arange(len(sorted_counts)) + 1) / len(sorted_counts) * 100\n \n # Plot cumulative distribution\n ax.plot(species_percentage, cumulative_percentage, \n linewidth=3, color='purple', alpha=0.8)\n ax.fill_between(species_percentage, cumulative_percentage, alpha=0.3, color='purple')\n \n # Add reference lines\n ax.axhline(80, color='red', linestyle='--', alpha=0.7, label='80% of proteins')\n ax.axvline(20, color='orange', linestyle=':', alpha=0.7, label='20% of species')\n \n # Find intersection points\n idx_80_proteins = np.argmax(cumulative_percentage >= 80)\n idx_20_species = int(0.2 * len(species_percentage))\n \n if idx_80_proteins < len(species_percentage):\n ax.plot(species_percentage[idx_80_proteins], 80, 'ro', markersize=8)\n ax.text(species_percentage[idx_80_proteins], 85, \n f'{species_percentage[idx_80_proteins]:.1f}% species\\nhave 80% proteins',\n ha='center', va='bottom', fontsize=9)\n \n ax.set_xlabel('Percentage of Species (%)', fontsize=12)\n ax.set_ylabel('Cumulative Percentage of Proteins (%)', fontsize=12)\n ax.set_title('Cumulative Distribution of Proteins', fontsize=14, fontweight='bold')\n ax.legend()\n ax.grid(alpha=0.3)\n ax.set_xlim(0, 100)\n ax.set_ylim(0, 100)\n\ndef plot_taxonomy_stats(stats: Dict, ax):\n \"\"\"Plot taxonomy statistics overview.\"\"\"\n ax.axis('off')\n \n # Create comprehensive statistics table\n stats_data = [\n ['Metric', 'Value'],\n ['Total Species', f\"{stats['num_species']:,}\"],\n ['Total Proteins', f\"{stats['num_proteins']:,}\"],\n ['Mean Proteins/Species', f\"{stats['proteins_per_species']['mean']:.1f}\"],\n ['Median Proteins/Species', f\"{stats['proteins_per_species']['median']:.1f}\"],\n ['Min Proteins/Species', f\"{stats['proteins_per_species']['min']:,}\"],\n ['Max Proteins/Species', f\"{stats['proteins_per_species']['max']:,}\"],\n ]\n \n # Add coverage statistics\n taxon_counts = list(stats['top_species'].values())\n if taxon_counts:\n top_10_coverage = sum(list(stats['top_species'].values())[:10]) / stats['num_proteins'] * 100\n top_20_coverage = sum(list(stats['top_species'].values())[:20]) / stats['num_proteins'] * 100\n \n stats_data.extend([\n ['', ''],\n ['Coverage Statistics', ''],\n ['Top 10 Species Coverage', f'{top_10_coverage:.1f}%'],\n ['Top 20 Species Coverage', f'{top_20_coverage:.1f}%'],\n ])\n \n # Create table\n table = ax.table(cellText=stats_data, \n cellLoc='center',\n loc='center',\n bbox=[0.1, 0.1, 0.8, 0.8])\n \n table.auto_set_font_size(False)\n table.set_fontsize(10)\n table.scale(1, 1.2)\n \n # Style table\n for i in range(len(stats_data)):\n for j in range(len(stats_data[0])):\n if i == 0: # Header row\n table[(i, j)].set_facecolor('#4B5563')\n table[(i, j)].set_text_props(weight='bold', color='white')\n elif i >= len(stats_data) - 4: # Coverage statistics\n table[(i, j)].set_facecolor('#E5E7EB')\n \n ax.set_title('Taxonomy Statistics Overview', fontsize=14, fontweight='bold')\n\ndef plot_protein_coverage(taxonomy: pd.DataFrame, ax):\n \"\"\"Plot protein coverage by species rank.\"\"\"\n taxon_counts = taxonomy['taxon_id'].value_counts()\n \n if len(taxon_counts) == 0:\n ax.text(0.5, 0.5, 'No species data available', \n ha='center', va='center', transform=ax.transAxes)\n ax.set_title('Protein Coverage', fontsize=14, fontweight='bold')\n return\n \n sorted_counts = np.sort(taxon_counts.values)[::-1]\n cumulative_coverage = np.cumsum(sorted_counts) / np.sum(sorted_counts) * 100\n \n # Create area plot\n species_ranks = np.arange(1, len(sorted_counts) + 1)\n \n ax.fill_between(species_ranks, cumulative_coverage, alpha=0.6, color='green')\n ax.plot(species_ranks, cumulative_coverage, linewidth=2, color='darkgreen')\n \n # Add reference lines and annotations\n for coverage in [50, 80, 90, 95]:\n idx = np.argmax(cumulative_coverage >= coverage)\n if idx < len(species_ranks):\n ax.axvline(species_ranks[idx], color='red', linestyle='--', alpha=0.5)\n ax.text(species_ranks[idx], coverage + 2, \n f'{species_ranks[idx]} species\\n({coverage}% coverage)',\n ha='center', va='bottom', fontsize=8,\n bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))\n \n ax.set_xlabel('Species Rank', fontsize=12)\n ax.set_ylabel('Cumulative Protein Coverage (%)', fontsize=12)\n ax.set_title('Protein Coverage by Species Rank', fontsize=14, fontweight='bold')\n ax.grid(alpha=0.3)\n ax.set_xscale('log') # Log scale for better visualization\n\ndef plot_species_rank(taxonomy: pd.DataFrame, ax):\n \"\"\"Plot species rank vs protein count (Zipf-like distribution).\"\"\"\n taxon_counts = taxonomy['taxon_id'].value_counts()\n \n if len(taxon_counts) == 0:\n ax.text(0.5, 0.5, 'No species data available', \n ha='center', va='center', transform=ax.transAxes)\n ax.set_title('Species Rank Distribution', fontsize=14, fontweight='bold')\n return\n \n sorted_counts = np.sort(taxon_counts.values)[::-1]\n ranks = np.arange(1, len(sorted_counts) + 1)\n \n # Create rank plot\n ax.loglog(ranks, sorted_counts, 'o-', alpha=0.7, markersize=4, linewidth=1)\n \n # Add power law reference line if applicable\n if len(sorted_counts) > 10:\n try:\n # Fit a simple power law for reference\n log_ranks = np.log(ranks[:100]) # Use first 100 points for fit\n log_counts = np.log(sorted_counts[:100])\n slope, intercept = np.polyfit(log_ranks, log_counts, 1)\n \n # Plot fitted line\n fitted_counts = np.exp(intercept) * ranks**slope\n ax.loglog(ranks, fitted_counts, 'r--', alpha=0.8, \n label=f'Power law fit (Ī±={-slope:.2f})')\n ax.legend()\n except:\n pass\n \n ax.set_xlabel('Species Rank (log scale)', fontsize=12)\n ax.set_ylabel('Protein Count (log scale)', fontsize=12)\n ax.set_title('Species Rank vs Protein Count', fontsize=14, fontweight='bold')\n ax.grid(alpha=0.3, which='both')\n\n# Simplified version for quick analysis\ndef quick_taxonomy_visualization(taxonomy: pd.DataFrame, \n save_path: str = None) -> Dict:\n \"\"\"\n Quick taxonomy visualization with essential plots only.\n \n Args:\n taxonomy: DataFrame with protein_id and taxon_id\n save_path: Optional path to save figure\n \n Returns:\n Dictionary with taxonomy statistics\n \"\"\"\n stats = analyze_taxonomy_distribution(taxonomy)\n \n # Create simplified visualization\n fig, axes = plt.subplots(1, 2, figsize=(15, 6))\n \n # 1. Top species\n plot_top_species(stats, 15, axes[0])\n \n # 2. Distribution\n plot_species_distribution(taxonomy, axes[1])\n \n plt.suptitle('Taxonomic Distribution Overview', fontsize=16, fontweight='bold')\n plt.tight_layout()\n \n if save_path:\n plt.savefig(save_path, dpi=300, bbox_inches='tight')\n print(f\"Saved quick taxonomy visualization to {save_path}\")\n \n plt.show()\n \n return stats\n\n# Enhanced analysis function with visualization integration\ndef analyze_taxonomy_distribution_enhanced(taxonomy: pd.DataFrame, \n visualize: bool = True,\n save_path: str = None) -> Dict:\n \"\"\"\n Enhanced taxonomy analysis with optional visualization.\n \n Args:\n taxonomy: DataFrame with protein_id and taxon_id\n visualize: Whether to create visualizations\n save_path: Optional path to save figure\n \n Returns:\n Dictionary with taxonomy statistics\n \"\"\"\n \n if visualize:\n if len(taxonomy) > 1000: # Large dataset - use quick visualization\n quick_taxonomy_visualization(taxonomy, save_path)\n else: # Small dataset - use comprehensive visualization\n visualize_taxonomy_distribution(taxonomy, save_path=save_path)\n \n # return stats","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T20:14:26.810607Z","iopub.execute_input":"2025-10-19T20:14:26.811438Z","iopub.status.idle":"2025-10-19T20:14:26.843902Z","shell.execute_reply.started":"2025-10-19T20:14:26.811396Z","shell.execute_reply":"2025-10-19T20:14:26.843323Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"analyze_taxonomy_distribution_enhanced(taxonomy, visualize=True, save_path=\"/kaggle/working/taxonomy_distribution.png\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T20:14:26.844963Z","iopub.execute_input":"2025-10-19T20:14:26.845246Z","iopub.status.idle":"2025-10-19T20:14:29.191588Z","shell.execute_reply.started":"2025-10-19T20:14:26.845230Z","shell.execute_reply":"2025-10-19T20:14:29.190896Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def analyze_annotations_by_species(annotations: pd.DataFrame, taxonomy: pd.DataFrame) -> pd.DataFrame:\n \"\"\"\n Analyze annotation patterns across species.\n \n Returns:\n DataFrame with per-species annotation statistics\n \"\"\"\n # Merge annotations with taxonomy\n merged = annotations.merge(taxonomy, on='protein_id', how='left')\n \n # Group by species\n species_stats = merged.groupby('taxon_id').agg({\n 'protein_id': 'nunique',\n 'go_term': 'count'\n }).reset_index()\n \n species_stats.columns = ['taxon_id', 'num_proteins', 'num_annotations']\n species_stats['annotations_per_protein'] = (\n species_stats['num_annotations'] / species_stats['num_proteins']\n )\n \n species_stats = species_stats.sort_values('num_proteins', ascending=False)\n \n return species_stats","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T20:14:32.875255Z","iopub.execute_input":"2025-10-19T20:14:32.875518Z","iopub.status.idle":"2025-10-19T20:14:32.880455Z","shell.execute_reply.started":"2025-10-19T20:14:32.875497Z","shell.execute_reply":"2025-10-19T20:14:32.879766Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport seaborn as sns\nimport numpy as np\nfrom typing import Dict, List\nimport pandas as pd\n\ndef visualize_annotations_by_species(annotations: pd.DataFrame, \n taxonomy: pd.DataFrame,\n top_n: int = 20,\n save_path: str = None) -> pd.DataFrame:\n \"\"\"\n Visualize annotation patterns across species with comprehensive plots.\n \n Args:\n annotations: DataFrame with protein_id, go_term, ontology\n taxonomy: DataFrame with protein_id and taxon_id\n top_n: Number of top species to display in detail\n save_path: Optional path to save figure\n \n Returns:\n DataFrame with per-species annotation statistics\n \"\"\"\n # Get statistics from the analysis function\n species_stats = analyze_annotations_by_species(annotations, taxonomy)\n \n # Create comprehensive visualization\n fig, axes = plt.subplots(2, 3, figsize=(22, 14))\n \n # 1. Top species by protein count\n ax1 = axes[0, 0]\n plot_species_protein_counts(species_stats, top_n, ax1)\n \n # 2. Top species by annotation count\n ax2 = axes[0, 1]\n plot_species_annotation_counts(species_stats, top_n, ax2)\n \n # 3. Annotations per protein distribution\n ax3 = axes[0, 2]\n plot_annotations_per_protein(species_stats, ax3)\n \n # 4. Correlation analysis\n ax4 = axes[1, 0]\n plot_correlation_analysis(species_stats, ax4)\n \n # 5. Species statistics overview\n ax5 = axes[1, 1]\n plot_species_statistics(species_stats, ax5)\n \n # 6. Annotation density analysis\n ax6 = axes[1, 2]\n plot_annotation_density(annotations, taxonomy, species_stats, ax6)\n \n plt.suptitle('Annotation Patterns Across Species', fontsize=18, fontweight='bold')\n plt.tight_layout()\n \n if save_path:\n plt.savefig(save_path, dpi=300, bbox_inches='tight')\n print(f\"Saved species annotation visualization to {save_path}\")\n \n plt.show()\n \n return species_stats\n\ndef plot_species_protein_counts(species_stats: pd.DataFrame, top_n: int, ax):\n \"\"\"Plot top N species by protein count.\"\"\"\n top_species = species_stats.head(top_n)\n \n if top_species.empty:\n ax.text(0.5, 0.5, 'No species data available', \n ha='center', va='center', transform=ax.transAxes)\n ax.set_title(f'Top {top_n} Species by Protein Count', fontsize=14, fontweight='bold')\n return\n \n species_ids = top_species['taxon_id'].astype(str)\n protein_counts = top_species['num_proteins']\n \n # Create horizontal bar chart for better readability\n bars = ax.barh(range(len(species_ids)), protein_counts,\n color=plt.cm.plasma(np.linspace(0, 1, len(species_ids))),\n alpha=0.8, edgecolor='black')\n \n ax.set_xlabel('Number of Proteins', fontsize=12)\n ax.set_ylabel('Species (Taxon ID)', fontsize=12)\n ax.set_title(f'Top {top_n} Species by Protein Count', fontsize=14, fontweight='bold')\n ax.set_yticks(range(len(species_ids)))\n ax.set_yticklabels(species_ids, fontsize=9)\n ax.grid(axis='x', alpha=0.3)\n ax.invert_yaxis() # Highest count at top\n \n # Add value labels on bars\n for i, (bar, count) in enumerate(zip(bars, protein_counts)):\n width = bar.get_width()\n ax.text(width + width * 0.01, bar.get_y() + bar.get_height()/2,\n f'{count:,}', ha='left', va='center', \n fontsize=9, fontweight='bold')\n \n # Add percentage of total\n total_proteins = species_stats['num_proteins'].sum()\n for i, (bar, count) in enumerate(zip(bars, protein_counts)):\n percentage = (count / total_proteins) * 100\n ax.text(width * 0.3, bar.get_y() + bar.get_height()/2,\n f'{percentage:.1f}%', ha='center', va='center',\n fontsize=8, color='white', fontweight='bold')\n\ndef plot_species_annotation_counts(species_stats: pd.DataFrame, top_n: int, ax):\n \"\"\"Plot top N species by annotation count.\"\"\"\n top_species = species_stats.head(top_n)\n \n if top_species.empty:\n ax.text(0.5, 0.5, 'No species data available', \n ha='center', va='center', transform=ax.transAxes)\n ax.set_title(f'Top {top_n} Species by Annotation Count', fontsize=14, fontweight='bold')\n return\n \n species_ids = top_species['taxon_id'].astype(str)\n annotation_counts = top_species['num_annotations']\n \n # Create bar chart\n bars = ax.bar(range(len(species_ids)), annotation_counts,\n color=plt.cm.viridis(np.linspace(0, 1, len(species_ids))),\n alpha=0.8, edgecolor='black')\n \n ax.set_xlabel('Species (Taxon ID)', fontsize=12)\n ax.set_ylabel('Number of Annotations', fontsize=12)\n ax.set_title(f'Top {top_n} Species by Annotation Count', fontsize=14, fontweight='bold')\n ax.set_xticks(range(len(species_ids)))\n ax.set_xticklabels(species_ids, rotation=45, ha='right', fontsize=9)\n ax.grid(axis='y', alpha=0.3)\n \n # Add value labels on bars\n for bar, count in zip(bars, annotation_counts):\n height = bar.get_height()\n ax.text(bar.get_x() + bar.get_width()/2, height,\n f'{count:,}', ha='center', va='bottom', \n fontsize=9, fontweight='bold')\n \n # Add annotation density\n for i, (bar, count, proteins) in enumerate(zip(bars, annotation_counts, top_species['num_proteins'])):\n density = count / proteins\n ax.text(bar.get_x() + bar.get_width()/2, height * 0.3,\n f'{density:.1f} annot/prot', ha='center', va='center',\n fontsize=7, color='white', fontweight='bold',\n bbox=dict(boxstyle='round', facecolor='black', alpha=0.7))\n\ndef plot_annotations_per_protein(species_stats: pd.DataFrame, ax):\n \"\"\"Plot distribution of annotations per protein across species.\"\"\"\n if species_stats.empty:\n ax.text(0.5, 0.5, 'No species data available', \n ha='center', va='center', transform=ax.transAxes)\n ax.set_title('Annotations per Protein Distribution', fontsize=14, fontweight='bold')\n return\n \n annotations_per_protein = species_stats['annotations_per_protein']\n \n # Create histogram\n n, bins, patches = ax.hist(annotations_per_protein, bins=30, \n alpha=0.7, color='lightcoral',\n edgecolor='black', linewidth=0.5)\n \n # Add statistics lines\n mean_density = annotations_per_protein.mean()\n median_density = annotations_per_protein.median()\n \n ax.axvline(mean_density, color='red', linestyle='--', linewidth=2,\n label=f'Mean: {mean_density:.2f}')\n ax.axvline(median_density, color='orange', linestyle=':', linewidth=2,\n label=f'Median: {median_density:.2f}')\n \n ax.set_xlabel('Annotations per Protein', fontsize=12)\n ax.set_ylabel('Number of Species', fontsize=12)\n ax.set_title('Distribution of Annotation Density Across Species', \n fontsize=14, fontweight='bold')\n ax.legend()\n ax.grid(alpha=0.3)\n \n # Add statistics box\n stats_text = (f'Species: {len(annotations_per_protein):,}\\n'\n f'Min: {annotations_per_protein.min():.2f}\\n'\n f'Max: {annotations_per_protein.max():.2f}\\n'\n f'Std: {annotations_per_protein.std():.2f}')\n ax.text(0.95, 0.95, stats_text, transform=ax.transAxes,\n ha='right', va='top', fontsize=10,\n bbox=dict(boxstyle='round', facecolor='white', alpha=0.9))\n\ndef plot_correlation_analysis(species_stats: pd.DataFrame, ax):\n \"\"\"Plot correlation between protein count and annotation count.\"\"\"\n if len(species_stats) < 2:\n ax.text(0.5, 0.5, 'Insufficient data for correlation analysis', \n ha='center', va='center', transform=ax.transAxes)\n ax.set_title('Protein vs Annotation Correlation', fontsize=14, fontweight='bold')\n return\n \n # Create scatter plot\n scatter = ax.scatter(species_stats['num_proteins'], \n species_stats['num_annotations'],\n c=species_stats['annotations_per_protein'],\n cmap='viridis', alpha=0.7, s=50,\n edgecolors='black', linewidth=0.5)\n \n # Add colorbar\n cbar = plt.colorbar(scatter, ax=ax)\n cbar.set_label('Annotations per Protein', fontsize=10)\n \n # Calculate correlation\n correlation = species_stats['num_proteins'].corr(species_stats['num_annotations'])\n \n # Add trend line\n if not species_stats['num_proteins'].isna().any() and not species_stats['num_annotations'].isna().any():\n z = np.polyfit(species_stats['num_proteins'], species_stats['num_annotations'], 1)\n p = np.poly1d(z)\n ax.plot(species_stats['num_proteins'], p(species_stats['num_proteins']), \n \"r--\", alpha=0.8, linewidth=2,\n label=f'Correlation: {correlation:.3f}')\n \n ax.set_xlabel('Number of Proteins', fontsize=12)\n ax.set_ylabel('Number of Annotations', fontsize=12)\n ax.set_title(f'Protein vs Annotation Correlation\\n(r = {correlation:.3f})', \n fontsize=14, fontweight='bold')\n ax.legend()\n ax.grid(alpha=0.3)\n \n # Use log scale if data range is large\n if (species_stats['num_proteins'].max() / species_stats['num_proteins'].min() > 100 or\n species_stats['num_annotations'].max() / species_stats['num_annotations'].min() > 100):\n ax.set_xscale('log')\n ax.set_yscale('log')\n ax.set_xlabel('Number of Proteins (log scale)')\n ax.set_ylabel('Number of Annotations (log scale)')\n\ndef plot_species_statistics(species_stats: pd.DataFrame, ax):\n \"\"\"Plot comprehensive species statistics overview.\"\"\"\n ax.axis('off')\n \n if species_stats.empty:\n ax.text(0.5, 0.5, 'No species statistics available', \n ha='center', va='center', transform=ax.transAxes)\n return\n \n # Calculate comprehensive statistics\n total_species = len(species_stats)\n total_proteins = species_stats['num_proteins'].sum()\n total_annotations = species_stats['num_annotations'].sum()\n \n # Coverage statistics\n top_10_protein_coverage = species_stats.head(10)['num_proteins'].sum() / total_proteins * 100\n top_10_annotation_coverage = species_stats.head(10)['num_annotations'].sum() / total_annotations * 100\n \n # Density statistics\n density_stats = species_stats['annotations_per_protein'].describe()\n \n # Create statistics table\n stats_data = [\n ['Metric', 'Value'],\n ['Total Species', f'{total_species:,}'],\n ['Total Proteins', f'{total_proteins:,}'],\n ['Total Annotations', f'{total_annotations:,}'],\n ['', ''],\n ['Coverage Statistics', ''],\n ['Top 10 Species - Protein Coverage', f'{top_10_protein_coverage:.1f}%'],\n ['Top 10 Species - Annotation Coverage', f'{top_10_annotation_coverage:.1f}%'],\n ['', ''],\n ['Annotation Density', ''],\n ['Mean Annotations/Protein', f'{density_stats[\"mean\"]:.2f}'],\n ['Median Annotations/Protein', f'{density_stats[\"50%\"]:.2f}'],\n ['Min Annotations/Protein', f'{density_stats[\"min\"]:.2f}'],\n ['Max Annotations/Protein', f'{density_stats[\"max\"]:.2f}'],\n ]\n \n # Create table\n table = ax.table(cellText=stats_data, \n cellLoc='left',\n loc='center',\n bbox=[0.1, 0.1, 0.8, 0.8])\n \n table.auto_set_font_size(False)\n table.set_fontsize(9)\n table.scale(1, 1.1)\n \n # Style table\n for i in range(len(stats_data)):\n for j in range(len(stats_data[0])):\n if i == 0: # Header row\n table[(i, j)].set_facecolor('#4B5563')\n table[(i, j)].set_text_props(weight='bold', color='white')\n elif i in [5, 9]: # Section headers\n table[(i, j)].set_facecolor('#E5E7EB')\n table[(i, j)].set_text_props(weight='bold')\n \n ax.set_title('Species Annotation Statistics Overview', fontsize=14, fontweight='bold')\n\ndef plot_annotation_density(annotations: pd.DataFrame, taxonomy: pd.DataFrame, \n species_stats: pd.DataFrame, ax):\n \"\"\"Plot annotation density analysis across species.\"\"\"\n if species_stats.empty:\n ax.text(0.5, 0.5, 'No species data available', \n ha='center', va='center', transform=ax.transAxes)\n ax.set_title('Annotation Density Analysis', fontsize=14, fontweight='bold')\n return\n \n # Calculate density percentiles\n density_percentiles = np.percentile(species_stats['annotations_per_protein'], \n [25, 50, 75, 90, 95])\n \n # Create density categories\n low_density = species_stats[species_stats['annotations_per_protein'] <= density_percentiles[0]]\n medium_density = species_stats[(species_stats['annotations_per_protein'] > density_percentiles[0]) & \n (species_stats['annotations_per_protein'] <= density_percentiles[2])]\n high_density = species_stats[species_stats['annotations_per_protein'] > density_percentiles[2]]\n \n # Create stacked bar chart\n categories = ['Low Density', 'Medium Density', 'High Density']\n counts = [len(low_density), len(medium_density), len(high_density)]\n colors = ['#EF4444', '#F59E0B', '#10B981']\n \n bars = ax.bar(categories, counts, color=colors, alpha=0.8, edgecolor='black')\n \n ax.set_ylabel('Number of Species', fontsize=12)\n ax.set_title('Species by Annotation Density Category', fontsize=14, fontweight='bold')\n ax.grid(axis='y', alpha=0.3)\n \n # Add value labels and percentages\n total_species = len(species_stats)\n for bar, count in zip(bars, counts):\n height = bar.get_height()\n percentage = (count / total_species) * 100\n ax.text(bar.get_x() + bar.get_width()/2, height,\n f'{count:,}\\n({percentage:.1f}%)', ha='center', va='bottom',\n fontsize=10, fontweight='bold')\n \n # Add density range annotations\n ax.text(0.02, 0.98, f'Density Ranges:\\n'\n f'Low: ā‰¤{density_percentiles[0]:.2f}\\n'\n f'Medium: {density_percentiles[0]:.2f}-{density_percentiles[2]:.2f}\\n'\n f'High: >{density_percentiles[2]:.2f}',\n transform=ax.transAxes, ha='left', va='top', fontsize=9,\n bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))\n\n# Quick visualization version\ndef quick_species_annotation_visualization(annotations: pd.DataFrame,\n taxonomy: pd.DataFrame,\n top_n: int = 15,\n save_path: str = None) -> pd.DataFrame:\n \"\"\"\n Quick visualization of species annotation patterns.\n \n Args:\n annotations: DataFrame with protein_id, go_term, ontology\n taxonomy: DataFrame with protein_id and taxon_id\n top_n: Number of top species to display\n save_path: Optional path to save figure\n \n Returns:\n DataFrame with per-species annotation statistics\n \"\"\"\n species_stats = analyze_annotations_by_species(annotations, taxonomy)\n \n # Create simplified visualization\n fig, axes = plt.subplots(1, 2, figsize=(16, 6))\n \n # 1. Top species by protein count\n plot_species_protein_counts(species_stats, top_n, axes[0])\n \n # 2. Correlation analysis\n plot_correlation_analysis(species_stats, axes[1])\n \n plt.suptitle('Species Annotation Patterns - Quick Overview', fontsize=16, fontweight='bold')\n plt.tight_layout()\n \n if save_path:\n plt.savefig(save_path, dpi=300, bbox_inches='tight')\n print(f\"Saved quick species annotation visualization to {save_path}\")\n \n plt.show()\n \n return species_stats\n\n# Enhanced analysis with visualization\ndef analyze_annotations_by_species_enhanced(annotations: pd.DataFrame,\n taxonomy: pd.DataFrame,\n visualize: bool = True,\n top_n: int = 20,\n save_path: str = None) -> pd.DataFrame:\n \"\"\"\n Enhanced species annotation analysis with optional visualization.\n \n Args:\n annotations: DataFrame with protein_id, go_term, ontology\n taxonomy: DataFrame with protein_id and taxon_id\n visualize: Whether to create visualizations\n top_n: Number of top species to display\n save_path: Optional path to save figure\n \n Returns:\n DataFrame with per-species annotation statistics\n \"\"\"\n species_stats = analyze_annotations_by_species(annotations, taxonomy)\n \n if visualize:\n if len(species_stats) > 50: # Large dataset - use quick visualization\n quick_species_annotation_visualization(annotations, taxonomy, top_n, save_path)\n else: # Small dataset - use comprehensive visualization\n visualize_annotations_by_species(annotations, taxonomy, top_n, save_path)\n \n return species_stats","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T20:14:48.946506Z","iopub.execute_input":"2025-10-19T20:14:48.946983Z","iopub.status.idle":"2025-10-19T20:14:48.982422Z","shell.execute_reply.started":"2025-10-19T20:14:48.946960Z","shell.execute_reply":"2025-10-19T20:14:48.981928Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"analyze_annotations_by_species_enhanced(annotations, taxonomy, visualize=True, save_path=\"/kaggle/working/annotations_by_species.png\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-19T20:14:48.983753Z","iopub.execute_input":"2025-10-19T20:14:48.983943Z","iopub.status.idle":"2025-10-19T20:14:52.643261Z","shell.execute_reply.started":"2025-10-19T20:14:48.983928Z","shell.execute_reply":"2025-10-19T20:14:52.642581Z"}},"outputs":[],"execution_count":null}]}