CNV-robust phasing, DBSCAN strategy, and expansion-sensitivity QC flags

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "STRdust"
-version = "0.18.0"
+version = "0.19.0"
 edition = "2024"
 
 

README.md

@@ -51,9 +51,11 @@ OPTIONS:
     -t, --threads <THREADS>            Number of parallel threads to use [default: 1]
         --sample <SAMPLE>              Sample name to use in VCF header, if not provided, the bam file name is used
         --somatic                      Print information on somatic variability
-        --unphased                     Reads are not phased, will use hierarchical clustering to phase expansions
+        --unphased                     Reads are not phased, will cluster the reads to phase expansions
         --consensus-reads              Maximum number of reads to use to build the consensus sequence [default: 20]
         --find-outliers                Identify poorly supported outlier expansions (only with --unphased)
+        --min-haplotype-fraction <F>   Minimum fraction of reads for a cluster to be a haplotype (only with --unphased) [default: 0.1]
+        --phasing <STRATEGY>           How to split unphased reads into haplotypes: 'ward', 'dbscan' or 'both' (only with --unphased) [default: ward]
         --haploid <HAPLOID>            comma-separated list of haploid (sex) chromosomes
         --alignment-all                Always use full alignment (disable fast reference check via CIGAR)
     -h, --help                         Print help information
@@ -65,14 +67,24 @@ OPTIONS:
 - BED files can be provided in plain text or gzipped format (`.bed` or `.bed.gz`)
 - Lowering the number of consensus reads may lead to lesser accurate alternative allele sequences (selecting randomly from the reads), but may greatly improve speed. Note that in the case of somatic length variation, a small number of randomly selected reads may lead to a bias and not be representative of the true repeat length.
 - Genotyping known pathogenic repeats with the `--pathogenic` flag will return a VCF with the pathogenic STRs from STRchive, but currently only for the GRCh38 reference.
