{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.12.12","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":31153,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## **CAFA 6 PROTEIN FUNCTION PREDICTION - COMPREHENSIVE FEATURE ENGINEERING**\n\n>**Competition:** Critical Assessment of Functional Annotation 6\n\n>**Task:** Predict Gene Ontology (GO) terms for proteins based on amino acid sequences\n\n>**Evaluation:** Weighted F1-score (Information Accretion) across MFO, BPO, CCO ontologies\n\n>**Author:** Dr. Hany Ghazal (PhD)\n\n>**Date:** December 2025\n\n>**Version: 3.0 - Production Ready**\n\n### ***NOTEBOOK STRUCTURE:***\n1. Setup & Data Loading\n2. GO Ontology Analysis\n3. Basic Sequence Features (Length, Composition, Molecular Properties)\n4. Advanced Physicochemical Features\n5. K-mer Features (Dipeptides, Tripeptides)\n6. Positional & Regional Features\n7. GO Hierarchy Features (IC, Depth, Ancestors)\n8. Sequence Pattern Detection\n9. Feature Engineering Pipeline\n10. Feature Selection & Optimization\n11. Export for Modeling\n\n### ***KEY FEATURES:***\n>- GO hierarchy-aware features using ontology structure\n>- Information Accretion (IA) weighted features\n>- Multi-ontology support (MFO, BPO, CCO)\n>- Memory-efficient batch processing\n>- Comprehensive feature validation\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"# !pip install obonet","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-31T00:21:47.659685Z","iopub.execute_input":"2025-12-31T00:21:47.660381Z","iopub.status.idle":"2025-12-31T00:21:47.663737Z","shell.execute_reply.started":"2025-12-31T00:21:47.660348Z","shell.execute_reply":"2025-12-31T00:21:47.663000Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport networkx as nx\nfrom collections import Counter, defaultdict\nfrom typing import Dict, List, Tuple, Set\nimport warnings\nfrom pathlib import Path\nimport gc\nimport obonet\n\nwarnings.filterwarnings('ignore')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-31T00:27:03.227152Z","iopub.execute_input":"2025-12-31T00:27:03.227827Z","iopub.status.idle":"2025-12-31T00:27:03.232058Z","shell.execute_reply.started":"2025-12-31T00:27:03.227799Z","shell.execute_reply":"2025-12-31T00:27:03.231274Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================================\n# 1. CONFIGURATION & SETUP\n# ============================================================================\n\nclass Config:\n    \"\"\"Competition-specific configuration.\"\"\"\n    \n    # File paths\n    DATA_DIR = Path('/kaggle/input/cafa-6-protein-function-prediction')\n    TRAIN_SEQUENCES = DATA_DIR / 'Train/train_sequences.fasta'\n    TRAIN_TERMS = DATA_DIR / 'Train/train_terms.tsv'\n    TRAIN_TAXONOMY = DATA_DIR / 'Train/train_taxonomy.tsv'\n    GO_GRAPH = DATA_DIR / 'Train/go-basic.obo'\n    TEST_SEQUENCES = DATA_DIR / 'Test/testsuperset.fasta'\n    IA_WEIGHTS = DATA_DIR / 'IA.tsv'\n    \n    # Ontology roots (from dataset description)\n    ONTOLOGY_ROOTS = {\n        'BPO': 'GO:0008150',  # Biological Process\n        'MFO': 'GO:0003674',  # Molecular Function\n        'CCO': 'GO:0005575'   # Cellular Component\n    }\n    \n    # Feature engineering settings\n    INCLUDE_DIPEPTIDES = True      # 400 features\n    INCLUDE_TRIPEPTIDES = False    # Too many features, use sparse\n    TRIPEPTIDE_TOP_N = 100         # Only top N if enabled\n    TERMINAL_WINDOW = 25           # N/C terminal window size\n    N_REGIONS = 5                  # Number of sequence regions\n    \n    # Memory management\n    BATCH_SIZE = 5000\n    VERBOSE = True\n\nprint(\"=\"*80)\nprint(\" \"*20 + \"CAFA 6 FEATURE ENGINEERING PIPELINE\")\nprint(\"=\"*80)\nprint(f\"\\nConfiguration loaded:\")\nprint(f\"  - Dipeptides: {Config.INCLUDE_DIPEPTIDES}\")\nprint(f\"  - Tripeptides: {Config.INCLUDE_TRIPEPTIDES}\")\nprint(f\"  - Batch size: {Config.BATCH_SIZE}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-31T00:27:03.635378Z","iopub.execute_input":"2025-12-31T00:27:03.635627Z","iopub.status.idle":"2025-12-31T00:27:03.642295Z","shell.execute_reply.started":"2025-12-31T00:27:03.635604Z","shell.execute_reply":"2025-12-31T00:27:03.641547Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================================\n# 2. DATA LOADING UTILITIES\n# ============================================================================\n\ndef load_fasta_sequences(fasta_path: Path) -> Dict[str, str]:\n    \"\"\"\n    Load protein sequences from FASTA file.\n    Format: >sp|PROTEIN_ID|NAME or >PROTEIN_ID\n    \"\"\"\n    sequences = {}\n    current_id = None\n    current_seq = []\n    \n    with open(fasta_path, 'r') as f:\n        for line in f:\n            line = line.strip()\n            if line.startswith('>'):\n                # Save previous sequence\n                if current_id:\n                    sequences[current_id] = ''.join(current_seq)\n                \n                # Parse new header\n                # Format: sp|P9WHI7|RECN_MYCT or just PROTEIN_ID\n                header = line[1:]\n                if '|' in header:\n                    parts = header.split('|')\n                    current_id = parts[1]  # Extract UniProt ID\n                else:\n                    current_id = header.split()[0]\n                \n                current_seq = []\n            else:\n                current_seq.append(line)\n        \n        # Save last sequence\n        if current_id:\n            sequences[current_id] = ''.join(current_seq)\n    \n    return sequences\n\n\n\ndef load_go_annotations(terms_path: Path) -> pd.DataFrame:\n    \"\"\"Load GO term annotations.\"\"\"\n    df = pd.read_csv(terms_path, sep='\\t', \n                     names=['protein_id', 'go_term', 'ontology'])\n    return df\n\ndef load_taxonomy(taxonomy_path: Path) -> pd.DataFrame:\n    \"\"\"Load taxonomy information.\"\"\"\n    df = pd.read_csv(taxonomy_path, sep='\\t',\n                     names=['protein_id', 'taxon_id'])\n    return df\n\ndef load_ia_weights(ia_path: Path) -> Dict[str, float]:\n    \"\"\"Load Information Accretion weights.\"\"\"\n    df = pd.read_csv(ia_path, sep='\\t',\n                     names=['go_term', 'ia_weight'])\n    return df.set_index('go_term')['ia_weight'].to_dict()\n\ndef load_go_graph(obo_path: Path) -> nx.DiGraph:\n    \"\"\"\n    Load GO ontology graph from OBO file.\n    Note: Requires obonet library (pip install obonet)\n    \"\"\"\n    try:\n        import obonet\n        graph = obonet.read_obo(str(obo_path))\n        return graph\n    except ImportError:\n        print(\"Warning: obonet not installed. Using simplified graph loading.\")\n        return create_simple_go_graph(obo_path)\n\ndef create_simple_go_graph(obo_path: Path) -> nx.DiGraph:\n    \"\"\"Simplified GO graph creation without obonet.\"\"\"\n    graph = nx.DiGraph()\n    current_term = None\n    \n    with open(obo_path, 'r') as f:\n        for line in f:\n            line = line.strip()\n            if line.startswith('id: GO:'):\n                current_term = line.split('id: ')[1]\n                graph.add_node(current_term)\n            elif line.startswith('is_a: GO:') and current_term:\n                parent = line.split('is_a: ')[1].split(' !')[0]\n                graph.add_edge(current_term, parent)\n    \n    return graph\n\nprint(\"\\n\" + \"=\"*80)\nprint(\"LOADING DATA...\")\nprint(\"=\"*80)\n# Load all data\ntrain_sequences = load_fasta_sequences(Config.TRAIN_SEQUENCES)\ntrain_annotations = load_go_annotations(Config.TRAIN_TERMS)\ntrain_taxonomy = load_taxonomy(Config.TRAIN_TAXONOMY)\nia_weights = load_ia_weights(Config.IA_WEIGHTS)\ngo_graph = load_go_graph(Config.GO_GRAPH)\n\nprint(f\"\\n✓ Data loaded successfully:\")\nprint(f\"  - Training sequences: {len(train_sequences):,}\")\nprint(f\"  - Unique proteins: {train_annotations['protein_id'].nunique():,}\")\nprint(f\"  - Total annotations: {len(train_annotations):,}\")\nprint(f\"  - Unique GO terms: {train_annotations['go_term'].nunique():,}\")\nprint(f\"  - GO graph nodes: {go_graph.number_of_nodes():,}\")\nprint(f\"  - GO graph edges: {go_graph.number_of_edges():,}\")\nprint(f\"  - IA weights: {len(ia_weights):,}\")\n\n# Ontology distribution\nprint(f\"\\n  Ontology distribution:\")\nfor ont in ['BPO', 'MFO', 'CCO']:\n    count = (train_annotations['ontology'] == ont).sum()\n    print(f\"    {ont}: {count:,} ({count/len(train_annotations)*100:.1f}%)\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-31T00:27:04.082357Z","iopub.execute_input":"2025-12-31T00:27:04.083026Z","iopub.status.idle":"2025-12-31T00:27:09.770577Z","shell.execute_reply.started":"2025-12-31T00:27:04.083001Z","shell.execute_reply":"2025-12-31T00:27:09.769962Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================================\n# 3. AMINO ACID PROPERTIES DATABASE\n# ============================================================================\n\nclass AAProperties:\n    \"\"\"Comprehensive amino acid physicochemical properties database.\"\"\"\n    \n    # All 20 standard amino acids\n    AMINO_ACIDS = ['A', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'K', 'L',\n                   'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'V', 'W', 'Y']\n    \n    # Molecular weights (Da)\n    MW = {'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        # Kyte-Doolittle hydrophobicity scale\n    HYDRO = {'A': 1.8, 'C': 2.5, 'D': -3.5, 'E': -3.5, 'F': 2.8,\n             'G': -0.4, 'H': -3.2, 'I': 4.5, 'K': -3.9, 'L': 3.8,\n             'M': 1.9, 'N': -3.5, 'P': -1.6, 'Q': -3.5, 'R': -4.5,\n             'S': -0.8, 'T': -0.7, 'V': 4.2, 'W': -0.9, 'Y': -1.3}\n    \n    # Charge at pH 7.0\n    CHARGE = {'A': 0, 'C': 0, 'D': -1, 'E': -1, 'F': 0, 'G': 0,\n              'H': 0.1, 'I': 0, 'K': 1, 'L': 0, 'M': 0, 'N': 0,\n              'P': 0, 'Q': 0, 'R': 1, 'S': 0, 'T': 0, 'V': 0,\n              'W': 0, 'Y': 0}\n    \n    # Polarity (1=polar, 0=nonpolar)\n    POLAR = {'A': 0, 'C': 1, 'D': 1, 'E': 1, 'F': 0, 'G': 0,\n             'H': 1, 'I': 0, 'K': 1, 'L': 0, 'M': 0, 'N': 1,\n             'P': 0, 'Q': 1, 'R': 1, 'S': 1, 'T': 1, 'V': 0,\n             'W': 0, 'Y': 1}\n\n      \n    # Van der Waals volume (ų)\n    VOLUME = {'A': 67, 'C': 86, 'D': 91, 'E': 109, 'F': 135, 'G': 48,\n              'H': 118, 'I': 124, 'K': 135, 'L': 124, 'M': 124, 'N': 96,\n              'P': 90, 'Q': 114, 'R': 148, 'S': 73, 'T': 93, 'V': 105,\n              'W': 163, 'Y': 141}\n    \n    # Aromaticity\n    AROMATIC = {'F': 1, 'W': 1, 'Y': 1, 'H': 1}\n    \n    # Functional groups\n    GROUPS = {\n        'aliphatic': ['A', 'V', 'L', 'I', 'M'],\n        'aromatic': ['F', 'Y', 'W', 'H'],\n        'polar_uncharged': ['S', 'T', 'N', 'Q', 'C'],\n        'positive': ['K', 'R', 'H'],\n        'negative': ['D', 'E'],\n        'small': ['A', 'G', 'S', 'C'],\n        'tiny': ['G', 'A', 'S'],\n        'sulfur': ['C', 'M'],\n        'hydroxyl': ['S', 'T', 'Y'],\n        'amide': ['N', 'Q']\n    }\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-31T00:27:57.297537Z","iopub.execute_input":"2025-12-31T00:27:57.298265Z","iopub.status.idle":"2025-12-31T00:27:57.307554Z","shell.execute_reply.started":"2025-12-31T00:27:57.298236Z","shell.execute_reply":"2025-12-31T00:27:57.306904Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================================\n# 4. FEATURE EXTRACTION: BASIC SEQUENCE FEATURES\n# ============================================================================\n\ndef extract_basic_features(seq: str) -> Dict[str, float]:\n    \"\"\"Extract basic sequence-level features.\"\"\"\n    features = {}\n    length = len(seq)\n    \n    # Length features\n    features['length'] = length\n    features['log_length'] = np.log1p(length)\n    features['sqrt_length'] = np.sqrt(length)\n    \n    # Molecular weight\n    if length > 0:\n        mw = sum(AAProperties.MW.get(aa, 0) for aa in seq)\n        mw -= (length - 1) * 18.015  # Subtract water from peptide bonds\n        features['molecular_weight'] = mw\n        features['mw_per_aa'] = mw / length\n    else:\n        features['molecular_weight'] = 0\n        features['mw_per_aa'] = 0\n    \n    # Amino acid composition (20 features)\n    counter = Counter(seq)\n    for aa in AAProperties.AMINO_ACIDS:\n        features[f'comp_{aa}'] = counter.get(aa, 0) / length if length > 0 else 0\n    \n    # Group composition (10 features)\n    for group_name, group_aas in AAProperties.GROUPS.items():\n        count = sum(counter.get(aa, 0) for aa in group_aas)\n        features[f'group_{group_name}'] = count / length if length > 0 else 0\n    \n    return features","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-31T00:27:58.545806Z","iopub.execute_input":"2025-12-31T00:27:58.546088Z","iopub.status.idle":"2025-12-31T00:27:58.553154Z","shell.execute_reply.started":"2025-12-31T00:27:58.546065Z","shell.execute_reply":"2025-12-31T00:27:58.552305Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================================\n# 5. FEATURE EXTRACTION: PHYSICOCHEMICAL PROPERTIES\n# ============================================================================\n\ndef extract_physicochemical_features(seq: str) -> Dict[str, float]:\n    \"\"\"Extract physicochemical property statistics.\"\"\"\n    features = {}\n    \n    if len(seq) == 0:\n        return {f'phys_{prop}_{stat}': 0.0 \n                for prop in ['hydro', 'charge', 'volume', 'polar']\n                for stat in ['mean', 'std', 'min', 'max']}\n    \n    # Hydrophobicity statistics\n    hydro_values = [AAProperties.HYDRO.get(aa, 0) for aa in seq]\n    features['phys_hydro_mean'] = np.mean(hydro_values)\n    features['phys_hydro_std'] = np.std(hydro_values)\n    features['phys_hydro_min'] = np.min(hydro_values)\n    features['phys_hydro_max'] = np.max(hydro_values)\n    \n    # Charge statistics\n    charge_values = [AAProperties.CHARGE.get(aa, 0) for aa in seq]\n    features['phys_charge_mean'] = np.mean(charge_values)\n    features['phys_charge_std'] = np.std(charge_values)\n    features['phys_charge_min'] = np.min(charge_values)\n    features['phys_charge_max'] = np.max(charge_values)\n    features['phys_net_charge'] = sum(charge_values)\n    features['phys_charge_density'] = sum(charge_values) / len(seq)\n\n    # Volume statistics\n    volume_values = [AAProperties.VOLUME.get(aa, 0) for aa in seq]\n    features['phys_volume_mean'] = np.mean(volume_values)\n    features['phys_volume_std'] = np.std(volume_values)\n    \n    # Polarity statistics\n    polar_values = [AAProperties.POLAR.get(aa, 0) for aa in seq]\n    features['phys_polarity'] = np.mean(polar_values)\n    \n    # Aromaticity\n    aromatic_count = sum(1 for aa in seq if aa in AAProperties.AROMATIC)\n    features['phys_aromaticity'] = aromatic_count / len(seq)\n\n   # Isoelectric point (simplified estimate)\n    pos_count = sum(1 for aa in seq if aa in ['K', 'R', 'H'])\n    neg_count = sum(1 for aa in seq if aa in ['D', 'E'])\n    if pos_count + neg_count > 0:\n        pI = 6.5 + 0.5 * (pos_count - neg_count) / len(seq) * 10\n        features['phys_pI'] = np.clip(pI, 3.0, 11.0)\n    else:\n        features['phys_pI'] = 7.0\n    \n    return features\n    ","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-31T00:28:05.047655Z","iopub.execute_input":"2025-12-31T00:28:05.048053Z","iopub.status.idle":"2025-12-31T00:28:05.056682Z","shell.execute_reply.started":"2025-12-31T00:28:05.048021Z","shell.execute_reply":"2025-12-31T00:28:05.055914Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================================\n# 6. FEATURE EXTRACTION: K-MER FEATURES\n# ============================================================================\n\ndef extract_kmer_features(seq: str, k: int = 2) -> Dict[str, float]:\n    \"\"\"Extract k-mer frequency features.\"\"\"\n    features = {}\n    \n    if len(seq) < k:\n        return features\n    \n    # Count k-mers\n    kmer_counts = Counter()\n    for i in range(len(seq) - k + 1):\n        kmer = seq[i:i+k]\n        if all(aa in AAProperties.AMINO_ACIDS for aa in kmer):\n            kmer_counts[kmer] += 1\n    \n    # Normalize\n    total = sum(kmer_counts.values())\n    if total > 0:\n        for kmer, count in kmer_counts.items():\n            features[f'kmer{k}_{kmer}'] = count / total\n    \n    return features","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-31T00:28:05.363614Z","iopub.execute_input":"2025-12-31T00:28:05.364138Z","iopub.status.idle":"2025-12-31T00:28:05.369228Z","shell.execute_reply.started":"2025-12-31T00:28:05.364111Z","shell.execute_reply":"2025-12-31T00:28:05.368581Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================================\n# 7. FEATURE EXTRACTION: POSITIONAL FEATURES\n# ============================================================================\n\ndef extract_positional_features(seq: str) -> Dict[str, float]:\n    \"\"\"Extract position-specific features.\"\"\"\n    features = {}\n    length = len(seq)\n    \n    if length == 0:\n        return features\n    \n    # N-terminal composition (first 25 AA)\n    n_term_size = min(Config.TERMINAL_WINDOW, length)\n    n_term = seq[:n_term_size]\n    n_counter = Counter(n_term)\n    for aa in AAProperties.AMINO_ACIDS:\n        features[f'nterm_{aa}'] = n_counter.get(aa, 0) / n_term_size\n\n    # C-terminal composition (last 25 AA)\n    c_term_size = min(Config.TERMINAL_WINDOW, length)\n    c_term = seq[-c_term_size:]\n    c_counter = Counter(c_term)\n    for aa in AAProperties.AMINO_ACIDS:\n        features[f'cterm_{aa}'] = c_counter.get(aa, 0) / c_term_size\n    \n    # Regional properties (divide sequence into regions)\n    region_size = length // Config.N_REGIONS\n    for i in range(Config.N_REGIONS):\n        start = i * region_size\n        end = start + region_size if i < Config.N_REGIONS - 1 else length\n        region = seq[start:end]\n        \n        if len(region) > 0:\n            # Hydrophobicity in region\n            hydro = np.mean([AAProperties.HYDRO.get(aa, 0) for aa in region])\n            features[f'region{i}_hydro'] = hydro\n            \n            # Charge in region\n            charge = np.sum([AAProperties.CHARGE.get(aa, 0) for aa in region])\n            features[f'region{i}_charge'] = charge\n    \n    return features","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-31T00:28:12.130205Z","iopub.execute_input":"2025-12-31T00:28:12.130729Z","iopub.status.idle":"2025-12-31T00:28:12.137342Z","shell.execute_reply.started":"2025-12-31T00:28:12.130702Z","shell.execute_reply":"2025-12-31T00:28:12.136619Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================================\n# 8. FEATURE EXTRACTION: GO HIERARCHY FEATURES\n# ============================================================================\n\ndef extract_go_features(protein_id: str, \n                       annotations: pd.DataFrame,\n                       go_graph: nx.DiGraph,\n                       ia_weights: Dict[str, float]) -> Dict[str, float]:\n    \"\"\"Extract GO ontology-based features (CRITICAL for CAFA).\"\"\"\n    features = {}\n    \n    # Get protein's GO terms\n    protein_annot = annotations[annotations['protein_id'] == protein_id]\n    go_terms = protein_annot['go_term'].tolist()\n    \n    if len(go_terms) == 0:\n        return {\n            'go_num_terms': 0, 'go_num_terms_mfo': 0, 'go_num_terms_bpo': 0, \n            'go_num_terms_cco': 0, 'go_avg_depth': 0, 'go_max_depth': 0,\n            'go_avg_ia': 0, 'go_max_ia': 0, 'go_sum_ia': 0,\n            'go_num_ancestors': 0, 'go_specificity_score': 0\n        }\n    \n    # Basic counts\n    features['go_num_terms'] = len(go_terms)\n\n    # Per-ontology counts\n    for ont in ['MFO', 'BPO', 'CCO']:\n        ont_count = (protein_annot['ontology'] == ont).sum()\n        features[f'go_num_terms_{ont.lower()}'] = ont_count\n        features[f'go_ratio_{ont.lower()}'] = ont_count / len(go_terms)\n    \n    # Calculate term depths from roots\n    depths = []\n    for term in go_terms:\n        for root_id, root_name in Config.ONTOLOGY_ROOTS.items():\n            try:\n                if nx.has_path(go_graph, term, root_id):\n                    depth = nx.shortest_path_length(go_graph, term, root_id)\n                    depths.append(depth)\n                    break\n            except:\n                continue\n\n    if depths:\n        features['go_avg_depth'] = np.mean(depths)\n        features['go_max_depth'] = np.max(depths)\n        features['go_min_depth'] = np.min(depths)\n        features['go_std_depth'] = np.std(depths)\n    else:\n        features['go_avg_depth'] = 0\n        features['go_max_depth'] = 0\n        features['go_min_depth'] = 0\n        features['go_std_depth'] = 0\n    \n    # Information Accretion statistics (CRITICAL for CAFA evaluation)\n    ia_values = [ia_weights.get(term, 0) for term in go_terms]\n    features['go_avg_ia'] = np.mean(ia_values) if ia_values else 0\n    features['go_max_ia'] = np.max(ia_values) if ia_values else 0\n    features['go_sum_ia'] = np.sum(ia_values) if ia_values else 0\n    features['go_std_ia'] = np.std(ia_values) if ia_values else 0\n    \n    # Ancestor count (reflects term specificity)\n    all_ancestors = set()\n    for term in go_terms:\n        try:\n            ancestors = nx.ancestors(go_graph, term)\n            all_ancestors.update(ancestors)\n        except:\n            continue\n    \n    features['go_num_ancestors'] = len(all_ancestors)\n\n    # Specificity score (deeper terms are more specific)\n    if len(go_terms) > 0 and depths:\n        features['go_specificity_score'] = np.mean(depths) * np.mean(ia_values) if ia_values else 0\n    else:\n        features['go_specificity_score'] = 0\n    \n    return features","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-31T00:28:31.591245Z","iopub.execute_input":"2025-12-31T00:28:31.591757Z","iopub.status.idle":"2025-12-31T00:28:31.601762Z","shell.execute_reply.started":"2025-12-31T00:28:31.591732Z","shell.execute_reply":"2025-12-31T00:28:31.601101Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================================\n# 9. FEATURE EXTRACTION: SEQUENCE PATTERNS\n# ============================================================================\n\ndef extract_pattern_features(seq: str) -> Dict[str, float]:\n    \"\"\"Extract sequence pattern and motif features.\"\"\"\n    features = {}\n    \n    if len(seq) == 0:\n        return features\n    \n    # Detect repeating patterns\n    repeat_counts = {3: 0, 4: 0, 5: 0, 6: 0}\n    for length in repeat_counts.keys():\n        seen = set()\n        for i in range(len(seq) - 2*length + 1):\n            pattern = seq[i:i+length]\n            next_pattern = seq[i+length:i+2*length] if i+2*length <= len(seq) else ''\n            if pattern == next_pattern and pattern not in seen:\n                repeat_counts[length] += 1\n                seen.add(pattern)\n    \n    for length, count in repeat_counts.items():\n        features[f'pattern_repeat{length}'] = count\n    features['pattern_total_repeats'] = sum(repeat_counts.values())\n\n    # Signal peptide indicator (hydrophobic N-terminus)\n    if len(seq) >= 30:\n        n_term_hydro = np.mean([AAProperties.HYDRO.get(aa, 0) for aa in seq[:30]])\n        features['pattern_signal_peptide'] = 1 if n_term_hydro > 1.5 else 0\n    else:\n        features['pattern_signal_peptide'] = 0\n    \n    # Transmembrane helices (long hydrophobic stretches)\n    hydrophobic = {'A', 'V', 'L', 'I', 'M', 'F', 'W', 'P'}\n    max_stretch = current_stretch = 0\n    for aa in seq:\n        if aa in hydrophobic:\n            current_stretch += 1\n            max_stretch = max(max_stretch, current_stretch)\n        else:\n            current_stretch = 0\n    features['pattern_max_hydrophobic_stretch'] = max_stretch\n    features['pattern_has_transmembrane'] = 1 if max_stretch >= 20 else 0\n\n    # Nuclear localization signal (K/R rich regions)\n    has_nls = 0\n    for i in range(len(seq) - 4):\n        window = seq[i:i+5]\n        if sum(1 for aa in window if aa in ['K', 'R']) >= 4:\n            has_nls = 1\n            break\n    features['pattern_nuclear_localization'] = has_nls\n    \n    # Glycosylation sites (N-X-S/T motif)\n    glyco_sites = 0\n    for i in range(len(seq) - 2):\n        if seq[i] == 'N' and seq[i+1] != 'P' and seq[i+2] in ['S', 'T']:\n            glyco_sites += 1\n    features['pattern_glycosylation_sites'] = glyco_sites\n    \n    # Low complexity regions (repeating single AA)\n    max_single_repeat = max((Counter(seq).values()), default=0)\n    features['pattern_max_aa_repeat'] = max_single_repeat\n    features['pattern_complexity_ratio'] = max_single_repeat / len(seq)\n    \n    return features","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-31T00:28:37.758354Z","iopub.execute_input":"2025-12-31T00:28:37.759038Z","iopub.status.idle":"2025-12-31T00:28:37.768462Z","shell.execute_reply.started":"2025-12-31T00:28:37.759008Z","shell.execute_reply":"2025-12-31T00:28:37.767716Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================================\n# 10. MAIN FEATURE EXTRACTION PIPELINE\n# ============================================================================\n\nclass CAFAFeatureExtractor:\n    \"\"\"Main feature extraction pipeline for CAFA 6 competition.\"\"\"\n    \n    def __init__(self, \n                 annotations: pd.DataFrame,\n                 go_graph: nx.DiGraph,\n                 ia_weights: Dict[str, float],\n                 include_dipeptides: bool = True):\n        self.annotations = annotations\n        self.go_graph = go_graph\n        self.ia_weights = ia_weights\n        self.include_dipeptides = include_dipeptides\n    \n    def extract_all_features(self, protein_id: str, sequence: str) -> Dict[str, float]:\n        \"\"\"Extract all features for a single protein.\"\"\"\n        features = {}\n        \n        # 1. Basic features (~35 features)\n        features.update(extract_basic_features(sequence))\n        \n        # 2. Physicochemical features (~18 features)\n        features.update(extract_physicochemical_features(sequence))\n        \n        # 3. Dipeptide features (400 features) - OPTIONAL\n        if self.include_dipeptides:\n            features.update(extract_kmer_features(sequence, k=2))\n        \n        # 4. Positional features (~50 features)\n        features.update(extract_positional_features(sequence))\n        \n        # 5. GO hierarchy features (~20 features) - CRITICAL\n        features.update(extract_go_features(\n            protein_id, self.annotations, self.go_graph, self.ia_weights))\n        \n        # 6. Pattern features (~15 features)\n        features.update(extract_pattern_features(sequence))\n        \n        return features\n    \n    def extract_batch(self, \n                     sequences: Dict[str, str],\n                     verbose: bool = True) -> pd.DataFrame:\n        \"\"\"Extract features for multiple proteins with progress tracking.\"\"\"\n        \n        if verbose:\n            print(f\"\\n{'='*80}\")\n            print(\" \"*20 + \"EXTRACTING FEATURES\")\n            print(f\"{'='*80}\")\n            print(f\"\\nProcessing {len(sequences):,} proteins...\")\n        \n        features_list = []\n        protein_ids = []\n        \n        for i, (protein_id, sequence) in enumerate(sequences.items()):\n            if verbose and (i + 1) % 5000 == 0:\n                print(f\"  Progress: {i+1:,} / {len(sequences):,} ({(i+1)/len(sequences)*100:.1f}%)\")\n            \n            try:\n                features = self.extract_all_features(protein_id, sequence)\n                features_list.append(features)\n                protein_ids.append(protein_id)\n            except Exception as e:\n                if verbose:\n                    print(f\"  Error processing {protein_id}: {e}\")\n                continue\n        \n        # Create DataFrame\n        df = pd.DataFrame(features_list, index=protein_ids)\n        \n        if verbose:\n            print(f\"\\n{'='*80}\")\n            print(\"✓ FEATURE EXTRACTION COMPLETE\")\n            print(f\"{'='*80}\")\n            print(f\"\\nResults:\")\n            print(f\"  - Total proteins: {len(df):,}\")\n            print(f\"  - Total features: {len(df.columns):,}\")\n            print(f\"\\nFeature breakdown:\")\n            self._print_feature_summary(df)\n        \n        return df\n    \n    def _print_feature_summary(self, df: pd.DataFrame):\n        \"\"\"Print detailed feature summary.\"\"\"\n        categories = {\n            'Basic (length, MW, composition)': sum(1 for c in df.columns \n                if any(x in c for x in ['length', 'molecular', 'comp_', 'group_'])),\n            'Physicochemical (hydro, charge, etc)': sum(1 for c in df.columns \n                if c.startswith('phys_')),\n            'K-mers (dipeptides)': sum(1 for c in df.columns \n                if c.startswith('kmer')),\n            'Positional (N/C-term, regions)': sum(1 for c in df.columns \n                if any(x in c for x in ['nterm_', 'cterm_', 'region'])),\n            'GO Hierarchy (CRITICAL for CAFA)': sum(1 for c in df.columns \n                if c.startswith('go_')),\n            'Patterns (motifs, repeats)': sum(1 for c in df.columns \n                if c.startswith('pattern_'))\n        }\n        \n        for category, count in categories.items():\n            print(f\"  - {category}: {count}\")\n\n# ============================================================================\n# 11. EXECUTE FEATURE EXTRACTION\n# ============================================================================\n\nprint(\"\\n\" + \"=\"*80)\nprint(\" \"*15 + \"STARTING FEATURE EXTRACTION PIPELINE\")\nprint(\"=\"*80)\n\n# Initialize extractor\nextractor = CAFAFeatureExtractor(\n    annotations=train_annotations,\n    go_graph=go_graph,\n    ia_weights=ia_weights,\n    include_dipeptides=Config.INCLUDE_DIPEPTIDES\n)\n\n# Extract features for training set\ntrain_features = extractor.extract_batch(train_sequences, verbose=True)\n\n# ============================================================================\n# 12. FEATURE VALIDATION & QUALITY CHECKS\n# ============================================================================\n\nprint(\"\\n\" + \"=\"*80)\nprint(\" \"*20 + \"FEATURE VALIDATION\")\nprint(\"=\"*80)\n\n# Check for missing values\nmissing_counts = train_features.isnull().sum()\nif missing_counts.sum() > 0:\n    print(f\"\\n⚠️  Found {missing_counts.sum()} missing values\")\n    print(f\"   Columns with missing data: {(missing_counts > 0).sum()}\")\n    # Fill missing values\n    train_features = train_features.fillna(0)\n    print(\"   ✓ Filled with zeros\")\nelse:\n    print(\"\\n✓ No missing values detected\")\n\n# Check for infinite values\ninf_mask = np.isinf(train_features.values)\nif inf_mask.any():\n    print(f\"\\n⚠️  Found {inf_mask.sum()} infinite values\")\n    train_features = train_features.replace([np.inf, -np.inf], 0)\n    print(\"   ✓ Replaced with zeros\")\nelse:\n    print(\"✓ No infinite values detected\")\n\n# Feature statistics\nprint(f\"\\n✓ Feature statistics:\")\nprint(f\"  - Mean feature variance: {train_features.var().mean():.4f}\")\nprint(f\"  - Features with zero variance: {(train_features.var() == 0).sum()}\")\nprint(f\"  - Memory usage: {train_features.memory_usage(deep=True).sum() / 1024**2:.2f} MB\")\n\n# ============================================================================\n# 13. FEATURE SELECTION & OPTIMIZATION\n# ============================================================================\n\nprint(\"\\n\" + \"=\"*80)\nprint(\" \"*18 + \"FEATURE OPTIMIZATION\")\nprint(\"=\"*80)\n\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.feature_selection import VarianceThreshold\n\n# Remove zero-variance features\nprint(\"\\n1️⃣  Removing zero-variance features...\")\ninitial_features = len(train_features.columns)\nselector = VarianceThreshold(threshold=0.0)\ntrain_features_selected = pd.DataFrame(\n    selector.fit_transform(train_features),\n    index=train_features.index,\n    columns=train_features.columns[selector.get_support()]\n)\nprint(f\"   Removed: {initial_features - len(train_features_selected.columns)}\")\nprint(f\"   Remaining: {len(train_features_selected.columns)}\")\n\n# Standardize features\nprint(\"\\n2️⃣  Standardizing features...\")\nscaler = StandardScaler()\ntrain_features_scaled = pd.DataFrame(\n    scaler.fit_transform(train_features_selected),\n    index=train_features_selected.index,\n    columns=train_features_selected.columns\n)\nprint(\"   ✓ Features standardized (mean=0, std=1)\")\n\n# ============================================================================\n# 14. PREPARE TARGET LABELS\n# ============================================================================\n\nprint(\"\\n\" + \"=\"*80)\nprint(\" \"*18 + \"PREPARING TARGET LABELS\")\nprint(\"=\"*80)\n\n# Create multi-label binary matrix for each ontology\ndef create_label_matrix(annotations: pd.DataFrame, \n                       protein_ids: List[str],\n                       ontology: str = None) -> pd.DataFrame:\n    \"\"\"Create binary label matrix for GO terms.\"\"\"\n    \n    # Filter by ontology if specified\n    if ontology:\n        annot_filtered = annotations[annotations['ontology'] == ontology]\n    else:\n        annot_filtered = annotations\n    \n    # Get unique GO terms\n    unique_terms = sorted(annot_filtered['go_term'].unique())\n    \n    # Create binary matrix\n    label_matrix = pd.DataFrame(0, index=protein_ids, columns=unique_terms)\n    \n    for _, row in annot_filtered.iterrows():\n        protein_id = row['protein_id']\n        go_term = row['go_term']\n        if protein_id in label_matrix.index:\n            label_matrix.loc[protein_id, go_term] = 1\n    \n    return label_matrix\n\n# Create label matrices for each ontology\nprotein_ids = train_features_scaled.index.tolist()\n\nprint(\"\\nCreating label matrices...\")\nlabels_mfo = create_label_matrix(train_annotations, protein_ids, 'MFO')\nlabels_bpo = create_label_matrix(train_annotations, protein_ids, 'BPO')\nlabels_cco = create_label_matrix(train_annotations, protein_ids, 'CCO')\n\nprint(f\"\\n✓ Label matrices created:\")\nprint(f\"  - MFO: {labels_mfo.shape[0]:,} proteins × {labels_mfo.shape[1]:,} terms\")\nprint(f\"  - BPO: {labels_bpo.shape[0]:,} proteins × {labels_bpo.shape[1]:,} terms\")\nprint(f\"  - CCO: {labels_cco.shape[0]:,} proteins × {labels_cco.shape[1]:,} terms\")\nprint(f\"\\n  Label statistics:\")\nprint(f\"  - MFO avg labels per protein: {labels_mfo.sum(axis=1).mean():.2f}\")\nprint(f\"  - BPO avg labels per protein: {labels_bpo.sum(axis=1).mean():.2f}\")\nprint(f\"  - CCO avg labels per protein: {labels_cco.sum(axis=1).mean():.2f}\")\n\n# ============================================================================\n# 15. SAVE PROCESSED DATA\n# ============================================================================\n\nprint(\"\\n\" + \"=\"*80)\nprint(\" \"*20 + \"SAVING PROCESSED DATA\")\nprint(\"=\"*80)\n\n# Save features\nprint(\"\\nSaving feature matrices...\")\ntrain_features_scaled.to_csv('train_features_scaled.csv')\nprint(\"  ✓ train_features_scaled.csv\")\n\n# Save labels\nlabels_mfo.to_csv('train_labels_mfo.csv')\nlabels_bpo.to_csv('train_labels_bpo.csv')\nlabels_cco.to_csv('train_labels_cco.csv')\nprint(\"  ✓ train_labels_mfo.csv\")\nprint(\"  ✓ train_labels_bpo.csv\")\nprint(\"  ✓ train_labels_cco.csv\")\n\n# Save feature names for later use\nfeature_info = pd.DataFrame({\n    'feature_name': train_features_scaled.columns,\n    'feature_index': range(len(train_features_scaled.columns))\n})\nfeature_info.to_csv('feature_info.csv', index=False)\nprint(\"  ✓ feature_info.csv\")\n\n# ============================================================================\n# 16. SUMMARY & NEXT STEPS\n# ============================================================================\n\nprint(\"\\n\" + \"=\"*80)\nprint(\" \"*15 + \"✓ FEATURE ENGINEERING COMPLETE!\")\nprint(\"=\"*80)\n\nprint(f\"\"\"\nSUMMARY:\n--------\n✓ Processed {len(train_features_scaled):,} proteins\n✓ Generated {len(train_features_scaled.columns)} features per protein\n✓ Created multi-label matrices for 3 ontologies\n✓ All data standardized and ready for modeling\n\nFEATURE CATEGORIES:\n-------------------\n- Basic sequence features (length, composition, MW)\n- Physicochemical properties (hydrophobicity, charge, volume)\n- K-mer features (dipeptide frequencies)\n- Positional features (N/C-terminal, regional properties)\n- GO hierarchy features (depth, IA weights, specificity) ⭐ CRITICAL\n- Sequence patterns (motifs, repeats, signals)\n\nOUTPUT FILES:\n-------------\n1. train_features_scaled.csv     - Standardized feature matrix\n2. train_labels_mfo.csv          - MFO labels\n3. train_labels_bpo.csv          - BPO labels  \n4. train_labels_cco.csv          - CCO labels\n5. feature_info.csv              - Feature metadata\n\nNEXT STEPS:\n-----------\n1. Train multi-label classifiers (one per ontology)\n2. Consider ensemble methods (Random Forest, XGBoost, Neural Networks)\n3. Implement hierarchical prediction (respect GO graph structure)\n4. Use Information Accretion weights in loss function\n5. Validate on hold-out set before submission\n6. Generate predictions for test superset\n\n\"\"\")\n\nprint(\"=\"*80)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-31T00:28:52.615770Z","iopub.execute_input":"2025-12-31T00:28:52.616383Z","iopub.status.idle":"2025-12-31T01:27:02.593018Z","shell.execute_reply.started":"2025-12-31T00:28:52.616356Z","shell.execute_reply":"2025-12-31T01:27:02.592314Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}