From e5f10cc676984a07fa99eafda609d058acbdc9cf Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Fri, 16 Oct 2026 04:34:25 -0700 Subject: [PATCH 1/2] Add cross-species candidates and evidence cards: aree build-cross-species-cards (WP6) Tiers orthogroup-level results behind five named gates (replication, direction, significance, phenotype relevance, quality) into cross_species_convergence / emerging / discordant / unsupported, with a transparent weighted score and a between-species I-squared penalty. Cards show every species side by side. Single-species cards now name their orthogroup. Fixes the dashboard pool summary reading orthogroup results as gene-level pools. Co-Authored-By: Claude Opus 5.5 --- README.md | 1 + docs/implementation_status.md | 6 +- docs/interpreting_candidate_scores.md | 69 ++++++ docs/interpreting_meta_analysis.md | 4 + docs/mussel_clam_integration_plan.md | 19 +- src/aree/cli.py | 37 +++ src/meta_analysis/orthogroups.py | 16 ++ src/prioritize/cross_species.py | 198 ++++++++++++++++ src/prioritize/rank.py | 3 + src/reporting/cross_species_cards.py | 313 +++++++++++++++++++++++++ src/reporting/dashboard.py | 5 +- src/reporting/evidence_cards.py | 20 +- tests/test_cross_species_candidates.py | 171 ++++++++++++++ 13 files changed, 848 insertions(+), 14 deletions(-) create mode 100644 src/prioritize/cross_species.py create mode 100644 src/reporting/cross_species_cards.py create mode 100644 tests/test_cross_species_candidates.py diff --git a/README.md b/README.md index 129599e..e3a1f5b 100644 --- a/README.md +++ b/README.md @@ -147,6 +147,7 @@ AREE/ | `aree check-study-scope --branch ` | Fail if a study branch changes files its study does not own (run by CI) | | `aree harmonize --study [--input ]` | Harmonize a study (or one processed table) into the evidence table | | `aree meta-analyze [--phenotype

] [--feature-type ] [--level gene\|orthogroup] [--include-multi-copy]` | Random-effects meta-analysis over the evidence table, with BH-adjusted p-values per phenotype / feature-type family; `--level orthogroup` pools species' gene-level estimates across species through single-copy orthogroups | +| `aree build-cross-species-cards [--phenotype

] [--feature-type ] [--max-adjusted-p ] [--all-cards] [--include-multi-copy]` | Tier orthogroup-level (cross-species) results into `reports/evidence_cards/cross_species/candidates.tsv` and write a card per significant row showing every species side by side | | `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 | diff --git a/docs/implementation_status.md b/docs/implementation_status.md index 3ca1f8c..e2d29f9 100644 --- a/docs/implementation_status.md +++ b/docs/implementation_status.md @@ -76,8 +76,10 @@ worth less than one that says what it got wrong. - Cross-species pooling. Orthogroups exist (OrthoFinder run `orthofinder3.1.5_2026-10-13`, attached to evidence at harmonization) and orthogroup-level meta-analysis exists (`aree meta-analyze --level - orthogroup`), but cross-species candidates, their evidence cards and the - dashboard gate are not built; see WP6 in mussel_clam_integration_plan.md. + orthogroup`), as do cross-species candidates and their evidence cards + (`aree build-cross-species-cards`), but the dashboard gate is not built and + no real cross-species row exists until a second species' real evidence is + on `main`; see WP6 in mussel_clam_integration_plan.md. - Improving UniProt coverage beyond the 8.4% of genes that UniProtKB itself cross-references to NCBI GeneID — most likely by mapping RefSeq protein accessions through the NCBI annotation rather than relying on UniProt xrefs. diff --git a/docs/interpreting_candidate_scores.md b/docs/interpreting_candidate_scores.md index 956bc9e..2b07c95 100644 --- a/docs/interpreting_candidate_scores.md +++ b/docs/interpreting_candidate_scores.md @@ -209,6 +209,75 @@ BH-adjusted p-values with the family size, and a forest-style table of per-study effect estimates with each study's own adjusted p, so a reviewer can see exactly which inputs drove the score rather than trusting a single number. +## Cross-species candidates + +`aree build-cross-species-cards` tiers the orthogroup-level meta-analysis +(see [interpreting_meta_analysis.md](interpreting_meta_analysis.md#cross-species-results---level-orthogroup)): +one row per orthogroup, phenotype and feature type, pooling each species' own +gene-level estimate. The rule is the same as for single-species candidates: +**hard gates first, score second.** It writes, under +`reports/evidence_cards/cross_species/`, a `candidates.tsv` with every row and +its component breakdown, an `index.json`, and a card for each row whose pooled +adjusted p, or any one species' gene-level adjusted p, is at or below +`--max-adjusted-p`. + +The five convergence gates, each named in the `gates_failed` column when it +fails: + +| Gate | Passes when | +|---|---| +| `replication` | at least two species each have k ≥ 2 studies of their own | +| `direction` | every species' estimate has the same sign | +| `significance` | the pooled cross-species BH-adjusted p ≤ 0.05 | +| `phenotype_relevance` | the weakest evidence class across all species' contributing records is more than exposure-only | +| `quality` | quality flags across contributing records are within the gene-level limit | + +The tiers: + +| Tier | Meaning | +|---|---| +| `cross_species_convergence` | every gate passes | +| `cross_species_emerging` | species agree in sign and the pooled estimate is significant, but another gate fails, usually replication | +| `cross_species_discordant` | species' estimates disagree in sign; shown, never dropped | +| `cross_species_unsupported` | species agree in sign but the pooled estimate is not significant | + +Sign agreement is judged on point estimates, so a near-zero estimate of the +opposite sign makes a row discordant. Each card prints every species' +confidence interval so a reader can tell a reversal from noise around zero. + +The cross-species score uses its own weights, which sum to 100: + +| Component | Weight | What it measures | +|---|---|---| +| `per_species_replication_score` | 20 | the smallest per-species study count, saturating at 2 = 0.5 | +| `n_species_score` | 15 | number of species, saturating at 2 = 0.5 | +| `direction_consistency_score` | 15 | share of species agreeing with the majority sign | +| `effect_magnitude_score` | 10 | absolute pooled cross-species effect | +| `significance_score` | 10 | pooled BH-adjusted p, on the −log10 scale | +| `phenotype_relevance_score` | 10 | weakest evidence class across all species' records | +| `orthology_score` | 8 | copy class: `single_copy_all` 1.0, `single_copy` and `demo_crosswalk` 0.9, `multi_copy` 0.3 | +| `sample_size_score` | 5 | total biological replicates across species | +| `mapping_confidence_score` | 4 | worst identifier mapping behind any species' estimate | +| `quality_score` | 3 | quality flags across all contributing records | + +Between-species I² is subtracted as a penalty of up to 15 points. The +smallest per-species study count carries the most weight because a +cross-species claim is only as strong as its weakest species. + +A cross-species card shows each species' gene, original identifiers, studies, +tissue, life stage, stressor, sample size, effect with 95% CI and gene-level +adjusted p side by side, then the pooled estimate and between-species +heterogeneity. It names the orthogroup run, or says the orthogroup is the demo +crosswalk's illustrative one. It also warns when the species' studies differ +in stressor, tissue or life stage, because species and study design are then +confounded. Agreement across species shows an association recurs in +orthologous genes. It does not show the genes do the same thing in each +species. + +On the demo, oyster and mussel larval studies give five emerging rows, one +discordant row and two unsupported rows. None converges, because the demo +mussel has one study per gene. + ## Related documentation - [interpreting_meta_analysis.md](interpreting_meta_analysis.md) diff --git a/docs/interpreting_meta_analysis.md b/docs/interpreting_meta_analysis.md index 8e38321..99d3493 100644 --- a/docs/interpreting_meta_analysis.md +++ b/docs/interpreting_meta_analysis.md @@ -196,6 +196,10 @@ As of 2026-10-14, the only cross-species rows come from the demo: oyster and mussel larval studies sharing demo orthogroups, none past the gate. Real cross-species rows need real evidence from a second species on `main`. +`aree build-cross-species-cards` tiers these rows and writes a card for each +one that shows every species side by side; see +[interpreting_candidate_scores.md](interpreting_candidate_scores.md#cross-species-candidates). + ## Related documentation - [interpreting_candidate_scores.md](interpreting_candidate_scores.md) diff --git a/docs/mussel_clam_integration_plan.md b/docs/mussel_clam_integration_plan.md index 303c823..56b98d4 100644 --- a/docs/mussel_clam_integration_plan.md +++ b/docs/mussel_clam_integration_plan.md @@ -2,8 +2,8 @@ **Status, 2026-10-13:** WP0 to WP4 are implemented. WP5 has its backlog and a first study (TUCCI2025_HEAT_PRIMING, Manila clam, reanalyzed). WP6 has its -orthogroup run and orthogroup-level meta-analysis; cross-species candidates -and the dashboard gate are not built yet (see the entries in section 5). Everything else is a proposal. This plan describes the work needed +orthogroup run, orthogroup-level meta-analysis and cross-species candidates +with their evidence cards; the dashboard gate is not built yet (see the entries in section 5). Everything else is a proposal. This plan describes the work needed to take AREE from a single-species (*Magallana gigas*) evidence engine to one that holds mussel and clam evidence, first species-by-species and then across species through orthology. @@ -442,7 +442,7 @@ k ≥ 2 per taxon. ### WP6. Cross-species orthology and pooling (the new science) -**Status: orthogroups done, 2026-10-13; pooling not started.** What exists: +**Status: orthogroups done, 2026-10-13; pooling and cross-species candidates done, 2026-10-14; dashboard gate not started.** What exists: - `aree fetch-reference --proteins` writes one protein per gene (the longest) for each species, keyed by NCBI GeneID; all seven reference assemblies have @@ -478,8 +478,17 @@ from the active run. On the demo, 8 orthogroups pool oyster with mussel Building it found that the demo's illustrative orthogroup IDs collide with real OrthoFinder IDs, so simulated evidence never consults the run (tested). -Still to build: the cross-species candidate category, evidence cards that -show both species, and the dashboard gate. The +**Cross-species candidates and cards: done, 2026-10-14.** `aree +build-cross-species-cards` tiers each orthogroup row behind five named gates +(replication, direction, significance, phenotype relevance, quality) into +`cross_species_convergence`, `cross_species_emerging`, +`cross_species_discordant` or `cross_species_unsupported`, with a transparent +weighted score and a between-species I-squared penalty. Each card shows every +species' estimate side by side, the pooled estimate, the orthogroup run and +the gates failed. Single-species cards now name their orthogroup. On the +demo: 5 emerging, 1 discordant, 2 unsupported, none convergent. + +Still to build: the dashboard gate. The OrthoFinder run used `-M dendroblast` (OrthoFinder 3 has no orthogroups-only mode) and only its `Orthogroups.tsv` is used; hierarchical orthogroups are a possible refinement. diff --git a/src/aree/cli.py b/src/aree/cli.py index c911251..ad039fb 100644 --- a/src/aree/cli.py +++ b/src/aree/cli.py @@ -469,6 +469,43 @@ def build_evidence_cards_cmd(phenotype, feature_type, max_adjusted_p, all_cards) click.echo(f" {tier}: {tiers[tier]}") +@main.command("build-cross-species-cards") +@click.option("--phenotype", default=None, help="Phenotype ontology term id to filter on (default: all).") +@click.option("--feature-type", default=None, help="Feature type to filter on (default: all).") +@click.option("--max-adjusted-p", default=DEFAULT_MAX_ADJUSTED_P, show_default=True, type=float, + help="Render a card only when the pooled cross-species adjusted p, or any one species' " + "gene-level adjusted p, is at or below this value. Every row is still listed in " + "reports/evidence_cards/cross_species/candidates.tsv.") +@click.option("--all-cards", is_flag=True, help="Render a card for every orthogroup row.") +@click.option("--include-multi-copy", is_flag=True, + help="Also pool multi-copy orthogroups where each species has evidence for exactly one gene.") +def build_cross_species_cards_cmd(phenotype, feature_type, max_adjusted_p, all_cards, include_multi_copy): + """Tier orthogroup-level (cross-species) results and write cross-species evidence cards.""" + from reporting.cross_species_cards import TIER_ORDER as CROSS_TIERS + from reporting.cross_species_cards import build_cross_species_cards_report + + try: + report = build_cross_species_cards_report( + phenotype=phenotype, feature_type=feature_type, max_adjusted_p=max_adjusted_p, + all_cards=all_cards, include_multi_copy=include_multi_copy) + except (FileNotFoundError, ValueError) as exc: + click.echo(click.style(str(exc), fg="red")) + sys.exit(1) + if report.n_candidates == 0: + click.echo("No orthogroup has gene-level estimates from two or more species, so there are no " + "cross-species candidates. Harmonize a study of a second species first.") + return + click.echo(click.style( + f"Tiered {report.n_candidates} cross-species rows ({report.candidates_path}); " + f"wrote {len(report.index)} cards to {report.candidates_path.parent}/", fg="green")) + counts = {} + for row in report.index: + counts[row["tier"]] = counts.get(row["tier"], 0) + 1 + for tier in CROSS_TIERS: + if tier in counts: + click.echo(f" {tier}: {counts[tier]}") + + @main.command("top-candidates") @click.option("--n", "n", default=10, show_default=True, type=int, help="Number of top-ranked candidates to show per phenotype.") diff --git a/src/meta_analysis/orthogroups.py b/src/meta_analysis/orthogroups.py index d1b027d..ffea3e9 100644 --- a/src/meta_analysis/orthogroups.py +++ b/src/meta_analysis/orthogroups.py @@ -55,7 +55,15 @@ "orthogroup_id", "phenotype", "feature_type", "simulated", "species_taxid", "feature_id_standardized", "k_studies", "studies", "pooled_effect", "pooled_se", "p_value", "adjusted_p_value", + # Carried from the gene-level row so cross-species prioritization can judge + # relevance, quality and context without re-running the gene-level pool. + "contributing_evidence_ids", "total_sample_size", "distinct_tissues", + "distinct_life_stages", "distinct_stressors", "mapping_confidences", "quality_flags_union", ] +# Defaults for gene-level rows built by hand (tests) without these columns. +_CARRIED_DEFAULTS = {"contributing_evidence_ids": "", "total_sample_size": 0, "distinct_tissues": 0, + "distinct_life_stages": 0, "distinct_stressors": "", "mapping_confidences": "", + "quality_flags_union": ""} def gene_to_orthogroup(evidence: pd.DataFrame) -> dict[tuple, str]: @@ -87,6 +95,13 @@ def _key(value): return text +def _carried(row: pd.Series, column: str): + value = row.get(column) + if value is None or (isinstance(value, float) and value != value): + return _CARRIED_DEFAULTS[column] + return value + + def run_orthogroup_meta_analysis( gene_results: pd.DataFrame, evidence: pd.DataFrame, @@ -166,6 +181,7 @@ def run_orthogroup_meta_analysis( "k_studies": int(r["k_studies"]), "studies": r["studies"], "pooled_effect": r["pooled_effect"], "pooled_se": r["pooled_se"], "p_value": r["p_value"], "adjusted_p_value": r.get("adjusted_p_value"), + **{c: _carried(r, c) for c in _CARRIED_DEFAULTS}, }) result = pd.DataFrame(rows) diff --git a/src/prioritize/cross_species.py b/src/prioritize/cross_species.py new file mode 100644 index 0000000..f2e4935 --- /dev/null +++ b/src/prioritize/cross_species.py @@ -0,0 +1,198 @@ +"""Cross-species candidates: orthogroups whose evidence converges across species. WP6. + +Input is the orthogroup-level meta-analysis (meta_analysis/orthogroups.py): one +row per orthogroup, phenotype and feature type, pooling each species' own +gene-level estimate. This module tiers and scores those rows the same way +prioritize/rank.py does for single-species candidates: **hard gates first, +score second**. A high score never promotes a row past a gate it fails. + +Tiers, strongest first: + +``cross_species_convergence`` + Every gate passes: the WP6 replication gate (at least two species, each + with k >= 2 studies of its own), every species' estimate has the same sign, + the pooled cross-species estimate survives BH adjustment (q <= 0.05), the + weakest contributing evidence class is more than exposure-only, and the + quality flags are within the gene-level limit. +``cross_species_emerging`` + Species agree in sign and the pooled estimate is significant, but at least + one other gate fails. Usually the replication gate: one species rests on a + single study. ``gates_failed`` names each failed gate. +``cross_species_discordant`` + The species' estimates disagree in sign. Surfaced, not dropped: an ortholog + responding in opposite directions in two species is a finding about the + biology, or about the studies, and either way a reader must see it. +``cross_species_unsupported`` + Species agree in sign but the pooled estimate is not significant. + +Agreement in sign is judged on point estimates, as at gene level, so a +near-zero estimate of opposite sign makes a row discordant. The card shows +each species' estimate and confidence interval so a reader can tell a real +reversal from noise around zero. +""" +from __future__ import annotations + +import math + +import pandas as pd + +from .rank import EVIDENCE_CLASS_BY_RELEVANCE, HIGH_PRIORITY_MIN_QUALITY_SCORE, _partition_value +from .scoring import ( + MAPPING_CONFIDENCE_SCORE, + PHENOTYPE_RELEVANCE_SCORE, + _clip01, + _saturating, + conservative_relevance, + significance_p_value, +) + +TIER_ORDER = [ + "cross_species_convergence", + "cross_species_emerging", + "cross_species_discordant", + "cross_species_unsupported", +] +CONVERGENCE_MAX_ADJUSTED_P = 0.05 +CONVERGENCE_MIN_RELEVANCE_SCORE = PHENOTYPE_RELEVANCE_SCORE["exposure_only"] # must exceed + +# Weights sum to 100; between-species heterogeneity is a separate penalty. +WEIGHTS = { + "n_species_score": 15, # more species agreeing is stronger conservation evidence + "per_species_replication_score": 20, # the weakest species' study count: a chain is its weakest link + "sample_size_score": 5, # total biological replication across species + "effect_magnitude_score": 10, # pooled cross-species effect + "significance_score": 10, # BH-adjusted pooled p, within its family + "direction_consistency_score": 15, # share of species agreeing with the majority sign + "phenotype_relevance_score": 10, # weakest evidence class among all contributing records + "orthology_score": 8, # how like-for-like the orthogroup is + "mapping_confidence_score": 4, # worst identifier mapping behind any species' estimate + "quality_score": 3, # quality flags across all contributing records +} +HETEROGENEITY_PENALTY_WEIGHT = 15 # subtracted, scaled by between-species I^2 / 100 + +# Single-copy in every species is the cleanest one-to-one comparison; single-copy +# among the species present is nearly as good; the demo crosswalk is one gene +# per species by construction; multi-copy (opt-in only) compares one paralog +# with another species' different paralog. +ORTHOLOGY_SCORE = { + "single_copy_all": 1.0, + "single_copy": 0.9, + "demo_crosswalk": 0.9, + "multi_copy": 0.3, +} + +GATES = ["replication", "direction", "significance", "phenotype_relevance", "quality"] + + +def _split(value) -> list[str]: + if value is None or (isinstance(value, float) and value != value): + return [] + return [x for x in str(value).split("|") if x] + + +def species_estimates_for(species_df: pd.DataFrame, row: dict) -> pd.DataFrame: + """The per-species estimates behind one orthogroup result row.""" + sim = _partition_value(row["simulated"]) + mask = ( + (species_df["orthogroup_id"].astype(str) == str(row["orthogroup_id"])) + & (species_df["phenotype"] == row["phenotype"]) + & (species_df["feature_type"] == row["feature_type"]) + & (species_df["simulated"].map(_partition_value) == sim) + ) + return species_df[mask] + + +def compute_components(row: dict, species: pd.DataFrame, relevance: str) -> dict: + """0-1 component scores for one orthogroup row. Pure: same inputs, same output.""" + k = pd.to_numeric(species["k_studies"], errors="coerce").fillna(0) + mapping = {m for v in species.get("mapping_confidences", []) for m in _split(v)} + flags = {f for v in species.get("quality_flags_union", []) for f in _split(v)} + n = pd.to_numeric(species.get("total_sample_size", pd.Series(dtype=float)), errors="coerce").fillna(0) + i2 = row.get("i_squared_between_species") + return { + "n_species_score": _saturating(row["n_species"], 2.0), + "per_species_replication_score": _saturating(k.min() if len(k) else 0, 2.0), + "sample_size_score": _saturating(n.sum(), 100.0), + "effect_magnitude_score": _saturating(abs(row["pooled_effect"]), 2.0), + "significance_score": _saturating(-math.log10(significance_p_value(row)), 5.0), + "direction_consistency_score": _clip01(row["direction_consistency_species"]), + "phenotype_relevance_score": PHENOTYPE_RELEVANCE_SCORE.get(relevance, 0.1), + "orthology_score": ORTHOLOGY_SCORE.get(row["copy_class"], 0.0), + "mapping_confidence_score": min((MAPPING_CONFIDENCE_SCORE.get(m, 0.0) for m in mapping), default=0.0), + "quality_score": _clip01(1 - len(flags) / 5.0), + "heterogeneity_penalty": _clip01((i2 if i2 == i2 and i2 is not None else 0.0) / 100.0), + } + + +def cross_species_score(components: dict) -> float: + positive = sum(WEIGHTS[name] * components[name] for name in WEIGHTS) + penalty = HETEROGENEITY_PENALTY_WEIGHT * components.get("heterogeneity_penalty", 0.0) + return round(max(0.0, min(100.0, positive - penalty)), 2) + + +def gates_failed(row: dict, components: dict) -> list[str]: + """Each convergence gate this row fails, in GATES order.""" + failed = [] + if not bool(row["meets_replication_gate"]): + failed.append("replication") + if row["direction_consistency_species"] < 1.0: + failed.append("direction") + if significance_p_value(row) > CONVERGENCE_MAX_ADJUSTED_P: + failed.append("significance") + if components["phenotype_relevance_score"] <= CONVERGENCE_MIN_RELEVANCE_SCORE: + failed.append("phenotype_relevance") + if components["quality_score"] < HIGH_PRIORITY_MIN_QUALITY_SCORE: + failed.append("quality") + return failed + + +def tier_for(failed: list[str]) -> str: + if not failed: + return "cross_species_convergence" + if "direction" in failed: + return "cross_species_discordant" + if "significance" in failed: + return "cross_species_unsupported" + return "cross_species_emerging" + + +def _relevance(row: dict, species: pd.DataFrame, classification_by_id: dict) -> str: + ids = [i for v in species.get("contributing_evidence_ids", []) for i in _split(v)] + return conservative_relevance(row["phenotype"], [classification_by_id.get(i) for i in ids]) + + +def build_cross_species_candidates( + og_results: pd.DataFrame, species_df: pd.DataFrame, evidence: pd.DataFrame +) -> pd.DataFrame: + """Tier and score every orthogroup result row. Every row is returned.""" + if og_results is None or len(og_results) == 0: + return pd.DataFrame() + classification_by_id = {} + if "resilience_classification" in evidence.columns: + classification_by_id = dict(zip(evidence["evidence_id"].astype(str), evidence["resilience_classification"])) + + rows = [] + for row in og_results.to_dict("records"): + species = species_estimates_for(species_df, row) + relevance = _relevance(row, species, classification_by_id) + components = compute_components(row, species, relevance) + failed = gates_failed(row, components) + rows.append({ + **row, + "tier": tier_for(failed), + "score": cross_species_score(components), + "gates_failed": "|".join(failed), + "evidence_class": EVIDENCE_CLASS_BY_RELEVANCE.get(relevance, "exposure_only"), + "min_k_studies_per_species": int(pd.to_numeric(species["k_studies"]).min()) if len(species) else 0, + "total_sample_size": int(pd.to_numeric(species["total_sample_size"], errors="coerce").fillna(0).sum()) + if "total_sample_size" in species.columns else 0, + "mapping_confidences": "|".join(sorted({m for v in species.get("mapping_confidences", []) + for m in _split(v)})), + "quality_flags_union": "|".join(sorted({f for v in species.get("quality_flags_union", []) + for f in _split(v)})), + **{f"component_{k}": v for k, v in components.items()}, + }) + out = pd.DataFrame(rows) + rank = out["tier"].map({t: i for i, t in enumerate(TIER_ORDER)}) + return (out.assign(_rank=rank).sort_values(["_rank", "score"], ascending=[True, False]) + .drop(columns="_rank").reset_index(drop=True)) diff --git a/src/prioritize/rank.py b/src/prioritize/rank.py index 40ff50c..f0f2f38 100644 --- a/src/prioritize/rank.py +++ b/src/prioritize/rank.py @@ -234,6 +234,8 @@ def build_candidates(meta_df: pd.DataFrame, evidence_df: pd.DataFrame) -> pd.Dat relevance = conservative_relevance(meta_dict["phenotype"], classifications) meta_dict["n_supporting_layers"] = n_supporting_layers + orthogroups = (sorted(subset["orthogroup_id"].dropna().astype(str).unique()) + if "orthogroup_id" in subset.columns else []) components = compute_components(meta_dict, mapping_confidences, quality_flags, relevance) score = candidate_score(components) @@ -265,6 +267,7 @@ def build_candidates(meta_df: pd.DataFrame, evidence_df: pd.DataFrame) -> pd.Dat "pooled_stressors": "|".join(pooled_stressors(meta_dict)), "is_high_priority": is_high_priority, "is_multi_omics_convergence": is_multi_omics, + "orthogroup_id": "|".join(orthogroups), "molecular_layer": own_layer, "n_supporting_layers": n_supporting_layers, "supporting_layers": "|".join(sorted(supporting)), diff --git a/src/reporting/cross_species_cards.py b/src/reporting/cross_species_cards.py new file mode 100644 index 0000000..37fd554 --- /dev/null +++ b/src/reporting/cross_species_cards.py @@ -0,0 +1,313 @@ +"""Cross-species evidence cards: one card per orthogroup, phenotype and feature type. WP6. + +`aree build-cross-species-cards` runs the orthogroup-level meta-analysis, tiers +and scores every row (prioritize/cross_species.py), and writes, under +reports/evidence_cards/cross_species/: + +* ``candidates.tsv``: every row, every tier, with its component breakdown and + the gates it failed, so the full ranking is auditable; +* ``index.json``: the rows that got a card; +* one markdown card per row with a significant signal: the pooled + cross-species adjusted p, or any one species' gene-level adjusted p, at or + below the threshold. The second clause keeps discordant rows on the page: + two species each significant in opposite directions pool to nothing. + +Each card shows both (or all) species side by side: each species' own +gene-level estimate, then the pooled cross-species estimate and how much the +species disagree. +""" +from __future__ import annotations + +import json +from dataclasses import dataclass +from pathlib import Path + +import pandas as pd + +from common import load_vocab, species_term_for_taxid +from harmonize.core import load_evidence_table +from prioritize.cross_species import ( + TIER_ORDER, + build_cross_species_candidates, + species_estimates_for, +) +from prioritize.rank import _partition_value + +from . import evidence_cards as _cards +from .evidence_cards import ( + EVIDENCE_CLASS_STATEMENT, + NEXT_STEP_BY_EVIDENCE_CLASS, + _fmt, + _join, + _slug, +) + +CROSS_SPECIES_SUBDIR = "cross_species" + +TIER_STATEMENT = { + "cross_species_convergence": ( + "Every convergence gate passes: at least two species each replicate the association in two " + "or more studies of their own, every species' estimate has the same sign, and the pooled " + "cross-species estimate survives false-discovery control." + ), + "cross_species_emerging": ( + "The species agree in sign and the pooled estimate is significant, but at least one " + "convergence gate fails (listed below). Treat it as a lead that needs replication, not as a " + "conserved association." + ), + "cross_species_discordant": ( + "The species' estimates **disagree in sign**. This is shown, not dropped: it may be a real " + "difference between species, or a difference between the studies (stressor, tissue, life " + "stage, phenotype definition). The per-species table is the place to tell which." + ), + "cross_species_unsupported": ( + "The species agree in sign, but pooled across species the estimate is not significant." + ), +} + +NEXT_STEP_BY_TIER = { + "cross_species_convergence": ( + "Test the association in an independent cohort of each species, then take the species with " + "the strongest within-species evidence into functional follow-up." + ), + "cross_species_emerging": ( + "Replicate within each species that rests on a single study before treating the association " + "as conserved." + ), + "cross_species_discordant": ( + "Compare the species' study designs before reading the reversal as biology. A study of both " + "species under one matched design is what would resolve it." + ), + "cross_species_unsupported": ( + "No validation step is warranted from this evidence." + ), +} + +GATE_EXPLANATION = { + "replication": "fewer than two species have k >= 2 studies of their own", + "direction": "the species' estimates do not all share one sign", + "significance": "the pooled cross-species adjusted p is above 0.05", + "phenotype_relevance": "the weakest contributing evidence class is exposure-only", + "quality": "quality flags across contributing records exceed the limit", +} + + +@dataclass +class CrossSpeciesCardResult: + index: list + n_candidates: int + candidates_path: Path | None + summary: dict + + +def cards_dir() -> Path: + """Resolved at call time so tests can redirect EVIDENCE_CARDS_DIR.""" + return _cards.EVIDENCE_CARDS_DIR / CROSS_SPECIES_SUBDIR + + +def _species_label(taxid) -> str: + term = species_term_for_taxid(taxid) + return f"*{term['scientific_name']}*" if term else f"taxid {taxid}" + + +def _evidence_context(evidence: pd.DataFrame, ids: list[str]) -> dict: + subset = evidence[evidence["evidence_id"].astype(str).isin(set(ids))] + return {col: sorted({str(v) for v in subset[col].dropna()}) if col in subset.columns else [] + for col in ("feature_id_original", "tissue", "life_stage", "stressor", "study_id")} + + +def species_rows_for_card(candidate: dict, species_df: pd.DataFrame, evidence: pd.DataFrame) -> list[dict]: + """One dict per species, with the evidence context a reader needs to compare them.""" + rows = [] + for r in species_estimates_for(species_df, candidate).to_dict("records"): + ids = [i for i in str(r.get("contributing_evidence_ids") or "").split("|") if i] + ctx = _evidence_context(evidence, ids) + effect, se = float(r["pooled_effect"]), float(r["pooled_se"]) + rows.append({ + "species_taxid": r["species_taxid"], "species": _species_label(r["species_taxid"]), + "gene": r["feature_id_standardized"], "original_ids": ctx["feature_id_original"], + "k_studies": int(r["k_studies"]), "studies": str(r["studies"]), + "tissues": ctx["tissue"], "life_stages": ctx["life_stage"], "stressors": ctx["stressor"], + "sample_size": r.get("total_sample_size"), + "effect": effect, "se": se, "ci_lower": effect - 1.959964 * se, "ci_upper": effect + 1.959964 * se, + "adjusted_p_value": r.get("adjusted_p_value"), + }) + return rows + + +def has_significant_signal(candidate: dict, species_rows: list[dict], max_adjusted_p: float) -> bool: + q = candidate.get("adjusted_p_value") + if q is not None and q == q and q <= max_adjusted_p: + return True + return any(r["adjusted_p_value"] is not None and r["adjusted_p_value"] == r["adjusted_p_value"] + and r["adjusted_p_value"] <= max_adjusted_p for r in species_rows) + + +def _orthology_line(candidate: dict, run_ref: dict | None) -> str: + if candidate["copy_class"] == "demo_crosswalk": + return ("Orthology source: the demo crosswalk's illustrative orthogroup, one gene per species by " + "construction. It is not from an OrthoFinder run.") + run = run_ref["run"] if run_ref else "unknown (no active run)" + copy = { + "single_copy_all": "single-copy in every species in the run", + "single_copy": "single-copy in each species that has it, absent from at least one", + "multi_copy": ("**multi-copy**: at least one species has paralogs. Pooled only because each " + "species has evidence for exactly one gene; that gene may not be the true ortholog"), + }.get(candidate["copy_class"], candidate["copy_class"]) + return f"Orthology source: OrthoFinder run `{run}`. Copy class: {copy}." + + +def _count(value) -> str: + number = pd.to_numeric(value, errors="coerce") + return "n/a" if pd.isna(number) or number <= 0 else str(int(number)) + + +def _list(values: list) -> str: + return ", ".join(values) if values else "not recorded" + + +def card_markdown(candidate: dict, species_rows: list[dict], phenotype_label: str, run_ref: dict | None) -> str: + tier = candidate["tier"] + evidence_class = candidate.get("evidence_class", "exposure_only") + simulated = _partition_value(candidate["simulated"]) == "true" + failed = [g for g in str(candidate.get("gates_failed") or "").split("|") if g] + L = [] + L.append(f"# Cross-species evidence card: {candidate['orthogroup_id']} — {phenotype_label}") + L.append("") + L.append(f"**Tier:** `{tier}` **Score:** {candidate['score']} / 100") + L.append(f"**Evidence class:** `{evidence_class}` **Copy class:** `{candidate['copy_class']}` " + f"**Species:** {candidate['n_species']}") + L.append("") + L.append("> A cross-species estimate says an association recurs in orthologous genes of different " + "species. It is not a validated biomarker, and it does not show the genes do the same thing " + "in each species. See Limitations below.") + L.append("") + L.append("## What this evidence does and does not show") + L.append(TIER_STATEMENT[tier]) + L.append("") + L.append(EVIDENCE_CLASS_STATEMENT.get(evidence_class, EVIDENCE_CLASS_STATEMENT["exposure_only"])) + L.append("") + L.append("## Orthology") + L.append(f"- Orthogroup: `{candidate['orthogroup_id']}`") + L.append(f"- {_orthology_line(candidate, run_ref)}") + for r in species_rows: + L.append(f"- {r['species']} (taxid `{r['species_taxid']}`): gene `{r['gene']}`, " + f"original ID(s) {_list(r['original_ids'])}") + L.append(f"- Mapping confidence across contributing records: {candidate.get('mapping_confidences') or 'n/a'}") + L.append(f"- Data origin: `{'simulated' if simulated else 'real'}`") + L.append("") + L.append("## Per-species estimates") + L.append("Each species' estimate is its own gene-level random-effects pool, exactly as on its " + "single-species card.") + L.append("") + L.append("| species | gene | studies (k) | tissue | life stage | stressor | n | effect | 95% CI | gene-level adj. p |") + L.append("|---|---|---|---|---|---|---|---|---|---|") + for r in species_rows: + L.append(f"| {r['species']} | {r['gene']} | {r['studies'].replace('|', ', ')} ({r['k_studies']}) | " + f"{_list(r['tissues'])} | {_list(r['life_stages'])} | {_list(r['stressors'])} | " + f"{_count(r['sample_size'])} | " + f"{r['effect']:.3f} | {r['ci_lower']:.3f} to {r['ci_upper']:.3f} | " + f"{_fmt(r['adjusted_p_value'])} |") + L.append("") + L.append("## Cross-species pooled estimate") + L.append(f"- Pooled effect across species: {candidate['pooled_effect']:.3f} " + f"(95% CI {candidate['ci_lower']:.3f} to {candidate['ci_upper']:.3f})") + L.append(f"- Pooled p-value: {_fmt(candidate['p_value'])}; BH-adjusted across the " + f"{_fmt(candidate.get('n_tests_in_family'), 'd')} orthogroups tested for " + f"`{candidate['phenotype']}` / `{candidate['feature_type']}`: {_fmt(candidate.get('adjusted_p_value'))}") + L.append(f"- Heterogeneity between species: Q = {_fmt(candidate['q_between_species'])}, " + f"I² = {_fmt(candidate['i_squared_between_species'], '.1f')}%, " + f"τ² = {_fmt(candidate['tau_squared_between_species'])}") + L.append(f"- Direction consistency across species: {candidate['direction_consistency_species']:.2f}") + L.append(f"- Studies per species: {candidate['k_studies_by_species']}. Replication gate " + f"(two or more species each with k ≥ 2): " + f"{'met' if bool(candidate['meets_replication_gate']) else '**not met**'}") + if failed: + L.append(f"- Convergence gates failed: {_join([f'{g} ({GATE_EXPLANATION[g]})' for g in failed])}") + if tier == "cross_species_discordant": + L.append("- ⚠️ **Species disagree in direction.** Compare each species' confidence interval: an " + "interval spanning zero is noise around no effect, not a reversal.") + L.append("") + L.append("## Limitations") + L.append("- Association, not mechanism, in every species shown.") + if candidate["copy_class"] == "demo_crosswalk": + L.append("- The orthogroup is illustrative, assigned by the demo crosswalk rather than inferred.") + else: + L.append("- Orthology is inferred from sequence. Shared ancestry does not guarantee shared function " + "or regulation.") + designs = {k: {v for r in species_rows for v in r[k]} for k in ("tissues", "life_stages", "stressors")} + differ = [label for key, label in (("stressors", "stressor"), ("tissues", "tissue"), + ("life_stages", "life stage")) if len(designs[key]) > 1] + if differ: + L.append(f"- The species' studies differ in {_join(differ)}, so species and design are confounded " + "in the between-species comparison.") + L.append(f"- Quality flags across contributing records: {candidate.get('quality_flags_union') or 'none recorded'}") + if evidence_class != "resilience_associated": + L.append(f"- At least one contributing study did not measure a resilience outcome, so the evidence " + f"class is `{evidence_class}`. A strong tier here is not a resilience claim.") + if simulated: + L.append("- This card contains simulated demo evidence; do not cite it as real evidence.") + L.append("") + L.append("## Recommended next validation step") + phenotype_gap = NEXT_STEP_BY_EVIDENCE_CLASS.get(evidence_class) + if phenotype_gap and tier in ("cross_species_convergence", "cross_species_emerging"): + L.append(phenotype_gap) + L.append("") + L.append(f"Once a resilience outcome is available, the tier-appropriate step is: {NEXT_STEP_BY_TIER[tier]}") + else: + L.append(NEXT_STEP_BY_TIER[tier]) + L.append("") + return "\n".join(L) + + +def card_slug(candidate: dict) -> str: + origin = "simulated" if _partition_value(candidate["simulated"]) == "true" else "real" + return _slug(f"{candidate['orthogroup_id']}_{candidate['phenotype']}_{candidate['feature_type']}_{origin}") + + +def build_cross_species_cards_report( + phenotype: str | None = None, + feature_type: str | None = None, + max_adjusted_p: float = _cards.DEFAULT_MAX_ADJUSTED_P, + all_cards: bool = False, + include_multi_copy: bool = False, +) -> CrossSpeciesCardResult: + from mappings.orthogroups import active_run_ref + from meta_analysis.run import write_orthogroup_meta_analysis + + og_results, species_df, summary, _ = write_orthogroup_meta_analysis( + phenotype, feature_type, include_multi_copy=include_multi_copy) + evidence = load_evidence_table() + candidates = build_cross_species_candidates(og_results, species_df, evidence) + if len(candidates) == 0: + return CrossSpeciesCardResult(index=[], n_candidates=0, candidates_path=None, summary=summary) + + out_dir = cards_dir() + out_dir.mkdir(parents=True, exist_ok=True) + phenotype_vocab = load_vocab("phenotype_ontology") + run_ref = active_run_ref() + + card_files, written = [], [] + for candidate in candidates.to_dict("records"): + species_rows = species_rows_for_card(candidate, species_df, evidence) + if not all_cards and not has_significant_signal(candidate, species_rows, max_adjusted_p): + card_files.append("") + continue + label = phenotype_vocab.get(candidate["phenotype"], {}).get("label", candidate["phenotype"]) + path = out_dir / f"{card_slug(candidate)}.md" + path.write_text(card_markdown(candidate, species_rows, label, run_ref)) + card_files.append(str(path)) + written.append({k: candidate.get(k) for k in ( + "orthogroup_id", "phenotype", "feature_type", "simulated", "copy_class", "species_taxids", + "tier", "evidence_class", "score", "k_studies_by_species", "adjusted_p_value", "gates_failed")} + | {"card_file": str(path)}) + + candidates = candidates.assign(card_file=card_files) + candidates_path = out_dir / "candidates.tsv" + candidates.to_csv(candidates_path, sep="\t", index=False) + (out_dir / "index.json").write_text(json.dumps(written, indent=2, default=str)) + return CrossSpeciesCardResult(index=written, n_candidates=len(candidates), + candidates_path=candidates_path, summary=summary) + + +__all__ = ["TIER_ORDER", "build_cross_species_cards_report", "card_markdown", "cards_dir"] diff --git a/src/reporting/dashboard.py b/src/reporting/dashboard.py index 12e63cb..c03e573 100644 --- a/src/reporting/dashboard.py +++ b/src/reporting/dashboard.py @@ -359,7 +359,10 @@ def _pool_summary() -> list[dict]: seen = set() # A per-phenotype run and the all-phenotypes run describe the same pools; # read the specific files first so `seen` credits each pool to one of them. - paths = sorted(META_ANALYSIS_DIR.glob("*_meta_analysis.tsv"), + # Orthogroup-level (cross-species) results share the suffix but not the + # columns; they are summarized separately, not as gene-level pools. + paths = sorted((p for p in META_ANALYSIS_DIR.glob("*_meta_analysis.tsv") + if not p.name.endswith("_orthogroup_meta_analysis.tsv")), key=lambda p: (p.name.startswith("all_phenotypes"), p.name)) for path in paths: df = pd.read_csv(path, sep="\t", low_memory=False) diff --git a/src/reporting/evidence_cards.py b/src/reporting/evidence_cards.py index ed541b6..e1c2277 100644 --- a/src/reporting/evidence_cards.py +++ b/src/reporting/evidence_cards.py @@ -239,12 +239,20 @@ def _card_markdown(candidate: dict, forest_points: list, other_omics: list, phen if term is not None: lines.append(f"- Species: *{term['scientific_name']}* ({term['label']})") lines.append(f"- Species taxid: `{candidate['species_taxid']}`") - lines.append( - "- Ortholog context: none assigned. This card pools evidence from one species " - "only; the same feature ID in another species has its own card and is never " - "pooled with this one. Cross-species comparison through orthogroups is not " - "yet implemented." - ) + orthogroup = candidate.get("orthogroup_id") + if isinstance(orthogroup, str) and orthogroup: + lines.append( + f"- Ortholog context: orthogroup `{orthogroup}`. This card pools evidence from one " + "species only; evidence from another species is never pooled with this one. Where " + "another species has evidence for the same orthogroup, the cross-species estimate is " + "on its own card under `evidence_cards/cross_species/` (`aree build-cross-species-cards`)." + ) + else: + lines.append( + "- Ortholog context: no orthogroup assigned. This card pools evidence from one species " + "only; evidence from another species is never pooled with this one, and this feature " + "cannot be compared across species until it is in an orthogroup." + ) simulated = _partition_value(candidate["simulated"]) == "true" lines.append(f"- Data origin: `{'simulated' if simulated else 'real'}`") lines.append(f"- Mapping confidence across contributing studies: {candidate['mapping_confidences']}") diff --git a/tests/test_cross_species_candidates.py b/tests/test_cross_species_candidates.py new file mode 100644 index 0000000..12ac683 --- /dev/null +++ b/tests/test_cross_species_candidates.py @@ -0,0 +1,171 @@ +"""Cross-species candidates and evidence cards, WP6 (src/prioritize/cross_species.py, +src/reporting/cross_species_cards.py).""" +from __future__ import annotations + +from pathlib import Path + +import pandas as pd +import pytest +from conftest import harmonize_all_demo_studies + +from meta_analysis.orthogroups import run_orthogroup_meta_analysis +from prioritize.cross_species import ( + WEIGHTS, + build_cross_species_candidates, + compute_components, + cross_species_score, + species_estimates_for, +) +from reporting.cross_species_cards import ( + card_markdown, + has_significant_signal, + species_rows_for_card, +) + +OYSTER, MUSSEL = 29159, 29158 +CLASSES = {"OG1": "single_copy_all"} + + +def _gene(feature, taxid, effect, se, k, p, classification="resilience"): + ids = "|".join(f"E{feature}_{i}" for i in range(k)) + return {"feature_id_standardized": feature, "phenotype": "thermal_tolerance", "feature_type": "gene", + "simulated": False, "species_taxid": taxid, "k_studies": k, + "studies": "|".join(f"S{feature}_{i}" for i in range(k)), + "pooled_effect": effect, "pooled_se": se, "p_value": p, "adjusted_p_value": p, + "contributing_evidence_ids": ids, "total_sample_size": 12 * k, "distinct_tissues": 1, + "distinct_life_stages": 1, "distinct_stressors": "temperature", + "mapping_confidences": "exact", "quality_flags_union": "", "_class": classification} + + +def _rank(*genes): + genes = pd.DataFrame(genes) + evidence = pd.DataFrame([ + {"evidence_id": eid, "feature_id_standardized": g["feature_id_standardized"], + "feature_id_original": f"LOC{g['feature_id_standardized']}", "simulated": "False", + "species_taxid": str(g["species_taxid"]), "orthogroup_id": "OG1", + "resilience_classification": g["_class"], "study_id": eid.split("_")[0], + "tissue": "gill", "life_stage": "adult", "stressor": "temperature"} + for g in genes.to_dict("records") for eid in g["contributing_evidence_ids"].split("|")]) + result, species, _ = run_orthogroup_meta_analysis(genes.drop(columns="_class"), evidence, CLASSES) + return build_cross_species_candidates(result, species, evidence), species, evidence + + +def test_two_replicated_concordant_species_converge(): + ranked, _, _ = _rank(_gene("101", OYSTER, 1.2, 0.2, k=2, p=1e-4), _gene("201", MUSSEL, 0.9, 0.25, k=2, p=1e-3)) + row = ranked.iloc[0] + assert row["tier"] == "cross_species_convergence" and row["gates_failed"] == "" + assert row["evidence_class"] == "resilience_associated" + assert row["min_k_studies_per_species"] == 2 and row["total_sample_size"] == 48 + + +def test_one_study_in_one_species_is_emerging_with_the_gate_named(): + ranked, _, _ = _rank(_gene("101", OYSTER, 1.2, 0.2, k=2, p=1e-4), _gene("201", MUSSEL, 0.9, 0.25, k=1, p=1e-3)) + row = ranked.iloc[0] + assert row["tier"] == "cross_species_emerging" and row["gates_failed"] == "replication" + + +def test_exposure_only_evidence_in_one_species_blocks_convergence(): + """The weakest evidence class across all species' records decides, as at gene level.""" + ranked, _, _ = _rank(_gene("101", OYSTER, 1.2, 0.2, k=2, p=1e-4), + _gene("201", MUSSEL, 0.9, 0.25, k=2, p=1e-3, classification="exposure_only")) + row = ranked.iloc[0] + assert row["tier"] == "cross_species_emerging" and row["gates_failed"] == "phenotype_relevance" + assert row["evidence_class"] == "exposure_only" + + +def test_opposite_signs_are_discordant_and_still_get_a_card(): + ranked, species, evidence = _rank(_gene("101", OYSTER, 1.5, 0.2, k=2, p=1e-6), + _gene("201", MUSSEL, -1.5, 0.2, k=2, p=1e-6)) + row = ranked.iloc[0].to_dict() + assert row["tier"] == "cross_species_discordant" + assert "direction" in row["gates_failed"].split("|") + rows = species_rows_for_card(row, species, evidence) + # Pooled, the opposite effects cancel; the species' own significance keeps the card. + assert row["adjusted_p_value"] > 0.05 and has_significant_signal(row, rows, 0.05) + text = card_markdown(row, rows, "Thermal tolerance", {"run": "testrun"}) + assert "disagree in sign" in text and "Species disagree in direction" in text + + +def test_concordant_but_not_significant_is_unsupported(): + ranked, _, _ = _rank(_gene("101", OYSTER, 0.1, 0.5, k=2, p=0.8), _gene("201", MUSSEL, 0.05, 0.5, k=2, p=0.9)) + assert ranked.iloc[0]["tier"] == "cross_species_unsupported" + + +def test_score_is_reproducible_and_penalizes_between_species_heterogeneity(): + assert sum(WEIGHTS.values()) == 100 + agree, species, _ = _rank(_gene("101", OYSTER, 1.0, 0.1, k=2, p=1e-4), _gene("201", MUSSEL, 1.0, 0.1, k=2, p=1e-4)) + differ, _, _ = _rank(_gene("101", OYSTER, 0.4, 0.1, k=2, p=1e-4), _gene("201", MUSSEL, 1.6, 0.1, k=2, p=1e-4)) + row = agree.iloc[0].to_dict() + components = compute_components(row, species_estimates_for(species, row), "resilience") + assert cross_species_score(components) == row["score"] + assert differ.iloc[0]["component_heterogeneity_penalty"] > 0.5 + assert differ.iloc[0]["score"] < agree.iloc[0]["score"] + + +def test_card_shows_each_species_side_by_side(): + ranked, species, evidence = _rank(_gene("101", OYSTER, 1.2, 0.2, k=2, p=1e-4), + _gene("201", MUSSEL, 0.9, 0.25, k=1, p=1e-3)) + row = ranked.iloc[0].to_dict() + text = card_markdown(row, species_rows_for_card(row, species, evidence), "Thermal tolerance", + {"run": "orthofinder_testrun"}) + assert "*Magallana gigas*" in text and "*Mytilus galloprovincialis*" in text + assert "LOC101" in text and "LOC201" in text + assert "orthofinder_testrun" in text + assert "replication (fewer than two species have k >= 2" in text + assert "not met" in text + + +def test_demo_builds_cross_species_cards_without_claiming_convergence(isolated_reports): + from click.testing import CliRunner + + from aree.cli import main + + harmonize_all_demo_studies() + out = CliRunner().invoke(main, ["build-cross-species-cards", "--feature-type", "gene"]) + assert out.exit_code == 0, out.output + cards_dir = isolated_reports["cards_dir"] / "cross_species" + ranked = pd.read_csv(cards_dir / "candidates.tsv", sep="\t") + assert len(ranked) > 0 + # The demo mussel contributes one study per gene, so nothing converges. + assert "cross_species_convergence" not in set(ranked["tier"]) + assert set(ranked["copy_class"]) == {"demo_crosswalk"} + cards = sorted(cards_dir.glob("*.md")) + assert cards and len(cards) == (ranked["card_file"].fillna("") != "").sum() + for card in cards: + text = card.read_text() + assert "*Magallana gigas*" in text and "*Mytilus" in text + assert "simulated demo evidence" in text and "not from an OrthoFinder run" in text + + +def test_demo_gene_cards_name_their_orthogroup(isolated_reports): + from reporting.evidence_cards import build_evidence_cards_report + + harmonize_all_demo_studies() + report = build_evidence_cards_report(phenotype="larval_viability", feature_type="gene") + hsp70 = [r for r in report.index if r["feature_id_standardized"] == "LOC105333935"] + assert hsp70, "demo HSP70 card missing" + text = Path(hsp70[0]["card_file"]).read_text() + assert "orthogroup `OG0000001`" in text + + +def test_dashboard_skips_orthogroup_results(isolated_reports, monkeypatch): + """Regression: the gene-level pool summary globbed *_meta_analysis.tsv, which + also matches the orthogroup output, whose columns differ.""" + import reporting.dashboard as dashboard + from meta_analysis.run import write_meta_analysis, write_orthogroup_meta_analysis + + monkeypatch.setattr(dashboard, "META_ANALYSIS_DIR", isolated_reports["meta_dir"]) + harmonize_all_demo_studies() + write_meta_analysis(None, "gene") + write_orthogroup_meta_analysis(None, "gene") + pools = dashboard._pool_summary() + assert pools and all("k_studies" not in str(p.get("error", "")) for p in pools) + + +@pytest.mark.parametrize("copy_class,expected", [("single_copy_all", 1.0), ("multi_copy", 0.3)]) +def test_orthology_score_reflects_copy_class(copy_class, expected): + row = {"n_species": 2, "pooled_effect": 1.0, "p_value": 0.01, "adjusted_p_value": 0.01, + "direction_consistency_species": 1.0, "copy_class": copy_class, "i_squared_between_species": 0.0} + species = pd.DataFrame({"k_studies": [2, 2], "total_sample_size": [10, 10], + "mapping_confidences": ["exact", "exact"], "quality_flags_union": ["", ""]}) + assert compute_components(row, species, "resilience")["orthology_score"] == expected From e150a1e3ccdb15a449c0dbc62c5c97724760edf3 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Fri, 16 Oct 2026 05:02:13 -0700 Subject: [PATCH 2/2] Cross-species tiers: discordant needs significant species; judge quality per species On the first real oyster-clam pool (5,681 single-copy orthogroups), 2,880 of 2,886 sign disagreements involved a non-significant species, and the union of quality flags across species failed the quality gate for every row. Discordant now requires species each significant at gene level with opposite signs; other sign disagreements are unsupported. Quality is the weakest species' own gene-level quality. Docs record the first real cross-species results. Co-Authored-By: Claude Opus 5.5 --- docs/implementation_status.md | 7 +++-- docs/interpreting_candidate_scores.md | 32 +++++++++++++------ docs/interpreting_meta_analysis.md | 8 +++-- docs/mussel_clam_integration_plan.md | 16 ++++++++-- src/prioritize/cross_species.py | 43 +++++++++++++++++++------- src/reporting/cross_species_cards.py | 14 +++++---- tests/test_cross_species_candidates.py | 18 +++++++++++ 7 files changed, 102 insertions(+), 36 deletions(-) diff --git a/docs/implementation_status.md b/docs/implementation_status.md index e2d29f9..b292004 100644 --- a/docs/implementation_status.md +++ b/docs/implementation_status.md @@ -77,9 +77,10 @@ worth less than one that says what it got wrong. `orthofinder3.1.5_2026-10-13`, attached to evidence at harmonization) and orthogroup-level meta-analysis exists (`aree meta-analyze --level orthogroup`), as do cross-species candidates and their evidence cards - (`aree build-cross-species-cards`), but the dashboard gate is not built and - no real cross-species row exists until a second species' real evidence is - on `main`; see WP6 in mussel_clam_integration_plan.md. + (`aree build-cross-species-cards`). On real data the oyster and the first + Manila clam study pool through 5,681 single-copy orthogroups, none past the + replication gate yet (one clam study). The dashboard gate is not built; see + WP6 in mussel_clam_integration_plan.md. - Improving UniProt coverage beyond the 8.4% of genes that UniProtKB itself cross-references to NCBI GeneID — most likely by mapping RefSeq protein accessions through the NCBI annotation rather than relying on UniProt xrefs. diff --git a/docs/interpreting_candidate_scores.md b/docs/interpreting_candidate_scores.md index 2b07c95..13c5207 100644 --- a/docs/interpreting_candidate_scores.md +++ b/docs/interpreting_candidate_scores.md @@ -230,7 +230,7 @@ fails: | `direction` | every species' estimate has the same sign | | `significance` | the pooled cross-species BH-adjusted p ≤ 0.05 | | `phenotype_relevance` | the weakest evidence class across all species' contributing records is more than exposure-only | -| `quality` | quality flags across contributing records are within the gene-level limit | +| `quality` | each species' own quality flags are within the gene-level limit (the weakest species decides) | The tiers: @@ -238,12 +238,15 @@ The tiers: |---|---| | `cross_species_convergence` | every gate passes | | `cross_species_emerging` | species agree in sign and the pooled estimate is significant, but another gate fails, usually replication | -| `cross_species_discordant` | species' estimates disagree in sign; shown, never dropped | -| `cross_species_unsupported` | species agree in sign but the pooled estimate is not significant | +| `cross_species_discordant` | species that are each significant at gene level disagree in sign; shown, never dropped | +| `cross_species_unsupported` | the pooled estimate is not significant, or species disagree in sign only because one of them is noise | -Sign agreement is judged on point estimates, so a near-zero estimate of the -opposite sign makes a row discordant. Each card prints every species' -confidence interval so a reader can tell a reversal from noise around zero. +The direction gate is judged on point estimates, so any sign disagreement +fails it. Only a disagreement between species that are each significant at +gene level is called discordant. Two null estimates disagree in sign about +half the time: on the first real oyster-clam pool, 2,880 of 2,886 sign +disagreements involved a non-significant species, and only 6 were reversals +between significant ones. The cross-species score uses its own weights, which sum to 100: @@ -258,7 +261,7 @@ The cross-species score uses its own weights, which sum to 100: | `orthology_score` | 8 | copy class: `single_copy_all` 1.0, `single_copy` and `demo_crosswalk` 0.9, `multi_copy` 0.3 | | `sample_size_score` | 5 | total biological replicates across species | | `mapping_confidence_score` | 4 | worst identifier mapping behind any species' estimate | -| `quality_score` | 3 | quality flags across all contributing records | +| `quality_score` | 3 | the weakest species' own gene-level quality (`1 - n_flags/5`) | Between-species I² is subtracted as a penalty of up to 15 points. The smallest per-species study count carries the most weight because a @@ -274,9 +277,18 @@ confounded. Agreement across species shows an association recurs in orthologous genes. It does not show the genes do the same thing in each species. -On the demo, oyster and mussel larval studies give five emerging rows, one -discordant row and two unsupported rows. None converges, because the demo -mussel has one study per gene. +On the demo, oyster and mussel larval studies give five emerging rows and +three unsupported rows. None converges, because the demo mussel has one study +per gene. The demo's one sign disagreement involves a non-significant mussel +estimate, so it is unsupported, not discordant. + +On real data (2026-10-16), the three harmonized oyster thermal-tolerance studies and TUCCI2025_HEAT_PRIMING pool through 5,681 single-copy orthogroups: 64 emerging, 6 discordant, 5,611 unsupported, none convergent. Every emerging row fails three gates: +replication (one clam study), phenotype relevance (the clam comparison is +exposure-only), and quality (each species' own flags exceed the gene-level +limit, as they do for its single-species candidates). The 6 discordant rows +are oyster genes with small significant effects against clam genes with large +significant effects of opposite sign. The two studies differ in stressor +timing, so they are leads to check, not established species differences. ## Related documentation diff --git a/docs/interpreting_meta_analysis.md b/docs/interpreting_meta_analysis.md index 99d3493..c18510a 100644 --- a/docs/interpreting_meta_analysis.md +++ b/docs/interpreting_meta_analysis.md @@ -192,9 +192,11 @@ How to read a row: - `adjusted_p_value` is Benjamini-Hochberg within each phenotype, feature type and origin family of orthogroup-level tests. -As of 2026-10-14, the only cross-species rows come from the demo: oyster and -mussel larval studies sharing demo orthogroups, none past the gate. Real -cross-species rows need real evidence from a second species on `main`. +The first real cross-species rows came on 2026-10-16, once the first real +clam study reached `main`: oyster and Manila clam thermal-tolerance evidence +pools through 5,681 single-copy orthogroups (4,890 multi-copy orthogroups set +aside). None is past the gate, because the clam has one study. The demo adds +oyster and mussel larval rows sharing demo orthogroups, also below the gate. `aree build-cross-species-cards` tiers these rows and writes a card for each one that shows every species side by side; see diff --git a/docs/mussel_clam_integration_plan.md b/docs/mussel_clam_integration_plan.md index 56b98d4..bd37d10 100644 --- a/docs/mussel_clam_integration_plan.md +++ b/docs/mussel_clam_integration_plan.md @@ -474,7 +474,8 @@ Q, I-squared and tau-squared, per-species estimates and study counts, and a Multi-copy orthogroups are set aside by default, as are orthogroups missing from the active run. On the demo, 8 orthogroups pool oyster with mussel (larval viability), none past the gate; on all committed real evidence -(465,531 oyster records), nothing pools across species yet, correctly. +(465,531 oyster records), nothing pooled across species yet, correctly, until +the first clam study merged (see below). Building it found that the demo's illustrative orthogroup IDs collide with real OrthoFinder IDs, so simulated evidence never consults the run (tested). @@ -486,7 +487,18 @@ build-cross-species-cards` tiers each orthogroup row behind five named gates weighted score and a between-species I-squared penalty. Each card shows every species' estimate side by side, the pooled estimate, the orthogroup run and the gates failed. Single-species cards now name their orthogroup. On the -demo: 5 emerging, 1 discordant, 2 unsupported, none convergent. +demo: 5 emerging, 3 unsupported, none convergent. + +First real cross-species rows, 2026-10-16, after TUCCI2025_HEAT_PRIMING +merged: oyster and Manila clam thermal-tolerance evidence pools through 5,681 +single-copy orthogroups, giving 64 emerging, 6 discordant and 5,611 +unsupported rows, none convergent, with 603 cards. Real data changed two +rules. "Discordant" now requires species that are each significant at gene +level and disagree in sign: 2,880 of 2,886 sign disagreements involved a +non-significant species, which is noise, not a reversal. Quality is now +judged per species, because the union of flags across species always +exceeded the limit. PERUZZA2023_HEATWAVE (claimed, PR #64) is the second clam +study, the one that can bring the clam to k = 2. Still to build: the dashboard gate. The OrthoFinder run used `-M dendroblast` (OrthoFinder 3 has no orthogroups-only diff --git a/src/prioritize/cross_species.py b/src/prioritize/cross_species.py index f2e4935..289d0f8 100644 --- a/src/prioritize/cross_species.py +++ b/src/prioritize/cross_species.py @@ -19,16 +19,17 @@ one other gate fails. Usually the replication gate: one species rests on a single study. ``gates_failed`` names each failed gate. ``cross_species_discordant`` - The species' estimates disagree in sign. Surfaced, not dropped: an ortholog + At least two species are each significant at gene level (adjusted + p <= 0.05) and disagree in sign. Surfaced, not dropped: an ortholog responding in opposite directions in two species is a finding about the biology, or about the studies, and either way a reader must see it. ``cross_species_unsupported`` - Species agree in sign but the pooled estimate is not significant. + Everything else: the pooled estimate is not significant, or the species + disagree in sign only because at least one of them is noise around zero. -Agreement in sign is judged on point estimates, as at gene level, so a -near-zero estimate of opposite sign makes a row discordant. The card shows -each species' estimate and confidence interval so a reader can tell a real -reversal from noise around zero. +The direction gate itself is judged on point estimates, as at gene level, so +any sign disagreement fails it and is named in ``gates_failed``. Only a +disagreement between significant species is called discordant. """ from __future__ import annotations @@ -106,7 +107,10 @@ def compute_components(row: dict, species: pd.DataFrame, relevance: str) -> dict """0-1 component scores for one orthogroup row. Pure: same inputs, same output.""" k = pd.to_numeric(species["k_studies"], errors="coerce").fillna(0) mapping = {m for v in species.get("mapping_confidences", []) for m in _split(v)} - flags = {f for v in species.get("quality_flags_union", []) for f in _split(v)} + # Quality is judged per species, on that species' own gene-level flags: a + # union across species grows with every study added and would fail the + # gate however good each species' evidence was. The weakest species decides. + per_species_flags = [len(set(_split(v))) for v in species.get("quality_flags_union", [])] n = pd.to_numeric(species.get("total_sample_size", pd.Series(dtype=float)), errors="coerce").fillna(0) i2 = row.get("i_squared_between_species") return { @@ -119,7 +123,7 @@ def compute_components(row: dict, species: pd.DataFrame, relevance: str) -> dict "phenotype_relevance_score": PHENOTYPE_RELEVANCE_SCORE.get(relevance, 0.1), "orthology_score": ORTHOLOGY_SCORE.get(row["copy_class"], 0.0), "mapping_confidence_score": min((MAPPING_CONFIDENCE_SCORE.get(m, 0.0) for m in mapping), default=0.0), - "quality_score": _clip01(1 - len(flags) / 5.0), + "quality_score": min((_clip01(1 - n / 5.0) for n in per_species_flags), default=1.0), "heterogeneity_penalty": _clip01((i2 if i2 == i2 and i2 is not None else 0.0) / 100.0), } @@ -130,6 +134,13 @@ def cross_species_score(components: dict) -> float: return round(max(0.0, min(100.0, positive - penalty)), 2) +def significant_signs(species: pd.DataFrame, max_adjusted_p: float = CONVERGENCE_MAX_ADJUSTED_P) -> set[int]: + """Signs of the species estimates that are significant at gene level.""" + q = pd.to_numeric(species.get("adjusted_p_value", pd.Series(dtype=float)), errors="coerce") + effects = pd.to_numeric(species["pooled_effect"], errors="coerce") + return {1 if e > 0 else -1 for e, qq in zip(effects, q) if qq == qq and qq <= max_adjusted_p and e != 0} + + def gates_failed(row: dict, components: dict) -> list[str]: """Each convergence gate this row fails, in GATES order.""" failed = [] @@ -146,12 +157,17 @@ def gates_failed(row: dict, components: dict) -> list[str]: return failed -def tier_for(failed: list[str]) -> str: +def tier_for(failed: list[str], significant_discordance: bool) -> str: + """Discordant requires species that are each significant at gene level and + disagree in sign. A sign disagreement involving a non-significant species + is noise around zero: on the first real oyster-clam pool, 2,880 of 2,886 + sign disagreements were of that kind, and labelling them "discordant" + would bury the 6 real reversals.""" if not failed: return "cross_species_convergence" - if "direction" in failed: + if significant_discordance: return "cross_species_discordant" - if "significance" in failed: + if "direction" in failed or "significance" in failed: return "cross_species_unsupported" return "cross_species_emerging" @@ -177,9 +193,12 @@ def build_cross_species_candidates( relevance = _relevance(row, species, classification_by_id) components = compute_components(row, species, relevance) failed = gates_failed(row, components) + signs = significant_signs(species) rows.append({ **row, - "tier": tier_for(failed), + "tier": tier_for(failed, len(signs) > 1), + "n_species_significant": int((pd.to_numeric(species["adjusted_p_value"], errors="coerce") + <= CONVERGENCE_MAX_ADJUSTED_P).sum()), "score": cross_species_score(components), "gates_failed": "|".join(failed), "evidence_class": EVIDENCE_CLASS_BY_RELEVANCE.get(relevance, "exposure_only"), diff --git a/src/reporting/cross_species_cards.py b/src/reporting/cross_species_cards.py index 37fd554..955e7e6 100644 --- a/src/reporting/cross_species_cards.py +++ b/src/reporting/cross_species_cards.py @@ -56,12 +56,14 @@ "conserved association." ), "cross_species_discordant": ( - "The species' estimates **disagree in sign**. This is shown, not dropped: it may be a real " - "difference between species, or a difference between the studies (stressor, tissue, life " - "stage, phenotype definition). The per-species table is the place to tell which." + "Species that are each significant at gene level **disagree in sign**. This is shown, not " + "dropped: it may be a real difference between species, or a difference between the studies " + "(stressor, tissue, life stage, phenotype definition). The per-species table is the place to " + "tell which." ), "cross_species_unsupported": ( - "The species agree in sign, but pooled across species the estimate is not significant." + "Pooled across species the estimate is not significant, or the species disagree in sign only " + "because at least one of them is not significant on its own. Neither is cross-species support." ), } @@ -225,8 +227,8 @@ def card_markdown(candidate: dict, species_rows: list[dict], phenotype_label: st if failed: L.append(f"- Convergence gates failed: {_join([f'{g} ({GATE_EXPLANATION[g]})' for g in failed])}") if tier == "cross_species_discordant": - L.append("- ⚠️ **Species disagree in direction.** Compare each species' confidence interval: an " - "interval spanning zero is noise around no effect, not a reversal.") + L.append("- ⚠️ **Species disagree in direction**, and each disagreeing species is significant " + "on its own. Check the designs below before reading it as a species difference.") L.append("") L.append("## Limitations") L.append("- Association, not mechanism, in every species shown.") diff --git a/tests/test_cross_species_candidates.py b/tests/test_cross_species_candidates.py index 12ac683..106af8a 100644 --- a/tests/test_cross_species_candidates.py +++ b/tests/test_cross_species_candidates.py @@ -169,3 +169,21 @@ def test_orthology_score_reflects_copy_class(copy_class, expected): species = pd.DataFrame({"k_studies": [2, 2], "total_sample_size": [10, 10], "mapping_confidences": ["exact", "exact"], "quality_flags_union": ["", ""]}) assert compute_components(row, species, "resilience")["orthology_score"] == expected + + +def test_sign_disagreement_with_a_non_significant_species_is_unsupported_not_discordant(): + """Two null estimates disagree in sign half the time; that is not a reversal.""" + ranked, _, _ = _rank(_gene("101", OYSTER, 1.5, 0.2, k=2, p=1e-6), _gene("201", MUSSEL, -0.1, 0.5, k=2, p=0.8)) + row = ranked.iloc[0] + assert row["tier"] == "cross_species_unsupported" + assert "direction" in row["gates_failed"].split("|") and row["n_species_significant"] == 1 + + +def test_quality_is_judged_per_species_not_on_the_union(): + """Each species carries two flags (four distinct in all): every species passes on its own.""" + a = _gene("101", OYSTER, 1.2, 0.2, k=2, p=1e-4) | {"quality_flags_union": "low_replication|single_timepoint"} + b = _gene("201", MUSSEL, 0.9, 0.25, k=2, p=1e-3) | {"quality_flags_union": "tank_level_replication|unbalanced_design"} + ranked, _, _ = _rank(a, b) + row = ranked.iloc[0] + assert row["component_quality_score"] == pytest.approx(0.6) + assert row["tier"] == "cross_species_convergence"