Metabolite → gene links: metabolite_linked_convergence tier

.gitignore

@@ -58,3 +58,10 @@ registry/study_registry.csv
 # OrthoFinder intermediates (staged proteomes, DIAMOND results); the run's
 # outputs under data/reference/orthogroups/<run_id>/ 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/

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
 

README.md

@@ -151,6 +151,7 @@ AREE/
 | `aree build-evidence-cards [--phenotype <p>] [--feature-type <t>] [--max-adjusted-p <q>] [--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 <N>] [--candidates <path>] [--out <path>]` | 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 <n> [--slug <s>]` | 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 <acc> [--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 <path> [--threads N] [--activate]` | Infer orthogroups across every species with a crosswalk (OrthoFinder on RefSeq proteomes); writes `data/reference/orthogroups/<run_id>/` |
 | `aree build-dashboard [--out <path>]` | Summarize registry, pipeline, and findings state into `docs/dashboard/data.json` (plus the top candidates' evidence cards) for the Quarto dashboard page |

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**."
             )
 

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

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"}

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

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:<id> 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	

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."
+  ]
+}

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

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

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"

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)

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 |

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`<span title="${v}">${v.slice(0, n - 1).trimEnd()}…</span>` : v ?? "")
 // "Calla B, Thompson NF, Burge CA (2026)" -> "Calla et al. (2026)"
 shortCite = (c) => {
@@ -129,7 +129,8 @@ html`<div class="kpi-grid">
 
   const tierText = html`<p>${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(", ")}.</p>`;
 
   const caveats = [];

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.