{"cells":[{"cell_type":"markdown","id":"title","metadata":{},"source":"# CAFA-6 GOA + BLAST Baseline (v8 - Complete v3 Reproduction)\n\n**Purpose**: v8 is an exact reproduction of v3 to confirm the 0.036 score and validate that v3's parameters are the correct baseline.\n\nThis notebook implements an optimized baseline for the CAFA-6 protein function prediction competition:\n\n1. **Diamond BLAST** - Find similar proteins in training set with optimized parameters\n2. **GOA Transfer** - Transfer GO annotations using adaptive multi-score fusion\n3. **GO Propagation** - Apply True Path Rule with Noisy-OR probabilistic combination\n\n**Why v8 (Reproducing v3)**:\n- v7 failed with score 0.023 (expected 0.32-0.36)\n- Analysis revealed v7 used different parameters than v3:\n  - v7: propagation_method=\"max\", v3: \"noisy_or\"\n  - v7: n_hits=50, v3: 75\n  - v7: min_identity=30.0, v3: 25.0\n  - v7: threshold=0.01, v3: 0.005\n  - v7: score_method=\"bitscore_weighted\", v3: \"multi_fusion\"\n- v8 reproduces v3 exactly to confirm 0.036 score\n\n**v3 Key Parameters (Preserved in v8)**:\n- **Fixed submission format**: Headerless TSV as required by Kaggle\n- **Adaptive multi-score fusion**: 40% bitscore + 40% identity + 20% e-value for better ranking\n- **Optimized threshold**: 0.005 (lower than v2) for better recall\n- **Increased coverage**: n_hits=75, min_identity=25% for more matches\n- **Noisy-OR propagation**: Probabilistic evidence combination\n- **Memory optimization**: Batch processing to handle large prediction sets\n\n**Expected Score**: 0.036 (matching v3)"},{"cell_type":"markdown","id":"setup-header","metadata":{},"source":"## 1. Environment Setup\n\nInstall Diamond BLAST from precompiled binary (no conda required)."},{"cell_type":"code","execution_count":null,"id":"install-diamond","metadata":{},"outputs":[],"source":"# Install Diamond BLAST from precompiled binary\n!wget -q https://github.com/bbuchfink/diamond/releases/download/v2.1.9/diamond-linux64.tar.gz\n!tar xzf diamond-linux64.tar.gz\n!mv diamond /usr/local/bin/\n!chmod +x /usr/local/bin/diamond\n!diamond version"},{"cell_type":"code","execution_count":null,"id":"install-deps","metadata":{},"outputs":[],"source":"# Install Python dependencies\n!pip install -q obonet networkx biopython tqdm"},{"cell_type":"markdown","id":"download-header","metadata":{},"source":"## 2. Download External Data\n\nDownload GO ontology and GOA annotations from UniProt."},{"cell_type":"code","execution_count":null,"id":"download-data","metadata":{},"outputs":[],"source":"import os\n\n# Create working directories\nos.makedirs('/kaggle/working/external', exist_ok=True)\n\n# Download GO Ontology\n!wget -q http://purl.obolibrary.org/obo/go/go-basic.obo -O /kaggle/working/external/go-basic.obo\nprint('Downloaded GO ontology')\n\n# Copy GOA data from Kaggle dataset (experimental evidence only)\n# Note: Kaggle automatically extracts .gz files, so we use the uncompressed filename\nGOA_INPUT_PATH = '/kaggle/input/cafa6-goa-data/goa_experimental.gaf'\n!cp {GOA_INPUT_PATH} /kaggle/working/external/goa.gaf\nprint('Copied GOA data from Kaggle dataset')"},{"cell_type":"code","execution_count":null,"id":"verify-files","metadata":{},"outputs":[],"source":"# Verify downloaded files\n!ls -lh /kaggle/working/external/\n!head -5 /kaggle/working/external/go-basic.obo\n!head -5 /kaggle/working/external/goa.gaf"},{"cell_type":"markdown","id":"imports-header","metadata":{},"source":"## 3. Import Libraries"},{"cell_type":"code","execution_count":null,"id":"imports","metadata":{},"outputs":[],"source":"import gzip\nimport logging\nimport subprocess\nimport tempfile\nfrom collections import defaultdict\nfrom dataclasses import dataclass\nfrom pathlib import Path\nfrom typing import Optional\n\nimport networkx as nx\nimport numpy as np\nimport obonet\nimport pandas as pd\nfrom tqdm.auto import tqdm\n\n# Setup logging\nlogging.basicConfig(level=logging.INFO, format='%(asctime)s - %(levelname)s - %(message)s')\nlogger = logging.getLogger(__name__)\n\nprint('Libraries imported successfully')"},{"cell_type":"markdown","id":"gograph-header","metadata":{},"source":"## 4. GO Graph Class\n\nLoad and navigate the Gene Ontology hierarchy."},{"cell_type":"code","execution_count":null,"id":"gograph","metadata":{},"outputs":[],"source":"class GOGraph:\n    \"\"\"Gene Ontology graph for term relationships.\"\"\"\n\n    def __init__(self, obo_path: str):\n        \"\"\"Load GO graph from OBO file.\"\"\"\n        self.obo_path = Path(obo_path)\n        self._load_graph()\n\n    def _load_graph(self):\n        \"\"\"Load and process the GO graph.\"\"\"\n        logger.info(f\"Loading GO graph from {self.obo_path}\")\n        self.graph = obonet.read_obo(str(self.obo_path))\n        logger.info(f\"Loaded {len(self.graph.nodes)} GO terms\")\n\n        # Create mapping from GO ID to term info\n        self.term_info = {}\n        for node, data in self.graph.nodes(data=True):\n            self.term_info[node] = {\n                \"name\": data.get(\"name\", \"\"),\n                \"namespace\": data.get(\"namespace\", \"\"),\n                \"is_obsolete\": data.get(\"is_obsolete\", False),\n            }\n\n        # Build caches\n        self._ancestor_cache: dict = {}\n        self._descendant_cache: dict = {}\n\n    def get_namespace(self, go_id: str) -> str:\n        \"\"\"Get the namespace of a GO term.\"\"\"\n        if go_id in self.term_info:\n            return self.term_info[go_id].get(\"namespace\", \"\")\n        return \"\"\n\n    def get_ontology(self, go_id: str) -> str:\n        \"\"\"Get the ontology code (BPO, MFO, CCO) of a GO term.\"\"\"\n        namespace = self.get_namespace(go_id)\n        namespace_map = {\n            \"biological_process\": \"BPO\",\n            \"molecular_function\": \"MFO\",\n            \"cellular_component\": \"CCO\",\n        }\n        return namespace_map.get(namespace, \"\")\n\n    def is_valid_term(self, go_id: str) -> bool:\n        \"\"\"Check if a GO term is valid (exists and not obsolete).\"\"\"\n        if go_id not in self.term_info:\n            return False\n        return not self.term_info[go_id].get(\"is_obsolete\", False)\n\n    def get_ancestors(self, go_id: str, include_self: bool = False) -> set:\n        \"\"\"Get all ancestor terms of a GO term.\"\"\"\n        if go_id in self._ancestor_cache:\n            ancestors = self._ancestor_cache[go_id].copy()\n        else:\n            ancestors = set()\n            if go_id in self.graph:\n                ancestors = nx.ancestors(self.graph, go_id)\n            self._ancestor_cache[go_id] = ancestors.copy()\n\n        if include_self:\n            ancestors.add(go_id)\n        return ancestors\n\n    def get_descendants(self, go_id: str, include_self: bool = False) -> set:\n        \"\"\"Get all descendant terms of a GO term.\"\"\"\n        if go_id in self._descendant_cache:\n            descendants = self._descendant_cache[go_id].copy()\n        else:\n            descendants = set()\n            if go_id in self.graph:\n                descendants = nx.descendants(self.graph, go_id)\n            self._descendant_cache[go_id] = descendants.copy()\n\n        if include_self:\n            descendants.add(go_id)\n        return descendants"},{"cell_type":"markdown","id":"propagator-header","metadata":{},"source":"## 5. GO Propagator Class\n\nPropagate GO term predictions following the True Path Rule."},{"cell_type":"code","execution_count":null,"id":"propagator","metadata":{},"outputs":[],"source":"import gc\n\nclass GOPropagator:\n    \"\"\"Propagates GO term predictions according to True Path Rule.\"\"\"\n\n    def __init__(self, go_graph: GOGraph):\n        self.go_graph = go_graph\n\n    def propagate_predictions(\n        self,\n        predictions: dict,\n        propagation_method: str = \"max\",\n        batch_size: int = 10000,\n    ) -> dict:\n        \"\"\"Propagate predictions to ancestor terms with batch processing for memory efficiency.\"\"\"\n        logger.info(f\"Propagating predictions using method: {propagation_method} (batch_size={batch_size})\")\n\n        protein_ids = list(predictions.keys())\n        propagated = {}\n        \n        for i in range(0, len(protein_ids), batch_size):\n            batch_ids = protein_ids[i:i + batch_size]\n            logger.info(f\"Processing batch {i // batch_size + 1}/{(len(protein_ids) + batch_size - 1) // batch_size}\")\n            \n            for protein_id in tqdm(batch_ids, desc=f\"Propagating batch\"):\n                propagated[protein_id] = self._propagate_single(predictions[protein_id], propagation_method)\n            \n            # Free memory after each batch\n            gc.collect()\n\n        return propagated\n\n    def _propagate_single(self, term_scores: dict, method: str) -> dict:\n        \"\"\"Propagate scores for a single protein.\"\"\"\n        result = term_scores.copy()\n        ancestor_scores: dict = defaultdict(list)\n\n        for go_term, score in term_scores.items():\n            if not self.go_graph.is_valid_term(go_term):\n                continue\n            ancestors = self.go_graph.get_ancestors(go_term)\n            for ancestor in ancestors:\n                if self.go_graph.is_valid_term(ancestor):\n                    ancestor_scores[ancestor].append(score)\n\n        for ancestor, scores in ancestor_scores.items():\n            if ancestor in result:\n                if method == \"max\":\n                    result[ancestor] = max(result[ancestor], max(scores))\n                elif method == \"mean\":\n                    all_scores = scores + [result[ancestor]]\n                    result[ancestor] = sum(all_scores) / len(all_scores)\n                elif method == \"noisy_or\":\n                    all_scores = scores + [result[ancestor]]\n                    prob_all_false = 1.0\n                    for s in all_scores:\n                        prob_all_false *= 1 - s\n                    result[ancestor] = 1 - prob_all_false\n            else:\n                if method == \"max\":\n                    result[ancestor] = max(scores)\n                elif method == \"mean\":\n                    result[ancestor] = sum(scores) / len(scores)\n                elif method == \"noisy_or\":\n                    prob_all_false = 1.0\n                    for s in scores:\n                        prob_all_false *= 1 - s\n                    result[ancestor] = 1 - prob_all_false\n\n        return result\n\n    def apply_negative_propagation(\n        self,\n        predictions: dict,\n        negative_threshold: float = 0.01,\n        batch_size: int = 10000,\n    ) -> dict:\n        \"\"\"Apply negative propagation with batch processing.\"\"\"\n        logger.info(f\"Applying negative propagation with threshold: {negative_threshold}\")\n\n        protein_ids = list(predictions.keys())\n        result = {}\n        \n        for i in range(0, len(protein_ids), batch_size):\n            batch_ids = protein_ids[i:i + batch_size]\n            \n            for protein_id in batch_ids:\n                result[protein_id] = self._negative_propagate_single(predictions[protein_id], negative_threshold)\n            \n            gc.collect()\n\n        return result\n\n    def _negative_propagate_single(self, term_scores: dict, threshold: float) -> dict:\n        \"\"\"Apply negative propagation for a single protein.\"\"\"\n        result = term_scores.copy()\n        negative_terms = {t for t, s in term_scores.items() if s < threshold}\n\n        for neg_term in negative_terms:\n            descendants = self.go_graph.get_descendants(neg_term)\n            for desc in descendants:\n                if desc in result:\n                    result[desc] = min(result[desc], term_scores[neg_term])\n\n        return result"},{"cell_type":"markdown","id":"goa-header","metadata":{},"source":"## 6. GOA Parser\n\nParse Gene Ontology Annotation (GOA) database."},{"cell_type":"code","execution_count":null,"id":"goa-parser","metadata":{},"outputs":[],"source":"# Evidence codes that indicate experimental validation\nEXPERIMENTAL_EVIDENCE_CODES = frozenset([\n    \"EXP\", \"IDA\", \"IPI\", \"IMP\", \"IGI\", \"IEP\",\n    \"TAS\", \"IC\", \"HTP\", \"HDA\", \"HMP\", \"HGI\", \"HEP\"\n])\n\nASPECT_TO_ONTOLOGY = {\"P\": \"BPO\", \"F\": \"MFO\", \"C\": \"CCO\"}\n\n\ndef parse_goa(\n    gaf_path: str,\n    experimental_only: bool = True,\n    chunk_size: int = 100000,\n) -> pd.DataFrame:\n    \"\"\"Parse GOA GAF file.\"\"\"\n    gaf_path = Path(gaf_path)\n    annotations = []\n\n    logger.info(f\"Parsing GOA file: {gaf_path}\")\n\n    open_fn = gzip.open if str(gaf_path).endswith(\".gz\") else open\n    mode = \"rt\" if str(gaf_path).endswith(\".gz\") else \"r\"\n\n    with open_fn(gaf_path, mode) as f:\n        batch = []\n        for line in tqdm(f, desc=\"Parsing GAF\"):\n            if line.startswith(\"!\"):\n                continue\n\n            fields = line.strip().split(\"\\t\")\n            if len(fields) < 15:\n                continue\n\n            protein_id = fields[1]\n            go_id = fields[4]\n            evidence_code = fields[6]\n            aspect = fields[8]\n\n            if experimental_only and evidence_code not in EXPERIMENTAL_EVIDENCE_CODES:\n                continue\n\n            batch.append({\n                \"protein_id\": protein_id,\n                \"go_id\": go_id,\n                \"aspect\": aspect,\n                \"evidence_code\": evidence_code,\n                \"ontology\": ASPECT_TO_ONTOLOGY.get(aspect, \"\"),\n            })\n\n            if len(batch) >= chunk_size:\n                annotations.extend(batch)\n                batch = []\n\n        if batch:\n            annotations.extend(batch)\n\n    df = pd.DataFrame(annotations)\n    logger.info(f\"Parsed {len(df)} annotations\")\n\n    # Remove duplicates\n    df = df.drop_duplicates(subset=[\"protein_id\", \"go_id\"])\n    logger.info(f\"After deduplication: {len(df)} annotations\")\n\n    return df"},{"cell_type":"markdown","id":"diamond-header","metadata":{},"source":"## 7. Diamond BLAST Wrapper\n\nFast protein sequence alignment using Diamond BLAST."},{"cell_type":"code","execution_count":null,"id":"diamond","metadata":{},"outputs":[],"source":"class DiamondBlast:\n    \"\"\"Wrapper for Diamond BLAST sequence alignment.\"\"\"\n\n    OUTFMT_COLUMNS = [\n        \"qseqid\", \"sseqid\", \"pident\", \"length\", \"mismatch\", \"gapopen\",\n        \"qstart\", \"qend\", \"sstart\", \"send\", \"evalue\", \"bitscore\"\n    ]\n\n    def __init__(self, db_path=None, threads: int = 4, sensitivity: str = \"very-sensitive\"):\n        self.db_path = db_path\n        self.threads = threads\n        self.sensitivity = sensitivity\n\n    def make_database(self, fasta_path: str, output_path: str):\n        \"\"\"Create Diamond database from FASTA file.\"\"\"\n        output_path = Path(output_path)\n        output_path.parent.mkdir(parents=True, exist_ok=True)\n\n        cmd = [\n            \"diamond\", \"makedb\",\n            \"--in\", str(fasta_path),\n            \"--db\", str(output_path),\n            \"--threads\", str(self.threads),\n        ]\n\n        logger.info(f\"Creating Diamond database: {' '.join(cmd)}\")\n        result = subprocess.run(cmd, capture_output=True, text=True)\n\n        if result.returncode != 0:\n            raise RuntimeError(f\"Diamond makedb failed: {result.stderr}\")\n\n        db_file = Path(str(output_path) + \".dmnd\")\n        logger.info(f\"Database created: {db_file}\")\n        self.db_path = db_file\n        return db_file\n\n    def search(\n        self,\n        query_fasta: str,\n        output_path: str = None,\n        max_target_seqs: int = 50,\n        evalue: float = 1e-3,\n        min_identity: float = 0.0,\n    ) -> pd.DataFrame:\n        \"\"\"Search query sequences against database.\"\"\"\n        if self.db_path is None:\n            raise ValueError(\"Database path not set. Call make_database() first.\")\n\n        if output_path is None:\n            tmp = tempfile.NamedTemporaryFile(mode=\"w\", suffix=\".tsv\", delete=False)\n            output_path = Path(tmp.name)\n            tmp.close()\n        else:\n            output_path = Path(output_path)\n            output_path.parent.mkdir(parents=True, exist_ok=True)\n\n        cmd = [\n            \"diamond\", \"blastp\",\n            \"--query\", str(query_fasta),\n            \"--db\", str(self.db_path),\n            \"--out\", str(output_path),\n            f\"--{self.sensitivity}\",\n            \"--threads\", str(self.threads),\n            \"--max-target-seqs\", str(max_target_seqs),\n            \"--evalue\", str(evalue),\n            \"--outfmt\", \"6\", *self.OUTFMT_COLUMNS,\n        ]\n\n        if min_identity > 0:\n            cmd.extend([\"--id\", str(min_identity)])\n\n        logger.info(f\"Running Diamond BLAST: {' '.join(cmd)}\")\n        result = subprocess.run(cmd, capture_output=True, text=True)\n\n        if result.returncode != 0:\n            raise RuntimeError(f\"Diamond blastp failed: {result.stderr}\")\n\n        df = pd.read_csv(\n            output_path, sep=\"\\t\", names=self.OUTFMT_COLUMNS,\n            dtype={\"qseqid\": str, \"sseqid\": str, \"pident\": float,\n                   \"length\": int, \"mismatch\": int, \"gapopen\": int,\n                   \"qstart\": int, \"qend\": int, \"sstart\": int,\n                   \"send\": int, \"evalue\": float, \"bitscore\": float}\n        )\n\n        logger.info(f\"Found {len(df)} hits for {df['qseqid'].nunique()} queries\")\n        return df\n\n    def get_best_hits(\n        self,\n        hits_df: pd.DataFrame,\n        n_best: int = 10,\n        min_identity: float = 30.0,\n        max_evalue: float = 1e-5,\n    ) -> pd.DataFrame:\n        \"\"\"Get best hits for each query.\"\"\"\n        filtered = hits_df[\n            (hits_df[\"pident\"] >= min_identity) & (hits_df[\"evalue\"] <= max_evalue)\n        ].copy()\n\n        filtered = filtered.sort_values(\"bitscore\", ascending=False)\n        best_hits = filtered.groupby(\"qseqid\").head(n_best)\n\n        logger.info(f\"Filtered to {len(best_hits)} best hits for {best_hits['qseqid'].nunique()} queries\")\n        return best_hits\n\n\n\n\ndef compute_blast_scores(hits_df: pd.DataFrame, method: str = \"multi_fusion\") -> pd.DataFrame:\n    \"\"\"Compute similarity scores from BLAST hits with enhanced scoring methods.\"\"\"\n    df = hits_df.copy()\n\n    if method == \"multi_fusion\":\n        # Multi-score fusion: combine bitscore, identity, and e-value\n        # Normalize bitscore per query\n        df[\"bitscore_norm\"] = df.groupby(\"qseqid\")[\"bitscore\"].transform(\n            lambda x: x / x.max() if x.max() > 0 else 0\n        )\n        # Normalize identity (already 0-100)\n        df[\"identity_norm\"] = df[\"pident\"] / 100.0\n        # Normalize e-value (log scale)\n        df[\"evalue_norm\"] = -np.log10(df[\"evalue\"].clip(lower=1e-180))\n        df[\"evalue_norm\"] = df.groupby(\"qseqid\")[\"evalue_norm\"].transform(\n            lambda x: x / x.max() if x.max() > 0 else 0\n        )\n        # Weighted combination: bitscore (40%), identity (40%), e-value (20%)\n        df[\"score\"] = (\n            0.4 * df[\"bitscore_norm\"] +\n            0.4 * df[\"identity_norm\"] +\n            0.2 * df[\"evalue_norm\"]\n        )\n    elif method == \"bitscore_weighted\":\n        df[\"score\"] = df.groupby(\"qseqid\")[\"bitscore\"].transform(\n            lambda x: x / x.max() if x.max() > 0 else 0\n        )\n    elif method == \"identity_weighted\":\n        df[\"score\"] = df[\"pident\"] / 100.0\n    elif method == \"evalue_weighted\":\n        df[\"score\"] = -np.log10(df[\"evalue\"].clip(lower=1e-180))\n        df[\"score\"] = df.groupby(\"qseqid\")[\"score\"].transform(\n            lambda x: x / x.max() if x.max() > 0 else 0\n        )\n    else:\n        raise ValueError(f\"Unknown method: {method}\")\n\n    return df[[\"qseqid\", \"sseqid\", \"score\"]].rename(\n        columns={\"qseqid\": \"query_id\", \"sseqid\": \"subject_id\"}\n    )\n    )"},{"cell_type":"markdown","id":"baseline-header","metadata":{},"source":"## 8. GOA + BLAST Baseline Class\n\nMain prediction pipeline that combines BLAST search with GOA transfer."},{"cell_type":"code","execution_count":null,"id":"baseline","metadata":{},"outputs":[],"source":"class GOABlastBaseline:\n    \"\"\"GOA + BLAST baseline predictor.\"\"\"\n\n    def __init__(\n        self,\n        go_graph: GOGraph,\n        goa_df: pd.DataFrame,\n        blast: DiamondBlast,\n        propagation_method: str = \"max\",\n        use_negative_propagation: bool = True,\n        negative_threshold: float = 0.05,\n    ):\n        self.go_graph = go_graph\n        self.blast = blast\n        self.propagation_method = propagation_method\n        self.use_negative_propagation = use_negative_propagation\n        self.negative_threshold = negative_threshold\n\n        # Build protein -> GO terms mapping from GOA\n        self.protein_annotations = self._build_annotation_map(goa_df)\n        logger.info(f\"Loaded annotations for {len(self.protein_annotations)} proteins\")\n\n        # Initialize propagator\n        self.propagator = GOPropagator(go_graph)\n\n    def _build_annotation_map(self, goa_df: pd.DataFrame) -> dict:\n        \"\"\"Build protein -> {go_term: confidence} mapping.\"\"\"\n        annotations = defaultdict(dict)\n        for _, row in goa_df.iterrows():\n            protein_id = row[\"protein_id\"]\n            go_id = row[\"go_id\"]\n            annotations[protein_id][go_id] = 1.0\n        return dict(annotations)\n\n    def _extract_accession(self, protein_id: str) -> str:\n        \"\"\"Extract UniProt accession from various ID formats.\"\"\"\n        if \"|\" in protein_id:\n            parts = protein_id.split(\"|\")\n            if len(parts) >= 2:\n                return parts[1]\n        return protein_id\n\n    def predict(\n        self,\n        query_fasta: str,\n        n_hits: int = 50,\n        min_identity: float = 30.0,\n        max_evalue: float = 1e-3,\n        score_method: str = \"bitscore_weighted\",\n    ) -> dict:\n        \"\"\"Predict GO terms for query proteins.\"\"\"\n        logger.info(\"Running BLAST search...\")\n        hits_df = self.blast.search(\n            query_fasta=query_fasta,\n            max_target_seqs=n_hits,\n            evalue=max_evalue,\n        )\n\n        if len(hits_df) == 0:\n            logger.warning(\"No BLAST hits found!\")\n            return {}\n\n        # Filter hits\n        logger.info(\"Filtering BLAST hits...\")\n        best_hits = self.blast.get_best_hits(\n            hits_df, n_best=n_hits, min_identity=min_identity, max_evalue=max_evalue\n        )\n        \n        # Free memory\n        del hits_df\n        gc.collect()\n\n        # Compute similarity scores\n        logger.info(f\"Computing similarity scores using {score_method}...\")\n        scores_df = compute_blast_scores(best_hits, method=score_method)\n        \n        del best_hits\n        gc.collect()\n\n        # Transfer annotations\n        logger.info(\"Transferring GO annotations...\")\n        predictions = self._transfer_annotations(scores_df)\n        \n        del scores_df\n        gc.collect()\n\n        # Apply propagation\n        logger.info(\"Applying GO propagation...\")\n        predictions = self.propagator.propagate_predictions(\n            predictions, self.propagation_method\n        )\n\n        # Apply negative propagation\n        if self.use_negative_propagation:\n            logger.info(\"Applying negative propagation...\")\n            predictions = self.propagator.apply_negative_propagation(\n                predictions, self.negative_threshold\n            )\n\n        return predictions\n\n    def _transfer_annotations(self, scores_df: pd.DataFrame) -> dict:\n        \"\"\"Transfer GO annotations from similar proteins.\"\"\"\n        predictions = defaultdict(dict)\n        matches_found = 0\n\n        for _, row in tqdm(scores_df.iterrows(), desc=\"Transferring annotations\"):\n            query_id = row[\"query_id\"]\n            subject_id = row[\"subject_id\"]\n            similarity_score = row[\"score\"]\n\n            subject_accession = self._extract_accession(subject_id)\n\n            if subject_accession not in self.protein_annotations:\n                continue\n\n            matches_found += 1\n            for go_term, annotation_conf in self.protein_annotations[subject_accession].items():\n                term_score = similarity_score * annotation_conf\n\n                if go_term in predictions[query_id]:\n                    predictions[query_id][go_term] = max(\n                        predictions[query_id][go_term], term_score\n                    )\n                else:\n                    predictions[query_id][go_term] = term_score\n\n        logger.info(f\"Found {matches_found} BLAST hits with GOA annotations\")\n        return dict(predictions)\n\n    def predict_to_tsv(\n        self,\n        query_fasta: str,\n        output_path: str,\n        threshold: float = 0.01,\n        max_terms_per_protein: int = 1500,\n        batch_size: int = 5000,\n        **predict_kwargs\n    ) -> tuple:\n        \"\"\"Predict and write directly to TSV file (memory efficient).\n        \n        Returns:\n            tuple: (total_predictions, unique_proteins)\n        \"\"\"\n        logger.info(\"Running BLAST search...\")\n        hits_df = self.blast.search(\n            query_fasta=query_fasta,\n            max_target_seqs=predict_kwargs.get('n_hits', 50),\n            evalue=predict_kwargs.get('max_evalue', 1e-3),\n        )\n\n        if len(hits_df) == 0:\n            logger.warning(\"No BLAST hits found!\")\n            return 0, 0\n\n        # Filter hits\n        logger.info(\"Filtering BLAST hits...\")\n        best_hits = self.blast.get_best_hits(\n            hits_df,\n            n_best=predict_kwargs.get('n_hits', 50),\n            min_identity=predict_kwargs.get('min_identity', 30.0),\n            max_evalue=predict_kwargs.get('max_evalue', 1e-3)\n        )\n        \n        del hits_df\n        gc.collect()\n\n        # Compute similarity scores\n        score_method = predict_kwargs.get('score_method', 'bitscore_weighted')\n        logger.info(f\"Computing similarity scores using {score_method}...\")\n        scores_df = compute_blast_scores(best_hits, method=score_method)\n        \n        del best_hits\n        gc.collect()\n\n        # Transfer annotations\n        logger.info(\"Transferring GO annotations...\")\n        predictions = self._transfer_annotations(scores_df)\n        \n        del scores_df\n        gc.collect()\n\n        # Process predictions in batches and write directly to TSV\n        logger.info(\"Processing predictions and writing to TSV...\")\n        protein_ids = list(predictions.keys())\n        total_predictions = 0\n        unique_proteins = set()\n        \n        # Initialize empty file (NO header - Kaggle expects headerless TSV)\n        with open(output_path, 'w') as f:\n            pass\n        \n        num_batches = (len(protein_ids) + batch_size - 1) // batch_size\n        \n        for batch_idx in range(0, len(protein_ids), batch_size):\n            batch_ids = protein_ids[batch_idx:batch_idx + batch_size]\n            batch_num = batch_idx // batch_size + 1\n            logger.info(f\"Processing batch {batch_num}/{num_batches}\")\n            \n            # Propagate this batch\n            batch_predictions = {pid: predictions[pid] for pid in batch_ids}\n            propagated_batch = self.propagator.propagate_predictions(\n                batch_predictions, self.propagation_method, batch_size=batch_size\n            )\n            \n            # Free original predictions for this batch\n            for pid in batch_ids:\n                del predictions[pid]\n            del batch_predictions\n            gc.collect()\n            \n            # Write batch to TSV\n            lines = []\n            for protein_id in batch_ids:\n                if protein_id not in propagated_batch:\n                    continue\n                    \n                term_scores = propagated_batch[protein_id]\n                \n                # Filter by threshold and sort by score\n                filtered_terms = [\n                    (go_term, score) for go_term, score in term_scores.items()\n                    if score > threshold and self.go_graph.is_valid_term(go_term)\n                ]\n                filtered_terms.sort(key=lambda x: x[1], reverse=True)\n                \n                # Limit terms per protein\n                for go_term, score in filtered_terms[:max_terms_per_protein]:\n                    lines.append(f\"{protein_id}\\t{go_term}\\t{score:.3f}\\n\")\n                    total_predictions += 1\n                \n                if filtered_terms:\n                    unique_proteins.add(protein_id)\n            \n            # Append to file\n            with open(output_path, 'a') as f:\n                f.writelines(lines)\n            \n            del propagated_batch\n            del lines\n            gc.collect()\n            \n            logger.info(f\"Batch {batch_num}: wrote predictions, total so far: {total_predictions}\")\n        \n        logger.info(f\"Finished writing {total_predictions} predictions for {len(unique_proteins)} proteins\")\n        return total_predictions, len(unique_proteins)"},{"cell_type":"markdown","id":"run-header","metadata":{},"source":"## 9. Load Data and Initialize\n\nLoad CAFA-6 competition data and initialize the baseline."},{"cell_type":"code","execution_count":null,"id":"paths","metadata":{},"outputs":[],"source":"# Paths\nDATA_DIR = Path('/kaggle/input/cafa-6-protein-function-prediction')\nWORKING_DIR = Path('/kaggle/working')\nEXTERNAL_DIR = WORKING_DIR / 'external'\n\n# CAFA-6 data paths\nTRAIN_FASTA = DATA_DIR / 'Train' / 'train_sequences.fasta'\nTEST_FASTA = DATA_DIR / 'Test' / 'testsuperset.fasta'\n\n# External data paths\nOBO_PATH = EXTERNAL_DIR / 'go-basic.obo'\nGOA_PATH = EXTERNAL_DIR / 'goa.gaf'\n\nprint(f\"Train FASTA: {TRAIN_FASTA}\")\nprint(f\"Test FASTA: {TEST_FASTA}\")\nprint(f\"GO OBO: {OBO_PATH}\")\nprint(f\"GOA data: {GOA_PATH}\")"},{"cell_type":"code","execution_count":null,"id":"load-go","metadata":{},"outputs":[],"source":"# Load GO graph\ngo_graph = GOGraph(str(OBO_PATH))\nprint(f\"Loaded GO graph with {len(go_graph.graph.nodes)} terms\")"},{"cell_type":"code","execution_count":null,"id":"load-goa","metadata":{},"outputs":[],"source":"# Load GOA data (experimental evidence only)\ngoa_df = parse_goa(str(GOA_PATH), experimental_only=True)\nprint(f\"Loaded {len(goa_df)} GOA annotations\")\nprint(f\"Unique proteins: {goa_df['protein_id'].nunique()}\")\nprint(f\"Unique GO terms: {goa_df['go_id'].nunique()}\")"},{"cell_type":"code","execution_count":null,"id":"create-db","metadata":{},"outputs":[],"source":"# Create BLAST database from training data\nblast = DiamondBlast(threads=4, sensitivity=\"very-sensitive\")\ndb_path = WORKING_DIR / 'diamond' / 'train_db'\nblast.make_database(str(TRAIN_FASTA), str(db_path))"},{"cell_type":"code","execution_count":null,"id":"init-baseline","metadata":{},"outputs":[],"source":"# Initialize baseline with optimized parameters\nbaseline = GOABlastBaseline(\n    go_graph=go_graph,\n    goa_df=goa_df,\n    blast=blast,\n    propagation_method=\"noisy_or\",  # Changed from \"max\" for probabilistic evidence combination\n    use_negative_propagation=False,  # Disabled to prevent memory overflow\n    negative_threshold=0.05,\n)\nprint(\"Baseline initialized with Noisy-OR propagation!\")"},{"cell_type":"markdown","id":"predict-header","metadata":{},"source":"## 10. Run Predictions\n\nGenerate predictions for the test set."},{"cell_type":"code","execution_count":null,"id":"predict","metadata":{},"outputs":[],"source":"# Run predictions with optimized parameters (v8 = v3)\noutput_path = WORKING_DIR / 'submission.tsv'\n\ntotal_preds, num_proteins = baseline.predict_to_tsv(\n    str(TEST_FASTA),\n    str(output_path),\n    threshold=0.005,  # v8 (same as v3): Lowered from 0.01 for better recall\n    max_terms_per_protein=1500,\n    batch_size=5000,\n    n_hits=75,  # v8 (same as v3): Increased from 50 for better coverage\n    min_identity=25.0,  # v8 (same as v3): Lowered from 30.0 for more matches\n    max_evalue=1e-3,  # Stricter from 1e-2\n    score_method=\"multi_fusion\",  # Enhanced: combines bitscore, identity, e-value\n)\n\nprint(f\"Generated {total_preds} predictions for {num_proteins} proteins\")"},{"cell_type":"code","execution_count":null,"id":"stats","metadata":{},"outputs":[],"source":"# Load submission and show statistics (headerless TSV)\nsubmission_df = pd.read_csv(output_path, sep='\\t', header=None, names=['Protein', 'GO_term', 'Confidence'])\nprint(f\"Submission file: {output_path}\")\nprint(f\"Total predictions: {len(submission_df)}\")\nprint(f\"Unique proteins: {submission_df['Protein'].nunique()}\")\nprint(f\"Unique GO terms: {submission_df['GO_term'].nunique()}\")\nprint(f\"\\nScore distribution:\")\nprint(submission_df['Confidence'].describe())"},{"cell_type":"markdown","id":"submission-header","metadata":{},"source":"## 11. Verify Submission File\n\nCheck the submission file format."},{"cell_type":"code","execution_count":null,"id":"submission","metadata":{},"outputs":[],"source":"# Show sample of submission\nprint(\"Sample predictions:\")\nsubmission_df.head(20)"},{"cell_type":"code","execution_count":null,"id":"save","metadata":{},"outputs":[],"source":"# Clean up DataFrame to free memory\ndel submission_df\ngc.collect()\nprint(\"Submission file ready at /kaggle/working/submission.tsv\")"},{"cell_type":"code","execution_count":null,"id":"verify","metadata":{},"outputs":[],"source":"# Verify submission file\n!head -20 /kaggle/working/submission.tsv\n!wc -l /kaggle/working/submission.tsv\n!ls -lh /kaggle/working/submission.tsv"},{"cell_type":"markdown","id":"end-header","metadata":{},"source":"## Done!\n\nThe submission file has been created at `/kaggle/working/submission.tsv`.\n\n**v8 Parameters (Exact v3 Reproduction)**:\n1. Fixed submission format (headerless TSV)\n2. Adaptive multi-score fusion (40% bitscore + 40% identity + 20% e-value)\n3. Optimized threshold (0.005) for better recall\n4. Increased coverage (75 hits, 25% min identity)\n5. Noisy-OR propagation (probabilistic evidence combination)\n\nNext steps:\n1. Submit to the competition\n2. Check the leaderboard score\n3. Expected score: 0.036 (matching v3)"}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.11.0"}},"nbformat":4,"nbformat_minor":5}