From 2ad2266855febe77234818b59a5f48c6225d49a0 Mon Sep 17 00:00:00 2001 From: James Kane Date: Mon, 29 Jun 2026 04:25:30 -0500 Subject: [PATCH 01/13] feat(age): FTDNA STR import, per-branch STR ASR, causality + sparse-node fixes MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Pragmatic (baseline) STR-aging integration: - FTDNA Y-STR import: genomics.biosample_str_profile (mig 0053), du_db::str_import parser/upsert, du-jobs `run-once ftdna-str` loader (kit→subject_id→sample_guid, manifest token match, FTDNA_STR_ALIAS override for own-genome accessions like WGS229=B5163). build_str_inputs UNIONs local + federated profiles. - Per-branch STR ASR: tree.haplogroup_str_asr (mig 0054) persists the parsimony-reconstructed ancestral motif; ystr::branch_str_asr diffs node vs parent; Y-tree node sidebar shows the per-marker mutations (en/es/fr). - Age correctness vs McDonald 2021: * §2.3 causality back-correction on COMBINED (parent ≥ child projection) — fixes parent-younger-than-child inversions introduced by the STR term. * §2.5.2 STR sparse-node gating (MIN_STR_TESTERS_FOR_COMBINE) — keeps reconstruction-collapsed (0–1 tester) nodes from dominating the SNP clock. * clear stale denormalised tmrca_ybp/formed_ybp on undatable nodes. - documents/proposals/aging-pipeline-audit-mcdonald2021.md: full audit of the pipeline against McDonald 2021 (faithful items, divergences, refinements P1–P6). Co-Authored-By: Claude Opus 4.8 (1M context) --- .../aging-pipeline-audit-mcdonald2021.md | 253 ++++++++++++++++++ rust/Cargo.lock | 1 + rust/crates/du-db/src/age.rs | 70 ++++- rust/crates/du-db/src/lib.rs | 1 + rust/crates/du-db/src/str_import.rs | 209 +++++++++++++++ rust/crates/du-db/src/ystr.rs | 141 +++++++++- rust/crates/du-jobs/Cargo.toml | 1 + rust/crates/du-jobs/src/ftdna_str.rs | 245 +++++++++++++++++ rust/crates/du-jobs/src/main.rs | 13 +- rust/crates/du-web/src/routes/tree.rs | 42 ++- .../du-web/templates/tree/snp_sidebar.html | 20 ++ rust/locales/en.txt | 4 + rust/locales/es.txt | 4 + rust/locales/fr.txt | 4 + rust/migrations/0053_biosample_str.sql | 29 ++ rust/migrations/0054_haplogroup_str_asr.sql | 24 ++ 16 files changed, 1049 insertions(+), 12 deletions(-) create mode 100644 documents/proposals/aging-pipeline-audit-mcdonald2021.md create mode 100644 rust/crates/du-db/src/str_import.rs create mode 100644 rust/crates/du-jobs/src/ftdna_str.rs create mode 100644 rust/migrations/0053_biosample_str.sql create mode 100644 rust/migrations/0054_haplogroup_str_asr.sql diff --git a/documents/proposals/aging-pipeline-audit-mcdonald2021.md b/documents/proposals/aging-pipeline-audit-mcdonald2021.md new file mode 100644 index 00000000..4eef0cb4 --- /dev/null +++ b/documents/proposals/aging-pipeline-audit-mcdonald2021.md @@ -0,0 +1,253 @@ +# Aging pipeline audit vs McDonald 2021 + +Audit of the DecodingUs Y-DNA branch-age pipeline against **McDonald, I. "Improved +Models of Coalescence Ages of Y-DNA Haplogroups." *Genes* 2021, 12, 862** +(doi:10.3390/genes12060862). The pipeline is intended to be a faithful +implementation of that paper; this document records where it is faithful, where it +diverges or approximates, what is missing, and candidate refinements (several of +which are under discussion by email with the author). + +Code under audit: +- `crates/du-db/src/pdf.rs` — discrete PDF machinery (Poisson/Gaussian/convolve/multiply/mixture/quantiles) +- `crates/du-db/src/age.rs` — SNP-Poisson propagation, genealogical anchors, Eq-1 combine, causality +- `crates/du-db/src/ystr.rs` — Y-STR `P(g|m)`, marker age PDFs, ancestral-motif reconstruction, STR propagation +- migrations `0013_str.sql`, `0014_str_age.sql`, `0051_y_callable_mask.sql` + +Paper reference frame: rates ~8×10⁻¹⁰ SNP·bp⁻¹·yr⁻¹; STR rates per-marker per-generation; +generation 33 yr; "present day" = tester/MRCA birth year. + +--- + +## 1. Verdict at a glance + +The **core probabilistic engine is faithful**: Eq 1 (product of evidence PDFs), +Eq 3 (Poisson SNP clock), Eq 11–14 + Table 1 (STR `P(g|m)` with multi-step ω), +ancestral-motif reconstruction (§2.5.2), and region-consistent SNP excision +(Appendix A.2/A.3) are all implemented and match the paper's numbers. + +The **largest gaps are in the uncertainty budget and the second-order priors**, +which the paper itself flags as dominant for real-world clades: + +1. **Mutation-rate uncertainty (σ_µ) is not propagated** (Appendix A.2.2/A.4). The + paper states this *dominates the error budget* for most real cases (~±67 yr + floor even at the well-tested R-S781). Our CIs are Poisson-only and therefore + too narrow for larger/older clades. **[highest priority]** +2. **STR saturation has no max-time taper** (Appendix A.5.2). The paper explicitly + recommends a Gaussian tapering of `P(g|m)` (or an exponential correction / revised + µ) to stop STRs underestimating deep ages via convergent mutations. We use a + crude tester-count gate instead. +3. **No NRR population-size prior** (Eq 25), **no shared-surname/autosomal historical + terms** (Eq 24, §2.6.3), **tester-birth offset omitted** (Appendix A.1), + **generational scatter `8√N` omitted** (Appendix A.5.1). + +There is also one **methodological divergence worth a decision**: our SNP TMRCA is a +*robust depth-quantile* estimator, **not** the paper's Eq 7/8 PDF-convolution build +(§3). The STR path *does* use Eq 8 convolution — so the two clocks are built by +different machinery (see §6). + +--- + +## 2. Section-by-section fidelity + +### §2.1 Driving principle — Eq 1 `P(t|e) = k ∏ P(t|eᵢ)` +**FAITHFUL.** `Pdf::multiply` (pointwise product + renormalise) is used per node to +combine SNP × STR × genealogical PDFs (`age.rs` combine loop). The `k` +normalisation is the renormalise step. Disjoint-product underflow falls back to an +inverse-variance Gaussian combine (`combine()`), which can't annihilate — a +reasonable engineering guard not in the paper. + +### §2.2 SNP PDFs from NGS — Eq 2–8 +- **Eq 3 `P(t|m)=Poisson(m,tbµ)`** — **FAITHFUL.** `Pdf::poisson_on` builds the + Gamma-shaped density `∝(tbµ)^m e^{-tbµ}` via log-sum-exp (overflow-safe for + deep backbone counts). Mode `m/bµ`, mean `(m+1)/bµ` verified by tests. +- **µ_SNP = 8.33×10⁻¹⁰** (`age.rs SNP_RATE`) — **FAITHFUL** to Helgason/[21] + (Icelandic, 33.5 yr/gen, MSY 21.3 Mbp). Single combined rate, matching the + paper's choice (§2.2.1/Appendix A.4); X-deg/ampliconic vs palindromic rate split + is not modelled (paper also uses the combined rate by default). ✓ +- **Eq 4 `b̄` (effective callable bp)** — **PARTIAL.** Paper defines `b̄` as the + intersection of coverage callable in **≥2 sub-clades** (≈ second-highest + coverage; §3.4.2 "second-highest coverage among constituent sub-clades"). Our + preferred path is a **uniform joint-call mask** (`genomics.y_callable_interval`, + mig 0051) used region-consistently for numerator *and* denominator — a cleaner + equivalent. **But** the fallback (no mask) uses the **per-node mean** tester + callable bp, not the paper's second-highest. Divergent only when the mask is + absent; flag for the fallback path. +- **§2.2.2 µ uncertainty (σ_µ)** — **MISSING.** `poisson_on` takes a point `mu`. + The paper applies σ_µ as a multiplicative tree-wide scaling (or per-node + broadening) at the time-conversion step. `HALLAST_RATE{,_LO,_HI}` constants exist + but are not folded into the CI. **This is the dominant missing error term.** +- **Eq 5–8 (clade build by convolution)** — **DIVERGENT (deliberate).** See §6. + We replace the Eq-6 YFull child-averaging (which the paper rejects) with a + *robust depth-quantile* estimator (q90 of descendant root-to-tip SNP depths, with + a q97≥150-SNP deep-ladder gate and a 2-children corroboration floor), then emit a + single Poisson PDF from that count and convolve the branch time for `formed` + (`propagate()`). This fixes deep-clade over-aging but is not the paper's Eq 7/8 + PDF product. + +### §2.3 Causality — Eq 9/10 +- **IMPLEMENTED (simplified).** Two layers now enforce parent ≥ child: + (a) the SNP path clamps child counts to parents top-down (`propagate`); (b) the + COMBINED path projects medians bottom-up, raising each parent to its oldest child + with a CI shift (`recompute_combined_ages`, added 2026-06). This *guarantees* the + constraint (verified: 0 inversions tree-wide). +- **DIVERGENT from the exact method.** The paper's Eq 9/10 *derive the child from + the parent* by reverse-convolution `P(t_c)=P(t_p) ⊛ P(t_{c→p}|−m)` and take a + **√(kits)-weighted average of Eq 9 and Eq 10** with `P(t<0)=0`. We instead raise + the parent to the child (a median projection), because our inversions came from + *under-aged parents* (a young STR term), not over-aged children. Note + `Pdf::convolve_sub` (the Eq 9/10 reverse convolution) already exists but is + currently unused — the building block for a faithful Eq 9/10 is in place. + +### §2.4 Ancient DNA +- **PARTIAL.** aDNA can be entered as a `tree.genealogical_anchor` with + `carbon_date_bp`. The paper models aDNA as a **lower limit** — the *cumulative* + calibrated ¹⁴C distribution (novel variants often uncallable) — and notes σ_µ + must be applied before mixing aDNA calendar dates with SNP ages. We model anchors + as symmetric Gaussians, not cumulative/one-sided, and don't apply σ_µ. Adequate + for proven-genealogy anchors; not yet a faithful aDNA constraint. + +### §2.5 Y-STR PDFs — Eq 11–23, Table 1 +- **Eq 11–12 `m_s=tµ_s`, `P(t|m_s)=Poisson⊗P(µ_s)`** — **PARTIAL.** `marker_age_pdf` + builds `Poisson(m, t·µ_year)` with `µ_year = mutation_rate/33`. The **⊗P(µ_s) + convolution over STR rate uncertainty is omitted** (we have + `mutation_rate_lower/upper` columns but use the point rate) — same σ_µ gap as SNPs. +- **Eq 13–14, Table 1 `P(g|m)`** — **FAITHFUL.** McDonald's Table 1 is embedded + verbatim (`ystr.rs TABLE1`, g,m≤10) and extended beyond its range by an + all-orders signed-step ω convolution; `marker_age_pdf` forms the mixture + `P(t|g)=Σ_m P(g|m)P(t|m)` (`Pdf::mixture`). Row/column sums ≈1 verified. +- **§2.5.3 multi-step ω** — **FAITHFUL.** `ω±1=0.96217, ω±2=0.032, ω±3=0.004`, ÷√10 + per further repeat, `w₊=w₋=0.5` default; per-marker `omega_plus/minus` drive + asymmetric tables when present. Exact match to Eq 15. +- **§2.5.1 "twice the TMRCA" for two modern tests** — **HANDLED CORRECTLY.** We + measure tester→reconstructed-ancestral-motif distance (1× TMRCA), not pairwise + (2×), so no double-count. +- **§2.5.2 ancestral-motif reconstruction** — **FAITHFUL.** `reconstruct_motifs` + does the up-pass (modal of sub-clades + direct testers) then down-pass (parent + fill), simple markers only. The paper warns this **underestimates** when mutations + are missed (null/ambiguous ancestral alleles bias young) — which is exactly the + collapse we observed and mitigated with the ≥2-tester combine gate + (`MIN_STR_TESTERS_FOR_COMBINE`). See §3 / §5. +- **Eq 23 / §2.5.4 multi-copy markers** — **DIVERGENT.** Paper recommends including + multi-copy markers (DYS464/CDY/YCAII/DYS385…) as **binary** "any mutation?" + (`P(g=0)` vs `Σ P(g>0)`). We **exclude multi-copy entirely** from scoring and + propagation (simple markers only). This drops signal from the fastest markers, + most useful for recent clades. Candidate refinement. +- **Eq 22 mixed multi-step (`f₂`+)** — acceptably neglected (paper says terms + `f₃`+ and mixed `[+²,−²]` are <1/900th weight). + +### §2.6 Historical / ancillary +- **§2.6.1 paper genealogies** — **PARTIAL.** Point/Gaussian anchors supported + (`genealogical_anchor` + 10%-or-25yr σ default). δ-function, boxcar, log-normal + shapes from the paper are not distinctly modelled (all collapse to Gaussian). +- **§2.6.2 shared surnames (Eq 24)** — **MISSING** (`ψ₂ + ½(1−ψ₂)[1+erf(...)]`). +- **§2.6.3 autosomal DNA constraints** — **MISSING** (log-normal / cumulative + Gaussian on close relatedness). +- **§2.6.4 NRR population-size prior (Eq 25 `P(t|NRR)=NRR^{t/G}`)** — **MISSING.** + The paper notes this skews TMRCAs older (many more distant than close cousins; + ~30% for NRR=1.3) and can dominate the budget for very-well-tested recent + families. Notable omission for the genealogical (recent) timescale. + +### §2.7 / Appendix A — calibration constants +| Item | Paper | Code | Verdict | +|---|---|---|---| +| Present day (A.1) | tester/MRCA birth; add ~64 yr (95% CI 35–91) to TMRCA | `PRESENT_YEAR=1950` for ¹⁴C; **tester-birth offset omitted** ("≈present, negligible") | **MINOR GAP** — ages ~1 generation too young; lose that CI term | +| SNP rate (A.4) | ~8×10⁻¹⁰, σ≈±8–10% | `SNP_RATE=8.33e-10`; σ_µ not used | rate ✓, **σ_µ MISSING** | +| Recurrent/errant SNPs (A.2/A.3) | excise PAR1/PAR2/DYZ19/Yq12/centromere (~940 kb), retain palindromes; mask SNP clusters <~100 bp (MNPs) | `recurrent_region_mask_sql` masks heterochromatin + inverted_repeat; palindromes retained; joint-call mask excludes non-callable | **MOSTLY FAITHFUL**; palindrome-retention ✓; **MNP <100 bp clustering not removed** | +| Generation length (A.5.1) | 33 yr (95% CI 29–37) + per-gen scatter `8√N` yr | `GENERATION_YEARS=33` | mean ✓; **scatter & CI MISSING** | +| STR rates (A.5.2) | per-marker, homogeneous study (Willems 2016 [47]) | `genomics.str_mutation_rate` per-marker, `DEFAULT_STR_RATE=0.0025` fallback; 137 rows seeded | **FAITHFUL** (ensure full ~700-marker coverage) | +| STR max-time taper (A.5.2) | Gaussian taper on `P(g|m)` / exp correction / revised µ at large t | none (only ≥2-tester gate) | **MISSING** | +| Pop-growth bias (A.4.1) | larger clades gain ~2.3 extra SNPs; remedy = causality | causality projection | acceptable | + +--- + +## 3. What the worked examples confirm + +- **Example 3 (Fig 3, ~4 ky):** STR-only underestimates badly (2560 vs 4000 yr) + from convergent/hidden mutations; combined SNP+STR (3667) beats SNP-only (3873). + → validates both that STR *adds* precision **and** that it needs a deep-time + taper (our §1.2 gap). Our reconstruction-collapse is the same phenomenon. +- **Example 4 (R-S781):** σ_µ sets a ~±67 yr error floor even with good data; + YFull's narrower CI is called out as *not* carrying full Poisson + rate + uncertainty. → our Poisson-only CIs share YFull's under-coverage; §1.1 gap. +- **`b̄` = second-highest sub-clade coverage; palindromes retained; ~940 kb of + PAR/DYZ19/Yq12/centromere excised** — confirms our mask design intent. + +--- + +## 4. Faithful — keep as is +- Eq 1 product combine (`multiply`) ✓ +- Eq 3 Poisson clock, log-safe (`poisson_on`) ✓ +- Table 1 `P(g|m)` verbatim + ω convolution extension ✓ +- Multi-step ω (0.96217/0.032/0.004, ÷√10) ✓ +- Per-marker STR rates from `str_mutation_rate` ✓ +- Ancestral-motif up/down reconstruction (§2.5.2) ✓ +- Region-consistent SNP excision; palindromes retained (A.2/A.3) ✓ +- 33 yr/gen, 1950 ¹⁴C present ✓ +- Causality guarantee (parent ≥ child) ✓ (mechanism simplified — see §6) + +--- + +## 5. Recently fixed (2026-06, this work) +- **Causality on COMBINED** (§2.3): bottom-up parent≥child projection — was wholly + absent from the combined term; caused the R-A804/R-A802 inversion. +- **STR collapse gate** (§2.5.2): `MIN_STR_TESTERS_FOR_COMBINE=2` keeps + reconstruction-only nodes (≈0 genetic distance ⇒ spuriously young+narrow PDF) + from dominating the SNP clock. A stop-gap for the missing §A.5.2 taper. +- Undatable nodes no longer retain stale denormalised `tmrca_ybp`. + +--- + +## 6. The central methodological question: SNP propagation + +The paper builds a node's age as the **normalised product of its children's PDFs, +each convolved with its branch-time PDF** (Eq 7/8), with `P(t_p<0)=0` guaranteeing +causality *by construction*. Our **STR** path does exactly this (`str_tmrca`: +`child ⊛ branch`, then product). Our **SNP** path does **not** — it computes a +robust **q90 depth-quantile** of descendant SNP depths and emits one Poisson PDF +(`propagate`). This was a deliberate fix for deep-clade over-aging (a single +over-called deep tip was inflating ancestors under naive max/averaging), but it +means: + +- the two clocks are built by **different estimators** (inconsistent uncertainty + semantics), and +- the SNP side departs from the paper's probabilistic Eq 8. + +**Decision needed:** either (a) return the SNP side to the Eq 7/8 convolution build +and layer the robustness fixes (outlier-robust child weighting, deep-ladder gate) +on top of it — unifying with the STR path and recovering true PDF shape; or +(b) keep the quantile estimator and document it as an intentional, validated +divergence. This is a good email topic with the author. + +--- + +## 7. Prioritised refinements + +**P1 — Propagate mutation-rate uncertainty σ_µ (A.2.2/A.4, Eq 23).** +The paper's dominant error term. Convolve/scale the final age PDFs by the rate +uncertainty (SNP ±~8–10% via `HALLAST_RATE` CI already present; STR via +`mutation_rate_lower/upper`). Widens CIs to honest coverage; biggest accuracy win. + +**P2 — STR deep-time taper (A.5.2).** +Add a Gaussian taper on `P(g|m)` (mathematically preferred) or a max-time/exp +correction so STR stops underestimating old clades via convergent mutations. +Replaces the blunt ≥2-tester gate with a principled saturation model; lets STR +inform internal nodes again. + +**P3 — Faithful Eq 9/10 causality** (reverse-convolution `convolve_sub`, √kits +weighted Eq9/Eq10 average) replacing the median projection — redistributes PDF mass +instead of shifting medians. + +**P4 — Tester-birth offset + generational scatter** (A.1, A.5.1): add ~64 yr (95% +CI 35–91) at the tip and `8√N`-yr generational scatter to STR ages. Small, cheap, +improves both central value and CI. + +**P5 — NRR population-size prior (Eq 25)** for the recent/genealogical timescale; +**multi-copy STR binary inclusion (§2.5.4)** for fast-marker signal; **shared-surname +(Eq 24) / autosomal** historical terms; **MNP <100 bp cluster removal (A.3)**. + +**P6 — Unify SNP propagation with Eq 8** (see §6) — larger, do after P1–P2. + +--- + +*Generated 2026-06-29. Pipeline state at audit: 10,614 COMBINED estimates, 2,013 +STR_VARIANCE, 0 parent Result, Option)> = sqlx::query_as( - "SELECT haplogroup_id, estimate_ybp, ci_low_ybp, ci_high_ybp \ + // Reconstruct a Gaussian for any STR_VARIANCE row with no fresh PDF (a curated + // value, or one predating profile data). Gated on the same within-clade + // diversity floor as the propagated PDFs (sample_count ≥ MIN_STR_TESTERS_FOR_COMBINE) + // so the fallback can't re-admit the reconstruction-collapsed nodes the + // propagation path excluded — see ystr::MIN_STR_TESTERS_FOR_COMBINE. + let str_rows: Vec<(i64, i32, Option, Option, Option)> = sqlx::query_as( + "SELECT haplogroup_id, estimate_ybp, ci_low_ybp, ci_high_ybp, sample_count \ FROM tree.haplogroup_age_estimate WHERE method='STR_VARIANCE' AND estimate_ybp IS NOT NULL", ) .fetch_all(pool) .await?; - for (hg, est, lo, hi) in str_rows { - if str_pdf.contains_key(&hg) { + for (hg, est, lo, hi, n) in str_rows { + if str_pdf.contains_key(&hg) + || (n.unwrap_or(0) as usize) < crate::ystr::MIN_STR_TESTERS_FOR_COMBINE + { continue; } let mean = est as f64; @@ -657,6 +664,45 @@ pub async fn recompute_combined_ages(pool: &PgPool) -> Result = HashMap::with_capacity(combined_writes.len()); + for (idx, w) in combined_writes.iter().enumerate() { + by_id.insert(w.0, idx); + } + let mut corrected = 0usize; + for &i in post_order(&clades).iter() { + let Some(&pw) = by_id.get(&ids[i]) else { continue }; + let parent_med = combined_writes[pw].1; + let mut target = parent_med; + for &c in &clades[i].children { + if let Some(&cw) = by_id.get(&ids[c]) { + target = target.max(combined_writes[cw].1); + } + } + if target > parent_med { + let delta = target - parent_med; + combined_writes[pw].1 = target; // median + combined_writes[pw].2 += delta; // ci low + combined_writes[pw].3 += delta; // ci high + corrected += 1; + } + } + if corrected > 0 { + tracing::info!(corrected, "combined-age causality back-correction (parent ≥ child)"); + } + } + // ── PHASE 2: write everything in one short transaction. The writes are tiny and // back-to-back, so the connection never idles long enough to be reaped (no CPU // work happens between them). @@ -690,6 +736,22 @@ pub async fn recompute_combined_ages(pool: &PgPool) -> Result &'static str { + match self { + SkipReason::Empty => "empty", + SkipReason::Html => "html", + SkipReason::LongFormat => "long-format", + SkipReason::AllMissing => "all-missing", + } + } +} + +/// Split a CSV line into trimmed, unquoted fields. The vendor exports never +/// embed a comma inside a value (marker names use `-`, not `,`), so a plain +/// comma split — then strip surrounding quotes + whitespace — is sufficient and +/// handles the quoted, unquoted, and space-padded variants alike. +fn fields(line: &str) -> Vec { + line.split(',') + .map(|f| f.trim().trim_matches('"').trim().to_string()) + .collect() +} + +/// Parse one vendor token into a lexicon `strValue`, or `None` if missing. +pub fn parse_token(raw: &str) -> Option { + let t = raw.trim().trim_matches('"').trim(); + if t.is_empty() || t == "-" { + return None; + } + let parts: Vec<&str> = t.split('-').map(str::trim).filter(|p| !p.is_empty()).collect(); + let ints: Option> = parts.iter().map(|p| p.parse::().ok()).collect(); + Some(match ints { + Some(v) if v.len() == 1 => json!({ "type": "simple", "repeats": v[0] }), + Some(v) if v.len() > 1 => json!({ "type": "multiCopy", "copies": v }), + // Unparseable (partial repeats, footnote marks): keep losslessly, unscored. + _ => json!({ "type": "complex", "raw": t }), + }) +} + +/// Parse a wide-format vendor CSV (header row of markers + one value row) into a +/// lexicon profile. `Err(SkipReason)` for inputs that aren't a wide export. +pub fn parse_wide_csv(text: &str) -> Result { + let lower = text.to_ascii_lowercase(); + if lower.contains(" = text.lines().filter(|l| !l.trim().is_empty()).collect(); + if lines.len() < 2 { + return Err(SkipReason::Empty); + } + // The `marker;value;` long layout has one marker per line and no comma header. + if lines[0].contains(';') && !lines[0].contains(',') { + return Err(SkipReason::LongFormat); + } + let header = fields(lines[0]); + let values = fields(lines[1]); + let mut markers = Vec::new(); + for (name, raw) in header.iter().zip(values.iter()) { + if name.is_empty() { + continue; + } + if let Some(value) = parse_token(raw) { + markers.push(json!({ "marker": name, "value": value })); + } + } + if markers.is_empty() { + return Err(SkipReason::AllMissing); + } + let total_markers = markers.len() as i32; + Ok(ParsedProfile { markers: Value::Array(markers), total_markers }) +} + +/// Map the given accessions to their `core.biosample.sample_guid` (absent +/// accessions are simply missing from the result — e.g. a NAS sample not loaded +/// into this database). The cohort `subject_id` is stored as the accession. +pub async fn guids_by_accessions( + pool: &PgPool, + accessions: &[String], +) -> Result, DbError> { + let rows: Vec<(String, Uuid)> = sqlx::query_as( + "SELECT accession, sample_guid FROM core.biosample \ + WHERE accession = ANY($1) AND deleted = false", + ) + .bind(accessions) + .fetch_all(pool) + .await?; + Ok(rows.into_iter().collect()) +} + +/// Upsert a sample's imported STR profile. A sample with several vendor kits +/// keeps the most complete one (highest marker count); re-imports of the same +/// kit are idempotent. +pub async fn upsert_profile( + pool: &PgPool, + sample_guid: Uuid, + imported_from: &str, + panel: &str, + profile: &ParsedProfile, + source_file: &str, +) -> Result<(), DbError> { + sqlx::query( + "INSERT INTO genomics.biosample_str_profile \ + (sample_guid, source, imported_from, panel, total_markers, markers, source_file) \ + VALUES ($1, 'IMPORTED', $2, $3, $4, $5, $6) \ + ON CONFLICT (sample_guid) DO UPDATE SET \ + imported_from = EXCLUDED.imported_from, panel = EXCLUDED.panel, \ + total_markers = EXCLUDED.total_markers, markers = EXCLUDED.markers, \ + source_file = EXCLUDED.source_file, imported_at = now() \ + WHERE EXCLUDED.total_markers >= genomics.biosample_str_profile.total_markers", + ) + .bind(sample_guid) + .bind(imported_from) + .bind(panel) + .bind(profile.total_markers) + .bind(&profile.markers) + .bind(source_file) + .execute(pool) + .await?; + Ok(()) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn parses_tokens() { + assert_eq!(parse_token(" 13"), Some(json!({"type":"simple","repeats":13}))); + assert_eq!(parse_token("\"13\""), Some(json!({"type":"simple","repeats":13}))); + assert_eq!(parse_token(" 11-15"), Some(json!({"type":"multiCopy","copies":[11,15]}))); + assert_eq!( + parse_token(" 15-15-16-17"), + Some(json!({"type":"multiCopy","copies":[15,15,16,17]})) + ); + assert_eq!(parse_token(" -"), None); + assert_eq!(parse_token(" "), None); + // partial-repeat / footnoted values are preserved as complex, not dropped + assert_eq!(parse_token("10.2"), Some(json!({"type":"complex","raw":"10.2"}))); + } + + #[test] + fn parses_wide_csv_quoted() { + let csv = "DYS393,DYS390,DYS385,DYS459\n\" 13\",\" 24\",\" 11-15\",\" 9-10\"\n"; + let p = parse_wide_csv(csv).unwrap(); + assert_eq!(p.total_markers, 4); + let arr = p.markers.as_array().unwrap(); + assert_eq!(arr[0], json!({"marker":"DYS393","value":{"type":"simple","repeats":13}})); + assert_eq!(arr[2], json!({"marker":"DYS385","value":{"type":"multiCopy","copies":[11,15]}})); + } + + #[test] + fn parses_wide_csv_unquoted_and_skips_missing() { + let csv = "DYS393,DYS390,DYS19\n13,24,-\n"; + let p = parse_wide_csv(csv).unwrap(); + assert_eq!(p.total_markers, 2); // DYS19 = '-' omitted + } + + #[test] + fn classifies_non_wide_inputs() { + assert_eq!(parse_wide_csv(" \n"), Err(SkipReason::Empty)); + assert_eq!(parse_wide_csv("\n"), Err(SkipReason::Html)); + assert_eq!(parse_wide_csv("DYS393;13;\nDYS390;24;\n"), Err(SkipReason::LongFormat)); + } +} diff --git a/rust/crates/du-db/src/ystr.rs b/rust/crates/du-db/src/ystr.rs index e581229d..f0a38305 100644 --- a/rust/crates/du-db/src/ystr.rs +++ b/rust/crates/du-db/src/ystr.rs @@ -1,6 +1,7 @@ -//! Y-STR per-branch modal signatures: aggregate the mirrored STR profiles -//! (`fed.str_profile`) by SNP-defined branch (haplogroup) into a modal haplotype -//! stored in `tree.haplogroup_ancestral_str`. This is the catalog-side half of +//! Y-STR per-branch modal signatures: aggregate the STR profiles — federated +//! (`fed.str_profile`) plus locally-imported (`genomics.biosample_str_profile`, +//! see [`crate::str_import`]) — by SNP-defined branch (haplogroup) into a modal +//! haplotype stored in `tree.haplogroup_ancestral_str`. This is the catalog-side half of //! the STR feature (the mirror is `fed::str_profile`); prediction over these //! signatures is Phase 2. //! @@ -640,13 +641,25 @@ struct StrInputs { /// motifs. Shared by [`recompute_signatures`] (which writes the signatures + ages) /// and [`str_tmrca_pdfs`] (which the combined-age step consumes). async fn build_str_inputs(pool: &PgPool) -> Result { + // Two sources, one population per branch: federated profiles (Navigator's + // mirrored com.decodingus.atmosphere.strProfile, keyed by the published + // y_haplogroup name) UNION locally-imported vendor profiles (academic / D2C + // cohorts in core.biosample, joined live through their Y tree placement — + // tree.haplogroup_sample — so the branch is the current one, not a name copy). let rows: Vec<(i64, Value)> = sqlx::query_as( "SELECT h.id, sp.markers \ FROM fed.str_profile sp \ JOIN fed.biosample b ON b.at_uri = sp.biosample_ref \ JOIN tree.haplogroup h ON h.name = b.y_haplogroup \ AND h.haplogroup_type::text = 'Y_DNA' AND h.valid_until IS NULL \ - WHERE b.y_haplogroup IS NOT NULL", + WHERE b.y_haplogroup IS NOT NULL \ + UNION ALL \ + SELECT hs.haplogroup_id, lp.markers \ + FROM genomics.biosample_str_profile lp \ + JOIN tree.haplogroup_sample hs ON hs.sample_guid = lp.sample_guid \ + AND hs.dna_type = 'Y_DNA'::core.dna_type AND hs.status = 'PLACED' \ + JOIN tree.haplogroup h ON h.id = hs.haplogroup_id AND h.valid_until IS NULL \ + WHERE hs.haplogroup_id IS NOT NULL", ) .fetch_all(pool) .await?; @@ -711,12 +724,31 @@ async fn build_str_inputs(pool: &PgPool) -> Result { Ok(StrInputs { ids, children, motif, testers, models, modal_by_hg }) } +/// Minimum direct STR testers on a node for its propagated TMRCA to enter the +/// combined-age product. STR ages are only trustworthy where there is genuine +/// within-clade genetic *diversity*: with 0 testers a node's age is pure motif +/// reconstruction, and with 1 the reconstructed ancestral motif equals that lone +/// tester — both make the genetic distance ≈0, so the STR PDF collapses to a +/// spuriously young, spuriously narrow estimate (McDonald §2.5.2: "TMRCAs become +/// underestimated if mutations are missed") that then dominates the sharp SNP +/// term. Requiring ≥2 testers keeps STR where the paper says it is most useful +/// (well-populated clades) and lets SNP date the rest; the COMBINED causality +/// pass ([`crate::age`]) is the hard backstop. +pub const MIN_STR_TESTERS_FOR_COMBINE: usize = 2; + /// Per-node tree-propagated STR TMRCA PDFs (keyed by haplogroup id) on the given -/// grid — consumed by the combined-age PDF product in [`crate::age`]. +/// grid — consumed by the combined-age PDF product in [`crate::age`]. Restricted +/// to nodes with real within-clade STR diversity (≥[`MIN_STR_TESTERS_FOR_COMBINE`] +/// direct testers); reconstruction-only nodes are left to the SNP clock. The +/// per-branch `STR_VARIANCE` rows (written by [`recompute_signatures`]) are not +/// gated — they remain a displayed contributing factor. pub async fn str_tmrca_pdfs(pool: &PgPool, res: f64, max_age: f64) -> Result, DbError> { let inp = build_str_inputs(pool).await?; let ages = propagate_str(&inp.children, &inp.motif, &inp.testers, &inp.models, res, max_age); - Ok(inp.ids.iter().zip(ages).filter_map(|(&id, a)| a.map(|p| (id, p))).collect()) + Ok((0..inp.ids.len()) + .filter(|&i| inp.testers[i].len() >= MIN_STR_TESTERS_FOR_COMBINE) + .filter_map(|i| ages[i].clone().map(|p| (inp.ids[i], p))) + .collect()) } // ── queries ─────────────────────────────────────────────────────────────────── @@ -766,6 +798,8 @@ pub struct RecomputeStats { pub haplogroups: usize, pub markers: usize, pub age_estimates: usize, + /// Reconstructed-ancestral-motif marker rows written to tree.haplogroup_str_asr. + pub asr_markers: usize, } /// Recompute every Y branch's modal signature from the mirrored profiles, then the @@ -855,6 +889,38 @@ pub async fn recompute_signatures(pool: &PgPool) -> Result Result, +} + +impl StrAsrMarker { + /// True when this marker mutated on the branch leading into the node. + pub fn changed(&self) -> bool { + matches!(self.parent_value, Some(p) if p != self.ancestral_value) + } +} + +/// A Y haplogroup's reconstructed ancestral STR motif (by name), each marker +/// paired with the parent's value. Markers sorted; empty when no ASR is stored +/// (no STR evidence in the subtree). The parent is the node's current valid +/// Y-tree parent — `tree.haplogroup_str_asr` is keyed to live nodes, so the diff +/// always reflects the present topology. +pub async fn branch_str_asr(pool: &PgPool, haplogroup: &str) -> Result, DbError> { + let rows: Vec<(String, i32, Option)> = sqlx::query_as( + "SELECT a.marker_name, a.ancestral_value, p.ancestral_value AS parent_value \ + FROM tree.haplogroup h \ + JOIN tree.haplogroup_str_asr a ON a.haplogroup_id = h.id \ + LEFT JOIN tree.haplogroup_relationship r \ + ON r.child_haplogroup_id = h.id AND r.valid_until IS NULL \ + LEFT JOIN tree.haplogroup_str_asr p \ + ON p.haplogroup_id = r.parent_haplogroup_id AND p.marker_name = a.marker_name \ + WHERE h.name = $1 AND h.haplogroup_type::text = 'Y_DNA' AND h.valid_until IS NULL \ + ORDER BY a.marker_name", + ) + .bind(haplogroup) + .fetch_all(pool) + .await?; + Ok(rows + .into_iter() + .map(|(marker_name, ancestral_value, parent_value)| StrAsrMarker { + marker_name, + ancestral_value, + parent_value, + }) + .collect()) +} + /// A per-branch age estimate (one contributing factor; method-labeled). #[derive(Debug, sqlx::FromRow)] pub struct AgeEstimate { @@ -990,6 +1107,18 @@ mod tests { (name.into(), StrValue::Simple(n)) } + #[test] + fn str_asr_marker_changed() { + let mk = |p: Option, v: i32| StrAsrMarker { + marker_name: "DYS393".into(), + ancestral_value: v, + parent_value: p, + }; + assert!(mk(Some(12), 13).changed()); // mutated on this branch + assert!(!mk(Some(13), 13).changed()); // unchanged from parent + assert!(!mk(None, 13).changed()); // root / no parent reconstruction + } + #[test] fn parses_marker_union() { let markers = json!([ diff --git a/rust/crates/du-jobs/Cargo.toml b/rust/crates/du-jobs/Cargo.toml index 48656256..6645ab53 100644 --- a/rust/crates/du-jobs/Cargo.toml +++ b/rust/crates/du-jobs/Cargo.toml @@ -24,6 +24,7 @@ anyhow = { workspace = true } serde = { workspace = true } serde_json = { workspace = true } chrono = { workspace = true } +uuid = { workspace = true } tokio-tungstenite = { workspace = true } futures-util = { workspace = true } diff --git a/rust/crates/du-jobs/src/ftdna_str.rs b/rust/crates/du-jobs/src/ftdna_str.rs new file mode 100644 index 00000000..69d068ef --- /dev/null +++ b/rust/crates/du-jobs/src/ftdna_str.rs @@ -0,0 +1,245 @@ +//! FTDNA Y-STR import. Bulk-loads the vendor "DYS Results" exports on the NAS +//! into `genomics.biosample_str_profile`, matched to the academic / D2C cohort. +//! +//! Layout: `FTDNA_STR_DIR//DYS_Results.csv` (or `strs.csv`) — one kit per +//! subdirectory. The directory name is the FTDNA kit number; a few are already a +//! cohort `subject_id` (UUID). Resolution is two hops: +//! kit ─(cohort_manifest.tsv: `kit` column, comma-split)→ subject_id +//! subject_id ─(= core.biosample.accession)→ sample_guid +//! The manifest's `kit` column is a comma-separated list (multi-vendor donors), +//! so every token is indexed and matched, not just single-lab rows. +//! +//! Idempotent + safe to re-run: parsing is pure ([`du_db::str_import`]) and the +//! upsert keeps each sample's most complete profile. Defaults to a preview that +//! rolls nothing back because it simply doesn't write — pass apply to commit. +//! +//! FTDNA_STR_DIR=/Volumes/nas/Genomics/FTDNA_STR \ +//! COHORT_MANIFEST=/Volumes/nas/Genomics/d2c/cohort_manifest.tsv \ +//! decodingus-jobs run-once ftdna-str --apply + +use std::collections::HashMap; +use std::fs; +use std::path::{Path, PathBuf}; + +use anyhow::{Context, Result}; +use du_db::str_import; +use du_db::PgPool; + +const IMPORTED_FROM: &str = "FTDNA"; +const PANEL: &str = "FTDNA Y-DNA DYS"; + +pub struct Config { + pub str_dir: PathBuf, + pub manifest: PathBuf, + pub apply: bool, + pub limit: usize, + /// Direct kit-token → biosample-accession overrides, for samples whose + /// biosample accession isn't the cohort `subject_id` (own/de-novo genomes + /// loaded under a friendly accession — e.g. FTDNA kit `B5163` is `WGS229`). + /// Bypasses the manifest kit→subject_id→accession chain. + pub aliases: HashMap, +} + +impl Config { + /// Build from env (`FTDNA_STR_DIR`, `COHORT_MANIFEST`, optional + /// `FTDNA_STR_ALIAS="KIT=ACCESSION,KIT=ACCESSION"`) + the trailing CLI flags + /// (`--apply`, `--limit N`). Returns `None` if a required path is unset. + pub fn from_env(args: &[String]) -> Option { + let str_dir = env_path("FTDNA_STR_DIR")?; + let manifest = env_path("COHORT_MANIFEST")?; + let apply = args.iter().any(|a| a == "--apply"); + let limit = args + .iter() + .position(|a| a == "--limit") + .and_then(|i| args.get(i + 1)) + .and_then(|n| n.parse().ok()) + .unwrap_or(0); + let aliases = std::env::var("FTDNA_STR_ALIAS") + .ok() + .map(|s| { + s.split(',') + .filter_map(|pair| pair.split_once('=')) + .map(|(k, v)| (k.trim().to_string(), v.trim().to_string())) + .filter(|(k, v)| !k.is_empty() && !v.is_empty()) + .collect() + }) + .unwrap_or_default(); + Some(Config { str_dir, manifest, apply, limit, aliases }) + } +} + +fn env_path(key: &str) -> Option { + std::env::var(key).ok().filter(|s| !s.is_empty()).map(PathBuf::from) +} + +/// kit-token → subject_id from the cohort manifest. The `kit` column is a +/// comma-separated list aligned (loosely) to the `lab` list; every token is +/// indexed so an FTDNA kit resolves regardless of the donor's other vendors. +fn load_kit_map(manifest: &Path) -> Result> { + let text = fs::read_to_string(manifest) + .with_context(|| format!("read manifest {}", manifest.display()))?; + let mut lines = text.lines(); + let header = lines.next().context("empty manifest")?; + let cols: Vec<&str> = header.split('\t').collect(); + let subj_i = cols.iter().position(|c| *c == "subject_id").context("no subject_id column")?; + let kit_i = cols.iter().position(|c| *c == "kit").context("no kit column")?; + + let mut map = HashMap::new(); + for line in lines { + let f: Vec<&str> = line.split('\t').collect(); + let (Some(subj), Some(kits)) = (f.get(subj_i), f.get(kit_i)) else { continue }; + let subj = subj.trim(); + if subj.is_empty() { + continue; + } + for kit in kits.split(',') { + let kit = kit.trim(); + if !kit.is_empty() { + map.insert(kit.to_string(), subj.to_string()); + } + } + } + Ok(map) +} + +/// Normalize a kit directory name to its kit token: strip the loader-suffixes +/// some dirs carry (`MK74582_str`, `N256757_YDNA_DYS`) and recover the kit from +/// the one dir saved under a full download path (`-Users-…-169776_YDNA_DYS…`). +fn kit_token(dir_name: &str) -> String { + if let Some(rest) = dir_name.strip_prefix("-Users") { + if let Some(pos) = rest.find("_YDNA_DYS") { + if let Some(kit) = rest[..pos].rsplit('-').next() { + return kit.to_string(); + } + } + } + for suf in ["_str", "_YDNA_DYS"] { + if let Some(base) = dir_name.strip_suffix(suf) { + return base.to_string(); + } + } + dir_name.to_string() +} + +fn looks_like_uuid(s: &str) -> bool { + s.len() == 36 && s.as_bytes().iter().filter(|&&b| b == b'-').count() == 4 +} + +/// The first `*.csv` in a kit directory (the export name varies: +/// `DYS_Results.csv` / `strs.csv`). +fn find_csv(dir: &Path) -> Option { + fs::read_dir(dir) + .ok()? + .flatten() + .map(|e| e.path()) + .find(|p| p.extension().and_then(|e| e.to_str()) == Some("csv")) +} + +#[derive(Default)] +struct Report { + dirs: usize, + no_csv: usize, + unresolved_kit: usize, // dir didn't map to a subject_id + not_in_db: usize, // subject_id not a loaded biosample + skipped: HashMap<&'static str, usize>, // by SkipReason + loaded: usize, +} + +pub async fn run(pool: &PgPool, cfg: &Config) -> Result<()> { + let kit_map = load_kit_map(&cfg.manifest)?; + + // Discover kit dirs and resolve each to a subject_id (kit token, or a dir + // that is itself a subject_id UUID). + let mut dirs: Vec<(String, PathBuf, String)> = Vec::new(); // (subject_id, csv, dir_name) + let mut rep = Report::default(); + let mut entries: Vec = fs::read_dir(&cfg.str_dir) + .with_context(|| format!("read {}", cfg.str_dir.display()))? + .flatten() + .map(|e| e.path()) + .filter(|p| p.is_dir()) + .collect(); + entries.sort(); + for dir in entries { + rep.dirs += 1; + let name = dir.file_name().and_then(|s| s.to_str()).unwrap_or("").to_string(); + let token = kit_token(&name); + // Alias overrides win over the manifest (own/de-novo genomes keyed by a + // friendly accession rather than the cohort subject_id). + let subject_id = match cfg.aliases.get(&token).or_else(|| kit_map.get(&token)) { + Some(s) => s.clone(), + None if looks_like_uuid(&name) => name.clone(), + None => { + rep.unresolved_kit += 1; + continue; + } + }; + let Some(csv) = find_csv(&dir) else { + rep.no_csv += 1; + continue; + }; + dirs.push((subject_id, csv, name)); + } + + // Resolve subject_ids to biosample guids in one query. + let accs: Vec = dirs.iter().map(|(s, _, _)| s.clone()).collect(); + let guids = str_import::guids_by_accessions(pool, &accs).await?; + + let mut sample = Vec::new(); + for (subject_id, csv, name) in &dirs { + if cfg.limit > 0 && rep.loaded >= cfg.limit { + break; + } + let Some(&guid) = guids.get(subject_id) else { + rep.not_in_db += 1; + continue; + }; + let text = match fs::read_to_string(csv) { + Ok(t) => t, + Err(e) => { + tracing::warn!(dir = %name, error = %e, "read failed"); + *rep.skipped.entry("read-error").or_default() += 1; + continue; + } + }; + let profile = match str_import::parse_wide_csv(&text) { + Ok(p) => p, + Err(reason) => { + *rep.skipped.entry(reason.as_str()).or_default() += 1; + continue; + } + }; + let source_file = csv.strip_prefix(&cfg.str_dir).unwrap_or(csv).to_string_lossy().into_owned(); + if cfg.apply { + str_import::upsert_profile(pool, guid, IMPORTED_FROM, PANEL, &profile, &source_file).await?; + } + rep.loaded += 1; + if sample.len() < 8 { + sample.push((name.clone(), profile.total_markers)); + } + } + + for (name, markers) in &sample { + println!("{name:<20} markers={markers}"); + } + let skipped_total: usize = rep.skipped.values().sum(); + let skipped_detail: Vec = { + let mut s: Vec<_> = rep.skipped.iter().map(|(k, v)| format!("{k}={v}")).collect(); + s.sort(); + s + }; + println!( + "\n{} dirs | {} loaded | {} kit-unresolved | {} not-in-db | {} no-csv | {} skipped [{}] | {}", + rep.dirs, + rep.loaded, + rep.unresolved_kit, + rep.not_in_db, + rep.no_csv, + skipped_total, + skipped_detail.join(", "), + if cfg.apply { "COMMITTED" } else { "preview (no writes) — pass --apply to load" } + ); + if cfg.apply { + println!("Run `decodingus-jobs run-once branch-age` to fold the STRs into branch ages."); + } + Ok(()) +} diff --git a/rust/crates/du-jobs/src/main.rs b/rust/crates/du-jobs/src/main.rs index ada35414..78be72a0 100644 --- a/rust/crates/du-jobs/src/main.rs +++ b/rust/crates/du-jobs/src/main.rs @@ -10,6 +10,7 @@ use std::sync::Arc; use std::time::Duration; mod ena; +mod ftdna_str; mod jetstream; mod publications; mod scheduler; @@ -122,8 +123,18 @@ async fn main() -> anyhow::Result<()> { let rep = du_db::tree_sample::recompute_placements(&pool, du_domain::enums::DnaType::YDna).await?; tracing::info!(placed = rep.placed, unplaced = rep.unplaced, "tree-samples-recompute complete (Y)"); } + // Bulk-load FTDNA Y-STR exports into genomics.biosample_str_profile, + // matched to the cohort by kit→subject_id. Feeds the STR-variance age + // model (run `branch-age` after). Previews unless `--apply`. + "ftdna-str" => { + let args: Vec = argv.collect(); + let cfg = ftdna_str::Config::from_env(&args).ok_or_else(|| { + anyhow::anyhow!("set FTDNA_STR_DIR + COHORT_MANIFEST") + })?; + ftdna_str::run(&pool, &cfg).await?; + } other => anyhow::bail!( - "unknown run-once job '{other}' (known: ybrowse, reconcile, yregions, branch-age, sequencer-consensus, discovery-consensus, coverage-norms, ibd-discovery-recompute, exchange-expire, tree-samples-recompute, dedup-candidates)" + "unknown run-once job '{other}' (known: ybrowse, reconcile, yregions, branch-age, ftdna-str, sequencer-consensus, discovery-consensus, coverage-norms, ibd-discovery-recompute, exchange-expire, tree-samples-recompute, dedup-candidates)" ), } return Ok(()); diff --git a/rust/crates/du-web/src/routes/tree.rs b/rust/crates/du-web/src/routes/tree.rs index a7890a14..41969bee 100644 --- a/rust/crates/du-web/src/routes/tree.rs +++ b/rust/crates/du-web/src/routes/tree.rs @@ -148,12 +148,25 @@ struct SnpSidebar { name: String, provenance: Option, variants: Vec, + /// Reconstructed ancestral Y-STR motif (phylogenetic ASR): the markers that + /// mutated on the branch leading into this node (parent→node), headlined. + str_changes: Vec, + /// Total markers in the node's reconstructed ancestral haplotype (context for + /// the change count; 0 ⇒ no STR evidence in the subtree). + str_motif_markers: usize, /// Placed non-D2C sample leaves at or below this node (capped for the sidebar). samples: Vec, /// How many more placed samples exist beyond the shown `samples` (0 ⇒ all shown). samples_more: i64, } +/// One Y-STR mutation along the branch leading into this node, e.g. DYS393 12→13. +struct StrChangeRow { + marker: String, + from: i32, + to: i32, +} + /// One placed sample row in the sidebar (label + optional paper citation). struct LeafRow { label: String, @@ -383,6 +396,24 @@ async fn snp_sidebar( None => None, }; + // Reconstructed ancestral Y-STR motif (Y only) — the per-marker mutations on + // the branch leading into this node (parent→node), plus the motif size. + let (str_changes, str_motif_markers) = if matches!(dna_type, DnaType::YDna) { + let asr = du_db::ystr::branch_str_asr(&st.pool, &name).await?; + let changes = asr + .iter() + .filter(|m| m.changed()) + .map(|m| StrChangeRow { + marker: m.marker_name.clone(), + from: m.parent_value.unwrap_or(m.ancestral_value), + to: m.ancestral_value, + }) + .collect(); + (changes, asr.len()) + } else { + (Vec::new(), 0) + }; + // Placed non-D2C sample leaves at or below this node (capped for the sidebar). let mut leaves = du_db::tree_sample::samples_under(&st.pool, &name, dna_type).await?; let samples_more = (leaves.len() as i64 - SIDEBAR_SAMPLE_CAP as i64).max(0); @@ -399,7 +430,16 @@ async fn snp_sidebar( }) .collect(); - Ok(html(&SnpSidebar { t: locale.t, name, provenance, variants, samples, samples_more })) + Ok(html(&SnpSidebar { + t: locale.t, + name, + provenance, + variants, + str_changes, + str_motif_markers, + samples, + samples_more, + })) } // ── Helpers ────────────────────────────────────────────────────────────────── diff --git a/rust/crates/du-web/templates/tree/snp_sidebar.html b/rust/crates/du-web/templates/tree/snp_sidebar.html index cdd676b0..0c14c59e 100644 --- a/rust/crates/du-web/templates/tree/snp_sidebar.html +++ b/rust/crates/du-web/templates/tree/snp_sidebar.html @@ -69,6 +69,26 @@ {% endif %} + {# Reconstructed ancestral Y-STR motif (phylogenetic ASR): the markers that + mutated on the branch into this node — the STR analogue of the SNP transitions. #} + {% if str_motif_markers > 0 %} +
+
+ {{ t.get("tree.strasr.title") }} +
+ {% if str_changes.is_empty() %} +

