{"cells":[{"cell_type":"markdown","id":"title","metadata":{},"source":"# CAFA-6 GOA + BLAST Baseline v7 (External GOA Reversion)\n\nThis notebook implements the v7 strategy: **Reverting to external GOA data** after v6's catastrophic failure.\n\n## Critical Fix: External GOA Data\n\n**v6 mistake**: Used competition training labels (537K annotations, 82K proteins) ❌  \n**v7 strategy**: Use external GOA data (613K annotations, 111K proteins) ✅\n\n**Why this matters**:\n- External GOA has 30%+ better coverage than competition labels\n- Reference baselines (0.35+ score) all use external GOA\n- v6 score: 0.017 (failed), v5 score: 0.032 (better), v3 score: 0.036\n\n## Strategy: Reference Baseline Reproduction\n\n1. **Diamond BLAST** - Find similar proteins with reference parameters\n2. **GOA Transfer** - Transfer GO annotations from external GOA database\n3. **GO Propagation** - Apply True Path Rule with MAX method (reference standard)\n\n**Key Parameters** (from reference baselines):\n- `propagation_method=\"max\"` - Reference standard (not noisy_or)\n- `use_negative_propagation=False` - Disabled initially for stability\n- `score_method=\"bitscore_weighted\"` - Simple, proven scoring\n- `n_hits=50` - Reference standard\n- `max_evalue=1e-3` - Reference threshold\n- `min_identity=30.0` - Reference minimum\n- `threshold=0.01` - Reference output threshold\n\n**Expected Score**: 0.32-0.36 (recovering to v3 level and approaching reference baselines)\n\n**Success Metrics**:\n- GOA proteins loaded: >= 100,000 (target: 111,775)\n- GOA annotations: >= 600,000 (target: 613,613)\n- BLAST GOA matches: >= 40,000\n- Total predictions: >= 150M (target: 250M+)\n- Unique proteins: >= 180,000 (target: 195,000)"},{"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":"imports-header","metadata":{},"source":"## 2. 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":"## 3. 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":"## 4. 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":"## 5. GOA Parser\n\nParse Gene Ontology Annotation (GOA) database from **EXTERNAL GOA DATA** (critical for v7)."},{"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":"## 6. 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\ndef compute_blast_scores(hits_df: pd.DataFrame, method: str = \"bitscore_weighted\") -> pd.DataFrame:\n    \"\"\"Compute similarity scores from BLAST hits.\"\"\"\n    df = hits_df.copy()\n\n    if 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    )"},{"cell_type":"markdown","id":"baseline-header","metadata":{},"source":"## 7. 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 = False,\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":"## 8. Load Data and Initialize\n\n**CRITICAL**: Load EXTERNAL GOA data, NOT competition training labels."},{"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'\nEXTERNAL_DIR.mkdir(exist_ok=True)\n\n# CAFA-6 data paths\nTRAIN_FASTA = DATA_DIR / 'Train' / 'train_sequences.fasta'\nTEST_FASTA = DATA_DIR / 'Test' / 'testsuperset.fasta'\nGO_OBO = DATA_DIR / 'Train' / 'go-basic.obo'\n\n# CRITICAL: Download EXTERNAL GOA data (not competition training labels)\nGOA_INPUT_PATH = EXTERNAL_DIR / 'goa_experimental.gaf'\n\nprint(\"=\" * 70)\nprint(\"v7 CONFIGURATION - EXTERNAL GOA DATA\")\nprint(\"=\" * 70)\nprint(f\"Train FASTA: {TRAIN_FASTA}\")\nprint(f\"Test FASTA: {TEST_FASTA}\")\nprint(f\"GO OBO: {GO_OBO}\")\nprint(f\"GOA data (EXTERNAL): {GOA_INPUT_PATH}\")\nprint(\"=\" * 70)\n"},{"cell_type":"code","execution_count":null,"id":"load-go","metadata":{},"outputs":[],"source":"# Load GO graph from competition data\ngo_graph = GOGraph(str(GO_OBO))\nprint(f\"✓ Loaded GO graph with {len(go_graph.graph.nodes)} terms\")"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# Download EXTERNAL GOA data from UniProt FTP\n# This is the CRITICAL fix - using external GOA instead of competition labels\nimport urllib.request\nimport gzip\nimport shutil\n\nprint(\"Downloading GOA data from UniProt FTP...\")\nGOA_URL = 'https://ftp.ebi.ac.uk/pub/databases/GO/goa/UNIPROT/goa_uniprot_gcrp.gaf.gz'\nGOA_GZ = EXTERNAL_DIR / 'goa.gaf.gz'\n\n# Download with progress\nurllib.request.urlretrieve(GOA_URL, GOA_GZ)\nprint(f\"Downloaded: {GOA_GZ.stat().st_size / 1024 / 1024:.1f} MB\")\n\n# Decompress\nprint(\"Decompressing GOA data...\")\nwith gzip.open(GOA_GZ, 'rb') as f_in:\n    with open(GOA_INPUT_PATH, 'wb') as f_out:\n        shutil.copyfileobj(f_in, f_out)\n\nprint(f\"GOA data ready: {GOA_INPUT_PATH}\")\nprint(f\"  Size: {GOA_INPUT_PATH.stat().st_size / 1024 / 1024:.1f} MB\")\n"},{"cell_type":"code","execution_count":null,"id":"load-goa","metadata":{},"outputs":[],"source":"# Load EXTERNAL GOA data (experimental evidence only)\n# This is the CRITICAL fix from v6's mistake\ngoa_df = parse_goa(GOA_INPUT_PATH, experimental_only=True)\n\n# VERIFICATION PRINTS - Critical Success Metrics\nprint(\"\\n\" + \"=\" * 70)\nprint(\"EXTERNAL GOA DATA VERIFICATION\")\nprint(\"=\" * 70)\nprint(f\"✓ GOA proteins loaded: {len(goa_df['protein_id'].nunique()):,}\")\nprint(f\"✓ GOA annotations: {len(goa_df):,}\")\nprint(f\"✓ Unique GO terms: {goa_df['go_id'].nunique():,}\")\nprint(f\"\\nTarget metrics:\")\nprint(f\"  - Expected proteins: 111,775\")\nprint(f\"  - Expected annotations: 613,613\")\nprint(\"=\" * 70)\n\n# Verify we have enough coverage\nassert len(goa_df['protein_id'].nunique()) >= 100000, \"GOA protein count too low!\"\nassert len(goa_df) >= 600000, \"GOA annotation count too low!\"\nprint(\"\\n✓ GOA data verification PASSED\")"},{"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 REFERENCE parameters\nbaseline = GOABlastBaseline(\n    go_graph=go_graph,\n    goa_df=goa_df,  # EXTERNAL GOA data\n    blast=blast,\n    propagation_method=\"max\",  # Reference standard (not noisy_or)\n    use_negative_propagation=False,  # Disabled initially for stability\n    negative_threshold=0.05,\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"BASELINE CONFIGURATION\")\nprint(\"=\" * 70)\nprint(f\"✓ Propagation method: max (reference standard)\")\nprint(f\"✓ Negative propagation: False (disabled for stability)\")\nprint(f\"✓ GOA proteins in baseline: {len(baseline.protein_annotations):,}\")\nprint(\"=\" * 70)"},{"cell_type":"markdown","id":"predict-header","metadata":{},"source":"## 9. Run Predictions\n\nGenerate predictions using REFERENCE baseline parameters."},{"cell_type":"code","execution_count":null,"id":"predict","metadata":{},"outputs":[],"source":"# Run predictions with REFERENCE parameters\noutput_path = WORKING_DIR / 'submission.tsv'\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"PREDICTION PARAMETERS (Reference Baseline)\")\nprint(\"=\" * 70)\nprint(f\"threshold: 0.01\")\nprint(f\"max_terms_per_protein: 1500\")\nprint(f\"n_hits: 50\")\nprint(f\"min_identity: 30.0\")\nprint(f\"max_evalue: 1e-3\")\nprint(f\"score_method: bitscore_weighted\")\nprint(\"=\" * 70 + \"\\n\")\n\ntotal_preds, num_proteins = baseline.predict_to_tsv(\n    str(TEST_FASTA),\n    str(output_path),\n    threshold=0.01,              # Reference threshold\n    max_terms_per_protein=1500,\n    batch_size=5000,\n    n_hits=50,                   # Reference standard\n    min_identity=30.0,           # Reference minimum\n    max_evalue=1e-3,             # Reference threshold\n    score_method=\"bitscore_weighted\",  # Simple, proven scoring\n)\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"PREDICTION RESULTS - SUCCESS METRICS\")\nprint(\"=\" * 70)\nprint(f\"✓ Total predictions: {total_preds:,}\")\nprint(f\"✓ Unique proteins with predictions: {num_proteins:,}\")\nprint(f\"\\nTarget metrics:\")\nprint(f\"  - Minimum predictions: 150M\")\nprint(f\"  - Target predictions: 250M+\")\nprint(f\"  - Minimum proteins: 180K\")\nprint(f\"  - Target proteins: 195K\")\nprint(\"=\" * 70)"},{"cell_type":"markdown","id":"stats-header","metadata":{},"source":"## 10. Verify Submission File"},{"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'])\n\nprint(\"\\n\" + \"=\" * 70)\nprint(\"SUBMISSION FILE STATISTICS\")\nprint(\"=\" * 70)\nprint(f\"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())\nprint(\"=\" * 70)"},{"cell_type":"code","execution_count":null,"id":"submission","metadata":{},"outputs":[],"source":"# Show sample of submission\nprint(\"\\nSample 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(\"\\n✓ Submission file ready at /kaggle/working/submission.tsv\")"},{"cell_type":"code","execution_count":null,"id":"verify","metadata":{},"outputs":[],"source":"# Verify submission file format\n!head -20 /kaggle/working/submission.tsv\n!echo \"\"\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## Expected Results\n\n**Expected Score**: 0.32-0.36\n\n**Success Criteria**:\n- **Minimum**: Score >= 0.32 (validates approach, recovers to v3 level)\n- **Target**: Score >= 0.35 (approaches reference baselines)\n- **Ideal**: Score >= 0.36 (matches/exceeds reference baselines)\n\n## Key Changes from v6 (The Failed Version)\n\n1. **External GOA Data**: Reverted to external GOA (111K proteins, 613K annotations) instead of competition training labels (82K proteins, 537K annotations)\n2. **Reference Parameters**: Using proven parameters from reference baselines (propagation_method=\"max\", n_hits=50, etc.)\n3. **Verification Prints**: Added comprehensive verification to confirm correct data loading\n\n## Next Steps\n\n1. **Submit to competition**\n2. **Check leaderboard score**\n3. **If score >= 0.32**: Proceed with v7.1 (enable negative propagation)\n4. **If score >= 0.35**: Proceed with parameter optimization and ensemble methods"}],"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.7.12"}},"nbformat":4,"nbformat_minor":5}