From d476bca527cb12c3313a77352203f90d34636047 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Fri, 16 Oct 2026 18:35:10 -0700 Subject: [PATCH] Link metabolite evidence to genes through KEGG reactions (metabolite_linked_convergence tier) Multi-omics convergence linked layers only by identical feature_id_standardized, so metabolite evidence (standardized to the metabolite name) could never converge with gene-level evidence. Add a second, explicitly weaker kind of convergence: metabolite name -> KEGG compound -> reaction -> KO -> crg GeneID. - src/mappings/metabolite_links.py: per-partition link tables (simulated demo vs real, never merged), currency-metabolite and hub (>25 genes) exclusions kept visible, link confidence never exact / one-to-one. - src/mappings/build_metabolite_links.py + `aree build-metabolite-links`: queries KEGG REST only for curated compounds, caches every response with sha256, resolves GeneIDs against the RS_2024_06 crosswalk and retired-ID sidecar. Output is git-ignored (KEGG redistribution terms). - rank.py: new tier metabolite_linked_convergence between multi_omics_convergence and emerging; linked layers never enter n_supporting_layers or the score; new convergence_basis / linked_* / link_table_* candidate fields. - Evidence cards: "Metabolite-gene links" section with the compound / reaction / KO chain, link confidence, table checksum, stressor caveat and exclusions. - Demo: fabricated ldh (significant) and sdha (not) in the demo crosswalk and GIGAS_HEAT01 chronic table, so GIGAS_GROW06 lactate <-> ldh converges by link and succinate/sdha and ATP do not. - Design note docs/metabolite_gene_links.md, provenance JSON schema, docs, app legend, dashboard KPI, Makefile target, tests. Co-Authored-By: Claude Opus 5.5 --- .gitignore | 7 + Makefile | 6 +- README.md | 1 + app/main.py | 5 + ...EAT01_chronic_heat_vs_control_dge_demo.tsv | 2 + data/mappings/currency_metabolites.yaml | 23 ++ data/mappings/gene_id_crosswalk.tsv | 2 + data/mappings/metabolite_compound_map.tsv | 19 + ...metabolite_gene_links_demo.provenance.json | 23 ++ data/mappings/metabolite_gene_links_demo.tsv | 9 + data/reference/genome_assemblies.yaml | 9 + docs/_quarto.yml | 2 + docs/identifier_mapping.md | 11 + docs/implementation_status.md | 1 + docs/index.qmd | 5 +- docs/interpreting_candidate_scores.md | 17 +- docs/metabolite_gene_links.md | 219 ++++++++++ docs/methods.md | 13 +- ...tabolite_gene_links_provenance.schema.json | 79 ++++ src/aree/cli.py | 35 ++ src/mappings/build_metabolite_links.py | 291 +++++++++++++ src/mappings/metabolite_links.py | 300 ++++++++++++++ src/prioritize/rank.py | 94 ++++- src/reporting/dashboard.py | 1 + src/reporting/evidence_cards.py | 136 +++++- tests/test_metabolite_links.py | 388 ++++++++++++++++++ 26 files changed, 1690 insertions(+), 8 deletions(-) create mode 100644 data/mappings/currency_metabolites.yaml create mode 100644 data/mappings/metabolite_compound_map.tsv create mode 100644 data/mappings/metabolite_gene_links_demo.provenance.json create mode 100644 data/mappings/metabolite_gene_links_demo.tsv create mode 100644 docs/metabolite_gene_links.md create mode 100644 schemas/metabolite_gene_links_provenance.schema.json create mode 100644 src/mappings/build_metabolite_links.py create mode 100644 src/mappings/metabolite_links.py create mode 100644 tests/test_metabolite_links.py diff --git a/.gitignore b/.gitignore index 8a9cf3e..93f3fcd 100644 --- a/.gitignore +++ b/.gitignore @@ -58,3 +58,10 @@ registry/study_registry.csv # OrthoFinder intermediates (staged proteomes, DIAMOND results); the run's # outputs under data/reference/orthogroups// ARE committed. data/reference/orthogroups/_work/ + +# KEGG-derived metabolite -> gene link table for real C. gigas evidence. Built +# locally by `aree build-metabolite-links`; NOT committed, because KEGG does not +# clearly grant redistribution of derived tables (docs/metabolite_gene_links.md). +# The simulated demo table (metabolite_gene_links_demo.*) IS committed. +data/mappings/metabolite_gene_links_mgigas.* +data/mappings/_sources/ diff --git a/Makefile b/Makefile index 6eeaf46..b6b1a80 100644 --- a/Makefile +++ b/Makefile @@ -1,5 +1,5 @@ # AREE task runner. Run `make help` for the list. -.PHONY: help registry install test lint demo demo-clean register harmonize meta cards crosswalk intake intake-check real-study real-pool dashboard docs app clean +.PHONY: help registry install test lint demo demo-clean register harmonize meta cards crosswalk intake intake-check real-study real-pool dashboard docs app clean metabolite-links # Every simulated study in registry/studies/, so a new demo study needs no edit here. DEMO_STUDIES = $(shell aree study-ids --simulated) @@ -12,6 +12,7 @@ help: @echo " lint run ruff over src and tests" @echo " demo run the full demo pipeline (register -> harmonize -> meta -> cards)" @echo " demo-clean remove generated reports/ and the registry index, then run demo" + @echo " metabolite-links build the real C. gigas metabolite->gene link table from KEGG (git-ignored)" @echo " crosswalk rebuild the real NCBI/UniProt identifier crosswalk (~380 MB download)" @echo " intake regenerate the real study's result tables from the published source" @echo " intake-check verify the committed result tables still reproduce (no writes)" @@ -62,6 +63,9 @@ demo-clean: clean demo crosswalk: aree build-crosswalk +metabolite-links: + aree build-metabolite-links + REAL_CROSSWALK := data/reference/crosswalk/mgigas_gene_id_crosswalk.tsv REAL_INTAKE := data/studies/HESSER2024_VCOR/intake.yaml diff --git a/README.md b/README.md index e3a1f5b..94de269 100644 --- a/README.md +++ b/README.md @@ -151,6 +151,7 @@ AREE/ | `aree build-evidence-cards [--phenotype

] [--feature-type ] [--max-adjusted-p ] [--all-cards]` | Rank every candidate into `reports/evidence_cards/candidates.tsv` and write an evidence card for each one with a BH-adjusted p ≤ `q` (default 0.05, pooled or in any contributing study) | | `aree top-candidates [--n ] [--candidates ] [--out ]` | Regroup an existing `candidates.tsv` by phenotype and write the top `N` (default 10) per phenotype, ranked by tier then score, to `reports/top_candidates_summary.md` | | `aree build-crosswalk --taxid [--slug ]` | Build a real identifier crosswalk from NCBI Gene + UniProtKB for a registered species | +| `aree build-metabolite-links [--offline]` | Build the real *C. gigas* metabolite → enzyme gene link table from the KEGG REST API (curated compounds only; output is git-ignored, not redistributed). See [docs/metabolite_gene_links.md](docs/metabolite_gene_links.md) | | `aree fetch-reference --assembly [--genome] [--proteins] [--plan] [--check]` | Fetch, verify and record a RefSeq reference bundle (transcriptome, GTF, tx2gene, optional genome, optional one-protein-per-gene proteome) | | `aree build-orthogroups --orthofinder [--threads N] [--activate]` | Infer orthogroups across every species with a crosswalk (OrthoFinder on RefSeq proteomes); writes `data/reference/orthogroups//` | | `aree build-dashboard [--out ]` | Summarize registry, pipeline, and findings state into `docs/dashboard/data.json` (plus the top candidates' evidence cards) for the Quarto dashboard page | diff --git a/app/main.py b/app/main.py index 304821f..cf8d332 100644 --- a/app/main.py +++ b/app/main.py @@ -149,6 +149,11 @@ "- **multi_omics_convergence** — ≥2 molecular layers (e.g. transcriptomics + DNA methylation) " "each with a significant record for the same standardized feature under the **same phenotype**, " "including the candidate's own layer.\n" + "- **metabolite_linked_convergence** — the candidate's own layer is significant and a " + "metabolite and a gene are linked only through KEGG reaction annotation " + "(metabolite → reaction → KO → gene), under the same phenotype. Weaker than shared gene " + "identity, ranked below it, and not counted in the score " + "(see docs/metabolite_gene_links.md).\n" "- **emerging** — single-study or otherwise unconfirmed; **requires replication**." ) diff --git a/data/demo/rnaseq/GIGAS_HEAT01_chronic_heat_vs_control_dge_demo.tsv b/data/demo/rnaseq/GIGAS_HEAT01_chronic_heat_vs_control_dge_demo.tsv index 8388108..b445075 100644 --- a/data/demo/rnaseq/GIGAS_HEAT01_chronic_heat_vs_control_dge_demo.tsv +++ b/data/demo/rnaseq/GIGAS_HEAT01_chronic_heat_vs_control_dge_demo.tsv @@ -7,3 +7,5 @@ LOC105336010 410.8 -1.40 0.28 -5.00 5.7e-07 0.0050 LOC105336200 325.9 0.80 0.33 2.42 0.0155 0.0600 LOC105331241 910.2 0.30 0.30 1.00 0.3170 0.4000 LOC105329871 701.3 0.20 0.31 0.65 0.5160 0.6000 +LOC105342001 1520.3 1.05 0.30 3.50 0.0005 0.0100 +LOC105342002 840.7 0.25 0.31 0.81 0.4200 0.5200 diff --git a/data/mappings/currency_metabolites.yaml b/data/mappings/currency_metabolites.yaml new file mode 100644 index 0000000..29addca --- /dev/null +++ b/data/mappings/currency_metabolites.yaml @@ -0,0 +1,23 @@ +# Currency (cofactor / inorganic) compounds that are never used as a metabolite -> gene link. +# +# Each takes part in hundreds of KEGG reactions (ATP, C00002, in 810 on 2026-10-16), +# so a change in one would "link" to most of the metabolic genome. They stay in the +# evidence table and on evidence cards as metabolites; they are only excluded as LINKS, +# and every exclusion is shown on the card. See docs/metabolite_gene_links.md. +version: "0.1.0" +compounds: + - {kegg_compound: C00001, name: H2O, reason: "solvent; participates in most hydrolytic reactions"} + - {kegg_compound: C00080, name: "H+", reason: "proton; balancing species"} + - {kegg_compound: C00002, name: ATP, reason: "energy currency; ~800 KEGG reactions"} + - {kegg_compound: C00008, name: ADP, reason: "energy currency"} + - {kegg_compound: C00020, name: AMP, reason: "energy currency"} + - {kegg_compound: C00003, name: "NAD+", reason: "redox cofactor"} + - {kegg_compound: C00004, name: NADH, reason: "redox cofactor"} + - {kegg_compound: C00006, name: "NADP+", reason: "redox cofactor"} + - {kegg_compound: C00005, name: NADPH, reason: "redox cofactor"} + - {kegg_compound: C00009, name: Orthophosphate, reason: "inorganic cofactor"} + - {kegg_compound: C00013, name: Diphosphate, reason: "inorganic cofactor"} + - {kegg_compound: C00010, name: CoA, reason: "acyl-carrier cofactor"} + - {kegg_compound: C00011, name: CO2, reason: "inorganic; carboxylation/decarboxylation"} + - {kegg_compound: C00007, name: Oxygen, reason: "inorganic; oxidase/oxygenase co-substrate"} + - {kegg_compound: C00014, name: Ammonia, reason: "inorganic nitrogen"} diff --git a/data/mappings/gene_id_crosswalk.tsv b/data/mappings/gene_id_crosswalk.tsv index db8607a..3f38563 100644 --- a/data/mappings/gene_id_crosswalk.tsv +++ b/data/mappings/gene_id_crosswalk.tsv @@ -38,3 +38,5 @@ LOC105335777 ENSCGRG00000018 K1QQ18 LOC105335777 novel2 OG0000018 LOC105340010 ENSCGRG00000019 K1QQ10 LOC105340010 aqp_like OG0000019 LOC105340020 ENSCGRG00000020 K1QQ20 LOC105340020 atp1a_nka OG0000020 LOC105341000 ENSCGRG00000021 K1QQ21 LOC105341000 actin OG0000021 +LOC105342001 ENSCGRG00000022 K1QQ22 LOC105342001 ldh OG0000022 +LOC105342002 ENSCGRG00000023 K1QQ23 LOC105342002 sdha OG0000023 diff --git a/data/mappings/metabolite_compound_map.tsv b/data/mappings/metabolite_compound_map.tsv new file mode 100644 index 0000000..41c50e2 --- /dev/null +++ b/data/mappings/metabolite_compound_map.tsv @@ -0,0 +1,19 @@ +# Curated metabolite name -> KEGG compound / ChEBI map. Used by src/mappings/metabolite_links.py. +# One row per (name, compound). A reported name that cannot be resolved to a single stereoisomer +# gets one row per candidate compound and match_basis=curator_judgement (links become 'inferred'). +# KEGG compound IDs and ChEBI cross-references checked against https://rest.kegg.jp/get/cpd: on 2026-10-16. +# Names are matched after normalization (lowercase; spaces/hyphens -> underscores). See docs/metabolite_gene_links.md. +metabolite_name kegg_compound chebi_id match_basis curator date_curated note +lactate C00186 CHEBI:16651 curator_judgement aree_curator 2026-10-16 LC-MS does not resolve stereoisomers; (S)-lactate. Molluscs carry both L- and D-LDH, so both are linked. +lactate C00256 CHEBI:16004 curator_judgement aree_curator 2026-10-16 LC-MS does not resolve stereoisomers; (R)-lactate (D-lactate). +l_lactate C00186 CHEBI:16651 exact_name aree_curator 2026-10-16 +d_lactate C00256 CHEBI:16004 exact_name aree_curator 2026-10-16 +succinate C00042 CHEBI:15741 exact_name aree_curator 2026-10-16 +glucose C00031 CHEBI:4167 synonym aree_curator 2026-10-16 D-glucose assumed (the biological enantiomer); anomer not resolved. +d_glucose C00031 CHEBI:4167 exact_name aree_curator 2026-10-16 +atp C00002 CHEBI:15422 exact_name aree_curator 2026-10-16 Currency metabolite; never credited as a link. +alanine C00041 CHEBI:16977 synonym aree_curator 2026-10-16 L-alanine assumed; D-alanine is minor in bivalve tissue but not excluded by the assay. +l_alanine C00041 CHEBI:16977 exact_name aree_curator 2026-10-16 +glycine_betaine C00719 CHEBI:17750 synonym aree_curator 2026-10-16 KEGG name 'Betaine'. +betaine C00719 CHEBI:17750 exact_name aree_curator 2026-10-16 +pyruvate C00022 CHEBI:15361 exact_name aree_curator 2026-10-16 diff --git a/data/mappings/metabolite_gene_links_demo.provenance.json b/data/mappings/metabolite_gene_links_demo.provenance.json new file mode 100644 index 0000000..e1ea420 --- /dev/null +++ b/data/mappings/metabolite_gene_links_demo.provenance.json @@ -0,0 +1,23 @@ +{ + "artifact": "metabolite_gene_links_demo.tsv", + "artifact_sha256": "2156a4428d9a46ea18f979cccff155ed280c5ca666414f2143c9181b9850d3f0", + "simulated": true, + "builder": "hand-written demo fixture", + "date_generated": "2026-10-16", + "ncbi_taxid": 29159, + "organism": "Magallana gigas (Crassostrea gigas) \u2014 SIMULATED demo partition only", + "source_summary": "SIMULATED demo links; real KEGG reaction/KO IDs, fabricated gene IDs", + "reference_annotation": "cgigas_annotation_v1_demo (fabricated demo crosswalk, data/mappings/gene_id_crosswalk.tsv)", + "sources": [ + { + "name": "KEGG REST (identifiers only, checked by hand)", + "url": "https://rest.kegg.jp/get/R00703+R02164", + "accessed": "2026-10-16" + } + ], + "max_genes_per_compound": 25, + "known_limitations": [ + "Simulated: gene IDs are fabricated demo LOC numbers that collide with real NCBI GeneIDs denoting different genes.", + "Used only for evidence with simulated=true and species_taxid=29159." + ] +} diff --git a/data/mappings/metabolite_gene_links_demo.tsv b/data/mappings/metabolite_gene_links_demo.tsv new file mode 100644 index 0000000..e4a1ffd --- /dev/null +++ b/data/mappings/metabolite_gene_links_demo.tsv @@ -0,0 +1,9 @@ +# SIMULATED DEMO metabolite -> gene link table. The gene IDs are the FABRICATED demo LOC numbers +# from data/mappings/gene_id_crosswalk.tsv and must never be read as real C. gigas genes. The KEGG +# compound, reaction and KO identifiers are real KEGG identifiers (checked 2026-10-16); only the +# KO -> gene assignment is fabricated. Never merge with a real table (metabolite_gene_links_mgigas.tsv). +# Columns are the same contract `aree build-metabolite-links` writes. See docs/metabolite_gene_links.md. +kegg_compound compound_name kegg_reaction equation_side ko ncbi_gene_id link_type +C00186 (S)-Lactate R00703 left K00016 LOC105342001 enzyme +C00022 Pyruvate R00703 right K00016 LOC105342001 enzyme +C00042 Succinate R02164 left K00234 LOC105342002 enzyme diff --git a/data/reference/genome_assemblies.yaml b/data/reference/genome_assemblies.yaml index 937dce2..4b8f82a 100644 --- a/data/reference/genome_assemblies.yaml +++ b/data/reference/genome_assemblies.yaml @@ -5,6 +5,13 @@ # that, so a typo or an undocumented assembly is caught at registration rather # than surfacing later as evidence that cannot be compared. # +# `metabolite_gene_links` (optional) is the metabolite -> enzyme gene link table +# that metabolomics evidence on this species is linked through for the +# `metabolite_linked_convergence` tier. It is built locally by +# `aree build-metabolite-links` and NOT committed (KEGG terms); until it is built, +# real candidates report `link_table_status: not_built`. See +# docs/metabolite_gene_links.md. +# # `identifier_crosswalk` is the crosswalk that evidence from studies on this # assembly is standardized against (repo-relative path, or null if none has been # built yet). Harmonizing a real study reads it from here, so a mussel study and @@ -34,6 +41,7 @@ assemblies: annotation_release: "GCF_963853765.1-RS_2024_06" annotation_release_date: "2024-06-20" identifier_crosswalk: "data/reference/crosswalk/mgigas_gene_id_crosswalk.tsv" + metabolite_gene_links: "data/mappings/metabolite_gene_links_mgigas.tsv" verified: "2026-08-28" notes: >- Current NCBI reference genome for the Pacific oyster. AREE's identifier @@ -52,6 +60,7 @@ assemblies: is_ncbi_reference: false annotation_release: null identifier_crosswalk: "data/reference/crosswalk/mgigas_gene_id_crosswalk.tsv" + metabolite_gene_links: "data/mappings/metabolite_gene_links_mgigas.tsv" verified: "2026-08-28" notes: >- The assembly most C. gigas studies published before ~2024 were mapped diff --git a/docs/_quarto.yml b/docs/_quarto.yml index 60eb0af..d717201 100644 --- a/docs/_quarto.yml +++ b/docs/_quarto.yml @@ -55,6 +55,8 @@ website: href: raw_vs_processed.md - text: "Identifier mapping" href: identifier_mapping.md + - text: "Metabolite–gene links" + href: metabolite_gene_links.md - section: "Interpreting results" contents: - text: "Meta-analysis outputs" diff --git a/docs/identifier_mapping.md b/docs/identifier_mapping.md index 1b7e1f6..ec790d5 100644 --- a/docs/identifier_mapping.md +++ b/docs/identifier_mapping.md @@ -361,6 +361,17 @@ annotations. (`src/prioritize/scoring.py`, `worst_mapping = min(...)`), so a candidate is only as trustworthy, identifier-wise, as its weakest link. +## Metabolites are linked, not mapped + +Metabolite evidence keeps the reported metabolite name as its +`feature_id_standardized`; it is never forced into the gene hierarchy above. +Its connection to genes is a separate, weaker relationship (metabolite → +KEGG compound → reaction → KO → crg GeneID) resolved at ranking time from +`data/mappings/metabolite_gene_links_*.tsv`, chosen per partition exactly as +the crosswalk is: the simulated demo table for simulated evidence, the table +registered under `metabolite_gene_links` in `genome_assemblies.yaml` for real +evidence. See [metabolite_gene_links.md](metabolite_gene_links.md). + ## Related documentation - [interpreting_candidate_scores.md](interpreting_candidate_scores.md) diff --git a/docs/implementation_status.md b/docs/implementation_status.md index d9df2d3..a4962f3 100644 --- a/docs/implementation_status.md +++ b/docs/implementation_status.md @@ -22,6 +22,7 @@ demo pipeline runs from the committed header-only registry). | Provenance manifests | `src/harmonize/core.py` (`_write_harmonize_manifest`) | checksum, params, versions, warnings per comparison | | Random-effects meta-analysis | `src/meta_analysis/` | DerSimonian–Laird, I²/τ², direction-consistency; Benjamini–Hochberg `adjusted_p_value` per phenotype/feature-type/origin/species family; tail p-values via the survival function (no underflow to 0); standard errors derived once for the whole table and groups sliced as flat arrays, so a 54k-record, 30k-feature pool runs in about 5 s (was ~95 s); poolability-aware replication counts; lower-confidence alias de-duplication; fails closed on repeated within-study contrasts unless the study prespecifies one comparison as `meta_analysis_primary`, with the others counted in `n_excluded_non_primary` | | Transparent candidate scoring + tiers | `src/prioritize/` | named weighted components scored on adjusted p, each unbounded component saturating rather than clipping so the top of a genome-wide ranking does not tie (126 candidates shared one score before 2026-09-05; the largest tie group is now 10); every candidate also carries `evidence_class` and `context_replication`, which say what the evidence is *of* and whether it replicates across biological contexts, independent of tier; hard gates that a score cannot override, including adjusted p ≤ 0.05 for the top tier; multi-omics convergence counts molecular layers (from the feature_types vocabulary) with significant same-phenotype support, not feature types anywhere; origin/species partitions preserved; tiers ordered strongest-first | +| Metabolite → gene links (`metabolite_linked_convergence` tier) | `src/mappings/metabolite_links.py`, `src/mappings/build_metabolite_links.py`, `data/mappings/metabolite_compound_map.tsv`, `data/mappings/currency_metabolites.yaml` | metabolite name → KEGG compound → reaction → KO → crg GeneID; currency and hub (> 25 genes) compounds shown but never credited; links never `exact`; separate tier below shared-identity convergence, no score contribution; demo GIGAS_GROW06 lactate ↔ simulated `ldh` converges, `sdha`/succinate and ATP do not. The real KEGG-derived table is built locally (`aree build-metabolite-links`, 191 links on 2026-10-16) and **not committed** (KEGG terms), so CI exercises the simulated table only; no real metabolomics study has passed QC yet (DUPERRET2025_POMS_STARVE failed), so no real candidate is linked | | Evidence cards | `src/reporting/evidence_cards.py` | every card opens with a "What this evidence does and does not show" section naming the evidence class, the contexts that do and do not vary, and — where no contributing study measured a resilience outcome — the phenotype gap ahead of any validation advice; every candidate incl. emerging ranked into `candidates.tsv`; a markdown card for each with a significant signal (BH-adjusted p ≤ 0.05 pooled or in a contributing study; `--all-cards` for all); single indexed pass, so a 30k-candidate real pool renders in about a minute; partition-safe filenames, forest/multi-omics context, limitations and next step | | `aree` CLI | `src/aree/cli.py` | all 6 required commands run on a clean machine, plus `build-crosswalk` and `intake-supplementary` | | Streamlit interface | `app/main.py` | browse/filter studies, evidence, candidates, search, CSV/TSV export; verified serving HTTP 200 | diff --git a/docs/index.qmd b/docs/index.qmd index 2d5d193..ff2b195 100644 --- a/docs/index.qmd +++ b/docs/index.qmd @@ -26,7 +26,7 @@ fmt = (n) => (n == null ? "—" : n.toLocaleString("en-US")) fmtP = (p) => (p == null ? "—" : p < 1e-3 ? p.toExponential(1) : p.toFixed(3)) phen = (id) => data.vocab.phenotype_labels[id] ?? id stress = (id) => data.vocab.stressor_labels[id] ?? id -TIER = new Map([["high_priority_cross_study", "High priority"], ["multi_omics_convergence", "Multi-omics"], ["emerging", "Emerging"]]) +TIER = new Map([["high_priority_cross_study", "High priority"], ["multi_omics_convergence", "Multi-omics"], ["metabolite_linked_convergence", "Metabolite-linked"], ["emerging", "Emerging"]]) clip = (v, n) => (v && v.length > n ? htl.html`${v.slice(0, n - 1).trimEnd()}…` : v ?? "") // "Calla B, Thompson NF, Burge CA (2026)" -> "Calla et al. (2026)" shortCite = (c) => { @@ -129,7 +129,8 @@ html`

