From db8881dc54263e9a86f752f45ad2dd3cfc0e7826 Mon Sep 17 00:00:00 2001 From: James Kane Date: Sun, 26 Jul 2026 07:25:23 -0500 Subject: [PATCH 1/4] feat(crawl): add CNCB/GSA linking support for raw seq files - Implemented GSA linking for CNCB studies, extending dual-archive coverage to East/Central-Asian BioProjects (`PRJCA`, `CRA`). - Added `GsaClient` for GSA open-tier crawls (`getRunInfoByCra` CSV parsing, `md5sum.txt` checksums). - Parametrized `ingest_libraries` to support GSA source alongside ENA. - Updated `crawl_project` to dispatch by accession prefix (`PRJCA/CRA` to GSA, `PRJEB/ENA` to ENA). - Introduced `pubs.study_source` migration for `"CNCB_GSA"` support. - Handled controlled-tier studies (`HRA`, `GVM`) by recording accession + DAR marker, skipping file ingest. - Verified prefix-agnostic `SAMC` biosample ingestion within existing infrastructure. Scoped to GSA open-tier only; controlled-tier workflows (GSA-Human/aspersa) deferred. --- rust/docs/cncb-gsa-sequence-linking.md | 237 +++++++++++++++++++++++++ 1 file changed, 237 insertions(+) create mode 100644 rust/docs/cncb-gsa-sequence-linking.md diff --git a/rust/docs/cncb-gsa-sequence-linking.md b/rust/docs/cncb-gsa-sequence-linking.md new file mode 100644 index 0000000..4a09708 --- /dev/null +++ b/rust/docs/cncb-gsa-sequence-linking.md @@ -0,0 +1,237 @@ +# CNCB / GSA sequence-file linking — design + +**Status:** proposal · **Date:** 2026-07-10 · **Author:** recon by Claude w/ James + +Adds a second external archive alongside ENA/NCBI for enumerating and linking raw +sequencing files to academic biosamples. The China National Center for +Bioinformation (CNCB) National Genomics Data Center (NGDC) is now the deposition +target for a growing body of East/Central-Asian popgen papers (the trigger here: +`GVM000900`, West China Hospital, 166 Central Asians, BioProject `PRJCA032194`). + +This doc scopes what we can link, how it maps onto our existing +`crawl-project` / `sequence_library` machinery, and the one structural change the +CNCB source forces (source is no longer implicitly "ENA"). + +--- + +## 1. Background: how CNCB is shaped + +CNCB/NGDC is a family of sub-archives bridged by a single **BioProject** +accession (`PRJCA…`, their `PRJEB/PRJNA` analog). The tiers that matter to us: + +| CNCB resource | Accession | Holds | Our analog | Linkable? | +|---|---|---|---|---| +| BioProject | `PRJCA…` | project umbrella | PRJEB/PRJNA | bridge only | +| **GSA** | `CRA`(study)/`CRX`(exp)/**`CRR`(run)**/`SAMC`(sample) | **open raw reads** (FASTQ/BAM) | ENA/SRA | **yes, anonymous** | +| GSA-Human | `HRA…` | controlled human raw reads | dbGaP/EGA | no — DAR + aspera key | +| GVM | `GVM…` | merged variant VCF/FASTA | dbSNP / joint VCF | open ones only; not reads | +| OMIX / GWH / OBIA | — | processed / assemblies | — | out of scope | + +**Two-tier reality.** Only the **open GSA (`CRA`)** tier is automatically +linkable. GSA-Human (`HRA`) and controlled-access GVM are DAR-gated (like +dbGaP/EGA): we can record the accession and a "controlled, DAR required" status, +but cannot fetch files. The linker must handle both — the happy path *and* the +metadata-only case — not assume every CNCB record yields files. + +### The `PRJCA032194` probe (why the two-tier handling is mandatory) + +Probed 2026-07-09. As of now this specific project exposes **nothing linkable**: + +- `getRunInfoByCra?searchTerm=PRJCA032194` and `getRunInfo` → **header-only CSV, + zero run rows** (no released GSA runs). +- GSA-Human "matches" (`HRA000087/150`, `CRA000112`) were **template-boilerplate + example accessions** — a nonsense control search returns the same ones. No real + `HRA` study is tied to the project. +- Only artifact on file is the controlled, embargoed **GVM VCF `GVM000900`** + (genotypes, not reads; release **2026-11-11**). + +So `PRJCA032194` is a controlled + embargoed + VCF-only record — the exact class +where the linker must degrade to "record accession + status," not file ingest. +Worth a re-probe after 2026-11-11. + +--- + +## 2. The programmatic surface (GSA open tier) + +No clean REST/JSON like ENA's `filereport`; it's POST-form → CSV, plus a +deterministic file host. Reference implementation: +[`BioOmics/iSeq`](https://github.com/BioOmics/iSeq) (`bin/iseq`). + +**Metadata (the `filereport` analog).** POST form-encoded, returns a ~25-column +CSV, **one row per file** (not per run — runs span multiple rows, grouped by the +`Run` column): + +``` +POST https://ngdc.cncb.ac.cn/gsa/search/getRunInfoByCra + searchTerm=&totalDatas=9999&downLoadCount=9999 +GET https://ngdc.cncb.ac.cn/gsa/search?searchTerm= # resolve run/exp → parent CRA +``` + +CSV columns (header verified live): + +``` +Run, Center, ReleaseDate, FileType, FileName, FileSize, Download_path, +Experiment, Title, LibraryName, LibraryStrategy, LibrarySelection, +LibrarySource, LibraryLayout, InsertSize, InsertDev, Platform, +BioProject, BioSample, SampleType, TaxID, ScientificName, SampleName, +Submission, Organization +``` + +Everything we need is in-band: `Download_path` gives the file URL directly (no +HTML scrape of `browse//` needed, unlike iSeq), plus `BioSample` +(`SAMC…`), `Platform`, `LibraryLayout`, `FileType`, `FileSize`. + +**File URLs** — deterministic, multi-mirror: + +``` +https://download.cncb.ac.cn/gsa///_f1.fq.gz # HTTPS (store this) +ftp://download.big.ac.cn/gsa///_f1.fq.gz # FTP mirror ++ HuaweiCloud / Qtrans accelerator mirrors when present +``` + +`` is a volume segment (`gsa`, `gsa2`, `gsa3`, …) — **do not hardcode**; take +it from the CSV `Download_path`. + +**Checksums** — no per-row MD5. One study-level manifest: +`https://download.cncb.ac.cn/gsa//md5sum.txt` (filename→md5). One fetch +per study, index by `FileName`. + +### Caveats vs ENA + +- POST-form + CSV + a separate md5 manifest — no single JSON call. +- Endpoint already renamed once (`getRunInfo` → `getRunInfoByCra`, Sep 2024) — + pin endpoint paths in config. +- One row **per file**; paired FASTQ = two rows sharing a `Run`. (ENA gives one + row per run with `;`-joined files.) Grouping is `Run → Sample`. +- Download hosts are CN-resident; expect slow HTTPS. We only *store* URLs, we + don't fetch bytes, so this is a curator/end-user concern, not ours. +- No credentials for open GSA. GSA-Human needs an aspera key + approved DAR. + +--- + +## 3. How it maps onto our existing machinery + +We already do exactly this for ENA in `crawl_project.rs` +(`du_external::ena` → group runs by sample → `biosample::upsert_by_accession` +→ `publication::link_biosample` → `sequence::ingest_libraries`). CNCB reuses the +whole write path; only the *fetch + parse* and the *source label* change. + +### 3.1 New: `du-external/src/cncb.rs` + +Mirror `ena.rs`. A `GsaClient` with: + +- `study_runs(cra: &str) -> Vec` — POST `getRunInfoByCra`, parse CSV + by header (robust to column reorder, same as `ena::parse_run_report`). +- `resolve_cra(run_or_exp: &str) -> Vec` — for `CRR`/`CRX`/`PRJCA` inputs. +- `md5_manifest(cra, volume) -> HashMap` — fetch + parse + `md5sum.txt` once per study. +- `GsaRunRow` fields mirror the CSV columns above. Parsing pure + unit-tested + against a captured CSV fixture; HTTP is a thin wrapper. (One row per file → + the client can pre-group into runs, or leave grouping to the job as ENA does.) + +### 3.2 Changed: `sequence::ingest_libraries` — parametrize the source + +Today the JSONB provenance is **hardcoded `"source": "ENA"`** +(`sequence.rs`, the `atproto` slot): + +```rust +.bind(json!({ "source": "ENA", "run_accession": lib.external_run_ref })) +``` + +Add a `source: &str` (or a small `enum SeqSource { Ena, Gsa }`) to +`NewSeqLibrary` (or the `ingest_libraries` call) so GSA runs record +`{"source": "GSA", "run_accession": "CRR…"}`. This is the one required change to +shared code; everything downstream (`biosample::report` read path, the sidecar +`http_locations`/`checksums` JSONB shapes) is source-agnostic already. + +### 3.3 Changed: `crawl_project.rs` — dispatch by accession shape + +`crawl_one_accession` / `crawl_pending` currently assume ENA. Options: + +- **(A, preferred) Dispatch inside the existing job.** Detect archive by + accession prefix — `PRJEB/PRJNA/ERP/SRP` → ENA path; `PRJCA/CRA` → GSA path — + and route to the matching client. Keeps one `crawl-project` command and one + pending-study drain; both clients feed the same `group-by-sample → upsert → + link → ingest` core. Factor that core to take a + `Vec<(sample_acc, Vec)>` so ENA and GSA just produce it + differently. +- (B) A parallel `crawl-project-cncb` run-once command. Less code churn now, + duplicates the drain loop and the study-table plumbing. Reject unless (A) + proves awkward. + +### 3.4 Required migration: extend `pubs.study_source` + +`study::upsert_by_accession` casts its source `::pubs.study_source`, and that +enum is `('ENA', 'NCBI_BIOPROJECT', 'NCBI_GENBANK')` (mig `0006_pubs.sql`) — **no +GSA value**. A `CRA`/`PRJCA` study insert would fail the cast. Add a migration: + +```sql +ALTER TYPE pubs.study_source ADD VALUE IF NOT EXISTS 'CNCB_GSA'; +``` + +(`ADD VALUE` can't run inside the same txn that uses it, and can't be dropped — +standard enum-extension caveats; add it in its own migration.) The accession +*string* itself is fine — `genomic_study.accession` is verbatim + unique, no +prefix logic. This plus the `sequence` source param (§3.2) are the two required +data-layer changes; `sequence_library` provenance is free JSONB (no migration). + +### 3.5 Biosample source — verified prefix-agnostic ✓ + +Crawl-created samples stay `SAMPLE_SOURCE = "EXTERNAL"` (public/academic), same +as ENA. The CNCB `SAMC…` BioSample accession goes in via +`upsert_by_accession(pool, samc_acc, "EXTERNAL", None)`. **Verified** +(`biosample.rs`): the upsert stores `accession.trim()` verbatim, unique on +`accession`, with no `SAM[END]`/prefix assumption anywhere in the write or report +read path; `core.biosample_source` already includes `EXTERNAL`. `SAMC…` needs no +change. (Open question #3 — resolved.) + +--- + +## 4. Controlled tier (GSA-Human / GVM-controlled) + +Out of scope for automated file ingest. Minimum useful handling: + +- When a crawl resolves a BioProject to only an `HRA` study or a controlled + `GVM`, **record the accession + a `controlled: true, dar_required: true` + marker** on the study/publication so a curator knows a manual DAR is the gate + (analogous to how we'd treat dbGaP/EGA). +- Do **not** attempt aspera/credentialed fetch. If we ever pursue it, it's a + separate manual-curation workflow (register → DAR to the study's DAC → aspera + key), not a job. + +--- + +## 5. Scope / cut list + +**In scope (first slice):** +- `du-external::cncb::GsaClient` (getRunInfoByCra CSV + md5 manifest), unit-tested + against a captured fixture. +- Migration: `ALTER TYPE pubs.study_source ADD VALUE 'CNCB_GSA'` (§3.4). +- `sequence::ingest_libraries` source parametrization (`ENA` | `GSA`). +- `crawl_project` dispatch-by-prefix, GSA open path reaching `ingest_libraries`. +- Controlled-tier: record accession + DAR marker, no fetch. + +**Deferred / out:** +- GSA-Human aspera/DAR fetch (manual curation, not a job). +- GVM variant-VCF ingestion (variant tier — separate from read linking; overlaps + ybrowse/reconcile territory, not this doc). +- OMIX/GWH/OBIA. +- A standalone ops python driver — the Rust `crawl-project` job is the driver; + no `link-gsa-sequence-files.py` needed (the ENA python script predated the job). + +--- + +## 6. Open questions / verify before building + +1. **Re-probe `PRJCA032194` after 2026-11-11** — does an open `CRA` study appear, + or does it stay HRA/GVM-controlled? Decides whether the trigger paper is ever + auto-linkable or is purely a controlled-tier example. +2. **Paired-FASTQ grouping** — confirm on a real multi-file `CRA` that rows share + the `Run` value and that `_f1/_r2` (or `_1/_2`) is the reliable mate signal. +3. ~~**`biosample` accession namespace** — verify `SAMC…` stored verbatim.~~ + **Resolved (2026-07-10):** prefix-agnostic, `EXTERNAL` source exists, no change + (§3.5). Verifying this also surfaced the required `study_source` migration (§3.4). +4. **Rate/politeness** — reuse the ENA `REQUEST_GAP` (150ms) + bounded batch; + confirm CNCB tolerates it (CN latency may want a longer gap). +5. **Endpoint stability** — the `getRunInfo*` rename history argues for + config-pinned endpoint paths + a fixture-refresh check in CI. From 3979d9337a3296ec2650779eae9e01a4db491d1a Mon Sep 17 00:00:00 2001 From: James Kane Date: Wed, 29 Jul 2026 05:08:20 -0500 Subject: [PATCH 2/4] feat(str): corpus-wide Y-STR marker report + mutation-rate provenance MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds /str-markers and /api/v1/reports/str-markers: the observed repeat range (min / modal / max) of every Y-STR marker across all profiles we hold, joined to its repeat motif and mutation rate. Until now the STR data was only ever read per-haplogroup (modal signatures, ASR, branch ages), so there was no way to ask what the corpus says about a marker. Corpus-wide only — the clade filter can be added later without a schema change. Migration 0069 adds genomics.str_marker, the first STR marker *reference* in the schema; marker identity was previously a bare marker_name TEXT everywhere. Seeded (142 rows, 125 with a motif) by parsing the free-text str_mutation_rate.source, the only motif information this database holds. `coordinates` mirrors the core.genome_region shape so a locus source can be loaded in place later; it is empty today. Marker-name aliasing is load-bearing, not defensive: profiles say DYS19 and Y-GATA-H4 where the rate table says DYS394 and YGATAH4, so a join on name alone silently drops the two most recognizable markers in the panel. DYS19/DYS394 are also duplicate rate rows for one locus, so the seed dedupes them. Every lookup folds through str_marker.aliases. Migration 0070 re-keys str_mutation_rate on (marker_name, method) with a partial unique index on the single active rate per marker, so a rate derived from our own tree can coexist with the published one instead of overwriting it; load_marker_models now reads WHERE is_active, making promotion a data operation rather than a code change. It also records, per branch age, which rates produced it. That split is worth having: only 114 of 829 observed markers have a published rate, so across all 1,963 STR_VARIANCE branches 1,393,730 marker-scorings fall back to DEFAULT_STR_RATE against 221,819 measured — previously invisible. A reported 0 is FTDNA's null-allele (deleted locus) convention rather than a zero-repeat count, so it is excluded from min/modal/max and counted separately; without that DYS448 reads min 0 instead of 17 and CDY reads min '0-0'. The same rule applies per copy of a multi-copy vector. Deriving rates from accumulating observations is designed but NOT built here. The estimator must take branch lengths from method='SNP_POISSON' ages only — COMBINED or STR_VARIANCE would feed STR-derived ages back into STR rate estimation and let the model converge on itself. Co-Authored-By: Claude Opus 5 (1M context) --- rust/crates/du-db/src/ystr.rs | 191 +++++++++++- rust/crates/du-db/tests/str_marker_stats.rs | 197 ++++++++++++ rust/crates/du-web/src/api/dto.rs | 78 +++++ rust/crates/du-web/src/api/mod.rs | 11 +- rust/crates/du-web/src/routes/mod.rs | 2 + rust/crates/du-web/src/routes/str_markers.rs | 285 ++++++++++++++++++ rust/crates/du-web/templates/base.html | 1 + rust/crates/du-web/templates/static/page.html | 1 + rust/crates/du-web/templates/str/markers.html | 87 ++++++ rust/locales/en.txt | 43 +++ rust/locales/es.txt | 43 +++ rust/locales/fr.txt | 43 +++ rust/migrations/0069_str_marker.sql | 96 ++++++ rust/migrations/0070_str_rate_provenance.sql | 57 ++++ rust/scripts/seed-str-mutation-rates.sql | 4 +- 15 files changed, 1127 insertions(+), 12 deletions(-) create mode 100644 rust/crates/du-db/tests/str_marker_stats.rs create mode 100644 rust/crates/du-web/src/routes/str_markers.rs create mode 100644 rust/crates/du-web/templates/str/markers.html create mode 100644 rust/migrations/0069_str_marker.sql create mode 100644 rust/migrations/0070_str_rate_provenance.sql diff --git a/rust/crates/du-db/src/ystr.rs b/rust/crates/du-db/src/ystr.rs index 1d5607f..9387a26 100644 --- a/rust/crates/du-db/src/ystr.rs +++ b/rust/crates/du-db/src/ystr.rs @@ -315,6 +315,9 @@ fn p_g_given_m(g: i64, m: i64) -> f64 { #[derive(Clone)] pub struct MarkerModel { pub mu_per_gen: f64, + /// Provenance of `mu_per_gen`: `PUBLISHED`/`DERIVED` from the rate row, or + /// `DEFAULT` for the [`DEFAULT_STR_RATE`] fallback used where no rate exists. + pub method: String, omega_pgm: Option>>>, } @@ -394,7 +397,7 @@ pub fn compute_str_age( if modal_simple.is_empty() { return None; } - let default = MarkerModel { mu_per_gen: DEFAULT_STR_RATE, omega_pgm: None }; + let default = MarkerModel { mu_per_gen: DEFAULT_STR_RATE, method: "DEFAULT".into(), omega_pgm: None }; // Per-marker age PDFs depend only on (marker, g) — memoize across testers. let mut cache: HashMap<(&str, i64), Option> = HashMap::new(); let mut acc: Option = None; @@ -624,7 +627,7 @@ fn propagate_str( max_age: f64, ) -> Vec> { let n = children.len(); - let default = MarkerModel { mu_per_gen: DEFAULT_STR_RATE, omega_pgm: None }; + let default = MarkerModel { mu_per_gen: DEFAULT_STR_RATE, method: "DEFAULT".into(), omega_pgm: None }; let mut cache: HashMap<(String, i64), Option> = HashMap::new(); let mut memo: Vec>> = vec![None; n]; for i in 0..n { @@ -767,11 +770,18 @@ pub async fn str_tmrca_pdfs(pool: &PgPool, res: f64, max_age: f64) -> Result Result, DbError> { let rows: Vec = sqlx::query_as( "SELECT marker_name, mutation_rate::float8 AS mutation_rate, \ omega_plus::float8 AS omega_plus, omega_minus::float8 AS omega_minus, \ - multi_step_rate::float8 AS multi_step_rate FROM genomics.str_mutation_rate", + multi_step_rate::float8 AS multi_step_rate, method \ + FROM genomics.str_mutation_rate WHERE is_active", ) .fetch_all(pool) .await?; @@ -790,7 +800,7 @@ pub async fn load_marker_models(pool: &PgPool) -> Result None, }; - (r.marker_name, MarkerModel { mu_per_gen: r.mutation_rate, omega_pgm }) + (r.marker_name, MarkerModel { mu_per_gen: r.mutation_rate, method: r.method, omega_pgm }) }) .collect()) } @@ -802,6 +812,7 @@ struct MarkerRateRow { omega_plus: Option, omega_minus: Option, multi_step_rate: Option, + method: String, } #[derive(Debug, Default)] @@ -878,16 +889,37 @@ pub async fn recompute_signatures(pool: &PgPool) -> Result = BTreeSet::new(); + let measured = motif[i] + .keys() + .filter_map(|m| models.get(m)) + .inspect(|m| { + methods.insert(m.method.as_str()); + }) + .count() as i32; + let rate_method = match methods.len() { + 0 => None, + 1 => methods.iter().next().map(|s| s.to_string()), + _ => Some("MIXED".to_string()), + }; sqlx::query( "INSERT INTO tree.haplogroup_age_estimate \ (haplogroup_id, method, estimate_ybp, ci_low_ybp, ci_high_ybp, \ - sample_count, marker_count, generation_years, computed_at) \ - VALUES ($1, 'STR_VARIANCE', $2, $3, $4, $5, $6, $7, now()) \ + sample_count, marker_count, generation_years, computed_at, \ + rate_method, measured_rate_markers, default_rate_markers) \ + VALUES ($1, 'STR_VARIANCE', $2, $3, $4, $5, $6, $7, now(), $8, $9, $10) \ ON CONFLICT (haplogroup_id, method) DO UPDATE SET \ estimate_ybp = EXCLUDED.estimate_ybp, ci_low_ybp = EXCLUDED.ci_low_ybp, \ ci_high_ybp = EXCLUDED.ci_high_ybp, sample_count = EXCLUDED.sample_count, \ marker_count = EXCLUDED.marker_count, generation_years = EXCLUDED.generation_years, \ - computed_at = now()", + computed_at = now(), rate_method = EXCLUDED.rate_method, \ + measured_rate_markers = EXCLUDED.measured_rate_markers, \ + default_rate_markers = EXCLUDED.default_rate_markers", ) .bind(ids[i]) .bind(med.round() as i32) @@ -896,6 +928,9 @@ pub async fn recompute_signatures(pool: &PgPool) -> Result, pub marker_count: Option, pub generation_years: Option, + /// Provenance of the mutation rates behind an `STR_VARIANCE` estimate: + /// `PUBLISHED`, `DERIVED`, or `MIXED`. None for other methods, and for an + /// estimate whose markers all fell back to [`DEFAULT_STR_RATE`]. + pub rate_method: Option, + pub measured_rate_markers: Option, + /// Markers scored at [`DEFAULT_STR_RATE`] for want of a rate row — a large + /// share means the estimate is weakly grounded in measured rates. + pub default_rate_markers: Option, } /// Contributing age estimates for a Y haplogroup (e.g. STR_VARIANCE). pub async fn branch_age_estimates(pool: &PgPool, haplogroup: &str) -> Result, DbError> { let rows = sqlx::query_as::<_, AgeEstimate>( "SELECT e.method, e.estimate_ybp, e.ci_low_ybp, e.ci_high_ybp, e.sample_count, \ - e.marker_count, e.generation_years::float8 AS generation_years \ + e.marker_count, e.generation_years::float8 AS generation_years, \ + e.rate_method, e.measured_rate_markers, e.default_rate_markers \ FROM tree.haplogroup_age_estimate e \ JOIN tree.haplogroup h ON h.id = e.haplogroup_id \ WHERE h.name = $1 AND h.haplogroup_type::text = 'Y_DNA' AND h.valid_until IS NULL \ @@ -1042,6 +1086,133 @@ pub async fn branch_age_estimates(pool: &PgPool, haplogroup: &str) -> Result, + pub modal_value: Option, + pub max_value: Option, + // Multi-copy markers: a rendered copy vector, e.g. "11-15". None for simple. + pub min_combination: Option, + pub modal_combination: Option, + pub max_combination: Option, + pub distinct_values: i64, + /// Observations reporting a deleted locus (a `0` repeat count, or any `0` + /// copy) — excluded from min/modal/max. + pub null_alleles: i64, + /// Partial-repeat / footnoted values (`{type:complex}`), counted but unscored. + pub complex_count: i64, + pub motif: Option, + pub period: Option, + pub coordinates: Option, + pub mutation_rate: Option, + pub rate_ci_low: Option, + pub rate_ci_high: Option, + pub rate_method: Option, + pub rate_source: Option, + /// How this marker enters the STR age model: `MEASURED_RATE` (an active rate + /// row), `DEFAULT_RATE` (falls back to [`DEFAULT_STR_RATE`]), or `EXCLUDED` + /// (multi-copy — never scored). + pub age_model_status: String, +} + +/// Corpus-wide observed range per Y-STR marker: min, modal and max across every +/// profile we hold, joined to the marker registry (motif/period) and the active +/// mutation rate. Ordered by observation count, descending. +/// +/// Both profile sources feed it — the federated mirror and the locally-imported +/// vendor exports — the same union [`build_str_inputs`] uses, minus the tree join +/// (this is a corpus summary, not a per-branch one). +/// +/// A reported `0` is FTDNA's null-allele (deleted locus) convention rather than a +/// zero-repeat count, so it is excluded from min/modal/max and counted in +/// `null_alleles`; without that, DYS448 reads min 0 instead of 17. The same rule +/// applies per copy of a multi-copy vector. +pub async fn marker_stats(pool: &PgPool) -> Result, DbError> { + let rows = sqlx::query_as::<_, MarkerStat>( + "WITH obs AS ( \ + SELECT p.sample_guid::text AS sid, e->>'marker' AS marker, e->'value' AS v \ + FROM genomics.biosample_str_profile p, jsonb_array_elements(p.markers) e \ + UNION ALL \ + SELECT sp.did || '/' || sp.rkey, e->>'marker', e->'value' \ + FROM fed.str_profile sp, jsonb_array_elements(sp.markers) e \ + ), typed AS ( \ + SELECT sid, marker, v->>'type' AS vtype, \ + CASE WHEN v->>'type' = 'simple' THEN (v->>'repeats')::int END AS repeats, \ + CASE WHEN v->>'type' = 'multiCopy' \ + THEN ARRAY(SELECT jsonb_array_elements_text(v->'copies')::int) END AS copies, \ + v->>'raw' AS raw \ + FROM obs \ + ), scored AS ( \ + SELECT *, (repeats IS NOT NULL AND repeats = 0) \ + OR (copies IS NOT NULL AND 0 = ANY(copies)) AS is_null_allele \ + FROM typed \ + ), agg AS ( \ + SELECT marker, \ + count(*) AS observations, \ + count(DISTINCT sid) AS samples, \ + COALESCE(bool_or(vtype = 'multiCopy'), false) AS observed_multi, \ + min(repeats) FILTER (WHERE NOT is_null_allele) AS min_value, \ + max(repeats) FILTER (WHERE NOT is_null_allele) AS max_value, \ + mode() WITHIN GROUP (ORDER BY repeats) \ + FILTER (WHERE NOT is_null_allele) AS modal_value, \ + (array_agg(array_to_string(copies, '-') ORDER BY copies) \ + FILTER (WHERE copies IS NOT NULL AND NOT is_null_allele))[1] AS min_copies, \ + (array_agg(array_to_string(copies, '-') ORDER BY copies DESC) \ + FILTER (WHERE copies IS NOT NULL AND NOT is_null_allele))[1] AS max_copies, \ + mode() WITHIN GROUP (ORDER BY copies) \ + FILTER (WHERE copies IS NOT NULL AND NOT is_null_allele) AS modal_copies, \ + count(DISTINCT COALESCE(repeats::text, array_to_string(copies, '-'), raw)) \ + AS distinct_values, \ + count(*) FILTER (WHERE is_null_allele) AS null_alleles, \ + count(*) FILTER (WHERE vtype = 'complex') AS complex_count \ + FROM scored \ + GROUP BY marker \ + ) \ + SELECT a.marker AS marker_name, \ + (a.observed_multi OR COALESCE(m.multi_copy, false)) AS multi_copy, \ + a.observations, a.samples, \ + a.min_value, a.modal_value, a.max_value, \ + a.min_copies AS min_combination, \ + array_to_string(a.modal_copies, '-') AS modal_combination, \ + a.max_copies AS max_combination, \ + a.distinct_values, a.null_alleles, a.complex_count, \ + m.motif, m.period, NULLIF(m.coordinates, '{}'::jsonb) AS coordinates, \ + r.mutation_rate::float8 AS mutation_rate, \ + r.mutation_rate_lower::float8 AS rate_ci_low, \ + r.mutation_rate_upper::float8 AS rate_ci_high, \ + r.method AS rate_method, r.source AS rate_source, \ + CASE WHEN a.observed_multi OR COALESCE(m.multi_copy, false) THEN 'EXCLUDED' \ + WHEN r.id IS NOT NULL THEN 'MEASURED_RATE' \ + ELSE 'DEFAULT_RATE' END AS age_model_status \ + FROM agg a \ + LEFT JOIN LATERAL ( \ + SELECT sm.* FROM genomics.str_marker sm \ + WHERE sm.marker_name = a.marker OR a.marker = ANY(sm.aliases) \ + ORDER BY (sm.marker_name <> a.marker) LIMIT 1 \ + ) m ON true \ + LEFT JOIN LATERAL ( \ + SELECT mr.* FROM genomics.str_mutation_rate mr \ + WHERE mr.is_active \ + AND (mr.marker_name = a.marker \ + OR mr.marker_name = ANY(COALESCE(m.aliases, '{}'))) \ + ORDER BY (mr.marker_name <> a.marker) LIMIT 1 \ + ) r ON true \ + ORDER BY a.observations DESC, a.marker", + ) + .fetch_all(pool) + .await?; + Ok(rows) +} + /// A ranked STR→branch prediction. #[derive(Debug, Clone, PartialEq)] pub struct Prediction { @@ -1164,7 +1335,7 @@ mod tests { } fn model(mu: f64) -> MarkerModel { - MarkerModel { mu_per_gen: mu, omega_pgm: None } + MarkerModel { mu_per_gen: mu, method: "PUBLISHED".into(), omega_pgm: None } } #[test] @@ -1266,7 +1437,7 @@ mod tests { let motifs = vec![flat(10), flat(11), flat(12)]; // root, mid, leaf let testers = vec![vec![], vec![], vec![flat(12)]]; // leaf tip on its motif (g=0) let models: HashMap = (1..=5) - .map(|i| (format!("DYS{i}"), MarkerModel { mu_per_gen: 0.01, omega_pgm: None })) + .map(|i| (format!("DYS{i}"), MarkerModel { mu_per_gen: 0.01, method: "PUBLISHED".into(), omega_pgm: None })) .collect(); let ages = propagate_str( &children, diff --git a/rust/crates/du-db/tests/str_marker_stats.rs b/rust/crates/du-db/tests/str_marker_stats.rs new file mode 100644 index 0000000..ca5af03 --- /dev/null +++ b/rust/crates/du-db/tests/str_marker_stats.rs @@ -0,0 +1,197 @@ +//! Live-DB test for the corpus-wide Y-STR marker report (`du_db::ystr::marker_stats`, +//! migrations 0069/0070). +//! +//! Runs against an isolated ephemeral database rather than the shared dev one: the +//! report aggregates *every* profile in the catalog, so assertions about exact +//! min/modal/max values are only meaningful when this test owns the whole corpus. +//! A fresh database also has an empty `genomics.str_marker` — the 0069 seed derives +//! from `genomics.str_mutation_rate`, which is populated by a script, not a +//! migration — so the reference rows are seeded here explicitly. +//! Re-runnable; skips (passes) when DATABASE_URL is unset. +//! +//! eval "$(./scripts/test-db.sh up)" && cargo test -p du-db --test str_marker_stats + +use du_db::ystr::{self, MarkerStat}; +use sqlx::PgPool; + +fn database_url() -> Option { + std::env::var("DATABASE_URL").ok().filter(|s| !s.is_empty()) +} + +/// One profile: a biosample plus its `strMarkerValue[]`, in the lexicon shape both +/// the federated mirror and the vendor importer write. +async fn seed_profile(pool: &PgPool, accession: &str, markers: serde_json::Value) { + let guid: uuid::Uuid = sqlx::query_scalar( + "INSERT INTO core.biosample (source, accession) \ + VALUES ('EXTERNAL'::core.biosample_source, $1) RETURNING sample_guid", + ) + .bind(accession) + .fetch_one(pool) + .await + .unwrap(); + let total = markers.as_array().map(|a| a.len()).unwrap_or(0) as i32; + sqlx::query( + "INSERT INTO genomics.biosample_str_profile (sample_guid, total_markers, markers) \ + VALUES ($1, $2, $3)", + ) + .bind(guid) + .bind(total) + .bind(&markers) + .execute(pool) + .await + .unwrap(); +} + +fn simple(marker: &str, repeats: i32) -> serde_json::Value { + serde_json::json!({ "marker": marker, "value": { "type": "simple", "repeats": repeats } }) +} + +fn multi(marker: &str, copies: &[i32]) -> serde_json::Value { + serde_json::json!({ "marker": marker, "value": { "type": "multiCopy", "copies": copies } }) +} + +fn find<'a>(stats: &'a [MarkerStat], marker: &str) -> &'a MarkerStat { + stats.iter().find(|m| m.marker_name == marker).unwrap_or_else(|| panic!("no row for {marker}")) +} + +#[tokio::test] +async fn marker_stats_reports_range_motif_and_rate_basis() { + let Some(url) = database_url() else { + eprintln!("DATABASE_URL unset — skipping marker_stats test"); + return; + }; + let db = du_db::testing::ephemeral_db(&url).await.expect("ephemeral db"); + let pool = db.pool().clone(); + + // Reference rows. DYS19 is registered under the name profiles use, carrying the + // literature's DYS394 designation as an alias — the rate below is filed under + // DYS394, so only alias folding can connect the two. + sqlx::query( + "INSERT INTO genomics.str_marker (marker_name, motif, period, aliases) \ + VALUES ('DYS19', 'AGAT', 4, ARRAY['DYS394'])", + ) + .execute(&pool) + .await + .unwrap(); + sqlx::query( + "INSERT INTO genomics.str_mutation_rate (marker_name, mutation_rate, method) \ + VALUES ('DYS394', 0.0028513662, 'PUBLISHED')", + ) + .execute(&pool) + .await + .unwrap(); + + // Three testers. DYS19 spans 12–14 with 14 twice (the mode). DYS448 carries a + // reported 0 — a deleted locus, not a zero repeat count. DYS385 is multi-copy. + // DYS447 has no rate row at all. DYS442 carries a partial repeat. + seed_profile( + &pool, + "STRTEST-1", + serde_json::json!([ + simple("DYS19", 12), + simple("DYS448", 0), + simple("DYS447", 24), + multi("DYS385", &[11, 15]), + serde_json::json!({ "marker": "DYS442", "value": { "type": "complex", "raw": "10.2" } }), + ]), + ) + .await; + seed_profile( + &pool, + "STRTEST-2", + serde_json::json!([ + simple("DYS19", 14), + simple("DYS448", 20), + simple("DYS447", 25), + multi("DYS385", &[11, 14]), + ]), + ) + .await; + seed_profile( + &pool, + "STRTEST-3", + serde_json::json!([ + simple("DYS19", 14), + simple("DYS448", 22), + simple("DYS447", 25), + multi("DYS385", &[11, 14]), + ]), + ) + .await; + + let stats = ystr::marker_stats(&pool).await.expect("marker_stats"); + + // Simple marker: min/modal/max, and the motif + rate reached through the alias. + let dys19 = find(&stats, "DYS19"); + assert_eq!((dys19.min_value, dys19.modal_value, dys19.max_value), (Some(12), Some(14), Some(14))); + assert_eq!(dys19.observations, 3); + assert_eq!(dys19.samples, 3); + assert_eq!(dys19.distinct_values, 2); + assert_eq!(dys19.motif.as_deref(), Some("AGAT")); + assert_eq!(dys19.period, Some(4)); + assert_eq!(dys19.age_model_status, "MEASURED_RATE"); + assert_eq!(dys19.rate_method.as_deref(), Some("PUBLISHED")); + + // A reported 0 is excluded from the range and counted as a null allele — + // without this the minimum reads 0 rather than 20. + let dys448 = find(&stats, "DYS448"); + assert_eq!(dys448.min_value, Some(20), "null allele must not become the minimum"); + assert_eq!(dys448.max_value, Some(22)); + assert_eq!(dys448.null_alleles, 1); + assert_eq!(dys448.observations, 3, "the null allele is still an observation"); + + // Multi-copy: rendered copy vectors, no scalar values, never age-scored. + let dys385 = find(&stats, "DYS385"); + assert!(dys385.multi_copy); + assert_eq!(dys385.min_value, None); + assert_eq!(dys385.min_combination.as_deref(), Some("11-14")); + assert_eq!(dys385.modal_combination.as_deref(), Some("11-14")); + assert_eq!(dys385.max_combination.as_deref(), Some("11-15")); + assert_eq!(dys385.age_model_status, "EXCLUDED"); + + // No rate row ⇒ the age model falls back to DEFAULT_STR_RATE, and says so. + let dys447 = find(&stats, "DYS447"); + assert_eq!(dys447.age_model_status, "DEFAULT_RATE"); + assert_eq!(dys447.mutation_rate, None); + assert_eq!((dys447.min_value, dys447.modal_value, dys447.max_value), (Some(24), Some(25), Some(25))); + + // Partial repeats are preserved and counted, but never scored. + let dys442 = find(&stats, "DYS442"); + assert_eq!(dys442.complex_count, 1); + assert_eq!(dys442.min_value, None); + assert_eq!(dys442.modal_value, None); +} + +#[tokio::test] +async fn load_marker_models_honors_the_active_rate() { + let Some(url) = database_url() else { + eprintln!("DATABASE_URL unset — skipping active-rate test"); + return; + }; + let db = du_db::testing::ephemeral_db(&url).await.expect("ephemeral db"); + let pool = db.pool().clone(); + + // A published rate alongside one derived from our own tree, as the re-keyed + // table now allows. Only the active row may reach the age model — promoting a + // derived rate is a data change, not a code change. + sqlx::query( + "INSERT INTO genomics.str_mutation_rate \ + (marker_name, mutation_rate, method, is_active) VALUES \ + ('DYSTEST1', 0.0030, 'PUBLISHED', true), \ + ('DYSTEST1', 0.0090, 'DERIVED', false), \ + ('DYSTEST2', 0.0040, 'DERIVED', true)", + ) + .execute(&pool) + .await + .unwrap(); + + let models = ystr::load_marker_models(&pool).await.expect("load_marker_models"); + + let m1 = models.get("DYSTEST1").expect("DYSTEST1 model"); + assert_eq!(m1.mu_per_gen, 0.0030, "the inactive DERIVED rate must not win"); + assert_eq!(m1.method, "PUBLISHED"); + + let m2 = models.get("DYSTEST2").expect("DYSTEST2 model"); + assert_eq!(m2.mu_per_gen, 0.0040); + assert_eq!(m2.method, "DERIVED", "an active derived rate is used, and labeled as such"); +} diff --git a/rust/crates/du-web/src/api/dto.rs b/rust/crates/du-web/src/api/dto.rs index 4f068da..c5afd98 100644 --- a/rust/crates/du-web/src/api/dto.rs +++ b/rust/crates/du-web/src/api/dto.rs @@ -645,6 +645,13 @@ pub struct AgeEstimateDto { pub sample_count: Option, pub marker_count: Option, pub generation_years: Option, + /// Provenance of the mutation rates behind an STR_VARIANCE estimate: + /// PUBLISHED, DERIVED or MIXED (null for other methods). + pub rate_method: Option, + pub measured_rate_markers: Option, + /// Markers scored at the default rate for want of a published one — a large + /// share means the estimate is only weakly grounded in measured rates. + pub default_rate_markers: Option, } impl From for AgeEstimateDto { @@ -657,6 +664,77 @@ impl From for AgeEstimateDto { sample_count: e.sample_count, marker_count: e.marker_count, generation_years: e.generation_years, + rate_method: e.rate_method, + measured_rate_markers: e.measured_rate_markers, + default_rate_markers: e.default_rate_markers, + } + } +} + +/// One Y-STR marker's corpus-wide observed range, with the reference and rate +/// information we hold for it (the `/str-markers` report). +#[derive(Serialize, ToSchema)] +pub struct StrMarkerStatDto { + pub marker: String, + /// Palindromic/duplicated locus: reported as a copy vector, never age-scored. + pub multi_copy: bool, + pub observations: i64, + pub samples: i64, + /// Simple markers — repeat counts. Null for multi-copy markers. + pub min_value: Option, + pub modal_value: Option, + pub max_value: Option, + /// Multi-copy markers — rendered copy vectors, e.g. "11-15". Null for simple. + pub min_combination: Option, + pub modal_combination: Option, + pub max_combination: Option, + pub distinct_values: i64, + /// Observations reporting a deleted locus (a `0` repeat count, or any `0` + /// copy). Excluded from min/modal/max. + pub null_alleles: i64, + /// Partial-repeat / footnoted values, preserved but unscored. + pub complex_count: i64, + /// Repeat unit, where known — we hold one for a minority of markers. + pub motif: Option, + pub period: Option, + /// Per-build locus coordinates. Always null today: no source loaded yet. + #[schema(value_type = Object)] + pub coordinates: Option, + pub mutation_rate: Option, + pub rate_ci_low: Option, + pub rate_ci_high: Option, + pub rate_method: Option, + pub rate_source: Option, + /// How this marker enters the STR age model: MEASURED_RATE, DEFAULT_RATE + /// (no rate row — falls back to the default), or EXCLUDED (multi-copy). + pub age_model_status: String, +} + +impl From for StrMarkerStatDto { + fn from(m: du_db::ystr::MarkerStat) -> Self { + StrMarkerStatDto { + marker: m.marker_name, + multi_copy: m.multi_copy, + observations: m.observations, + samples: m.samples, + min_value: m.min_value, + modal_value: m.modal_value, + max_value: m.max_value, + min_combination: m.min_combination, + modal_combination: m.modal_combination, + max_combination: m.max_combination, + distinct_values: m.distinct_values, + null_alleles: m.null_alleles, + complex_count: m.complex_count, + motif: m.motif, + period: m.period, + coordinates: m.coordinates, + mutation_rate: m.mutation_rate, + rate_ci_low: m.rate_ci_low, + rate_ci_high: m.rate_ci_high, + rate_method: m.rate_method, + rate_source: m.rate_source, + age_model_status: m.age_model_status, } } } diff --git a/rust/crates/du-web/src/api/mod.rs b/rust/crates/du-web/src/api/mod.rs index 78d6e13..9316b08 100644 --- a/rust/crates/du-web/src/api/mod.rs +++ b/rust/crates/du-web/src/api/mod.rs @@ -177,6 +177,13 @@ async fn reports_ancestry(State(st): State) -> Result) -> Result>, AppError> { + let rows = du_db::ystr::marker_stats(&st.pool).await?; + Ok(Json(rows.into_iter().map(StrMarkerStatDto::from).collect())) +} + #[utoipa::path(get, path = "/api/v1/reports/haplogroups", tag = "reports", responses((status = 200, description = "Y/MT haplogroup distribution across mirrored biosamples", body = [HaplogroupCountDto])))] async fn reports_haplogroups(State(st): State) -> Result>, AppError> { @@ -514,7 +521,7 @@ fn csv_field(s: &str) -> String { y_tree, mt_tree, y_tree_full, mt_tree_full, y_tree_version, mt_tree_version, y_node_samples, mt_node_samples, coverage_benchmarks, sequencer_lab, sequencer_lab_instruments, discovery_proposals, discovery_proposal, test_types, test_type_by_code, references_details, biosample_report, sample_report, biosample_studies, list_variants, get_variant, variants_by_haplogroup, export_metadata, export_variants, export_variants_gff, list_region_builds, regions_by_build, - reports_coverage, reports_ancestry, reports_haplogroups, + reports_coverage, reports_ancestry, reports_haplogroups, reports_str_markers, haplogroup_str_signature, haplogroup_age, str_predict, ), components(schemas( @@ -524,6 +531,7 @@ fn csv_field(s: &str) -> String { GenomeRegionDto, StudyDto, ExportMetadataDto, Page, Page, Page, FedCoverageByBuildDto, AncestryShareDto, HaplogroupCountDto, StrSignatureMarkerDto, StrPredictRequest, StrPredictionDto, StrPredictResponseDto, AgeEstimateDto, + StrMarkerStatDto, )), tags( (name = "tree", description = "Y/MT haplogroup trees"), @@ -559,6 +567,7 @@ pub fn router() -> Router { .route("/api/v1/reports/coverage", get(reports_coverage)) .route("/api/v1/reports/ancestry", get(reports_ancestry)) .route("/api/v1/reports/haplogroups", get(reports_haplogroups)) + .route("/api/v1/reports/str-markers", get(reports_str_markers)) .route("/api/v1/references/details", get(references_details)) .route("/api/v1/references/details/:publication_id/biosamples", get(biosample_report)) .route("/api/v1/samples/:slug", get(sample_report)) diff --git a/rust/crates/du-web/src/routes/mod.rs b/rust/crates/du-web/src/routes/mod.rs index 6bf657d..6c4978b 100644 --- a/rust/crates/du-web/src/routes/mod.rs +++ b/rust/crates/du-web/src/routes/mod.rs @@ -41,6 +41,7 @@ pub mod samples; pub mod sequencer; pub mod social; pub mod social_edge; +pub mod str_markers; pub mod tree; pub mod variants; pub mod versioning; @@ -65,6 +66,7 @@ pub fn app(state: AppState) -> Router { .merge(samples::router()) .merge(maps::router()) .merge(coverage::router()) + .merge(str_markers::router()) .merge(pages::router()) .merge(auth_routes::router()) .merge(curator::router()) diff --git a/rust/crates/du-web/src/routes/str_markers.rs b/rust/crates/du-web/src/routes/str_markers.rs new file mode 100644 index 0000000..791e466 --- /dev/null +++ b/rust/crates/du-web/src/routes/str_markers.rs @@ -0,0 +1,285 @@ +//! Public Y-STR marker report: the observed range (min / modal / max) of every +//! marker across the whole corpus, with its repeat motif and mutation rate where +//! we hold them. +//! +//! This is the horizontal view of the STR data. `du_db::ystr` otherwise only ever +//! reads it per-haplogroup (modal signatures, ASR, branch ages), so there was no +//! way to ask what the corpus says about a marker. The shape mirrors the MIN/MAX/ +//! MODE summary a vendor project page publishes, over a larger cohort. +//! +//! Two columns exist to keep the numbers honest rather than merely tidy: +//! null alleles (a reported `0`, which is a deleted locus, not a zero repeat +//! count) are excluded from the range and counted separately, and each marker +//! states how it enters the age model — most markers have no published rate and +//! are scored at `du_db::ystr::DEFAULT_STR_RATE`. + +use crate::error::AppError; +use crate::i18n::{Locale, T}; +use crate::render::html; +use crate::state::AppState; +use axum::extract::{Query, State}; +use axum::response::Response; +use axum::routing::get; +use axum::Router; +use serde::Deserialize; + +pub fn router() -> Router { + Router::new().route("/str-markers", get(page)) +} + +/// How a marker enters the STR age model, as a rendered badge. Kept as a small +/// enum-ish trio of flags because askama can't compare strings in the template. +struct AgeStatus { + label: String, + /// Bootstrap contextual class for the badge. + css: &'static str, +} + +/// One marker's row, fully pre-formatted — the template does no formatting or +/// comparison (this crate's convention: display logic lives in the route module). +struct MarkerRow { + marker: String, + multi_copy: bool, + motif: String, + period: String, + /// Min / modal / max, already rendered: repeat counts for simple markers, + /// copy vectors (e.g. `11-15`) for multi-copy ones. + min: String, + modal: String, + max: String, + observations: String, + samples: String, + distinct_values: String, + /// Compact caveat cell, e.g. `2 null · 1 partial`; empty when clean. + notes: String, + rate: String, + rate_ci: String, + rate_source: String, + age_status: AgeStatus, +} + +/// Corpus totals for the summary line above the table. +struct Totals { + markers: String, + observations: String, + with_motif: String, + with_rate: String, +} + +#[derive(Deserialize)] +struct MarkerQuery { + /// Case-insensitive substring filter on the marker name. + q: Option, + /// `observations` (default), `name`, `motif`, or `range`. + sort: Option, +} + +#[derive(askama::Template)] +#[template(path = "str/markers.html")] +struct MarkersTemplate { + t: T, + next: String, + user: Option, + rows: Vec, + totals: Totals, + q: String, + sort: String, + /// The rate used where a marker has none, shown in the legend. + default_rate: String, +} + +/// Integer with thousands separators (e.g. `12,345`). +fn fmt_count(n: i64) -> String { + let s = n.abs().to_string(); + let mut out = String::new(); + for (i, ch) in s.chars().enumerate() { + if i > 0 && (s.len() - i) % 3 == 0 { + out.push(','); + } + out.push(ch); + } + if n < 0 { + format!("-{out}") + } else { + out + } +} + +/// A mutation rate in the scale it is actually read at — per generation, and small +/// enough that plain decimal notation is unreadable below ~1e-4. +fn fmt_rate(v: Option) -> String { + match v { + Some(r) if r >= 0.001 => format!("{r:.5}"), + Some(r) if r > 0.0 => format!("{r:.2e}"), + _ => "—".into(), + } +} + +/// Pick the rendered value for one end of the range: simple markers carry a +/// repeat count, multi-copy markers a copy vector. Exactly one is ever set. +fn range_cell(value: Option, combination: &Option) -> String { + match (value, combination) { + (Some(v), _) => v.to_string(), + (None, Some(c)) => c.clone(), + _ => "—".into(), + } +} + +/// Null-allele and partial-repeat counts folded into one cell, so two columns +/// that are empty for nearly every marker don't widen the table. +fn notes(null_alleles: i64, complex: i64, t: &T) -> String { + let mut parts = Vec::new(); + if null_alleles > 0 { + parts.push(format!("{null_alleles} {}", t.get("str.markers.note.null"))); + } + if complex > 0 { + parts.push(format!("{complex} {}", t.get("str.markers.note.partial"))); + } + parts.join(" · ") +} + +fn age_status(status: &str, t: &T) -> AgeStatus { + match status { + "MEASURED_RATE" => { + AgeStatus { label: t.get("str.markers.status.measured").to_string(), css: "text-bg-success" } + } + "EXCLUDED" => { + AgeStatus { label: t.get("str.markers.status.excluded").to_string(), css: "text-bg-secondary" } + } + _ => AgeStatus { label: t.get("str.markers.status.default").to_string(), css: "text-bg-warning" }, + } +} + +fn to_row(m: du_db::ystr::MarkerStat, t: &T) -> MarkerRow { + let ci = match (m.rate_ci_low, m.rate_ci_high) { + (Some(lo), Some(hi)) => format!("{}–{}", fmt_rate(Some(lo)), fmt_rate(Some(hi))), + _ => "—".into(), + }; + MarkerRow { + marker: m.marker_name, + multi_copy: m.multi_copy, + motif: m.motif.unwrap_or_else(|| "—".into()), + period: m.period.map(|p| p.to_string()).unwrap_or_else(|| "—".into()), + min: range_cell(m.min_value, &m.min_combination), + modal: range_cell(m.modal_value, &m.modal_combination), + max: range_cell(m.max_value, &m.max_combination), + observations: fmt_count(m.observations), + samples: fmt_count(m.samples), + distinct_values: fmt_count(m.distinct_values), + notes: notes(m.null_alleles, m.complex_count, t), + rate: fmt_rate(m.mutation_rate), + rate_ci: ci, + rate_source: m.rate_source.unwrap_or_default(), + age_status: age_status(&m.age_model_status, t), + } +} + +/// Sort key for a marker name: `DYS` numerically (so DYS19 precedes +/// DYS385 rather than following DYS1xx alphabetically), everything else after it +/// by name. +fn name_key(marker: &str) -> (u8, u32, String) { + if let Some(rest) = marker.strip_prefix("DYS") { + let digits: String = rest.chars().take_while(|c| c.is_ascii_digit()).collect(); + if !digits.is_empty() { + if let Ok(n) = digits.parse::() { + // Suffixed variants (DYS389I/II) sort after the bare number. + return (0, n, rest[digits.len()..].to_string()); + } + } + } + (1, 0, marker.to_string()) +} + +async fn page( + State(st): State, + locale: Locale, + user: crate::auth::MaybeUser, + Query(query): Query, +) -> Result { + let all = du_db::ystr::marker_stats(&st.pool).await?; + + // Totals describe the whole corpus, not the filtered view — they are the + // reference-coverage headline (how much of what we observe we can explain). + let totals = Totals { + markers: fmt_count(all.len() as i64), + observations: fmt_count(all.iter().map(|m| m.observations).sum()), + with_motif: fmt_count(all.iter().filter(|m| m.motif.is_some()).count() as i64), + with_rate: fmt_count( + all.iter().filter(|m| m.age_model_status == "MEASURED_RATE").count() as i64 + ), + }; + + let q = query.q.unwrap_or_default(); + let needle = q.trim().to_ascii_lowercase(); + let mut rows: Vec = all + .into_iter() + .filter(|m| needle.is_empty() || m.marker_name.to_ascii_lowercase().contains(&needle)) + .collect(); + + // SQL already returns observations-descending; re-sort only for other keys. + let sort = query.sort.unwrap_or_default(); + match sort.as_str() { + "name" => rows.sort_by(|a, b| name_key(&a.marker_name).cmp(&name_key(&b.marker_name))), + // Motif-less markers last, so the sort surfaces what we can explain. + "motif" => rows.sort_by(|a, b| { + (a.motif.is_none(), a.period, name_key(&a.marker_name)).cmp(&( + b.motif.is_none(), + b.period, + name_key(&b.marker_name), + )) + }), + // Widest observed spread first — the most variable markers in the corpus. + "range" => rows.sort_by(|a, b| b.distinct_values.cmp(&a.distinct_values)), + _ => {} + } + + let t = &locale.t; + let rows = rows.into_iter().map(|m| to_row(m, t)).collect(); + + Ok(html(&MarkersTemplate { + t: locale.t, + next: locale.next, + user: user.nav(), + rows, + totals, + q, + sort, + default_rate: fmt_rate(Some(du_db::ystr::DEFAULT_STR_RATE)), + })) +} + +#[cfg(test)] +mod tests { + use super::{fmt_count, fmt_rate, name_key, range_cell}; + + #[test] + fn markers_sort_numerically_within_dys() { + let mut got = vec!["DYS448", "DYS19", "CDY", "DYS389II", "DYS389I", "DYS385", "Y-GATA-H4"]; + got.sort_by_key(|m| name_key(m)); + // DYS19 before DYS385 (numeric, not lexical); non-DYS names after. + assert_eq!( + got, + vec!["DYS19", "DYS385", "DYS389I", "DYS389II", "DYS448", "CDY", "Y-GATA-H4"] + ); + } + + #[test] + fn range_cell_picks_whichever_shape_the_marker_has() { + assert_eq!(range_cell(Some(13), &None), "13"); + assert_eq!(range_cell(None, &Some("11-15".to_string())), "11-15"); + assert_eq!(range_cell(None, &None), "—"); + } + + #[test] + fn rates_stay_readable_at_both_scales() { + assert_eq!(fmt_rate(Some(0.00278)), "0.00278"); + assert_eq!(fmt_rate(Some(9.74635e-05)), "9.75e-5"); + assert_eq!(fmt_rate(None), "—"); + } + + #[test] + fn counts_get_thousands_separators() { + assert_eq!(fmt_count(875), "875"); + assert_eq!(fmt_count(508227), "508,227"); + } +} diff --git a/rust/crates/du-web/templates/base.html b/rust/crates/du-web/templates/base.html index f5bbfd2..7c64fb6 100644 --- a/rust/crates/du-web/templates/base.html +++ b/rust/crates/du-web/templates/base.html @@ -34,6 +34,7 @@ diff --git a/rust/crates/du-web/templates/static/page.html b/rust/crates/du-web/templates/static/page.html index adccee6..0e839ea 100644 --- a/rust/crates/du-web/templates/static/page.html +++ b/rust/crates/du-web/templates/static/page.html @@ -17,6 +17,7 @@

