{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":116062,"databundleVersionId":14084779,"sourceType":"competition"}],"dockerImageVersionId":31153,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nfrom pathlib import Path\nimport re\nfrom collections import defaultdict, Counter\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.preprocessing import MultiLabelBinarizer\nfrom sklearn.metrics import f1_score\nimport warnings\nwarnings.filterwarnings('ignore')\n\n\n# CONFIGURATION\n\nDATA_PATH = Path('/kaggle/input/cafa-6-protein-function-prediction')\nMAX_SEQ_LENGTH = 1000  # Truncate long sequences\nKMER_SIZE = 3  # For k-mer features\nTOP_N_TERMS = 1000  # Limit to most frequent GO terms for faster training\n\n\n# 1. DATA LOADING\n\nprint(\"Loading data...\")\n\ndef parse_fasta(fasta_file):\n    \"\"\"Parse FASTA file and return dictionary of sequences\"\"\"\n    sequences = {}\n    current_id = None\n    current_seq = []\n    \n    with open(fasta_file, 'r') as f:\n        for line in f:\n            line = line.strip()\n            if line.startswith('>'):\n                if current_id:\n                    sequences[current_id] = ''.join(current_seq)\n                # Extract protein ID from header (e.g., sp|P9WHI7|RECN_MYCT)\n                parts = line[1:].split('|')\n                if len(parts) >= 2:\n                    current_id = parts[1]\n                else:\n                    current_id = line[1:].split()[0]\n                current_seq = []\n            else:\n                current_seq.append(line)\n        \n        if current_id:\n            sequences[current_id] = ''.join(current_seq)\n    \n    return sequences\n\n# Load sequences\ntrain_sequences = parse_fasta(DATA_PATH / 'Train/train_sequences.fasta')\ntest_sequences = parse_fasta(DATA_PATH / 'Test/testsuperset.fasta')\n\nprint(f\"Loaded {len(train_sequences)} training sequences\")\nprint(f\"Loaded {len(test_sequences)} test sequences\")\n\n# Load annotations\ntrain_terms_df = pd.read_csv(DATA_PATH / 'Train/train_terms.tsv', sep='\\t', \n                             names=['protein_id', 'go_term', 'ontology'])\n\n# Load IA weights\nia_df = pd.read_csv(DATA_PATH / 'IA.tsv', sep='\\t', \n                    names=['go_term', 'ia_weight'])\nia_weights = dict(zip(ia_df['go_term'], ia_df['ia_weight']))\n\nprint(f\"Loaded {len(train_terms_df)} annotations\")\nprint(f\"Unique GO terms: {train_terms_df['go_term'].nunique()}\")\n\n\n# 2. FEATURE ENGINEERING\n\nprint(\"\\nGenerating features...\")\n\ndef get_amino_acid_composition(seq):\n    \"\"\"Calculate amino acid composition\"\"\"\n    amino_acids = 'ACDEFGHIKLMNPQRSTVWY'\n    total = len(seq)\n    if total == 0:\n        return np.zeros(20)\n    \n    composition = np.array([seq.count(aa) / total for aa in amino_acids])\n    return composition\n\ndef get_dipeptide_composition(seq):\n    \"\"\"Calculate dipeptide composition\"\"\"\n    amino_acids = 'ACDEFGHIKLMNPQRSTVWY'\n    dipeptides = [aa1 + aa2 for aa1 in amino_acids for aa2 in amino_acids]\n    total = len(seq) - 1\n    if total <= 0:\n        return np.zeros(400)\n    \n    composition = np.array([seq.count(dp) / total for dp in dipeptides])\n    return composition\n\ndef get_kmer_features(seq, k=3, max_features=100):\n    \"\"\"Get k-mer frequency features (limited for efficiency)\"\"\"\n    if len(seq) < k:\n        return np.zeros(max_features)\n    \n    kmers = [seq[i:i+k] for i in range(len(seq) - k + 1)]\n    kmer_counts = Counter(kmers)\n    \n    # Get top k-mers\n    top_kmers = [kmer for kmer, _ in kmer_counts.most_common(max_features)]\n    features = np.array([kmer_counts.get(kmer, 0) / len(kmers) for kmer in top_kmers])\n    \n    # Pad if needed\n    if len(features) < max_features:\n        features = np.pad(features, (0, max_features - len(features)))\n    \n    return features\n\ndef get_sequence_features(seq):\n    \"\"\"Extract comprehensive sequence features\"\"\"\n    seq = seq[:MAX_SEQ_LENGTH]  # Truncate if too long\n    \n    features = []\n    \n    # Basic properties\n    features.append(len(seq))\n    features.append(seq.count('C') / len(seq) if len(seq) > 0 else 0)  # Cysteine content\n    features.append(seq.count('M') / len(seq) if len(seq) > 0 else 0)  # Methionine content\n    \n    # Charged residues\n    positive = sum(seq.count(aa) for aa in 'RHK')\n    negative = sum(seq.count(aa) for aa in 'DE')\n    features.append(positive / len(seq) if len(seq) > 0 else 0)\n    features.append(negative / len(seq) if len(seq) > 0 else 0)\n    \n    # Hydrophobicity\n    hydrophobic = sum(seq.count(aa) for aa in 'AILMFWYV')\n    features.append(hydrophobic / len(seq) if len(seq) > 0 else 0)\n    \n    # Amino acid composition\n    aa_comp = get_amino_acid_composition(seq)\n    features.extend(aa_comp)\n    \n    # Dipeptide composition (sample subset for efficiency)\n    dipep_comp = get_dipeptide_composition(seq)\n    features.extend(dipep_comp[:100])  # Use first 100 dipeptides\n    \n    return np.array(features)\n\n# Generate features for training set\nprint(\"Extracting training features...\")\nX_train_list = []\ny_train_proteins = []\n\nfor protein_id in train_sequences.keys():\n    if protein_id in train_terms_df['protein_id'].values:\n        X_train_list.append(get_sequence_features(train_sequences[protein_id]))\n        y_train_proteins.append(protein_id)\n\nX_train = np.array(X_train_list)\nprint(f\"Training feature matrix shape: {X_train.shape}\")\n\n\n# 3. LABEL PREPARATION\n\nprint(\"\\nPreparing labels...\")\n\n# Get most frequent GO terms to limit complexity\nterm_counts = train_terms_df['go_term'].value_counts()\ntop_terms = term_counts.head(TOP_N_TERMS).index.tolist()\n\nprint(f\"Using top {len(top_terms)} GO terms\")\n\n# Create label matrix\nprotein_terms = defaultdict(list)\nfor _, row in train_terms_df.iterrows():\n    if row['go_term'] in top_terms:\n        protein_terms[row['protein_id']].append(row['go_term'])\n\ny_train_labels = [protein_terms[pid] for pid in y_train_proteins]\n\n# Binarize labels\nmlb = MultiLabelBinarizer(classes=top_terms)\ny_train = mlb.fit_transform(y_train_labels)\n\nprint(f\"Label matrix shape: {y_train.shape}\")\nprint(f\"Average labels per protein: {y_train.sum(axis=1).mean():.2f}\")\n\n\n# 4. MODEL TRAINING\n\nprint(\"\\nTraining model...\")\n\nfrom sklearn.multioutput import MultiOutputClassifier\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.ensemble import RandomForestClassifier\n\n# Use Logistic Regression for efficiency (faster than RF on large datasets)\n# For better performance, consider using XGBoost or neural networks\nbase_model = LogisticRegression(max_iter=100, random_state=42, n_jobs=-1)\nmodel = MultiOutputClassifier(base_model, n_jobs=-1)\n\n# Train on subset for faster execution (remove this for full training)\nsample_size = min(10000, len(X_train))\nindices = np.random.choice(len(X_train), sample_size, replace=False)\nX_train_sample = X_train[indices]\ny_train_sample = y_train[indices]\n\nprint(f\"Training on {sample_size} samples...\")\nmodel.fit(X_train_sample, y_train_sample)\n\nprint(\"Training complete!\")\n\n\n# 5. GENERATE PREDICTIONS\n\nprint(\"\\nGenerating predictions for test set...\")\n\n# Process test sequences in batches\nbatch_size = 1000\npredictions = []\n\ntest_protein_ids = list(test_sequences.keys())\nn_batches = (len(test_protein_ids) + batch_size - 1) // batch_size\n\nfor batch_idx in range(n_batches):\n    start_idx = batch_idx * batch_size\n    end_idx = min((batch_idx + 1) * batch_size, len(test_protein_ids))\n    batch_ids = test_protein_ids[start_idx:end_idx]\n    \n    # Extract features\n    X_test_batch = np.array([\n        get_sequence_features(test_sequences[pid]) \n        for pid in batch_ids\n    ])\n    \n    # Predict probabilities\n    y_pred_proba = model.predict_proba(X_test_batch)\n    \n    # Convert to predictions\n    for i, protein_id in enumerate(batch_ids):\n        for j, go_term in enumerate(mlb.classes_):\n            # Get probability (handle different sklearn versions)\n            if isinstance(y_pred_proba[j], np.ndarray):\n                prob = y_pred_proba[j][i, 1] if y_pred_proba[j].shape[1] > 1 else y_pred_proba[j][i]\n            else:\n                prob = y_pred_proba[j]\n            \n            # Only include predictions above threshold\n            if prob > 0.01:\n                predictions.append({\n                    'protein_id': protein_id,\n                    'go_term': go_term,\n                    'probability': min(prob, 0.999)  # Cap at 0.999\n                })\n    \n    if (batch_idx + 1) % 10 == 0:\n        print(f\"Processed {end_idx}/{len(test_protein_ids)} proteins...\")\n\n# 6. CREATE SUBMISSION FILE\n\nprint(\"\\nCreating submission file...\")\n\nsubmission_df = pd.DataFrame(predictions)\n\n# Sort by protein_id and probability\nsubmission_df = submission_df.sort_values(['protein_id', 'probability'], \n                                          ascending=[True, False])\n\n# Limit to 1500 terms per protein\nsubmission_df = submission_df.groupby('protein_id').head(1500)\n\n# Format probabilities to 3 significant figures\nsubmission_df['probability'] = submission_df['probability'].apply(\n    lambda x: f\"{x:.3g}\"\n)\n\n# Save submission\nsubmission_df.to_csv('submission.tsv', sep='\\t', header=False, index=False)\n\nprint(f\"\\nSubmission file created with {len(submission_df)} predictions\")\nprint(f\"Unique proteins: {submission_df['protein_id'].nunique()}\")\nprint(f\"Average predictions per protein: {len(submission_df) / submission_df['protein_id'].nunique():.1f}\")\n\n# Show sample\nprint(\"\\nSample predictions:\")\nprint(submission_df.head(10).to_string(index=False))\n\nprint(\"\\n✓ Complete! Submission file saved as 'submission.tsv'\")","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-11-04T06:06:30.938603Z","iopub.execute_input":"2025-11-04T06:06:30.939603Z","iopub.status.idle":"2025-11-04T06:32:12.177973Z","shell.execute_reply.started":"2025-11-04T06:06:30.939569Z","shell.execute_reply":"2025-11-04T06:32:12.176831Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"read_sub = pd.read_csv('/kaggle/working/submission.tsv')\nread_sub.head(10)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-04T06:36:18.041582Z","iopub.execute_input":"2025-11-04T06:36:18.042033Z","iopub.status.idle":"2025-11-04T06:36:27.785451Z","shell.execute_reply.started":"2025-11-04T06:36:18.042004Z","shell.execute_reply":"2025-11-04T06:36:27.784664Z"}},"outputs":[],"execution_count":null}]}