const tierText = html`

${fmt(nRanked)} gene × phenotype candidates are ranked: ${fmt(k.n_candidates_real_high_priority)} high priority (replicated across studies), - ${fmt(k.n_candidates_real_multi_omics)} multi-omics, ${fmt(k.n_candidates_real_emerging)} emerging. + ${fmt(k.n_candidates_real_multi_omics)} multi-omics, ${fmt(k.n_candidates_real_metabolite_linked)} metabolite-linked, + ${fmt(k.n_candidates_real_emerging)} emerging. By evidence class: ${classes.map(d => `${fmt(d.n)} ${d.evidence_class.replace(/_/g, " ")}`).join(", ")}.

`; const caveats = []; diff --git a/docs/interpreting_candidate_scores.md b/docs/interpreting_candidate_scores.md index 13c5207..a826f4a 100644 --- a/docs/interpreting_candidate_scores.md +++ b/docs/interpreting_candidate_scores.md @@ -123,7 +123,22 @@ Tier assignment, in order: qualifies, because every cross-layer overlap in the demo is cross-phenotype. The card still lists that adjacent evidence and marks each row with whether it counts. -3. **`emerging`** — everything else: single-study support, or support that +3. **`metabolite_linked_convergence`** — the candidate's own layer is + significant, and another layer carries a significant same-phenotype record + for a feature **linked** to it biochemically rather than by a shared ID: a + metabolite and a gene whose product KEGG annotates (by orthology) to + catalyse a reaction involving that metabolite. Checked only if neither tier + above applies. This is a deliberately weaker, separately named kind of + convergence: it does not add to `n_supporting_layers` or to the score; + currency metabolites (ATP, NAD⁺, …) and hub compounds reaching more than + 25 genes are never credited; and every link is `inferred` or + `one_to_many_ortholog`, never `exact`. The card lists each link's + compound, reaction, KO, confidence and source checksum, plus the + exclusions. Candidates carry `convergence_basis` + (`shared_identifier` / `metabolite_gene_link` / `none`), `linked_layers`, + `linked_features`, `link_confidences` and `link_table_status`. See + [metabolite_gene_links.md](metabolite_gene_links.md). +4. **`emerging`** — everything else: single-study support, or support that fails one or more of the high-priority gates. Always labeled as requiring replication. diff --git a/docs/metabolite_gene_links.md b/docs/metabolite_gene_links.md new file mode 100644 index 0000000..68afcbd --- /dev/null +++ b/docs/metabolite_gene_links.md @@ -0,0 +1,219 @@ +# Design note: linking metabolite evidence to genes + +Status: **implemented (MVP).** Decisions in section 7 were confirmed on +2026-10-16. + +## 1. Problem + +The `multi_omics_convergence` tier (`src/prioritize/rank.py`, +`EvidenceIndex.supporting_layers` / `other_layers`) links molecular layers only +by an identical `feature_id_standardized`. RNA-seq, methylation and proteomics +evidence all standardize to an NCBI GeneID, so they can meet. Metabolite +evidence (`src/harmonize/metabolomics.py`) standardizes to the putative +metabolite *name* (`lactate`, `succinate`, ...), so it never shares an ID with +gene-level evidence. As a result, metabolomics can never take part in +convergence, whatever the data show. + +Metabolites are not gene products, so a shared ID would be the wrong fix. The +honest link is biochemical: *this metabolite is a substrate or product of a +reaction that this gene's product catalyses*. That is a weaker claim than "the +same gene was measured twice". It must be labelled weaker, scored as weaker, +and explained on the card. + +## 2. The link chain + +``` +metabolite name (as reported, preserved) + └─ curated name → compound (KEGG C-number, ChEBI ID) [curation step] + └─ compound → reaction (KEGG REACTION, substrate/product side) + └─ reaction → KO (KEGG Orthology) + └─ KO → crg gene (KEGG GENES, organism code `crg`) + └─ crg gene = NCBI GeneID → AREE crosswalk / orthogroup +``` + +Facts checked on 2026-10-16 against the KEGG REST API: + +- `crg` (T03920, *Magallana gigas*) has 26,085 proteins, and 10,575 of them + have a KO. `link/ko/crg` returns 10,878 gene→KO pairs. +- `crg` gene IDs are NCBI GeneIDs. **10,874 of the 10,878 are current GeneIDs + in AREE's RS_2024_06 crosswalk** (`data/reference/crosswalk/mgigas_gene_id_crosswalk.tsv`). + None are retired, and 4 are unknown. The builder still routes every ID through + the retired-ID sidecar and records the 4 that don't resolve. It does not drop + them silently. +- Currency metabolites fan out enormously. For example, ATP (C00002) appears in + 810 KEGG reactions. Without a hub filter, ATP alone would "link" most of the + metabolic genome to any ATP change. + +## 3. Proposed files (`data/mappings/`) + +| File | Contents | +|---|---| +| `metabolite_compound_map.tsv` | Curated `metabolite_name` (normalized, lowercase), `synonyms`, `kegg_compound`, `chebi_id`, `match_basis` (`exact_name` / `synonym` / `curator_judgement`), `curator`, `date`, `note`. Hand-curated and small. It holds only the names that appear in registered studies. | +| `currency_metabolites.yaml` | The explicit list of compounds that are never used as a link (ATP, ADP, AMP, NAD⁺/NADH, NADP⁺/NADPH, H₂O, H⁺, CO₂, O₂, Pi, PPi, CoA, NH₃), with the reason. These metabolites still appear as evidence. They just don't link. | +| `metabolite_gene_links_mgigas.tsv` | Built locally, **not committed** (decision A), one row per (compound, reaction, KO, gene): `kegg_compound`, `kegg_reaction`, `equation_side` (`left`/`right`), `ko`, `ncbi_gene_id`, `link_type`. The per-compound gene count and `link_confidence` are derived when the table is loaded, so they cannot drift from the rows. | +| `metabolite_gene_links_mgigas.provenance.json` | Same shape as the crosswalk provenance. It records the KEGG release (`info/kegg`), the access date, every REST URL fetched with its sha256, the builder version, the sha256 of the output, coverage counts, unresolved GeneIDs, and limitations. | +| `metabolite_gene_links_demo.tsv` (+ provenance) | **Simulated**, against the fabricated demo LOC IDs. It obeys the same rule as the demo crosswalk: it is never merged with the real table. | + +A builder, `src/mappings/build_metabolite_links.py`, fetches **only** the +compounds listed in `metabolite_compound_map.tsv`. It makes one +`link/reaction/` call per compound, then the reaction→KO and +reaction-equation lookups, plus a single `link/ko/crg` call. It never mirrors +KEGG. `genome_assemblies.yaml` gains a `metabolite_gene_links:` key next to +`identifier_crosswalk:`, so a real study resolves links through its assembly +and a simulated one through the demo table, exactly as identifiers do today. + +**No change to the evidence schema.** Metabolite records keep +`feature_id_original` (the source feature ID) and `feature_id_standardized` +(the name). The link is resolved at ranking time from the table, so the link +table can be rebuilt without re-harmonizing any study. + +## 4. Mapping confidence + +The existing `mapping_confidence` vocabulary is reused for `link_confidence`. +A metabolite→gene link is **never** `exact` or `one_to_one_ortholog`: + +| Situation | `link_confidence` | +|---|---| +| Compound reaches exactly one crg gene, through a KO transferred from KEGG's ortholog assignment | `inferred` | +| Compound reaches 2 or more crg genes (the usual case: paralogs, or several reactions) | `one_to_many_ortholog` | +| Name→compound match was `curator_judgement`, not an exact name or synonym | downgraded to `inferred` at most, and flagged | +| Compound reaches more than `MAX_GENES_PER_COMPOUND` genes (proposed: 25) | link **not credited**. Recorded as `hub_compound` | +| Compound in `currency_metabolites.yaml` | link **not credited** | + +Only metabolites at MSI annotation level 1–2 have a standardized name today +(level 3–4 → `unresolved`), so only those can link. That gate stays as it is. + +`link_type` is `enzyme` for MVP. Transporter links (for example, KEGG BRITE +ko02000 / TCDB) are listed as future work. KEGG REACTION rarely encodes +transport, and inventing a transporter link would overstate what is known. + +## 5. Ranking changes (`src/prioritize/rank.py`) + +A second, **separately named** kind of convergence. It never merges into the +first. + +- `EvidenceIndex` loads the partition's link table and builds both directions: + gene → linked metabolite names, and metabolite name → linked GeneIDs. +- New per-candidate fields: + - `linked_layers`: layers with a significant record (adj. p ≤ 0.05) under the + **same phenotype and partition** for a *linked* feature, excluding layers + already in `supporting_layers`. + - `n_linked_layers`, `linked_features` (for example, + `lactate→LOC… via R00703/K00016`), `link_confidences`, + `link_table_sha256`. + - `convergence_basis`: `shared_identifier`, `metabolite_gene_link`, or `none`. +- New tier `metabolite_linked_convergence`. It is ordered **after** + `multi_omics_convergence` and **before** `emerging`. The gate is the + candidate's own layer is significant AND at least one linked layer is + significant. High-priority and shared-identity gates are evaluated first, so + a link never promotes a candidate past a stronger tier and never demotes one. +- **Score:** linked layers do **not** enter `n_supporting_layers` or + `assay_diversity_score`. The transparent score stays a function of + shared-identity evidence only. The tier and the card carry the link. (See + decision B for the alternative.) +- Direction across the link is reported but not required to agree. That is the + same policy as shared-identity convergence. A product rising while its enzyme + falls has several valid explanations. + +## 6. Evidence card + +There is a new subsection under "Multi-omics context": +**"Metabolite–gene links (biochemical, weaker than shared gene identity)"**. +It contains: + +- A plain statement: "Linked through KEGG reaction annotation, not measured. + The enzyme assignment is inferred by KEGG orthology. The link means the gene's + product is annotated to catalyse a reaction involving this metabolite. It does + not mean the gene regulates or controls the metabolite's level." +- A table with study, layer, linked feature, compound, reaction, KO, side, + `link_confidence`, n genes for the compound, direction, adj. p, same + phenotype, significant, and counts. +- The source line: link table file, KEGG release, access date, sha256. +- Excluded links listed with their reason (`currency_metabolite`, + `hub_compound`), so the filter is visible rather than silent. +- A stressor caveat when the linked records come from a different stressor + (for example, a hypoxia metabolite and a heat transcript under + `growth_under_stress`). + +## 7. Decisions (confirmed 2026-10-16) + +**A. Source and licensing: builder only, no committed KEGG table.** KEGG +REST is provided for academic use, and KEGG does not clearly grant +redistribution of derived tables. AREE therefore commits the builder +(`aree build-metabolite-links`), the curated compound map, the currency list +and the provenance schema (`schemas/metabolite_gene_links_provenance.schema.json`). +It does not commit the KEGG-derived link table. `data/mappings/metabolite_gene_links_mgigas.*` +is git-ignored. When it has not been built, real-data candidates carry +`link_table_status = not_built` and no linked convergence: the absence is +reported, never silent. CI uses only the simulated demo table. An open-licence +route (Rhea + ChEBI + an eggNOG-mapper run on the RS_2024_06 proteome) remains +future work, because it would allow the table to be committed. + +**B. Scoring: tier and card only.** Linked layers never enter +`n_supporting_layers` or `assay_diversity_score`. + +**C. A fourth tier**, `metabolite_linked_convergence`, after +`multi_omics_convergence` and before `emerging`. + +**D. Demo:** two fabricated demo genes were added to the demo crosswalk and to +GIGAS_HEAT01 `chronic_heat_vs_control`: + +- `ldh` (LOC105342001): significant, linked to lactate through R00703/K00016. +- `sdha` (LOC105342002): not significant, linked to succinate through R02164/K00234. + +Both are simulated. The KEGG reaction and KO IDs in the demo table are real +KEGG identifiers, but the gene IDs they point to are fabricated. + +**E. Cap:** `MAX_GENES_PER_COMPOUND = 25` (`src/mappings/metabolite_links.py`). + +### Confidence rule as implemented + +A metabolite→gene link is never `exact` or `one_to_one_ortholog`. A compound +that reaches more than one crg gene is `one_to_many_ortholog`. A compound that +reaches exactly one gene is `inferred`, because the enzyme assignment is still +KEGG orthology transfer, not oyster biochemistry. A name→compound match by +`curator_judgement` makes the link `inferred`, which is the weaker rank in +`mapping_confidence.yaml`. The same applies when one reported name maps to +more than one compound. For example, `lactate` maps to both L-lactate +(C00186) and D-lactate (C00256), because LC-MS does not resolve +stereoisomers, and molluscs carry both L- and D-lactate dehydrogenases. + +KEGG writes every reaction equation with `<=>`. The table therefore records +the compound's `equation_side` (`left`/`right`), not substrate or product: the +equation does not encode the physiological direction. + +Pumps and transporters have no KEGG REACTION, so they get no link. For +example, SERCA (K05853) has none. A significant ATP change and a significant +`atp2a` change are therefore not linked by this mechanism, and ATP is a +currency metabolite in any case. + +### First local build (2026-10-16, not committed) + +`aree build-metabolite-links` against KEGG REST produced 191 links for 6 of the +7 queried compounds: + +- ATP was skipped as a currency compound. +- Pyruvate (47 genes) exceeds the cap and is a hub, so it is not credited. +- 10,875 of 10,879 KEGG `crg` genes with a KO resolved to current RS_2024_06 + GeneIDs. The 4 unresolved IDs are rRNA/tRNA pseudo-IDs (`s-rRNA-1`, + `tRNA-Arg`, ...), and they are listed in the provenance. +- One biologically meaningful result: **no `crg` gene carries the + L-lactate dehydrogenase KO (K00016)**, so (S)-lactate has no link, and + `lactate` links only through D-lactate (C00256, 2 genes). This is consistent + with bivalve anaerobic metabolism relying on D-LDH and opine dehydrogenases. + It is also why reported `lactate` is mapped to both stereoisomers rather + than assumed to be L. The simulated demo, by contrast, links lactate through + K00016 to a fabricated gene. That is a demo device, not oyster biology. + +## 8. Assumptions + +- KEGG's crg KO assignments are adequate as *inferred* enzyme annotation for + *C. gigas*. They come from orthology, not oyster biochemistry. +- Reaction participation is the only link type in the MVP. Pathway + co-membership is deliberately **not** a link, because it would connect + almost everything. +- The real table is built for RS_2024_06 only. Other assemblies resolve + through NCBI GeneIDs, as identifiers already do. +- The only real metabolomics study (DUPERRET2025_POMS_STARVE) failed QC + (sr320/AREE#70). Real-data convergence is therefore not expected yet, and + tests use the demo. diff --git a/docs/methods.md b/docs/methods.md index ef5b89a..76d41e5 100644 --- a/docs/methods.md +++ b/docs/methods.md @@ -120,7 +120,18 @@ adjustment and each measured under the same phenotype, is assigned to a multi-omics convergence tier, with the specific cross-layer linkage shown explicitly; the link is shared gene identity, not a measured regulatory relationship, and evidence for the same feature under a different phenotype -is displayed but not credited. All remaining candidates — including those with a strong effect +is displayed but not credited. Because a metabolite cannot share a gene +identifier, metabolomics evidence is linked to gene-level evidence +biochemically instead: a curated metabolite name is mapped to KEGG compounds, +then to the KEGG reactions in which the compound participates, the KEGG +orthologs annotated to catalyse those reactions, and the *C. gigas* genes +KEGG assigns to those orthologs (NCBI GeneIDs, checked against the +RS_2024_06 crosswalk). Currency cofactors and compounds linked to more than +25 genes are excluded, and links are never assigned higher than +one-to-many-ortholog confidence. A candidate whose own layer is significant +and whose linked feature is significant in another layer under the same +phenotype is placed in a separate metabolite-linked convergence tier, ranked +below shared-identity convergence and not credited in the numeric score. All remaining candidates — including those with a strong effect in a single study — are assigned to an emerging tier and explicitly labeled as requiring independent replication. This gating structure is the mechanism by which the framework prevents a single study's statistical diff --git a/schemas/metabolite_gene_links_provenance.schema.json b/schemas/metabolite_gene_links_provenance.schema.json new file mode 100644 index 0000000..6d941aa --- /dev/null +++ b/schemas/metabolite_gene_links_provenance.schema.json @@ -0,0 +1,79 @@ +{ + "$schema": "http://json-schema.org/draft-07/schema#", + "$id": "https://github.com/sr320/AREE/schemas/metabolite_gene_links_provenance.schema.json", + "title": "AREE metabolite -> gene link table provenance sidecar", + "description": "Provenance for data/mappings/metabolite_gene_links_.tsv. The real (KEGG-derived) table and its sidecar are built locally by `aree build-metabolite-links` and are not committed; the simulated demo sidecar is. See docs/metabolite_gene_links.md.", + "type": "object", + "required": [ + "artifact", "artifact_sha256", "simulated", "date_generated", "ncbi_taxid", + "source_summary", "sources", "max_genes_per_compound", "known_limitations" + ], + "properties": { + "artifact": {"type": "string", "pattern": "\\.tsv$"}, + "artifact_sha256": {"type": "string", "pattern": "^[0-9a-f]{64}$"}, + "simulated": {"type": "boolean"}, + "builder": {"type": "string"}, + "builder_version": {"type": "string"}, + "date_generated": {"type": "string", "format": "date"}, + "ncbi_taxid": {"type": "integer"}, + "kegg_organism": {"type": "string"}, + "kegg_release": { + "type": "object", + "required": ["summary"], + "properties": {"summary": {"type": "string"}, "ko": {"type": "string"}} + }, + "source_summary": {"type": "string", "minLength": 1}, + "licence_note": {"type": "string"}, + "identifier_crosswalk": {"type": "string"}, + "compound_map": {"$ref": "#/definitions/checksummed_file"}, + "currency_metabolites": { + "allOf": [ + {"$ref": "#/definitions/checksummed_file"}, + {"type": "object", "required": ["skipped"], + "properties": {"skipped": {"type": "array", "items": {"type": "string", "pattern": "^C\\d{5}$"}}}} + ] + }, + "max_genes_per_compound": {"type": "integer", "minimum": 1}, + "sources": { + "type": "array", + "minItems": 1, + "items": { + "type": "object", + "required": ["url"], + "properties": { + "url": {"type": "string"}, + "local_copy": {"type": "string"}, + "sha256": {"type": "string", "pattern": "^[0-9a-f]{64}$"}, + "accessed": {"type": "string"} + } + } + }, + "coverage": { + "type": "object", + "required": ["kegg_genes_with_ko", "kegg_genes_unresolved", "compounds_queried", "hub_compounds", "rows"] + }, + "unresolved_kegg_gene_ids": {"type": "array", "items": {"type": "string"}}, + "known_limitations": {"type": "array", "items": {"type": "string"}} + }, + "if": {"properties": {"simulated": {"const": false}}}, + "then": { + "required": [ + "builder", "builder_version", "kegg_organism", "kegg_release", "licence_note", + "identifier_crosswalk", "compound_map", "currency_metabolites", "coverage", + "unresolved_kegg_gene_ids" + ], + "properties": { + "sources": {"items": {"required": ["url", "local_copy", "sha256"]}} + } + }, + "definitions": { + "checksummed_file": { + "type": "object", + "required": ["path", "sha256"], + "properties": { + "path": {"type": "string"}, + "sha256": {"type": "string", "pattern": "^[0-9a-f]{64}$"} + } + } + } +} diff --git a/src/aree/cli.py b/src/aree/cli.py index ad039fb..c6507bf 100644 --- a/src/aree/cli.py +++ b/src/aree/cli.py @@ -5,6 +5,7 @@ from __future__ import annotations import datetime +import json import subprocess import sys from pathlib import Path @@ -780,5 +781,39 @@ def fail(message): click.echo(click.style(f" identifier_crosswalk: \"{shown}\"", fg="cyan")) + +@main.command("build-metabolite-links") +@click.option("--taxid", default=29159, show_default=True, type=int, + help="NCBI taxid; must have a KEGG organism code in mappings.build_metabolite_links.") +@click.option("--offline", is_flag=True, + help="Rebuild only from cached KEGG responses in data/mappings/_sources/kegg/.") +def build_metabolite_links_cmd(taxid, offline): + """Build the real metabolite -> enzyme gene link table from the KEGG REST API. + + Queries KEGG only for the compounds curated in + data/mappings/metabolite_compound_map.tsv and writes + data/mappings/metabolite_gene_links_.tsv plus a .provenance.json + sidecar. The output is git-ignored and must not be committed (KEGG terms). + See docs/metabolite_gene_links.md. + """ + from mappings.build_metabolite_links import build + from mappings.metabolite_links import clear_caches + + try: + path = build(taxid=taxid, offline=offline) + except (OSError, ValueError) as exc: + click.echo(click.style(str(exc), fg="red")) + sys.exit(1) + clear_caches() + prov = json.loads(path.with_suffix(".provenance.json").read_text()) + cov = prov["coverage"] + click.echo(click.style(f"Wrote {path} ({cov['rows']} links)", fg="green")) + click.echo(f" compounds with links: {cov['compounds_with_links']} of {cov['compounds_queried']} queried" + f" (currency skipped: {', '.join(prov['currency_metabolites']['skipped']) or 'none'})") + click.echo(f" hub compounds (not credited): {', '.join(cov['hub_compounds']) or 'none'}") + click.echo(f" KEGG genes: {cov['kegg_genes_current']} current, {cov['kegg_genes_retired_remapped']} " + f"remapped from retired IDs, {cov['kegg_genes_unresolved']} unresolved") + click.echo(" Not for redistribution: do not commit this file (docs/metabolite_gene_links.md).") + if __name__ == "__main__": main() diff --git a/src/mappings/build_metabolite_links.py b/src/mappings/build_metabolite_links.py new file mode 100644 index 0000000..72694d6 --- /dev/null +++ b/src/mappings/build_metabolite_links.py @@ -0,0 +1,291 @@ +"""Build the real C. gigas metabolite -> gene link table from the KEGG REST API. + + aree build-metabolite-links # fetch (or reuse cached) KEGG responses + aree build-metabolite-links --offline # rebuild from the cached responses only + +Only the compounds named in data/mappings/metabolite_compound_map.tsv are +queried — never a mirror of KEGG — and currency compounds +(data/mappings/currency_metabolites.yaml) are not queried at all. The chain is: + + link/ko/crg every C. gigas gene -> KO assignment (one call) + link/reaction/cpd:C... reactions each curated compound takes part in + link/ko/R...+R... KOs annotated to catalyse those reactions + get/R...+R... reaction equations, for the compound's side + get/cpd:C...+C... compound names + +KEGG gene IDs for `crg` are NCBI GeneIDs. Each one is checked against the +species' AREE crosswalk; a retired GeneID is remapped through the +*_retired_gene_ids.tsv sidecar and counted, and one that resolves to nothing is +dropped and listed in the provenance, never silently. + +The output (metabolite_gene_links_mgigas.tsv + .provenance.json) is NOT +committed: KEGG REST is provided for academic use and KEGG does not clearly +grant redistribution of derived tables. Every raw response is cached under +data/mappings/_sources/kegg/ (also git-ignored) with its sha256 recorded, so a +build can be reproduced offline and audited. See docs/metabolite_gene_links.md. +""" +from __future__ import annotations + +import csv +import hashlib +import json +import re +import urllib.request +from collections import defaultdict +from datetime import date +from pathlib import Path +from typing import Callable + +from common import MAPPINGS_DIR, REPO_ROOT, load_assemblies, sha256_file + +from .metabolite_links import ( + COMPOUND_MAP_PATH, + CURRENCY_PATH, + LINK_COLUMNS, + MAX_GENES_PER_COMPOUND, + load_compound_map, + load_currency_compounds, + provenance_path, +) + +BUILDER_VERSION = "1.0.0" +KEGG_REST = "https://rest.kegg.jp" +KEGG_ORGANISM = {29159: "crg"} +DEFAULT_SLUG = {29159: "mgigas"} +SOURCES_DIR = MAPPINGS_DIR / "_sources" / "kegg" +# KEGG REST accepts at most 10 entries per get/link call. +BATCH = 10 + +Fetcher = Callable[[str], str] + + +def output_path(taxid: int) -> Path: + return MAPPINGS_DIR / f"metabolite_gene_links_{DEFAULT_SLUG[taxid]}.tsv" + + +class CachedKegg: + """Fetch KEGG REST operations, caching each raw response and its sha256.""" + + def __init__(self, cache_dir: Path = SOURCES_DIR, offline: bool = False, fetch: Fetcher | None = None): + self.cache_dir = cache_dir + self.offline = offline + self._fetch = fetch or self._http_get + self.log = [] + + @staticmethod + def _http_get(operation: str) -> str: + with urllib.request.urlopen(f"{KEGG_REST}/{operation}", timeout=60) as response: # noqa: S310 + return response.read().decode("utf-8") + + def get(self, operation: str) -> str: + path = self.cache_dir / (re.sub(r"[^A-Za-z0-9_.+-]+", "_", operation) + ".txt") + if path.exists(): + text = path.read_text() + elif self.offline: + raise FileNotFoundError(f"--offline, and no cached KEGG response for {operation!r} at {path}") + else: + text = self._fetch(operation) + path.parent.mkdir(parents=True, exist_ok=True) + path.write_text(text) + self.log.append({ + "url": f"{KEGG_REST}/{operation}", + "local_copy": str(path.relative_to(REPO_ROOT)) if path.is_relative_to(REPO_ROOT) else str(path), + "sha256": hashlib.sha256(text.encode("utf-8")).hexdigest(), + }) + return text + + +def _batches(items: list, size: int = BATCH): + for i in range(0, len(items), size): + yield items[i:i + size] + + +def _pairs(text: str) -> list: + """Parse a KEGG `link` response into (source, target) pairs without db prefixes.""" + out = [] + for line in text.splitlines(): + if "\t" in line: + a, b = line.split("\t", 1) + out.append((a.split(":", 1)[-1].strip(), b.split(":", 1)[-1].strip())) + return out + + +def parse_entries(text: str, fields: tuple) -> dict: + """Parse a flat-file `get` response into {ENTRY id: {field: first-line value}}.""" + entries, current = {}, None + for line in text.splitlines(): + if line.startswith("ENTRY"): + current = line.split()[1] + entries[current] = {} + elif current and line[:12].strip() in fields and line[:12].strip() not in entries[current]: + entries[current][line[:12].strip()] = line[12:].strip() + return entries + + +def equation_side(equation: str, compound: str) -> str | None: + """'left' / 'right' / 'both' for a compound in a KEGG equation ('A + B <=> C').""" + if "<=>" not in equation: + return None + left, right = equation.split("<=>", 1) + on_left = re.search(rf"\b{compound}\b", left) is not None + on_right = re.search(rf"\b{compound}\b", right) is not None + if on_left and on_right: + return "both" + return "left" if on_left else "right" if on_right else None + + +def load_gene_resolver(taxid: int) -> tuple: + """(current GeneIDs, retired -> current map, crosswalk path) for a taxid.""" + crosswalk = next( + (a["identifier_crosswalk"] for a in load_assemblies().values() + if a.get("ncbi_taxid") == taxid and a.get("is_ncbi_reference") and a.get("identifier_crosswalk")), + None, + ) + if crosswalk is None: + raise ValueError(f"no identifier crosswalk registered for the reference assembly of taxid {taxid}") + path = REPO_ROOT / crosswalk + with open(path) as fh: + current = {r["ncbi_gene_id"] for r in csv.DictReader( + (line for line in fh if not line.startswith("#")), delimiter="\t")} + retired = {} + retired_path = path.with_name(path.name.replace("_gene_id_crosswalk.tsv", "_retired_gene_ids.tsv")) + if retired_path.exists(): + with open(retired_path) as fh: + for r in csv.DictReader((line for line in fh if not line.startswith("#")), delimiter="\t"): + if r.get("current_gene_id") and r["current_gene_id"] != "-": + retired[r["discontinued_gene_id"]] = r["current_gene_id"] + return current, retired, path + + +def build( + taxid: int = 29159, + out_path: Path | None = None, + offline: bool = False, + fetch: Fetcher | None = None, + cache_dir: Path = SOURCES_DIR, + gene_resolver: tuple | None = None, +) -> Path: + """Build the link table and its provenance sidecar; return the table path.""" + if taxid not in KEGG_ORGANISM: + raise ValueError(f"no KEGG organism code recorded for taxid {taxid} (known: {sorted(KEGG_ORGANISM)})") + org = KEGG_ORGANISM[taxid] + out_path = out_path or output_path(taxid) + kegg = CachedKegg(cache_dir=cache_dir, offline=offline, fetch=fetch) + current, retired, crosswalk_path = gene_resolver or load_gene_resolver(taxid) + + release = kegg.get("info/kegg").strip().splitlines() + genes_by_ko = defaultdict(set) + raw_genes = set() + for gene, ko in _pairs(kegg.get(f"link/ko/{org}")): + raw_genes.add(gene) + genes_by_ko[ko].add(gene) + + resolved, n_retired_remapped, unresolved = {}, 0, [] + for gene in sorted(raw_genes): + if gene in current: + resolved[gene] = gene + elif gene in retired and retired[gene] in current: + resolved[gene] = retired[gene] + n_retired_remapped += 1 + else: + unresolved.append(gene) + + currency = load_currency_compounds(CURRENCY_PATH) + compound_map = load_compound_map(COMPOUND_MAP_PATH) + curated = sorted({c["kegg_compound"] for cs in compound_map.values() for c in cs}) + queried = [c for c in curated if c not in currency] + + reactions_by_compound = { + cpd: sorted({r for _, r in _pairs(kegg.get(f"link/reaction/cpd:{cpd}"))}) for cpd in queried + } + all_reactions = sorted({r for rs in reactions_by_compound.values() for r in rs}) + kos_by_reaction = defaultdict(set) + equations = {} + for batch in _batches(all_reactions): + for reaction, ko in _pairs(kegg.get("link/ko/" + "+".join(f"rn:{r}" for r in batch))): + kos_by_reaction[reaction].add(ko) + for entry, fields in parse_entries(kegg.get("get/" + "+".join(f"rn:{r}" for r in batch)), + ("EQUATION",)).items(): + equations[entry] = fields.get("EQUATION", "") + names = {} + for batch in _batches(queried): + for entry, fields in parse_entries(kegg.get("get/" + "+".join(f"cpd:{c}" for c in batch)), + ("NAME",)).items(): + names[entry] = fields.get("NAME", "").rstrip(";") + + rows = [] + for cpd in queried: + for reaction in reactions_by_compound[cpd]: + side = equation_side(equations.get(reaction, ""), cpd) + for ko in sorted(kos_by_reaction.get(reaction, ())): + for gene in sorted(genes_by_ko.get(ko, ())): + if gene not in resolved: + continue + rows.append({ + "kegg_compound": cpd, "compound_name": names.get(cpd, ""), + "kegg_reaction": reaction, "equation_side": side or "unknown", "ko": ko, + "ncbi_gene_id": resolved[gene], "link_type": "enzyme", + }) + rows = [dict(t) for t in sorted({tuple(r.items()) for r in rows}, + key=lambda t: tuple(v for _, v in t))] + + out_path.parent.mkdir(parents=True, exist_ok=True) + with open(out_path, "w") as fh: + fh.write(f"# AREE metabolite -> gene link table, taxid {taxid} (KEGG organism {org}).\n") + fh.write(f"# Built by src/mappings/build_metabolite_links.py v{BUILDER_VERSION} on " + f"{date.today().isoformat()} from the KEGG REST API. NOT for redistribution; see\n") + fh.write("# the .provenance.json sidecar and docs/metabolite_gene_links.md.\n") + fh.write("\t".join(LINK_COLUMNS) + "\n") + for r in rows: + fh.write("\t".join(str(r[c]) for c in LINK_COLUMNS) + "\n") + + genes_per_compound = defaultdict(set) + for r in rows: + genes_per_compound[r["kegg_compound"]].add(r["ncbi_gene_id"]) + release_line = next((line.strip() for line in release if line.strip().startswith("kegg")), release[0]) + ko_line = next((" ".join(line.split()) for line in release if line.strip().startswith("ko ")), "") + provenance = { + "artifact": out_path.name, + "artifact_sha256": sha256_file(out_path), + "simulated": False, + "builder": "src/mappings/build_metabolite_links.py", + "builder_version": BUILDER_VERSION, + "date_generated": date.today().isoformat(), + "ncbi_taxid": taxid, + "kegg_organism": org, + "kegg_release": {"summary": release_line, "ko": ko_line}, + "source_summary": f"KEGG REST {org}, accessed {date.today().isoformat()} ({ko_line or release_line})", + "licence_note": ("KEGG REST is provided for academic use by academic users. This derived table " + "is not committed or redistributed by AREE; rebuild it locally."), + "identifier_crosswalk": str(crosswalk_path.relative_to(REPO_ROOT)) + if crosswalk_path.is_relative_to(REPO_ROOT) else str(crosswalk_path), + "compound_map": {"path": str(COMPOUND_MAP_PATH.relative_to(REPO_ROOT)), + "sha256": sha256_file(COMPOUND_MAP_PATH)}, + "currency_metabolites": {"path": str(CURRENCY_PATH.relative_to(REPO_ROOT)), + "sha256": sha256_file(CURRENCY_PATH), + "skipped": sorted(c for c in curated if c in currency)}, + "max_genes_per_compound": MAX_GENES_PER_COMPOUND, + "sources": kegg.log, + "coverage": { + "kegg_genes_with_ko": len(raw_genes), + "kegg_genes_current": len(raw_genes) - n_retired_remapped - len(unresolved), + "kegg_genes_retired_remapped": n_retired_remapped, + "kegg_genes_unresolved": len(unresolved), + "compounds_curated": len(curated), + "compounds_queried": len(queried), + "compounds_with_links": len(genes_per_compound), + "hub_compounds": sorted(c for c, g in genes_per_compound.items() if len(g) > MAX_GENES_PER_COMPOUND), + "genes_per_compound": {c: len(g) for c, g in sorted(genes_per_compound.items())}, + "rows": len(rows), + }, + "unresolved_kegg_gene_ids": unresolved, + "known_limitations": [ + "KO assignments are KEGG's orthology transfer to C. gigas, not oyster biochemistry; every " + "link is 'inferred' or 'one_to_many_ortholog', never 'exact'.", + "KEGG equations are written reversibly; equation_side is not substrate/product.", + "Pumps and transporters without a KEGG REACTION (e.g. SERCA, K05853) get no link.", + "Only compounds in metabolite_compound_map.tsv are queried; add names there first.", + ], + } + provenance_path(out_path).write_text(json.dumps(provenance, indent=2) + "\n") + return out_path diff --git a/src/mappings/metabolite_links.py b/src/mappings/metabolite_links.py new file mode 100644 index 0000000..6ea152b --- /dev/null +++ b/src/mappings/metabolite_links.py @@ -0,0 +1,300 @@ +"""Metabolite -> enzyme gene links, the weaker second kind of multi-omics convergence. + +Gene-level layers (transcriptomics, proteomics, DNA methylation) converge by +sharing a standardized gene identifier. A metabolite is not a gene product, so +it can only be *linked* to genes biochemically: + + reported metabolite name + -> KEGG compound (curated: data/mappings/metabolite_compound_map.tsv) + -> KEGG reaction (the compound appears on either side of the equation) + -> KEGG orthology (the KO annotated to catalyse that reaction) + -> crg gene (KEGG's KO assignment for C. gigas = an NCBI GeneID) + +Every step is annotation, not measurement, so a link is never `exact` and never +`one_to_one_ortholog` (see docs/metabolite_gene_links.md section 7). Currency +metabolites (ATP, NAD+, H2O, ...) and hub compounds reaching more than +MAX_GENES_PER_COMPOUND genes are kept visible but never credited: without that +filter a change in ATP would "link" to most of the metabolic genome. + +Link tables are chosen per evidence partition, like the identifier crosswalk: +simulated evidence resolves against the simulated demo table and real evidence +against the table registered for its species' assembly in +data/reference/genome_assemblies.yaml (`metabolite_gene_links`). The two are +never merged. The real KEGG-derived table is built locally by +`aree build-metabolite-links` and is not committed (KEGG redistribution terms). +""" +from __future__ import annotations + +import csv +import re +from collections import defaultdict +from dataclasses import dataclass, field +from functools import lru_cache +from pathlib import Path + +from common import ( + MAPPINGS_DIR, + REPO_ROOT, + load_assemblies, + load_json, + load_yaml, + sha256_file, +) + +COMPOUND_MAP_PATH = MAPPINGS_DIR / "metabolite_compound_map.tsv" +CURRENCY_PATH = MAPPINGS_DIR / "currency_metabolites.yaml" +DEMO_LINKS_PATH = MAPPINGS_DIR / "metabolite_gene_links_demo.tsv" + +# A compound annotated to more crg genes than this is a hub: crediting it would +# turn "this metabolite changed" into "much of the metabolic genome is linked". +MAX_GENES_PER_COMPOUND = 25 + +LINK_COLUMNS = [ + "kegg_compound", "compound_name", "kegg_reaction", "equation_side", "ko", "ncbi_gene_id", "link_type", +] + +# Link-table availability for a partition, as reported on every candidate. +STATUS_AVAILABLE = "available" +STATUS_NOT_BUILT = "not_built" # registered for the assembly, not built on this machine +STATUS_NOT_REGISTERED = "not_registered" # no table exists for this species / data origin + + +def normalize_metabolite_name(name) -> str: + """Lowercase, trim, and fold spaces/hyphens to underscores: 'Glycine betaine' -> 'glycine_betaine'.""" + return re.sub(r"[\s\-]+", "_", str(name).strip().casefold()) + + +def _read_tsv(path: Path) -> list: + with open(path) as fh: + return list(csv.DictReader((line for line in fh if not line.startswith("#")), delimiter="\t")) + + +@lru_cache(maxsize=4) +def load_compound_map(path: Path = COMPOUND_MAP_PATH) -> dict: + """normalized metabolite name -> list of {kegg_compound, chebi_id, match_basis, note}.""" + out = defaultdict(list) + for row in _read_tsv(path): + out[normalize_metabolite_name(row["metabolite_name"])].append({ + "kegg_compound": row["kegg_compound"], + "chebi_id": row.get("chebi_id") or None, + "match_basis": row["match_basis"], + "note": row.get("note") or "", + }) + return dict(out) + + +@lru_cache(maxsize=4) +def load_currency_compounds(path: Path = CURRENCY_PATH) -> dict: + """KEGG compound -> reason it is never used as a link.""" + doc = load_yaml(path) + return {c["kegg_compound"]: c["reason"] for c in doc["compounds"]} + + +@dataclass(frozen=True) +class MetaboliteGeneLink: + """One metabolite -> compound -> reaction -> KO -> gene path.""" + + metabolite: str + kegg_compound: str + compound_match_basis: str + kegg_reaction: str + equation_side: str + ko: str + ncbi_gene_id: str + link_type: str + n_genes_for_compound: int + link_confidence: str + credited: bool + exclusion_reason: str | None = None + + def describe(self) -> str: + return (f"{self.metabolite}→{self.ncbi_gene_id} via {self.kegg_compound}/" + f"{self.kegg_reaction}/{self.ko}") + + +@dataclass(frozen=True) +class CompoundExclusion: + """A compound a metabolite resolves to that contributes no credited link.""" + + metabolite: str + kegg_compound: str + reason: str # currency_metabolite | hub_compound | no_reaction_in_table + + +def link_confidence(match_basis: str, n_compounds_for_name: int, n_genes_for_compound: int) -> str: + """Mapping confidence of a metabolite -> gene link (never 'exact' or one-to-one). + + A curator's judgement call on the name, or a name that resolves to more + than one compound (e.g. 'lactate' -> L- and D-lactate), is the weakest + rung, `inferred`. Otherwise a compound reaching several genes is + `one_to_many_ortholog`, and one reaching a single gene is still `inferred`: + the enzyme assignment is KEGG orthology transfer, not oyster biochemistry. + """ + if match_basis == "curator_judgement" or n_compounds_for_name > 1: + return "inferred" + if n_genes_for_compound > 1: + return "one_to_many_ortholog" + return "inferred" + + +@dataclass +class MetaboliteGeneLinks: + """A loaded link table, indexed in both directions.""" + + path: Path | None + status: str + sha256: str | None = None + source: str | None = None + links: list = field(default_factory=list) + exclusions: list = field(default_factory=list) + + def __post_init__(self): + self._by_metabolite = defaultdict(list) + self._by_gene = defaultdict(list) + self._exclusions_by_metabolite = defaultdict(list) + for link in self.links: + self._by_metabolite[link.metabolite].append(link) + self._by_gene[link.ncbi_gene_id].append(link) + for exclusion in self.exclusions: + self._exclusions_by_metabolite[exclusion.metabolite].append(exclusion) + + @property + def available(self) -> bool: + return self.status == STATUS_AVAILABLE + + def for_metabolite(self, name) -> list: + return self._by_metabolite.get(normalize_metabolite_name(name), []) + + def for_gene(self, gene_id) -> list: + return self._by_gene.get(str(gene_id), []) + + def excluded_for_metabolite(self, name) -> list: + return self._exclusions_by_metabolite.get(normalize_metabolite_name(name), []) + + def links_for_feature(self, feature_id, is_metabolite: bool) -> list: + """Links touching this standardized feature, from whichever side it is on.""" + return self.for_metabolite(feature_id) if is_metabolite else self.for_gene(feature_id) + + def label(self) -> str: + if not self.available: + return self.status + rel = self.path.relative_to(REPO_ROOT) if self.path.is_relative_to(REPO_ROOT) else self.path + return f"{rel} (sha256 {self.sha256[:12]}…, {self.source})" + + +def build_links( + link_rows: list, + compound_map: dict, + currency: dict, + max_genes_per_compound: int = MAX_GENES_PER_COMPOUND, +) -> tuple: + """Expand link-table rows into per-metabolite links plus visible exclusions. + + Returns (links, exclusions). Pure: tests build tables from plain rows. + """ + rows_by_compound = defaultdict(list) + for row in link_rows: + rows_by_compound[row["kegg_compound"]].append(row) + genes_by_compound = { + cpd: {str(r["ncbi_gene_id"]) for r in rows} for cpd, rows in rows_by_compound.items() + } + + links, exclusions = [], [] + for metabolite, compounds in compound_map.items(): + n_compounds = len(compounds) + for compound in compounds: + cpd = compound["kegg_compound"] + if cpd in currency: + exclusions.append(CompoundExclusion(metabolite, cpd, "currency_metabolite")) + continue + rows = rows_by_compound.get(cpd) + if not rows: + exclusions.append(CompoundExclusion(metabolite, cpd, "no_reaction_in_table")) + continue + n_genes = len(genes_by_compound[cpd]) + hub = n_genes > max_genes_per_compound + if hub: + exclusions.append(CompoundExclusion(metabolite, cpd, "hub_compound")) + confidence = link_confidence(compound["match_basis"], n_compounds, n_genes) + seen = set() + for row in rows: + key = (row["kegg_reaction"], row["ko"], str(row["ncbi_gene_id"])) + if key in seen: + continue + seen.add(key) + links.append(MetaboliteGeneLink( + metabolite=metabolite, + kegg_compound=cpd, + compound_match_basis=compound["match_basis"], + kegg_reaction=row["kegg_reaction"], + equation_side=row["equation_side"], + ko=row["ko"], + ncbi_gene_id=str(row["ncbi_gene_id"]), + link_type=row.get("link_type") or "enzyme", + n_genes_for_compound=n_genes, + link_confidence=confidence, + credited=not hub, + exclusion_reason="hub_compound" if hub else None, + )) + return links, exclusions + + +def provenance_path(table_path: Path) -> Path: + return table_path.with_suffix(".provenance.json") + + +def load_link_table(path: Path) -> MetaboliteGeneLinks: + """Load and index a built link table, or report why there is none.""" + if not path.exists(): + return MetaboliteGeneLinks(path=path, status=STATUS_NOT_BUILT) + provenance = load_json(provenance_path(path)) if provenance_path(path).exists() else {} + links, exclusions = build_links(_read_tsv(path), load_compound_map(), load_currency_compounds()) + return MetaboliteGeneLinks( + path=path, + status=STATUS_AVAILABLE, + sha256=sha256_file(path), + source=provenance.get("source_summary", "no provenance sidecar"), + links=links, + exclusions=exclusions, + ) + + +def _demo_taxid() -> str | None: + sidecar = provenance_path(DEMO_LINKS_PATH) + return str(load_json(sidecar)["ncbi_taxid"]) if sidecar.exists() else None + + +def registered_link_path(taxid) -> Path | None: + """The real link table registered for any assembly of this taxid, or None.""" + for record in load_assemblies().values(): + if str(record.get("ncbi_taxid")) != str(taxid): + continue + value = record.get("metabolite_gene_links") + if value: + path = Path(value) + return path if path.is_absolute() else REPO_ROOT / path + return None + + +@lru_cache(maxsize=16) +def links_for_partition(simulated: str, taxid: str) -> MetaboliteGeneLinks: + """The link table for one evidence partition (`_simulated`, `_species_taxid`). + + Simulated evidence only ever sees the simulated demo table, and only for the + taxid it was fabricated for. Real evidence only ever sees a real table. + """ + if simulated == "true": + if taxid == _demo_taxid(): + return load_link_table(DEMO_LINKS_PATH) + return MetaboliteGeneLinks(path=None, status=STATUS_NOT_REGISTERED) + path = registered_link_path(taxid) + if path is None: + return MetaboliteGeneLinks(path=None, status=STATUS_NOT_REGISTERED) + return load_link_table(path) + + +def clear_caches() -> None: + """Forget loaded tables (after a rebuild, or in tests that swap files).""" + load_compound_map.cache_clear() + load_currency_compounds.cache_clear() + links_for_partition.cache_clear() diff --git a/src/prioritize/rank.py b/src/prioritize/rank.py index f0f2f38..f787b36 100644 --- a/src/prioritize/rank.py +++ b/src/prioritize/rank.py @@ -17,6 +17,7 @@ import pandas as pd from common import load_vocab +from mappings.metabolite_links import links_for_partition from .scoring import ( candidate_score, @@ -41,8 +42,22 @@ MULTI_OMICS_MIN_LAYERS = 2 LAYER_SUPPORT_MAX_ADJUSTED_P = 0.05 +# A metabolite cannot share a gene identifier, so metabolomics is linked to the +# gene-level layers biochemically instead (metabolite -> KEGG reaction -> KO -> +# C. gigas gene; src/mappings/metabolite_links.py). That is annotation, not a +# second measurement of the same molecule, so it is a separate, weaker tier: +# never merged into `multi_omics_convergence`, never counted in +# n_supporting_layers, and never part of the numeric score +# (docs/metabolite_gene_links.md). +METABOLITE_LAYER = "metabolomics" + # Order in which tiers are presented: strongest evidence first. -TIER_ORDER = ["high_priority_cross_study", "multi_omics_convergence", "emerging"] +TIER_ORDER = [ + "high_priority_cross_study", + "multi_omics_convergence", + "metabolite_linked_convergence", + "emerging", +] # What KIND of evidence supports a candidate, independent of how strong it is. # The tier answers "how well replicated?"; this answers "replicated evidence of @@ -213,6 +228,60 @@ def supporting_layers(self, meta_row: dict) -> set: layers.add(layer_of(meta_row["feature_type"])) return layers + def links_for(self, meta_row: dict) -> list: + """Metabolite <-> gene links touching this candidate's feature, from the + link table of its own partition (simulated demo or real).""" + key = candidate_key(meta_row) + table = links_for_partition(key[3], key[4]) + if not table.available: + return [] + return table.links_for_feature(key[0], layer_of(meta_row["feature_type"]) == METABOLITE_LAYER) + + @staticmethod + def linked_feature_id(link, from_metabolite: bool) -> str: + return link.ncbi_gene_id if from_metabolite else link.metabolite + + def linked_records(self, meta_row: dict) -> pd.DataFrame: + """Evidence for every feature linked to this candidate, under any + phenotype, from layers other than the candidate's own, with the link + that connects each record in `_link`. Callers use `phenotype`, + `_significant` and `_link.credited` to tell support from context.""" + key = candidate_key(meta_row) + own_layer = layer_of(meta_row["feature_type"]) + from_metabolite = own_layer == METABOLITE_LAYER + frames = [] + for link in self.links_for(meta_row): + linked_id = self.linked_feature_id(link, from_metabolite) + records = self._slice(self._by_feature, (linked_id, key[3], key[4])) + records = records[records["_layer"] != own_layer] + if len(records): + frames.append(records.assign(_link=[link] * len(records))) + if not frames: + return self.evidence.iloc[0:0].assign(_link=None) + return pd.concat(frames) + + def linked_support(self, meta_row: dict, supporting: set) -> tuple: + """(linked_layers, credited links behind them) for this candidate. + + A linked layer is one with a significant record for a CREDITED linked + feature under this candidate's own phenotype and partition, that is not + already a shared-identifier supporting layer — shared identity is the + stronger link and is never double-counted as a biochemical one. + """ + key = candidate_key(meta_row) + from_metabolite = layer_of(meta_row["feature_type"]) == METABOLITE_LAYER + layers, links = set(), [] + for link in self.links_for(meta_row): + if not link.credited: + continue + linked_id = self.linked_feature_id(link, from_metabolite) + found = self._supported_layers.get((linked_id, key[1], key[3], key[4]), set()) - supporting + found.discard(layer_of(meta_row["feature_type"])) + if found: + layers |= found + links.append(link) + return layers, links + def build_candidates(meta_df: pd.DataFrame, evidence_df: pd.DataFrame) -> pd.DataFrame: if len(meta_df) == 0: @@ -228,6 +297,9 @@ def build_candidates(meta_df: pd.DataFrame, evidence_df: pd.DataFrame) -> pd.Dat own_layer = layer_of(meta_dict["feature_type"]) supporting = index.supporting_layers(meta_dict) n_supporting_layers = len(supporting) + linked_layers, linked_links = index.linked_support(meta_dict, supporting) + link_table = links_for_partition(_partition_value(meta_dict["simulated"]), + _partition_value(meta_dict["species_taxid"])) classifications = (subset["resilience_classification"] if "resilience_classification" in subset.columns else []) @@ -250,14 +322,26 @@ def build_candidates(meta_df: pd.DataFrame, evidence_df: pd.DataFrame) -> pd.Dat # layers converging on a gene this candidate itself shows no signal for # is not evidence for this candidate. is_multi_omics = own_layer in supporting and n_supporting_layers >= MULTI_OMICS_MIN_LAYERS + # Same own-layer rule for the biochemical link: the candidate itself must + # carry signal, and some other layer must carry it for a linked feature. + is_metabolite_linked = own_layer in supporting and bool(linked_layers) if is_high_priority: tier = "high_priority_cross_study" elif is_multi_omics: tier = "multi_omics_convergence" + elif is_metabolite_linked: + tier = "metabolite_linked_convergence" else: tier = "emerging" + if is_multi_omics: + convergence_basis = "shared_identifier" + elif is_metabolite_linked: + convergence_basis = "metabolite_gene_link" + else: + convergence_basis = "none" + rows.append({ **meta_dict, "score": score, @@ -271,6 +355,14 @@ def build_candidates(meta_df: pd.DataFrame, evidence_df: pd.DataFrame) -> pd.Dat "molecular_layer": own_layer, "n_supporting_layers": n_supporting_layers, "supporting_layers": "|".join(sorted(supporting)), + "is_metabolite_linked_convergence": is_metabolite_linked, + "convergence_basis": convergence_basis, + "n_linked_layers": len(linked_layers), + "linked_layers": "|".join(sorted(linked_layers)), + "linked_features": "|".join(sorted({link.describe() for link in linked_links})), + "link_confidences": "|".join(sorted({link.link_confidence for link in linked_links})), + "link_table_status": link_table.status, + "link_table_sha256": link_table.sha256 or "", "mapping_confidences": "|".join(mapping_confidences), "quality_flags_union": "|".join(quality_flags), **{f"component_{k}": v for k, v in components.items()}, diff --git a/src/reporting/dashboard.py b/src/reporting/dashboard.py index 7cf1aa2..5562876 100644 --- a/src/reporting/dashboard.py +++ b/src/reporting/dashboard.py @@ -518,6 +518,7 @@ def _tier_count(simulated: bool, tier: str) -> int: "n_evidence_records_simulated": int(evidence["simulated"].sum()) if len(evidence) else 0, "n_candidates_real_high_priority": _tier_count(False, "high_priority_cross_study"), "n_candidates_real_multi_omics": _tier_count(False, "multi_omics_convergence"), + "n_candidates_real_metabolite_linked": _tier_count(False, "metabolite_linked_convergence"), "n_candidates_real_emerging": _tier_count(False, "emerging"), "n_real_stressor_classes": len(real_stressors), "real_stressor_classes": real_stressors, diff --git a/src/reporting/evidence_cards.py b/src/reporting/evidence_cards.py index e1c2277..62b8121 100644 --- a/src/reporting/evidence_cards.py +++ b/src/reporting/evidence_cards.py @@ -26,8 +26,10 @@ from common import REPORTS_DIR, load_vocab, species_term_for_taxid from harmonize.core import load_evidence_table +from mappings.metabolite_links import MAX_GENES_PER_COMPOUND, links_for_partition from meta_analysis.run import run_meta_analysis from prioritize.rank import ( + METABOLITE_LAYER, EvidenceIndex, _partition_value, add_partition_columns, @@ -47,6 +49,11 @@ "Investigate the regulatory link across molecular layers (e.g. paired methylation/expression " "measurement in the same individuals) before committing to validation." ), + "metabolite_linked_convergence": ( + "Measure the metabolite and the linked enzyme's transcript or activity in the same animals " + "and the same stressor before treating the link as biological; the current link is KEGG " + "reaction annotation transferred by orthology, not a measured relationship." + ), "emerging": ( "Requires replication in an independent study before further validation investment." ), @@ -195,6 +202,118 @@ def _other_omics_context(index: EvidenceIndex, candidate: dict) -> list: return rows +EXCLUSION_REASON_TEXT = { + "currency_metabolite": "currency metabolite (data/mappings/currency_metabolites.yaml); never used as a link", + "hub_compound": f"hub compound reaching more than {MAX_GENES_PER_COMPOUND} genes; not credited", + "no_reaction_in_table": "no reaction for this compound in the link table", +} + + +def _linked_omics_context(index: EvidenceIndex, candidate: dict) -> dict: + """Evidence reached through a metabolite <-> gene link, for the card. + + Every linked record is shown, under any phenotype. A row counts toward + `metabolite_linked_convergence` only when its link is credited, it is for + this candidate's own phenotype, its study reported adjusted p <= 0.05, and + its layer is not already a shared-identifier supporting layer.""" + table = links_for_partition(_partition_value(candidate["simulated"]), + _partition_value(candidate["species_taxid"])) + is_metabolite = candidate.get("molecular_layer") == METABOLITE_LAYER + supporting = {x for x in str(candidate.get("supporting_layers") or "").split("|") if x} + rows = [] + for record in index.linked_records(candidate).to_dict("records"): + link = record["_link"] + same_phenotype = record["phenotype"] == candidate["phenotype"] + significant = bool(record["_significant"]) + rows.append({ + "study_id": record["study_id"], + "layer": record["_layer"], + "linked_feature": index.linked_feature_id(link, is_metabolite), + "kegg_compound": link.kegg_compound, + "kegg_reaction": link.kegg_reaction, + "ko": link.ko, + "equation_side": link.equation_side, + "link_confidence": link.link_confidence, + "n_genes_for_compound": link.n_genes_for_compound, + "credited": link.credited, + "phenotype": record["phenotype"], + "stressor": record["stressor"], + "molecular_direction": record["molecular_direction"], + "effect_size": record["effect_size"], + "adjusted_p_value": record["adjusted_p_value"], + "same_phenotype": same_phenotype, + "significant": significant, + "counts": (link.credited and same_phenotype and significant + and record["_layer"] not in supporting), + }) + exclusions = ( + [{"kegg_compound": e.kegg_compound, "reason": e.reason} + for e in table.excluded_for_metabolite(candidate["feature_id_standardized"])] + if is_metabolite and table.available else [] + ) + return {"table": table, "rows": rows, "exclusions": exclusions, "is_metabolite": is_metabolite} + + +def _linked_section(candidate: dict, linked: dict) -> list: + # A gene with no linked metabolite evidence gets no section: the table's + # status is reported on every metabolite card and in candidates.tsv instead. + if not linked["is_metabolite"] and not linked["rows"]: + return [] + table = linked["table"] + lines = ["", "### Metabolite–gene links (biochemical; weaker than shared gene identity)"] + if not table.available: + lines.append( + f"No metabolite–gene link table for this partition (`{table.status}`). " + + ("Build it with `aree build-metabolite-links`; until then metabolomics evidence cannot " + "be linked to gene-level evidence for this species." if table.status == "not_built" + else "Links are only defined for species with a registered table.") + ) + return lines + lines.append( + "These links come from KEGG reaction annotation, not from measurement: the linked gene's " + "product is *annotated* (by KEGG orthology transfer) to catalyse a reaction in which the " + "metabolite takes part. A link does not mean the gene regulates, or controls the level of, " + "the metabolite. KEGG writes reactions as reversible equations, so the side the compound is " + "on does not say whether it is consumed or produced. Linked layers do not add to this " + "candidate's score." + ) + lines.append("") + lines.append(f"Link table: {table.label()}.") + stressors = {r["stressor"] for r in linked["rows"] if r["counts"]} + own_stressors = {x for x in str(candidate.get("pooled_stressors") or "").split("|") if x} + if stressors and own_stressors and not stressors <= own_stressors: + lines.append("") + lines.append( + f"⚠️ The linked support comes from a different stressor (`{'|'.join(sorted(stressors))}`) " + f"than this candidate's own evidence (`{'|'.join(sorted(own_stressors))}`). Same " + "phenotype, different exposure: treat as a hypothesis connecting two experiments." + ) + if linked["rows"]: + lines.append("") + lines.append("Only rows marked **yes** in the last column count toward " + "`metabolite_linked_convergence`:") + lines.append("| study | layer | linked feature | compound | reaction | KO | side | link confidence | " + "genes for compound | credited | phenotype | stressor | direction | effect | adj. p | counts |") + lines.append("|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|") + for r in linked["rows"]: + lines.append( + f"| {r['study_id']} | {r['layer']} | {r['linked_feature']} | {r['kegg_compound']} | " + f"{r['kegg_reaction']} | {r['ko']} | {r['equation_side']} | {r['link_confidence']} | " + f"{r['n_genes_for_compound']} | {'yes' if r['credited'] else 'no (hub)'} | {r['phenotype']} | " + f"{r['stressor']} | {r['molecular_direction']} | {_fmt(r['effect_size'], '.3f')} | " + f"{_fmt(r['adjusted_p_value'])} | {'**yes**' if r['counts'] else 'no'} |" + ) + else: + lines.append("") + lines.append("No evidence in other layers for any feature linked to this one.") + if linked["exclusions"]: + lines.append("") + lines.append("Compounds this metabolite resolves to that contribute no credited link:") + for e in linked["exclusions"]: + lines.append(f"- `{e['kegg_compound']}`: {EXCLUSION_REASON_TEXT.get(e['reason'], e['reason'])}") + return lines + + def has_significant_signal(candidate: dict, forest_points: list, max_adjusted_p: float) -> bool: """Does anything about this candidate clear the significance threshold? @@ -212,7 +331,8 @@ def has_significant_signal(candidate: dict, forest_points: list, max_adjusted_p: ) -def _card_markdown(candidate: dict, forest_points: list, other_omics: list, phenotype_label: str) -> str: +def _card_markdown(candidate: dict, forest_points: list, other_omics: list, phenotype_label: str, + linked: dict | None = None) -> str: lines = [] lines.append(f"# Evidence card: {candidate['feature_id_standardized']} — {phenotype_label}") lines.append("") @@ -321,6 +441,17 @@ def _card_markdown(candidate: dict, forest_points: list, other_omics: list, phen else: lines.append("") lines.append("No evidence from other molecular layers for this standardized feature yet.") + if candidate.get("is_metabolite_linked_convergence"): + lines.append("") + lines.append( + "**How the layers were linked:** not by a shared identifier. Layer(s) " + f"{', '.join(f'`{x}`' for x in str(candidate.get('linked_layers')).split('|'))} were linked " + f"through metabolite–enzyme annotation ({candidate.get('linked_features')}; link confidence: " + f"{candidate.get('link_confidences')}). This is a weaker kind of convergence than shared gene " + "identity and is ranked below it." + ) + if linked is not None: + lines.extend(_linked_section(candidate, linked)) lines.append("") lines.append("## Limitations") lines.append(f"- Quality flags present across contributing evidence: {candidate['quality_flags_union'] or 'none recorded'}") @@ -419,7 +550,8 @@ def build_evidence_cards_report( other_omics = _other_omics_context(index, candidate) phenotype_label = phenotype_vocab.get(candidate["phenotype"], {}).get("label", candidate["phenotype"]) card_path = EVIDENCE_CARDS_DIR / f"{_card_slug(candidate)}.md" - card_path.write_text(_card_markdown(candidate, forest_points, other_omics, phenotype_label)) + linked = _linked_omics_context(index, candidate) + card_path.write_text(_card_markdown(candidate, forest_points, other_omics, phenotype_label, linked)) card_files.append(str(card_path)) written.append({ "feature_id_standardized": candidate["feature_id_standardized"], diff --git a/tests/test_metabolite_links.py b/tests/test_metabolite_links.py new file mode 100644 index 0000000..63ff7d9 --- /dev/null +++ b/tests/test_metabolite_links.py @@ -0,0 +1,388 @@ +"""Metabolite -> gene links: the weaker, separately labelled kind of convergence. + +See docs/metabolite_gene_links.md. +""" +import json +import re +from pathlib import Path + +import jsonschema +import pandas as pd +import pytest +from conftest import harmonize_all_demo_studies, oyster + +import mappings.metabolite_links as ml +import prioritize.rank as rank +import reporting.evidence_cards as evidence_cards +from common import REPO_ROOT, SCHEMAS_DIR, load_json, load_vocab, sha256_file +from harmonize.core import load_evidence_table +from mappings.build_metabolite_links import build, equation_side, parse_entries +from prioritize.rank import TIER_ORDER, build_candidates + +LDH, SDHA = "LOC105342001", "LOC105342002" + + +@pytest.fixture(autouse=True) +def _fresh_link_caches(): + ml.clear_caches() + yield + ml.clear_caches() + + +# --------------------------------------------------------------------------- # +# Curated inputs +# --------------------------------------------------------------------------- # + + +def test_compound_map_is_well_formed(): + compound_map = ml.load_compound_map() + assert compound_map, "compound map is empty" + for name, compounds in compound_map.items(): + assert name == ml.normalize_metabolite_name(name) + for c in compounds: + assert re.fullmatch(r"C\d{5}", c["kegg_compound"]), c + assert c["match_basis"] in {"exact_name", "synonym", "curator_judgement"} + assert c["chebi_id"] is None or c["chebi_id"].startswith("CHEBI:") + + +def test_every_demo_metabolite_name_with_a_standardized_id_is_curated(): + """A level 1-2 demo metabolite missing from the map would silently never link.""" + table = pd.read_csv( + REPO_ROOT / "data/demo/metabolomics/GIGAS_GROW06_hypoxia_vs_control_metabolite_features_demo.tsv", + sep="\t", + ) + names = table.loc[table["annotation_confidence_level"] <= 2, "putative_metabolite_name"] + missing = {n for n in names if ml.normalize_metabolite_name(n) not in ml.load_compound_map()} + assert not missing + + +def test_atp_is_a_currency_metabolite(): + assert "C00002" in ml.load_currency_compounds() + + +def test_demo_link_table_provenance_matches_the_file_and_is_simulated(): + prov = load_json(ml.provenance_path(ml.DEMO_LINKS_PATH)) + assert prov["simulated"] is True + assert prov["artifact_sha256"] == sha256_file(ml.DEMO_LINKS_PATH) + jsonschema.validate(prov, load_json(SCHEMAS_DIR / "metabolite_gene_links_provenance.schema.json")) + + +def test_real_kegg_link_table_is_not_tracked_by_git(): + """KEGG redistribution terms: the real table is built locally, never committed.""" + ignored = (REPO_ROOT / ".gitignore").read_text() + assert "data/mappings/metabolite_gene_links_mgigas.*" in ignored + + +# --------------------------------------------------------------------------- # +# Mapping confidence +# --------------------------------------------------------------------------- # + + +@pytest.mark.parametrize("basis,n_compounds,n_genes,expected", [ + ("exact_name", 1, 1, "inferred"), + ("exact_name", 1, 4, "one_to_many_ortholog"), + ("synonym", 1, 4, "one_to_many_ortholog"), + ("curator_judgement", 1, 4, "inferred"), + ("exact_name", 2, 4, "inferred"), +]) +def test_link_confidence_is_never_exact(basis, n_compounds, n_genes, expected): + confidence = ml.link_confidence(basis, n_compounds, n_genes) + assert confidence == expected + assert confidence not in {"exact", "one_to_one_ortholog"} + assert confidence in load_vocab("mapping_confidence") + + +def _row(cpd, gene, reaction="R1", ko="K1", side="left"): + return {"kegg_compound": cpd, "compound_name": cpd, "kegg_reaction": reaction, + "equation_side": side, "ko": ko, "ncbi_gene_id": gene, "link_type": "enzyme"} + + +def test_hub_and_currency_compounds_are_shown_but_never_credited(): + compound_map = { + "hubby": [{"kegg_compound": "C00099", "match_basis": "exact_name"}], + "atp": [{"kegg_compound": "C00002", "match_basis": "exact_name"}], + "rare": [{"kegg_compound": "C00050", "match_basis": "exact_name"}], + } + rows = [_row("C00099", f"G{i}") for i in range(4)] + [_row("C00002", "G1"), _row("C00050", "G9")] + links, exclusions = ml.build_links(rows, compound_map, {"C00002": "currency"}, max_genes_per_compound=3) + + hub = [link for link in links if link.metabolite == "hubby"] + assert len(hub) == 4 and not any(link.credited for link in hub) + assert {(e.metabolite, e.reason) for e in exclusions} == { + ("hubby", "hub_compound"), ("atp", "currency_metabolite")} + assert not [link for link in links if link.metabolite == "atp"] + (rare,) = [link for link in links if link.metabolite == "rare"] + assert rare.credited and rare.link_confidence == "inferred" + + +def test_ambiguous_name_links_every_candidate_compound_as_inferred(): + """'lactate' from LC-MS cannot be told apart from its stereoisomer.""" + compound_map = {"lactate": [{"kegg_compound": "C00186", "match_basis": "curator_judgement"}, + {"kegg_compound": "C00256", "match_basis": "curator_judgement"}]} + links, exclusions = ml.build_links([_row("C00256", "G1"), _row("C00256", "G2", ko="K2")], + compound_map, {}) + assert {link.link_confidence for link in links} == {"inferred"} + assert [(e.kegg_compound, e.reason) for e in exclusions] == [("C00186", "no_reaction_in_table")] + + +# --------------------------------------------------------------------------- # +# Partition selection: demo and real tables are never mixed +# --------------------------------------------------------------------------- # + + +def test_simulated_oyster_partition_uses_the_demo_table_only(): + table = ml.links_for_partition("true", "29159") + assert table.available and table.path == ml.DEMO_LINKS_PATH + assert {link.metabolite for link in table.for_gene(LDH)} == {"lactate", "l_lactate", "pyruvate"} + + +def test_simulated_other_species_has_no_link_table(): + """The demo mussel reuses oyster demo IDs; it must not inherit oyster links.""" + assert ml.links_for_partition("true", "29158").status == ml.STATUS_NOT_REGISTERED + + +def test_real_partition_reports_an_unbuilt_table(monkeypatch, tmp_path): + monkeypatch.setattr(ml, "registered_link_path", lambda taxid: tmp_path / "missing.tsv") + table = ml.links_for_partition("false", "29159") + assert table.status == ml.STATUS_NOT_BUILT + assert table.label() == "not_built" + + +# --------------------------------------------------------------------------- # +# Ranking +# --------------------------------------------------------------------------- # + + +def test_linked_tier_sits_between_shared_identity_and_emerging(): + assert TIER_ORDER.index("multi_omics_convergence") < TIER_ORDER.index( + "metabolite_linked_convergence") < TIER_ORDER.index("emerging") + + +def _candidate(feature, feature_type, evidence_id, study): + return pd.DataFrame([{ + "feature_id_standardized": feature, "phenotype": "growth_under_stress", + "feature_type": feature_type, "simulated": False, "species_taxid": 29159, + "k_studies": 1, "studies": study, "n_evidence_records": 1, "total_sample_size": 10, + "pooled_effect": 1.2, "pooled_se": 0.3, "ci_lower": 0.6, "ci_upper": 1.8, "z": 4.0, + "p_value": 1e-4, "adjusted_p_value": 0.01, "n_tests_in_family": 10, + "q_statistic": 0.0, "i_squared": 0.0, "tau_squared": 0.0, "direction_consistency": 1.0, + "distinct_tissues": 1, "distinct_life_stages": 1, "distinct_stressors": "hypoxia", + "contributing_evidence_ids": evidence_id, + }]) + + +def _evidence(rows): + return pd.DataFrame([{ + "evidence_id": eid, "study_id": study, "feature_id_standardized": feature, + "phenotype": phenotype, "feature_type": ftype, "simulated": False, "species_taxid": 29159, + "mapping_confidence": "exact", "quality_flags": [], "adjusted_p_value": adj_p, + } for eid, study, feature, ftype, phenotype, adj_p in rows]) + + +@pytest.fixture +def synthetic_links(monkeypatch): + """Real partition, taxid 29159, linked through an in-memory table.""" + links, exclusions = ml.build_links( + [_row("C00186", "GENE_LDH", reaction="R00703", ko="K00016")], + {"lactate": [{"kegg_compound": "C00186", "match_basis": "exact_name"}]}, {}, + ) + table = ml.MetaboliteGeneLinks(path=Path("synthetic.tsv"), status=ml.STATUS_AVAILABLE, + sha256="0" * 64, source="test", links=links, exclusions=exclusions) + monkeypatch.setattr(rank, "links_for_partition", lambda simulated, taxid: table) + monkeypatch.setattr(evidence_cards, "links_for_partition", lambda simulated, taxid: table) + return table + + +def test_gene_and_metabolite_converge_through_a_credited_link(synthetic_links): + evidence = _evidence([ + ("rna", "RNA", "GENE_LDH", "gene", "growth_under_stress", 0.01), + ("met", "MET", "lactate", "metabolite_feature", "growth_under_stress", 0.02), + ]) + gene = build_candidates(_candidate("GENE_LDH", "gene", "rna", "RNA"), evidence).iloc[0] + metabolite = build_candidates(_candidate("lactate", "metabolite_feature", "met", "MET"), evidence).iloc[0] + + for candidate, linked_layer in [(gene, "metabolomics"), (metabolite, "transcriptomics")]: + assert candidate["tier"] == "metabolite_linked_convergence" + assert candidate["convergence_basis"] == "metabolite_gene_link" + assert candidate["linked_layers"] == linked_layer + # A link is never counted as a shared-identifier layer. + assert candidate["n_supporting_layers"] == 1 + assert candidate["linked_features"] == "lactate→GENE_LDH via C00186/R00703/K00016" + assert candidate["link_table_sha256"] == "0" * 64 + + +def test_link_does_not_change_the_score(synthetic_links, monkeypatch): + evidence = _evidence([ + ("rna", "RNA", "GENE_LDH", "gene", "growth_under_stress", 0.01), + ("met", "MET", "lactate", "metabolite_feature", "growth_under_stress", 0.02), + ]) + linked = build_candidates(_candidate("GENE_LDH", "gene", "rna", "RNA"), evidence).iloc[0] + monkeypatch.setattr(rank, "links_for_partition", + lambda simulated, taxid: ml.MetaboliteGeneLinks(path=None, status="not_registered")) + unlinked = build_candidates(_candidate("GENE_LDH", "gene", "rna", "RNA"), evidence).iloc[0] + assert unlinked["tier"] == "emerging" + assert linked["score"] == unlinked["score"] + + +@pytest.mark.parametrize("phenotype,adj_p", [("thermal_tolerance", 0.01), ("growth_under_stress", 0.3)]) +def test_linked_evidence_must_be_significant_and_same_phenotype(synthetic_links, phenotype, adj_p): + evidence = _evidence([ + ("rna", "RNA", "GENE_LDH", "gene", "growth_under_stress", 0.01), + ("met", "MET", "lactate", "metabolite_feature", phenotype, adj_p), + ]) + candidate = build_candidates(_candidate("GENE_LDH", "gene", "rna", "RNA"), evidence).iloc[0] + assert candidate["tier"] == "emerging" + assert candidate["convergence_basis"] == "none" + + +def test_candidate_must_carry_its_own_signal(synthetic_links): + evidence = _evidence([ + ("rna", "RNA", "GENE_LDH", "gene", "growth_under_stress", 0.4), + ("met", "MET", "lactate", "metabolite_feature", "growth_under_stress", 0.01), + ]) + meta = _candidate("GENE_LDH", "gene", "rna", "RNA").assign(p_value=0.3, adjusted_p_value=0.4) + candidate = build_candidates(meta, evidence).iloc[0] + assert candidate["tier"] == "emerging" + + +def test_shared_identity_outranks_a_link(synthetic_links): + evidence = _evidence([ + ("rna", "RNA", "GENE_LDH", "gene", "growth_under_stress", 0.01), + ("meth", "METH", "GENE_LDH", "methylation_region", "growth_under_stress", 0.01), + ("met", "MET", "lactate", "metabolite_feature", "growth_under_stress", 0.01), + ]) + candidate = build_candidates(_candidate("GENE_LDH", "gene", "rna", "RNA"), evidence).iloc[0] + assert candidate["tier"] == "multi_omics_convergence" + assert candidate["convergence_basis"] == "shared_identifier" + # The metabolite layer is still reported as linked context, not as shared identity. + assert candidate["supporting_layers"] == "dna_methylation|transcriptomics" + assert candidate["linked_layers"] == "metabolomics" + + +# --------------------------------------------------------------------------- # +# Demo end to end: GIGAS_GROW06 (metabolomics) x GIGAS_HEAT01 chronic (RNA-seq) +# --------------------------------------------------------------------------- # + + +def _demo_growth_candidates(): + harmonize_all_demo_studies() + from reporting.evidence_cards import rank_all_candidates + candidates = rank_all_candidates(load_evidence_table(), "growth_under_stress", None) + return oyster(candidates).set_index("feature_id_standardized") + + +def test_demo_lactate_and_ldh_converge_by_link_only(isolated_reports): + candidates = _demo_growth_candidates() + for feature in (LDH, "lactate"): + assert candidates.loc[feature, "tier"] == "metabolite_linked_convergence" + assert candidates.loc[feature, "link_confidences"] == "inferred" + # sdha is linked to a significant succinate change, but is not significant itself. + assert candidates.loc[SDHA, "tier"] == "emerging" + assert candidates.loc["succinate", "tier"] == "emerging" + assert candidates.loc["ATP", "tier"] == "emerging" + assert not (candidates["tier"] == "multi_omics_convergence").any() + + +def test_demo_cards_explain_the_link_and_its_exclusions(isolated_reports): + harmonize_all_demo_studies() + index = evidence_cards.build_evidence_cards(phenotype="growth_under_stress") + cards = {(e["feature_id_standardized"], int(e["species_taxid"])): Path(e["card_file"]).read_text() + for e in index} + lactate = cards[("lactate", 29159)] + assert "not by a shared identifier" in lactate + assert "Metabolite–gene links (biochemical; weaker than shared gene identity)" in lactate + assert "metabolite_gene_links_demo.tsv (sha256 " in lactate + assert "| GIGAS_HEAT01 | transcriptomics | LOC105342001 | C00186 | R00703 | K00016 |" in lactate + assert "different stressor" in lactate + assert "`C00256`: no reaction" in lactate + + atp = cards[("ATP", 29159)] + assert "`C00002`: currency metabolite" in atp + assert "not by a shared identifier" not in atp + + candidates = pd.read_csv(isolated_reports["cards_dir"] / "candidates.tsv", sep="\t") + for column in ("convergence_basis", "linked_layers", "linked_features", "link_confidences", + "link_table_status", "link_table_sha256"): + assert column in candidates.columns + + +# --------------------------------------------------------------------------- # +# Builder (no network: a fake KEGG) +# --------------------------------------------------------------------------- # + + +def test_equation_side(): + assert equation_side("C00186 + C00003 <=> C00022 + C00004", "C00186") == "left" + assert equation_side("C00186 + C00003 <=> C00022 + C00004", "C00022") == "right" + assert equation_side("C00031 <=> C00031", "C00031") == "both" + assert equation_side("C00031 <=> C00022", "C00186") is None + + +def test_parse_entries_reads_first_value_per_entry(): + text = ("ENTRY R00703 Reaction\nNAME (S)-lactate:NAD+ oxidoreductase\n" + "EQUATION C00186 + C00003 <=> C00022 + C00004 + C00080\n///\n" + "ENTRY C00186 Compound\nNAME (S)-Lactate;\n" + " L-Lactate;\n///\n") + parsed = parse_entries(text, ("EQUATION", "NAME")) + assert parsed["R00703"]["EQUATION"].startswith("C00186 + C00003") + assert parsed["C00186"]["NAME"] == "(S)-Lactate;" + + +FAKE_KEGG = { + "info/kegg": "kegg\tKEGG (Kyoto Encyclopedia of Genes and Genomes)\n\tko 28,545 2026/10/16\n", + "link/ko/crg": ("crg:1001\tko:K03778\ncrg:1002\tko:K03778\ncrg:9999\tko:K00234\n" + "crg:tRNA-Arg\tko:K00234\ncrg:1003\tko:K00873\n"), + "link/reaction/cpd:C00256": "cpd:C00256\trn:R00704\n", + "link/reaction/cpd:C00042": "cpd:C00042\trn:R02164\n", +} + + +def _fake_fetch(calls): + def fetch(operation): + calls.append(operation) + if operation in FAKE_KEGG: + return FAKE_KEGG[operation] + if operation.startswith("link/reaction/"): + return "" + if operation.startswith("link/ko/"): + return "rn:R00704\tko:K03778\nrn:R02164\tko:K00234\n" + if operation.startswith("get/rn:"): + return ("ENTRY R00704\nEQUATION C00256 + C00003 <=> C00022 + C00004\n///\n" + "ENTRY R02164\nEQUATION C15602 + C00042 <=> C15603 + C00122\n///\n") + if operation.startswith("get/cpd:"): + return "ENTRY C00256\nNAME (R)-Lactate;\n///\nENTRY C00042\nNAME Succinate;\n///\n" + raise AssertionError(f"unexpected KEGG call {operation}") + return fetch + + +def test_builder_writes_a_resolved_table_with_provenance(tmp_path): + calls = [] + out = tmp_path / "links.tsv" + # 9999 is a retired GeneID replaced by 1004; "tRNA-Arg" resolves to nothing. + resolver = ({"1001", "1002", "1003", "1004"}, {"9999": "1004"}, tmp_path / "crosswalk.tsv") + build(taxid=29159, out_path=out, fetch=_fake_fetch(calls), cache_dir=tmp_path / "cache", + gene_resolver=resolver) + + table = pd.read_csv(out, sep="\t", comment="#", dtype=str) + assert set(table.columns) == set(ml.LINK_COLUMNS) + assert set(zip(table["kegg_compound"], table["ncbi_gene_id"])) == { + ("C00256", "1001"), ("C00256", "1002"), ("C00042", "1004")} + assert set(table.loc[table["kegg_compound"] == "C00256", "equation_side"]) == {"left"} + # Currency compounds are never queried. + assert "link/reaction/cpd:C00002" not in calls + + prov = load_json(ml.provenance_path(out)) + jsonschema.validate(prov, load_json(SCHEMAS_DIR / "metabolite_gene_links_provenance.schema.json")) + assert prov["artifact_sha256"] == sha256_file(out) + assert prov["coverage"]["kegg_genes_retired_remapped"] == 1 + assert prov["unresolved_kegg_gene_ids"] == ["tRNA-Arg"] + assert "C00002" in prov["currency_metabolites"]["skipped"] + assert all(len(s["sha256"]) == 64 for s in prov["sources"]) + + # Offline rebuild from the cache reproduces the same table. + out2 = tmp_path / "links2.tsv" + build(taxid=29159, out_path=out2, offline=True, cache_dir=tmp_path / "cache", gene_resolver=resolver, + fetch=lambda op: (_ for _ in ()).throw(AssertionError("network used offline"))) + body = [line for line in out.read_text().splitlines() if not line.startswith("#")] + assert body == [line for line in out2.read_text().splitlines() if not line.startswith("#")] + assert json.loads(ml.provenance_path(out2).read_text())["artifact_sha256"] == sha256_file(out2)