{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"}},"cells":[{"cell_type":"markdown","source":"# Competition: cafa-6-protein-function-prediction\n\n**Generated by Alexandria Research Assistant**\n\n**Dataset:** cafa-6-protein-function-prediction\n\n**Task:** competition - kaggle-competition\n\n---\n\n⚠️ **Note:** This notebook contains Alexandria markers (lines starting with `# ⚠️ ALEXANDRIA MARKER`) at the top of each code cell. These markers enable the 'Sync from Kaggle' feature to track cell outputs. Please do not delete them.","metadata":{}},{"cell_type":"markdown","source":"## Setup & Imports\n\nInstall and import all necessary libraries for protein sequence analysis, data handling, and modeling.","metadata":{}},{"cell_type":"code","source":"# ⚠️ ALEXANDRIA MARKER - DO NOT DELETE (used for syncing outputs from Kaggle)\nprint(\"===ALEXANDRIA_CELL_1_START===\")\n\n# Setup & Imports for Protein Sequence Analysis\n\n# Install required packages (if not already available)\ntry:\n    import Bio\nexcept ImportError:\n    print(\"Installing biopython...\")\n    !pip install biopython\n\ntry:\n    import tqdm\nexcept ImportError:\n    print(\"Installing tqdm...\")\n    !pip install tqdm\n\n# Import libraries for sequence parsing, data handling, and modeling\nimport os\nimport random\n\nimport numpy as np\nimport pandas as pd\nfrom tqdm import tqdm\n\nfrom Bio import SeqIO\n\nimport torch\nimport torch.nn as nn\nimport torch.optim as optim\n\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import accuracy_score\n\n# Set random seeds for reproducibility\nSEED = 42\nrandom.seed(SEED)\nnp.random.seed(SEED)\ntorch.manual_seed(SEED)\ntorch.cuda.manual_seed_all(SEED)\nos.environ[\"PYTHONHASHSEED\"] = str(SEED)\n\n# Define dataset path\nDATA_PATH = \"/kaggle/input/cafa-6-protein-function-prediction\"\n\nprint(\"Setup complete. Libraries imported and random seeds set.\")\nprint(f\"Dataset path: {DATA_PATH}\")\nprint(\"Directory contents:\")\nprint(os.listdir(DATA_PATH))\nprint(\"Train directory contents:\")\nprint(os.listdir(os.path.join(DATA_PATH, \"Train\")))\nprint(\"Test directory contents:\")\nprint(os.listdir(os.path.join(DATA_PATH, \"Test\")))","outputs":[],"metadata":{},"execution_count":null},{"cell_type":"markdown","source":"## Data Loading\n\nLoad all relevant CAFA-6 files, including FASTA sequences, GO terms, and taxonomy files, into memory.","metadata":{}},{"cell_type":"code","source":"# ⚠️ ALEXANDRIA MARKER - DO NOT DELETE (used for syncing outputs from Kaggle)\nprint(\"===ALEXANDRIA_CELL_2_START===\")\n\n# ⚠️ ALEXANDRIA MARKER - DO NOT DELETE (used for syncing outputs from Kaggle)\nprint(\"===ALEXANDRIA_CELL_2_START===\")\n\n# Cell 2: Data Loading\n# Load all relevant CAFA-6 files including FASTA sequences, GO terms, and taxonomy files\n\nprint(\"\\n\" + \"=\"*80)\nprint(\"CELL 2: DATA LOADING\")\nprint(\"=\"*80)\n\n# Define file paths\ntrain_sequences_path = os.path.join(DATA_PATH, \"Train\", \"train_sequences.fasta\")\ntest_sequences_path = os.path.join(DATA_PATH, \"Test\", \"testsuperset.fasta\")\ntrain_terms_path = os.path.join(DATA_PATH, \"Train\", \"train_terms.tsv\")\ntrain_taxonomy_path = os.path.join(DATA_PATH, \"Train\", \"train_taxonomy.tsv\")\ngo_basic_path = os.path.join(DATA_PATH, \"Train\", \"go-basic.obo\")\ntest_taxon_path = os.path.join(DATA_PATH, \"Test\", \"testsuperset-taxon-list.tsv\")\nia_path = os.path.join(DATA_PATH, \"IA.tsv\")\nsample_submission_path = os.path.join(DATA_PATH, \"sample_submission.tsv\")\n\n# Verify all files exist\nprint(\"\\nVerifying file presence...\")\nrequired_files = {\n    \"train_sequences.fasta\": train_sequences_path,\n    \"testsuperset.fasta\": test_sequences_path,\n    \"train_terms.tsv\": train_terms_path,\n    \"train_taxonomy.tsv\": train_taxonomy_path,\n    \"go-basic.obo\": go_basic_path,\n    \"testsuperset-taxon-list.tsv\": test_taxon_path,\n    \"IA.tsv\": ia_path,\n    \"sample_submission.tsv\": sample_submission_path\n}\n\nmissing_files = []\nfor file_name, file_path in required_files.items():\n    if os.path.exists(file_path):\n        print(f\"✓ {file_name}\")\n    else:\n        print(f\"✗ {file_name} - NOT FOUND\")\n        missing_files.append(file_name)\n\nif missing_files:\n    print(f\"\\nWarning: {len(missing_files)} file(s) missing: {missing_files}\")\nelse:\n    print(\"\\n✓ All required files found!\")\n\n# Load FASTA sequences using Bio.SeqIO\nprint(\"\\n\" + \"-\"*80)\nprint(\"Loading FASTA sequences...\")\nprint(\"-\"*80)\n\ntry:\n    print(\"\\nParsing train_sequences.fasta...\")\n    train_sequences = {}\n    train_seq_count = 0\n    for record in SeqIO.parse(train_sequences_path, \"fasta\"):\n        train_sequences[record.id] = record\n        train_seq_count += 1\n    print(f\"✓ Loaded {train_seq_count} training sequences\")\n    \n    # Display sample training sequence info\n    if train_seq_count > 0:\n        sample_id = list(train_sequences.keys())[0]\n        sample_record = train_sequences[sample_id]\n        print(f\"  Sample sequence ID: {sample_id}\")\n        print(f\"  Sample sequence length: {len(sample_record.seq)}\")\n        print(f\"  Sample sequence (first 50 bp): {str(sample_record.seq)[:50]}\")\nexcept Exception as e:\n    print(f\"✗ Error loading train_sequences.fasta: {e}\")\n    train_sequences = {}\n\ntry:\n    print(\"\\nParsing testsuperset.fasta...\")\n    test_sequences = {}\n    test_seq_count = 0\n    for record in SeqIO.parse(test_sequences_path, \"fasta\"):\n        test_sequences[record.id] = record\n        test_seq_count += 1\n    print(f\"✓ Loaded {test_seq_count} test sequences\")\n    \n    # Display sample test sequence info\n    if test_seq_count > 0:\n        sample_id = list(test_sequences.keys())[0]\n        sample_record = test_sequences[sample_id]\n        print(f\"  Sample sequence ID: {sample_id}\")\n        print(f\"  Sample sequence length: {len(sample_record.seq)}\")\n        print(f\"  Sample sequence (first 50 bp): {str(sample_record.seq)[:50]}\")\nexcept Exception as e:\n    print(f\"✗ Error loading testsuperset.fasta: {e}\")\n    test_sequences = {}\n\n# Load TSV files using pandas\nprint(\"\\n\" + \"-\"*80)\nprint(\"Loading TSV files...\")\nprint(\"-\"*80)\n\ntry:\n    print(\"\\nReading train_terms.tsv...\")\n    train_terms = pd.read_csv(train_terms_path, sep=\"\\t\")\n    print(f\"✓ Loaded train_terms: shape {train_terms.shape}\")\n    print(f\"  Columns: {list(train_terms.columns)}\")\n    print(f\"  First few rows:\")\n    print(train_terms.head())\nexcept Exception as e:\n    print(f\"✗ Error loading train_terms.tsv: {e}\")\n    train_terms = pd.DataFrame()\n\ntry:\n    print(\"\\nReading train_taxonomy.tsv...\")\n    train_taxonomy = pd.read_csv(train_taxonomy_path, sep=\"\\t\")\n    print(f\"✓ Loaded train_taxonomy: shape {train_taxonomy.shape}\")\n    print(f\"  Columns: {list(train_taxonomy.columns)}\")\n    print(f\"  First few rows:\")\n    print(train_taxonomy.head())\nexcept Exception as e:\n    print(f\"✗ Error loading train_taxonomy.tsv: {e}\")\n    train_taxonomy = pd.DataFrame()\n\ntry:\n    print(\"\\nReading testsuperset","outputs":[],"metadata":{},"execution_count":null},{"cell_type":"markdown","source":"## Exploratory Data Analysis (EDA)\n\nExplore the distribution of sequence lengths, GO term frequencies, and taxonomic diversity in the training data.","metadata":{}},{"cell_type":"code","source":"# ⚠️ ALEXANDRIA MARKER - DO NOT DELETE (used for syncing outputs from Kaggle)\nprint(\"===ALEXANDRIA_CELL_3_START===\")\n\n# ⚠️ ALEXANDRIA MARKER - DO NOT DELETE (used for syncing outputs from Kaggle)\nprint(\"===ALEXANDRIA_CELL_3_START===\")\n\n# Exploratory Data Analysis (EDA) for CAFA-6 Training Data\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom collections import Counter, defaultdict\n\nprint(\"\\n\" + \"=\"*80)\nprint(\"CELL 3: EXPLORATORY DATA ANALYSIS (EDA)\")\nprint(\"=\"*80)\n\n# 1. Sequence Length Distribution\nprint(\"\\n[1] Analyzing sequence length distribution...\")\n\ntry:\n    seq_lengths = [len(record.seq) for record in train_sequences.values()]\n    print(f\"  Number of training sequences: {len(seq_lengths)}\")\n    print(f\"  Min length: {np.min(seq_lengths)}, Max length: {np.max(seq_lengths)}, Median: {np.median(seq_lengths)}\")\n    \n    plt.figure(figsize=(10,5))\n    sns.histplot(seq_lengths, bins=50, kde=True, color='skyblue')\n    plt.title(\"Training Sequence Length Distribution\")\n    plt.xlabel(\"Sequence Length (aa)\")\n    plt.ylabel(\"Count\")\n    plt.grid(True, axis='y', linestyle='--', alpha=0.5)\n    plt.show()\n    \n    # Log scale plot for long tail\n    plt.figure(figsize=(10,5))\n    sns.histplot(seq_lengths, bins=50, log_scale=(False, True), color='salmon')\n    plt.title(\"Training Sequence Length Distribution (Log Count)\")\n    plt.xlabel(\"Sequence Length (aa)\")\n    plt.ylabel(\"Log Count\")\n    plt.grid(True, axis='y', linestyle='--', alpha=0.5)\n    plt.show()\nexcept Exception as e:\n    print(f\"✗ Error analyzing sequence lengths: {e}\")\n\n# 2. GO Term Frequency Analysis\nprint(\"\\n[2] Analyzing GO term frequencies...\")\n\ntry:\n    # train_terms: columns = ['EntryID', 'GO_ID']\n    go_counts = train_terms['GO_ID'].value_counts()\n    print(f\"  Number of unique GO terms: {go_counts.shape}\")\n    print(\"  Top 10 most frequent GO terms:\")\n    print(go_counts.head(10))\n    \n    plt.figure(figsize=(12,6))\n    sns.barplot(x=go_counts.head(20).index, y=go_counts.head(20).values, palette='viridis')\n    plt.title(\"Top 20 Most Frequent GO Terms in Training Data\")\n    plt.xlabel(\"GO Term\")\n    plt.ylabel(\"Frequency\")\n    plt.xticks(rotation=75, ha='right')\n    plt.tight_layout()\n    plt.show()\nexcept Exception as e:\n    print(f\"✗ Error analyzing GO term frequencies: {e}\")\n\n# 2b. GO Term Hierarchy (optional: show parent terms if possible)\nprint(\"\\n[2b] Parsing GO hierarchy for top terms...\")\n\ndef parse_go_obo(obo_path, go_ids):\n    \"\"\"Parse OBO file for names and namespaces of selected GO IDs.\"\"\"\n    go_info = {}\n    try:\n        with open(obo_path, \"r\") as f:\n            current_id = None\n            current_name = None\n            current_ns = None\n            for line in f:\n                if line.startswith(\"[Term]\"):\n                    current_id = None\n                    current_name = None\n                    current_ns = None\n                elif line.startswith(\"id: GO:\"):\n                    current_id = line.strip().split(\"id: \")[1]\n                elif line.startswith(\"name: \"):\n                    current_name = line.strip().split(\"name: \")[1]\n                elif line.startswith(\"namespace: \"):\n                    current_ns = line.strip().split(\"namespace: \")[1]\n                elif line.strip() == \"\" and current_id and current_id in go_ids:\n                    go_info[current_id] = {\"name\": current_name, \"namespace\": current_ns}\n        return go_info\n    except Exception as e:\n        print(f\"✗ Error parsing OBO file: {e}\")\n        return {}\n\ntry:\n    top_go_ids = go_counts.head(10).index.tolist()\n    go_obo_info = parse_go_obo(go_basic_path, set(top_go_ids))\n    print(\"  Top 10 GO term details:\")\n    for go_id in top_go_ids:\n        info = go_obo_info.get(go_id, {})\n        print(f\"    {go_id}: {info.get('name','?')} [{info.get('namespace','?')}]\")\nexcept Exception as e:\n    print(f\"✗ Error extracting GO hierarchy: {e}\")\n\n# 3. Taxonomic Diversity\nprint(\"\\n[3] Summarizing taxonomic coverage...\")\n\ntry:\n    # train_taxonomy: columns = ['EntryID', 'Taxon']\n    taxon_counts = train_taxonomy['Taxon'].value_counts()\n    print(f\"  Number of unique taxa: {taxon_counts.shape}\")\n    print(\"  Top 10 most frequent taxa:\")\n    print(taxon_counts.head(10))\n    \n    plt.figure(figsize=(12,6))\n    sns.barplot(x=taxon_counts.head(20).index.astype(str), y=taxon_counts.head(20).values, palette='mako')\n    plt.title(\"Top 20 Most Frequent Taxa in Training Data\")\n    plt.xlabel(\"Taxon ID\")\n    plt.ylabel(\"Frequency\")\n    plt.xticks(rotation=75, ha='right')\n    plt.tight_layout()\n    plt.show()\nexcept Exception as e:\n    print(f\"✗ Error analyzing taxonomic diversity: {e}\")\n\n# 4. Missing or Anomalous Entries\nprint(\"\\n[4] Identifying missing or anomalous entries...\")\n\ntry:\n    # Check for missing sequence records referenced in train_terms or train_taxonomy\n    seq_ids = set(train_sequences.keys())\n    terms_ids = set(train_terms['EntryID'])\n    taxon_ids = set(train_taxonomy['EntryID'])\n    \n    missing_in_seq_terms = terms_ids - seq_ids\n    missing_in_seq_taxon = taxon_ids - seq_ids\n    missing_in_terms = seq_ids - terms_ids\n    missing_in_taxon = seq_ids - taxon_ids\n    \n    print(f\"  Entries in train_terms missing sequence: {len(missing_in_seq_terms)}\")\n    print(f\"  Entries in train_taxonomy missing sequence: {len(missing_in_seq_taxon)}\")\n    print(f\"  Sequences missing GO annotation: {len(missing_in_terms)}\")\n    print(f\"  Sequences missing taxonomy: {len(missing_in_taxon)}\")\n    \n    # Check for anomalous sequence lengths (e.g., <20 or >5000 aa)\n    short_seqs = [sid for sid, rec in train_sequences.items() if len(rec.seq) < 20]\n    long_seqs = [sid for sid, rec in train_sequences.items() if len(rec.seq) > 5000]\n    print(f\"  Sequences shorter than 20 aa: {len(short_seqs)}\")\n    print(f\"  Sequences longer than 5000 aa: {len(long_seqs)}\")\nexcept Exception as e:\n    print(f\"✗ Error checking","outputs":[],"metadata":{},"execution_count":null},{"cell_type":"markdown","source":"## Preprocessing\n\nClean and preprocess sequences and labels, preparing them for feature extraction and modeling.","metadata":{}},{"cell_type":"code","source":"# ⚠️ ALEXANDRIA MARKER - DO NOT DELETE (used for syncing outputs from Kaggle)\nprint(\"===ALEXANDRIA_CELL_4_START===\")\n\n# ⚠️ ALEXANDRIA MARKER - DO NOT DELETE (used for syncing outputs from Kaggle)\nprint(\"===ALEXANDRIA_CELL_4_START===\")\n\n# CELL 4: PREPROCESSING\n# Clean and preprocess sequences and labels for feature extraction and modeling\n\nprint(\"\\n\" + \"=\"*80)\nprint(\"CELL 4: PREPROCESSING\")\nprint(\"=\"*80)\n\n# Define valid amino acids (IUPAC standard, 20 canonical + 'X' for unknown)\nVALID_AMINO_ACIDS = set(\"ACDEFGHIKLMNPQRSTVWYX\")\n\n# 1. Clean and standardize training sequences\nprint(\"\\n[1] Cleaning and standardizing training sequences...\")\n\nclean_train_sequences = {}\ninvalid_seq_ids = []\nfor seq_id, record in tqdm(train_sequences.items(), desc=\"Processing train sequences\"):\n    seq_str = str(record.seq).upper().replace(\" \", \"\").replace(\"\\n\", \"\")\n    # Remove non-amino acid characters\n    seq_str_clean = \"\".join([aa if aa in VALID_AMINO_ACIDS else \"X\" for aa in seq_str])\n    # Check for missing/invalid sequence (e.g., empty or mostly invalid)\n    if len(seq_str_clean) < 20 or seq_str_clean.count(\"X\") > 0.5 * len(seq_str_clean):\n        invalid_seq_ids.append(seq_id)\n        continue\n    # Update record with cleaned sequence\n    record.seq = record.seq.__class__(seq_str_clean)\n    clean_train_sequences[seq_id] = record\n\nprint(f\"✓ Cleaned training sequences: {len(clean_train_sequences)} valid, {len(invalid_seq_ids)} removed (short/invalid)\")\n\n# 2. Clean and standardize test sequences\nprint(\"\\n[2] Cleaning and standardizing test sequences...\")\n\nclean_test_sequences = {}\ninvalid_test_seq_ids = []\nfor seq_id, record in tqdm(test_sequences.items(), desc=\"Processing test sequences\"):\n    seq_str = str(record.seq).upper().replace(\" \", \"\").replace(\"\\n\", \"\")\n    seq_str_clean = \"\".join([aa if aa in VALID_AMINO_ACIDS else \"X\" for aa in seq_str])\n    if len(seq_str_clean) < 20 or seq_str_clean.count(\"X\") > 0.5 * len(seq_str_clean):\n        invalid_test_seq_ids.append(seq_id)\n        continue\n    record.seq = record.seq.__class__(seq_str_clean)\n    clean_test_sequences[seq_id] = record\n\nprint(f\"✓ Cleaned test sequences: {len(clean_test_sequences)} valid, {len(invalid_test_seq_ids)} removed (short/invalid)\")\n\n# 3. Map protein IDs to GO terms (labels)\nprint(\"\\n[3] Mapping protein IDs to GO terms...\")\n\n# Remove GO annotations for invalid sequences\nvalid_train_ids = set(clean_train_sequences.keys())\nfiltered_train_terms = train_terms[train_terms['EntryID'].isin(valid_train_ids)].copy()\n\n# Optionally filter out rare GO terms (e.g., <5 occurrences)\ngo_term_counts = filtered_train_terms['GO_ID'].value_counts()\nrare_go_terms = go_term_counts[go_term_counts < 5].index.tolist()\nfiltered_train_terms = filtered_train_terms[~filtered_train_terms['GO_ID'].isin(rare_go_terms)]\n\nprint(f\"✓ Filtered GO terms: {filtered_train_terms['GO_ID'].nunique()} unique (rare terms removed)\")\nprint(f\"  Remaining train-term pairs: {filtered_train_terms.shape}\")\n\n# 4. Map protein IDs to taxonomy\nprint(\"\\n[4] Mapping protein IDs to taxonomy...\")\n\nfiltered_train_taxonomy = train_taxonomy[train_taxonomy['EntryID'].isin(valid_train_ids)].copy()\nprint(f\"✓ Filtered taxonomy: {filtered_train_taxonomy['Taxon'].nunique()} unique taxa\")\nprint(f\"  Remaining train-taxonomy pairs: {filtered_train_taxonomy.shape}\")\n\n# 5. Prepare final training set (intersection of valid sequence, GO, and taxonomy)\nprint(\"\\n[5] Preparing final training set...\")\n\nfinal_train_ids = set(filtered_train_terms['EntryID']) & set(filtered_train_taxonomy['EntryID']) & valid_train_ids\nfinal_train_sequences = {sid: clean_train_sequences[sid] for sid in final_train_ids}\nfinal_train_terms = filtered_train_terms[filtered_train_terms['EntryID'].isin(final_train_ids)].copy()\nfinal_train_taxonomy = filtered_train_taxonomy[filtered_train_taxonomy['EntryID'].isin(final_train_ids)].copy()\n\nprint(f\"✓ Final training set: {len(final_train_sequences)} sequences, {final_train_terms.shape} GO annotations, {final_train_taxonomy.shape} taxonomy entries\")\n\n# 6. Summary statistics\nprint(\"\\n[6] Preprocessing summary:\")\nprint(f\"  Training sequences: {len(train_sequences)} original, {len(final_train_sequences)} usable\")\nprint(f\"  Test sequences: {len(test_sequences)} original, {len(clean_test_sequences)} usable\")\nprint(f\"  Unique GO terms: {final_train_terms['GO_ID'].nunique()}\")\nprint(f\"  Unique taxa: {final_train_taxonomy['Taxon'].nunique()}\")\n\n# 7. Save cleaned/preprocessed data for downstream steps (optional, in-memory for now)\npreprocessed_data = {\n    \"train_sequences\": final_train_sequences,\n    \"train_terms\": final_train_terms,\n    \"train_taxonomy\": final_train_taxonomy,\n    \"test_sequences\": clean_test_sequences\n}\n\nprint(\"\\n✓ Preprocessing complete. Data ready for feature extraction and modeling.\")","outputs":[],"metadata":{},"execution_count":null},{"cell_type":"markdown","source":"## Feature Engineering\n\nExtract biologically relevant features from protein sequences for model input.","metadata":{}},{"cell_type":"code","source":"# ⚠️ ALEXANDRIA MARKER - DO NOT DELETE (used for syncing outputs from Kaggle)\nprint(\"===ALEXANDRIA_CELL_5_START===\")\n\nprint(\"===ALEXANDRIA_CELL_5_START===\")\n\n# CELL 5: FEATURE ENGINEERING\n# Extract biologically relevant features from protein sequences for model input\n\nimport itertools\nfrom collections import Counter, defaultdict\n\nprint(\"\\n\" + \"=\"*80)\nprint(\"CELL 5: FEATURE ENGINEERING\")\nprint(\"=\"*80)\n\n# 1. K-mer (n-gram) frequency extraction for protein sequences\ndef compute_kmer_frequencies(seq, k=3, alphabet=\"ACDEFGHIKLMNPQRSTVWYX\"):\n    \"\"\"Compute normalized k-mer frequencies for a protein sequence.\"\"\"\n    seq = str(seq)\n    total_kmers = max(len(seq) - k + 1, 1)\n    kmer_counts = Counter([seq[i:i+k] for i in range(len(seq) - k + 1)])\n    # Generate all possible kmers for the alphabet\n    all_kmers = [''.join(p) for p in itertools.product(alphabet, repeat=k)]\n    freq_vec = np.array([kmer_counts[kmer]/total_kmers if total_kmers > 0 else 0.0 for kmer in all_kmers], dtype=np.float32)\n    return freq_vec\n\nprint(\"\\n[1] Computing k-mer frequencies for training and test sequences (k=3)...\")\nk = 3\nalphabet = \"ACDEFGHIKLMNPQRSTVWYX\"\nall_kmers = [''.join(p) for p in itertools.product(alphabet, repeat=k)]\nkmer_dim = len(all_kmers)\nprint(f\"  Using {k}-mers: {kmer_dim} possible combinations\")\n\ntrain_kmer_features = {}\nfor seq_id, record in tqdm(final_train_sequences.items(), desc=\"Train k-mer features\"):\n    train_kmer_features[seq_id] = compute_kmer_frequencies(record.seq, k=k, alphabet=alphabet)\nprint(f\"✓ Computed k-mer features for {len(train_kmer_features)} training sequences\")\n\ntest_kmer_features = {}\nfor seq_id, record in tqdm(clean_test_sequences.items(), desc=\"Test k-mer features\"):\n    test_kmer_features[seq_id] = compute_kmer_frequencies(record.seq, k=k, alphabet=alphabet)\nprint(f\"✓ Computed k-mer features for {len(test_kmer_features)} test sequences\")\n\n# 2. Sequence embeddings (ESM-1b) using torch.hub (if available)\nprint(\"\\n[2] Generating sequence embeddings using pretrained ESM-1b model (if available)...\")\ntry:\n    import esm\n    esm_model, alphabet_obj = esm.pretrained.esm1b_t33_650M_UR50S()\n    batch_converter = alphabet_obj.get_batch_converter()\n    esm_model.eval()\n    print(\"✓ ESM-1b model loaded\")\n    \n    def get_esm_embedding(seq):\n        \"\"\"Get mean-pooled ESM-1b embedding for a single sequence.\"\"\"\n        batch_labels, batch_strs, batch_tokens = batch_converter([(\"protein\", seq)])\n        with torch.no_grad():\n            results = esm_model(batch_tokens, repr_layers=, return_contacts=False)\n        token_representations = results[\"representations\"]\n        # Exclude start/end tokens ([CLS], [EOS])\n        emb = token_representations[0, 1:len(seq)+1].mean(0).cpu().numpy()\n        return emb.astype(np.float32)\n    \n    train_esm_embeddings = {}\n    for seq_id, record in tqdm(final_train_sequences.items(), desc=\"Train ESM embeddings\"):\n        try:\n            train_esm_embeddings[seq_id] = get_esm_embedding(str(record.seq))\n        except Exception as e:\n            print(f\"  [Train] {seq_id}: {e}\")\n    print(f\"✓ Generated ESM embeddings for {len(train_esm_embeddings)} training sequences\")\n    \n    test_esm_embeddings = {}\n    for seq_id, record in tqdm(clean_test_sequences.items(), desc=\"Test ESM embeddings\"):\n        try:\n            test_esm_embeddings[seq_id] = get_esm_embedding(str(record.seq))\n        except Exception as e:\n            print(f\"  [Test] {seq_id}: {e}\")\n    print(f\"✓ Generated ESM embeddings for {len(test_esm_embeddings)} test sequences\")\nexcept Exception as e:\n    print(f\"✗ ESM embedding generation failed: {e}\")\n    train_esm_embeddings = {}\n    test_esm_embeddings = {}\n\n# 3. Incorporate taxonomy/domain features (one-hot encoding for taxonomy)\nprint(\"\\n[3] Encoding taxonomy features (one-hot)...\")\ntry:\n    unique_taxa = sorted(final_train_taxonomy['Taxon'].unique())\n    taxon2idx = {taxon: idx for idx, taxon in enumerate(unique_taxa)}\n    taxon_dim = len(taxon2idx)\n    print(f\"  {taxon_dim} unique taxa\")\n\n    train_taxon_features = {}\n    for _, row in final_train_taxonomy.iterrows():\n        vec = np.zeros(taxon_dim, dtype=np.float32)\n        idx = taxon2idx[row['Taxon']]\n        vec[idx] = 1.0\n        train_taxon_features[row['EntryID']] = vec\n    print(f\"✓ Encoded taxonomy features for {len(train_taxon_features)} training sequences\")\nexcept Exception as e:\n    print(f\"✗ Taxonomy encoding failed: {e}\")\n    train_taxon_features = {}\n\n# 4. Prepare multi-hot encoded GO term targets\nprint(\"\\n[4] Preparing multi-hot encoded GO term targets...\")\ntry:\n    unique_go_terms = sorted(final_train_terms['GO_ID'].unique())\n    go2idx = {go: idx for idx, go in enumerate(unique_go_terms)}\n    go_dim = len(go2idx)\n    print(f\"  {go_dim} unique GO terms\")\n\n    # Map EntryID to set of GO term indices\n    entryid2goidx = defaultdict(set)\n    for _, row in final_train_terms.iterrows():\n        entryid2goidx[row['EntryID']].add(go2idx[row['GO_ID']])\n    \n    train_go_targets = {}\n    for seq_id in final_train_sequences.keys():\n        vec = np.zeros(go_dim, dtype=np.float32)\n        for idx in entryid2goidx.get(seq_id, []):\n            vec[idx] = 1.0\n        train_go_targets[seq_id] = vec\n    print(f\"✓ Multi-hot GO term targets prepared for {len(train_go_targets)} training sequences\")\nexcept Exception as e:\n    print(f\"✗ GO term multi-hot encoding failed: {e}\")\n    train_go_targets = {}\n\n# 5. Assemble final feature matrices for model input\nprint(\"\\n[5] Assembling feature matrices for model input...\")\n\ndef stack_features(feature_dict, ids, dim):\n    X = np.zeros((len(ids), dim), dtype=np.float32)\n    for i, sid in enumerate(ids):\n        if sid in feature_dict:\n            X[i] = feature_dict[sid]\n    return X\n\ntrain_ids = list(final_train_sequences.keys())\nX_kmer = stack_features(train_kme","outputs":[],"metadata":{},"execution_count":null},{"cell_type":"markdown","source":"## Model Training\n\nTrain a multi-label classifier (e.g., deep neural network) to predict GO terms from sequence features.","metadata":{}},{"cell_type":"code","source":"# ⚠️ ALEXANDRIA MARKER - DO NOT DELETE (used for syncing outputs from Kaggle)\nprint(\"===ALEXANDRIA_CELL_6_START===\")\n\nprint(\"===ALEXANDRIA_CELL_6_START===\")\n\n# CELL 6: MODEL TRAINING\n# Train a multi-label classifier to predict GO terms from sequence features\n\nimport torch\nimport torch.nn as nn\nimport torch.optim as optim\nfrom torch.utils.data import Dataset, DataLoader\nfrom sklearn.model_selection import train_test_split\n\nprint(\"\\n\" + \"=\"*80)\nprint(\"CELL 6: MODEL TRAINING\")\nprint(\"=\"*80)\n\n# 1. Prepare feature matrices and targets\nprint(\"\\n[1] Preparing feature matrices and targets for modeling...\")\n\ntry:\n    # Use k-mer features, taxonomy features, and optionally ESM embeddings if available\n    use_esm = len(train_esm_embeddings) == len(train_ids)\n    feature_list = []\n    feature_dims = []\n    if use_esm:\n        print(\"  Using ESM embeddings + taxonomy features\")\n        X_esm = stack_features(train_esm_embeddings, train_ids, len(next(iter(train_esm_embeddings.values()))))\n        feature_list.append(X_esm)\n        feature_dims.append(X_esm.shape[1])\n    else:\n        print(\"  Using k-mer features + taxonomy features\")\n        X_kmer = stack_features(train_kmer_features, train_ids, kmer_dim)\n        feature_list.append(X_kmer)\n        feature_dims.append(X_kmer.shape[1])\n    X_taxon = stack_features(train_taxon_features, train_ids, taxon_dim)\n    feature_list.append(X_taxon)\n    feature_dims.append(X_taxon.shape[1])\n    X = np.concatenate(feature_list, axis=1)\n    y = stack_features(train_go_targets, train_ids, go_dim)\n    print(f\"  Feature matrix shape: {X.shape}\")\n    print(f\"  Target matrix shape: {y.shape}\")\nexcept Exception as e:\n    print(f\"✗ Error preparing feature matrices: {e}\")\n    raise\n\n# 2. Train/validation split\nprint(\"\\n[2] Splitting data into training and validation sets...\")\ntry:\n    X_train, X_val, y_train, y_val = train_test_split(\n        X, y, test_size=0.15, random_state=SEED, stratify=None\n    )\n    print(f\"  Training set: {X_train.shape} samples\")\n    print(f\"  Validation set: {X_val.shape} samples\")\nexcept Exception as e:\n    print(f\"✗ Error splitting data: {e}\")\n    raise\n\n# 3. Define PyTorch Dataset\nclass ProteinDataset(Dataset):\n    def __init__(self, X, y):\n        self.X = torch.tensor(X, dtype=torch.float32)\n        self.y = torch.tensor(y, dtype=torch.float32)\n    def __len__(self):\n        return self.X.shape\n    def __getitem__(self, idx):\n        return self.X[idx], self.y[idx]\n\ntrain_dataset = ProteinDataset(X_train, y_train)\nval_dataset = ProteinDataset(X_val, y_val)\n\n# 4. Model architecture (Feed-forward DNN)\nprint(\"\\n[3] Defining model architecture...\")\n\nclass MultiLabelClassifier(nn.Module):\n    def __init__(self, input_dim, output_dim, hidden_dims=[512, 256, 128], dropout=0.3):\n        super().__init__()\n        layers = []\n        dims = [input_dim] + hidden_dims\n        for i in range(len(hidden_dims)):\n            layers.append(nn.Linear(dims[i], dims[i+1]))\n            layers.append(nn.BatchNorm1d(dims[i+1]))\n            layers.append(nn.ReLU())\n            layers.append(nn.Dropout(dropout))\n        layers.append(nn.Linear(hidden_dims[-1], output_dim))\n        self.net = nn.Sequential(*layers)\n    def forward(self, x):\n        return torch.sigmoid(self.net(x))\n\ninput_dim = X.shape[1]\noutput_dim = y.shape[1]\nmodel = MultiLabelClassifier(input_dim, output_dim)\nprint(f\"  Model input dim: {input_dim}, output dim: {output_dim}\")\n\n# 5. Training setup\nbatch_size = 128\nnum_epochs = 50\npatience = 7\nmin_delta = 0.001\nlr = 1e-3\ndevice = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\nmodel = model.to(device)\n\ntrain_loader = DataLoader(train_dataset, batch_size=batch_size, shuffle=True, num_workers=2)\nval_loader = DataLoader(val_dataset, batch_size=batch_size, shuffle=False, num_workers=2)\n\ncriterion = nn.BCELoss()\noptimizer = optim.Adam(model.parameters(), lr=lr)\nscheduler = optim.lr_scheduler.ReduceLROnPlateau(optimizer, mode='max', factor=0.5, patience=3, verbose=True)\n\n# 6. Early stopping and checkpointing\nbest_val_score = -np.inf\nepochs_no_improve = 0\ncheckpoint_path = \"best_model.pt\"\n\ndef compute_f1(y_true, y_pred, threshold=0.5):\n    y_pred_bin = (y_pred > threshold).astype(np.float32)\n    tp = (y_true * y_pred_bin).sum(axis=1)\n    fp = ((1 - y_true) * y_pred_bin).sum(axis=1)\n    fn = (y_true * (1 - y_pred_bin)).sum(axis=1)\n    f1 = 2 * tp / (2 * tp + fp + fn + 1e-8)\n    return np.mean(f1)\n\nprint(\"\\n[4] Starting model training...\")\nfor epoch in range(num_epochs):\n    model.train()\n    train_loss = 0.0\n    for Xb, yb in train_loader:\n        Xb, yb = Xb.to(device), yb.to(device)\n        optimizer.zero_grad()\n        outputs = model(Xb)\n        loss = criterion(outputs, yb)\n        loss.backward()\n        optimizer.step()\n        train_loss += loss.item() * Xb.size(0)\n    train_loss /= len(train_loader.dataset)\n\n    # Validation\n    model.eval()\n    val_loss = 0.0\n    val_preds = []\n    val_targets = []\n    with torch.no_grad():\n        for Xb, yb in val_loader:\n            Xb, yb = Xb.to(device), yb.to(device)\n            outputs = model(Xb)\n            loss = criterion(outputs, yb)\n            val_loss += loss.item() * Xb.size(0)\n            val_preds.append(outputs.cpu().numpy())\n            val_targets.append(yb.cpu().numpy())\n    val_loss /= len(val_loader.dataset)\n    val_preds = np.concatenate(val_preds, axis=0)\n    val_targets = np.concatenate(val_targets, axis=0)\n    val_f1 = compute_f1(val_targets, val_preds, threshold=0.5)\n\n    scheduler.step(val_f1)\n\n    print(f\"Epoch {epoch+1:02d}/{num_epochs} | Train Loss: {train_loss:.4f} | Val Loss: {val_loss:.4f} | Val F1: {val_f1:.4f}\")\n\n    # Early stopping\n    if val_f1 > best_val_score + min_delta:","outputs":[],"metadata":{},"execution_count":null},{"cell_type":"markdown","source":"## Evaluation\n\nEvaluate model performance using CAFA-relevant metrics and visualize results.","metadata":{}},{"cell_type":"code","source":"# ⚠️ ALEXANDRIA MARKER - DO NOT DELETE (used for syncing outputs from Kaggle)\nprint(\"===ALEXANDRIA_CELL_7_START===\")\n\n# ⚠️ ALEXANDRIA MARKER - DO NOT DELETE (used for syncing outputs from Kaggle)\nprint(\"===ALEXANDRIA_CELL_7_START===\")\n\n# CELL 7: EVALUATION\n# Evaluate model performance using CAFA-relevant metrics and visualize results\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom sklearn.metrics import (\n    precision_score, recall_score, f1_score, roc_auc_score,\n    hamming_loss, accuracy_score, precision_recall_curve, roc_curve, auc\n)\nimport warnings\nwarnings.filterwarnings('ignore')\n\nprint(\"\\n\" + \"=\"*80)\nprint(\"CELL 7: EVALUATION\")\nprint(\"=\"*80)\n\n# 1. Load best model checkpoint and generate predictions\nprint(\"\\n[1] Loading best model and generating predictions...\")\ntry:\n    model.load_state_dict(torch.load(checkpoint_path, map_location=device))\n    model.eval()\n    print(\"✓ Best model loaded from checkpoint\")\nexcept Exception as e:\n    print(f\"⚠ Could not load checkpoint: {e}. Using current model.\")\n\n# Generate predictions on validation set\nprint(\"\\n[2] Generating predictions on validation set...\")\nval_preds_all = []\nval_targets_all = []\nwith torch.no_grad():\n    for Xb, yb in val_loader:\n        Xb, yb = Xb.to(device), yb.to(device)\n        outputs = model(Xb)\n        val_preds_all.append(outputs.cpu().numpy())\n        val_targets_all.append(yb.cpu().numpy())\n\nval_preds_proba = np.concatenate(val_preds_all, axis=0)\nval_targets = np.concatenate(val_targets_all, axis=0)\nprint(f\"✓ Predictions shape: {val_preds_proba.shape}\")\n\n# 3. Compute CAFA-relevant metrics\nprint(\"\\n[3] Computing CAFA-relevant evaluation metrics...\")\n\n# Threshold for binary predictions\nthreshold = 0.5\nval_preds_binary = (val_preds_proba > threshold).astype(np.float32)\n\n# Example-based metrics\nprint(\"\\n  [Example-Based Metrics]\")\ntry:\n    # Exact Match Ratio (EMR) / Subset Accuracy\n    emr = accuracy_score(val_targets, val_preds_binary)\n    print(f\"    Exact Match Ratio (EMR): {emr:.4f}\")\nexcept Exception as e:\n    print(f\"    EMR computation failed: {e}\")\n    emr = None\n\ntry:\n    # 0/1 Loss\n    zero_one_loss = 1.0 - emr if emr is not None else None\n    print(f\"    0/1 Loss: {zero_one_loss:.4f}\")\nexcept Exception as e:\n    print(f\"    0/1 Loss computation failed: {e}\")\n\ntry:\n    # Hamming Loss\n    h_loss = hamming_loss(val_targets, val_preds_binary)\n    print(f\"    Hamming Loss: {h_loss:.4f}\")\nexcept Exception as e:\n    print(f\"    Hamming Loss computation failed: {e}\")\n    h_loss = None\n\n# Label-based metrics (macro and micro averaging)\nprint(\"\\n  [Label-Based Metrics]\")\ntry:\n    # Macro-averaged precision, recall, F1\n    macro_precision = precision_score(val_targets, val_preds_binary, average='macro', zero_division=0)\n    macro_recall = recall_score(val_targets, val_preds_binary, average='macro', zero_division=0)\n    macro_f1 = f1_score(val_targets, val_preds_binary, average='macro', zero_division=0)\n    print(f\"    Macro Precision: {macro_precision:.4f}\")\n    print(f\"    Macro Recall: {macro_recall:.4f}\")\n    print(f\"    Macro F1: {macro_f1:.4f}\")\nexcept Exception as e:\n    print(f\"    Macro metrics computation failed: {e}\")\n\ntry:\n    # Micro-averaged precision, recall, F1\n    micro_precision = precision_score(val_targets, val_preds_binary, average='micro', zero_division=0)\n    micro_recall = recall_score(val_targets, val_preds_binary, average='micro', zero_division=0)\n    micro_f1 = f1_score(val_targets, val_preds_binary, average='micro', zero_division=0)\n    print(f\"    Micro Precision: {micro_precision:.4f}\")\n    print(f\"    Micro Recall: {micro_recall:.4f}\")\n    print(f\"    Micro F1: {micro_f1:.4f}\")\nexcept Exception as e:\n    print(f\"    Micro metrics computation failed: {e}\")\n\ntry:\n    # Weighted-averaged F1\n    weighted_f1 = f1_score(val_targets, val_preds_binary, average='weighted', zero_division=0)\n    print(f\"    Weighted F1: {weighted_f1:.4f}\")\nexcept Exception as e:\n    print(f\"    Weighted F1 computation failed: {e}\")\n\n# 4. Compute Fmax (F-measure at optimal threshold)\nprint(\"\\n[4] Computing Fmax (F-measure at optimal threshold)...\")\ntry:\n    fmax = 0.0\n    fmax_threshold = 0.5\n    thresholds_to_test = np.arange(0.1, 1.0, 0.05)\n    \n    for t in thresholds_to_test:\n        preds_t = (val_preds_proba > t).astype(np.float32)\n        f1_t = f1_score(val_targets, preds_t, average='micro', zero_division=0)\n        if f1_t > fmax","outputs":[],"metadata":{},"execution_count":null},{"cell_type":"markdown","source":"## Submission/Results\n\nFormat predictions for Kaggle submission and display key results.","metadata":{}},{"cell_type":"code","source":"# ⚠️ ALEXANDRIA MARKER - DO NOT DELETE (used for syncing outputs from Kaggle)\nprint(\"===ALEXANDRIA_CELL_8_START===\")\n\nprint(\"===ALEXANDRIA_CELL_8_START===\")\n\n# CELL 8: SUBMISSION/RESULTS\n# Format predictions for Kaggle submission and display key results\n\nimport os\nimport pandas as pd\nimport numpy as np\n\nprint(\"\\n\" + \"=\"*80)\nprint(\"CELL 8: SUBMISSION/RESULTS\")\nprint(\"=\"*80)\n\n# 1. Load sample submission and test protein IDs\nsample_submission_path = os.path.join(DATA_PATH, \"sample_submission.tsv\")\ntry:\n    sample_submission = pd.read_csv(sample_submission_path, sep=\"\\t\", header=None)\n    print(f\"✓ Loaded sample_submission.tsv: shape {sample_submission.shape}\")\nexcept Exception as e:\n    print(f\"✗ Error loading sample_submission.tsv: {e}\")\n    sample_submission = pd.DataFrame()\n\n# Columns: Protein (EntryID), GO (GO_ID), Prediction (float)\nif sample_submission.shape[1] == 3:\n    sample_submission.columns = [\"Protein\", \"GO\", \"Prediction\"]\n\n# 2. Prepare test protein IDs and GO terms for submission\ntest_protein_ids = list(clean_test_sequences.keys())\nsubmission_go_terms = sorted(go2idx.keys())  # Use trained GO terms only\n\nprint(f\"  Number of test proteins: {len(test_protein_ids)}\")\nprint(f\"  Number of GO terms for prediction: {len(submission_go_terms)}\")\n\n# 3. Generate predictions for test set proteins\nprint(\"\\n[1] Generating predictions for test set proteins...\")\n\n# Prepare test features\ntry:\n    use_esm_test = len(test_esm_embeddings) == len(test_protein_ids)\n    test_feature_list = []\n    if use_esm_test:\n        print(\"  Using ESM embeddings + taxonomy features for test set\")\n        X_test_esm = stack_features(test_esm_embeddings, test_protein_ids, len(next(iter(test_esm_embeddings.values()))))\n        test_feature_list.append(X_test_esm)\n    else:\n        print(\"  Using k-mer features + taxonomy features for test set\")\n        X_test_kmer = stack_features(test_kmer_features, test_protein_ids, kmer_dim)\n        test_feature_list.append(X_test_kmer)\n    # For taxonomy, fill zeros for missing test taxa\n    X_test_taxon = np.zeros((len(test_protein_ids), taxon_dim), dtype=np.float32)\n    test_feature_list.append(X_test_taxon)\n    X_test = np.concatenate(test_feature_list, axis=1)\n    print(f\"  Test feature matrix shape: {X_test.shape}\")\nexcept Exception as e:\n    print(f\"✗ Error preparing test features: {e}\")\n    X_test = np.zeros((len(test_protein_ids), input_dim), dtype=np.float32)\n\n# Generate predictions\nmodel.eval()\nwith torch.no_grad():\n    X_test_tensor = torch.tensor(X_test, dtype=torch.float32).to(device)\n    test_preds_proba = model(X_test_tensor).cpu().numpy()\nprint(\"✓ Test predictions generated\")\n\n# 4. Format output to match sample_submission.tsv requirements\nprint(\"\\n[2] Formatting predictions for Kaggle submission...\")\n\n# Threshold predictions (clip to [0,1], optionally apply cutoff)\ncutoff_threshold_low = 0.1\ntest_preds_proba = np.clip(test_preds_proba, 0, 1)\n\nsubmission_rows = []\nfor i, protein_id in enumerate(test_protein_ids):\n    for j, go_term in enumerate(submission_go_terms):\n        pred_score = test_preds_proba[i, j]\n        if pred_score >= cutoff_threshold_low:\n            submission_rows.append([protein_id, go_term, round(float(pred_score), 6)])\n\nsubmission_df = pd.DataFrame(submission_rows, columns=[\"Protein\", \"GO\", \"Prediction\"])\nprint(f\"✓ Submission dataframe created: {submission_df.shape} rows\")\n\n# 5. Save submission file\nsubmission_file = \"submission.tsv\"\ntry:\n    submission_df.to_csv(submission_file, sep=\"\\t\", header=False, index=False)\n    print(f\"✓ Submission file saved: {submission_file}\")\nexcept Exception as e:\n    print(f\"✗ Error saving submission file: {e}\")\n\n# 6. Summarize key findings and next steps\nprint(\"\\n[3] Key Results & Next Steps\")\nprint(\"-\"*80)\ntry:\n    print(f\"  Validation Macro F1: {macro_f1:.4f}\")\n    print(f\"  Validation Micro F1: {micro_f1:.4f}\")\n    print(f\"  Validation Weighted F1: {weighted_f1:.4f}\")\n    print(f\"  Validation Exact Match Ratio (EMR): {emr:.4f}\")\n    print(f\"  Validation Hamming Loss: {h_loss:.4f}\")\nexcept Exception as e:\n    print(f\"  Could not display validation metrics: {e}\")\n\nprint(\"\\n  Next Steps:\")\nprint(\"  - Experiment with alternative sequence embeddings (e.g., ProtT5, ProtBERT)\")\nprint(\"  - Incorporate ontology graph constraints in postprocessing\")\nprint(\"  - Tune model architecture and hyperparameters\")\nprint(\"  - Analyze error cases and improve feature engineering\")\n\nprint(\"\\n✓ Submission and results summary complete. Ready for Kaggle upload.\")","outputs":[],"metadata":{},"execution_count":null}],"nbformat":4,"nbformat_minor":4}