+- For unphased data (`--unphased`), STRdust splits the reads into (at most) two haplotypes before building consensus. Three strategies are available via `--phasing`:
+  - `ward` (default): hierarchical (Ward) clustering on a **length-weighted Levenshtein distance**. The distance is dominated by how much two reads differ in length, so reads are grouped primarily by length. Robust for the common case where alleles differ mainly in length, but it tends to fragment a single length-variable expansion across several length bins.
+  - `dbscan` (experimental): DBSCAN on **length-invariant k-mer composition feature vectors**, i.e. it groups reads by their *sequence composition* (which motif they are made of) rather than primarily by length. This is the key difference from `ward`: two reads of the same motif cluster together even when their lengths differ a lot, so a length-variable expansion is kept together as one allele. The trade-off is that the reference and expanded alleles must differ in composition for DBSCAN to separate them. The two largest clusters become the haplotypes; remaining clusters and noise reads are reported as `OUTLIERS`, and the total number of clusters is reported in `NCLUSTERS` (so loci with `NCLUSTERS > 2`, i.e. complex/multi-population loci, can be flagged downstream). Its internal parameters (neighbourhood radius and length weight) are hardcoded — see [Tuning and hardcoded parameters](#tuning-and-hardcoded-parameters).
+  - `both` (QC mode): report the `ward` call as usual, but additionally run `dbscan` and, when the two disagree by more than 2x on the longer allele, raise a `DISCORDANT_LENGTH` flag and report the DBSCAN allele lengths in `DBSCAN_RB`. This deliberately over-flags (it prioritises sensitivity) and is intended as a triage signal: a worklist of loci worth reviewing where the two orthogonal clustering approaches disagree. The reported genotype is unchanged from `ward`.
 - By default, STRdust uses a fast reference check (QUICKREF) to skip full alignment at loci that appear to be homozygous reference. It inspects the CIGAR strings of the first 25 reads spanning a locus, and if at least 5 are found and none show a length difference from the reference of more than 3 bp, the locus is called 0|0 immediately. Loci called this way are marked with a `QUICKREF` flag in the VCF INFO field. This substantially speeds up runs on samples with many reference-like loci. To disable this optimisation and always perform full alignment, use `--alignment-all`.
 
 ## Output format
 
 STRdust produces a VCF file per sample. The consensus sequence is in the ALT field, with sequences from each read in the SEQS INFO field (when running with --somatic).
 The FRB FORMAT field is the total repeat length, of the two alleles, in nucleotides. The RB field is the difference between the indidiual allele lengths and the reference length.
+The MRL FORMAT field is the median read length per allele (relative to the reference). Because it is a median rather than the POA consensus length, it is more robust to a long tail of length-variable reads, and comparing it to RB highlights alleles whose consensus length is being pulled by a skewed read distribution.
 The SC FORMAT field is a measure of accuracy of the consensus sequence compared to the overlap graph from the individual reads, which could be influenced by the presence of sequencing errors or somatic variation.
 
+Two locus-level quality flags are raised automatically when running with `--unphased` (no extra options needed):
+
+- `EXPANSION_OUTLIER`: at least 2 reads are more than 2x longer than the longer called allele. This catches a larger expansion that was missed by both alleles because its few reads fell to noise/outliers during clustering, or were dropped as length outliers when building the consensus. The thresholds (2 reads, 2x) are currently fixed; they were chosen to suppress single-read artifacts while still surfacing genuine multi-read expansions.
+- `IMPRECISE_LENGTH`: a called allele's reads have a wide length spread (coefficient of variation, std_dev/mean, above 0.2). At such loci the underlying length distribution is continuous or long-tailed, so a single consensus length is not a faithful summary — compare RB and MRL, and treat the call with caution. The 0.2 threshold is currently fixed.
+
 Example output:
 
 ```text
@@ -81,18 +93,38 @@ Example output:
 ##INFO=<ID=STDEV,Number=2,Type=Integer,Description="Standard deviation of the repeat length">
 ##INFO=<ID=SEQS,Number=1,Type=String,Description="Sequences supporting the two alleles">
 ##INFO=<ID=OUTLIERS,Number=1,Type=String,Description="Outlier sequences much longer than the alleles">
+##INFO=<ID=NCLUSTERS,Number=1,Type=Integer,Description="Number of read clusters found by the DBSCAN phasing strategy (>2 indicates a complex multi-population locus)">
 ##INFO=<ID=CLUSTERFAILURE,Number=0,Type=Flag,Description="If unphased input failed to cluster in two haplotype">
+##INFO=<ID=EXPANSION_OUTLIER,Number=0,Type=Flag,Description="At least 2 reads are more than 2x longer than the longer called allele, suggesting a larger expansion missed by both alleles (only with --unphased)">
+##INFO=<ID=IMPRECISE_LENGTH,Number=0,Type=Flag,Description="A called allele has a wide read-length spread (coefficient of variation > 0.2); the reported consensus length may not be representative">
+##INFO=<ID=DISCORDANT_LENGTH,Number=0,Type=Flag,Description="With --phasing both: the reported Ward call and the DBSCAN call differ by more than 2x on the longer allele; see DBSCAN_RB (QC only, fires liberally)">
+##INFO=<ID=DBSCAN_RB,Number=2,Type=Integer,Description="With --phasing both: DBSCAN allele lengths relative to reference, reported when DISCORDANT_LENGTH is set">
 ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
 ##FORMAT=<ID=RB,Number=2,Type=Integer,Description="Repeat length of the two alleles in bases relative to reference">
 ##FORMAT=<ID=FRB,Number=2,Type=Integer,Description="Full repeat length of the two alleles in bases">
+##FORMAT=<ID=MRL,Number=2,Type=Integer,Description="Median read length of the two alleles' clusters in bases relative to reference, robust to a long length tail">
 ##FORMAT=<ID=PS,Number=1,Type=Integer,Description="Phase set identifier">
 ##FORMAT=<ID=SUP,Number=2,Type=Integer,Description="Read support per allele">
 ##FORMAT=<ID=SC,Number=2,Type=Integer,Description="Consensus score per allele">
 #CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  HG00271.hg38
-chr1    1435798 .       TGGCGCGGAGCGGCGCGGAGCG  GCTGGCGCGGAGCGGCGCGGA,GCGGGCGCGCGCAGGA  .       .       END=1435818;STDEV=1,2     GT:RB:FRB:SUP:SC        1|2:1,-4:21,16:18,6:63,41
+chr1    1435798 .       TGGCGCGGAGCGGCGCGGAGCG  GCTGGCGCGGAGCGGCGCGGA,GCGGGCGCGCGCAGGA  .       .       END=1435818;STDEV=1,2     GT:RB:FRB:MRL:SUP:SC        1|2:1,-4:21,16:1,-4:18,6:63,41
 chr1    57367044        .       AAAATAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATAAAT     AAAATAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATAAA,AAAATAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATAAA        .       .       END=57367125;STDEV=3,0 GT:RB:FRB:SUP:SC 1|2:-9,-34:72,47:17,12:216,141
 ```
 
+## Tuning and hardcoded parameters
+
+STRdust was developed while investigating the [pathogenic GOLGA8A repeat expansion](https://www.nature.com/articles/s41588-026-02537-7), and a [benchmark of repeat genotypers](https://www.biorxiv.org/content/10.64898/2026.02.28.708646v1) found it to be the most sensitive of the tested tools for detecting actual pathogenic expansions. Maximising that sensitivity is a priority, and several of the heuristics that support it are intentionally kept as **hardcoded constants** rather than command-line arguments. A parameter is only worth exposing if a user can tell how to tune it, and exposing all of these would inflate the option list without helping most users. Each is defined as a documented constant in the source so it can be changed (and promoted to a CLI argument) if real cohort experience shows a need:
+
+| Parameter | Value | Where (source) | Tuning intuition |
+|---|---|---|---|
+| `EXPANSION_OUTLIER` min reads / ratio | ≥ 2 reads, > 2x longer allele | `src/genotype.rs` | lower either for more sensitivity to few-read / single-read expansions, at the cost of more artifact flags |
+| `IMPRECISE_LENGTH` length CV | > 0.2 | `src/consensus.rs` | lower flags more loci as length-imprecise; raise to flag only the most spread-out |
+| `DISCORDANT_LENGTH` ratio (`--phasing both`) | > 2x longer allele | `src/genotype.rs` | lower to surface smaller Ward/DBSCAN disagreements for QC review |
+| DBSCAN neighbourhood radius (`eps`) | 0.4 | `src/phase_insertions.rs` | smaller = tighter/more clusters (more splitting, more reads dropped to noise → can undercall length-variable expansions); larger = looser clusters that can merge distinct alleles. Adjust in ±0.1 steps and watch `NCLUSTERS` |
+| DBSCAN length weight | 0.3 | `src/phase_insertions.rs` | lower = composition dominates (keeps length-variable expansions together); higher = length matters more (separates same-motif alleles that differ only in length, behaving more like `ward`) |
+
+If you have a specifically challenging repeat expansion where the defaults do not work well, we would be happy to work with you — please [open an issue](https://github.com/wdecoster/STRdust/issues) with details.
+
 ## Development
 
 ### Getting Started for Contributors

src/consensus.rs

@@ -1,19 +1,49 @@
 use bio::alignment::{pairwise::Scoring, poa::Aligner};
 use log::debug;
+use rand::SeedableRng;
+use rand::rngs::StdRng;
 use rand::seq::IteratorRandom;
 use std::fmt;
 
+/// Fixed seed for the read downsampling RNG. Downsampling is purely a
+/// performance/memory measure, so the choice of subset should not make the
+/// genotype non-reproducible: seeding with a constant makes every run (and
+/// every comparison between runs) deterministic and independent of locus or
+/// thread ordering.
+const DOWNSAMPLE_SEED: u64 = 42;
+
+/// A called allele is flagged IMPRECISE_LENGTH when the coefficient of variation
+/// (std_dev / mean) of its cluster's read lengths exceeds this threshold. 0.2 (20%
+/// spread) separates clean, near-fixed-length alleles from continuous/long-tailed
+/// length distributions (e.g. the GOLGA8A repeat) where a single consensus length
+/// is not a faithful representation of the underlying reads. Hardcoded for now;
+/// can be promoted to a CLI argument later if cohort experience warrants it.
+const IMPRECISE_LENGTH_CV: f64 = 0.2;
+
 #[derive(Clone)]
 pub struct Consensus {
     pub seq: Option<String>,
     pub support: usize,
     pub std_dev: usize,
     pub score: i32,
+    /// Median full length of the cluster reads (before length-outlier removal),
+    /// reported as MRL. More robust to a long tail than the POA consensus length.
+    pub median_length: usize,
+    /// True when the cluster's read-length coefficient of variation exceeds
+    /// [`IMPRECISE_LENGTH_CV`]; surfaced as the IMPRECISE_LENGTH locus flag.
+    pub imprecise: bool,
 }
 
 impl Default for Consensus {
     fn default() -> Consensus {
-        Consensus { seq: None, support: 0, std_dev: 0, score: -1 }
+        Consensus {
+            seq: None,
+            support: 0,
+            std_dev: 0,
+            score: -1,
+            median_length: 0,
+            imprecise: false,
+        }
     }
 }
 
@@ -37,25 +67,39 @@ pub fn consensus(
     repeat: &crate::repeats::RepeatInterval,
 ) -> Consensus {
     if seqs.is_empty() {
-        return Consensus { seq: None, support: 0, std_dev: 0, score: -1 };
+        return Consensus::default();
     }
     let num_reads_ = seqs.len();
-    let (seqs, std_dev) = remove_outliers(seqs, repeat);
+    // Length statistics over all cluster reads (before length-outlier removal), so the
+    // reported median and the imprecision verdict reflect the true read distribution.
+    let (mean, std_dev, median_length) = length_stats(seqs);
+    let imprecise = mean > 0 && (std_dev as f64 / mean as f64) > IMPRECISE_LENGTH_CV;
+    let seqs = remove_outliers(seqs, mean, std_dev, repeat);
     let num_reads = seqs.len();
     debug!("{repeat}: Kept {}/{} reads after dropping outliers", num_reads, num_reads_);
     if num_reads < support {
-        Consensus { seq: None, support: num_reads, std_dev, score: -1 }
+        Consensus { seq: None, support: num_reads, std_dev, score: -1, median_length, imprecise }
     } else if consensus_reads == 1 {
         // if only on read should be used to generate the consensus, the consensus is a randomly selected read
-        let seq = seqs.into_iter().choose(&mut rand::rng()).unwrap();
-        Consensus { seq: Some(seq.to_string()), support: num_reads, std_dev, score: 0 }
+        let seq = seqs
+            .into_iter()
+            .choose(&mut StdRng::seed_from_u64(DOWNSAMPLE_SEED))
+            .unwrap();
+        Consensus {
+            seq: Some(seq.to_string()),
+            support: num_reads,
+            std_dev,
+            score: 0,
+            median_length,
+            imprecise,
+        }
     } else {
         // if there are more than <consensus_reads> reads, downsample before taking the consensus
         // for performance and memory reasons
         let seqs_bytes = if num_reads > consensus_reads {
             debug!("{repeat}: Too many reads, downsampling to {consensus_reads}");
             seqs.into_iter()
-                .sample(&mut rand::rng(), consensus_reads)
+                .sample(&mut StdRng::seed_from_u64(DOWNSAMPLE_SEED), consensus_reads)
                 .into_iter()
                 .map(|seq| seq.bytes().collect::<Vec<u8>>())
                 .collect::<Vec<Vec<u8>>>()
@@ -84,52 +128,90 @@ pub fn consensus(
             support: num_reads,
             std_dev,
             score,
+            median_length,
+            imprecise,
         }
     }
 }
 
-fn remove_outliers<'a>(
-    seqs: &'a [String],
-    repeat: &crate::repeats::RepeatInterval,
-) -> (Vec<&'a String>, usize) {
-    // remove sequences that are shorter or longer than two standard deviations from the mean
-    // except if the stdev is small
-    let lengths = seqs.iter().map(|x| x.len()).collect::<Vec<usize>>();
-    debug!("{repeat}: lengths: {:?}", lengths);
-
+/// Mean, (floored) standard deviation, and median of the sequence lengths.
+fn length_stats(seqs: &[String]) -> (usize, usize, usize) {
+    let mut lengths = seqs.iter().map(|x| x.len()).collect::<Vec<usize>>();
     let mean = lengths.iter().sum::<usize>() / lengths.len();
     let variance = lengths
         .iter()
         .map(|x| (*x as isize - mean as isize).pow(2) as usize)
         .sum::<usize>()
         / lengths.len();
-    // note that the casting to usize will floor the std_dev to the integer below, rather than properly rounding it to the nearest integer
-    // so this keeps the std_dev smaller than it really is, but not by a lot
-    // the places where this matter are probably negligible
+    // casting to usize floors the std_dev to the integer below, rather than rounding to nearest,
+    // so this keeps the std_dev slightly smaller than it really is, but not by a lot
     let std_dev = (variance as f64).sqrt() as usize;
-    debug!("mean: {}, std_dev: {}", mean, std_dev);
+    lengths.sort_unstable();
+    let median = if lengths.len() % 2 == 0 {
+        (lengths[lengths.len() / 2] + lengths[lengths.len() / 2 - 1]) / 2
+    } else {
+        lengths[lengths.len() / 2]
+    };
+    (mean, std_dev, median)
+}
+
+fn remove_outliers<'a>(
+    seqs: &'a [String],
+    mean: usize,
+    std_dev: usize,
+    repeat: &crate::repeats::RepeatInterval,
+) -> Vec<&'a String> {
+    // remove sequences that are shorter or longer than two standard deviations from the mean
+    // except if the stdev is small
+    debug!("{repeat}: mean: {}, std_dev: {}", mean, std_dev);
     if std_dev < 5 {
         debug!("std_dev < 5, not removing any outliers");
-        (seqs.iter().collect::<Vec<&String>>(), std_dev)
+        seqs.iter().collect::<Vec<&String>>()
     } else {
         // avoid underflowing usize
         let min_val = mean.saturating_sub(2 * std_dev);
         let max_val = mean + 2 * std_dev;
         debug!("Removing outliers outside of [{},{}]", min_val, max_val);
-        let filtered_seqs = seqs
-            .iter()
-            .zip(lengths.iter())
-            .filter(|(_, len)| **len > min_val && **len < max_val)
-            .map(|(seq, _)| seq)
-            .collect::<Vec<&String>>();
-        (filtered_seqs, std_dev)
+        seqs.iter()
+            .filter(|seq| seq.len() > min_val && seq.len() < max_val)
+            .collect::<Vec<&String>>()
     }
 }
 
 #[cfg(test)]
 mod tests {
     use super::*;
     use bio::alignment::pairwise::Scoring;
+
+    fn dummy_repeat() -> crate::repeats::RepeatInterval {
+        crate::repeats::RepeatInterval {
+            chrom: "chr1".to_string(),
+            start: 1,
+            end: 100,
+            created: None,
+        }
+    }
+
+    #[test]
+    fn test_consensus_reports_median_and_not_imprecise_when_tight() {
+        // near-uniform lengths (~30 bp): low CV -> not imprecise, median ~30
+        let seqs: Vec<String> = (0..10).map(|i| "A".repeat(30 + i % 2)).collect();
+        let cons = consensus(&seqs, 2, 1, &dummy_repeat());
+        assert!(!cons.imprecise, "tight length distribution should not be imprecise");
+        assert!((29..=31).contains(&cons.median_length), "median {}", cons.median_length);
+    }
+
+    #[test]
+    fn test_consensus_flags_imprecise_when_length_spread_is_wide() {
+        // continuous/long-tailed lengths: CV well above 0.2 -> imprecise
+        let seqs: Vec<String> = vec![30, 50, 120, 300, 900, 2700]
+            .into_iter()
+            .map(|l| "A".repeat(l))
+            .collect();
+        let cons = consensus(&seqs, 2, 1, &dummy_repeat());
+        assert!(cons.imprecise, "wide length distribution should be flagged imprecise");
+    }
+
     #[test]
     fn test_consensus() {
         // I created this test because these sequences segfaulted on bianca