{{ t.get("tree.strasr.nochange") }}

+ {% else %} +
+ {% for c in str_changes %} + {{ c.marker }} {{ c.from }}→{{ c.to }} + {% endfor %} +
+ {% endif %} +
{{ str_motif_markers }} {{ t.get("tree.strasr.markers") }}
+
+ {% endif %} + {# Placed non-D2C sample leaves at or below this node (YFull-style). #} {% if !samples.is_empty() %}
diff --git a/rust/locales/en.txt b/rust/locales/en.txt index 97c57fb2..95c58168 100644 --- a/rust/locales/en.txt +++ b/rust/locales/en.txt @@ -131,6 +131,10 @@ tree.prov.aka=Also known as tree.prov.updated=Updated tree.prov.ci=95% CI tree.prov.samples=Samples (n) +tree.strasr.title=Y-STR ancestral motif (ASR) +tree.strasr.help=Reconstructed Y-STR mutation on the branch into this node (parent→node) +tree.strasr.nochange=No Y-STR change reconstructed on this branch. +tree.strasr.markers=markers reconstructed references.title=References references.search.placeholder=Search by title, journal, or DOI… diff --git a/rust/locales/es.txt b/rust/locales/es.txt index 078ee9e9..b0372a1e 100644 --- a/rust/locales/es.txt +++ b/rust/locales/es.txt @@ -83,6 +83,10 @@ tree.prov.aka=También conocido como tree.prov.updated=Actualizado tree.prov.ci=IC 95% tree.prov.samples=Muestras (n) +tree.strasr.title=Motivo Y-STR ancestral (ASR) +tree.strasr.help=Mutación Y-STR reconstruida en la rama hacia este nodo (ancestro→nodo) +tree.strasr.nochange=No se reconstruyó ningún cambio Y-STR en esta rama. +tree.strasr.markers=marcadores reconstruidos references.title=Referencias references.search.placeholder=Buscar por título, revista o DOI… diff --git a/rust/locales/fr.txt b/rust/locales/fr.txt index 6bf22c2b..fe2c5f73 100644 --- a/rust/locales/fr.txt +++ b/rust/locales/fr.txt @@ -83,6 +83,10 @@ tree.prov.aka=Aussi connu sous tree.prov.updated=Mis à jour tree.prov.ci=IC 95 % tree.prov.samples=Échantillons (n) +tree.strasr.title=Motif Y-STR ancestral (ASR) +tree.strasr.help=Mutation Y-STR reconstruite sur la branche menant à ce nœud (parent→nœud) +tree.strasr.nochange=Aucun changement Y-STR reconstruit sur cette branche. +tree.strasr.markers=marqueurs reconstruits references.title=Références references.search.placeholder=Rechercher par titre, revue ou DOI… diff --git a/rust/migrations/0053_biosample_str.sql b/rust/migrations/0053_biosample_str.sql new file mode 100644 index 00000000..e53b6848 --- /dev/null +++ b/rust/migrations/0053_biosample_str.sql @@ -0,0 +1,29 @@ +-- Locally-imported Y-STR profiles, keyed to a core.biosample (sample_guid). +-- +-- The federation path (fed.str_profile → fed.biosample, mirrored from Navigator's +-- com.decodingus.atmosphere.strProfile) covers citizen samples that publish their +-- own STRs. It can't carry the academic / D2C cohorts: those live in +-- core.biosample (accession = the cohort subject_id) and are placed in the Y tree +-- via tree.haplogroup_sample, with no atproto identity. This table is their STR +-- home — bulk-loaded from vendor exports (FTDNA Y-DNA DYS results, …) and matched +-- to the cohort by accession. +-- +-- Both sources feed the SAME STR-variance age model: du_db::ystr::build_str_inputs +-- UNIONs this table (joined live through tree.haplogroup_sample, so no stale +-- haplogroup-name copy) with fed.str_profile before grouping profiles per branch. +-- +-- Markers are stored lossless as the lexicon's strMarkerValue[] JSONB — the same +-- shape parse_markers() reads — so simple + multi-copy score and complex is +-- preserved but unscored, exactly as for the federated profiles. + +CREATE TABLE genomics.biosample_str_profile ( + sample_guid UUID PRIMARY KEY REFERENCES core.biosample(sample_guid) ON DELETE CASCADE, + source TEXT NOT NULL DEFAULT 'IMPORTED', -- DIRECT_TEST / WGS_DERIVED / IMPORTED + imported_from TEXT, -- FTDNA / YSEQ / YFULL / ... + panel TEXT, -- vendor panel label (e.g. FTDNA Y-DNA DYS) + total_markers INTEGER NOT NULL DEFAULT 0, -- non-null marker count in `markers` + markers JSONB NOT NULL DEFAULT '[]'::jsonb,-- lossless strMarkerValue[] + source_file TEXT, -- provenance: NAS path / kit dir + imported_at TIMESTAMPTZ NOT NULL DEFAULT now() +); +CREATE INDEX biosample_str_profile_source_idx ON genomics.biosample_str_profile (imported_from); diff --git a/rust/migrations/0054_haplogroup_str_asr.sql b/rust/migrations/0054_haplogroup_str_asr.sql new file mode 100644 index 00000000..3340a2b5 --- /dev/null +++ b/rust/migrations/0054_haplogroup_str_asr.sql @@ -0,0 +1,24 @@ +-- Per-branch reconstructed ancestral Y-STR motif (phylogenetic ASR). +-- +-- du_db::ystr already reconstructs each node's ancestral simple-marker motif +-- (McDonald §2.5.2 parsimony up-/down-pass over the SNP-defined tree) as an +-- in-memory input to the STR-variance age model. This table persists that motif +-- so it can be surfaced per branch: the inferred ancestral haplotype, and — by +-- diffing a node against its parent — the per-marker STR mutations along the +-- branch leading into it (the STR analogue of the SNP anc→der transitions in the +-- node sidebar). +-- +-- Distinct from tree.haplogroup_ancestral_str (which holds the per-branch MODAL +-- of *direct testers* + curator MANUAL overrides, for the signature display and +-- STR→branch prediction). This is the reconstructed ancestral state, simple +-- markers only (multi-copy/complex aren't reconstructed in v1). Full-refreshed +-- by ystr::recompute_signatures on each `branch-age` run. + +CREATE TABLE tree.haplogroup_str_asr ( + haplogroup_id BIGINT NOT NULL REFERENCES tree.haplogroup(id) ON DELETE CASCADE, + marker_name TEXT NOT NULL, + ancestral_value INTEGER NOT NULL, -- reconstructed repeat count (simple markers) + recomputed_at TIMESTAMPTZ NOT NULL DEFAULT now(), + PRIMARY KEY (haplogroup_id, marker_name) +); +CREATE INDEX haplogroup_str_asr_hg_idx ON tree.haplogroup_str_asr (haplogroup_id); From 72891e943a784e8e30e45f698d38f4cddbb07cb6 Mon Sep 17 00:00:00 2001 From: James Kane Date: Mon, 29 Jun 2026 11:53:06 -0500 Subject: [PATCH 02/13] =?UTF-8?q?feat(age):=20faithful=20McDonald=202021?= =?UTF-8?q?=20implementation=20(Eq=207/8=20build=20+=20P1=E2=80=93P5)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Restores the paper's probabilistic age model for side-by-side comparison with the pragmatic depth-quantile branch (feat/ftdna-str-aging): - SNP Eq 7/8: replace the depth-quantile estimator with the bottom-up PDF-convolution build — a node's TMRCA is the normalised product over children of (child TMRCA ⊛ branch-time) and over tester tips; P(t<0)=0 keeps a parent older than its children by construction. - Eq 9 causality: top-down reverse-convolution constraint (Pdf::convolve_sub) at nodes with ≥2 informative children (§3.4 multiplicative propagation). - P1 σ_µ (Appendix A.2.2/A.4): fold the Helgason 95% rate band into each CI in quadrature — the paper's dominant error term. - P2 STR max-reliable-age cap (A.5.2): drop STR past STR_MAX_RELIABLE_YBP where it saturates/underestimates; SNP dates the deep clades. - P4 tester-birth offset (A.1): +63±14 yr once at each tip. - P5 NRR population prior (Eq 25) mechanism, default off pending calibration. - pdf::from_weights for custom priors. Validated on the corrected FTDNA-refined tree: faithful ages reproduce accepted YFull TMRCAs within ~5% for most backbone clades (CT-M168, R-L23/L51/P310/U106/ L21/DF13/M222). See documents/proposals/aging-pipeline-audit-mcdonald2021.md. Co-Authored-By: Claude Opus 4.8 (1M context) --- rust/crates/du-db/src/age.rs | 262 +++++++++++++++++----------------- rust/crates/du-db/src/pdf.rs | 9 ++ rust/crates/du-db/src/ystr.rs | 29 ++-- 3 files changed, 159 insertions(+), 141 deletions(-) diff --git a/rust/crates/du-db/src/age.rs b/rust/crates/du-db/src/age.rs index 17787fb5..f1f7ee73 100644 --- a/rust/crates/du-db/src/age.rs +++ b/rust/crates/du-db/src/age.rs @@ -41,46 +41,11 @@ const AGE_REFRESH_GUARD: &str = "NOT COALESCE((provenance->>'age_curated')::bool /// 5 Mbp (~35% of callable MSY) floor drops only the broken ones. pub const MIN_TESTER_CALLABLE_BP: f64 = 5_000_000.0; -/// A node's TMRCA is the [`TMRCA_DEEP_QUANTILE`] of the SNP depths of its descendant -/// tips (root-to-tip path length, in callable defining + private SNPs). The raw -/// *maximum* descendant depth — the old coalescent floor's deepest-lineage rule — is a -/// max order-statistic over N tips: it grows ~log N with sample size and latches onto a -/// single overdispersed / over-called lineage (U106's deepest tip is a ~6σ academic -/// outlier), so the large R clades read ~1.65× canonical. A high quantile is -/// sample-size-robust and also recovers the central depth the confidence filter pulls -/// down (q90-on-filtered ≈ median-on-unfiltered): on the filtered counts it reproduces -/// YFull within ~10% for U106/P312/L21/M222. See -/// `~/Genomics/ytree/results/deep_clade_age_residual.md`. -pub const TMRCA_DEEP_QUANTILE: f64 = 0.90; - -/// Deep-ladder gate. A caller over-call adds a roughly *fixed* number of excess SNPs to -/// a lineage (tens), so it inflates a shallow young clade by a large fraction but a deep -/// backbone lineage (hundreds of SNPs of real coalescence) only negligibly. Where a -/// substantial fraction of a node's tips ([`DEEP_GATE_QUANTILE`]) is itself deeper than -/// [`DEEP_LADDER_SNPS`], the deep tail is real signal corroborated by many independent -/// lineages (the unary-ladder backbone, e.g. CT-M168), not a lone over-call — so the -/// node is floored to its deepest descendant (max), restoring the validated deep-backbone -/// ages where a spurious tip is a negligible fraction. The gate *condition* uses a high -/// quantile so a lone mis-placed deep tip in an otherwise-shallow clade cannot trip it; -/// the threshold sits at the empirical young/deep crossover (~R-L23 / R-M343): a young -/// clade's q97 depth (≤ ~110 SNPs) never reaches it, the backbone's always does. -pub const DEEP_LADDER_SNPS: f64 = 150.0; -pub const DEEP_GATE_QUANTILE: f64 = 0.97; - -/// Linear-interpolated quantile of an already-ascending slice. Empty → 0; singleton → -/// itself; `q` in [0,1]. -fn quantile_sorted(sorted: &[f64], q: f64) -> f64 { - match sorted.len() { - 0 => 0.0, - 1 => sorted[0], - n => { - let pos = q * (n as f64 - 1.0); - let lo = pos.floor() as usize; - let hi = pos.ceil() as usize; - sorted[lo] + (pos - lo as f64) * (sorted[hi] - sorted[lo]) - } - } -} +// NOTE (feat/faithful-mcdonald-age): the depth-quantile estimator +// (TMRCA_DEEP_QUANTILE / DEEP_LADDER_SNPS / DEEP_GATE_QUANTILE / quantile_sorted / +// global_depths) of the pragmatic branch is intentionally dropped here — this +// branch restores the paper's Eq 7/8 PDF-convolution build instead. See +// `propagate` and documents/proposals/aging-pipeline-audit-mcdonald2021.md §6. /// Independent cross-check clock from Hallast et al. 2026 (142 population-scale Y /// assemblies, BEAST v1.10.4 strict molecular clock on the X-degenerate mask): @@ -96,6 +61,44 @@ pub const HALLAST_RATE_HI: f64 = 0.86e-9; /// "Before present" reference year (radiocarbon convention) for calendar anchors. pub const PRESENT_YEAR: i32 = 1950; +/// Tester-birth offset (McDonald Appendix A.1): commercial testers average ~64 yr +/// old (95% CI 35–91), and a TMRCA should be referenced to the tester's *birth*, +/// not the sampling date. Added once at each tip (mean ≈63 yr, σ≈14 from the CI) +/// so the whole tree is in the birth frame; it does not compound up the tree +/// because the branch-time convolutions carry it upward unchanged. +pub const TESTER_BIRTH_MEAN_YBP: f64 = 63.0; +pub const TESTER_BIRTH_SIGMA: f64 = 14.0; + +/// Male net reproduction rate for the population-size prior (McDonald Eq 25: +/// `P(t|NRR) = NRR^{t/G}`). Because a person has far more distant than close +/// cousins, a uniform-prior TMRCA is skewed too young; this prior up-weights older +/// `t`. NRR varies by population/epoch (~1.0012 deep-time average; ~1.3 for some +/// recent families) and can dominate the budget for very-well-tested recent +/// lineages, so a single global value is only a first approximation — **default +/// 1.0 (off)** until per-population calibration. The mechanism is wired in +/// `recompute_combined_ages`; set >1.0 to enable. G = [`crate::ystr::GENERATION_YEARS`]. +pub const NRR_DEFAULT: f64 = 1.0; + +/// 95% relative half-width of the SNP mutation rate (Helgason CI 7.57–9.17×10⁻¹⁰ +/// about 8.33×10⁻¹⁰ ⇒ ±9.6%; Appendix A.4.2). McDonald §2.2.2/Appendix A.4: the +/// rate uncertainty is a near-common multiplicative scaling of the whole tree and +/// **dominates the error budget** for larger/older clades (~±67 yr floor even at +/// R-S781). We fold it into each age's CI in quadrature at the time-conversion +/// step — faithful to "compute in nominal mutation timescales, convert to physical +/// time afterwards" — rather than re-propagating at perturbed rates. +pub const SNP_RATE_REL_95: f64 = 0.096; + +/// Widen a `(median, lo, hi)` ybp CI by the rate uncertainty, added in quadrature +/// to each Poisson half-width (median unchanged: a common scaling shifts every node +/// together, so it is an uncertainty, not a bias, on the relative tree). +fn with_rate_uncertainty(med: i32, lo: i32, hi: i32) -> (i32, i32, i32) { + let m = med as f64; + let r = SNP_RATE_REL_95 * m; // 95% rate half-width in years at this age + let lo2 = m - (((m - lo as f64).max(0.0)).powi(2) + r * r).sqrt(); + let hi2 = m + (((hi as f64 - m).max(0.0)).powi(2) + r * r).sqrt(); + (med, lo2.max(0.0).round() as i32, hi2.round() as i32) +} + // ── PDF-based tree propagation (McDonald 2021 §2.2, Eq 5–8) ─────────────────── // // The SNP age of a clade is built bottom-up: a node's TMRCA is the product over @@ -188,111 +191,94 @@ fn post_order(clades: &[Clade]) -> Vec { order } -/// Root-to-node cumulative branch SNPs (`gd[i]` = Σ branch_snps from a root down to -/// `i`), filled parents-before-children by walking `post` in reverse. -fn global_depths(clades: &[Clade], post: &[usize]) -> Vec { - let mut gd = vec![0.0f64; clades.len()]; - for &i in post.iter().rev() { - for &ch in &clades[i].children { - gd[ch] = gd[i] + clades[ch].branch_snps as f64; - } - } - gd -} - -/// Compute every clade's TMRCA + formed-age PDFs. A node's TMRCA (in SNPs, relative to -/// the node) is the larger of two terms: -/// • its own [`TMRCA_DEEP_QUANTILE`] over the pooled SNP depths of all descendant -/// tips (each tip's root-to-tip path length minus the node's own depth); and -/// • a **corroboration floor**: the *second*-deepest of its children's -/// (q90-robustified) ages. ≥2 independent deep sub-clades — a real early split, as -/// on the deep backbone — floor the node deep; a *single* deep child does NOT (it's -/// a lone over-deep lineage the q90 already discounts). Because each child's age is -/// itself robustified, a lone deep tip is discounted at every level before it can -/// spuriously corroborate. -/// The result is wrapped in a Poisson age PDF over the callable denominator; the formed -/// age is that TMRCA convolved with the node's own branch time. A node with no -/// descendant tip yields `None`. +/// Compute every clade's TMRCA + formed-age PDFs — **faithful McDonald Eq 7/8**. +/// +/// Bottom-up (Eq 8): a node's TMRCA PDF is the normalised **product** of its +/// independent lines of evidence, each lifted into the node's frame — +/// • each child sub-clade: `P(t_child) ⊛ P(t_{child→node})` (the child's TMRCA +/// convolved with the parent→child branch-time PDF, Eq 7); and +/// • each direct tester tip: a Poisson age `P(t | m_private, b̄)` over the node's +/// callable bp (tester birth ≈ present; the A.1 offset is added in `recompute`). +/// `Pdf::convolve` keeps `P(t<0)=0` (the branch time is non-negative), so a parent +/// is older than its children *by construction* in the common case. /// -/// This replaces the bottom-up product + coalescent-floor propagation, whose -/// max-descendant floor grew ~log N with sample size and whose `truncate_below` -/// renormalisation compounded that inflation up each spine node (see -/// [`TMRCA_DEEP_QUANTILE`]). +/// Top-down (Eq 9, §3.4 multiplicative propagation): where the stochastic SNP +/// counts still leave a child older than its parent, the parent acts as a +/// semi-independent constraint — `P(t_child | parent) = P(t_parent) ⊟ P(branch)` +/// (`convolve_sub`, P(t<0)=0) — combined multiplicatively into the child. Only +/// applied at nodes with ≥2 informative children, so the parent is not dominated +/// by the single child it would constrain (the paper's circularity caveat). +/// +/// This is the paper-faithful counterpart to the pragmatic depth-quantile +/// estimator on `feat/ftdna-str-aging`; the robustness trade-offs differ (the +/// product is sharper but more sensitive to a single over-called sub-clade — the +/// very effect §2.3/Appendix A.4.1 and the causality pass address). pub fn propagate(clades: &[Clade], mu: f64, res: f64, max_age: f64) -> Vec> { let post = post_order(clades); - let gd = global_depths(clades, &post); - // Per node, the ascending multiset of descendant-tip absolute depths (gd + the tip's - // private SNPs). Built bottom-up: each node drains its children's lists and appends - // its own direct testers, keeping the merged list for its parent (each list moves up - // exactly once). - let mut tip_depths: Vec> = vec![Vec::new(); clades.len()]; - let mut tmrca_snps: Vec> = vec![None; clades.len()]; + let branch = |c: usize| branch_time(clades, c, mu, res, max_age); + // Tester-birth offset (A.1) — convolved once into each tip, not per level. + let birth = Pdf::gaussian_on(TESTER_BIRTH_MEAN_YBP, TESTER_BIRTH_SIGMA, res, max_age); + let mut tmrca: Vec> = vec![None; clades.len()]; + + // ── Eq 8: bottom-up product build ──────────────────────────────────────── for &i in &post { - let mut depths = Vec::new(); + let mut factors: Vec = Vec::new(); for &ch in &clades[i].children { - depths.append(&mut tip_depths[ch]); - } - for &s in &clades[i].tester_snps { - depths.push(gd[i] + s as f64); + if let Some(ct) = &tmrca[ch] { + factors.push(ct.convolve(&branch(ch))); // child TMRCA ⊛ branch (Eq 7) + } } - if depths.is_empty() { - continue; + for &m in &clades[i].tester_snps { + // Tester age (Eq 3) referenced to the tester's birth (Eq 8 single-tester + // case p_k = P(t_c|m_c) ⊛ P(t_b); A.1 birth PDF). + factors.push(Pdf::poisson_on(m, clades[i].callable_bp, mu, res, max_age).convolve(&birth)); } - depths.sort_by(f64::total_cmp); - let local = (quantile_sorted(&depths, TMRCA_DEEP_QUANTILE) - gd[i]).max(0.0); - // Children's q90-robustified ages, lifted to this node's frame (+ branch). - let mut child_ages: Vec = clades[i] - .children - .iter() - .filter_map(|&ch| tmrca_snps[ch].map(|t| t + clades[ch].branch_snps as f64)) - .collect(); - // Corroboration floor = the 2nd-deepest child (depth reached by ≥2 independent - // sub-clades); a single deep child is not enough. - let corroborated = if child_ages.len() >= 2 { - child_ages.sort_by(|a, b| b.total_cmp(a)); - child_ages[1] - } else { - 0.0 - }; - // Deep-ladder gate (see [`DEEP_LADDER_SNPS`]): where a substantial fraction of - // tips is itself deep — real backbone, not a lone over-call — floor to the deepest - // descendant. Condition on q97 (robust to a single mis-placed deep tip), value = - // max depth (the validated deep-backbone estimate). - let gate_q = (quantile_sorted(&depths, DEEP_GATE_QUANTILE) - gd[i]).max(0.0); - let deep = if gate_q >= DEEP_LADDER_SNPS { - (depths[depths.len() - 1] - gd[i]).max(0.0) - } else { - 0.0 - }; - tmrca_snps[i] = Some(local.max(corroborated).max(deep)); - tip_depths[i] = depths; + tmrca[i] = pdf_product(&factors); } - // Top-down monotonicity: a child's TMRCA cannot exceed its parent's (its MRCA is more - // recent than the parent's). Where the robust estimate left an over-deep child under a - // corrected younger parent — an over-called subclade the parent's pooled q90 already - // discounted — clamp it to the parent, propagating the de-inflation down the subtree - // and removing the parent/child age inversions. Backbone order (deep parent > deep - // child) is untouched. Reverse post-order visits parents before children. + + // ── Eq 9: top-down causality constraint (parents before children) ───────── for &i in post.iter().rev() { - let Some(pt) = tmrca_snps[i] else { continue }; - for &ch in &clades[i].children { - if let Some(ct) = tmrca_snps[ch] { - if ct > pt { - tmrca_snps[ch] = Some(pt); - } - } + let Some(parent) = tmrca[i].clone() else { continue }; + let informative: Vec = + clades[i].children.iter().copied().filter(|&c| tmrca[c].is_some()).collect(); + if informative.len() < 2 { + continue; // can't constrain a child the parent's age depends entirely on + } + for &ch in &informative { + let implied = parent.convolve_sub(&branch(ch)); // P(t_parent) ⊟ P(branch) + let child = tmrca[ch].as_ref().unwrap(); + let refined = child.multiply(&implied); + // If the child wholly violated causality (disjoint from the parent + // constraint), pin it to the constraint rather than annihilating it. + tmrca[ch] = Some(if refined.total() > 0.0 { refined } else { implied }); } } + (0..clades.len()) .map(|i| { - let snps = tmrca_snps[i]?; - let tmrca = Pdf::poisson_on(snps.round() as i64, clades[i].callable_bp, mu, res, max_age); - let formed = tmrca.convolve(&branch_time(clades, i, mu, res, max_age)); - Some(CladeAge { tmrca, formed }) + let t = tmrca[i].clone()?; + let formed = t.convolve(&branch(i)); + Some(CladeAge { tmrca: t, formed }) }) .collect() } +/// Normalised product of independent age PDFs (Eq 8 / Eq 1). `None` if empty. +fn pdf_product(factors: &[Pdf]) -> Option { + factors.split_first().map(|(first, rest)| rest.iter().fold(first.clone(), |a, f| a.multiply(f))) +} + +/// McDonald Eq 25 population-size prior as a PDF: `mass[i] ∝ nrr^{age_i / g}` +/// (normalised in `multiply`). With `nrr > 1` older ages carry more prior weight, +/// correcting the young skew from there being more distant than close cousins. +fn nrr_prior_pdf(nrr: f64, g: f64, res: f64, max_age: f64) -> Pdf { + let bins = (max_age / res).round() as usize + 1; + let ln = nrr.ln(); + // exp((t/g)·ln nrr); kept on the same grid as the age PDFs. + let mass: Vec = (0..bins).map(|i| ((i as f64 * res / g) * ln).exp()).collect(); + Pdf::from_weights(res, mass) +} + /// SNPs in recurrent / FP-prone sequence are masked from age counting — they sit /// outside the callable denominator (`y_xdegen+y_ampliconic+y_palindromic`) and the /// paper excises recurrent regions self-consistently (Appendix A.2/A.3). Ampliconic @@ -627,6 +613,7 @@ pub async fn recompute_combined_ages(pool: &PgPool) -> Result crate::ystr::STR_MAX_RELIABLE_YBP { continue; } @@ -645,13 +632,17 @@ pub async fn recompute_combined_ages(pool: &PgPool) -> Result = BTreeSet::new(); node_set.extend(snp_pdf.keys().chain(gen_pdf.keys()).chain(str_pdf.keys()).copied()); + // Eq 25 population-size prior (NRR^{t/G}); None when disabled (NRR=1.0 ⇒ no-op). + let nrr_prior = (NRR_DEFAULT != 1.0).then(|| { + nrr_prior_pdf(NRR_DEFAULT, crate::ystr::GENERATION_YEARS, TREE_RESOLUTION_YEARS, TREE_MAX_AGE_YEARS) + }); let mut combined_writes: Vec<(i64, i32, i32, i32, i32)> = Vec::new(); // (hg, med, lo, hi, n_terms) for hg in node_set { let factors: Vec<&Pdf> = [snp_pdf.get(&hg), gen_pdf.get(&hg), str_pdf.get(&hg)].into_iter().flatten().collect(); let Some((first, rest)) = factors.split_first() else { continue }; let product = rest.iter().fold((*first).clone(), |acc, f| acc.multiply(f)); - let combined = if product.total() > 0.0 { + let mut combined = if product.total() > 0.0 { product } else { let params: Vec<(f64, f64)> = factors.iter().map(|p| pdf_gaussian_params(p)).collect(); @@ -660,6 +651,9 @@ pub async fn recompute_combined_ages(pool: &PgPool) -> Result (*first).clone(), } }; + if let Some(prior) = &nrr_prior { + combined = combined.multiply(prior); // Eq 25: up-weight older t + } let (med, lo, hi) = combined.ci95(); combined_writes.push((hg, med.round() as i32, lo.round() as i32, hi.round() as i32, factors.len() as i32)); } @@ -711,7 +705,8 @@ pub async fn recompute_combined_ages(pool: &PgPool) -> Result Result) -> Pdf { + let mut p = Pdf { res, mass: weights }; + p.normalize(); + p + } + /// A near-delta PDF at a precisely known age (e.g. a proven MRCA birth year). pub fn point(years: f64) -> Pdf { let mut pdf = Pdf::zeros(RESOLUTION_YEARS, MAX_AGE_YEARS); diff --git a/rust/crates/du-db/src/ystr.rs b/rust/crates/du-db/src/ystr.rs index f0a38305..f4cbdd44 100644 --- a/rust/crates/du-db/src/ystr.rs +++ b/rust/crates/du-db/src/ystr.rs @@ -174,6 +174,17 @@ pub const GENERATION_YEARS: f64 = 33.0; /// TMRCAs is negligible. const STR_M_MAX: i64 = 30; +/// Maximum TMRCA (years) at which a Y-STR estimate is trusted in the combined age +/// (McDonald Appendix A.5.2: "a maximum time limit should be associated with the +/// inclusion of Y-STRs"). Beyond a few thousand years a single marker's observed +/// distance saturates (back/parallel/multi-step mutations), and reconstructed +/// ancestral motifs miss events, so STR systematically *underestimates* deep +/// clades (Example 3: 2560 vs 4000 yr). Rather than the paper's harder-to-calibrate +/// Gaussian taper on `P(g|m)`, we apply the simpler stated option: a node whose STR +/// median exceeds this limit does not contribute its STR term to the Eq-1 product +/// (the SNP clock dates it); STR informs the recent clades where it is most useful. +pub const STR_MAX_RELIABLE_YBP: f64 = 10_000.0; + // ── Multi-step P(g|m): McDonald 2021 §2.5.3–2.5.4 (Eq 14–23, Table 1) ────────── // // A Y-STR's observed genetic distance g (|obs − ancestral|) underdetermines the @@ -725,15 +736,12 @@ async fn build_str_inputs(pool: &PgPool) -> Result { } /// Minimum direct STR testers on a node for its propagated TMRCA to enter the -/// combined-age product. STR ages are only trustworthy where there is genuine -/// within-clade genetic *diversity*: with 0 testers a node's age is pure motif -/// reconstruction, and with 1 the reconstructed ancestral motif equals that lone -/// tester — both make the genetic distance ≈0, so the STR PDF collapses to a -/// spuriously young, spuriously narrow estimate (McDonald §2.5.2: "TMRCAs become -/// underestimated if mutations are missed") that then dominates the sharp SNP -/// term. Requiring ≥2 testers keeps STR where the paper says it is most useful -/// (well-populated clades) and lets SNP date the rest; the COMBINED causality -/// pass ([`crate::age`]) is the hard backstop. +/// combined-age product. STR ages need genuine within-clade genetic *diversity*: +/// with 0–1 testers the reconstructed ancestral motif ≈ the tester(s), the genetic +/// distance ≈0, and the STR PDF collapses to a spuriously young, narrow estimate +/// (McDonald §2.5.2) that would dominate the sharp SNP term. ≥2 keeps STR where the +/// paper says it is most useful; [`STR_MAX_RELIABLE_YBP`] additionally caps it to +/// recent time, and the Eq 9 causality pass is the hard backstop. pub const MIN_STR_TESTERS_FOR_COMBINE: usize = 2; /// Per-node tree-propagated STR TMRCA PDFs (keyed by haplogroup id) on the given @@ -748,6 +756,9 @@ pub async fn str_tmrca_pdfs(pool: &PgPool, res: f64, max_age: f64) -> Result= MIN_STR_TESTERS_FOR_COMBINE) .filter_map(|i| ages[i].clone().map(|p| (inp.ids[i], p))) + // Max-reliable-age cap (A.5.2): drop STR where it saturates and would + // underestimate; the SNP clock dates the deep clades. + .filter(|(_, p)| p.median() <= STR_MAX_RELIABLE_YBP) .collect()) } From 5aa17b6b267dc73731edaa8eef7808ece5ee1524 Mon Sep 17 00:00:00 2001 From: James Kane Date: Mon, 29 Jun 2026 18:29:58 -0500 Subject: [PATCH 03/13] feat(age): per-sample tester Poisson (region-consistent coverage) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Each terminal tester now contributes a Poisson age over its OWN callable bp rather than a shared clade/joint denominator. A tester's private SNPs are its own calls, so by construction they lie inside its own mask — using a uniform denominator (or the Eq-4 second-highest coverage) instead is a global rescale that pushed every clade ~28% older, diverging from generally-accepted ages. `Clade.tester_snps` becomes `Vec<(count, callable_bp)>`. Defining/branch SNPs keep the joint-call mask they were ascertained over (region-consistent with the SNP-region positive filter); only the tester denominator goes per-sample. Net effect: deep branch-dominated nodes (BT, CT) are preserved within ~1%, while recent tester-driven clades (M222, U106) move a few percent older to reflect their real per-sample coverage. Co-Authored-By: Claude Opus 4.8 (1M context) --- rust/crates/du-db/src/age.rs | 66 +++++++++++++++++------------------- 1 file changed, 32 insertions(+), 34 deletions(-) diff --git a/rust/crates/du-db/src/age.rs b/rust/crates/du-db/src/age.rs index f1f7ee73..9701b705 100644 --- a/rust/crates/du-db/src/age.rs +++ b/rust/crates/du-db/src/age.rs @@ -120,14 +120,17 @@ pub struct Clade { /// the branch time when this node feeds its parent, and its own "formed" age. /// 0 for a root. pub branch_snps: i64, - /// Effective callable bp (`b̄`) over which this clade's SNPs are counted. + /// Callable bp (`b̄`) for this clade's **defining (branch) SNPs** — the joint-call + /// mask region those SNPs were ascertained over. Testers carry their own bp (below). pub callable_bp: f64, /// Child clade indices. pub children: Vec, - /// Private-SNP counts of testers sitting directly on this node (terminal tips); - /// each contributes a Poisson age factor (tester birth ≈ present is omitted as - /// a negligible offset). - pub tester_snps: Vec, + /// Testers sitting directly on this node (terminal tips), each as + /// `(private_SNP_count, that sample's callable bp)`. Per McDonald §2.2.3 + /// region-consistency, each contributes a Poisson age over its **own** coverage + /// — its private SNPs are its own calls, so they lie in its own mask — rather + /// than a shared clade/joint denominator (tester birth handled in `propagate`). + pub tester_snps: Vec<(i64, f64)>, } /// A clade's computed age PDFs. @@ -228,10 +231,11 @@ pub fn propagate(clades: &[Clade], mu: f64, res: f64, max_age: f64) -> Vec Result<(Vec, Vec), DbError> } // Per-build joint-call callable mask (genomics.y_callable_interval), if loaded. - // When present it makes the numerator and denominator region-consistent: - // • denominator — the UNIFORM mask size (the region ASR branch counts were - // ascertained over), replacing the per-sample coverage average; and - // • numerator — a POSITIVE filter so only SNPs INSIDE the mask are counted - // (the recurrent mask is negative-only and lets non-callable SNPs through). - // Absent (empty table) → falls back to the prior per-sample behaviour. + // Used here as the NUMERATOR's SNP-region positive filter — only SNPs INSIDE the + // mask are counted (the recurrent mask is negative-only and lets non-callable + // SNPs through). The DENOMINATOR (b̄) is the per-clade Eq-4 coverage computed + // below, not the uniform mask size. Absent (empty table) → no positive filter. let mask_bp: Option = sqlx::query_scalar( "SELECT sum(upper(span) - lower(span))::float8 FROM genomics.y_callable_interval", ) @@ -394,7 +396,6 @@ async fn build_clades(pool: &PgPool) -> Result<(Vec, Vec), DbError> )) .fetch_all(pool) .await?; - let (mut bp_sum, mut bp_cnt) = (vec![0.0f64; ids.len()], vec![0u32; ids.len()]); for (hg, snps, b) in testers { // Skip sliver-coverage samples: they divide their SNPs by a tiny callable // denominator and produce impossible (hundreds-of-ky) ages that floor the spine. @@ -402,14 +403,13 @@ async fn build_clades(pool: &PgPool) -> Result<(Vec, Vec), DbError> // deep/rare singleton lineages — A/B/Q/N — not caller over-calls, so a // clade-blind cap would wrongly youthen them. See branch_snp_qc_report.md.) if let (Some(&i), true) = (idx.get(&hg), b >= MIN_TESTER_CALLABLE_BP) { - clades[i].tester_snps.push(snps); - bp_sum[i] += b; - bp_cnt[i] += 1; + // Per-sample: keep the tester's own callable bp for its Poisson denominator. + clades[i].tester_snps.push((snps, b)); } } - // Representative b̄ per node: mean of its testers' callable bp, else the - // catalog-wide mean (so SNP-less internal branches still get a branch time). + // Catalog-wide mean coverage — fallback denominator for a defining-SNP branch + // with no joint mask loaded. let default_b: f64 = sqlx::query_scalar::<_, Option>(&format!( "SELECT avg({cbp})::float8 FROM genomics.biosample_callable_loci cl WHERE cl.chromosome IN ('chrY','Y')" )) @@ -417,14 +417,12 @@ async fn build_clades(pool: &PgPool) -> Result<(Vec, Vec), DbError> .await? .filter(|b| *b > 0.0) .unwrap_or(15_000_000.0); + + // Branch (defining-SNP) denominator: the joint-call mask the ASR branch counts + // were ascertained over (uniform region, region-consistent with the SNP-region + // positive filter above). Testers do NOT use this — each carries its own bp. for i in 0..ids.len() { - // Uniform joint-call mask size when loaded (ASR branch counts are ascertained - // over that fixed region); else the prior per-sample mean (or catalog default). - clades[i].callable_bp = match mask_bp { - Some(b) if b > 0.0 => b, - _ if bp_cnt[i] > 0 => bp_sum[i] / bp_cnt[i] as f64, - _ => default_b, - }; + clades[i].callable_bp = mask_bp.filter(|b| *b > 0.0).unwrap_or(default_b); } Ok((clades, ids)) @@ -846,7 +844,7 @@ mod tests { #[test] fn tmrca_of_single_tester_is_poisson_mode() { - let clades = vec![Clade { branch_snps: 0, callable_bp: B, children: vec![], tester_snps: vec![3] }]; + let clades = vec![Clade { branch_snps: 0, callable_bp: B, children: vec![], tester_snps: vec![(3, B)] }]; let ages = propagate(&clades, MU, RES, MAXA); let tmrca = &ages[0].as_ref().unwrap().tmrca; // Poisson mode m/(b·µ) = 300 yr, shifted by the tester-birth offset (A.1). @@ -859,7 +857,7 @@ mod tests { // parent(0) → child(1); child has 2 private SNPs and is 1 SNP below parent. let clades = vec![ Clade { branch_snps: 0, callable_bp: B, children: vec![1], tester_snps: vec![] }, - Clade { branch_snps: 1, callable_bp: B, children: vec![], tester_snps: vec![2] }, + Clade { branch_snps: 1, callable_bp: B, children: vec![], tester_snps: vec![(2, B)] }, ]; let ages = propagate(&clades, MU, RES, MAXA); let parent = ages[0].as_ref().unwrap(); @@ -877,7 +875,7 @@ mod tests { let mut clades = vec![Clade { branch_snps: 0, callable_bp: B, children: (1..=5).collect(), tester_snps: vec![] }]; for _ in 0..5 { - clades.push(Clade { branch_snps: 1, callable_bp: B, children: vec![], tester_snps: vec![50] }); + clades.push(Clade { branch_snps: 1, callable_bp: B, children: vec![], tester_snps: vec![(50, B)] }); } let parent = propagate(&clades, MU, RES, MAXA)[0].as_ref().unwrap().tmrca.median(); assert!(parent > 4500.0, "consistent deep clade parent {parent} stays deep"); @@ -896,9 +894,9 @@ mod tests { tester_snps: vec![], }]; for _ in 0..9 { - clades.push(Clade { branch_snps: 1, callable_bp: B, children: vec![], tester_snps: vec![2] }); + clades.push(Clade { branch_snps: 1, callable_bp: B, children: vec![], tester_snps: vec![(2, B)] }); } - clades.push(Clade { branch_snps: 1, callable_bp: B, children: vec![], tester_snps: vec![60] }); + clades.push(Clade { branch_snps: 1, callable_bp: B, children: vec![], tester_snps: vec![(60, B)] }); let parent = propagate(&clades, MU, RES, MAXA)[0].as_ref().unwrap().tmrca.median(); assert!(parent < 3000.0, "robust parent {parent} not floored to the lone deep outlier"); assert!(parent > 200.0, "parent {parent} still older than the shallow tips"); @@ -979,7 +977,7 @@ mod tests { // (a) build_clades: het-masking + structure. let (clades, ids) = build_clades(&pool).await.unwrap(); let at = |id: i64| ids.iter().position(|&x| x == id).unwrap(); - assert_eq!(clades[at(leaf)].tester_snps, vec![5], "het private SNP masked → 5 counted"); + assert_eq!(clades[at(leaf)].tester_snps, vec![(5, 12_500_000.0)], "het private SNP masked → 5 counted, with its callable bp"); assert_eq!(clades[at(leaf)].branch_snps, 3, "het defining SNP masked → 3"); assert!(clades[at(mid)].children.contains(&at(leaf))); assert!(clades[at(root)].children.contains(&at(mid))); From 1fd1c072c63b97530aa6593be05dd3c54363daba Mon Sep 17 00:00:00 2001 From: James Kane Date: Mon, 29 Jun 2026 20:58:23 -0500 Subject: [PATCH 04/13] fix(age): zero age-countable-SNP nodes are age-transparent, not short MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A node can be validly defined for PLACEMENT by a marker the SNP clock can't count — a palindromic SNP that gene-converts onto both arms (e.g. ZZ11), or a recurrent/non-callable site the masks exclude. Such a node legitimately has branch_snps = 0, but it is NOT a short branch: we have no clock information about that edge. The bug: branch_time fed m=0 into Pdf::poisson_on, whose m=0 case is the exponential exp(-t·b·µ) (~1/(b·µ) ≈ 90 yr mean), not a point mass. So every stacked clock-unstable generation (Z46516 → ZZ11 → …) silently padded ~90 yr onto every ancestor above it, over-aging the deep backbone of richly sub-divided clades (P312 had two such phantom generations pinning it). Fix: when branch_snps == 0, the branch time is a point mass at t=0 — a zero-length, age-transparent edge. The node stays in the tree (placement intact); its descendants' ages propagate straight through to the parent with no per-generation padding. General: applies to any palindromic/recurrent/ non-callable-defined node, not just ZZ11. Effect on the corrected tree: deep clock-stable backbone unchanged (BT/CT move 0-2 yr), while the P312/U106 backbone shifts ~60-100 yr younger (P312 5090→4994, U106 5115→5023) as the phantom padding is removed. Z46516 and ZZ11 are now coincident with P312, as they should be. Co-Authored-By: Claude Opus 4.8 (1M context) --- rust/crates/du-db/src/age.rs | 55 ++++++++++++++++++++++++++++++++++-- 1 file changed, 52 insertions(+), 3 deletions(-) diff --git a/rust/crates/du-db/src/age.rs b/rust/crates/du-db/src/age.rs index 9701b705..95e9bd1b 100644 --- a/rust/crates/du-db/src/age.rs +++ b/rust/crates/du-db/src/age.rs @@ -116,9 +116,12 @@ fn with_rate_uncertainty(med: i32, lo: i32, hi: i32) -> (i32, i32, i32) { /// One clade (haplogroup node) of the propagation input. #[derive(Debug, Clone, Default)] pub struct Clade { - /// SNPs on the edge from this node's parent down to it (`m_{parent→node}`): - /// the branch time when this node feeds its parent, and its own "formed" age. - /// 0 for a root. + /// **Age-countable** SNPs on the edge from this node's parent down to it + /// (`m_{parent→node}`): the branch time when this node feeds its parent, and its + /// own "formed" age. 0 for a root — and 0 for any node defined only by markers the + /// clock can't count (palindromic/recurrent/non-callable, valid for placement but + /// excluded from counting), which [`branch_time`] treats as a zero-length, + /// age-transparent edge rather than a short branch. pub branch_snps: i64, /// Callable bp (`b̄`) for this clade's **defining (branch) SNPs** — the joint-call /// mask region those SNPs were ascertained over. Testers carry their own bp (below). @@ -153,7 +156,24 @@ pub const TREE_RESOLUTION_YEARS: f64 = 50.0; pub const TREE_MAX_AGE_YEARS: f64 = 500_000.0; /// Branch-time PDF for clade `x`: `P(t | m_branch)` over its callable bp. +/// +/// A node with **zero age-countable SNPs** (`branch_snps == 0`) is age-transparent, +/// NOT a short branch. Such a node is real and validly defined for *placement* — but +/// by a marker the SNP clock can't count: a palindromic SNP that gene-converts onto +/// both arms (e.g. ZZ11), or a recurrent/non-callable site the masks exclude. We have +/// no clock information about that edge, so its branch time is a point mass at t=0 (a +/// zero-length edge), making the node coincident with its parent for age propagation. +/// Its descendants' ages then flow straight through to the parent. +/// +/// The alternative — `poisson_on(0, …)` — is an *exponential* `exp(−t·b·µ)` with a +/// ~`1/(b·µ)` ≈ 90 yr mean (only the t=0 bin gets `ln 1`; every later bin keeps the +/// decaying tail), so each stacked clock-unstable generation (Z46516 → ZZ11 → …) +/// silently pads ~90 yr onto every ancestor above it. That is the "missing SNP +/// collapses the algorithm" smell: treat `m=0` as "no information," not "short." fn branch_time(clades: &[Clade], x: usize, mu: f64, res: f64, max_age: f64) -> Pdf { + if clades[x].branch_snps <= 0 { + return Pdf::from_weights(res, vec![1.0]); // δ(t=0): zero-length, age-transparent + } Pdf::poisson_on(clades[x].branch_snps, clades[x].callable_bp, mu, res, max_age) } @@ -868,6 +888,35 @@ mod tests { assert!(child.formed.median() > child.tmrca.median(), "formed > tmrca"); } + #[test] + fn zero_snp_node_is_age_transparent_not_padding() { + // A node defined only by a clock-unstable marker (a palindromic SNP like ZZ11, + // valid for placement but not age-countable) has branch_snps = 0. It must be a + // zero-length, age-transparent edge — coincident with its parent — NOT an + // exponential ~1/(b·µ) branch that pads ~90 yr onto every ancestor above it. + // root(0) → mid(branch) → tip(1, deep 50-SNP tester). + let chain = |mid_branch: i64| { + let clades = vec![ + Clade { branch_snps: 0, callable_bp: B, children: vec![1], tester_snps: vec![] }, + Clade { branch_snps: mid_branch, callable_bp: B, children: vec![2], tester_snps: vec![] }, + Clade { branch_snps: 1, callable_bp: B, children: vec![], tester_snps: vec![(50, B)] }, + ]; + propagate(&clades, MU, RES, MAXA) + }; + let transp = chain(0); + let (root_t, mid_t) = + (transp[0].as_ref().unwrap().tmrca.median(), transp[1].as_ref().unwrap().tmrca.median()); + // 0-SNP mid carries no branch time: root is coincident with mid (no padding). + assert!((root_t - mid_t).abs() <= RES, "0-SNP mid age-transparent: root {root_t} ≈ mid {mid_t}"); + // Give mid a real countable branch → the parent is now strictly older than mid… + let real = chain(5); + let (root_r, mid_r) = + (real[0].as_ref().unwrap().tmrca.median(), real[1].as_ref().unwrap().tmrca.median()); + assert!(root_r > mid_r + RES, "real branch ages the parent: root {root_r} > mid {mid_r}"); + // …and the transparent parent is younger than the padded one (no spurious age). + assert!(root_t < root_r, "transparent root {root_t} younger than real-branch root {root_r}"); + } + #[test] fn consistent_deep_clade_stays_deep() { // All five sub-clades genuinely deep (50 private SNPs ≈ 5000 yr): the deep age From 3904252892ff7f897a57153f0c309a6e8cbd9e0e Mon Sep 17 00:00:00 2001 From: James Kane Date: Tue, 30 Jun 2026 06:32:35 -0500 Subject: [PATCH 05/13] fix(denovo): exempt reverse-polarity SNPs from the monophyletic gate MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The regenerated FTDNA ingests carry per-variant confidence flags; the loader excludes !jointConfirmed || !monophyletic from branch-age counting. But the `monophyletic` flag is invalid for reverse-polarity SNPs: the reference (CHM13/HG002) is haplogroup J, deep inside CT, so every backbone SNP above J (BT/CT/F/CF) carries the derived allele AS the reference. Their derived carriers are reference-matching (invisible to variant calling) while only the ancestral A/B outgroup shows a variant — a paraphyletic set the monophyly test (computed from variant calls) falsely flags non-monophyletic. Gating on it zeroed the ENTIRE deep backbone (CT 289/289, F 163/163, BT 8/8 flagged), which — with the age engine's zero-SNP-node transparency — collapsed every TMRCA below CT to a single ~16.8 kya depth (the whole tree descends from CT, so the parent≥child constraint capped all of it). Fix: exempt reverse-polarity SNPs from the monophyletic clause (jointConfirmed still applies to all; monophyletic still excludes genuine forward homoplasy). The deep backbone counts again — BT 95.7 kya, CT 63.2 kya, A-M31 11.6 kya (matches YFull). The joint VCF confirms these SNPs are real (GQ=99, clean AC=5 synapomorphies); the flag was wrong, not the calls. The real fix belongs upstream: the ASR/Fitch monophyly computation should be polarity-aware. Documented for the tree team in documents/proposals/denovo-ingest-confidence-flags.md. The INDEL artifact remains unverified — staying on the SNP artifact. Co-Authored-By: Claude Opus 4.8 (1M context) --- .../denovo-ingest-confidence-flags.md | 83 +++++++++++++++++++ rust/crates/du-db/src/denovo.rs | 54 +++++++++++- 2 files changed, 136 insertions(+), 1 deletion(-) create mode 100644 documents/proposals/denovo-ingest-confidence-flags.md diff --git a/documents/proposals/denovo-ingest-confidence-flags.md b/documents/proposals/denovo-ingest-confidence-flags.md new file mode 100644 index 00000000..7effa603 --- /dev/null +++ b/documents/proposals/denovo-ingest-confidence-flags.md @@ -0,0 +1,83 @@ +# De-novo ingest confidence flags — issues for the tree team + +**Status:** findings from age-pipeline work (2026-06-30, `feat/faithful-mcdonald-age`). +**Audience:** ytree pipeline maintainers (`~/Genomics/ytree/bin`). +**TL;DR:** The regenerated FTDNA ingest artifacts carry `jointConfirmed` / `monophyletic` +per-variant confidence flags. The AppView loader excludes flagged SNPs from branch-age +counting. One of those flags (`monophyletic`) is computed in a way that is **invalid for +reverse-polarity SNPs**, which silently zeroes the **entire deep backbone** (BT/CT/F) and +collapses every TMRCA below CT. A loader-side workaround is in place; the real fix belongs +upstream in the ASR / flagging stages. + +## Where this bites + +The loader (`rust/crates/du-db/src/denovo.rs::DenovoVariant::low_confidence`) maps: + +``` +low_confidence = !jointConfirmed || !monophyletic (missing flag → confident) +``` + +and the age engine (`du-db/src/age.rs::build_clades`) excludes `low_confidence` defining +SNPs from a node's branch count. Combined with the age engine's *zero-age-countable-SNP +nodes are age-transparent* rule (a node with 0 countable SNPs is treated as a zero-length, +pass-through edge), any node whose defining SNPs are **all** flagged collapses onto its +parent. When that happens to the deep backbone, the whole tree below it inherits the +collapse via the parent≥child causality constraint. + +## Issue 1 — `monophyletic=false` on reverse-polarity backbone SNPs (the CHM13-is-J bug) + +The reference is **CHM13 / HG002, whose Y is haplogroup J** — which sits deep *inside* CT. +So every backbone SNP **above J** (CF, F, CT, BT, and most of the deep A spine) carries the +**derived** allele *as the reference* → `polarity: "reverse"`. + +For a reverse-polarity SNP the derived carriers are reference-matching and therefore +**invisible to variant calling**, while only the ancestral outgroup (A/B lineages) shows a +variant call. A monophyly test computed from variant calls then sees the *ancestral* +carriers — a paraphyletic set — and flags the SNP `monophyletic=false`. This is a systematic +artifact of the inversion, not real homoplasy. + +Observed in `chrY.ftdna.refined.ingest.json` and `chrY.ftdna.refined.indel.ingest.json`: + +| node | defining SNPs | polarity | monophyletic=false | would-count (pre-fix) | +|---------|---------------|------------|--------------------|-----------------------| +| CT-M168 | 289 | all reverse| 286 | **0** | +| F-M89 | 163 | 162 reverse| 163 | **0** | +| BT-M42 | 8 | all reverse| 8 | **0** | +| A-V168 | 205 | all reverse| 197 | 8 | + +These are foundational markers (M168, M89, M42 themselves). Zeroing them drove the whole +tree to a single collapsed depth (every TMRCA pinned at CT's collapsed ~16.8 kya). + +**Loader workaround (shipped):** exempt reverse-polarity SNPs from the `monophyletic` gate. +`jointConfirmed` still applies to all; `monophyletic` still excludes genuine *forward* +homoplasy. + +**Upstream fix wanted:** the monophyly computation (Fitch / ASR, `bin/82b_asr.py` + +`bin/86_flag_confidence.py`) should be **polarity-aware** — evaluate clade-membership on the +*derived* state in phylogenetic (ancestral/derived) space, not on raw reference-relative +variant calls. Then reverse-polarity backbone SNPs come out `monophyletic=true` and no +loader exemption is needed. + +## Issue 2 — `jointConfirmed=false` excludes ~27% of forward SNPs + +~80,400 forward-polarity defining SNPs (27% of all defining calls) are `jointConfirmed=false` +and excluded from counting. This is **kept** for now — it plausibly reflects real +"AEngine positive-only / not confirmed in joint genotyping" calls and may legitimately tame +the deep-backbone over-counting seen earlier. Flagging here only so the tree team is aware of +the magnitude; if ages come out systematically young, this gate is the first suspect. + +## Issue 3 — the INDEL artifact is unverified + +`chrY.ftdna.refined.indel.ingest.json` (ASR + indels + flags) has **not** been confirmed to +be generated correctly. Loading it surfaced the same confidence-flag behavior as the SNP +artifact (the indels themselves load fine: +10,576 indel variants), but until the indel +build is validated we are **staying on the SNP artifact** (`chrY.ftdna.refined.ingest.json`). +The SNP artifact should track the latest export format (it does — same flag schema). + +## Provenance + +- Deep-backbone SNP reality was spot-checked in the joint VCF: e.g. A1a (A-M31) stem SNPs + are GQ=99, normal DP, ~95% clean AC=5 synapomorphies — i.e. the data is real; the *flag* + was wrong, not the calls. +- A-M31 ages to formed ~119 kya / TMRCA ~12 kya, matching YFull — confirming the deep + spine is correctly calibrated once the backbone SNPs are allowed to count. diff --git a/rust/crates/du-db/src/denovo.rs b/rust/crates/du-db/src/denovo.rs index 1c5ade4b..2877a664 100644 --- a/rust/crates/du-db/src/denovo.rs +++ b/rust/crates/du-db/src/denovo.rs @@ -117,8 +117,19 @@ pub struct DenovoVariant { impl DenovoVariant { /// A defining call is low-confidence if it is not joint-confirmed (AEngine /// positive-only) or not monophyletic (homoplasic). Missing flags → confident. + /// + /// EXCEPTION: the `monophyletic` flag is unreliable for **reverse-polarity** SNPs + /// and must not gate them. The reference (CHM13/HG002) is haplogroup J, which sits + /// deep inside CT, so every backbone SNP above J (BT/CT/F/CF) carries the *derived* + /// allele as the reference. Their derived carriers are therefore reference-matching + /// (invisible to variant calling) while the ancestral A/B outgroup shows as the + /// variant — a paraphyletic set the monophyly test (computed from variant calls) + /// falsely flags non-monophyletic. Gating on it zeroes the entire deep backbone + /// (M168/M89/M42 themselves) and collapses the tree. See + /// documents/proposals/denovo-ingest-confidence-flags.md. fn low_confidence(&self) -> bool { - !self.joint_confirmed.unwrap_or(true) || !self.monophyletic.unwrap_or(true) + let reverse = self.polarity.as_deref() == Some("reverse"); + !self.joint_confirmed.unwrap_or(true) || (!reverse && !self.monophyletic.unwrap_or(true)) } } @@ -726,3 +737,44 @@ pub async fn list_conflicts( .await?; Ok(Page { items, total, page: page.max(1), page_size: limit }) } + +#[cfg(test)] +mod tests { + use super::*; + + fn var(polarity: &str, joint: Option, mono: Option) -> DenovoVariant { + DenovoVariant { + chrom: "chrY".into(), + pos: 1, + ref_: "A".into(), + alt: "C".into(), + ancestral: "A".into(), + derived: "C".into(), + reversion: false, + polarity: Some(polarity.into()), + joint_confirmed: joint, + monophyletic: mono, + } + } + + #[test] + fn reverse_polarity_is_exempt_from_monophyletic_gate() { + // The deep backbone (BT/CT/F) above the J reference is all reverse-polarity and + // falsely flagged monophyletic=false; it must still COUNT. + assert!( + !var("reverse", Some(true), Some(false)).low_confidence(), + "reverse-polarity monophyletic=false must NOT be low-confidence (CHM13=J inversion)" + ); + // Forward-polarity homoplasy is genuine → still excluded. + assert!( + var("forward", Some(true), Some(false)).low_confidence(), + "forward-polarity monophyletic=false stays low-confidence (real homoplasy)" + ); + // jointConfirmed still gates regardless of polarity. + assert!(var("reverse", Some(false), Some(true)).low_confidence(), "reverse !jointConfirmed"); + assert!(var("forward", Some(false), Some(true)).low_confidence(), "forward !jointConfirmed"); + // Clean calls (and missing flags → confident) count. + assert!(!var("forward", Some(true), Some(true)).low_confidence()); + assert!(!var("forward", None, None).low_confidence(), "missing flags → confident"); + } +} From 7945ed217bf0e2738e0887e476eab182cc6b52bf Mon Sep 17 00:00:00 2001 From: James Kane Date: Tue, 30 Jun 2026 07:13:20 -0500 Subject: [PATCH 06/13] =?UTF-8?q?fix(age):=20derive=20formed=5Fybp=20from?= =?UTF-8?q?=20the=20combined=20tmrca=20(formed=20=E2=89=A5=20tmrca)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit formed_ybp was written from the SNP-only propagation (age.formed.median), while tmrca_ybp comes from the COMBINED product (SNP × STR × genealogical) and is then raised by the causality back-correction. Two different sources: when the STR term or the causality lift pushed the combined tmrca above the SNP-only formed, the node ended up with formed < tmrca — an impossible ordering (a clade can't split from its parent more recently than its own members coalesce). Observed at R-L21 (formed 2795 < tmrca 3779) and R-DF13 (2667 < 3779). Now formed = combined_tmrca + the node's own branch time, off the SAME causality-corrected tmrca, so formed ≥ tmrca by construction. The branch offset (formed − tmrca) is taken from the SNP propagation — a node's branch length is a SNP-count property, independent of the STR/genealogical TMRCA evidence — and is ≥0 (0 for an age-transparent zero-SNP node). Tree-wide formed --- rust/crates/du-db/src/age.rs | 43 +++++++++++++++++++++--------------- 1 file changed, 25 insertions(+), 18 deletions(-) diff --git a/rust/crates/du-db/src/age.rs b/rust/crates/du-db/src/age.rs index 95e9bd1b..bfcafb5b 100644 --- a/rust/crates/du-db/src/age.rs +++ b/rust/crates/du-db/src/age.rs @@ -560,19 +560,26 @@ pub async fn recompute_combined_ages(pool: &PgPool) -> Result = HashMap::new(); - // (id, med, lo, hi, tester_count, formed) — written in phase 2. - let mut snp_writes: Vec<(i64, i32, i32, i32, i32, i32)> = Vec::new(); + // (id, med, lo, hi, tester_count) — written in phase 2. + let mut snp_writes: Vec<(i64, i32, i32, i32, i32)> = Vec::new(); + // The node's own branch time in years (`formed − tmrca`, ≥0 — convolution only adds + // age; 0 for an age-transparent zero-SNP node). `formed_ybp` is derived from the + // *combined* tmrca plus this offset (phase 2) so it stays ≥ tmrca even when the STR + // term or the causality back-correction raises the combined tmrca above the SNP-only + // formed. Keying the offset off the SNP propagation is correct: a node's branch time + // is a SNP-count property, independent of the STR/genealogical TMRCA evidence. + let mut branch_offset: HashMap = HashMap::new(); for (i, age) in ages.iter().enumerate() { let Some(age) = age else { continue }; let (med, lo, hi) = age.tmrca.ci95(); snp_pdf.insert(ids[i], age.tmrca.clone()); + branch_offset.insert(ids[i], (age.formed.median() - age.tmrca.median()).max(0.0).round() as i32); snp_writes.push(( ids[i], med.round() as i32, lo.round() as i32, hi.round() as i32, clades[i].tester_snps.len() as i32, - age.formed.median().round() as i32, )); } @@ -722,17 +729,9 @@ pub async fn recompute_combined_ages(pool: &PgPool) -> Result Result Date: Tue, 30 Jun 2026 15:34:38 -0500 Subject: [PATCH 07/13] docs(denovo): record INDEL extraction placement-collapsing issue MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The indel path is used to debug sample placements visually in the local web view; its known problem manifests there as placement collapsing. Note this is upstream in extraction, NOT the loader/age engine — with the reverse-polarity fix the indel artifact loads cleanly (backbone counts, +10,560 indel links), so it'll drop in once extraction is corrected. Until then we stay on the SNP artifact; load indel only to eyeball placements, revert before trusting ages. Co-Authored-By: Claude Opus 4.8 (1M context) --- .../denovo-ingest-confidence-flags.md | 22 ++++++++++++++----- 1 file changed, 17 insertions(+), 5 deletions(-) diff --git a/documents/proposals/denovo-ingest-confidence-flags.md b/documents/proposals/denovo-ingest-confidence-flags.md index 7effa603..a77b9f23 100644 --- a/documents/proposals/denovo-ingest-confidence-flags.md +++ b/documents/proposals/denovo-ingest-confidence-flags.md @@ -66,13 +66,25 @@ and excluded from counting. This is **kept** for now — it plausibly reflects r the deep-backbone over-counting seen earlier. Flagging here only so the tree team is aware of the magnitude; if ages come out systematically young, this gate is the first suspect. -## Issue 3 — the INDEL artifact is unverified +## Issue 3 — the INDEL artifact is unverified (extraction has placement-collapsing issues) `chrY.ftdna.refined.indel.ingest.json` (ASR + indels + flags) has **not** been confirmed to -be generated correctly. Loading it surfaced the same confidence-flag behavior as the SNP -artifact (the indels themselves load fine: +10,576 indel variants), but until the indel -build is validated we are **staying on the SNP artifact** (`chrY.ftdna.refined.ingest.json`). -The SNP artifact should track the latest export format (it does — same flag schema). +be generated correctly. The **indel extraction itself has a known problem**, observed +tree-side as **placement collapsing** — samples/branches collapsing in the topology (the +indel path is currently used to debug sample placements visually in the local web view, and +that is where the collapsing shows up). + +Important: this is an **upstream extraction** problem, **not** an AppView loader or age-engine +problem. Mechanically the indel artifact loads cleanly now — with the reverse-polarity +monophyletic fix (Issue 1) the deep backbone counts (CT 286/293, F 162/163, BT 8/8) and the ++10,560 indel variant-links load without error. So once the extraction is corrected the indel +tree should drop in without further loader work. + +Until the extraction is validated we are **staying on the SNP artifact** +(`chrY.ftdna.refined.ingest.json`). The SNP artifact tracks the latest export format (same +flag schema). Workflow when iterating on the indel build: load the indel artifact only to +eyeball placements in the local web view; revert to the SNP artifact (reload + branch-age +recompute) before trusting any ages. ## Provenance From 922b26fc2992e4ae30f16fb5f03352d6d8e00dc9 Mon Sep 17 00:00:00 2001 From: James Kane Date: Wed, 1 Jul 2026 00:30:18 -0500 Subject: [PATCH 08/13] fix(age): count only SNPs in the Poisson clock, exclude indels MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit McDonald 2021 §2.2.1: the SNP clock (Eq 2/3, µ_SNP≈8e-10 SNP·bp⁻¹·yr⁻¹) applies to SNPs; indels/complex variants may be used "provided a mutation rate for them can be accurately defined" — i.e. as a SEPARATE Eq-1 evidence term with its own µ, never folded into the SNP count m. build_clades counted every haplogroup_variant / biosample_private_variant row (SNPs AND indels) into one Poisson m at the SNP rate, which on an indel-bearing tree biases every indel-dense branch older. Fix: both the branch and tester count queries filter mutation_type = 'SNP' (mutation_type cleanly separates SNP from INS/DEL/INDEL/MNP/STR). A node whose defining variants are all indels then has 0 age-countable SNPs and is handled as an age-transparent zero-length branch (existing branch_time δ(t=0) rule) — exactly like the palindromic/empty backbone nodes. Result: recompute over the indel tree now yields ages IDENTICAL to the SNP-only tree (BT 95712, CT 63255, F 47833, A-M31 11584, P312 4405 — matching to the year) while retaining the indel-resolved topology (11409 nodes vs 11363; +46 branches that break SNP-only polytomies). 0 formed --- rust/crates/du-db/src/age.rs | 40 ++++++++++++++++++++++++++++++++---- 1 file changed, 36 insertions(+), 4 deletions(-) diff --git a/rust/crates/du-db/src/age.rs b/rust/crates/du-db/src/age.rs index bfcafb5b..9d46da56 100644 --- a/rust/crates/du-db/src/age.rs +++ b/rust/crates/du-db/src/age.rs @@ -385,12 +385,23 @@ async fn build_clades(pool: &PgPool) -> Result<(Vec, Vec), DbError> "" }; + // Count ONLY single-nucleotide SNPs into the Poisson clock. The clock rate + // `SNP_RATE` (µ_SNP) is defined for SNPs; McDonald 2021 §2.2.1 permits indels/ + // complex variants only "provided a mutation rate for them can be accurately + // defined" — i.e. as a *separate* evidence term (Eq 1) with its own µ, never + // folded into the SNP `m`. So exclude INDEL/DEL/INS/MNP/STR here. A node whose + // defining variants are ALL indels then has 0 age-countable SNPs and is handled + // as an age-transparent zero-length branch (see `branch_time`), exactly like the + // palindromic/empty backbone nodes. (A dedicated indel clock — µ_indel + its own + // callable denominator — is future work: documents/proposals/…-confidence-flags.md.) + let snp_only = "AND v.mutation_type = 'SNP'::core.mutation_type"; + // Branch defining-SNP counts (recurrent-region-masked, callable-mask-intersected). let mask = recurrent_mask(); let branch: Vec<(i64, i64)> = sqlx::query_as(&format!( "SELECT hv.haplogroup_id, count(*)::bigint FROM tree.haplogroup_variant hv \ JOIN core.variant v ON v.id=hv.variant_id \ - WHERE hv.valid_until IS NULL AND NOT hv.low_confidence AND {mask} {callable_filter} GROUP BY hv.haplogroup_id" + WHERE hv.valid_until IS NULL AND NOT hv.low_confidence {snp_only} AND {mask} {callable_filter} GROUP BY hv.haplogroup_id" )) .fetch_all(pool) .await?; @@ -411,7 +422,7 @@ async fn build_clades(pool: &PgPool) -> Result<(Vec, Vec), DbError> LEFT JOIN genomics.biosample_callable_loci cl \ ON cl.sample_guid=pv.sample_guid AND cl.chromosome IN ('chrY','Y') \ WHERE pv.status='ACTIVE' AND pv.haplogroup_type='Y_DNA'::core.dna_type \ - AND pv.terminal_haplogroup_id IS NOT NULL AND {mask} {callable_filter} \ + AND pv.terminal_haplogroup_id IS NOT NULL {snp_only} AND {mask} {callable_filter} \ GROUP BY pv.terminal_haplogroup_id, pv.sample_guid" )) .fetch_all(pool) @@ -986,6 +997,18 @@ mod tests { .unwrap() } + /// An indel variant — must be EXCLUDED from the SNP Poisson clock (McDonald §2.2.1). + async fn ins_indel(pool: &PgPool, name: &str) -> i64 { + sqlx::query_scalar( + "INSERT INTO core.variant (canonical_name, mutation_type, naming_status, annotations) \ + VALUES ($1, 'INS'::core.mutation_type, 'NAMED'::core.naming_status, '{}'::jsonb) RETURNING id", + ) + .bind(name) + .fetch_one(pool) + .await + .unwrap() + } + /// Seed a 3-node chain with one tester, run the whole pipeline, and check the /// het-mask, causality (parent older), and formed > tmrca — against real PG. #[tokio::test] @@ -1019,6 +1042,12 @@ mod tests { } let hetdef = ins_var(&pool, "LEAFDEFHET", true).await; sqlx::query("INSERT INTO tree.haplogroup_variant (haplogroup_id, variant_id) VALUES ($1,$2)").bind(leaf).bind(hetdef).execute(&pool).await.unwrap(); + // Two INDEL defining variants on the leaf branch — must NOT be counted into the + // SNP clock (they need their own rate; McDonald §2.2.1). branch_snps stays 3. + for i in 0..2 { + let v = ins_indel(&pool, &format!("LEAFINDEL{i}")).await; + sqlx::query("INSERT INTO tree.haplogroup_variant (haplogroup_id, variant_id) VALUES ($1,$2)").bind(leaf).bind(v).execute(&pool).await.unwrap(); + } // One tester under leaf: 12.5 Mbp callable, 5 private SNPs + 1 het (masked). sqlx::query("INSERT INTO core.biosample (sample_guid, source) VALUES ($1::uuid, 'CITIZEN')").bind(GUID).execute(&pool).await.unwrap(); @@ -1029,12 +1058,15 @@ mod tests { } let hv = ins_var(&pool, "PRIVHET", true).await; sqlx::query("INSERT INTO tree.biosample_private_variant (sample_guid, variant_id, haplogroup_type, terminal_haplogroup_id) VALUES ($1::uuid,$2,'Y_DNA'::core.dna_type,$3)").bind(GUID).bind(hv).bind(leaf).execute(&pool).await.unwrap(); + // An INDEL private variant — excluded from the tester's SNP count (stays 5). + let piv = ins_indel(&pool, "PRIVINDEL").await; + sqlx::query("INSERT INTO tree.biosample_private_variant (sample_guid, variant_id, haplogroup_type, terminal_haplogroup_id) VALUES ($1::uuid,$2,'Y_DNA'::core.dna_type,$3)").bind(GUID).bind(piv).bind(leaf).execute(&pool).await.unwrap(); // (a) build_clades: het-masking + structure. let (clades, ids) = build_clades(&pool).await.unwrap(); let at = |id: i64| ids.iter().position(|&x| x == id).unwrap(); - assert_eq!(clades[at(leaf)].tester_snps, vec![(5, 12_500_000.0)], "het private SNP masked → 5 counted, with its callable bp"); - assert_eq!(clades[at(leaf)].branch_snps, 3, "het defining SNP masked → 3"); + assert_eq!(clades[at(leaf)].tester_snps, vec![(5, 12_500_000.0)], "het private SNP masked + INDEL private excluded → 5 counted, with its callable bp"); + assert_eq!(clades[at(leaf)].branch_snps, 3, "het defining SNP masked + 2 INDEL defining excluded → 3"); assert!(clades[at(mid)].children.contains(&at(leaf))); assert!(clades[at(root)].children.contains(&at(mid))); assert!((clades[at(leaf)].callable_bp - 12_500_000.0).abs() < 1.0); From 2cfe1771243e172ffe7f32b19a4ce6f6365eec1f Mon Sep 17 00:00:00 2001 From: James Kane Date: Wed, 1 Jul 2026 01:03:39 -0500 Subject: [PATCH 09/13] docs(denovo): note engine now excludes indels from the SNP clock MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Separates the two concerns for the tree team: the age engine handles indels correctly (SNP-only clock, indel-only nodes age-transparent → indel-tree ages == SNP-tree ages), so any remaining indel-view oddities are the upstream extraction/placement problem, not the clock. Co-Authored-By: Claude Opus 4.8 (1M context) --- .../proposals/denovo-ingest-confidence-flags.md | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/documents/proposals/denovo-ingest-confidence-flags.md b/documents/proposals/denovo-ingest-confidence-flags.md index a77b9f23..d16b108c 100644 --- a/documents/proposals/denovo-ingest-confidence-flags.md +++ b/documents/proposals/denovo-ingest-confidence-flags.md @@ -81,10 +81,19 @@ monophyletic fix (Issue 1) the deep backbone counts (CT 286/293, F 162/163, BT 8 tree should drop in without further loader work. Until the extraction is validated we are **staying on the SNP artifact** -(`chrY.ftdna.refined.ingest.json`). The SNP artifact tracks the latest export format (same -flag schema). Workflow when iterating on the indel build: load the indel artifact only to -eyeball placements in the local web view; revert to the SNP artifact (reload + branch-age -recompute) before trusting any ages. +(`chrY.ftdna.refined.ingest.json`) for published ages. The SNP artifact tracks the latest +export format (same flag schema). + +**Engine vs extraction are now cleanly separated.** The age engine no longer counts indels +into the SNP clock (`build_clades` filters `mutation_type = 'SNP'`; McDonald §2.2.1 — indels +need their own rate, as a separate Eq-1 term, never folded into the SNP `m`). Nodes defined +only by indels get 0 age-countable SNPs and go age-transparent (zero-length branch), like the +palindromic/empty backbone nodes. Consequence: a branch-age recompute on the **indel tree now +yields ages identical to the SNP tree** (BT 95712, CT 63255, A-M31 11584, P312 4405 — matching +to the year), while retaining the indel-resolved topology (+46 branches that break SNP-only +polytomies). So the indel artifact can be loaded and aged safely; any remaining oddities in the +indel view are the **extraction/placement** problem above, not the clock. A dedicated indel +clock (µ_indel + its own callable denominator) remains future work. ## Provenance From 92092668224dbfc2efdc01b2c45eaa8484923822 Mon Sep 17 00:00:00 2001 From: James Kane Date: Wed, 1 Jul 2026 06:52:26 -0500 Subject: [PATCH 10/13] fix(web): clear inverted-block reverse-complement display artifacts MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two Y-tree sidebar symptoms, both rooted in the same inverted-block liftover leaving affected catalog SNPs stored reverse-complemented (CHM13/HG002 Y is haplogroup J, deep inside CT — inverted blocks land RC): - Y18975 spurious "back-mutation" badge. The web reconstructs back-mutation from tree link.derived == catalog ancestral, but Y18975 is single-origin (recurrent=false); its link T>A vs catalog A>T is pure strand artifact. Gate back_mutation on `recurrent` — a non-recurrent SNP has nothing to revert from. Clears 1009 non-recurrent false positives incl. Y18975. - A9005 shown as "chrY:21905878C>A". The de-novo loader minted a coordinate name because the catalog's A9005 (G>T) is the RC of the tree link (C>A), so the exact-allele join at load missed it. Recover the real name for display via an RC-aware lateral: match a real-named catalog SNP at the same hs1 position whose allele set equals the link's OR its reverse- complement. GIN containment on position keeps it cheap; runs only for coordinate-named/unnamed links. No tree rebuild needed. See documents/proposals/denovo-ingest-confidence-flags.md. Co-Authored-By: Claude Opus 4.8 (1M context) --- rust/crates/du-db/src/haplogroup.rs | 25 +++++++++++++++++++++++-- rust/crates/du-web/src/routes/tree.rs | 11 ++++++++--- 2 files changed, 31 insertions(+), 5 deletions(-) diff --git a/rust/crates/du-db/src/haplogroup.rs b/rust/crates/du-db/src/haplogroup.rs index c089cddd..a87319c3 100644 --- a/rust/crates/du-db/src/haplogroup.rs +++ b/rust/crates/du-db/src/haplogroup.rs @@ -704,8 +704,16 @@ pub async fn variants_of( link_derived: Option, recurrent: bool, } + // A defining SNP the de-novo loader minted a coordinate name for (`chrY:pos…>…`) + // because the inverted-block liftover left the catalog's copy reverse-complemented, + // so the exact-allele match at load missed the real named entry. Recover the real + // name for display by finding a real-named catalog SNP at the same hs1 position whose + // alleles equal this link's set OR its reverse-complement. (GIN-indexed containment + // on the position keeps the per-row lookup cheap; it runs only for coordinate-named / + // unnamed links.) See documents/proposals/denovo-ingest-confidence-flags.md. let rows: Vec = sqlx::query_as( - "SELECT v.canonical_name, v.mutation_type::text AS mutation_type, v.aliases, v.coordinates, \ + "SELECT COALESCE(named.nm, v.canonical_name) AS canonical_name, \ + v.mutation_type::text AS mutation_type, v.aliases, v.coordinates, \ hv.ancestral_allele AS link_ancestral, hv.derived_allele AS link_derived, \ EXISTS (SELECT 1 FROM tree.haplogroup_variant hv2 \ WHERE hv2.variant_id = v.id AND hv2.valid_until IS NULL \ @@ -713,8 +721,21 @@ pub async fn variants_of( FROM core.variant v \ JOIN tree.haplogroup_variant hv ON hv.variant_id = v.id AND hv.valid_until IS NULL \ JOIN tree.haplogroup h ON h.id = hv.haplogroup_id \ + LEFT JOIN LATERAL ( \ + SELECT v2.canonical_name AS nm FROM core.variant v2 \ + WHERE v2.id <> v.id AND v2.mutation_type = v.mutation_type \ + AND v2.canonical_name IS NOT NULL AND v2.canonical_name !~ '^chr' \ + AND v2.coordinates @> jsonb_build_object('hs1', jsonb_build_object('position', v.coordinates->'hs1'->'position')) \ + AND least(v2.coordinates->'hs1'->>'ancestral', v2.coordinates->'hs1'->>'derived') \ + || greatest(v2.coordinates->'hs1'->>'ancestral', v2.coordinates->'hs1'->>'derived') IN ( \ + least(v.coordinates->'hs1'->>'ancestral', v.coordinates->'hs1'->>'derived') \ + || greatest(v.coordinates->'hs1'->>'ancestral', v.coordinates->'hs1'->>'derived'), \ + least(translate(v.coordinates->'hs1'->>'ancestral','ACGT','TGCA'), translate(v.coordinates->'hs1'->>'derived','ACGT','TGCA')) \ + || greatest(translate(v.coordinates->'hs1'->>'ancestral','ACGT','TGCA'), translate(v.coordinates->'hs1'->>'derived','ACGT','TGCA')) ) \ + ORDER BY core.ysnp_name_rank(v2.canonical_name), v2.canonical_name LIMIT 1 \ + ) named ON v.canonical_name IS NULL OR v.canonical_name ~ '^chr[0-9A-Za-z]*:' \ WHERE h.name = $1 AND h.haplogroup_type::text = $2 AND h.valid_until IS NULL \ - ORDER BY v.canonical_name", + ORDER BY COALESCE(named.nm, v.canonical_name)", ) .bind(name) .bind(pg_enum_label(&dna_type)?) diff --git a/rust/crates/du-web/src/routes/tree.rs b/rust/crates/du-web/src/routes/tree.rs index 41969bee..33b09f05 100644 --- a/rust/crates/du-web/src/routes/tree.rs +++ b/rust/crates/du-web/src/routes/tree.rs @@ -335,9 +335,14 @@ async fn snp_sidebar( (Some(a), Some(d)) => Some(format!("{a}>{d}")), _ => None, }; - // Back-mutation: this branch's derived state is the SNP's ancestral allele. - let back_mutation = - matches!((&v.link_derived, coord_ancestral(&v.coordinates)), (Some(d), Some(a)) if *d == a); + // Back-mutation: this branch's derived state is the SNP's ancestral allele — + // but only meaningful for a RECURRENT SNP (one that occurs on another branch, + // so there is something to revert). A single-origin SNP has nothing to back- + // mutate from, so a derived==ancestral coincidence there is a catalog strand + // artifact (e.g. Y18975, stored reverse-complemented by the inverted-block + // liftover: link T>A vs catalog A>T), not biology. Gate on `recurrent`. + let back_mutation = v.recurrent + && matches!((&v.link_derived, coord_ancestral(&v.coordinates)), (Some(d), Some(a)) if *d == a); let aliases = json_str_list(&v.aliases); // UNNAMED variants (homoplasy collisions) fall back to an alias. let name = v From 2378d2e503ec5c8c08f0b1b62b76a0a4bf6b6208 Mon Sep 17 00:00:00 2001 From: James Kane Date: Wed, 1 Jul 2026 09:15:34 -0500 Subject: [PATCH 11/13] fix(tree): protect de-novo topology tips from the call-based placement prune MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `tree_sample::recompute_placements` places non-D2C samples under the node their published paper call resolves to, then prunes every haplogroup_sample row not in its candidate set. But the de-novo loader (`denovo::load`) places tree tips by ML-tree *topology*, and those cohort samples (1000G/HGDP/PRJEB) carry no Y/mt call under the keys `pick_original_call` reads — so a single recompute, or the daily `tree-samples-recompute` cron, silently wiped all ~9.6k topology tips (observed: 0 Y_DNA rows after a load+recompute). Protect de-novo-origin biosamples (loader stamps `source_attrs->>'denovo'='true'` at creation) from both the prune and any overwrite, exactly like CURATED rows. No new status value — keeps every existing rollup/consumer (ystr, dedup, cladogram) unchanged. Regression test: a denovo-flagged PLACED tip with no call survives recompute and still counts toward its node rollup. Co-Authored-By: Claude Opus 4.8 (1M context) --- rust/crates/du-db/src/tree_sample.rs | 21 +++++++++- rust/crates/du-db/tests/tree_sample.rs | 54 ++++++++++++++++++++++++++ 2 files changed, 73 insertions(+), 2 deletions(-) diff --git a/rust/crates/du-db/src/tree_sample.rs b/rust/crates/du-db/src/tree_sample.rs index 8796f771..cafeff5f 100644 --- a/rust/crates/du-db/src/tree_sample.rs +++ b/rust/crates/du-db/src/tree_sample.rs @@ -74,14 +74,31 @@ async fn recompute_placements_locked(pool: &PgPool, dna_type: DnaType) -> Result .into_iter() .collect(); + // De-novo tree tips are placed by the loader (`denovo::load`) by tree *topology*, + // not by a published call — the sample sits under the node its lineage resolves to + // in the ML tree. This call-based recompute has no way to reproduce that placement + // (these cohort samples carry no Y/mt call under the keys `pick_original_call` reads), + // so it must not touch them: protect de-novo-origin biosamples from both the prune + // and any overwrite, exactly like CURATED. Without this, a single recompute (or the + // daily cron) silently wipes every topology tip. The loader stamps these biosamples + // with `source_attrs->>'denovo' = 'true'` at creation. + let denovo: std::collections::HashSet = sqlx::query_scalar::<_, Uuid>( + "SELECT sample_guid FROM core.biosample WHERE source_attrs->>'denovo' = 'true'", + ) + .fetch_all(pool) + .await? + .into_iter() + .collect(); + // Resolve each sample's call to a node id, caching by call text (calls repeat heavily). let mut cache: HashMap> = HashMap::new(); let mut processed: Vec = Vec::new(); let mut placements: Vec<(Uuid, Option, String)> = Vec::new(); let mut report = PlacementReport::default(); for (guid, arr) in rows { - // Preserve curator placements: keep the row (protect from prune), don't re-resolve. - if curated.contains(&guid) { + // Preserve curator + de-novo-topology placements: keep the row (protect from + // prune), don't re-resolve or overwrite. + if curated.contains(&guid) || denovo.contains(&guid) { processed.push(guid); continue; } diff --git a/rust/crates/du-db/tests/tree_sample.rs b/rust/crates/du-db/tests/tree_sample.rs index f2baf56b..57c45a50 100644 --- a/rust/crates/du-db/tests/tree_sample.rs +++ b/rust/crates/du-db/tests/tree_sample.rs @@ -164,3 +164,57 @@ async fn places_non_d2c_samples_and_records_unplaced() { .unwrap(); assert_eq!(curated, "CURATED", "curator placement preserved across recompute"); } + +/// A de-novo tree tip is placed by the loader by topology (no published call). The +/// call-based recompute must NOT prune or overwrite it — the biosample's +/// `source_attrs->>'denovo' = 'true'` marker protects it, like CURATED. Regression for +/// the "tree-samples-recompute wipes every de-novo tip" bug. +#[tokio::test] +async fn denovo_topology_tips_survive_recompute() { + let Some(url) = database_url() else { + eprintln!("DATABASE_URL unset — skipping denovo-tip test"); + return; + }; + let db = du_db::testing::ephemeral_db(&url).await.expect("ephemeral db"); + let pool = db.pool().clone(); + + let node_id = node(&pool, "R-M269", &[]).await; + + // A de-novo cohort tip: marked denovo, carries NO published Y call (empty + // original_haplogroups), and is placed at a node directly by the loader. + let guid: Uuid = sqlx::query_scalar( + "INSERT INTO core.biosample (source, accession, original_haplogroups, source_attrs) \ + VALUES ('STANDARD'::core.biosample_source, 'DENOVO-TIP', '[]'::jsonb, \ + jsonb_build_object('denovo', true)) RETURNING sample_guid", + ) + .fetch_one(&pool) + .await + .expect("insert denovo biosample"); + sqlx::query( + "INSERT INTO tree.haplogroup_sample (sample_guid, dna_type, haplogroup_id, call_text, status, refreshed_at) \ + VALUES ($1, 'Y_DNA'::core.dna_type, $2, 'R-M269', 'PLACED', now())", + ) + .bind(guid) + .bind(node_id) + .execute(&pool) + .await + .expect("loader-style tip placement"); + + // Recompute (which finds no call for this sample) must leave the tip intact. + tree_sample::recompute_placements(&pool, DnaType::YDna).await.unwrap(); + + let (status, hid): (String, Option) = + sqlx::query_as("SELECT status, haplogroup_id FROM tree.haplogroup_sample WHERE sample_guid = $1") + .bind(guid) + .fetch_optional(&pool) + .await + .unwrap() + .expect("de-novo tip must survive recompute, not be pruned"); + assert_eq!(status, "PLACED", "de-novo placement untouched"); + assert_eq!(hid, Some(node_id), "still placed at its topology node"); + assert_eq!( + tree_sample::counts_by_node(&pool, DnaType::YDna).await.unwrap().get(&node_id).copied(), + Some(1), + "de-novo tip still counts toward the node rollup", + ); +} From 783709672888cfcac9b31a63fa8ec7abe35765c5 Mon Sep 17 00:00:00 2001 From: James Kane Date: Wed, 1 Jul 2026 12:28:46 -0500 Subject: [PATCH 12/13] fix(denovo): collapse self-labeled private singletons, add --keep-private MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The de-novo pipeline never leaves a private singleton's `label` null — it names the branch after the single sample's own id. The old collapse predicate only matched `label.is_none()`, so it caught none of them: 6,549 UUID-named branches leaked into the tree as public nodes, each with a tip of the same id hanging off it (e.g. `7d7e9716…` under R-FT88981). Widen "no public name" to also mean `label == that tip's sample id`. On the real ftdna-indel ingest this collapses all 8,226 self-labeled singletons onto their public parent (seeding their SNPs into the discovery substrate) while preserving the 437 genuinely SNP-named single-sample terminals (R-FT49699, R-BY95127, …) as real public branches. Add a preprod-only `--keep-private` flag (threaded through load_denovo → denovo::load) that disables the collapse so private branches render as nodes for visual placement debugging. Production omits it; a forgotten flag is safe (default collapses). Co-Authored-By: Claude Opus 4.8 (1M context) --- rust/crates/du-db/src/denovo.rs | 50 ++++++++++++++------- rust/crates/du-db/tests/denovo.rs | 10 ++--- rust/crates/du-migrate/src/bin/tree_init.rs | 17 ++++--- 3 files changed, 50 insertions(+), 27 deletions(-) diff --git a/rust/crates/du-db/src/denovo.rs b/rust/crates/du-db/src/denovo.rs index 2877a664..ebfcce34 100644 --- a/rust/crates/du-db/src/denovo.rs +++ b/rust/crates/du-db/src/denovo.rs @@ -168,7 +168,7 @@ fn confidence(support: Option) -> &'static str { /// Load a de-novo tree document as the tree foundation. Assumes `tree.*` is /// already cleared (greenfield). Commits the topology, then recomputes backbone /// and bumps the tree revision. -pub async fn load(pool: &PgPool, doc: &DenovoTree) -> Result { +pub async fn load(pool: &PgPool, doc: &DenovoTree, keep_private: bool) -> Result { let dna: DnaType = match doc.haplogroup_type.as_str() { "Y_DNA" => DnaType::YDna, "MT_DNA" => DnaType::MtDna, @@ -210,12 +210,22 @@ pub async fn load(pool: &PgPool, doc: &DenovoTree) -> Result = HashMap::new(); let mut parent_of: HashMap<&str, &str> = HashMap::new(); for node in &doc.nodes { @@ -225,19 +235,25 @@ pub async fn load(pool: &PgPool, doc: &DenovoTree) -> Result = HashMap::new(); + let mut tip_sample: HashMap<&str, &str> = HashMap::new(); for tip in &doc.tips { *tip_count.entry(tip.parent_node.as_str()).or_default() += 1; + tip_sample.entry(tip.parent_node.as_str()).or_insert(tip.sample.as_str()); } - let private: HashSet<&str> = doc - .nodes - .iter() - .filter(|n| { - n.label.is_none() - && child_count.get(n.id.as_str()).copied().unwrap_or(0) == 0 - && tip_count.get(n.id.as_str()).copied().unwrap_or(0) == 1 - }) - .map(|n| n.id.as_str()) - .collect(); + let private: HashSet<&str> = if keep_private { + HashSet::new() + } else { + doc.nodes + .iter() + .filter(|n| { + child_count.get(n.id.as_str()).copied().unwrap_or(0) == 0 + && tip_count.get(n.id.as_str()).copied().unwrap_or(0) == 1 + && (n.label.is_none() + || n.label.as_deref() == tip_sample.get(n.id.as_str()).copied()) + }) + .map(|n| n.id.as_str()) + .collect() + }; rep.private_collapsed = private.len(); for node in &doc.nodes { diff --git a/rust/crates/du-db/tests/denovo.rs b/rust/crates/du-db/tests/denovo.rs index 52890abf..7e04253f 100644 --- a/rust/crates/du-db/tests/denovo.rs +++ b/rust/crates/du-db/tests/denovo.rs @@ -118,8 +118,8 @@ async fn denovo_lineages_coexist_and_clear_independently() { let db = du_db::testing::ephemeral_db(&url).await.expect("ephemeral db"); let pool = db.pool().clone(); - denovo::load(&pool, &doc()).await.expect("load Y"); // Y_DNA: 4 nodes, Node4 collapses → 3 - denovo::load(&pool, &mt_doc()).await.expect("load mt"); // MT_DNA, 2 nodes + denovo::load(&pool, &doc(), false).await.expect("load Y"); // Y_DNA: 4 nodes, Node4 collapses → 3 + denovo::load(&pool, &mt_doc(), false).await.expect("load mt"); // MT_DNA, 2 nodes assert_eq!(count_dna(&pool, "Y_DNA").await, 3); assert_eq!(count_dna(&pool, "MT_DNA").await, 2); @@ -145,7 +145,7 @@ async fn clear_dna_neutralizes_private_variant_refs_for_repeatable_reload() { seed_catalog_snp(&pool, "M269", 21452686, "T", "C").await; seed_catalog_snp(&pool, "L21", 13500000, "G", "A").await; - denovo::load(&pool, &doc()).await.expect("load Y"); + denovo::load(&pool, &doc(), false).await.expect("load Y"); // Seed one DENOVO + one FED private-variant row referencing a loaded node. let sample: uuid::Uuid = @@ -184,7 +184,7 @@ async fn clear_dna_neutralizes_private_variant_refs_for_repeatable_reload() { assert_eq!(fed_terminal, None, "FED row's stale node link nulled"); // The actual goal: a fresh load now succeeds end-to-end over the cleared lineage. - denovo::load(&pool, &doc()).await.expect("reload Y"); + denovo::load(&pool, &doc(), false).await.expect("reload Y"); assert_eq!(count_dna(&pool, "Y_DNA").await, 3, "reloaded cleanly"); } @@ -201,7 +201,7 @@ async fn loads_denovo_tree_with_catalog_reuse_and_mint() { seed_catalog_snp(&pool, "M269", 21452686, "T", "C").await; seed_catalog_snp(&pool, "L21", 13500000, "G", "A").await; - let rep = denovo::load(&pool, &doc()).await.expect("load"); + let rep = denovo::load(&pool, &doc(), false).await.expect("load"); // Node4 is an unlabeled single-tip leaf → collapsed as a private singleton: it is // not published as a node, its tip hangs on Node4's public parent (Node3), and its diff --git a/rust/crates/du-migrate/src/bin/tree_init.rs b/rust/crates/du-migrate/src/bin/tree_init.rs index 4db9090c..876861a9 100644 --- a/rust/crates/du-migrate/src/bin/tree_init.rs +++ b/rust/crates/du-migrate/src/bin/tree_init.rs @@ -38,6 +38,13 @@ struct Args { /// Apply the load (mutates the tree / mask; required). #[arg(long)] apply: bool, + /// Keep private singleton branches as visible nodes instead of collapsing them. + /// For PREPROD only (visual placement debugging). Production omits this: private + /// terminals (a leaf with one sample and no public SNP name) collapse onto their + /// public parent, and their SNPs are seeded into the discovery substrate to prime + /// branch discovery. Applies to `--denovo-y` / `--denovo-mt`. + #[arg(long)] + keep_private: bool, } /// Parse a BED (`chrom\tstart\tend`, half-open) into `(start, end)` spans, keeping @@ -59,14 +66,14 @@ fn parse_bed(path: &str) -> anyhow::Result> { /// Clear one lineage and load its de-novo ingest JSON. `expect` is the document's /// declared `haplogroup_type` (`Y_DNA`/`MT_DNA`) — a mismatch aborts before any write. -async fn load_denovo(pool: &PgPool, path: &str, dna: DnaType, expect: &str, apply: bool) -> anyhow::Result<()> { +async fn load_denovo(pool: &PgPool, path: &str, dna: DnaType, expect: &str, apply: bool, keep_private: bool) -> anyhow::Result<()> { anyhow::ensure!(apply, "--denovo-* mutates the tree; pass --apply"); let doc: du_db::denovo::DenovoTree = serde_json::from_str(&std::fs::read_to_string(path)?)?; anyhow::ensure!(doc.haplogroup_type == expect, "expected a {expect} document, got {}", doc.haplogroup_type); - tracing::info!(%path, hgtype = expect, nodes = doc.nodes.len(), tips = doc.tips.len(), root = %doc.root, "de-novo: loading foundation"); + tracing::info!(%path, hgtype = expect, nodes = doc.nodes.len(), tips = doc.tips.len(), root = %doc.root, keep_private, "de-novo: loading foundation"); let cleared = du_db::haplogroup::clear_dna(pool, dna).await?; tracing::info!(cleared_haplogroups = cleared, hgtype = expect, "de-novo: cleared lineage"); - let rep = du_db::denovo::load(pool, &doc).await?; + let rep = du_db::denovo::load(pool, &doc, keep_private).await?; tracing::info!( hgtype = expect, nodes = rep.nodes, edges = rep.edges, variant_links = rep.variant_links, variants_reused = rep.variants_reused, variants_created = rep.variants_created, @@ -95,8 +102,8 @@ async fn main() -> anyhow::Result<()> { du_db::run_migrations(&pool).await?; match (&args.denovo_y, &args.denovo_mt) { - (Some(path), _) => load_denovo(&pool, path, DnaType::YDna, "Y_DNA", args.apply).await?, - (None, Some(path)) => load_denovo(&pool, path, DnaType::MtDna, "MT_DNA", args.apply).await?, + (Some(path), _) => load_denovo(&pool, path, DnaType::YDna, "Y_DNA", args.apply, args.keep_private).await?, + (None, Some(path)) => load_denovo(&pool, path, DnaType::MtDna, "MT_DNA", args.apply, args.keep_private).await?, (None, None) if args.callable_mask.is_none() => { anyhow::bail!("pass --denovo-y , --denovo-mt , or --callable-mask ") } From 9f253ef8f339611b4c52534a09b640561061ef28 Mon Sep 17 00:00:00 2001 From: James Kane Date: Wed, 1 Jul 2026 15:14:31 -0500 Subject: [PATCH 13/13] feat(tree): show placed-sample count in sidebar instead of the list MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The SNP sidebar's "Placed samples" section listed every placed non-D2C sample leaf at or below the node (capped at 50 with a "+N more" note). Replace the list with just the count. - du-db: add `tree_sample::count_under` — same subtree CTE + PLACED/CURATED filters as `samples_under`, but `count(*)` with no row materialization or publication join. - du-web: `SnpSidebar` drops `samples`/`samples_more` (and the now-unused `LeafRow`/`SIDEBAR_SAMPLE_CAP`) for a single `sample_count`. - template: render the bold count under the existing "Placed samples" header. - locales: drop the orphaned `tree.samples.more` key (en/es/fr). Co-Authored-By: Claude Opus 4.8 (1M context) --- rust/crates/du-db/src/tree_sample.rs | 23 ++++++++++++ rust/crates/du-web/src/routes/tree.rs | 36 +++---------------- .../du-web/templates/tree/snp_sidebar.html | 19 ++-------- rust/locales/en.txt | 1 - rust/locales/es.txt | 1 - rust/locales/fr.txt | 1 - 6 files changed, 31 insertions(+), 50 deletions(-) diff --git a/rust/crates/du-db/src/tree_sample.rs b/rust/crates/du-db/src/tree_sample.rs index cafeff5f..ef6ad06d 100644 --- a/rust/crates/du-db/src/tree_sample.rs +++ b/rust/crates/du-db/src/tree_sample.rs @@ -304,3 +304,26 @@ pub async fn samples_under(pool: &PgPool, node_name: &str, dna_type: DnaType) -> .fetch_all(pool) .await?) } + +/// Count of placed samples at-or-below a node — the same subtree/filters as [`samples_under`] +/// but without materializing the rows or their publications (the sidebar shows just the count). +pub async fn count_under(pool: &PgPool, node_name: &str, dna_type: DnaType) -> Result { + Ok(sqlx::query_scalar( + "WITH RECURSIVE sub AS ( \ + SELECT id FROM tree.haplogroup \ + WHERE name = $1 AND haplogroup_type::text = $2 AND valid_until IS NULL \ + UNION ALL \ + SELECT r.child_haplogroup_id FROM tree.haplogroup_relationship r \ + JOIN sub ON r.parent_haplogroup_id = sub.id \ + WHERE r.valid_until IS NULL \ + ) \ + SELECT count(*)::bigint FROM tree.haplogroup_sample hs \ + JOIN sub ON sub.id = hs.haplogroup_id \ + JOIN core.biosample b ON b.sample_guid = hs.sample_guid \ + WHERE hs.dna_type::text = $2 AND hs.status IN ('PLACED','CURATED') AND b.deleted = false", + ) + .bind(node_name) + .bind(pg_enum_label(&dna_type)?) + .fetch_one(pool) + .await?) +} diff --git a/rust/crates/du-web/src/routes/tree.rs b/rust/crates/du-web/src/routes/tree.rs index 33b09f05..34a28490 100644 --- a/rust/crates/du-web/src/routes/tree.rs +++ b/rust/crates/du-web/src/routes/tree.rs @@ -154,10 +154,8 @@ struct SnpSidebar { /// Total markers in the node's reconstructed ancestral haplotype (context for /// the change count; 0 ⇒ no STR evidence in the subtree). str_motif_markers: usize, - /// Placed non-D2C sample leaves at or below this node (capped for the sidebar). - samples: Vec, - /// How many more placed samples exist beyond the shown `samples` (0 ⇒ all shown). - samples_more: i64, + /// Count of placed non-D2C sample leaves at or below this node (shown instead of the list). + sample_count: i64, } /// One Y-STR mutation along the branch leading into this node, e.g. DYS393 12→13. @@ -167,16 +165,6 @@ struct StrChangeRow { to: i32, } -/// One placed sample row in the sidebar (label + optional paper citation). -struct LeafRow { - label: String, - source: String, - citation: Option, -} - -/// Max leaf rows rendered in the sidebar before collapsing to an "+N more" note. -const SIDEBAR_SAMPLE_CAP: usize = 50; - struct VariantRow { name: String, mutation_type: String, @@ -419,21 +407,8 @@ async fn snp_sidebar( (Vec::new(), 0) }; - // Placed non-D2C sample leaves at or below this node (capped for the sidebar). - let mut leaves = du_db::tree_sample::samples_under(&st.pool, &name, dna_type).await?; - let samples_more = (leaves.len() as i64 - SIDEBAR_SAMPLE_CAP as i64).max(0); - leaves.truncate(SIDEBAR_SAMPLE_CAP); - let samples: Vec = leaves - .into_iter() - .map(|s| { - let label = s.accession.or(s.alias).unwrap_or_else(|| s.sample_guid.to_string()); - let citation = s.pub_title.map(|t| match s.pub_doi { - Some(doi) => format!("{t} ({doi})"), - None => t, - }); - LeafRow { label, source: s.source, citation } - }) - .collect(); + // Count of placed non-D2C sample leaves at or below this node (shown instead of the list). + let sample_count = du_db::tree_sample::count_under(&st.pool, &name, dna_type).await?; Ok(html(&SnpSidebar { t: locale.t, @@ -442,8 +417,7 @@ async fn snp_sidebar( variants, str_changes, str_motif_markers, - samples, - samples_more, + sample_count, })) } diff --git a/rust/crates/du-web/templates/tree/snp_sidebar.html b/rust/crates/du-web/templates/tree/snp_sidebar.html index 0c14c59e..05b1f2b6 100644 --- a/rust/crates/du-web/templates/tree/snp_sidebar.html +++ b/rust/crates/du-web/templates/tree/snp_sidebar.html @@ -89,24 +89,11 @@
{% endif %} - {# Placed non-D2C sample leaves at or below this node (YFull-style). #} - {% if !samples.is_empty() %} + {# Count of placed non-D2C sample leaves at or below this node (YFull-style). #} + {% if sample_count > 0 %}
{{ t.get("tree.samples.title") }}
-
    - {% for s in samples %} -
  • - {{ s.label }} - {{ s.source }} - {% if let Some(c) = s.citation %} -
    {{ c }}
    - {% endif %} -
    • - {% endfor %} -
- {% if samples_more > 0 %} -
+{{ samples_more }} {{ t.get("tree.samples.more") }}
- {% endif %} + {{ sample_count }}
{% endif %} diff --git a/rust/locales/en.txt b/rust/locales/en.txt index 95c58168..a0c6c354 100644 --- a/rust/locales/en.txt +++ b/rust/locales/en.txt @@ -117,7 +117,6 @@ tree.legend.default=Other tree.variants=variants tree.samples=samples tree.samples.title=Placed samples -tree.samples.more=more tree.snp.none=No defining variants recorded. tree.snp.aliases=aka tree.snp.recurrent=recurrent diff --git a/rust/locales/es.txt b/rust/locales/es.txt index b0372a1e..ae09bf79 100644 --- a/rust/locales/es.txt +++ b/rust/locales/es.txt @@ -69,7 +69,6 @@ tree.legend.default=Otro tree.variants=variantes tree.samples=muestras tree.samples.title=Muestras ubicadas -tree.samples.more=más tree.snp.none=No hay variantes definitorias registradas. tree.snp.aliases=alias tree.snp.recurrent=recurrente diff --git a/rust/locales/fr.txt b/rust/locales/fr.txt index fe2c5f73..892d7ef0 100644 --- a/rust/locales/fr.txt +++ b/rust/locales/fr.txt @@ -69,7 +69,6 @@ tree.legend.default=Autre tree.variants=variants tree.samples=échantillons tree.samples.title=Échantillons placés -tree.samples.more=de plus tree.snp.none=Aucun variant déterminant enregistré. tree.snp.aliases=alias tree.snp.recurrent=récurrent