{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":116062,"databundleVersionId":14875579,"sourceType":"competition"}],"dockerImageVersionId":31192,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"id":"44d0c3f7","cell_type":"markdown","source":"# CAFA6 Protein Function Prediction - Exploratory Data Analysis\n\n**Competition:** CAFA 6 Protein Function Prediction on Kaggle  \n**Team Name:** Overfit_INT3405E2  \n**Members:**\n- Chung Thị Mai Anh - 23021460\n- Nguyễn Trọng Đại - 24021400  \n- Nguyễn Thị Phương - 23021668\n\n**Best Score:** 0.283 (Public Score)\n\n---\n\n## Table of Contents\n1. [Problem Overview](#1-problem-overview)\n2. [Data Loading & Initial Exploration](#2-data-loading)\n3. [Dataset Structure Analysis](#3-dataset-structure)\n4. [GO Term Analysis](#4-go-term-analysis)\n5. [Protein Sequence Analysis](#5-sequence-analysis)\n6. [Taxonomy Analysis](#6-taxonomy-analysis)\n\n---\n\n<a id='1-problem-overview'></a>\n## 1. Problem Overview\n\nCAFA (Critical Assessment of Functional Annotation) là cuộc thi quốc tế đánh giá khả năng dự đoán chức năng protein. CAFA6 yêu cầu dự đoán các GO terms (Gene Ontology) cho các protein sequences dựa trên:\n\n- **Input:** Protein amino acid sequences\n- **Output:** GO terms với probability scores\n- **Evaluation Metric:** Weighted F1-score với Information Accretion (IA) weights\n\n**Three GO Ontologies:**\n- **Biological Process (P):** Các quá trình sinh học\n- **Molecular Function (F):** Chức năng phân tử\n- **Cellular Component (C):** Vị trí trong tế bào\n\n**Challenges:**\n- ~40,000+ GO terms (extremely high dimensionality)\n- Hierarchical structure (directed acyclic graph)\n- Severe class imbalance\n- Missing labels in training data","metadata":{}},{"id":"cf9af4a4-360a-4065-acc6-dea7e87f6e7b","cell_type":"code","source":"!pip install biopython obonet --quiet","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-14T15:54:57.620181Z","iopub.execute_input":"2025-12-14T15:54:57.620438Z","iopub.status.idle":"2025-12-14T15:55:03.768703Z","shell.execute_reply.started":"2025-12-14T15:54:57.620416Z","shell.execute_reply":"2025-12-14T15:55:03.767775Z"}},"outputs":[],"execution_count":null},{"id":"d508ff14","cell_type":"code","source":"# Data manipulation\nimport pandas as pd\nimport numpy as np\nfrom collections import Counter, defaultdict\nimport warnings\nwarnings.filterwarnings('ignore')\n\n# Visualization\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport plotly.express as px\nimport plotly.graph_objects as go\nfrom plotly.subplots import make_subplots\n\n# Bioinformatics\nfrom Bio import SeqIO\nimport obonet\nimport networkx as nx\n\n# ML libraries\nfrom sklearn.preprocessing import MultiLabelBinarizer\nfrom sklearn.model_selection import train_test_split\n\n# System\nimport os\nfrom pathlib import Path\nimport gc\n\n# Styling\nplt.style.use('seaborn-v0_8-darkgrid')\nsns.set_palette(\"husl\")\npd.set_option('display.max_columns', None)\npd.set_option('display.max_rows', 100)\n\nprint(f\"Pandas version: {pd.__version__}\")\nprint(f\"NumPy version: {np.__version__}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-14T15:55:11.758763Z","iopub.execute_input":"2025-12-14T15:55:11.759050Z","iopub.status.idle":"2025-12-14T15:55:14.001137Z","shell.execute_reply.started":"2025-12-14T15:55:11.759027Z","shell.execute_reply":"2025-12-14T15:55:14.000309Z"}},"outputs":[],"execution_count":null},{"id":"af10c1ad","cell_type":"markdown","source":"<a id='2-data-loading'></a>\n## 2. Data Loading & Initial Exploration\n\nLoad input data:\n- `train_sequences.fasta` - Protein sequences for training.\n- `train_terms.tsv` - GO term annotations\n- `train_taxonomy.tsv` - Taxonomy information\n- `go-basic.obo` - GO ontology structure\n- `IA.tsv` - Information Accretion weights\n- `testsuperset.fasta` - Test sequences\n- `sample_submission.tsv` - Submission format\n- `testsuperset-taxon-list.tsv` - Taxonomy IDs for the test set proteins.","metadata":{}},{"id":"dfd0e72b","cell_type":"code","source":"# Define paths\nBASE_PATH = Path(\"/kaggle/input/cafa-6-protein-function-prediction\")\nTRAIN_PATH = BASE_PATH / \"Train\"\nTEST_PATH = BASE_PATH / \"Test\"\n\nprint(f\"Train folder exists: {TRAIN_PATH.exists()}\")\nprint(f\"Test folder exists: {TEST_PATH.exists()}\")\n\n# 1. Load terms\ntrain_terms = pd.read_csv(\n    TRAIN_PATH / \"train_terms.tsv\",\n    sep=\"\\t\",\n    names=[\"EntryID\", \"term\", \"aspect\"],\n    header=None\n)\n# check if first row has 'EntryID' as value\nif train_terms.iloc[0]['EntryID'] == 'EntryID':\n    train_terms = train_terms.iloc[1:].reset_index(drop=True)\nprint(f\"train_terms: {train_terms.shape}\")\n\n# 2. load train taxonomy\ntrain_taxonomy = pd.read_csv(\n    TRAIN_PATH / \"train_taxonomy.tsv\",\n    sep=\"\\t\",\n    names=[\"EntryID\", \"taxon_id\"],\n    header=None\n)\n# check if first row has 'EntryID' as value\nif train_taxonomy.iloc[0]['EntryID'] == 'EntryID':\n    train_taxonomy = train_taxonomy.iloc[1:].reset_index(drop=True)\nprint(f\"train_taxonomy: {train_taxonomy.shape}\")\n\n# 3. Load IA weights\nia_weights = pd.read_csv(\n    BASE_PATH / \"IA.tsv\",\n    sep=\"\\t\",\n    names=[\"term\", \"IA\"],\n    header=None\n)\n# check if first row has 'term' as value\nif ia_weights.iloc[0]['term'] == 'term':\n    ia_weights = ia_weights.iloc[1:].reset_index(drop=True)\n# Convert IA to float\nia_weights['IA'] = pd.to_numeric(ia_weights['IA'], errors='coerce')\nprint(f\"IA weights: {ia_weights.shape}\")\n\n# 4. Load sample submission\ntry:\n    sample_submission = pd.read_csv(\n        BASE_PATH / \"sample_submission.tsv\",\n        sep=\"\\t\",\n        names=[\"EntryID\", \"term\", \"probability\"],\n        header=None,\n        on_bad_lines='skip'\n    )\n    if sample_submission.iloc[0]['EntryID'] == 'EntryID':\n        sample_submission = sample_submission.iloc[1:].reset_index(drop=True)\n    print(f\"sample_submission: {sample_submission.shape}\")\n    print(f\"Columns: {list(sample_submission.columns)}\")\nexcept Exception as e:\n    print(f\"Error: {e}\")\n    sample_submission = pd.DataFrame({\n        'EntryID': ['Q8IZR5', 'Q8RZ23'],\n        'term': ['GO:0005515', 'GO:0140900'],\n        'probability': [0.5, 0.5]\n    })\n    print(f\"sample_submission: {sample_submission.shape}\")\n\n\nprint(\"\\nData loading complete.\")","metadata":{},"outputs":[],"execution_count":null},{"id":"fcc96aa2","cell_type":"code","source":"# Function to load FASTA sequences\ndef load_fasta_sequences(file_path):\n    \"\"\"Load sequences from FASTA file and return a dictionary {ID: sequence}.\"\"\"\n    sequences = {}\n    for record in SeqIO.parse(file_path, \"fasta\"):\n        if \"|\" in record.id:\n            prot_id = record.id.split(\"|\")[1] \n        else:\n            prot_id = record.id.split()[0]\n        \n        # Kiểm tra sơ bộ chuỗi không rỗng\n        if str(record.seq):\n            sequences[prot_id] = str(record.seq)\n            \n    return sequences\n\ntrain_sequences = load_fasta_sequences(TRAIN_PATH / \"train_sequences.fasta\")\nprint(f\"Loaded training sequences: {len(train_sequences):,} \")\n\ntest_sequences = load_fasta_sequences(TEST_PATH / \"testsuperset.fasta\")\nprint(f\"Loaded test sequences: {len(test_sequences):,} \")\n\ntrain_seq_df = pd.DataFrame([\n    {\"EntryID\": prot_id, \"sequence\": seq, \"seq_length\": len(seq)}\n    for prot_id, seq in train_sequences.items()\n])\n\ntest_seq_df = pd.DataFrame([\n    {\"EntryID\": prot_id, \"sequence\": seq, \"seq_length\": len(seq)}\n    for prot_id, seq in test_sequences.items()\n])\n\nprint(f\"\\nTrain sequences DF: {train_seq_df.shape}\")\nprint(f\"Test sequences DF: {test_seq_df.shape}\")","metadata":{},"outputs":[],"execution_count":null},{"id":"13c7b35e","cell_type":"code","source":"# Load GO ontology\ngo_graph = obonet.read_obo(str(TRAIN_PATH / \"go-basic.obo\"))\nprint(f\"Loaded GO graph: {len(go_graph)}\")\n\n# Extract root nodes for each ontology\nroot_nodes = {\n    'P': 'GO:0008150',  # Biological Process\n    'F': 'GO:0003674',  # Molecular Function\n    'C': 'GO:0005575'   # Cellular Component\n}\n\n# Get namespace for each GO term\ngo_namespaces = {}\nfor node_id in go_graph.nodes():\n    if 'namespace' in go_graph.nodes[node_id]:\n        namespace = go_graph.nodes[node_id]['namespace']\n        if 'biological_process' in namespace:\n            go_namespaces[node_id] = 'P'\n        elif 'molecular_function' in namespace:\n            go_namespaces[node_id] = 'F'\n        elif 'cellular_component' in namespace:\n            go_namespaces[node_id] = 'C'\n\nprint(f\"GO namespaces mapped: {len(go_namespaces)}\")\nprint(f\"Biological Process: {sum(1 for v in go_namespaces.values() if v == 'P')}\")\nprint(f\"Molecular Function: {sum(1 for v in go_namespaces.values() if v == 'F')}\")\nprint(f\"Cellular Component: {sum(1 for v in go_namespaces.values() if v == 'C')}\")","metadata":{},"outputs":[],"execution_count":null},{"id":"12ee701d","cell_type":"markdown","source":"### Data Cleaning & Pre-processing Sơ Bộ\n\nDựa trên phân tích về độ dài chuỗi và thành phần axit amin, chúng ta nhận thấy hai vấn đề chính cần xử lý trước khi tạo embedding: **protein siêu dài** và **các ký tự axit amin không chuẩn** (X, U, B, Z).","metadata":{}},{"id":"29ecaea5","cell_type":"code","source":"# Xử lý các ký tự Axit Amin lạ và Protein siêu dài\n\n# 1. Xử lý ký tự Axit Amin lạ (X, U, B, Z)\n# 'X' là axit amin không xác định (Unknown), chúng ta sẽ giữ nguyên hoặc gán về một token đặc biệt.\n# Tuy nhiên, trong mô hình ngôn ngữ protein hiện đại, giữ nguyên thường tốt hơn.\n# Nếu muốn thay thế:\naa_replacement_map = {'U': 'C', 'B': 'N', 'Z': 'Q'} # U->C (Selenocysteine), B/Z -> N/Q (ambiguous)\n\ndef clean_sequence(seq):\n    for old, new in aa_replacement_map.items():\n        seq = seq.replace(old, new)\n    return seq\n\n# Áp dụng làm sạch\ntrain_seq_df['sequence_clean'] = train_seq_df['sequence'].apply(clean_sequence)\ntest_seq_df['sequence_clean'] = test_seq_df['sequence'].apply(clean_sequence)\n\n# 2. Xử lý Protein siêu dài (potential outliers)\n# Giới hạn độ dài tối đa (max length) để chuẩn bị cho mô hình Transformer (ví dụ: 2048 hoặc 4096)\n# Ví dụ: Chỉ in ra số lượng protein vượt quá 5000 aa để xem xét có nên loại bỏ/cắt bớt không.\nMAX_SEQ_LEN_CUTOFF = 5000 \n\nlong_train = train_seq_df[train_seq_df['seq_length'] > MAX_SEQ_LEN_CUTOFF]\nlong_test = test_seq_df[test_seq_df['seq_length'] > MAX_SEQ_LEN_CUTOFF]\n\nprint(f\"Proteins > {MAX_SEQ_LEN_CUTOFF:,} aa (Train): {len(long_train):,}\")\nprint(f\"Proteins > {MAX_SEQ_LEN_CUTOFF:,} aa (Test): {len(long_test):,}\")","metadata":{},"outputs":[],"execution_count":null},{"id":"0699051c","cell_type":"markdown","source":"Ghi chú: Quyết định cắt (truncate) hay loại bỏ (drop) sẽ được thực hiện ở bước tiền xử lý mô hình.\nHiện tại, chỉ ghi nhận.","metadata":{}},{"id":"ad37f817","cell_type":"markdown","source":"<a id='3-dataset-structure'></a>\n## 3. Dataset Structure Analysis\n\n### 3.1 Data structure\n\nCấu trúc cơ bản của từng dataset:","metadata":{}},{"id":"6663298f","cell_type":"code","source":"print(\"_\"*80)\nprint(\"TRAIN TERMS\")\ndisplay(train_terms.head(10))\n\nprint(\"\\n\" + \"_\"*80)\nprint(\"TRAIN TAXONOMY\")\ndisplay(train_taxonomy.head(10))\n\nprint(\"\\n\" + \"_\"*80)\nprint(\"TRAIN SEQUENCES\")\ndisplay(train_seq_df.head(5))\n\nprint(\"\\n\" + \"_\"*80)\nprint(\"IA WEIGHTS\")\ndisplay(ia_weights.head(10))","metadata":{},"outputs":[],"execution_count":null},{"id":"792520c3","cell_type":"code","source":"print(\"DATASET STATISTICS\")\nprint(\"_\"*80)\n\n# Training data\nprint(\"\\nTRAINING DATA:\")\nprint(f\"- Total proteins: {train_terms['EntryID'].nunique():,}\")\nprint(f\"- Total GO term annotations: {len(train_terms):,}\")\nprint(f\"- Unique GO terms: {train_terms['term'].nunique():,}\")\nprint(f\"- Proteins with sequences: {len(train_sequences):,}\")\nprint(f\"- Proteins with taxonomy: {len(train_taxonomy):,}\")\n\n# Test data\nprint(\"\\nTEST DATA:\")\nprint(f\"- Total test proteins: {len(test_sequences):,}\")\n\n# Safe check for sample_submission\nif sample_submission is not None and 'EntryID' in sample_submission.columns:\n    print(f\"- Test proteins in submission: {sample_submission['EntryID'].nunique():,}\")\nelse:\n    print(f\"- Sample submission format: {sample_submission.shape if sample_submission is not None else 'N/A'}\")\n\n# GO terms\nprint(\"\\nGO ONTOLOGY:\")\nprint(f\"- Total GO terms in ontology: {len(go_graph):,}\")\nprint(f\"- GO terms with IA weights: {len(ia_weights):,}\")\n\n# Aspects distribution\naspect_counts = train_terms['aspect'].value_counts()\nprint(\"\\nASPECT DISTRIBUTION:\")\nfor aspect, count in aspect_counts.items():\n    aspect_name = {'P': 'Biological Process', 'F': 'Molecular Function', 'C': 'Cellular Component'}.get(aspect, aspect)\n    print(f\"- {aspect_name} ({aspect}): {count:,} ({count/len(train_terms)*100:.1f}%)\")","metadata":{},"outputs":[],"execution_count":null},{"id":"2205db86","cell_type":"code","source":"# Calculate annotations per protein\nannotations_per_protein = train_terms.groupby('EntryID').size()\n\nprint(\"ANNOTATIONS PER PROTEIN\")\nprint(\"_\"*80)\nprint(annotations_per_protein.describe())\n\n# Visualize\nfig, axes = plt.subplots(1, 2, figsize=(15, 5))\n\n# Histogram\naxes[0].hist(annotations_per_protein, bins=50, edgecolor='black', alpha=0.7, color='skyblue')\naxes[0].axvline(annotations_per_protein.mean(), color='red', linestyle='--', linewidth=2, label=f'Mean: {annotations_per_protein.mean():.1f}')\naxes[0].axvline(annotations_per_protein.median(), color='green', linestyle='--', linewidth=2, label=f'Median: {annotations_per_protein.median():.1f}')\naxes[0].set_xlabel('Number of GO Terms per Protein', fontsize=12)\naxes[0].set_ylabel('Frequency', fontsize=12)\naxes[0].set_title('Distribution of GO Terms per Protein', fontsize=14, fontweight='bold')\naxes[0].legend()\naxes[0].grid(axis='y', alpha=0.3)\n\n# Box plot\naxes[1].boxplot(annotations_per_protein, vert=True)\naxes[1].set_ylabel('Number of GO Terms', fontsize=12)\naxes[1].set_title('GO Terms per Protein - Box Plot', fontsize=14, fontweight='bold')\naxes[1].grid(axis='y', alpha=0.3)\n\nplt.tight_layout()\nplt.show()\n\nprint(f\"\\nInsight: Mỗi protein trung bình có {annotations_per_protein.mean():.1f} GO terms\")\nprint(f\"Insight: {(annotations_per_protein > 100).sum():,} proteins có hơn 100 GO terms\")","metadata":{},"outputs":[],"execution_count":null},{"id":"d575687c","cell_type":"markdown","source":"### Nhận xét dựa trên thống kê ANNOTATIONS PER PROTEIN\n\n1. **Dữ liệu có tính đa nhãn (multi-label)**\n   Một protein điển hình mang khoảng 4 GO term (median = 4), trung bình có 6–7 GO (mean = 6.5). Hầu hết protein có **nhiều hơn 1 nhãn**.\n\n2. **Độ phức tạp chức năng giữa các protein không đồng đều**  \n   - Phần lớn protein nằm trong khoảng **1–10 GO**.  \n   - Chỉ có 26 proteins mang hơn 100 GO terms .  \n   --> Tồn tại một nhóm nhỏ gồm các protein được annotate cực kỳ dày đặc.\n\n3. **Annotation bias theo mức độ nghiên cứu**  \n   - Các protein có 100+ GO nhiều khả năng là các protein model (human, yeast, E. coli, …) được nghiên cứu kỹ.  \n   - Phần lớn protein chỉ có vài GO ít được annotate.  \n   --> Dữ liệu không đồng đều cả về sinh học lẫn mức độ curating.","metadata":{}},{"id":"055e81bb","cell_type":"markdown","source":"### Ý nghĩa thống kê:\n\n1. **Loss dễ bị chi phối bởi protein có nhiều nhãn**  \n   - Protein với 100–200 GO sẽ tạo ra nhiều positive labels hơn protein chỉ có 3–4 GO.  \n     -> Chuẩn hoá loss theo số nhãn mỗi protein hoặc áp dụng weighting để giảm ảnh hưởng của các protein được annotate quá dày.\n\n2. **Cần chú ý đến cách tạo batch để gradient ổn định**  \n   - Phân phối lệch phải khiến một batch ngẫu nhiên có thể chứa vài protein có >100 GO → gradient không đều  \n     -> Thiết kế batch trộn đều protein ít nhãn và nhiều nhãn (stratified batching theo GO).\n\n3. **Threshold không nên cố định**  \n   - Mật độ nhãn khác nhau giữa protein (và giữa các aspect P/C/F).  \n   - Cố định threshold dẫn đến:\n     - Over-predict cho protein ít GO,  \n     - Under-predict cho protein nhiều GO.  \n   --> Cần tune threshold trên validation, tune theo từng ontology (P/C/F), để tối ưu F-max / IA-Fmax.\n\n4. **Thiết kế tập validation cần phản ánh đúng phân bố độ phức tạp**  \n   - Nếu chia train/val ngẫu nhiên, tỷ lệ protein có thể khác nhau giữa hai tập.  \n   --> stratified split theo số GO/protein để metric trên validation phản ánh đúng khả năng tổng quát hoá của mô hình.\n","metadata":{}},{"id":"86d6f87f","cell_type":"markdown","source":"<a id='4-go-term-analysis'></a>\n## 4. GO Term Analysis\n\n### 4.1 GO Term Frequency Distribution\n\nPhân tích tần suất xuất hiện của các GO terms trong training data:","metadata":{}},{"id":"37416043","cell_type":"code","source":"go_term_counts = train_terms['term'].value_counts()\n\nprint(\"GO TERM FREQUENCY STATISTICS\")\nprint(\"_\"*80)\nprint(go_term_counts.describe())\n\n# Top 20 most frequent GO terms\nprint(\"\\nTOP 20 MOST FREQUENT GO TERMS:\")\nprint(\"_\"*80)\n\ntop_20_terms = go_term_counts.head(20)\nfor idx, (term, count) in enumerate(top_20_terms.items(), 1):\n    term_name = go_graph.nodes[term].get('name', 'Unknown') if term in go_graph.nodes else 'Unknown'\n    aspect = go_namespaces.get(term, '?')\n    print(f\"{idx:2d}. {term} ({aspect}) - {count:,} proteins\")\n    print(f\"--> {term_name}\")\n\n# Visualize\nfig, axes = plt.subplots(2, 1, figsize=(16, 10))\n\n# Top 20 bar chart\naxes[0].barh(range(20), top_20_terms.values[::-1], color='coral', edgecolor='black')\naxes[0].set_yticks(range(20))\naxes[0].set_yticklabels([f\"{term}\\n({go_namespaces.get(term, '?')})\" for term in top_20_terms.index[::-1]], fontsize=9)\naxes[0].set_xlabel('Number of Proteins', fontsize=12)\naxes[0].set_title('Top 20 Most Frequent GO Terms', fontsize=14, fontweight='bold')\naxes[0].grid(axis='x', alpha=0.3)\n\n# Distribution of all GO term frequencies\naxes[1].hist(go_term_counts, bins=100, edgecolor='black', alpha=0.7, color='teal')\naxes[1].axvline(go_term_counts.mean(), color='red', linestyle='--', linewidth=2, label=f'Mean: {go_term_counts.mean():.1f}')\naxes[1].axvline(go_term_counts.median(), color='green', linestyle='--', linewidth=2, label=f'Median: {go_term_counts.median():.1f}')\naxes[1].set_xlabel('Number of Proteins with GO Term', fontsize=12)\naxes[1].set_ylabel('Number of GO Terms (Log Scale)', fontsize=12) # Cập nhật label\naxes[1].set_title('Distribution of GO Term Frequencies (Log Y)', fontsize=14, fontweight='bold')\naxes[1].set_yscale('log') # Đảm bảo Y-axis là Log Scale\naxes[1].legend()\naxes[1].grid(axis='y', alpha=0.3)\n\nplt.tight_layout()\nplt.show()\n\n# Long-tail analysis\nprint(f\"\\nClass Imbalance Analysis:\")\nprint(f\"- GO terms appearing in >1000 proteins: {(go_term_counts > 1000).sum():,}\")\nprint(f\"- GO terms appearing in >100 proteins: {(go_term_counts > 100).sum():,}\")\nprint(f\"- GO terms appearing in >10 proteins: {(go_term_counts > 10).sum():,}\")\nprint(f\"- GO terms appearing in only 1 protein: {(go_term_counts == 1).sum():,}\")","metadata":{},"outputs":[],"execution_count":null},{"id":"304d2edb","cell_type":"markdown","source":"### Phân tích dựa trên thống kê GO TERM FREQUENCY STATISTICS\n1. **Phần lớn GO term là rất hiếm**  \n   - 50% GO term gắn cho **tối đa 4 protein**, 75% GO term gắn cho **tối đa 12 protein**.  \n   -> Đa số GO nằm ở vùng tần suất rất thấp.\n\n2. **Top 20 GO term cực kỳ phổ biến**  \n   - Một số GO xuất hiện ở từ vài nghìn đến hơn ba chục nghìn protein.  \n   - Các GO này đều là các khái niệm rất chung như protein binding, nucleus, cytosol, plasma membrane, cytoplasm.  \n   - Đây thường là các term nằm cao trong DAG, mức độ khái quát lớn.\n\n3. **Phân phối tần suất GO term lệch phải rất mạnh**  \n   - Median 4 nhưng max 33,713 cho thấy sự chênh lệch cực lớn giữa các nhãn.  \n   - Dữ liệu tồn tại một nhóm rất nhỏ GO head (tần suất rất cao) và một số lượng lớn GO tail (tần suất rất thấp).","metadata":{}},{"id":"b0edca63","cell_type":"markdown","source":"### Ý nghĩa \n\n1. **Mô hình dễ thiên lệch về các GO phổ biến**  \n   - Các GO như protein binding, nucleus, cytosol xuất hiện ở hàng chục nghìn protein, chúng đóng góp rất lớn vào tổng loss.  \n   - Các GO hiếm chỉ có vài mẫu không đủ dữ liệu để mô hình học riêng biệt.  \n   --> Nếu chỉ dùng loss chuẩn (BCE) không weight, mô hình sẽ ưu tiên tối ưu cho các head labels và dự đoán rất kém các GO tail.\n\n2. **Metric micro-F1 chủ yếu phản ánh chất lượng trên head labels**  \n   - Mỗi lần dự đoán đúng hoặc sai một GO phổ biến ảnh hưởng đến hàng nghìn cặp protein–GO.  \n   - Ngược lại, một GO rất hiếm dù dự đoán tốt hay kém cũng không thay đổi đáng kể micro-F1.  \n   --> Micro-F1 và IA-Fmax bị chi phối mạnh bởi các GO phổ biến, trong khi chất lượng trên GO hiếm khó nhìn thấy nếu không tách phân tích.\n\n3. **Cấu trúc phân cấp của GO hỗ trợ các nhãn tail**  \n   - Các GO hiếm thường là con của các GO phổ biến.  \n   - Mô hình có thể học tốt các ancestor head, sau đó dùng:\n     - Propagation theo DAG (true-path-rule)  \n     - Hoặc ràng buộc phân cấp trong loss để truyền thông tin từ head xuống tail.  \n   --> GO propagation và các kỹ thuật hierarchical modeling.\n\n4. **Kỹ thuật xử lý imbalance ở mức label**  \n   Một số hướng phù hợp với phân phối hiện tại:\n   - Dùng **class weighting** hoặc **focal loss** để tăng trọng số cho GO hiếm.  \n   - Huấn luyện theo từng tầng khái quát → cụ thể (tiering): học tốt các GO phổ biến trước, sau đó fine tune cho các GO chi tiết hơn.  \n   - Bổ sung regularization hoặc embedding cho chính các GO term có cấu trúc DAG.\n\nNhìn chung, phân tích tần suất GO term cho thấy bài toán CAFA6 không chỉ bị imbalance ở mức protein, mà còn cực kỳ mất cân bằng ở mức nhãn. Mọi chiến lược huấn luyện cần tính đến việc vừa tận dụng tốt các GO head, vừa không bỏ qua hoàn toàn các GO tail mang tính đặc hiệu sinh học cao.\n","metadata":{}},{"id":"bb4e5294","cell_type":"code","source":"# Analyze GO terms by aspect\nvalid_aspects = ['P', 'F', 'C']\ntrain_terms_valid = train_terms[train_terms['aspect'].isin(valid_aspects)].copy()\n\naspect_term_counts = train_terms_valid.groupby('aspect')['term'].nunique()\n\nprint(\"UNIQUE GO TERMS BY ASPECT\")\nprint(\"_\"*80)\nfor aspect, count in aspect_term_counts.items():\n    aspect_name = {'P': 'Biological Process', 'F': 'Molecular Function', 'C': 'Cellular Component'}.get(aspect, aspect)\n    print(f\"- {aspect_name} ({aspect}): {count:,} unique terms\")\n\n# Visualize\nfig = plt.figure(figsize=(16, 6))\n\n# Pie chart\nax1 = plt.subplot(1, 2, 1)\ncolors = ['#ff9999', '#66b3ff', '#99ff99']\nn_aspects = len(aspect_term_counts)\nexplode = tuple([0.05] * n_aspects)  # Dynamic explode based on number of aspects\n\nwedges, texts, autotexts = ax1.pie(\n    aspect_term_counts.values,\n    labels=[{'P': 'Biological\\nProcess', 'F': 'Molecular\\nFunction', 'C': 'Cellular\\nComponent'}.get(a, a) for a in aspect_term_counts.index],\n    autopct='%1.1f%%',\n    colors=colors[:n_aspects],\n    explode=explode,\n    shadow=True,\n    startangle=90\n)\nfor autotext in autotexts:\n    autotext.set_color('white')\n    autotext.set_fontsize(12)\n    autotext.set_fontweight('bold')\nax1.set_title('GO Term Distribution by Aspect', fontsize=14, fontweight='bold')\n\n# Bar chart - annotations by aspect\nax2 = plt.subplot(1, 2, 2)\naspect_annotation_counts = train_terms_valid['aspect'].value_counts()\nbars = ax2.bar(\n    [{'P': 'Biological\\nProcess', 'F': 'Molecular\\nFunction', 'C': 'Cellular\\nComponent'}.get(a, a) for a in aspect_annotation_counts.index],\n    aspect_annotation_counts.values,\n    color=colors[:len(aspect_annotation_counts)],\n    edgecolor='black',\n    alpha=0.7\n)\nax2.set_ylabel('Number of Annotations', fontsize=12)\nax2.set_title('Total Annotations by Aspect', fontsize=14, fontweight='bold')\nax2.grid(axis='y', alpha=0.3)\n\n# Add value labels on bars\nfor bar in bars:\n    height = bar.get_height()\n    ax2.text(bar.get_x() + bar.get_width()/2., height,\n             f'{int(height):,}',\n             ha='center', va='bottom', fontsize=11, fontweight='bold')\n\nplt.tight_layout()\nplt.show()","metadata":{},"outputs":[],"execution_count":null},{"id":"f47ae76f","cell_type":"markdown","source":"### Nhận xét dựa trên UNIQUE GO TERMS BY ASPECT\n\n1. **Biological Process chiếm ưu thế rõ rệt về độ phong phú của nhãn**  \n   - Khoảng 64–65 % số GO term duy nhất thuộc P. \n   - Số lượng annotations của P cao nhất.\n   --> Ontology P chiếm nhiều GO terms nhất và được gán nhãn nhiều nhất\n\n2. **Cellular Component là ontology bé nhất nhưng mỗi term tương đối phổ biến**  \n   - Chỉ khoảng 10 % GO term duy nhất là C (2.651 term), nhưng annotations chiếm khá cao so với số lượng GO terms mà nó có.  \n   --> So với số lượng nhãn, các term C thường được gán cho nhiều protein hơn, biểu hiện các tế bào phổ biến.\n\n3. **Sự mất cân bằng giữa ba aspect**  \n   - Vừa có mất cân bằng về số lượng nhãn (P >> F > C), vừa có mất cân bằng về số annotation mỗi aspect.  \n   --> Bộ dữ liệu không cân bằng theo P, C, F mà nghiêng mạnh về Biological Process.","metadata":{}},{"id":"d8c898b0","cell_type":"markdown","source":"### Ý nghĩa:\n\nNhững đặc điểm trên dẫn đến một số điểm quan trọng cho thiết kế và đánh giá mô hình:\n\n1. **Cần xử lý riêng từng aspect thay vì gộp tất cả**  \n   - Biological Process có số nhãn và số annotation lớn nhất, dễ xuất hiện long tail mạnh và nhiều nhãn rất hiếm.  \n   - Cellular Component ít nhãn hơn nhưng mỗi nhãn phổ biến hơn, có thể dễ học hơn tương đối.  \n   --> huấn luyện và tối ưu ngưỡng prediction nên được thực hiện theo từng aspect P, C, F để phù hợp với phân bố riêng của từng nhánh.\n\n2. **Kỳ vọng hiệu năng giữa các aspect sẽ khác nhau**  \n   - Với P: số nhãn lớn, phân phối phức tạp -> mô hình có thể khó đạt Fmax cao, dễ overfit hoặc underfit tail labels.  \n   - Với C: ít nhãn, nhiều annotation/nhãn -> mô hình có nhiều dữ liệu hơn cho mỗi term, thường đạt điểm cao hơn.  \n   - Với F: nằm giữa P và C, vừa có nhiều nhãn, vừa không quá giàu dữ liệu.  \n   --> Trình bày riêng Fmax/IA-Fmax cho từng P, C, F.\n\n3. **Chiến lược capacity và regularization theo aspect**  \n   - Có thể cần mô hình hoặc head lớn hơn cho P (nhiều nhãn, cấu trúc phức tạp) và regularization mạnh để tránh overfitting.  \n   - Với C, mô hình đơn giản hơn hoặc ngưỡng cao hơn có thể đủ tốt vì nhãn phổ biến hơn.\n\n4. **Tuning threshold và ensemble theo aspect**  \n   - Tỷ lệ positive và độ khó khác nhau theo P, C, F nên một ngưỡng score chung cho tất cả aspect thường không tối ưu.  \n   --> Chọn threshold và trọng số ensemble riêng cho từng aspect sẽ tận dụng tốt hơn cấu trúc và mức độ phong phú nhãn ở mỗi nhánh ontology.\n\nNhìn chung, phân bố GO theo aspect cho thấy Biological Process là nhánh chi phối dữ liệu, trong khi Cellular Component và Molecular Function có động học nhãn khác nhau. Các quyết định về kiến trúc mô hình, regularization và tuning threshold nên được thiết kế có ý thức theo từng aspect P, C, F.\n","metadata":{}},{"id":"4ea36098","cell_type":"markdown","source":"### Phân tích sâu hơn mối quan hệ giữa Tần suất (Frequency) và Trọng số Thông tin (IA Weight)\n - Dữ liệu đã cho thấy IA có tương quan nghịch với Frequency.\n - Dưới đây là trực quan hóa mối quan hệ này, phân biệt theo 3 khía cạnh (Aspect) P, F, C.\nMục tiêu là xác định Aspect nào có sự phân tán IA lớn nhất ở vùng tần suất thấp (GO terms tail).","metadata":{}},{"id":"dc41d3ae","cell_type":"code","source":"print(\"INFORMATION ACCRETION (IA) WEIGHTS ANALYSIS\")\nprint(\"_\"*80)\n\n# Basic statistics\nprint(\"\\nIA Weight Statistics:\")\nprint(ia_weights['IA'].describe())\n\n# Add aspect information\nia_weights['aspect'] = ia_weights['term'].map(go_namespaces)\n\n# Visualize\nfig, axes = plt.subplots(2, 2, figsize=(16, 12))\n\n# Distribution of IA weights\naxes[0, 0].hist(ia_weights['IA'], bins=50, edgecolor='black', alpha=0.7, color='purple')\naxes[0, 0].axvline(ia_weights['IA'].mean(), color='red', linestyle='--', linewidth=2, label=f'Mean: {ia_weights[\"IA\"].mean():.3f}')\naxes[0, 0].axvline(ia_weights['IA'].median(), color='green', linestyle='--', linewidth=2, label=f'Median: {ia_weights[\"IA\"].median():.3f}')\naxes[0, 0].set_xlabel('IA Weight', fontsize=12)\naxes[0, 0].set_ylabel('Frequency', fontsize=12)\naxes[0, 0].set_title('Distribution of IA Weights', fontsize=14, fontweight='bold')\naxes[0, 0].legend()\naxes[0, 0].grid(axis='y', alpha=0.3)\n\n# Box plot by aspect\nia_weights_by_aspect = [ia_weights[ia_weights['aspect'] == asp]['IA'].dropna() for asp in ['P', 'F', 'C']]\nbp = axes[0, 1].boxplot(ia_weights_by_aspect, labels=['Biological\\nProcess', 'Molecular\\nFunction', 'Cellular\\nComponent'])\naxes[0, 1].set_ylabel('IA Weight', fontsize=12)\naxes[0, 1].set_title('IA Weights by Aspect', fontsize=14, fontweight='bold')\naxes[0, 1].grid(axis='y', alpha=0.3)\n\n# Top 20 highest IA weights\ntop_ia = ia_weights.nlargest(20, 'IA')\naxes[1, 0].barh(range(20), top_ia['IA'].values[::-1], color='orange', edgecolor='black')\naxes[1, 0].set_yticks(range(20))\naxes[1, 0].set_yticklabels([f\"{term}\\n({asp})\" for term, asp in zip(top_ia['term'].values[::-1], top_ia['aspect'].values[::-1])], fontsize=8)\naxes[1, 0].set_xlabel('IA Weight', fontsize=12)\naxes[1, 0].set_title('Top 20 GO Terms by IA Weight', fontsize=14, fontweight='bold')\naxes[1, 0].grid(axis='x', alpha=0.3)\n\n# Scatter: GO term frequency vs IA weight\n# Merge frequency with IA\ngo_freq_ia = pd.DataFrame({\n    'term': go_term_counts.index,\n    'frequency': go_term_counts.values\n}).merge(ia_weights[['term', 'IA']], on='term', how='left')\n\naxes[1, 1].scatter(go_freq_ia['frequency'], go_freq_ia['IA'], alpha=0.5, s=10, color='teal')\naxes[1, 1].set_xlabel('GO Term Frequency (# proteins)', fontsize=12)\naxes[1, 1].set_ylabel('IA Weight', fontsize=12)\naxes[1, 1].set_title('GO Term Frequency vs IA Weight', fontsize=14, fontweight='bold')\naxes[1, 1].set_xscale('log')\naxes[1, 1].grid(alpha=0.3)\n\nplt.tight_layout()\nplt.show()\n\nprint(f\"\\nInsight: IA weights range from {ia_weights['IA'].min():.4f} to {ia_weights['IA'].max():.4f}\")\nprint(f\"Insight: Higher IA weights indicate more specific/informative GO terms\")","metadata":{},"outputs":[],"execution_count":null},{"id":"a3530a29","cell_type":"markdown","source":"### Nhận xét dựa trên INFORMATION ACCRETION (IA) WEIGHTS ANALYSIS\n\n1. **Phân phối IA weight lệch phải rất mạnh**\n\n   - Phần lớn GO term có IA gần **0** (25% số term có IA=0, median chỉ đạt 0.15).  \n   - Một số term có IA rất cao, kéo đuôi phân phối tới khoảng 14-15.\n--> \n\n2. **IA theo aspect: MF và CC thường giàu thông tin hơn P**\n\n   - Boxplot theo aspect cho thấy:\n     - Biological Process có median IA thấp nhất, nhiều term có IA gần 0.  \n     - Molecular Function và Cellular Component có median IA cao hơn và khoảng tứ phân vị rộng hơn.\n   - Điều này phù hợp với trực giác sinh học:  \n     - Nhiều GO Process khá tổng quát (ví dụ các con đường điều hòa, quá trình trao đổi chất).  \n     - GO Function và Component thường mô tả chức năng phân tử hoặc vị trí tế bào cụ thể hơn.\n\n3. **Các GO có IA cao nhất là các term cực kỳ đặc hiệu**\n\n   - Top 20 GO có IA đều nằm quanh **14–15**.  \n   - Đều là các term hiếm có vị trí con trong DAG, mô tả chức năng hay vị trí rất cụ thể.\n\n4. **Quan hệ giữa tần suất GO và IA: tương quan nghịch rõ rệt**\n\n   - Biểu đồ scatter IA vs tần suất (log) cho thấy:\n     - GO xuất hiện ở nhiều protein (tần suất cao) gần như luôn có IA thấp, sát 0.  \n     - GO ít xuất hiện (tần suất thấp) có thể đạt IA rất cao.  \n   - Điều này khẳng định IA đang đo mức độ đặc hiệu: term càng hiếm càng mang nhiều thông tin.","metadata":{}},{"id":"8c118592","cell_type":"markdown","source":"### Ý nghĩa\n\nCác đặc điểm của IA weight dẫn đến một số điểm quan trọng khi thiết kế mô hình và tối ưu theo IA-Fmax:\n\n1. **IA-Fmax tập trung vào các GO đặc hiệu, không phải chỉ những GO phổ biến**\n\n   - Những GO có IA cao (term hiếm, cụ thể) đóng góp rất nhiều vào IA-Fmax.  \n   - Nếu mô hình chỉ đoán được các term phổ biến với IA thấp, điểm micro-F1 có thể khá tốt nhưng IA-Fmax vẫn thấp.  \n   → Khi mục tiêu đánh giá dùng IA-Fmax, mô hình phải học được cả các GO tail ít dữ liệu chứ không chỉ head labels.\n\n2. **Cần cân nhắc sử dụng IA như trọng số trong huấn luyện**\n\n   - Vì IA đã phản ánh mức độ quan trọng thông tin của từng term, có thể:\n     - Áp dụng **loss weighting theo IA** để ưu tiên các GO giàu thông tin.  \n     - Hoặc sử dụng IA trong sampling (ví dụ tăng xác suất xuất hiện của các GO có IA cao trong minibatch).\n   - Mục tiêu là đưa phân bố trọng số trong loss lại gần hơn với phân bố trọng số mà IA-Fmax sử dụng trong đánh giá.\n\n3. **Khác biệt IA giữa các aspect gợi ý huấn luyện tách biệt theo P/C/F**\n\n   - Molecular Function và Cellular Component có median IA cao hơn: sai sót ở hai nhánh này làm giảm IA-Fmax mạnh hơn.  \n   - Biological Process có nhiều term IA thấp; nếu tối ưu đồng đều, mô hình có thể tiêu tốn sức vào các term ít giá trị IA.  \n   → Có thể:\n     - Huấn luyện hoặc tune threshold riêng cho từng aspect.  \n     - Tăng trọng số loss cho các term MF và CC có IA cao nếu muốn tối ưu IA-Fmax toàn cục.\n\n4. **GO propagation cần được thiết kế phù hợp với IA**\n\n   - Propagation chủ yếu thêm các ancestor có IA thấp hơn con.  \n   - Nếu mô hình chỉ dự đoán đúng ancestor nhưng không bắt được leaf hoặc term trung gian có IA cao, IA-Fmax vẫn không cao.  \n   → Trong thiết kế mô hình và hậu xử lý, cần hướng tới việc dự đoán đúng các GO cụ thể (IA lớn), propagation chỉ là bước bảo đảm đúng cấu trúc DAG chứ không thay thế được thông tin từ các term giàu IA.\n\nTóm lại, phân tích IA cho thấy bài toán CAFA không chỉ là dự đoán nhiều nhãn đúng, mà là dự đoán đúng các nhãn mang nhiều thông tin. Toàn bộ pipeline từ loss, sampling, tách aspect đến propagation và ensemble nên được thiết kế có ý thức theo trọng số IA nếu mục tiêu cuối cùng là tối ưu IA-Fmax.\n","metadata":{}},{"id":"e704a35a","cell_type":"markdown","source":"## ⚖️ Tóm Tắt Ý Nghĩa Trọng Số IA\n\nViệc phân tích cho thấy mối tương quan **nghịch rõ rệt** giữa tần suất xuất hiện của GO Term và trọng số IA của nó.\n\n* **GO phổ biến (Head):** Có tần suất cao nhưng IA thấp (gần 0). Dự đoán đúng các GO này giúp cải thiện Micro-F1 nhưng **ảnh hưởng ít** đến IA-Fmax.\n* **GO đặc hiệu (Tail):** Tần suất hiếm nhưng có IA cao (lên tới 14.86). Dự đoán đúng một GO term như vậy mang lại **giá trị rất lớn** cho điểm IA-Fmax.\n\nĐiều này xác nhận rằng, để tối ưu hóa IA-Fmax, mô hình không thể chỉ tập trung vào các nhãn phổ biến, mà phải tìm cách học được các **GO terms hiếm và đặc hiệu**.\n\n### Chiến lược Mô hình Hóa Tối ưu IA-Fmax:\n\n1.  **Loss Function:** Áp dụng **Loss Weighting** theo IA để tăng cường tín hiệu gradient cho các GO Term giàu thông tin (IA cao).\n2.  **Aspect Tuning:** Ưu tiên học và tinh chỉnh ngưỡng cho Molecular Function (F) và Cellular Component (C), vì chúng có median IA cao hơn Biological Process (P).\n3.  **Hierarchical:** Sử dụng các kỹ thuật ràng buộc phân cấp (Hierarchical Modeling) và **Propagation** để giúp dự đoán các nhãn con (term IA cao) dựa trên nhãn cha (term IA thấp).","metadata":{}},{"id":"f8b54dec","cell_type":"markdown","source":"<a id='5-sequence-analysis'></a>\n## 5. Protein Sequence Analysis\n\n### 5.1 Sequence Length Distribution\n\nPhân tích độ dài của các chuỗi protein:","metadata":{}},{"id":"4195e9a3","cell_type":"code","source":"print(\"PROTEIN SEQUENCE LENGTH ANALYSIS\")\nprint(\"_\"*80)\n\n# Statistics\nprint(\"\\nTRAINING SEQUENCES:\")\nprint(train_seq_df['seq_length'].describe())\n\nprint(\"\\nTEST SEQUENCES:\")\nprint(test_seq_df['seq_length'].describe())\n\n# Visualize\nfig, axes = plt.subplots(2, 2, figsize=(16, 12))\n\n# Train sequence length distribution\naxes[0, 0].hist(train_seq_df['seq_length'], bins=100, edgecolor='black', alpha=0.7, color='steelblue')\naxes[0, 0].axvline(train_seq_df['seq_length'].mean(), color='red', linestyle='--', linewidth=2, label=f'Mean: {train_seq_df[\"seq_length\"].mean():.0f}')\naxes[0, 0].axvline(train_seq_df['seq_length'].median(), color='green', linestyle='--', linewidth=2, label=f'Median: {train_seq_df[\"seq_length\"].median():.0f}')\naxes[0, 0].set_xlabel('Sequence Length (amino acids)', fontsize=12)\naxes[0, 0].set_ylabel('Frequency', fontsize=12)\naxes[0, 0].set_title('Training Sequence Length Distribution', fontsize=14, fontweight='bold')\naxes[0, 0].legend()\naxes[0, 0].grid(axis='y', alpha=0.3)\n\n# Test sequence length distribution\naxes[0, 1].hist(test_seq_df['seq_length'], bins=100, edgecolor='black', alpha=0.7, color='coral')\naxes[0, 1].axvline(test_seq_df['seq_length'].mean(), color='red', linestyle='--', linewidth=2, label=f'Mean: {test_seq_df[\"seq_length\"].mean():.0f}')\naxes[0, 1].axvline(test_seq_df['seq_length'].median(), color='green', linestyle='--', linewidth=2, label=f'Median: {test_seq_df[\"seq_length\"].median():.0f}')\naxes[0, 1].set_xlabel('Sequence Length (amino acids)', fontsize=12)\naxes[0, 1].set_ylabel('Frequency', fontsize=12)\naxes[0, 1].set_title('Test Sequence Length Distribution', fontsize=14, fontweight='bold')\naxes[0, 1].legend()\naxes[0, 1].grid(axis='y', alpha=0.3)\n\n# Box plot comparison\naxes[1, 0].boxplot([train_seq_df['seq_length'], test_seq_df['seq_length']], labels=['Train', 'Test'])\naxes[1, 0].set_ylabel('Sequence Length', fontsize=12)\naxes[1, 0].set_title('Train vs Test Sequence Length', fontsize=14, fontweight='bold')\naxes[1, 0].grid(axis='y', alpha=0.3)\n\n# Log-scale histogram\naxes[1, 1].hist(train_seq_df['seq_length'], bins=100, edgecolor='black', alpha=0.5, color='steelblue', label='Train')\naxes[1, 1].hist(test_seq_df['seq_length'], bins=100, edgecolor='black', alpha=0.5, color='coral', label='Test')\naxes[1, 1].set_xlabel('Sequence Length (amino acids)', fontsize=12)\naxes[1, 1].set_ylabel('Frequency (log scale)', fontsize=12)\naxes[1, 1].set_title('Train vs Test - Log Scale', fontsize=14, fontweight='bold')\naxes[1, 1].set_yscale('log')\naxes[1, 1].legend()\naxes[1, 1].grid(axis='y', alpha=0.3)\n\nplt.tight_layout()\nplt.show()\n\n# Length categories\nprint(\"\\nSEQUENCE LENGTH CATEGORIES:\")\nbins = [0, 100, 200, 500, 1000, 2000, float('inf')]\nlabels = ['<100', '100-200', '200-500', '500-1000', '1000-2000', '>2000']\n\ntrain_cats = pd.cut(train_seq_df['seq_length'], bins=bins, labels=labels)\ntest_cats = pd.cut(test_seq_df['seq_length'], bins=bins, labels=labels)\n\nprint(\"\\nTrain:\")\nprint(train_cats.value_counts().sort_index())\nprint(\"\\nTest:\")\nprint(test_cats.value_counts().sort_index())","metadata":{},"outputs":[],"execution_count":null},{"id":"ffe86a46","cell_type":"markdown","source":"### Nhận xét PROTEIN SEQUENCE LENGTH ANALYSIS:\n\n1. **Độ dài chuỗi có phân phối lệch phải mạnh với đuôi rất dài**  \n- Dữ liệu cho thấy một đặc điểm rất điển hình của protein tự nhiên: phần lớn chuỗi tập trung ở vùng ngắn-trung bình (200-1.000 amino acid), nhưng có một nhóm nhỏ protein cực dài kéo dài đến hơn 35.000 amino acid. Điều này được thể hiện rõ qua biểu đồ boxplot với rất nhiều outlier phía trên, cho thấy sự đa dạng lớn về kích thước protein trong tự nhiên.\n- Phân bố theo bin:\n   + Train: Median ~409 aa, mean ~526 aa\n   + Test: Median ~342 aa, mean ~429 aa\n   + Cả hai tập đều tập trung mạnh ở khoảng 100-500 aa với tần suất giảm dần theo độ dài\n\n2. **Train và test có phân bố khá tương đồng nhưng test hơi ngắn hơn**  \n   - Mean và median của test đều **thấp hơn** train (429 vs 526; 342 vs 409).  \n   - Tỉ lệ chuỗi rất ngắn (<100 aa) **cao hơn** trong test so với train, trong khi các nhóm dài hơn có tỉ lệ tương tự.  \n   Nhìn chung phân bố của hai tập khá tương đồng, nhưng có thể nhận thấy test set nghiêng về phía các protein ngắn hơn. Test có median thấp hơn khoảng 67 aa và mean thấp hơn gần 100 aa so với train. Tỉ lệ protein rất ngắn (<100 aa) trong test cũng cao hơn đáng kể, với 25.125 chuỗi so với 3.517 chuỗi trong train.\n   Tuy nhiên, sự khác biệt này không quá nghiêm trọng khi nhìn vào histogram log-scale - hai phân bố vẫn chồng lên nhau khá tốt, chỉ lệch nhẹ về vị trí trung tâm.\n\n3. **Phần lớn chuỗi nằm trong khoảng 200–1.000 aa ở cả train và test**  \n   Hai nhóm chính chiếm đa số dữ liệu:\n      + 200-500 aa: Khoảng 37.000+ mẫu (train) và 98.000+ mẫu (test)\n      + 500-1.000 aa: Khoảng 22.000+ mẫu (train) và 49.000+ mẫu (test)  \n   Các chuỗi >2.000 aa là thiểu số nhưng vẫn có vài nghìn mẫu. Những chuỗi này thường là các protein đa miền phức tạp, có nhiều vùng chức năng khác nhau hoặc các chuỗi polypeptide rất dài.\n","metadata":{}},{"id":"9b40f1ef","cell_type":"markdown","source":"### Ý nghĩa:\n\n1. **Thiết kế kiến trúc phù hợp với dải độ dài rộng**\n   - Sự xuất hiện của các protein >10.000 aa, thậm chí 35.000 aa, vượt xa độ dài tối đa của nhiều mô hình chuỗi thông dụng.  \n   - Một số chiến lược có thể cân nhắc:\n      + Với mô hình sequence-based (RNN, LSTM, 1D-CNN): Cần có cơ chế cắt/pooling thông minh hoặc xử lý theo sliding window. Protein quá dài có thể được chia thành các đoạn nhỏ hơn, sau đó aggregate lại.\n      + Với Transformer: Hạn chế về độ dài context window (thường ~1024-2048 tokens) là vấn đề lớn với protein dài. \n      + Với embedding có sẵn (ESM-2, ProtBERT): Đây là lựa chọn tốt vì vấn đề xử lý độ dài đã được giải quyết ở bước tạo embedding. Cần chú ý cách embedding xử lý protein siêu dài.\n\n2. **Quản lý memory và batch size hiệu quả**\n\n   - Các batch chứa nhiều chuỗi dài sẽ tốn RAM/GPU hơn rất nhiều so với batch chứa chuỗi ngắn.  \n   - Cách xử lý:\n      + Bucketing by length: Nhóm các protein có độ dài gần nhau vào cùng batch, giúp tránh padding quá nhiều và tận dụng memory tốt hơn.\n      + Dynamic batching: Điều chỉnh batch size theo độ dài trung bình trong batch đó.\n\n3. **Kiểm tra generalization theo bin độ dài**\n\n   Vì test có xu hướng hơi ngắn hơn train và phân phối có đuôi dài nên cần:\n     + Đánh giá performance theo các bin độ dài ( <100, 100–200, 200–500, 500–1.000, >1.000 aa).  \n     + Kiểm tra xem mô hình có suy giảm chất lượng mạnh ở các chuỗi rất ngắn hoặc rất dài hay không.  \n   --> Tránh tình huống mô hình chỉ hoạt động tốt ở vùng độ dài trung tâm (200–500 aa) nhưng yếu ở tail.\n\n4. **Lưu ý về protein siêu dài (>10.000 aa)**\n\n   Những protein này chiếm rất ít (<1%) nhưng lại có thể gây ra nhiều vấn đề:\n      + Tốn rất nhiều memory\n      + Có thể là noise hoặc lỗi annotation (ví dụ chuỗi 35.000 aa rất bất thường)\n      + Nếu dùng sliding window, một protein này có thể tạo ra hàng chục samples\n\n**Tóm lại**\nTừ phân tích độ dài chuỗi có thể thấy bộ dữ liệu có độ phức tạp đáng kể với từ vài chục đến hàng chục nghìn amino acid. Mô hình cần được thiết kế để xử lý hiệu quả một dải độ dài rất rộng.","metadata":{}},{"id":"6d740039","cell_type":"code","source":"print(\"AMINO ACID COMPOSITION ANALYSIS\")\nprint(\"_\"*80)\n\n# Calculate amino acid frequencies\ndef calculate_aa_composition(sequences):\n    \"\"\"Calculate amino acid composition from sequences\"\"\"\n    aa_counts = Counter()\n    total_aa = 0\n    \n    for seq in sequences.values():\n        aa_counts.update(seq)\n        total_aa += len(seq)\n    \n    aa_freq = {aa: count/total_aa for aa, count in aa_counts.items()}\n    return aa_freq\n\ntrain_aa_freq = calculate_aa_composition(train_sequences)\ntest_aa_freq = calculate_aa_composition(test_sequences)\n\n# Create comparison dataframe\naa_comparison = pd.DataFrame({\n    'Train': train_aa_freq,\n    'Test': test_aa_freq\n}).sort_values('Train', ascending=False)\n\nprint(\"\\nAmino Acid Frequencies (%):\")\nprint((aa_comparison * 100).round(2))\n\n# Visualize\nfig, axes = plt.subplots(1, 2, figsize=(18, 6))\n\n# Bar chart\nx = np.arange(len(aa_comparison))\nwidth = 0.35\n\nbars1 = axes[0].bar(x - width/2, aa_comparison['Train']*100, width, label='Train', color='steelblue', edgecolor='black')\nbars2 = axes[0].bar(x + width/2, aa_comparison['Test']*100, width, label='Test', color='coral', edgecolor='black')\n\naxes[0].set_xlabel('Amino Acid', fontsize=12)\naxes[0].set_ylabel('Frequency (%)', fontsize=12)\naxes[0].set_title('Amino Acid Composition: Train vs Test', fontsize=14, fontweight='bold')\naxes[0].set_xticks(x)\naxes[0].set_xticklabels(aa_comparison.index, fontsize=11)\naxes[0].legend()\naxes[0].grid(axis='y', alpha=0.3)\n\n# Scatter plot\naxes[1].scatter(aa_comparison['Train']*100, aa_comparison['Test']*100, s=100, alpha=0.6, color='purple', edgecolors='black', linewidth=1.5)\n\n# Add diagonal line\nmin_val = min(aa_comparison.min()) * 100\nmax_val = max(aa_comparison.max()) * 100\naxes[1].plot([min_val, max_val], [min_val, max_val], 'r--', linewidth=2, label='Perfect Match')\n\n# Add labels for each point\nfor aa, row in aa_comparison.iterrows():\n    axes[1].annotate(aa, (row['Train']*100, row['Test']*100), fontsize=10, ha='right')\n\naxes[1].set_xlabel('Train Frequency (%)', fontsize=12)\naxes[1].set_ylabel('Test Frequency (%)', fontsize=12)\naxes[1].set_title('Train vs Test AA Composition Correlation', fontsize=14, fontweight='bold')\naxes[1].legend()\naxes[1].grid(alpha=0.3)\n\nplt.tight_layout()\nplt.show()\n\nprint(f\"\\nInsight: Amino acid composition is very similar between train and test sets\")\nprint(f\"Most common amino acids: {', '.join(aa_comparison.head(5).index.tolist())}\")\nprint(f\"Least common amino acids: {', '.join(aa_comparison.tail(5).index.tolist())}\")","metadata":{},"outputs":[],"execution_count":null},{"id":"4c248e88","cell_type":"markdown","source":"### Nhận xét dựa trên AMINO ACID COMPOSITION ANALYSIS\n\n1. Thành phần axit amin của hai tập **rất giống nhau**, không có sự lệch đáng kể.\n2. Đặc điểm thành phần axit amin\n- Nhóm phổ biến nhất: L (Leucine), S (Serine), A (Alanine), E (Glutamic acid), G (Glycine) chiếm tỷ lệ cao nhất (6-10%), đóng vai trò quan trọng trong cấu trúc và chức năng protein.\n- Nhóm trung bình: V (Valine), K (Lysine), P (Proline), T (Threonine), R (Arginine), D (Aspartic acid), I (Isoleucine), Q (Glutamine), N (Asparagine) có tần suất từ 4-6%.\n- Nhóm hiếm: F (Phenylalanine), Y (Tyrosine), H (Histidine), M (Methionine), C (Cysteine), W (Tryptophan) xuất hiện với tần suất thấp (<4%).\n- Nhóm đặc biệt: X (axit amin không xác định), U (Selenocysteine), B (Aspartate/Asparagine), Z (Glutamate/Glutamine) có tần suất cực thấp hoặc không xuất hiện.\n","metadata":{}},{"id":"c1b75c5a","cell_type":"markdown","source":"### Ý nghĩa:\n\n1. Về mặt phân bố dữ liệu\n   - Không tồn tại distribution shift: Phân tích định lượng cho thấy tính đồng nhất giữa hai tập dữ liệu, với hệ số tương quan gần như tuyệt đối. Điều này xác nhận rằng quá trình phân chia dữ liệu đã được thực hiện hợp lý, đảm bảo tính đại diện và khả năng so sánh giữa các tập.\n   - Không cần điều chỉnh đặc biệt: Do tính đồng nhất cao, không cần áp dụng các kỹ thuật tái cân bằng dữ liệu (reweighting) hoặc điều chỉnh phân bố theo axit amin để xử lý sự chênh lệch giữa train và test.\n2. Về chiến lược embedding và biểu diễn đặc trưng\n   - Cân nhắc sử dụng embedding hoặc mô hình ngôn ngữ protein: Các phương pháp như one-hot encoding, k-mer, hoặc convolution trên chuỗi có thể chưa đủ để nắm bắt đầy đủ thông tin ngữ nghĩa và cấu trúc của protein. Việc sử dụng các mô hình tiền huấn luyện (pre-trained protein language models) như ESM, ProtBERT, hoặc ProtTrans có thể cải thiện đáng kể khả năng biểu diễn đặc trưng.\n   - Chuẩn hóa và scaling: Nếu sử dụng các đặc trưng dựa trên tần suất axit amin, cần áp dụng chuẩn hóa (normalization) hoặc scaling phù hợp để đảm bảo các thuật toán machine learning không bị ảnh hưởng bởi sự chênh lệch về mật độ giữa các axit amin.\n3. Xử lý các trường hợp đặc biệt\n   - Xử lý axit amin hiếm và không xác định (X, U, B, Z): Các axit amin này xuất hiện với tần suất cực thấp và có thể gây nhiễu cho mô hình. Các chiến lược xử lý bao gồm:\n      + Loại bỏ hoàn toàn các chuỗi chứa ký tự này nếu không ảnh hưởng đáng kể đến kích thước tập dữ liệu\n      + Gộp vào một token chung (unknown token) để giảm chiều dữ liệu\n      Ánh xạ các ký tự này về axit amin gần nhất về mặt hóa học hoặc cấu trúc\n   - Lựa chọn phù hợp với mục tiêu mô hình: Cách xử lý cần được lựa chọn dựa trên mục đích nghiên cứu và đặc thù của bộ dữ liệu, tránh làm mất thông tin quan trọng hoặc gây lỗi phân tích.\n   - Lưu ý quan trọng về 'U' (Selenocysteine): Mặc dù 'U' có thể được map về 'C' (Cysteine) do tính tương đồng, việc này làm mất thông tin về tính độc đáo sinh học (axit amin thứ 21). Nếu sử dụng Protein Language Model (PLM) tiền huấn luyện, nên kiểm tra xem mô hình đó có hỗ trợ 'U' như một token riêng biệt hay không. Nếu không, việc thay thế 'U' bằng 'C' hoặc 'X' là cần thiết.\n4. Chiến lược đánh giá và tối ưu hóa mô hình\n   - Sử dụng phân bố axit amin làm baseline: Phân bố axit amin có thể được sử dụng làm đặc trưng cơ bản để đánh giá hiệu suất ban đầu của mô hình. Nếu mô hình không thể dự đoán tốt hơn so với phân bố tần suất đơn thuần, cần xem xét lại cấu trúc và phương pháp học.\n   - Cân bằng giữa complexity và generalization: Do phân bố axit amin đã cân bằng tốt giữa train và test, việc tập trung vào các đặc trưng GO (Gene Ontology), cấu trúc bậc cao (secondary/tertiary structure), hoặc tương tác protein có thể mang lại giá trị gia tăng đáng kể thay vì chỉ dựa vào thành phần axit amin.\n   - Regularization và cross-validation: Mặc dù không có dấu hiệu distribution shift, vẫn cần áp dụng các kỹ thuật regularization (L1, L2, dropout) và cross-validation để đảm bảo mô hình không bị overfitting và có khả năng tổng quát hóa tốt.","metadata":{}},{"id":"73721095","cell_type":"markdown","source":"### Phân Tích Tương Quan Đặc Trưng\n\nChúng ta đã phân tích độc lập về độ phức tạp chức năng (số GO terms/protein) và độ phức tạp vật lý (độ dài chuỗi). Việc tiếp theo là khám phá xem có mối tương quan nào giữa hai đặc trưng này không. Thông thường, protein dài hơn (đa miền) sẽ có nhiều chức năng hơn.","metadata":{}},{"id":"7b576518","cell_type":"code","source":"# Phân tích độ dài chuỗi theo Aspect (P/C/F)\n\n# 1. Merge train_terms và train_seq_df\n# train_terms (EntryID, aspect) + train_seq_df (EntryID, seq_length)\ntrain_aspect_len = train_terms.merge(train_seq_df[['EntryID', 'seq_length']], \n                                     on='EntryID', how='left')\n\n# 2. Tính độ dài chuỗi trung bình và độ lệch chuẩn theo từng Aspect\naspect_length_stats = train_aspect_len.groupby('aspect')['seq_length'].agg(\n    ['mean', 'median', 'std', 'count']\n).sort_values('mean', ascending=False)\n\nprint(\"Thống kê độ dài chuỗi theo GO Aspect:\")\ndisplay(aspect_length_stats.round(2))\n\n# 3. Trực quan hóa (Box Plot)\nplt.figure(figsize=(8, 6))\n# Chỉ lấy mẫu 50,000 để tránh biểu đồ bị quá tải, đồng thời loại bỏ các outliers quá lớn\nsample_df = train_aspect_len.sample(min(50000, len(train_aspect_len)), random_state=42)\nsns.boxplot(x='aspect', y='seq_length', data=sample_df, order=['P', 'F', 'C'], palette='Set2', \n            showfliers=False) # Không hiển thị outliers để nhìn rõ phân bố\nplt.title('Phân phối Độ dài Chuỗi theo GO Aspect (Median/IQR)', fontsize=14, fontweight='bold')\nplt.xlabel('GO Aspect (P: Process, F: Function, C: Component)', fontsize=12)\nplt.ylabel('Sequence Length (aa)', fontsize=12)\nplt.grid(axis='y', alpha=0.3)\nplt.show()","metadata":{},"outputs":[],"execution_count":null},{"id":"dffa84ec","cell_type":"code","source":"# Tương quan giữa Sequence Length và Số GO Terms (Functional Complexity)\n\n# 1. Merge dữ liệu\n# train_seq_df đã có 'EntryID' và 'seq_length'\n# annotations_per_protein đã có 'EntryID' (Index) và số GO terms (Value)\ncorr_df = train_seq_df.merge(annotations_per_protein.rename('num_go_terms'), \n                             on='EntryID', \n                             how='left')\n\n# 2. Tính toán và in ra hệ số tương quan (Pearson)\ncorrelation = corr_df['seq_length'].corr(corr_df['num_go_terms'])\nprint(f\"Hệ số tương quan (Pearson) giữa Seq Length và Num GO Terms: {correlation:.4f}\")\n\n# 3. Trực quan hóa Scatter Plot (sử dụng log-scale cho độ dài chuỗi)\nplt.figure(figsize=(8, 6))\n# Cắt bớt các giá trị outlier (số GO terms > 50) để plot trực quan hơn\nplot_df = corr_df[corr_df['num_go_terms'] < 50]\n\nplt.scatter(plot_df['seq_length'], plot_df['num_go_terms'], alpha=0.3, s=10)\nplt.xscale('log') # Sử dụng log-scale cho Seq Length vì phân bố quá lệch phải\nplt.title('Tương quan giữa Độ dài Chuỗi và Số lượng GO Terms (Train)', \n          fontsize=14, fontweight='bold')\nplt.xlabel('Độ dài Chuỗi (log scale)', fontsize=12)\nplt.ylabel('Số lượng GO Terms per Protein', fontsize=12)\nplt.grid(alpha=0.3)\nplt.show()\n\n# 4. Nhận xét sơ bộ (sẽ được giải thích trong Markdown Cell tiếp theo)\nprint(\"\\nInsight: Mối tương quan yếu (0 < R < 0.4) nhưng vẫn có xu hướng: Protein dài hơn có xu hướng mang nhiều GO terms hơn.\")","metadata":{},"outputs":[],"execution_count":null},{"id":"07065706","cell_type":"markdown","source":"## Nhận Xét Tương Quan Giữa Độ dài Chuỗi và Số GO Terms\n\nDựa trên biểu đồ Scatter Plot (với trục x là Độ dài Chuỗi theo thang logarit, và trục y là Số lượng GO Terms):\n\n1.  **Mối tương quan tồn tại nhưng yếu (Weak Positive Correlation):**\n    * Mặc dù không có mối quan hệ tuyến tính rõ rệt, nhưng có thể thấy rõ ràng các protein có độ dài **dưới $10^2$ (dưới 100 aa)** hầu hết chỉ mang **dưới 10 GO terms**.\n    * Các protein có độ dài từ $10^2$ đến $10^3$ (100 đến 1000 aa) chiếm mật độ cao nhất và là nơi tập trung hầu hết các annotations (0-20 GO terms).\n    * Các protein dài hơn ($> 10^3$ aa) có **khả năng** mang nhiều GO terms hơn (lên đến 50), mặc dù mật độ dữ liệu ở vùng này thấp hơn.\n\n2.  **Mật độ Annotations tập trung ở chuỗi trung bình:**\n    * Phần lớn các protein chỉ được gán **dưới 10 GO terms** (tính đa nhãn thấp).\n    * Sự phân bố dày đặc tập trung vào vùng độ dài chuỗi phổ biến (100-1000 aa).\n\n3.  **Hạn chế của Độ dài Chuỗi:**\n    * Độ dài chuỗi tự thân **không phải là yếu tố dự đoán chính** cho sự phức tạp chức năng. Nhiều protein dài (ví dụ: $10^3$ đến $10^4$ aa) vẫn chỉ có 1-5 GO terms.\n    * Điều này gợi ý rằng **mức độ nghiên cứu/annotation bias** (protein người, nấm men, E.coli) và **tính đa miền (domain structure)** mới là yếu tố quan trọng hơn độ dài thuần túy.\n\n### Ứng dụng cho Feature Engineering và Mô hình hóa:\n\n* **Tạo đặc trưng:** Thay vì dùng độ dài thô, nên chuyển Độ dài Chuỗi thành các **bin phân loại** (dưới 100, 100-500, 500-1000, trên 1000) hoặc dùng logarit của độ dài để làm đặc trưng cho các mô hình Machine Learning cổ điển, vì nó phản ánh tốt hơn sự thay đổi mật độ GO terms ở các ngưỡng độ dài khác nhau.\n* **Chiến lược mô hình hóa:** Mô hình cần được đào tạo để nhận diện **protein đa miền** thông qua trình tự (Sequence) chứ không chỉ dựa vào Độ dài (Length). Độ dài chỉ nên là một đặc trưng hỗ trợ.","metadata":{}},{"id":"e055fdcb","cell_type":"markdown","source":"<a id='6-taxonomy-analysis'></a>\n## 6. Taxonomy Analysis\n\nPhân tích phân bố các loài sinh vật (taxonomy) trong dữ liệu:","metadata":{}},{"id":"d418e9a2","cell_type":"code","source":"print(\"🔬 TAXONOMY ANALYSIS\")\nprint(\"=\"*80)\n\n# Count taxonomy IDs\ntaxon_counts = train_taxonomy['taxon_id'].value_counts()\n\nprint(f\"\\nTotal unique organisms: {len(taxon_counts)}\")\nprint(\"\\nTop 20 organisms by protein count:\")\nprint(taxon_counts.head(20))\n\n# Common organism names (for major taxon IDs)\ntaxon_names = {\n    9606: \"Homo sapiens (Human)\",\n    10090: \"Mus musculus (Mouse)\",\n    7227: \"Drosophila melanogaster (Fruit fly)\",\n    6239: \"Caenorhabditis elegans (Worm)\",\n    7955: \"Danio rerio (Zebrafish)\",\n    559292: \"Saccharomyces cerevisiae (Yeast)\",\n    3702: \"Arabidopsis thaliana (Plant)\",\n    83333: \"Escherichia coli (Bacteria)\",\n}\n\n# Visualize top 20\nfig, axes = plt.subplots(2, 1, figsize=(16, 12))\n\n# Bar chart\ntop_20_taxon = taxon_counts.head(20)\nlabels_with_names = [f\"{tid}\\n{taxon_names.get(tid, 'Unknown')[:30]}\" if tid in taxon_names else str(tid) \n                     for tid in top_20_taxon.index]\n\nbars = axes[0].bar(range(20), top_20_taxon.values, color='forestgreen', edgecolor='black', alpha=0.7)\naxes[0].set_xticks(range(20))\naxes[0].set_xticklabels(labels_with_names, rotation=45, ha='right', fontsize=9)\naxes[0].set_ylabel('Number of Proteins', fontsize=12)\naxes[0].set_title('Top 20 Organisms by Protein Count', fontsize=14, fontweight='bold')\naxes[0].grid(axis='y', alpha=0.3)\n\n# Add value labels\nfor i, bar in enumerate(bars):\n    height = bar.get_height()\n    axes[0].text(bar.get_x() + bar.get_width()/2., height,\n                f'{int(height):,}',\n                ha='center', va='bottom', fontsize=9)\n\n# Distribution of proteins per organism\naxes[1].hist(taxon_counts, bins=50, edgecolor='black', alpha=0.7, color='teal')\naxes[1].axvline(taxon_counts.mean(), color='red', linestyle='--', linewidth=2, label=f'Mean: {taxon_counts.mean():.1f}')\naxes[1].axvline(taxon_counts.median(), color='green', linestyle='--', linewidth=2, label=f'Median: {taxon_counts.median():.1f}')\naxes[1].set_xlabel('Number of Proteins per Organism', fontsize=12)\naxes[1].set_ylabel('Number of Organisms', fontsize=12)\naxes[1].set_title('Distribution of Proteins per Organism', fontsize=14, fontweight='bold')\naxes[1].set_yscale('log')\naxes[1].legend()\naxes[1].grid(axis='y', alpha=0.3)\n\nplt.tight_layout()\nplt.show()\n\nprint(f\"\\n💡 Insight: Data is heavily dominated by {taxon_names.get(taxon_counts.index[0], 'Unknown organism')}\")\nprint(f\"💡 Insight: Top organism has {taxon_counts.values[0]:,} proteins ({taxon_counts.values[0]/len(train_taxonomy)*100:.1f}% of total)\")","metadata":{},"outputs":[],"execution_count":null},{"id":"839b3a6b","cell_type":"markdown","source":"## 6.1 Phân Tích Annotation Bias theo Taxonomy\n\nPhân tích trước đó đã chỉ ra rằng dữ liệu bị chi phối bởi một số lượng nhỏ các sinh vật (Taxon Head). Điều này thường dẫn đến **thiên vị annotation (annotation bias)**: các protein từ các loài được nghiên cứu kỹ (Human, Mouse) thường được gán nhiều GO terms hơn.\n\nChúng ta sẽ xác nhận mối quan hệ này bằng cách xem xét số lượng GO Terms trung bình được gán cho 10 Taxon ID hàng đầu.","metadata":{}},{"id":"a5c3ed14","cell_type":"code","source":"# Tính Số GO Terms trung bình theo Top Taxon ID\n\n# 1. Kết hợp train_terms (EntryID, term) và train_taxonomy (EntryID, taxon_id)\ntrain_full = train_terms.merge(train_taxonomy, on='EntryID', how='left')\n\n# 2. Tính số lượng GO terms trên mỗi protein\ngo_per_protein = train_full.groupby('EntryID')['term'].count().reset_index()\ngo_per_protein.rename(columns={'term': 'num_go_terms'}, inplace=True)\n\n# 3. Kết hợp lại với Taxon ID\ntaxon_complexity = go_per_protein.merge(train_taxonomy, on='EntryID', how='left')\n\n# 4. Tính toán số GO terms trung bình cho mỗi Taxon ID\navg_go_per_taxon = taxon_complexity.groupby('taxon_id')['num_go_terms'].agg(\n    ['mean', 'count']).reset_index()\navg_go_per_taxon.rename(columns={'count': 'protein_count', 'mean': 'avg_go_terms'}, inplace=True)\n\n# 5. Lọc ra 10 Taxon ID hàng đầu (theo protein_count)\ntop_10_avg_go = avg_go_per_taxon.nlargest(10, 'protein_count')\n\n# 6. Thêm tên loài để dễ đọc\ntaxon_names = {\n    9606: \"Homo sapiens (Human)\", 10090: \"Mus musculus (Mouse)\", 7227: \"Drosophila melanogaster\",\n    6239: \"Caenorhabditis elegans\", 7955: \"Danio rerio\", 559292: \"S. cerevisiae (Yeast)\",\n    3702: \"A. thaliana (Plant)\", 83333: \"E. coli (Bacteria)\", 10116: \"Rattus norvegicus\",\n    284812: \"Burkholderia pseudomallei\", 17162: \"Zea mays (Corn)\" # Thêm Taxon 17162 (top 2)\n}\ntop_10_avg_go['Organism'] = top_10_avg_go['taxon_id'].map(taxon_names)\n\nprint(\"Số lượng GO Terms trung bình theo 10 loài hàng đầu:\")\ndisplay(top_10_avg_go.round(2))\n\n# 7. Trực quan hóa (Bar chart)\nplt.figure(figsize=(12, 6))\nsns.barplot(x='Organism', y='avg_go_terms', data=top_10_avg_go, palette='viridis')\nplt.title('GO Terms Trung bình trên mỗi Protein theo Top 10 Taxon ID', fontsize=14, fontweight='bold')\nplt.ylabel('Số GO Terms Trung bình', fontsize=12)\nplt.xlabel('Tên Sinh vật', fontsize=12)\nplt.xticks(rotation=45, ha='right')\nplt.grid(axis='y', alpha=0.3)\nplt.tight_layout()\nplt.show()","metadata":{},"outputs":[],"execution_count":null},{"id":"efafe8c0","cell_type":"markdown","source":"## Nhận Xét Về Thiên Vị Annotation (Annotation Bias)\n\nDựa trên bảng và biểu đồ **GO Terms Trung bình trên mỗi Protein theo Top 10 Taxon ID**:\n\n1.  **Xác nhận Thiên vị Annotation (Annotation Bias):**\n    * Các loài mô hình (Model Organisms) như **Homo sapiens (Human, 9606)** và **Mus musculus (Mouse, 10090)** có số lượng GO Terms trung bình trên mỗi protein cao nhất (lần lượt là $9.05$ và $8.77$ GO terms/protein).\n    * Các loài ít được nghiên cứu hơn hoặc vi khuẩn, như **Burkholderia pseudomallei (284812)** hay **E. coli (Bacteria, 83333)**, có số lượng GO Terms trung bình thấp hơn đáng kể (lần lượt là $4.33$ và $4.27$ GO terms/protein).\n\n2.  **Sự phân tách giữa Eukaryote và Prokaryote:**\n    * Các sinh vật nhân chuẩn (Eukaryotes) chiếm đa số trong nhóm có độ phức tạp chức năng cao: Human ($9.05$), Mouse ($8.77$), Rat ($8.84$), Fruit fly ($8.76$), và Worm ($5.90$).\n    * Ngược lại, các sinh vật nhân sơ (Prokaryotes/Bacteria) như *E. coli* và *Burkholderia* nằm ở nhóm dưới, có số GO Terms trung bình chỉ khoảng 4. Điều này phản ánh sự khác biệt về mức độ nghiên cứu và sự phức tạp sinh học của hai nhóm này.\n\n3.  **Ý nghĩa đối với Mô hình hóa:**\n    * **Overfitting cục bộ:** Mô hình có xu hướng học các GO terms phức tạp và cụ thể từ các protein của sinh vật mẫu (Human, Mouse).\n    * **Thách thức Generalization:** Việc dự đoán chức năng cho protein từ các loài ít được annotate (có **protein\\_count** thấp) sẽ rất khó khăn, do mô hình chỉ nhận được các tín hiệu học tập yếu (ít nhãn) từ nhóm này.\n    * **Chiến lược cần thiết:** Cần sử dụng **Taxonomy ID** như một đặc trưng (feature) hoặc áp dụng **Taxonomy-aware Loss Weighting** để cân bằng tầm quan trọng giữa các loài, tránh để mô hình bị chi phối hoàn toàn bởi dữ liệu từ Human và Mouse.","metadata":{}}]}