{"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":"nvidiaTeslaT4","dataSources":[{"sourceId":116062,"databundleVersionId":14084779,"sourceType":"competition"}],"dockerImageVersionId":31154,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"markdown","source":"# **GoID Fit & Predict: RandomForest**","metadata":{}},{"cell_type":"code","source":"!pip install biopython","metadata":{"trusted":true,"_kg_hide-output":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom Bio import SeqIO\nimport os\nimport gc","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-18T04:12:46.279591Z","iopub.execute_input":"2025-10-18T04:12:46.279971Z","iopub.status.idle":"2025-10-18T04:12:46.284885Z","shell.execute_reply.started":"2025-10-18T04:12:46.279934Z","shell.execute_reply":"2025-10-18T04:12:46.283956Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Process Train Data","metadata":{}},{"cell_type":"code","source":"file_path = '/kaggle/input/cafa-6-protein-function-prediction/Train/train_sequences.fasta'\n\nif not os.path.exists(file_path):\n    print(f\"Error: File not found at {file_path}\")\nelse:\n    sequences = list(SeqIO.parse(file_path, 'fasta'))\n    n_sequences = len(sequences)\n    print(f\"✅ Total number of sequences loaded: {n_sequences}\")\n\n    lengths = [len(record.seq) for record in sequences]\n\n    df_lengths = pd.DataFrame(lengths, columns=['Length'])\n\nif 'sequences' not in locals() or not sequences:\n    print(\"Error: The 'sequences' list is not found or is empty. Please ensure the FASTA file loading was successful.\")\nelse:\n    # 1. Data Extraction\n    data = []\n    for record in sequences:\n        # Split the EntryID by the '|' (pipe) delimiter\n        parts = record.id.split('|')\n\n        if len(parts) == 3:\n            accession_type = parts[0]\n            accession_id = parts[1]\n            protein_name = parts[2]\n        else:\n            accession_type = 'N/A'\n            accession_id = record.id\n            protein_name = 'N/A'\n\n        data.append({\n            'Full_EntryID': record.id,\n            'Length': len(record.seq),\n            'Accession_Type': accession_type,\n            'Accession_ID': accession_id,\n            'Protein_Name': protein_name,\n            'Sequence': str(record.seq)[0:32] # The sequence itself is usually omitted due to its length\n        })\n\n    # 2. DataFrame Creation\n    df = pd.DataFrame(data)\n    display(df.head())\n    df.to_csv('fasta.csv',index=False)\n    \n    fasta1=df[['Accession_ID','Sequence']]\n    fasta1.columns=['EntryID','Sequence']\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"    fasta0=pd.read_csv('/kaggle/input/fasta-protein-data-length-and-sequence/fasta.csv')\n    fasta1=fasta0[['Accession_ID','Sequence']]\n    fasta1.columns=['EntryID','Sequence']\n    display(fasta1[0:3])","metadata":{"execution":{"iopub.status.busy":"2025-10-18T04:12:46.286045Z","iopub.execute_input":"2025-10-18T04:12:46.286362Z","iopub.status.idle":"2025-10-18T04:12:46.573582Z","shell.execute_reply.started":"2025-10-18T04:12:46.286338Z","shell.execute_reply":"2025-10-18T04:12:46.572801Z"}}},{"cell_type":"code","source":"df0 = pd.read_csv('/kaggle/input/cafa-6-protein-function-prediction/Train/train_terms.tsv', sep='\\t')\nprint(len(df0))\ndf1=df0[['EntryID','term']]\ndf1.columns=['EntryID','GoID']\ndf2=df1.merge(fasta1,on='EntryID',how='left')\ndisplay(df2)\n\ndf=df2[['EntryID','Sequence','GoID']][0:40000] #max=300000,#n=540000\ndf['value'] = 1\ndf_pivot = df.pivot_table(index=['EntryID','Sequence'], columns='GoID', values='value', fill_value=0)\ndisplay(df_pivot)\ngc.collect()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-18T04:12:46.575399Z","iopub.execute_input":"2025-10-18T04:12:46.575636Z","iopub.status.idle":"2025-10-18T04:12:46.987939Z","shell.execute_reply.started":"2025-10-18T04:12:46.575617Z","shell.execute_reply":"2025-10-18T04:12:46.987046Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Fit","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nfrom sklearn.feature_extraction.text import CountVectorizer\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.multioutput import MultiOutputClassifier\nfrom sklearn.ensemble import RandomForestClassifier\nfrom sklearn.metrics import jaccard_score, hamming_loss, classification_report\n\n# 1. Define the k-mer function (k=3)\ndef get_kmers(sequence, k=3):\n    \"\"\"Splits a sequence into overlapping k-mers.\"\"\"\n    return [sequence[i:i + k] for i in range(len(sequence) - k + 1)]\n\n# 2. Extract Sequences (X_text) and Target Labels (Y)\n# Use the Sequence index from the pivot table for alignment\nsequences = df_pivot.index.get_level_values('Sequence').astype(str)\nY = df_pivot.values  # Target matrix (0s and 1s for each GoID)\n\n# 3. Apply k-mer tokenization\n# We join the k-mers back into a string for CountVectorizer\nk = 3\nX_kmers = [' '.join(get_kmers(seq, k=k)) for seq in sequences]\ndel sequences  \ngc.collect()\n\n# 4. Use CountVectorizer to create the feature matrix X\nvectorizer = CountVectorizer(analyzer='word', ngram_range=(1, 1))\nX = vectorizer.fit_transform(X_kmers)\n\nprint(f\"✅ Feature Matrix (X) Shape: {X.shape}\")\nprint(f\"✅ Target Matrix (Y) Shape: {Y.shape}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-18T04:12:47.244629Z","iopub.execute_input":"2025-10-18T04:12:47.244889Z","iopub.status.idle":"2025-10-18T04:12:49.499242Z","shell.execute_reply.started":"2025-10-18T04:12:47.24487Z","shell.execute_reply":"2025-10-18T04:12:49.498355Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 5. Split the data (0.9:0.1 as requested)\nX_train, X_test, Y_train, Y_test = train_test_split(\n    X, Y,\n    test_size=0.1,      # 10% for testing\n    random_state=42     # for reproducibility\n)\ndel X, Y\ngc.collect()\nprint(f\"Data Split: Train samples={X_train.shape[0]}, Test samples={X_test.shape[0]}\")\nprint(\"-\" * 50)\n\nfrom sklearn.multioutput import MultiOutputClassifier\nfrom sklearn.ensemble import RandomForestClassifier\nfrom sklearn.metrics import jaccard_score, hamming_loss, classification_report\nimport numpy as np\n\n# --- 6. Initialize and Train the Multi-Output Classifier ---\nprint(\"## 6. Model Training (Binary Relevance / MultiOutputClassifier)\")\n\n# Base estimator: Random Forest Classifier\nrf_estimator = RandomForestClassifier(\n    n_estimators=40,      # Use 100 decision trees\n    random_state=42,\n    n_jobs=-1              # Utilize all CPU cores for parallel processing\n)\n\n# Multi-label model wrapper (Binary Relevance Strategy)\nmulti_rf_model = MultiOutputClassifier(rf_estimator)\n\nprint(\"Starting model training... (Training a separate RF for each GO Term)\")\nmulti_rf_model.fit(X_train, Y_train)\nprint(\"Training complete! 🎉\")\nprint(\"-\" * 50)\ngc.collect()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-18T04:12:49.500117Z","iopub.execute_input":"2025-10-18T04:12:49.500522Z","execution_failed":"2025-10-18T04:19:47.582Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# --- 7. Prediction (Predicting Probabilities) - SAFELY MODIFIED ---\nprint(\"## 7. Predicting Probabilities & Defining Y_pred_thresholded\")\n\nY_proba_list = []\n\n#for model in multi_rf_model.estimators_:\nfor idx, model in enumerate(multi_rf_model.estimators_):\n    # Predict probabilities for the current GO Term\n    proba_output = model.predict_proba(X_test)\n    \n    # Check the shape for single-class models\n    if proba_output.shape[1] == 2:\n        # Standard case: Two columns available. Extract P(class=1).\n        Y_proba_list.append(proba_output[:, 1])\n    else:\n        # Single-class case: Assume P(class=1) is 0.0 for safety.\n        # This handles the IndexError and provides a valid probability column.\n        n_samples = proba_output.shape[0]\n        Y_proba_list.append(np.zeros(n_samples))\n    \n    if (idx + 1) % 1000 == 0:\n        gc.collect()\n\n\n# Combine the results into the final (n_samples, n_labels) NumPy array\nimport numpy as np\nY_proba = np.stack(Y_proba_list, axis=1)\nprint(f\"Y_proba (Probabilities) successfully generated. Shape: {Y_proba.shape}\")\n\n# Convert probabilities to binary predictions using a 0.5 threshold\n# THIS LINE DEFINES THE MISSING VARIABLE: Y_pred_thresholded\nY_pred_thresholded = (Y_proba >= 0.5).astype(int)\nprint(\"Binary predictions (Y_pred_thresholded) created using Threshold = 0.5.\")\nprint(\"-\" * 50)\n\n\n# --- 8. Evaluation ---\nprint(\"## 8. Evaluation Metrics (Using Threshold 0.5)\")\nfrom sklearn.metrics import jaccard_score, hamming_loss, classification_report\n\n# Evaluate using the thresholded binary predictions\njaccard = jaccard_score(Y_test, Y_pred_thresholded, average='samples')\nhamming = hamming_loss(Y_test, Y_pred_thresholded)\n\nprint(f\"Multi-Label Classification Results (Threshold 0.5):\")\nprint(f\"-> Jaccard Index (Similarity): {jaccard:.4f}\")\nprint(f\"-> Hamming Loss (Error Rate): {hamming:.4f}\")\n\n# Example: Get a detailed report for the first GO term (column index 0)\nprint(\"\\nClassification Report for the first GO Term:\")\nprint(classification_report(Y_test[:, 0], Y_pred_thresholded[:, 0], zero_division=0))","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"markdown","source":"# Process Test Data and Predict","metadata":{}},{"cell_type":"code","source":"# ========== DISK-OPTIMIZED: Sparse predictions + Compression ==========\n\nimport numpy as np\nimport pandas as pd\nfrom Bio import SeqIO\nimport gc\nimport gzip\n\n# Load GO terms data\ngo_terms_raw = pd.read_csv('/kaggle/input/cafa-6-protein-function-prediction/Train/train_terms.tsv', sep='\\t')\n\n# Extract UNIQUE GO IDs from the pivot table columns (these are the trained models)\n# Use the GO IDs that were actually trained on\nunique_go_ids = df_pivot.columns.values  # These are the 10,045 GO terms with trained models\nprint(f\"Total trained GO Terms: {len(unique_go_ids)}\")\nprint(f\"Number of trained models: {len(multi_rf_model.estimators_)}\")\n\n# Verify they match\nassert len(unique_go_ids) == len(multi_rf_model.estimators_), \"Mismatch between GO terms and models!\"\n\ntsequences = list(SeqIO.parse('/kaggle/input/cafa-6-protein-function-prediction/Test/testsuperset.fasta', 'fasta'))\nprint(f'n_test_sequences: {len(tsequences)}')\n\n# Create DataFrame\ndata = {\n    'EntryID': [record.id for record in tsequences],\n    'Sequence': [str(record.seq)[0:320] for record in tsequences]\n}\n\ndf_TEST = pd.DataFrame(data)\nprint(f\"Test DataFrame shape: {df_TEST.shape}\")\nprint(\"-\" * 50)\n\n# ===== Apply preprocessing =====\nprint(\"Applying k-mer vectorization to test data...\")\n\ndef get_kmers(sequence, k=3):\n    return [sequence[i:i + k] for i in range(len(sequence) - k + 1)]\n\nk = 3\nX_kmers_test = [' '.join(get_kmers(seq, k=k)) for seq in df_TEST['Sequence'].astype(str)]\n\nX_TEST_vectorized = vectorizer.transform(X_kmers_test)\ndel X_kmers_test\ngc.collect()\n\nprint(f\"✅ Test Feature Matrix Shape: {X_TEST_vectorized.shape}\")\nprint(\"-\" * 50)\n\n# ===== DISK-EFFICIENT: Process + write TSV, filtering low scores =====\nprint(\"Processing predictions and writing to TSV file...\")\n\nsubmission_file_path = 'submission.tsv'\nbatch_size = 500\nn_go_terms = len(unique_go_ids)\nscore_threshold = 0.01  # Only write scores > threshold (reduces file size)\n\nwith open(submission_file_path, 'w') as f:\n    for batch_start in range(0, n_go_terms, batch_size):\n        batch_end = min(batch_start + batch_size, n_go_terms)\n        batch_size_actual = batch_end - batch_start\n        \n        batch_indices = range(batch_start, batch_end)\n        \n        # Process this batch of GO Terms\n        batch_predictions = []\n        \n        for idx in batch_indices:\n            model = multi_rf_model.estimators_[idx]\n            proba_output = model.predict_proba(X_TEST_vectorized)\n            \n            if proba_output.shape[1] == 2:\n                proba = proba_output[:, 1]\n            else:\n                proba = np.zeros(proba_output.shape[0])\n            \n            batch_predictions.append(proba)\n        \n        # Stack batch predictions\n        Y_batch = np.stack(batch_predictions, axis=1)\n        del batch_predictions\n        \n        # Get GO IDs for this batch\n        go_ids_batch = unique_go_ids[batch_start:batch_end]\n        \n        # Write results - ONLY scores above threshold (sparse output)\n        for i, entry_id in enumerate(df_TEST['EntryID'].values):\n            for j, go_id in enumerate(go_ids_batch):\n                score = Y_batch[i, j]\n                # Only write if score exceeds threshold\n                if score >= score_threshold:\n                    f.write(f\"{entry_id}\\t{go_id}\\t{score:.6f}\\n\")\n        \n        del Y_batch\n        gc.collect()\n        \n        print(f\"  Processed GO Terms {batch_end} / {n_go_terms}\")\n\nprint(f\"✅ Submission file saved: {submission_file_path}\")\nprint(f\"   Predictions written with score >= {score_threshold}\")\nprint(f\"   File format: Tab-separated values (EntryID, GO_ID, Score)\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}