{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.10.0"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":116062,"databundleVersionId":14875579,"sourceType":"competition"}],"dockerImageVersionId":31236,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# 🧬 CAFA 6: Baseline with Regularization\n\n## 📋 Overview\n\nThis notebook demonstrates overfitting prevention techniques for protein function prediction.\n\n**Key Features:**\n- ✅ **All 3 GO aspects** (Biological Process, Molecular Function, Cellular Component)\n- ✅ **1500 GO terms** across all aspects\n- ✅ **Regularization techniques** (Early stopping, Dropout, Weight Decay)\n- ✅ **Adaptive learning rate** (ReduceLROnPlateau)\n- ✅ **100% official competition data** - no external dependencies\n- ✅ **Internet OFF compatible**\n\n**Technical approach:**\n1. Parse official FASTA and GO annotation files\n2. Extract TF-IDF k-mer features (3-grams)\n3. Train multi-label classifier with regularization\n4. Monitor validation loss with early stopping\n5. Propagate predictions through GO hierarchy\n6. Generate submission file","metadata":{}},{"cell_type":"code","source":"# No external packages needed!\nimport os\nimport gc\nimport numpy as np\nimport pandas as pd\nimport torch\nimport torch.nn as nn\nimport torch.optim as optim\nfrom torch.utils.data import Dataset, DataLoader, random_split\nfrom sklearn.feature_extraction.text import TfidfVectorizer\nfrom sklearn.preprocessing import MultiLabelBinarizer\nfrom collections import defaultdict, Counter\nfrom tqdm.auto import tqdm\n\n# Device\nDEVICE = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\nprint(f\"🚀 Device: {DEVICE}\")\nprint(f\"📊 PyTorch: {torch.__version__}\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Configuration\nclass CFG:\n    # Official data paths\n    TRAIN_TERMS = '/kaggle/input/cafa-6-protein-function-prediction/Train/train_terms.tsv'\n    TRAIN_SEQS = '/kaggle/input/cafa-6-protein-function-prediction/Train/train_sequences.fasta'\n    TRAIN_TAX = '/kaggle/input/cafa-6-protein-function-prediction/Train/train_taxonomy.tsv'\n    TEST_SEQS = '/kaggle/input/cafa-6-protein-function-prediction/Test/testsuperset.fasta'\n    OBO_FILE = '/kaggle/input/cafa-6-protein-function-prediction/Train/go-basic.obo'\n    \n    # Model parameters\n    BATCH_SIZE = 128\n    EPOCHS = 50             # Early stopping will likely stop earlier\n    LEARNING_RATE = 1e-3\n    HIDDEN_DIM = 256        # Appropriately sized for dataset\n    DROPOUT = 0.5           # Strong regularization\n    WEIGHT_DECAY = 1e-4     # L2 regularization\n    \n    # Early stopping & learning rate scheduling\n    EARLY_STOP_PATIENCE = 7  # Stop if no improvement for 7 epochs\n    LR_PATIENCE = 3          # Reduce LR if no improvement for 3 epochs\n    LR_FACTOR = 0.5          # Reduce LR by half\n    MIN_LR = 1e-6\n    \n    # Feature extraction\n    KMER_SIZE = 3           # 3-mer amino acid k-mers\n    MAX_FEATURES = 10000    # TF-IDF feature limit\n    \n    # Post-processing\n    THRESHOLD = 0.1         # Minimum confidence for predictions\n    MIN_PREDS = 25          # Minimum predictions per protein\n    MAX_PREDS = 150         # Maximum predictions per protein\n    \n    # Target labels - all 3 GO aspects\n    TOP_K_LABELS = 1500         # Total GO terms to predict\n    TOP_K_PER_ASPECT = 500      # 500 per aspect (P, F, C)\n    ASPECTS = ['P', 'F', 'C']   # Biological Process, Molecular Function, Cellular Component\n    \n    DEVICE = DEVICE\n\nprint(\"✅ Configuration loaded\")\nprint(f\"   Model: {CFG.HIDDEN_DIM} → {CFG.HIDDEN_DIM//2} hidden units\")\nprint(f\"   Regularization: Dropout {CFG.DROPOUT}, Weight Decay {CFG.WEIGHT_DECAY}\")\nprint(f\"   GO Terms: {CFG.TOP_K_LABELS} across {len(CFG.ASPECTS)} aspects\")\nprint(f\"   Early stopping patience: {CFG.EARLY_STOP_PATIENCE} epochs\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 📖 Part 1: Parse Official Data\n\nAll parsing from official competition files - no external dependencies.","metadata":{}},{"cell_type":"code","source":"# Manual FASTA parser (no BioPython needed)\ndef parse_fasta_manual(fasta_path):\n    \"\"\"Parse FASTA without external libraries\"\"\"\n    sequences = {}\n    current_id = None\n    current_seq = []\n    \n    print(f\"📖 Parsing {os.path.basename(fasta_path)}...\")\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                # Parse new ID\n                parts = line[1:].split('|')\n                current_id = parts[1] if len(parts) > 1 else line[1:].split()[0]\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# Load sequences\ntrain_seqs = parse_fasta_manual(CFG.TRAIN_SEQS)\ntest_seqs = parse_fasta_manual(CFG.TEST_SEQS)\n\nprint(f\"✅ Train: {len(train_seqs):,} sequences\")\nprint(f\"✅ Test: {len(test_seqs):,} sequences\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Load training annotations - all 3 GO aspects\nprint(\"📥 Loading training annotations...\")\ntrain_terms_df = pd.read_csv(CFG.TRAIN_TERMS, sep='\\t', header=None, names=['protein', 'term', 'aspect'])\n\n# Include all aspects (P, F, C) for comprehensive function prediction\nprint(f\"✅ Total annotations: {len(train_terms_df):,}\")\nprint(f\"   Aspect P (Biological Process): {(train_terms_df['aspect'] == 'P').sum():,}\")\nprint(f\"   Aspect F (Molecular Function): {(train_terms_df['aspect'] == 'F').sum():,}\")\nprint(f\"   Aspect C (Cellular Component): {(train_terms_df['aspect'] == 'C').sum():,}\")\n\n# Create protein -> terms mapping\nprotein_to_terms = train_terms_df.groupby('protein')['term'].apply(list).to_dict()\n\nprint(f\"✅ {len(protein_to_terms):,} proteins with annotations\")\nprint(f\"✅ {train_terms_df['term'].nunique():,} unique GO terms (all aspects)\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Select top-K most common GO terms - balanced across aspects\nprint(\"📊 Selecting top GO terms per aspect...\")\n\ntop_terms = set()\n\n# Select top 500 from each aspect for balanced representation\nfor aspect in CFG.ASPECTS:\n    aspect_df = train_terms_df[train_terms_df['aspect'] == aspect]\n    term_counts = aspect_df['term'].value_counts()\n    top_aspect_terms = set(term_counts.head(CFG.TOP_K_PER_ASPECT).index)\n    top_terms.update(top_aspect_terms)\n    print(f\"   {aspect}: {len(top_aspect_terms)} terms\")\n\nprint(f\"✅ Selected total {len(top_terms)} GO terms (across all aspects)\")\n\n# Filter proteins that have at least one selected term\nvalid_proteins = []\nfor pid in protein_to_terms:\n    if pid in train_seqs:\n        filtered_terms = [t for t in protein_to_terms[pid] if t in top_terms]\n        if filtered_terms:\n            protein_to_terms[pid] = filtered_terms\n            valid_proteins.append(pid)\n\nprint(f\"✅ {len(valid_proteins):,} proteins with sequences and annotations\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 🧬 Part 2: Extract K-mer Features\n\nConvert amino acid sequences to k-mer features using TF-IDF.\n\n**What are k-mers?**\n- Sliding window of k amino acids\n- Example: \"ACDEFG\" → 3-mers: \"ACD\", \"CDE\", \"DEF\", \"EFG\"\n- Captures local sequence patterns","metadata":{}},{"cell_type":"code","source":"# K-mer extraction\ndef extract_kmers(sequence, k=3):\n    \"\"\"Extract k-mers from amino acid sequence\"\"\"\n    if len(sequence) < k:\n        return sequence\n    kmers = [sequence[i:i+k] for i in range(len(sequence) - k + 1)]\n    return ' '.join(kmers)\n\n# Extract k-mers for training sequences\nprint(f\"🔧 Extracting {CFG.KMER_SIZE}-mers...\")\ntrain_texts = [extract_kmers(train_seqs[pid], CFG.KMER_SIZE) for pid in tqdm(valid_proteins, desc=\"Train k-mers\")]\n\n# TF-IDF vectorization\nprint(f\"📐 Computing TF-IDF features...\")\ntfidf = TfidfVectorizer(\n    max_features=CFG.MAX_FEATURES,\n    analyzer='word',\n    token_pattern=r'(?u)\\b\\w+\\b'\n)\n\nX_train_sparse = tfidf.fit_transform(train_texts)\nX_train = X_train_sparse.toarray().astype(np.float32)\n\nprint(f\"✅ Feature matrix: {X_train.shape}\")\nprint(f\"   {X_train.shape[0]} proteins x {X_train.shape[1]} features\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Create multi-label targets\nprint(\"🎯 Creating multi-label targets...\")\n\nmlb = MultiLabelBinarizer(classes=sorted(top_terms))\ny_labels = [protein_to_terms[pid] for pid in valid_proteins]\nY_train = mlb.fit_transform(y_labels).astype(np.float32)\n\nprint(f\"✅ Label matrix: {Y_train.shape}\")\nprint(f\"   {Y_train.shape[0]} proteins x {Y_train.shape[1]} GO terms\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 🤖 Part 3: Train Model","metadata":{}},{"cell_type":"code","source":"# PyTorch Dataset\nclass ProteinDataset(Dataset):\n    def __init__(self, X, y=None):\n        self.X = torch.FloatTensor(X)\n        self.y = torch.FloatTensor(y) if y is not None else None\n    \n    def __len__(self):\n        return len(self.X)\n    \n    def __getitem__(self, idx):\n        if self.y is not None:\n            return self.X[idx], self.y[idx]\n        return self.X[idx]\n\n# Multi-layer Perceptron with strong regularization\nclass SimpleMLP(nn.Module):\n    def __init__(self, input_dim, output_dim, hidden_dim=256, dropout=0.5):\n        super().__init__()\n        # Architecture: input → 256 → 128 → output\n        # Moderate size prevents overfitting while maintaining capacity\n        self.net = nn.Sequential(\n            nn.Linear(input_dim, hidden_dim),\n            nn.BatchNorm1d(hidden_dim),\n            nn.ReLU(),\n            nn.Dropout(dropout),\n            nn.Linear(hidden_dim, hidden_dim // 2),\n            nn.BatchNorm1d(hidden_dim // 2),\n            nn.ReLU(),\n            nn.Dropout(dropout),\n            nn.Linear(hidden_dim // 2, output_dim)\n        )\n    \n    def forward(self, x):\n        return self.net(x)\n\nprint(\"✅ Model classes defined\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Prepare data\ndataset = ProteinDataset(X_train, Y_train)\ntrain_size = int(0.9 * len(dataset))\nval_size = len(dataset) - train_size\n\ntrain_ds, val_ds = random_split(dataset, [train_size, val_size])\ntrain_loader = DataLoader(train_ds, batch_size=CFG.BATCH_SIZE, shuffle=True)\nval_loader = DataLoader(val_ds, batch_size=CFG.BATCH_SIZE)\n\n# Initialize model\ninput_dim = X_train.shape[1]\noutput_dim = Y_train.shape[1]\n\nmodel = SimpleMLP(input_dim, output_dim, CFG.HIDDEN_DIM, CFG.DROPOUT).to(CFG.DEVICE)\nprint(f\"\\n🤖 Model: {input_dim} → {CFG.HIDDEN_DIM} → {CFG.HIDDEN_DIM//2} → {output_dim}\")\nprint(f\"Parameters: {sum(p.numel() for p in model.parameters()):,}\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Training with early stopping and adaptive learning rate\noptimizer = optim.Adam(model.parameters(), lr=CFG.LEARNING_RATE, weight_decay=CFG.WEIGHT_DECAY)\ncriterion = nn.BCEWithLogitsLoss()\n\n# Learning rate scheduler - reduces LR when validation loss plateaus\nscheduler = optim.lr_scheduler.ReduceLROnPlateau(\n    optimizer, \n    mode='min',\n    patience=CFG.LR_PATIENCE,\n    factor=CFG.LR_FACTOR,\n    min_lr=CFG.MIN_LR\n)\n\nprint(f\"\\n🚀 Training for up to {CFG.EPOCHS} epochs (with early stopping)...\\n\")\n\n# Early stopping variables\nbest_val_loss = float('inf')\nepochs_no_improve = 0\nbest_model_state = None\n\nfor epoch in range(CFG.EPOCHS):\n    # Training phase\n    model.train()\n    train_loss = 0\n    for X_batch, y_batch in train_loader:\n        X_batch, y_batch = X_batch.to(CFG.DEVICE), y_batch.to(CFG.DEVICE)\n        \n        optimizer.zero_grad()\n        outputs = model(X_batch)\n        loss = criterion(outputs, y_batch)\n        loss.backward()\n        optimizer.step()\n        \n        train_loss += loss.item()\n    \n    train_loss /= len(train_loader)\n    \n    # Validation phase\n    model.eval()\n    val_loss = 0\n    with torch.no_grad():\n        for X_batch, y_batch in val_loader:\n            X_batch, y_batch = X_batch.to(CFG.DEVICE), y_batch.to(CFG.DEVICE)\n            outputs = model(X_batch)\n            val_loss += criterion(outputs, y_batch).item()\n    \n    val_loss /= len(val_loader)\n    \n    # Update learning rate based on validation performance\n    old_lr = optimizer.param_groups[0]['lr']\n    scheduler.step(val_loss)\n    new_lr = optimizer.param_groups[0]['lr']\n    \n    # Track best model\n    if val_loss < best_val_loss:\n        best_val_loss = val_loss\n        epochs_no_improve = 0\n        best_model_state = model.state_dict().copy()\n        checkpoint_marker = \" ✅ [Best]\"\n    else:\n        epochs_no_improve += 1\n        checkpoint_marker = \"\"\n    \n    # Show LR change\n    lr_marker = f\" [LR: {old_lr:.2e} → {new_lr:.2e}]\" if old_lr != new_lr else \"\"\n    \n    # Progress logging\n    if (epoch + 1) % 5 == 0 or epoch == 0 or checkpoint_marker or lr_marker:\n        print(f\"Epoch {epoch+1:2d}/{CFG.EPOCHS}: Train {train_loss:.4f} | Val {val_loss:.4f}{lr_marker}{checkpoint_marker}\")\n    \n    # Early stopping check\n    if epochs_no_improve >= CFG.EARLY_STOP_PATIENCE:\n        print(f\"\\n⏹️  Early stopping triggered at epoch {epoch+1}\")\n        print(f\"   No improvement for {CFG.EARLY_STOP_PATIENCE} consecutive epochs\")\n        break\n\n# Restore best model weights\nif best_model_state is not None:\n    model.load_state_dict(best_model_state)\n    print(f\"\\n✅ Training complete! Best validation loss: {best_val_loss:.4f}\")\nelse:\n    print(\"\\n✅ Training complete!\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 🔮 Part 4: Generate Predictions","metadata":{}},{"cell_type":"code","source":"# Extract test features\nprint(\"🧪 Processing test sequences...\")\ntest_pids = list(test_seqs.keys())\ntest_texts = [extract_kmers(test_seqs[pid], CFG.KMER_SIZE) for pid in tqdm(test_pids, desc=\"Test k-mers\")]\n\nX_test_sparse = tfidf.transform(test_texts)\nX_test = X_test_sparse.toarray().astype(np.float32)\n\nprint(f\"✅ Test features: {X_test.shape}\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Predict\ntest_ds = ProteinDataset(X_test)\ntest_loader = DataLoader(test_ds, batch_size=64)\n\nprint(\"🔮 Generating predictions...\")\nmodel.eval()\nall_preds = []\n\nwith torch.no_grad():\n    for X_batch in tqdm(test_loader, desc=\"Predicting\"):\n        X_batch = X_batch.to(CFG.DEVICE)\n        outputs = torch.sigmoid(model(X_batch))\n        all_preds.append(outputs.cpu().numpy())\n\npredictions = np.vstack(all_preds)\nprint(f\"✅ Predictions shape: {predictions.shape}\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 🌳 Part 5: GO Hierarchy Propagation","metadata":{}},{"cell_type":"code","source":"# Parse GO hierarchy\ndef parse_go_hierarchy(obo_file):\n    parents = defaultdict(set)\n    current_id = None\n    \n    try:\n        with open(obo_file, 'r') as f:\n            for line in f:\n                line = line.strip()\n                if line == '[Term]':\n                    current_id = None\n                elif line.startswith('id: '):\n                    current_id = line[4:]\n                elif line.startswith('is_a: ') and current_id:\n                    parent = line.split()[1]\n                    parents[current_id].add(parent)\n                elif line.startswith('relationship: part_of ') and current_id:\n                    parts = line.split()\n                    if len(parts) >= 3:\n                        parent = parts[2]\n                        parents[current_id].add(parent)\n        return dict(parents)\n    except:\n        return {}\n\ngo_parents = parse_go_hierarchy(CFG.OBO_FILE)\nprint(f\"✅ Loaded GO hierarchy: {len(go_parents):,} terms\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Get ancestors\ndef get_ancestors(term, parents_map, cache={}):\n    if term in cache:\n        return cache[term]\n    \n    ancestors = set()\n    for parent in parents_map.get(term, []):\n        ancestors.add(parent)\n        ancestors.update(get_ancestors(parent, parents_map, cache))\n    \n    cache[term] = ancestors\n    return ancestors\n\n# Create submission rows\nprint(\"📝 Creating submission rows...\")\nsubmission_rows = []\n\nfor i, pid in enumerate(tqdm(test_pids, desc=\"Processing\")):\n    scores = predictions[i]\n    \n    # Get top predictions\n    top_indices = np.argsort(scores)[-CFG.MAX_PREDS:]\n    top_indices = [idx for idx in top_indices if scores[idx] >= CFG.THRESHOLD]\n    \n    # Ensure minimum predictions\n    if len(top_indices) < CFG.MIN_PREDS:\n        top_indices = np.argsort(scores)[-CFG.MIN_PREDS:]\n    \n    # Collect terms and propagate\n    term_scores = {}\n    for idx in top_indices:\n        term = mlb.classes_[idx]\n        score = float(scores[idx])\n        term_scores[term] = max(term_scores.get(term, 0), score)\n        \n        # Propagate to ancestors\n        if go_parents:\n            for ancestor in get_ancestors(term, go_parents):\n                if ancestor in mlb.classes_:\n                    term_scores[ancestor] = max(term_scores.get(ancestor, 0), score)\n    \n    # Add to submission\n    for term, score in term_scores.items():\n        if score > 0:\n            submission_rows.append({'pid': pid, 'term': term, 'p': score})\n\nprint(f\"✅ Generated {len(submission_rows):,} predictions\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 📤 Part 6: Create Submission","metadata":{}},{"cell_type":"code","source":"# Create submission DataFrame\nsubmission_df = pd.DataFrame(submission_rows)\nsubmission_df = submission_df.sort_values(['pid', 'p'], ascending=[True, False])\n\nprint(f\"📊 Submission Statistics:\")\nprint(f\"   Predictions: {len(submission_df):,}\")\nprint(f\"   Proteins: {submission_df['pid'].nunique():,}\")\nprint(f\"   GO terms: {submission_df['term'].nunique():,}\")\nprint(f\"   Avg per protein: {len(submission_df)/submission_df['pid'].nunique():.1f}\")\n\n# Save (NO HEADER!)\nsubmission_df.to_csv('submission.tsv', sep='\\t', index=False, header=False)\nprint(f\"\\n✅ Submission saved: submission.tsv\")\nprint(f\"   Size: {os.path.getsize('submission.tsv')/1024/1024:.1f} MB\")\n\nsubmission_df.head(10)","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 🎯 Technical Notes\n\n### Regularization Techniques\n\n**Overfitting Prevention:**\n- Early stopping (patience=7) - monitors validation loss\n- Dropout (0.5) - prevents co-adaptation\n- Weight Decay (L2, 1e-4) - penalizes large weights\n- Batch Normalization - stabilizes training\n- Model size: 2.3M parameters\n\n**Adaptive Learning:**\n- ReduceLROnPlateau scheduler\n- Patience: 3 epochs\n- Reduction factor: 0.5\n- Minimum LR: 1e-6\n\n**Training:**\n- Train/validation split (90/10)\n- Batch size: 128\n- Best model restoration\n\n### GO Term Coverage\n\n**All 3 aspects:**\n- P (Biological Process) - cellular processes\n- F (Molecular Function) - molecular activities\n- C (Cellular Component) - subcellular locations\n\n**1500 total terms** (500 per aspect)\n\n**GO Hierarchy Propagation:**\n- Leverages GO DAG structure\n- Predictions propagate to ancestors\n\n### Possible Extensions\n\n**Features:**\n- Multiple k-mer sizes (2, 3, 4)\n- Amino acid composition\n- Sequence properties (length, hydrophobicity)\n\n**Model:**\n- Deeper networks\n- Attention mechanisms\n- Protein Language Models (ProtT5, ESM2)\n\n**Training:**\n- Cross-validation\n- Focal loss for imbalance\n- Ensemble methods\n\n**Post-processing:**\n- Aspect-specific thresholds\n- Weighted hierarchy propagation\n- Confidence calibration","metadata":{}}]}