{{ title }}

  • {{ t.get("nav.variants") }}
  • {{ t.get("nav.references") }}
  • {{ t.get("nav.coverage") }}
    • +
  • {{ t.get("nav.strMarkers") }}
    • {% else if page == "faq" %} diff --git a/rust/crates/du-web/templates/str/markers.html b/rust/crates/du-web/templates/str/markers.html new file mode 100644 index 0000000..f23af9f --- /dev/null +++ b/rust/crates/du-web/templates/str/markers.html @@ -0,0 +1,87 @@ +{% extends "base.html" %} +{% block title %}{{ t.get("str.markers.title") }} — {{ t.get("app.name") }}{% endblock %} +{% block content %} +

    {{ t.get("str.markers.title") }}

    +

    {{ t.get("str.markers.intro") }}

    + +

    + {{ t.get("str.markers.total.markers") }}: {{ totals.markers }} + {{ t.get("str.markers.total.observations") }}: {{ totals.observations }} + {{ t.get("str.markers.total.withMotif") }}: {{ totals.with_motif }} + {{ t.get("str.markers.total.withRate") }}: {{ totals.with_rate }} +

    + +
    +
    + + +
    +
    + + +
    +
    + +
    +
    + +
    + + + + + + + + + + + + + + + + + + + + + + + + + {% for r in rows %} + + + + + + + + + + + + + + + + {% else %} + + {% endfor %} + +
    {{ t.get("str.markers.col.marker") }}{{ t.get("str.markers.col.repeat") }}{{ t.get("str.markers.col.observedRange") }}{{ t.get("str.markers.col.observations") }}{{ t.get("str.markers.col.mutationRate") }}{{ t.get("str.markers.col.ageModel") }}
    {{ t.get("str.markers.col.motif") }}{{ t.get("str.markers.col.period") }}{{ t.get("str.markers.col.min") }}{{ t.get("str.markers.col.modal") }}{{ t.get("str.markers.col.max") }}{{ t.get("str.markers.col.total") }}{{ t.get("str.markers.col.samples") }}{{ t.get("str.markers.col.distinct") }}{{ t.get("str.markers.col.notes") }}{{ t.get("str.markers.col.rate") }}{{ t.get("str.markers.col.ci") }}
    + {{ r.marker }} + {% if r.multi_copy %}{{ t.get("str.markers.multiCopy") }}{% endif %} + {{ r.motif }}{{ r.period }}{{ r.min }}{{ r.modal }}{{ r.max }}{{ r.observations }}{{ r.samples }}{{ r.distinct_values }}{{ r.notes }}{{ r.rate }}{{ r.rate_ci }}{{ r.age_status.label }}
    {{ t.get("str.markers.none") }}
    +
    + +

    {{ t.get("str.markers.legend.null") }}

    +

    {{ t.get("str.markers.legend.default") }} {{ default_rate }}.

    +

    {{ t.get("str.markers.legend.motif") }}

    +{% endblock %} diff --git a/rust/locales/en.txt b/rust/locales/en.txt index a682d44..e86dbf5 100644 --- a/rust/locales/en.txt +++ b/rust/locales/en.txt @@ -7,6 +7,7 @@ nav.variants=Variants nav.references=References nav.map=Map nav.coverage=Coverage +nav.strMarkers=Y-STR Markers nav.curator=Curator nav.tools=Tools nav.profile=Profile @@ -764,3 +765,45 @@ dc.col.magnitude=Magnitude dc.col.home=Home node dc.col.foreign=Foreign in dc.col.away=Members away + +str.markers.title=Y-STR Markers +str.markers.intro=Observed repeat range for every Y-STR marker across all profiles we hold — the minimum, modal and maximum value reported per marker, with its repeat motif and mutation rate where known. +str.markers.total.markers=Markers +str.markers.total.observations=Observations +str.markers.total.withMotif=With known motif +str.markers.total.withRate=With measured rate +str.markers.filter.marker=Marker +str.markers.filter.sort=Sort by +str.markers.filter.apply=Apply +str.markers.sort.observations=Most observed +str.markers.sort.name=Marker name +str.markers.sort.motif=Motif +str.markers.sort.range=Most variable +str.markers.col.marker=Marker +str.markers.col.repeat=Repeat unit +str.markers.col.motif=Motif +str.markers.col.period=Length +str.markers.col.observedRange=Observed range +str.markers.col.min=Min +str.markers.col.modal=Modal +str.markers.col.max=Max +str.markers.col.observations=Observations +str.markers.col.total=Total +str.markers.col.samples=Samples +str.markers.col.distinct=Distinct +str.markers.col.notes=Notes +str.markers.col.mutationRate=Mutation rate +str.markers.col.rate=Per generation +str.markers.col.ci=95% CI +str.markers.col.ageModel=Age model +str.markers.multiCopy=multi-copy +str.markers.multiCopy.help=Palindromic or duplicated locus, reported as a copy vector (e.g. 11-15) and not used for age estimation. +str.markers.note.null=null +str.markers.note.partial=partial +str.markers.status.measured=Measured rate +str.markers.status.default=Default rate +str.markers.status.excluded=Not scored +str.markers.none=No STR markers match. +str.markers.legend.null=A reported value of 0 is a null allele (a deleted locus), not a zero repeat count, so it is excluded from the range and counted under Notes. Partial repeats (e.g. 10.2) are preserved but not scored. +str.markers.legend.default=Markers without a published mutation rate are scored in the branch-age model at the default rate of +str.markers.legend.motif=Motifs and repeat lengths are shown only where a published source gives one; most extended markers have none. diff --git a/rust/locales/es.txt b/rust/locales/es.txt index 1a12bfe..cc7c1d1 100644 --- a/rust/locales/es.txt +++ b/rust/locales/es.txt @@ -7,6 +7,7 @@ nav.variants=Variantes nav.references=Referencias nav.map=Mapa nav.coverage=Cobertura +nav.strMarkers=Marcadores Y-STR nav.curator=Curador nav.tools=Herramientas nav.profile=Perfil @@ -758,3 +759,45 @@ announce.post=Publicar anuncio announce.recent=Anuncios recientes announce.none=Aún no hay anuncios. announce.replies=respuestas + +str.markers.title=Marcadores Y-STR +str.markers.intro=Rango de repeticiones observado para cada marcador Y-STR en todos los perfiles que conservamos: el valor mínimo, modal y máximo notificado por marcador, con su motivo de repetición y su tasa de mutación cuando se conocen. +str.markers.total.markers=Marcadores +str.markers.total.observations=Observaciones +str.markers.total.withMotif=Con motivo conocido +str.markers.total.withRate=Con tasa medida +str.markers.filter.marker=Marcador +str.markers.filter.sort=Ordenar por +str.markers.filter.apply=Aplicar +str.markers.sort.observations=Más observados +str.markers.sort.name=Nombre del marcador +str.markers.sort.motif=Motivo +str.markers.sort.range=Más variables +str.markers.col.marker=Marcador +str.markers.col.repeat=Unidad de repetición +str.markers.col.motif=Motivo +str.markers.col.period=Longitud +str.markers.col.observedRange=Rango observado +str.markers.col.min=Mín +str.markers.col.modal=Modal +str.markers.col.max=Máx +str.markers.col.observations=Observaciones +str.markers.col.total=Total +str.markers.col.samples=Muestras +str.markers.col.distinct=Distintos +str.markers.col.notes=Notas +str.markers.col.mutationRate=Tasa de mutación +str.markers.col.rate=Por generación +str.markers.col.ci=IC 95% +str.markers.col.ageModel=Modelo de edad +str.markers.multiCopy=multicopia +str.markers.multiCopy.help=Locus palindrómico o duplicado, notificado como vector de copias (p. ej. 11-15) y no utilizado para estimar la edad. +str.markers.note.null=nulo +str.markers.note.partial=parcial +str.markers.status.measured=Tasa medida +str.markers.status.default=Tasa por defecto +str.markers.status.excluded=No puntuado +str.markers.none=Ningún marcador STR coincide. +str.markers.legend.null=Un valor de 0 es un alelo nulo (un locus eliminado), no un recuento de cero repeticiones, por lo que se excluye del rango y se cuenta en Notas. Las repeticiones parciales (p. ej. 10,2) se conservan pero no se puntúan. +str.markers.legend.default=Los marcadores sin tasa de mutación publicada se puntúan en el modelo de edad de rama con la tasa por defecto de +str.markers.legend.motif=Los motivos y las longitudes de repetición solo se muestran cuando una fuente publicada los proporciona; la mayoría de los marcadores extendidos carecen de ellos. diff --git a/rust/locales/fr.txt b/rust/locales/fr.txt index 65e09dd..e7a1bba 100644 --- a/rust/locales/fr.txt +++ b/rust/locales/fr.txt @@ -7,6 +7,7 @@ nav.variants=Variants nav.references=Références nav.map=Carte nav.coverage=Couverture +nav.strMarkers=Marqueurs Y-STR nav.curator=Curateur nav.tools=Outils nav.profile=Profil @@ -758,3 +759,45 @@ announce.post=Publier l'annonce announce.recent=Annonces récentes announce.none=Aucune annonce pour l'instant. announce.replies=réponses + +str.markers.title=Marqueurs Y-STR +str.markers.intro=Plage de répétitions observée pour chaque marqueur Y-STR dans l'ensemble des profils que nous détenons : valeur minimale, modale et maximale rapportée par marqueur, avec son motif de répétition et son taux de mutation lorsqu'ils sont connus. +str.markers.total.markers=Marqueurs +str.markers.total.observations=Observations +str.markers.total.withMotif=Avec motif connu +str.markers.total.withRate=Avec taux mesuré +str.markers.filter.marker=Marqueur +str.markers.filter.sort=Trier par +str.markers.filter.apply=Appliquer +str.markers.sort.observations=Les plus observés +str.markers.sort.name=Nom du marqueur +str.markers.sort.motif=Motif +str.markers.sort.range=Les plus variables +str.markers.col.marker=Marqueur +str.markers.col.repeat=Unité de répétition +str.markers.col.motif=Motif +str.markers.col.period=Longueur +str.markers.col.observedRange=Plage observée +str.markers.col.min=Min +str.markers.col.modal=Modale +str.markers.col.max=Max +str.markers.col.observations=Observations +str.markers.col.total=Total +str.markers.col.samples=Échantillons +str.markers.col.distinct=Distincts +str.markers.col.notes=Notes +str.markers.col.mutationRate=Taux de mutation +str.markers.col.rate=Par génération +str.markers.col.ci=IC 95 % +str.markers.col.ageModel=Modèle d'âge +str.markers.multiCopy=multicopie +str.markers.multiCopy.help=Locus palindromique ou dupliqué, rapporté sous forme de vecteur de copies (p. ex. 11-15) et non utilisé pour l'estimation de l'âge. +str.markers.note.null=nul +str.markers.note.partial=partiel +str.markers.status.measured=Taux mesuré +str.markers.status.default=Taux par défaut +str.markers.status.excluded=Non évalué +str.markers.none=Aucun marqueur STR correspondant. +str.markers.legend.null=Une valeur de 0 correspond à un allèle nul (locus supprimé), et non à un compte de zéro répétition ; elle est donc exclue de la plage et comptée dans les Notes. Les répétitions partielles (p. ex. 10,2) sont conservées mais non évaluées. +str.markers.legend.default=Les marqueurs sans taux de mutation publié sont évalués dans le modèle d'âge de branche au taux par défaut de +str.markers.legend.motif=Les motifs et longueurs de répétition ne sont affichés que lorsqu'une source publiée les fournit ; la plupart des marqueurs étendus n'en ont pas. diff --git a/rust/migrations/0069_str_marker.sql b/rust/migrations/0069_str_marker.sql new file mode 100644 index 0000000..596ca5c --- /dev/null +++ b/rust/migrations/0069_str_marker.sql @@ -0,0 +1,96 @@ +-- Y-STR marker reference: the repeat unit (motif), its period, and eventually the +-- locus coordinates for each marker name we observe in profiles. +-- +-- Until now marker identity was a bare `marker_name TEXT` everywhere — in +-- genomics.biosample_str_profile.markers, fed.str_profile.markers, +-- tree.haplogroup_ancestral_str, tree.haplogroup_str_asr and +-- genomics.str_mutation_rate. Nothing said what DYS390 *is*. This table is that +-- missing dimension, and the join target for the corpus-wide marker report +-- (du_db::ystr::marker_stats → /str-markers, /api/v1/reports/str-markers). +-- +-- `coordinates` deliberately mirrors the core.genome_region shape +-- ({GRCh38:{contig,start,end}, hs1:{...}}) so the existing coordinate-lift +-- machinery can fill it in place once we have a locus source. It is empty today: +-- no table in this database has ever held an STR locus position. +-- +-- `aliases` is not cosmetic. The profiles say DYS19 and Y-GATA-H4; the mutation-rate +-- table (seeded from Willems 2016 / YHRD) says DYS394 and YGATAH4 for the same two +-- markers. A join on name alone silently drops the two most recognizable markers in +-- the panel, so every lookup folds through aliases. Canonical `marker_name` is +-- always the form that appears in profiles — that is what a tester sees on a report. +-- +-- Seeded from the only motif information this database holds: the free-text +-- `source` strings on genomics.str_mutation_rate, e.g. +-- 'Willems et al. 2016 (1000G MUTEA); motif AGAT, n=468' +-- which yields a motif for 126 of the 137 rate rows. That is ~132 markers against +-- the 829 we actually observe: the extended Big Y markers (DYS425, DYF395S1, +-- DYS520, ...) have no published motif in any source we currently hold, and the +-- report renders them blank rather than guessing. Idempotent (corrective re-run). + +CREATE TABLE genomics.str_marker ( + marker_name TEXT PRIMARY KEY, -- canonical name, as it appears in profiles + motif TEXT, -- repeat unit, e.g. AGAT + period SMALLINT, -- length(motif) + coordinates JSONB NOT NULL DEFAULT '{}'::jsonb, -- {GRCh38:{...}, hs1:{...}}; empty today + panel_names TEXT[], + aliases TEXT[], -- alternate names used by other sources + multi_copy BOOLEAN NOT NULL DEFAULT false, -- palindromic/duplicated; unscored by the age model + source TEXT, + created_at TIMESTAMPTZ NOT NULL DEFAULT now(), + updated_at TIMESTAMPTZ NOT NULL DEFAULT now() +); +CREATE INDEX str_marker_aliases_gin ON genomics.str_marker USING gin (aliases); + +WITH alias(canonical, alt) AS ( + -- canonical (the form profiles use) ← alternate (rate-table / literature form) + VALUES ('DYS19', 'DYS394'), -- same locus; Willems/YHRD use the DYS394 designation + ('Y-GATA-H4', 'YGATAH4'), + ('Y-GATA-A10', 'YGATAA10'), + ('DYS389I', 'DYS389i'), + ('DYS389II', 'DYS389ii') +), canon AS ( + SELECT COALESCE(a.canonical, r.marker_name) AS marker_name, + (regexp_match(r.source, 'motif ([ACGT]+)'))[1] AS motif, + r.panel_names, + r.source + FROM genomics.str_mutation_rate r + LEFT JOIN alias a ON a.alt = r.marker_name +), folded AS ( + -- The seed emitted DYS19 and DYS394 as separate rows for one locus (identical + -- rate and motif), so both fold to DYS19 here. Prefer the row that carries a + -- motif; ties break on name for determinism. + SELECT DISTINCT ON (marker_name) * + FROM canon + ORDER BY marker_name, (motif IS NULL), source +) +INSERT INTO genomics.str_marker (marker_name, motif, period, panel_names, aliases, multi_copy, source) +SELECT f.marker_name, + f.motif, + length(f.motif)::smallint, + f.panel_names, + NULLIF(ARRAY(SELECT alt FROM alias WHERE canonical = f.marker_name), '{}'), + -- Multi-copy (palindromic/duplicated) markers: the 7 observed as + -- {type:multiCopy} in our profiles, plus DYF387S1, which is multi-copy but + -- absent from the FTDNA DYS panel export. du_db::ystr scores none of them. + f.marker_name IN ('DYS385','DYS464','CDY','YCAII','DYS459','DYF395S1','DYS413','DYF387S1'), + 'derived from genomics.str_mutation_rate (Willems 2016 1000G MUTEA / YHRD)' + FROM folded f +ON CONFLICT (marker_name) DO UPDATE SET + motif = COALESCE(EXCLUDED.motif, genomics.str_marker.motif), + period = COALESCE(EXCLUDED.period, genomics.str_marker.period), + panel_names = EXCLUDED.panel_names, + aliases = EXCLUDED.aliases, + multi_copy = EXCLUDED.multi_copy, + source = EXCLUDED.source, + updated_at = now(); + +-- The multi-copy markers absent from the rate table (they are unscored, so they +-- were never given one) still belong in the registry — otherwise the flag is true +-- for DYS385 and silently false for the six equally palindromic loci beside it. +-- No motif: a duplicated locus has one per copy, and we hold none of them. +INSERT INTO genomics.str_marker (marker_name, multi_copy, source) +SELECT m, true, 'multi-copy locus; observed as {type:multiCopy} in vendor profiles' + FROM unnest(ARRAY['DYS464','CDY','YCAII','DYS459','DYF395S1','DYS413']) AS m +ON CONFLICT (marker_name) DO UPDATE SET + multi_copy = true, + updated_at = now(); diff --git a/rust/migrations/0070_str_rate_provenance.sql b/rust/migrations/0070_str_rate_provenance.sql new file mode 100644 index 0000000..a80a0cf --- /dev/null +++ b/rust/migrations/0070_str_rate_provenance.sql @@ -0,0 +1,57 @@ +-- Mutation-rate provenance: which rate produced which branch age, and room for +-- rates we derive ourselves. +-- +-- Two gaps this closes. +-- +-- (1) genomics.str_mutation_rate was one row per marker (marker_name UNIQUE), so a +-- rate we derive from our own tree could only land by overwriting the published +-- Willems/YHRD value — destroying the comparison that makes the derived number +-- trustworthy. Re-keyed on (marker_name, method) so PUBLISHED and DERIVED +-- coexist, with a partial unique index guaranteeing exactly one *active* rate +-- per marker. du_db::ystr::load_marker_models reads `WHERE is_active`, which +-- makes promoting a derived rate a data operation rather than a migration. +-- +-- (2) Nothing recorded which rates an age estimate actually used. Both +-- compute_str_age and propagate_str silently fall back to +-- ystr::DEFAULT_STR_RATE (0.0025) for any marker without a rate row — and only +-- 114 of the 829 markers we observe have one, so most of the STR age signal +-- rests on that placeholder. tree.haplogroup_age_estimate now carries the +-- split, so a branch age states how much of it is measured. +-- +-- The observation columns (mutations_observed / meioses_observed) exist so a +-- derived rate carries its own evidence and CI, and so `is_active` can stay false +-- until a marker has accumulated enough meioses to be worth trusting. The +-- derivation job itself is a follow-up; see docs — the estimator must take branch +-- lengths from method='SNP_POISSON' ages ONLY, never COMBINED or STR_VARIANCE, +-- which would feed STR-derived ages back into STR rate estimation. + +ALTER TABLE genomics.str_mutation_rate + ADD COLUMN method TEXT NOT NULL DEFAULT 'PUBLISHED', -- PUBLISHED | DERIVED + ADD COLUMN mutations_observed NUMERIC, -- fractional: multi-step branches contribute >1 + ADD COLUMN meioses_observed NUMERIC, -- Σ branch length in generations + ADD COLUMN derived_at TIMESTAMPTZ, + ADD COLUMN tree_revision BIGINT, -- tree.tree_revision the derivation ran against + ADD COLUMN is_active BOOLEAN NOT NULL DEFAULT true; + +-- The existing 137 rows are all PUBLISHED and active by virtue of the defaults. +ALTER TABLE genomics.str_mutation_rate + DROP CONSTRAINT str_mutation_rate_marker_name_key; +CREATE UNIQUE INDEX str_mutation_rate_marker_method_key + ON genomics.str_mutation_rate (marker_name, method); +-- At most one rate per marker may feed the age model. +CREATE UNIQUE INDEX str_mutation_rate_active_key + ON genomics.str_mutation_rate (marker_name) WHERE is_active; + +COMMENT ON COLUMN genomics.str_mutation_rate.method IS + 'PUBLISHED (literature: Willems 2016 / YHRD) or DERIVED (estimated from our own tree)'; +COMMENT ON COLUMN genomics.str_mutation_rate.is_active IS + 'The single rate per marker that du_db::ystr::load_marker_models feeds to the age model'; + +-- Which rates produced this age estimate. Written by ystr::recompute_signatures. +ALTER TABLE tree.haplogroup_age_estimate + ADD COLUMN rate_method TEXT, -- PUBLISHED / DERIVED / MIXED + ADD COLUMN measured_rate_markers INTEGER, -- markers scored with a real rate row + ADD COLUMN default_rate_markers INTEGER; -- markers that fell back to DEFAULT_STR_RATE + +COMMENT ON COLUMN tree.haplogroup_age_estimate.default_rate_markers IS + 'Markers scored at ystr::DEFAULT_STR_RATE for want of a rate row — high values mean the estimate is weakly grounded'; diff --git a/rust/scripts/seed-str-mutation-rates.sql b/rust/scripts/seed-str-mutation-rates.sql index ce1ae08..a14db51 100644 --- a/rust/scripts/seed-str-mutation-rates.sql +++ b/rust/scripts/seed-str-mutation-rates.sql @@ -165,7 +165,9 @@ VALUES ('DYS726', 0.0003711496, 0.0002828180, 0.0004870696, ARRAY['Willems2016-1kG'], 'Willems et al. 2016 (1000G MUTEA); motif AAGG, n=668'), ('Y-GATA-A10', 0.0035312625, 0.0030108908, 0.0041415699, ARRAY['Willems2016-1kG'], 'Willems et al. 2016 (1000G MUTEA); motif AGAT, n=847'), ('YGATAH4', 0.0025500000, 0.0017471797, 0.0033455338, ARRAY['YHRD'], 'YHRD combined (39 mutations / 15316 meioses)') -ON CONFLICT (marker_name) DO UPDATE SET +-- Keyed on (marker_name, method) since mig 0070 — these are all PUBLISHED rates and +-- must not disturb any DERIVED row estimated from our own tree. +ON CONFLICT (marker_name, method) DO UPDATE SET mutation_rate = EXCLUDED.mutation_rate, mutation_rate_lower = EXCLUDED.mutation_rate_lower, mutation_rate_upper = EXCLUDED.mutation_rate_upper, From 8d065a819326984f1097596d0e88440c85094d50 Mon Sep 17 00:00:00 2001 From: James Kane Date: Wed, 29 Jul 2026 05:37:21 -0500 Subject: [PATCH 3/4] fix(str): use sort_by_key for the marker-report sorts clippy::unnecessary_sort_by. The name and range arms rewrite directly; the motif arm follows for consistency, since a plain tuple key reads better than a hand-written comparison against a reconstructed tuple. Co-Authored-By: Claude Opus 5 (1M context) --- rust/crates/du-web/src/routes/str_markers.rs | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/rust/crates/du-web/src/routes/str_markers.rs b/rust/crates/du-web/src/routes/str_markers.rs index 791e466..3af0822 100644 --- a/rust/crates/du-web/src/routes/str_markers.rs +++ b/rust/crates/du-web/src/routes/str_markers.rs @@ -219,17 +219,13 @@ async fn page( // SQL already returns observations-descending; re-sort only for other keys. let sort = query.sort.unwrap_or_default(); match sort.as_str() { - "name" => rows.sort_by(|a, b| name_key(&a.marker_name).cmp(&name_key(&b.marker_name))), + "name" => rows.sort_by_key(|m| name_key(&m.marker_name)), // Motif-less markers last, so the sort surfaces what we can explain. - "motif" => rows.sort_by(|a, b| { - (a.motif.is_none(), a.period, name_key(&a.marker_name)).cmp(&( - b.motif.is_none(), - b.period, - name_key(&b.marker_name), - )) - }), + "motif" => { + rows.sort_by_key(|m| (m.motif.is_none(), m.period, name_key(&m.marker_name))) + } // Widest observed spread first — the most variable markers in the corpus. - "range" => rows.sort_by(|a, b| b.distinct_values.cmp(&a.distinct_values)), + "range" => rows.sort_by_key(|m| std::cmp::Reverse(m.distinct_values)), _ => {} } From 197d84f903c94038598e9bcb9d0a54b20ce7fe02 Mon Sep 17 00:00:00 2001 From: James Kane Date: Wed, 29 Jul 2026 06:04:05 -0500 Subject: [PATCH 4/4] fix(test): isolate deletes_split_covering_duplicate on an ephemeral DB MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The test connected to the shared DATABASE_URL database directly, unlike the two tests beside it in the same file, so it assumed that database was already migrated. That holds on a dev box and does not in CI, where the Postgres service starts empty — the insert failed with `relation "tree.haplogroup" does not exist`. Main has been red on this since 97f0853; it is unrelated to the STR work in this branch but blocks it, so it is fixed here. Switching to du_db::testing::ephemeral_db matches its siblings and makes the trailing manual cleanup redundant — the ephemeral database drops itself. Verified by pointing DATABASE_URL at a freshly created empty database and running the suite in parallel (CI's mode, not the --test-threads=1 that had been masking this locally alongside an already-migrated dev DB): 3/3 in variant_naming, and du-db green throughout. Co-Authored-By: Claude Opus 5 (1M context) --- rust/crates/du-db/tests/variant_naming.rs | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/rust/crates/du-db/tests/variant_naming.rs b/rust/crates/du-db/tests/variant_naming.rs index 2902456..8665c84 100644 --- a/rust/crates/du-db/tests/variant_naming.rs +++ b/rust/crates/du-db/tests/variant_naming.rs @@ -310,7 +310,11 @@ async fn variant_naming_authority_flow() { #[tokio::test] async fn deletes_split_covering_duplicate() { let Some(url) = database_url() else { return }; - let pool = PgPool::connect(&url).await.expect("connect"); + // Isolated + freshly migrated, like the two tests above. Connecting to the + // shared DATABASE_URL database directly assumes it is already migrated, which + // holds for a dev box but not for a CI Postgres service that starts empty. + let db = du_db::testing::ephemeral_db(&url).await.expect("ephemeral db"); + let pool = db.pool().clone(); let branch = mk_hg(&pool, "R-SPLITTEST").await; @@ -360,10 +364,5 @@ async fn deletes_split_covering_duplicate() { "a really-named row is never an erroneous duplicate" ); - // cleanup - for id in [same_site, on_branch] { - let _ = sqlx::query("DELETE FROM tree.haplogroup_variant WHERE variant_id = $1").bind(id).execute(&pool).await; - let _ = sqlx::query("DELETE FROM core.variant WHERE id = $1").bind(id).execute(&pool).await; - } - let _ = sqlx::query("DELETE FROM tree.haplogroup WHERE id = $1").bind(branch).execute(&pool).await; + // No cleanup: the ephemeral database drops itself. }