{"cells": [{"cell_type": "markdown", "metadata": {}, "source": "# CAFA-6 Diamond BLASTP Solution v1\n\n**Objetivo:** Usar Diamond BLASTP para encontrar proteinas similares y transferir GO terms\n\n**Tecnicas:**\n1. Diamond BLASTP alignment\n2. Transferencia de GO terms por similitud\n3. Ponderacion por % identidad\n4. GO propagation"}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "import os\nimport gc\nimport pandas as pd\nimport numpy as np\nfrom collections import defaultdict\nfrom tqdm.auto import tqdm\n\nprint('Diamond BLASTP Solution v1')"}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "# ==========================================\n# 1. SETUP DIAMOND\n# ==========================================\nprint('[1/7] Setting up Diamond...')\n\nif not os.path.exists('diamond'):\n    !wget -q https://github.com/bbuchfink/diamond/releases/download/v2.1.8/diamond-linux64.tar.gz\n    !tar xzf diamond-linux64.tar.gz\n    !chmod +x diamond\n    print('Diamond installed!')\nelse:\n    print('Diamond already exists')"}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "# ==========================================\n# 2. PATHS\n# ==========================================\nTRAIN_FASTA = '/kaggle/input/cafa-6-protein-function-prediction/Train/train_sequences.fasta'\nTEST_FASTA = '/kaggle/input/cafa-6-protein-function-prediction/Test/testsuperset.fasta'\nTRAIN_TERMS = '/kaggle/input/cafa-6-protein-function-prediction/Train/train_terms.tsv'\nGO_OBO = '/kaggle/input/cafa-6-protein-function-prediction/Train/go-basic.obo'\n\n# Check files exist\nfor f in [TRAIN_FASTA, TEST_FASTA, TRAIN_TERMS, GO_OBO]:\n    if os.path.exists(f):\n        print(f'OK: {os.path.basename(f)}')\n    else:\n        print(f'MISSING: {f}')"}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "# ==========================================\n# 3. CREATE DIAMOND DATABASE\n# ==========================================\nprint('[2/7] Creating Diamond database from training sequences...')\n\n!./diamond makedb --in {TRAIN_FASTA} -d train_db --quiet\n\nif os.path.exists('train_db.dmnd'):\n    print('Database created successfully!')\nelse:\n    print('ERROR: Database creation failed')"}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "# ==========================================\n# 4. RUN DIAMOND BLASTP\n# ==========================================\nprint('[3/7] Running Diamond BLASTP alignment...')\nprint('This may take 10-15 minutes...')\n\n# Run alignment - get top 10 hits per query with identity percentage\n!./diamond blastp \\\n    -d train_db \\\n    -q {TEST_FASTA} \\\n    --outfmt 6 qseqid sseqid pident evalue bitscore \\\n    --top 10 \\\n    --more-sensitive \\\n    --quiet \\\n    --out diamond_results.tsv\n\nif os.path.exists('diamond_results.tsv'):\n    results = pd.read_csv('diamond_results.tsv', sep='\\t', header=None,\n                          names=['query', 'target', 'pident', 'evalue', 'bitscore'])\n    print(f'Alignment complete! {len(results):,} hits found')\n    print(f'Unique test proteins with hits: {results[\"query\"].nunique():,}')\n    print(results.head())\nelse:\n    print('ERROR: Alignment failed')"}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "# ==========================================\n# 5. LOAD TRAINING GO TERMS\n# ==========================================\nprint('[4/7] Loading training GO terms...')\n\ntrain_terms = pd.read_csv(TRAIN_TERMS, sep='\\t', header=None, \n                          names=['protein', 'go_term', 'aspect'])\nprint(f'Loaded {len(train_terms):,} GO annotations')\nprint(f'Unique proteins: {train_terms[\"protein\"].nunique():,}')\nprint(f'Unique GO terms: {train_terms[\"go_term\"].nunique():,}')\n\n# Create protein -> GO terms mapping\nprotein_to_go = defaultdict(set)\nfor _, row in tqdm(train_terms.iterrows(), total=len(train_terms), desc='Building GO map'):\n    protein_to_go[row['protein']].add(row['go_term'])\n\nprint(f'GO map built for {len(protein_to_go):,} proteins')"}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "# ==========================================\n# 6. TRANSFER GO TERMS BY SIMILARITY\n# ==========================================\nprint('[5/7] Transferring GO terms based on sequence similarity...')\n\n# Load alignment results\nresults = pd.read_csv('diamond_results.tsv', sep='\\t', header=None,\n                      names=['query', 'target', 'pident', 'evalue', 'bitscore'])\n\n# Filter by minimum identity (e.g., 30%)\nMIN_IDENTITY = 30.0\nresults = results[results['pident'] >= MIN_IDENTITY]\nprint(f'Hits after filtering (>={MIN_IDENTITY}% identity): {len(results):,}')\n\n# Transfer GO terms weighted by identity\npredictions = defaultdict(lambda: defaultdict(float))\n\nfor _, row in tqdm(results.iterrows(), total=len(results), desc='Transferring GO terms'):\n    query = row['query']\n    target = row['target']\n    \n    # FIX: Parse UniProt ID from FASTA header format (sp|ID|NAME or tr|ID|NAME)\n    if '|' in target:\n        target = target.split('|')[1]\n    \n    pident = row['pident'] / 100.0  # Convert to 0-1 range\n    \n    # Get GO terms from target protein\n    if target in protein_to_go:\n        for go_term in protein_to_go[target]:\n            # Score = identity percentage (higher similarity = higher confidence)\n            predictions[query][go_term] = max(predictions[query][go_term], pident)\n\nprint(f'Predictions generated for {len(predictions):,} test proteins')"}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "# ==========================================\n# 7. PARSE GO ONTOLOGY FOR PROPAGATION\n# ==========================================\nprint('[6/7] Parsing GO ontology and propagating...')\n\ndef parse_obo(obo_path):\n    parents = defaultdict(set)\n    with open(obo_path, 'r') as f:\n        cur_id = None\n        for line in f:\n            line = line.strip()\n            if line == '[Term]':\n                cur_id = None\n            elif line.startswith('id: '):\n                cur_id = line.split('id: ')[1].strip()\n            elif line.startswith('is_a: '):\n                pid = line.split()[1].strip()\n                if cur_id:\n                    parents[cur_id].add(pid)\n            elif line.startswith('relationship: part_of '):\n                parts = line.split()\n                if len(parts) >= 3:\n                    pid = parts[2].strip()\n                    if cur_id:\n                        parents[cur_id].add(pid)\n    return parents\n\ndef get_ancestors(term, parents_map, cache):\n    if term in cache:\n        return cache[term]\n    ancestors = set()\n    stack = [term]\n    while stack:\n        cur = stack.pop()\n        for p in parents_map.get(cur, []):\n            if p not in ancestors:\n                ancestors.add(p)\n                stack.append(p)\n    cache[term] = ancestors\n    return ancestors\n\nparents_map = parse_obo(GO_OBO)\nprint(f'Parsed {len(parents_map)} GO terms')\n\n# Propagate predictions to ancestors\nROOTS = {'GO:0003674', 'GO:0008150', 'GO:0005575'}\nancestor_cache = {}\nfinal_rows = []\n\nfor protein, terms_dict in tqdm(predictions.items(), desc='Propagating'):\n    final_scores = dict(terms_dict)\n    \n    # Propagate to ancestors\n    for term, score in list(terms_dict.items()):\n        for ancestor in get_ancestors(term, parents_map, ancestor_cache):\n            final_scores[ancestor] = max(final_scores.get(ancestor, 0), score)\n    \n    # Set roots to 1.0\n    for root in ROOTS:\n        final_scores[root] = 1.0\n    \n    # Collect results (threshold 0.01)\n    for term, score in final_scores.items():\n        if score >= 0.01:\n            final_rows.append((protein, term, round(score, 6)))\n\nprint(f'Total predictions after propagation: {len(final_rows):,}')"}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "# ==========================================\n# 8. SAVE SUBMISSION\n# ==========================================\nprint('[7/7] Saving submission...')\n\nfinal_df = pd.DataFrame(final_rows, columns=['protein_id', 'go_term', 'score'])\nfinal_df = final_df.sort_values(['protein_id', 'score'], ascending=[True, False])\nfinal_df.to_csv('submission.tsv', sep='\\t', index=False, header=False)\n\nprint(f'\\n{\"=\"*60}')\nprint('DONE - Diamond BLASTP Solution v1')\nprint(f'Predictions: {len(final_df):,}')\nprint(f'Proteins: {final_df[\"protein_id\"].nunique():,}')\nprint(f'{\"=\"*60}')\nprint(final_df.head(10))"}], "metadata": {"kernelspec": {"display_name": "Python 3", "language": "python", "name": "python3"}, "language_info": {"name": "python", "version": "3.10.0"}}, "nbformat": 4, "nbformat_minor": 4}