CNV-robust phasing, DBSCAN strategy, and expansion-sensitivity QC flags#17
Merged
Conversation
The minimum cluster size required for a read cluster to be treated as a haplotype in unphased mode was computed purely as a fraction of the total read count (min_haplotype_fraction). Under copy-number gain the inflated read count raises this fraction-based threshold, so a genuine minority/ expanded allele supported by enough reads gets discarded as "too small". find_roots then descends into the majority cluster and over-splits it into two near-identical haplotypes, and find_outliers dumps the real expansion as length outliers (false-negative genotype). Cap the fraction-based threshold by the absolute --support count so that an allele backed by at least `support` reads can never be swallowed by high coverage. Behaviour is unchanged whenever fraction*n <= support (all normal coverage, incl. every existing test); it only relaxes the bar in the high- coverage/CNV regime that caused the bug. Regression test added from chr15:34419425-34419451 (GOLGA8A) where a length-variable expansion on 5/44 reads was lost. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Read downsampling for consensus used an unseeded rand::rng(), so any haplotype with more reads than --consensus-reads (default 20) drew a different random subset each run, making the consensus length (and thus the genotype) non-reproducible run-to-run. This dominated cohort A/B diffs, especially at high-variance expansion loci. Seed a StdRng with a fixed constant per consensus() call: deterministic and independent of locus/thread ordering. Not exposed as a CLI flag — the subset choice is purely a performance measure and shouldn't affect reproducibility. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Adds an opt-in alternative to the default Ward strategy for splitting unphased
reads into haplotypes. Instead of a length-dominated Levenshtein distance (which
fragments length-variable expansions and can lose them to the outlier bin), it
builds length-invariant k-mer composition feature vectors and clusters them with
DBSCAN, so reads of the same motif group together regardless of repeat length.
- src/features.rs: canonical k-mer feature vectors + tandem-period detection
(lifted from the trout cohort tool).
- src/dbscan.rs: custom DBSCAN over CSR adjacency returning cluster labels.
- split_dbscan(): min_samples = --support; the two largest clusters become the
haplotypes, remaining clusters + noise become outliers. Reports n_clusters.
- New INFO/NCLUSTERS field (DBSCAN mode only) so complex multi-population loci
(>2 clusters) can be flagged downstream.
- CLI: --phasing-strategy {ward,dbscan} (default ward), --dbscan-eps (0.2),
--dbscan-length-weight (1.0).
Default behaviour is unchanged (Ward remains the default). On the GOLGA8A
regression locus the DBSCAN path recovers the expansion as a distinct
haplotype (hap2 median ~606bp) instead of dropping it as outliers.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Change the startup info line from "Collected arguments" to log the entire parsed Cli struct, so each run's log records exactly which parameters were used (phasing strategy, dbscan eps / length-weight, support, etc.). This makes cohort run comparisons self-documenting. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
trout's defaults (eps=0.2, length_weight=1.0) were tuned for clean per-sample consensus points; on STRdust's raw, error-containing reads they over-fragment, dumping ~5 reads/locus to the outlier bin regardless of biology and losing real expansions. On the GOLGA8A cohort, eps=0.4/length_weight=0.3 restores discrimination (outlier rate 2.1 in expansion vs 1.3 in reference loci; flag rate 47% -> 10%) and fixes the catastrophic drops (e.g. UBC66_1 44 -> 1114, NHB18_068 27 -> 871). Updated the regression test to the new defaults. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Add the unphased phasing controls (--phasing-strategy, --dbscan-eps, --dbscan-length-weight, --min-haplotype-fraction) to the arguments list, a Notes section explaining the ward vs dbscan strategies and their trade-offs, and the NCLUSTERS INFO field to the output format example. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Pure formatting (line-wrap) of the StdRng::seed_from_u64 chain introduced by the deterministic-downsampling fix; no behaviour change. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Adds three sensitivity/QC outputs for --unphased calling, all on by default: - EXPANSION_OUTLIER: >=2 reads more than 2x the longer called allele, catching large expansions missed by both alleles (dropped to noise/outliers or removed as length outliers during consensus). - IMPRECISE_LENGTH: called-allele read-length coefficient of variation > 0.2, flagging continuous/long-tailed loci where a single consensus length is not representative. - MRL FORMAT field: median read length per allele relative to reference, robust to a long length tail; compare against RB. Renames --phasing-strategy to --phasing and adds a 'both' QC mode that reports the Ward call but also runs DBSCAN and flags DISCORDANT_LENGTH (+ DBSCAN_RB) when the two differ by more than 2x on the longer allele. The thresholds (EXPANSION_OUTLIER, IMPRECISE_LENGTH, DISCORDANT_LENGTH) and the DBSCAN eps/length-weight are kept as documented constants rather than CLI args; README documents each with tuning intuition. Breaking: removes --phasing-strategy (now --phasing), --dbscan-eps and --dbscan-length-weight. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
This branch improves STRdust's sensitivity to real expansions on
--unphaseddata and adds QC signals for hard loci (motivated by the GOLGA8A cohort analysis, where neither phasing method alone captures every case).Calling robustness
--supportso a copy-number gain (which inflates total read count, and thus the fraction-based threshold) cannot swallow a genuine minority allele.Phasing strategies (
--phasing)ward(default): length-weighted Levenshtein + Ward clustering (groups primarily by length).dbscan: DBSCAN on length-invariant k-mer composition features — keeps a length-variable expansion together as one allele where Ward fragments it. ReportsNCLUSTERS.both(QC): reports the Ward call but also runs DBSCAN and raisesDISCORDANT_LENGTH(+DBSCAN_RB) when the two differ by >2x on the longer allele. Deliberately over-flags for triage; genotype unchanged from Ward.New QC outputs (on by default for
--unphased)EXPANSION_OUTLIER: ≥2 reads more than 2x the longer called allele — catches large expansions missed by both alleles (dropped to noise/outliers or removed as length outliers during consensus).IMPRECISE_LENGTH: called-allele read-length CV > 0.2 — flags continuous/long-tailed loci where a single consensus length is not representative.MRLFORMAT field: median read length per allele (relative to ref), robust to a long tail; compare againstRB.Parameters
The QC thresholds and DBSCAN
eps/length-weightare kept as documented constants rather than CLI args (a parameter is only worth exposing if a user can tell how to tune it). The README "Tuning and hardcoded parameters" section documents each value, its source location, and tuning intuition, and offers to collaborate on challenging loci (referencing the bioRxiv benchmark).Breaking changes
--phasing-strategyrenamed to--phasing.--dbscan-epsand--dbscan-length-weightremoved (now hardcoded constants).Testing
cargo test: 67 unit + 3 integration pass; new tests coverEXPANSION_OUTLIER, length discordance, and the median/imprecision logic.cargo fmt --checkandcargo clippy --all-targetsclean.--phasing bothand the new FORMAT/INFO output on the test BAM.🤖 Generated with Claude Code