From 9a3970037976b06a51b5f933a8ac0798d76bd1b8 Mon Sep 17 00:00:00 2001
From: wdecoster <decosterwouter@gmail.com>
Date: Thu, 4 Jun 2026 10:02:23 +0200
Subject: [PATCH 1/8] fix: cap haplotype cluster-size threshold by --support
 (CNV-robust)

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>
---
 src/genotype.rs         |  2 +
 src/phase_insertions.rs | 90 +++++++++++++++++++++++++++++++++++++++--
 2 files changed, 89 insertions(+), 3 deletions(-)

diff --git a/src/genotype.rs b/src/genotype.rs
index 400853c..fe78d39 100644
--- a/src/genotype.rs
+++ b/src/genotype.rs
@@ -276,6 +276,7 @@ pub fn genotype_with_extracted_reads(
             repeat,
             args.find_outliers,
             args.min_haplotype_fraction,
+            args.support,
         );
         match phased.hap2 {
             Some(phase2) => {
@@ -556,6 +557,7 @@ fn genotype_repeat(
             repeat,
             args.find_outliers,
             args.min_haplotype_fraction,
+            args.support,
         );
         match phased.hap2 {
             Some(phase2) => {
diff --git a/src/phase_insertions.rs b/src/phase_insertions.rs
index 70f2f56..665b2e4 100644
--- a/src/phase_insertions.rs
+++ b/src/phase_insertions.rs
@@ -15,6 +15,7 @@ pub fn split(
     repeat: &crate::repeats::RepeatInterval,
     check_outliers: bool,
     min_haplotype_fraction: f32,
+    support: usize,
 ) -> SplitSequences {
     // the insertions are from an unphased experiment
     // and should be split in one (if homozygous) or two haplotypes
@@ -112,12 +113,19 @@ pub fn split(
     // create a vector with candidate haplotype-clusters
     let mut haplotype_clusters = vec![];
     // clusters have to represent at least min_haplotype_fraction of the reads, but never less than 1
-    let min_cluster_size = max((insertions.len() as f32 * min_haplotype_fraction) as usize, 1);
+    // the fraction-based threshold is capped by `support` so that copy-number gains (which inflate
+    // the total read count, and thus the fraction-based threshold) cannot swallow a genuine minority
+    // allele that is supported by at least `support` reads
+    let min_cluster_size = max(
+        ((insertions.len() as f32 * min_haplotype_fraction) as usize).min(support),
+        1,
+    );
     debug!(
-        "{repeat}: Minimum cluster size: {} reads ({}% of {} total)",
+        "{repeat}: Minimum cluster size: {} reads (min of {}% of {} total and support {})",
         min_cluster_size,
         min_haplotype_fraction * 100.0,
-        insertions.len()
+        insertions.len(),
+        support
     );
 
     for (index, step) in dend.steps().iter().enumerate() {
@@ -383,6 +391,76 @@ mod tests {
     #[allow(unused_imports)]
     use super::*;
 
+    #[test]
+    fn test_split_expansion_minority_allele() {
+        // Regression test for chr15:34419425-34419451 from rr_NA99_170.log
+        // A real, length-variable expansion is present on only 5/44 reads (the rest
+        // are the short reference allele). The expansion reads must be recovered as a
+        // haplotype and must NOT all be discarded as length outliers.
+        let r43 = "TCTTTCTTTCTTTCCTTTCCTTTCCTTTCCTTTCCTTTCCTTCCTTCCCTTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTTTTTCTTTCTTTCTTT".to_string();
+        let r2 = "TCTTTCTTTCTTCCTTTCCTTTCCTTTCCTTTCCTTTCCTTCCTTCCCTTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTTTCTTTCTTTCTTT".to_string();
+        let r33 = "CTTTCTTTCCTTTCCTTTCCTTTCCTTTCCTTTCCTTCCTTCCCCCTGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAAAAAAAAAAAAAAAAAAAAAAAATCTCTCTCTCTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTTTCTTTCTTTCTTTCTTTCTTTC".to_string();
+        let r8 = "TCTTTCTTTCTTTCCTTTCCTTTCCTTTCCTTTCCTTTCCTTCCTTCCCTTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTTTCTTCTTTCTTT".to_string();
+        let r23 = "TCTTTCTTTCCTTTCCTTTCCTTTCCTTTCCTTTCCTTCCTTCCCTTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTTCTCTCTCTCTCTCTCTCTCTTTCTTTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTTTCTTTCTTTCTTT".to_string();
+
+        let mut insertions = vec![r43, r2, r33, r8, r23];
+        // 39 short reference-allele reads (lengths matching the logged distribution)
+        for _ in 0..20 {
+            insertions.push("TTTCTTTCTTTCTTTCTTTCTTTCTTT".to_string()); // 27
+        }
+        for _ in 0..7 {
+            insertions.push("TTTCTTTCTTTCTTTCTTTCTTTCTTTT".to_string()); // 28
+        }
+        for _ in 0..2 {
+            insertions.push("TTTCTTTCTTTCTTTCTTTCTTTCTT".to_string()); // 26
+        }
+        for _ in 0..2 {
+            insertions.push("TTTCTTTCTTTCTTTCTTTCTTTCTTTTT".to_string()); // 29
+        }
+        insertions.push("CTTTTCTTTCTTTCTTTCTTTC".to_string()); // 22
+        insertions.push("TTTCTTTCTTTCTTTCTTTCTTT".to_string()); // 23
+        insertions.push("TTTTCTTTCTTTCTTTCTTTCTTT".to_string()); // 24
+        insertions.push("TTTCTTTCTTTCTTTCTTTCTTTCTTTCTTT".to_string()); // 31
+        insertions.push("AGTTTCTTTCTTTCTTTCTTTCTTTTTTTCTTT".to_string()); // 33
+        insertions.push("TTTCTTTCTTTCTTTCTTTTTTTCTTTTTTTCTTT".to_string()); // 35
+        for _ in 0..2 {
+            insertions.push("TTTCTTTCTTTCTTTCTTTCTTTCTTTTTTTTTT".to_string()); // 34
+        }
+        assert_eq!(insertions.len(), 44);
+
+        use rand::seq::SliceRandom;
+        let mut rng = rand::rng();
+        insertions.shuffle(&mut rng);
+
+        let splitseqs = split(
+            &insertions,
+            &crate::repeats::RepeatInterval {
+                chrom: "chr15".to_string(),
+                start: 34419425,
+                end: 34419451,
+                created: None,
+            },
+            true,
+            0.1,
+            2,
+        );
+
+        let hap2 = splitseqs
+            .hap2
+            .expect("expected a heterozygous call with two haplotypes");
+        // One of the two haplotypes should capture the expansion (median length >> reference)
+        let expansion_captured =
+            find_median(&splitseqs.hap1) > 100 || find_median(&hap2) > 100;
+        let n_outliers = splitseqs.outliers.as_ref().map_or(0, |o| o.len());
+        assert!(
+            expansion_captured,
+            "expansion not captured as a haplotype; hap1 median {}, hap2 median {}, {} outliers",
+            find_median(&splitseqs.hap1),
+            find_median(&hap2),
+            n_outliers
+        );
+    }
+
     #[test]
     fn test_split_length() {
         // test that the split function identifies two haplotype, mainly based on insertion length
@@ -409,6 +487,7 @@ mod tests {
             },
             false,
             0.1,
+            3,
         );
         assert!(splitseqs.hap1.len() == splitseqs.hap2.unwrap().len());
         // check that all sequences in hap1 are the same length
@@ -453,6 +532,7 @@ mod tests {
             },
             false,
             0.1,
+            3,
         );
         let mut hap1 = splitseqs.hap1;
         let mut hap2 = splitseqs.hap2.unwrap();
@@ -494,6 +574,7 @@ mod tests {
             },
             false,
             0.1,
+            3,
         );
         assert!(splitseqs.hap1.len() + splitseqs.hap2.unwrap().len() == insertions.len());
     }
@@ -538,6 +619,7 @@ mod tests {
             },
             false,
             0.1,
+            3,
         );
         assert!(splitseqs.hap2.is_none());
         println!("hap1: {:?}", splitseqs.hap1);
@@ -580,6 +662,7 @@ mod tests {
             },
             false,
             0.1,
+            3,
         );
         let mut hap1 = splitseqs.hap1;
         let mut hap2 = splitseqs.hap2.unwrap();
@@ -653,6 +736,7 @@ mod tests {
             },
             false,
             0.1,
+            3,
         );
         let mut hap1 = splitseqs.hap1;
         let mut hap2 = splitseqs.hap2.unwrap();

From 0d4e3eb4625d73c9af52eeef921757d3403aff40 Mon Sep 17 00:00:00 2001
From: wdecoster <decosterwouter@gmail.com>
Date: Thu, 4 Jun 2026 10:41:30 +0200
Subject: [PATCH 2/8] fix: seed downsampling RNG for deterministic consensus
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

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>
---
 src/consensus.rs | 13 +++++++++++--
 1 file changed, 11 insertions(+), 2 deletions(-)

diff --git a/src/consensus.rs b/src/consensus.rs
index d4e7672..ea71e13 100644
--- a/src/consensus.rs
+++ b/src/consensus.rs
@@ -1,8 +1,17 @@
 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;
+
 #[derive(Clone)]
 pub struct Consensus {
     pub seq: Option<String>,
@@ -47,7 +56,7 @@ pub fn consensus(
         Consensus { seq: None, support: num_reads, std_dev, score: -1 }
     } 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();
+        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 }
     } else {
         // if there are more than <consensus_reads> reads, downsample before taking the consensus
@@ -55,7 +64,7 @@ pub fn consensus(
         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>>>()

From 55761f6906dc0c3ab7c1c6b33d89a12987dfe19f Mon Sep 17 00:00:00 2001
From: wdecoster <decosterwouter@gmail.com>
Date: Thu, 4 Jun 2026 11:16:43 +0200
Subject: [PATCH 3/8] feat: optional DBSCAN phasing strategy
 (--phasing-strategy dbscan)

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>
---
 src/dbscan.rs           | 168 ++++++++++++++++++++++++++++++++++++
 src/features.rs         | 186 ++++++++++++++++++++++++++++++++++++++++
 src/genotype.rs         |  71 ++++++++++++---
 src/main.rs             |  30 ++++++-
 src/phase_insertions.rs | 167 +++++++++++++++++++++++++++++++++---
 src/vcf.rs              |  23 ++++-
 6 files changed, 617 insertions(+), 28 deletions(-)
 create mode 100644 src/dbscan.rs
 create mode 100644 src/features.rs

diff --git a/src/dbscan.rs b/src/dbscan.rs
new file mode 100644
index 0000000..d92b5e5
--- /dev/null
+++ b/src/dbscan.rs
@@ -0,0 +1,168 @@
+//! DBSCAN clustering for the optional `--phasing-strategy dbscan` path.
+//!
+//! Adapted from the `trout` cohort-outlier tool, which used a custom DBSCAN
+//! over a precomputed CSR adjacency for speed (one pairwise distance per pair,
+//! no per-range-query allocations). Here it returns full cluster *labels*
+//! (not just a noise mask) because phasing needs the read groupings.
+//!
+//! Distances are squared Euclidean over the feature vectors from
+//! [`crate::features`]; `min_samples` is the standard DBSCAN core-point
+//! threshold (a point counts itself, so a core point has ≥ `min_samples`
+//! members including itself).
+
+/// Cluster label for each point. `None` == noise (DBSCAN outlier).
+pub type Label = Option<usize>;
+
+/// Run DBSCAN and return a per-point label vector. Clustered points get a
+/// `Some(cluster_id)`; noise points get `None`. Cluster ids are assigned in
+/// the order clusters are discovered (0, 1, 2, …).
+pub fn cluster(points: &[Vec<f64>], eps: f64, min_samples: usize) -> Vec<Label> {
+    let n = points.len();
+    if n == 0 {
+        return Vec::new();
+    }
+    let eps_sq = eps * eps;
+
+    // Collect neighbour edges once, then build CSR adjacency from a single
+    // allocation pair (avoids growing n separate small Vecs).
+    let mut edges: Vec<(u32, u32)> = Vec::new();
+    for i in 0..n {
+        let pi = points[i].as_slice();
+        for (j, pj_vec) in points.iter().enumerate().skip(i + 1) {
+            let pj = pj_vec.as_slice();
+            // LLVM auto-vectorises this zipped sum on f64 slices.
+            let dsq: f64 = pi
+                .iter()
+                .zip(pj.iter())
+                .map(|(a, b)| {
+                    let d = a - b;
+                    d * d
+                })
+                .sum();
+            if dsq <= eps_sq {
+                edges.push((i as u32, j as u32));
+            }
+        }
+    }
+
+    // CSR build: row_starts[i..i+1] indexes col_indices for point i's neighbours.
+    let mut row_lens = vec![1usize; n]; // each point includes itself
+    for &(i, j) in &edges {
+        row_lens[i as usize] += 1;
+        row_lens[j as usize] += 1;
+    }
+    let mut row_starts = vec![0usize; n + 1];
+    for i in 0..n {
+        row_starts[i + 1] = row_starts[i] + row_lens[i];
+    }
+    let mut col_indices: Vec<u32> = vec![0u32; row_starts[n]];
+    let mut cursor = row_starts.clone();
+    for i in 0..n {
+        col_indices[cursor[i]] = i as u32; // self-loop
+        cursor[i] += 1;
+    }
+    for &(i, j) in &edges {
+        col_indices[cursor[i as usize]] = j;
+        cursor[i as usize] += 1;
+        col_indices[cursor[j as usize]] = i;
+        cursor[j as usize] += 1;
+    }
+    drop(edges);
+    drop(cursor);
+    drop(row_lens);
+
+    let neighbours_of = |i: usize| -> &[u32] { &col_indices[row_starts[i]..row_starts[i + 1]] };
+
+    #[derive(Clone, Copy, PartialEq, Eq)]
+    enum St {
+        Unvisited,
+        Noise,
+        Clustered(usize),
+    }
+    let mut status = vec![St::Unvisited; n];
+    let mut queue: Vec<usize> = Vec::new();
+    let mut next_cluster = 0usize;
+
+    for p in 0..n {
+        if status[p] != St::Unvisited {
+            continue;
+        }
+        if neighbours_of(p).len() < min_samples {
+            status[p] = St::Noise;
+            continue;
+        }
+        // p is a core point — start a new cluster and expand.
+        let cid = next_cluster;
+        next_cluster += 1;
+        status[p] = St::Clustered(cid);
+        queue.clear();
+        for &q in neighbours_of(p) {
+            let q = q as usize;
+            if q != p && matches!(status[q], St::Unvisited | St::Noise) {
+                status[q] = St::Clustered(cid);
+                queue.push(q);
+            }
+        }
+        let mut idx = 0;
+        while idx < queue.len() {
+            let q = queue[idx];
+            idx += 1;
+            // Only core points propagate the cluster.
+            if neighbours_of(q).len() >= min_samples {
+                for &r in neighbours_of(q) {
+                    let r = r as usize;
+                    if matches!(status[r], St::Unvisited | St::Noise) {
+                        status[r] = St::Clustered(cid);
+                        queue.push(r);
+                    }
+                }
+            }
+        }
+    }
+
+    status
+        .into_iter()
+        .map(|s| match s {
+            St::Clustered(cid) => Some(cid),
+            _ => None,
+        })
+        .collect()
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn empty() {
+        assert!(cluster(&[], 0.5, 3).is_empty());
+    }
+
+    #[test]
+    fn two_separated_clusters_plus_noise() {
+        // two dense 1-D blobs near 0.0 and 10.0, one far-away noise point
+        let mut pts: Vec<Vec<f64>> = Vec::new();
+        for i in 0..6 {
+            pts.push(vec![0.0 + i as f64 * 0.01]);
+        }
+        for i in 0..6 {
+            pts.push(vec![10.0 + i as f64 * 0.01]);
+        }
+        pts.push(vec![100.0]); // noise
+        let labels = cluster(&pts, 0.5, 3);
+        // first 6 share a label, next 6 share another, last is noise
+        let c0 = labels[0];
+        let c1 = labels[6];
+        assert!(c0.is_some() && c1.is_some() && c0 != c1);
+        assert!(labels[..6].iter().all(|l| *l == c0));
+        assert!(labels[6..12].iter().all(|l| *l == c1));
+        assert_eq!(labels[12], None);
+    }
+
+    #[test]
+    fn all_noise_when_sparse() {
+        let pts: Vec<Vec<f64>> = (0..10).map(|i| vec![(i * 100) as f64]).collect();
+        let labels = cluster(&pts, 1.0, 3);
+        assert!(labels.iter().all(|l| l.is_none()));
+    }
+}
diff --git a/src/features.rs b/src/features.rs
new file mode 100644
index 0000000..95f3d22
--- /dev/null
+++ b/src/features.rs
@@ -0,0 +1,186 @@
+//! Sequence feature vectors for the optional DBSCAN phasing strategy.
+//!
+//! Each insertion is mapped to a point `[length_weight * normalized_length,
+//! canonical_kmer_freq_0, ...]`. The k-mer frequencies are length-invariant:
+//! two reads of the same motif have near-identical k-mer vectors regardless of
+//! how long the repeat is, which is what lets a length-variable expansion
+//! cluster together (unlike the length-dominated Levenshtein distance used by
+//! the default Ward strategy). Lifted from the `trout` cohort-outlier tool.
+//!
+//! Canonical k-mers fold cyclic rotations together (CAG ≡ AGC ≡ GCA) so the
+//! representation is agnostic to where in the motif a read happens to start.
+
+use std::collections::BTreeMap;
+
+/// Bounds for `-k auto` period detection (k=1 is too crude; 4^6 already gives
+/// ~700 canonical features, beyond which DBSCAN distances become noise).
+pub const AUTO_K_MIN: usize = 2;
+pub const AUTO_K_MAX: usize = 6;
+/// Fallback when no clean period is detectable: trinucleotide composition has
+/// the most signal among the common pathogenic motif lengths.
+pub const AUTO_K_DEFAULT: usize = 3;
+/// Match-rate threshold for `detect_period` (~35% mismatch tolerance).
+const PERIOD_MATCH_THRESHOLD: f64 = 0.65;
+
+/// Detect the dominant tandem-repeat period of `seq` by counting positions
+/// where `seq[i] == seq[i+p]`. Returns the smallest p in `AUTO_K_MIN..=k_max`
+/// clearing `PERIOD_MATCH_THRESHOLD` (the fundamental period), or None.
+pub fn detect_period(seq: &[u8], k_max: usize) -> Option<usize> {
+    if seq.len() < 2 * k_max {
+        return None;
+    }
+    for p in AUTO_K_MIN..=k_max {
+        let total = seq.len() - p;
+        let matches = (0..total).filter(|&i| seq[i] == seq[i + p]).count();
+        if (matches as f64) / (total as f64) >= PERIOD_MATCH_THRESHOLD {
+            return Some(p);
+        }
+    }
+    None
+}
+
+/// Precomputed mapping from every raw k-mer index to its canonical (lex-first
+/// cyclic rotation) compact index.
+pub struct KmerTable {
+    /// raw_idx → compact canonical index (length: 4^k)
+    pub compact: Vec<usize>,
+    /// number of canonical classes
+    pub n: usize,
+}
+
+impl KmerTable {
+    pub fn new(k: usize) -> Self {
+        let num_raw = 4usize.pow(k as u32);
+
+        let decode = |idx: usize| -> Vec<usize> {
+            (0..k)
+                .map(|pos| (idx / 4usize.pow((k - 1 - pos) as u32)) % 4)
+                .collect()
+        };
+        let encode = |bases: &[usize]| -> usize { bases.iter().fold(0, |acc, &b| acc * 4 + b) };
+
+        // For each raw kmer, the raw index of its canonical (min) rotation.
+        let canonical_raw: Vec<usize> = (0..num_raw)
+            .map(|idx| {
+                let bases = decode(idx);
+                (0..k)
+                    .map(|r| {
+                        let rotated: Vec<usize> = (0..k).map(|i| bases[(i + r) % k]).collect();
+                        encode(&rotated)
+                    })
+                    .min()
+                    .unwrap_or(idx)
+            })
+            .collect();
+
+        // Assign compact indices to canonical raws in sorted order (deterministic).
+        let mut seen: BTreeMap<usize, usize> = BTreeMap::new();
+        let mut next = 0usize;
+        let compact: Vec<usize> = canonical_raw
+            .iter()
+            .map(|&cr| {
+                *seen.entry(cr).or_insert_with(|| {
+                    let i = next;
+                    next += 1;
+                    i
+                })
+            })
+            .collect();
+
+        KmerTable { compact, n: next }
+    }
+}
+
+/// Build a feature vector: `[length_weight * normalized_length, canonical_kmer_freqs...]`.
+pub fn build_feature_vector(
+    seq: &str,
+    table: &KmerTable,
+    k: usize,
+    max_len: usize,
+    length_weight: f64,
+) -> Vec<f64> {
+    let normalized_len = if max_len > 0 {
+        seq.len() as f64 / max_len as f64
+    } else {
+        0.0
+    };
+    let mut features = vec![normalized_len * length_weight];
+    features.extend(kmer_frequencies(seq, table, k));
+    features
+}
+
+fn kmer_frequencies(seq: &str, table: &KmerTable, k: usize) -> Vec<f64> {
+    let seq_bytes = seq.as_bytes();
+    if seq_bytes.len() < k {
+        return vec![0.0; table.n];
+    }
+    let mut counts = vec![0u32; table.n];
+    let total = seq_bytes.len() - k + 1;
+    for window in seq_bytes.windows(k) {
+        if let Some(raw) = kmer_raw_index(window) {
+            counts[table.compact[raw]] += 1;
+        }
+    }
+    let total_f = total as f64;
+    counts.iter().map(|&c| c as f64 / total_f).collect()
+}
+
+fn kmer_raw_index(kmer: &[u8]) -> Option<usize> {
+    let mut idx = 0usize;
+    for &b in kmer {
+        idx *= 4;
+        idx += match b {
+            b'A' | b'a' => 0,
+            b'C' | b'c' => 1,
+            b'G' | b'g' => 2,
+            b'T' | b't' => 3,
+            _ => return None,
+        };
+    }
+    Some(idx)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_kmer_frequencies_sum_to_one() {
+        let table = KmerTable::new(2);
+        let v = build_feature_vector("ACGTACGT", &table, 2, 8, 1.0);
+        let sum: f64 = v[1..].iter().sum();
+        assert!((sum - 1.0).abs() < 1e-9);
+    }
+
+    #[test]
+    fn test_canonical_rotations_merged() {
+        let table = KmerTable::new(3);
+        // CAG, AGC, GCA all canonicalize to the same class
+        let cag = kmer_raw_index(b"CAG").unwrap();
+        let agc = kmer_raw_index(b"AGC").unwrap();
+        let gca = kmer_raw_index(b"GCA").unwrap();
+        assert_eq!(table.compact[cag], table.compact[agc]);
+        assert_eq!(table.compact[cag], table.compact[gca]);
+    }
+
+    #[test]
+    fn test_detect_period_cag() {
+        assert_eq!(detect_period(b"CAGCAGCAGCAGCAGCAG", 6), Some(3));
+    }
+
+    #[test]
+    fn test_length_invariance_of_composition() {
+        // Same motif, very different lengths -> near-identical kmer part.
+        let table = KmerTable::new(3);
+        let short = "CAG".repeat(10);
+        let long = "CAG".repeat(200);
+        let vs = build_feature_vector(&short, &table, 3, long.len(), 1.0);
+        let vl = build_feature_vector(&long, &table, 3, long.len(), 1.0);
+        let comp_dist: f64 = vs[1..]
+            .iter()
+            .zip(vl[1..].iter())
+            .map(|(a, b)| (a - b).powi(2))
+            .sum();
+        assert!(comp_dist < 1e-3, "composition should be length-invariant");
+    }
+}
diff --git a/src/genotype.rs b/src/genotype.rs
index fe78d39..c93e794 100644
--- a/src/genotype.rs
+++ b/src/genotype.rs
@@ -209,6 +209,8 @@ pub fn genotype_with_extracted_reads(
     } else {
         None
     };
+    // number of clusters found by the DBSCAN phasing strategy (None for Ward / phased input)
+    let mut n_clusters: Option<usize> = None;
 
     // alignments can be extracted in an unphased manner, if the chromosome is --haploid or the --unphased is set
     let unphased = (args.haploid.is_some()
@@ -271,13 +273,23 @@ pub fn genotype_with_extracted_reads(
         }
 
         debug!("{repeat}: Phasing {} insertions", insertions.len());
-        let phased = crate::phase_insertions::split(
-            &insertions,
-            repeat,
-            args.find_outliers,
-            args.min_haplotype_fraction,
-            args.support,
-        );
+        let phased = match args.phasing_strategy {
+            crate::PhasingStrategy::Ward => crate::phase_insertions::split(
+                &insertions,
+                repeat,
+                args.find_outliers,
+                args.min_haplotype_fraction,
+                args.support,
+            ),
+            crate::PhasingStrategy::Dbscan => crate::phase_insertions::split_dbscan(
+                &insertions,
+                repeat,
+                args.find_outliers,
+                args.support,
+                args.dbscan_eps,
+                args.dbscan_length_weight,
+            ),
+        };
         match phased.hap2 {
             Some(phase2) => {
                 consenses.push(crate::consensus::consensus(
@@ -320,6 +332,7 @@ pub fn genotype_with_extracted_reads(
                 }
             }
         }
+        n_clusters = phased.n_clusters;
         outliers = phased.outliers;
     } else {
         // Phased
@@ -371,6 +384,7 @@ pub fn genotype_with_extracted_reads(
         repeat_ref_seq,
         all_insertions,
         outliers,
+        n_clusters,
         repeat,
         reads.ps,
         flags,
@@ -482,6 +496,8 @@ fn genotype_repeat(
     } else {
         None
     };
+    // number of clusters found by the DBSCAN phasing strategy (None for Ward / phased input)
+    let mut n_clusters: Option<usize> = None;
 
     // The rest of the function has three mutually exclusive options from here.
     // Either the reads are from a haploid chromosome, unphased or phased by a tool like WhatsHap/hiphase/...
@@ -552,13 +568,23 @@ fn genotype_repeat(
         }
 
         debug!("{repeat}: Phasing {} insertions", insertions.len(),);
-        let phased = crate::phase_insertions::split(
-            &insertions,
-            repeat,
-            args.find_outliers,
-            args.min_haplotype_fraction,
-            args.support,
-        );
+        let phased = match args.phasing_strategy {
+            crate::PhasingStrategy::Ward => crate::phase_insertions::split(
+                &insertions,
+                repeat,
+                args.find_outliers,
+                args.min_haplotype_fraction,
+                args.support,
+            ),
+            crate::PhasingStrategy::Dbscan => crate::phase_insertions::split_dbscan(
+                &insertions,
+                repeat,
+                args.find_outliers,
+                args.support,
+                args.dbscan_eps,
+                args.dbscan_length_weight,
+            ),
+        };
         match phased.hap2 {
             Some(phase2) => {
                 consenses.push(crate::consensus::consensus(
@@ -602,6 +628,7 @@ fn genotype_repeat(
                 }
             }
         }
+        n_clusters = phased.n_clusters;
         if let Some(ref mut outliers_vec) = outliers
             && let Some(outliers_found) = phased.outliers
         {
@@ -637,6 +664,7 @@ fn genotype_repeat(
         repeat_ref_seq,
         all_insertions,
         outliers,
+        n_clusters,
         repeat,
         reads.ps,
         flags,
@@ -806,6 +834,9 @@ mod tests {
             unphased: false,
             find_outliers: false,
             min_haplotype_fraction: 0.1,
+            phasing_strategy: crate::PhasingStrategy::Ward,
+            dbscan_eps: 0.2,
+            dbscan_length_weight: 1.0,
             threads: 1,
             sample: None,
             haploid: None,
@@ -841,6 +872,9 @@ mod tests {
             unphased: true,
             find_outliers: false,
             min_haplotype_fraction: 0.1,
+            phasing_strategy: crate::PhasingStrategy::Ward,
+            dbscan_eps: 0.2,
+            dbscan_length_weight: 1.0,
             threads: 1,
             sample: None,
             haploid: Some(String::from("chr7")),
@@ -870,6 +904,9 @@ mod tests {
             unphased: true,
             find_outliers: false,
             min_haplotype_fraction: 0.1,
+            phasing_strategy: crate::PhasingStrategy::Ward,
+            dbscan_eps: 0.2,
+            dbscan_length_weight: 1.0,
             threads: 1,
             sample: None,
             haploid: None,
@@ -905,6 +942,9 @@ mod tests {
             unphased: false,
             find_outliers: false,
             min_haplotype_fraction: 0.1,
+            phasing_strategy: crate::PhasingStrategy::Ward,
+            dbscan_eps: 0.2,
+            dbscan_length_weight: 1.0,
             threads: 1,
             sample: None,
             haploid: None,
@@ -942,6 +982,9 @@ mod tests {
             unphased: false,
             find_outliers: false,
             min_haplotype_fraction: 0.1,
+            phasing_strategy: crate::PhasingStrategy::Ward,
+            dbscan_eps: 0.2,
+            dbscan_length_weight: 1.0,
             threads: 1,
             sample: None,
             haploid: None,
diff --git a/src/main.rs b/src/main.rs
index 0b6a9c3..273b842 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -1,5 +1,5 @@
 #![allow(non_snake_case)]
-use clap::Parser;
+use clap::{Parser, ValueEnum};
 use log::{info, warn};
 use std::path::PathBuf;
 
@@ -7,6 +7,8 @@ pub mod bam_pool;
 pub mod batching;
 pub mod call;
 pub mod consensus;
+pub mod dbscan;
+pub mod features;
 pub mod genotype;
 pub mod motif;
 pub mod parse_bam;
@@ -15,6 +17,15 @@ pub mod repeats;
 pub mod utils;
 pub mod vcf;
 
+/// Strategy for splitting unphased reads into haplotypes.
+#[derive(ValueEnum, Clone, Copy, Debug, PartialEq, Eq)]
+pub enum PhasingStrategy {
+    /// Length-weighted Levenshtein distance + Ward hierarchical clustering.
+    Ward,
+    /// k-mer composition feature vectors + DBSCAN.
+    Dbscan,
+}
+
 // The arguments end up in the Cli struct
 #[derive(Parser, Debug)]
 #[command(author, version, about = "Tool to genotype STRs from long reads", long_about = None)]
@@ -71,6 +82,23 @@ pub struct Cli {
     #[arg(long, default_value_t = 0.1)]
     min_haplotype_fraction: f32,
 
+    /// Strategy for splitting unphased reads into haplotypes (only with --unphased).
+    /// 'ward': length-weighted Levenshtein + hierarchical clustering (default).
+    /// 'dbscan': k-mer composition features + DBSCAN (experimental, robust to length-variable expansions).
+    #[arg(long, value_enum, default_value_t = PhasingStrategy::Ward)]
+    phasing_strategy: PhasingStrategy,
+
+    /// DBSCAN neighbourhood radius in normalized [0,1] feature space (only with --phasing-strategy dbscan).
+    /// Smaller values split more readily. Tune together with --dbscan-length-weight.
+    #[arg(long, default_value_t = 0.2)]
+    dbscan_eps: f64,
+
+    /// Weight of the (normalized) length axis relative to k-mer composition for DBSCAN
+    /// (only with --phasing-strategy dbscan). Lower values let length-variable expansions of the
+    /// same motif cluster together; higher values better separate same-motif alleles that differ only in length.
+    #[arg(long, default_value_t = 1.0)]
+    dbscan_length_weight: f64,
+
     /// comma-separated list of haploid (sex) chromosomes
     #[arg(long)]
     haploid: Option<String>,
diff --git a/src/phase_insertions.rs b/src/phase_insertions.rs
index 665b2e4..7bf6535 100644
--- a/src/phase_insertions.rs
+++ b/src/phase_insertions.rs
@@ -8,6 +8,9 @@ pub struct SplitSequences {
     pub hap2: Option<Vec<String>>,
     pub flag: Option<String>,
     pub outliers: Option<Vec<String>>,
+    /// Number of clusters found, only set by the DBSCAN strategy (None for Ward).
+    /// Surfaced in the VCF as NCLUSTERS so complex multi-population loci can be flagged.
+    pub n_clusters: Option<usize>,
 }
 
 pub fn split(
@@ -116,10 +119,8 @@ pub fn split(
     // the fraction-based threshold is capped by `support` so that copy-number gains (which inflate
     // the total read count, and thus the fraction-based threshold) cannot swallow a genuine minority
     // allele that is supported by at least `support` reads
-    let min_cluster_size = max(
-        ((insertions.len() as f32 * min_haplotype_fraction) as usize).min(support),
-        1,
-    );
+    let min_cluster_size =
+        max(((insertions.len() as f32 * min_haplotype_fraction) as usize).min(support), 1);
     debug!(
         "{repeat}: Minimum cluster size: {} reads (min of {}% of {} total and support {})",
         min_cluster_size,
@@ -219,6 +220,7 @@ pub fn split(
                 } else {
                     None
                 },
+                n_clusters: None,
             }
         }
         1 => {
@@ -237,6 +239,7 @@ pub fn split(
                 } else {
                     None
                 },
+                n_clusters: None,
             }
         }
         2 => {
@@ -263,6 +266,7 @@ pub fn split(
                 } else {
                     None
                 },
+                n_clusters: None,
             }
         }
         _ => {
@@ -272,6 +276,111 @@ pub fn split(
     }
 }
 
+/// Alternative phasing strategy to [`split`] (Ward over length-weighted
+/// Levenshtein). Builds k-mer composition feature vectors (length-invariant)
+/// and clusters them with DBSCAN. Noise points and any clusters beyond the two
+/// largest become outliers. Returns `n_clusters` so complex multi-population
+/// loci can be flagged in the VCF.
+pub fn split_dbscan(
+    insertions: &[String],
+    repeat: &crate::repeats::RepeatInterval,
+    check_outliers: bool,
+    support: usize,
+    eps: f64,
+    length_weight: f64,
+) -> SplitSequences {
+    use crate::features::{self, KmerTable};
+
+    // Auto-detect the repeat period (k) from the median-length insertion.
+    let mut by_len: Vec<&String> = insertions.iter().collect();
+    by_len.sort_by_key(|s| s.len());
+    let median_seq = by_len[by_len.len() / 2];
+    let k = features::detect_period(median_seq.as_bytes(), features::AUTO_K_MAX)
+        .unwrap_or(features::AUTO_K_DEFAULT);
+    let table = KmerTable::new(k);
+    let max_len = insertions.iter().map(|s| s.len()).max().unwrap_or(0);
+
+    let points: Vec<Vec<f64>> = insertions
+        .iter()
+        .map(|s| features::build_feature_vector(s, &table, k, max_len, length_weight))
+        .collect();
+
+    // min_samples = support, so every DBSCAN cluster is backed by >= support reads.
+    let labels = crate::dbscan::cluster(&points, eps, support);
+
+    let mut clusters: HashMap<usize, Vec<String>> = HashMap::new();
+    let mut noise: Vec<String> = Vec::new();
+    for (seq, label) in insertions.iter().zip(labels.iter()) {
+        match label {
+            Some(cid) => clusters.entry(*cid).or_default().push(seq.clone()),
+            None => noise.push(seq.clone()),
+        }
+    }
+    let n_clusters = Some(clusters.len());
+    debug!(
+        "{repeat}: DBSCAN (k={k}, eps={eps}, min_samples={support}) found {} clusters and {} noise reads",
+        clusters.len(),
+        noise.len()
+    );
+
+    // Sort clusters by size (desc); tie-break on total bp then id for determinism.
+    let mut sized: Vec<(usize, Vec<String>)> = clusters.into_iter().collect();
+    sized.sort_by(|a, b| {
+        b.1.len()
+            .cmp(&a.1.len())
+            .then_with(|| {
+                b.1.iter()
+                    .map(String::len)
+                    .sum::<usize>()
+                    .cmp(&a.1.iter().map(String::len).sum::<usize>())
+            })
+            .then_with(|| a.0.cmp(&b.0))
+    });
+
+    let outliers_opt =
+        |extra: Vec<String>| -> Option<Vec<String>> { check_outliers.then_some(extra).filter(|e| !e.is_empty()) };
+
+    match sized.len() {
+        0 => {
+            debug!("{repeat}: DBSCAN found no clusters; treating as homozygous (CLUSTERFAILURE)");
+            SplitSequences {
+                hap1: insertions.to_vec(),
+                hap2: None,
+                flag: Some("CLUSTERFAILURE".to_string()),
+                outliers: None,
+                n_clusters,
+            }
+        }
+        1 => {
+            let hap1 = sized.into_iter().next().unwrap().1;
+            SplitSequences {
+                hap1,
+                hap2: None,
+                flag: None,
+                outliers: outliers_opt(noise),
+                n_clusters,
+            }
+        }
+        _ => {
+            let mut it = sized.into_iter();
+            let hap1 = it.next().unwrap().1;
+            let hap2 = it.next().unwrap().1;
+            // any clusters beyond the two largest join the noise as outliers
+            let mut extra = noise;
+            for (_, c) in it {
+                extra.extend(c);
+            }
+            SplitSequences {
+                hap1,
+                hap2: Some(hap2),
+                flag: None,
+                outliers: outliers_opt(extra),
+                n_clusters,
+            }
+        }
+    }
+}
+
 fn find_roots(
     top_root: usize,
     cluster_to_subclusters: &HashMap<usize, (usize, usize)>,
@@ -391,12 +500,10 @@ mod tests {
     #[allow(unused_imports)]
     use super::*;
 
-    #[test]
-    fn test_split_expansion_minority_allele() {
-        // Regression test for chr15:34419425-34419451 from rr_NA99_170.log
-        // A real, length-variable expansion is present on only 5/44 reads (the rest
-        // are the short reference allele). The expansion reads must be recovered as a
-        // haplotype and must NOT all be discarded as length outliers.
+    /// The 44 insertions of chr15:34419425-34419451 (GOLGA8A) from rr_NA99_170.log:
+    /// 5 length-variable, CT/CCTT-rich expansion reads + 39 short TTTC reference reads.
+    #[allow(dead_code)]
+    fn golga8a_insertions() -> Vec<String> {
         let r43 = "TCTTTCTTTCTTTCCTTTCCTTTCCTTTCCTTTCCTTTCCTTCCTTCCCTTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTTTTTCTTTCTTTCTTT".to_string();
         let r2 = "TCTTTCTTTCTTCCTTTCCTTTCCTTTCCTTTCCTTTCCTTCCTTCCCTTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTTTCTTTCTTTCTTT".to_string();
         let r33 = "CTTTCTTTCCTTTCCTTTCCTTTCCTTTCCTTTCCTTCCTTCCCCCTGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAAAAAAAAAAAAAAAAAAAAAAAATCTCTCTCTCTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTTTCTTTCTTTCTTTCTTTCTTTC".to_string();
@@ -427,6 +534,43 @@ mod tests {
             insertions.push("TTTCTTTCTTTCTTTCTTTCTTTCTTTTTTTTTT".to_string()); // 34
         }
         assert_eq!(insertions.len(), 44);
+        insertions
+    }
+
+    #[test]
+    fn test_split_dbscan_captures_expansion() {
+        // The DBSCAN strategy on the GOLGA8A data: with length-invariant k-mer
+        // composition features, the length-variable expansion is recovered as a
+        // distinct haplotype (rather than being fragmented by length and lost to
+        // the outlier bin, as happens with the length-dominated Ward distance).
+        let insertions = golga8a_insertions();
+        let repeat = crate::repeats::RepeatInterval {
+            chrom: "chr15".to_string(),
+            start: 34419425,
+            end: 34419451,
+            created: None,
+        };
+        // default CLI params (eps=0.2, length_weight=1.0), min_samples = support = 2
+        let s = split_dbscan(&insertions, &repeat, true, 2, 0.2, 1.0);
+        let hap2 = s.hap2.expect("expected a heterozygous call");
+        // exactly the reference + expansion alleles -> two clusters
+        assert_eq!(s.n_clusters, Some(2));
+        // one haplotype is the short reference, the other captures the expansion
+        let med1 = find_median(&s.hap1);
+        let med2 = find_median(&hap2);
+        assert!(
+            med1.min(med2) < 50 && med1.max(med2) > 100,
+            "expected one reference (~27bp) and one expansion (>100bp) haplotype, got {med1} and {med2}"
+        );
+    }
+
+    #[test]
+    fn test_split_expansion_minority_allele() {
+        // Regression test for chr15:34419425-34419451 from rr_NA99_170.log
+        // A real, length-variable expansion is present on only 5/44 reads (the rest
+        // are the short reference allele). The expansion reads must be recovered as a
+        // haplotype and must NOT all be discarded as length outliers.
+        let mut insertions = golga8a_insertions();
 
         use rand::seq::SliceRandom;
         let mut rng = rand::rng();
@@ -449,8 +593,7 @@ mod tests {
             .hap2
             .expect("expected a heterozygous call with two haplotypes");
         // One of the two haplotypes should capture the expansion (median length >> reference)
-        let expansion_captured =
-            find_median(&splitseqs.hap1) > 100 || find_median(&hap2) > 100;
+        let expansion_captured = find_median(&splitseqs.hap1) > 100 || find_median(&hap2) > 100;
         let n_outliers = splitseqs.outliers.as_ref().map_or(0, |o| o.len());
         assert!(
             expansion_captured,
diff --git a/src/vcf.rs b/src/vcf.rs
index 8811a00..cf0747e 100644
--- a/src/vcf.rs
+++ b/src/vcf.rs
@@ -52,6 +52,7 @@ pub struct VCFRecord {
     pub score: (String, String),
     pub somatic_info_field: String,
     pub outliers: String,
+    pub n_clusters: String,
     pub time_taken: String,
     pub ps: Option<u32>, // phase set identifier
     pub flags: String,
@@ -64,6 +65,7 @@ impl VCFRecord {
         repeat_ref_sequence: String,
         all_insertions: Option<Vec<String>>,
         outlier_insertions: Option<Vec<String>>,
+        n_clusters: Option<usize>,
         repeat: &crate::repeats::RepeatInterval,
         ps: Option<u32>,
         flag: Vec<String>,
@@ -126,6 +128,11 @@ impl VCFRecord {
             None => "".to_string(),
         };
 
+        let n_clusters = match n_clusters {
+            Some(n) => format!(";NCLUSTERS={n}"),
+            None => "".to_string(),
+        };
+
         let flags = if flag.is_empty() {
             "".to_string()
         } else {
@@ -153,6 +160,7 @@ impl VCFRecord {
             score: (allele1.score, allele2.score),
             somatic_info_field,
             outliers,
+            n_clusters,
             time_taken,
             ps,
             flags,
@@ -193,6 +201,7 @@ impl VCFRecord {
             score: (".".to_string(), ".".to_string()),
             somatic_info_field: "".to_string(),
             outliers: "".to_string(),
+            n_clusters: "".to_string(),
             time_taken,
             ps: None,
             flags: flags_str,
@@ -228,6 +237,7 @@ impl VCFRecord {
             score: (".".to_string(), ".".to_string()),
             somatic_info_field: "".to_string(),
             outliers: "".to_string(),
+            n_clusters: "".to_string(),
             time_taken,
             ps: None,
             flags: "".to_string(),
@@ -272,6 +282,7 @@ impl VCFRecord {
             score: (".".to_string(), ".".to_string()),
             somatic_info_field,
             outliers: "".to_string(),
+            n_clusters: "".to_string(),
             time_taken,
             ps,
             flags: if flag.is_empty() {
@@ -426,7 +437,7 @@ impl fmt::Display for VCFRecord {
                 };
                 write!(
                     f,
-                    "{chrom}\t{start}\t.\t{ref}\t{alt}\t.\t.\t{flags}END={end};STDEV={sd1},{sd2}{somatic}{outliers}{time_taken}\t{FORMAT}\t{genotype1}|{genotype2}:{l1},{l2}:{fl1},{fl2}:{sup1},{sup2}:{score1},{score2}{ps}",
+                    "{chrom}\t{start}\t.\t{ref}\t{alt}\t.\t.\t{flags}END={end};STDEV={sd1},{sd2}{somatic}{outliers}{n_clusters}{time_taken}\t{FORMAT}\t{genotype1}|{genotype2}:{l1},{l2}:{fl1},{fl2}:{sup1},{sup2}:{score1},{score2}{ps}",
                     chrom = self.chrom,
                     start = self.start,
                     flags = self.flags,
@@ -441,6 +452,7 @@ impl fmt::Display for VCFRecord {
                     sd2 = self.std_dev.1,
                     somatic = self.somatic_info_field,
                     outliers = self.outliers,
+                    n_clusters = self.n_clusters,
                     time_taken = self.time_taken,
                     genotype1 = self.allele.0,
                     genotype2 = self.allele.1,
@@ -530,6 +542,9 @@ pub fn write_vcf_header(args: &Cli) {
     println!(
         r#"##INFO=<ID=OUTLIERS,Number=1,Type=String,Description="Outlier sequences much longer than the alleles">"#
     );
+    println!(
+        r#"##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)">"#
+    );
     println!(
         r#"##INFO=<ID=CLUSTERFAILURE,Number=0,Type=Flag,Description="If unphased input failed to cluster in two haplotype">"#
     );
@@ -581,6 +596,9 @@ fn test_write_vcf_header_from_bam() {
         unphased: true,
         find_outliers: false,
         min_haplotype_fraction: 0.1,
+        phasing_strategy: crate::PhasingStrategy::Ward,
+        dbscan_eps: 0.2,
+        dbscan_length_weight: 1.0,
         threads: 1,
         sample: None,
         haploid: Some(String::from("chr7")),
@@ -608,6 +626,9 @@ fn test_write_vcf_header_from_name() {
         unphased: true,
         find_outliers: false,
         min_haplotype_fraction: 0.1,
+        phasing_strategy: crate::PhasingStrategy::Ward,
+        dbscan_eps: 0.2,
+        dbscan_length_weight: 1.0,
         threads: 1,
         sample: Some("test_sample".to_string()),
         haploid: Some(String::from("chr7")),

From 0970a8e10ad821607893d4dc2373b75b35cb4e05 Mon Sep 17 00:00:00 2001
From: wdecoster <decosterwouter@gmail.com>
Date: Thu, 4 Jun 2026 11:42:30 +0200
Subject: [PATCH 4/8] log: print the full parsed arguments at startup

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>
---
 src/main.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/main.rs b/src/main.rs
index 273b842..39b0a9c 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -150,7 +150,7 @@ fn main() {
     if args.find_outliers && !args.unphased {
         warn!("--find-outliers is only effective with --unphased");
     }
-    info!("Collected arguments");
+    info!("Collected arguments: {args:?}");
     call::genotype_repeats(args);
 }
 

From 528c7272a284a92388f2136a45ac26a5fc96a8da Mon Sep 17 00:00:00 2001
From: wdecoster <decosterwouter@gmail.com>
Date: Thu, 4 Jun 2026 11:51:08 +0200
Subject: [PATCH 5/8] tune: DBSCAN defaults to eps=0.4, length_weight=0.3

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>
---
 src/main.rs             | 4 ++--
 src/phase_insertions.rs | 9 +++++----
 2 files changed, 7 insertions(+), 6 deletions(-)

diff --git a/src/main.rs b/src/main.rs
index 39b0a9c..67e2298 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -90,13 +90,13 @@ pub struct Cli {
 
     /// DBSCAN neighbourhood radius in normalized [0,1] feature space (only with --phasing-strategy dbscan).
     /// Smaller values split more readily. Tune together with --dbscan-length-weight.
-    #[arg(long, default_value_t = 0.2)]
+    #[arg(long, default_value_t = 0.4)]
     dbscan_eps: f64,
 
     /// Weight of the (normalized) length axis relative to k-mer composition for DBSCAN
     /// (only with --phasing-strategy dbscan). Lower values let length-variable expansions of the
     /// same motif cluster together; higher values better separate same-motif alleles that differ only in length.
-    #[arg(long, default_value_t = 1.0)]
+    #[arg(long, default_value_t = 0.3)]
     dbscan_length_weight: f64,
 
     /// comma-separated list of haploid (sex) chromosomes
diff --git a/src/phase_insertions.rs b/src/phase_insertions.rs
index 7bf6535..937200a 100644
--- a/src/phase_insertions.rs
+++ b/src/phase_insertions.rs
@@ -337,8 +337,9 @@ pub fn split_dbscan(
             .then_with(|| a.0.cmp(&b.0))
     });
 
-    let outliers_opt =
-        |extra: Vec<String>| -> Option<Vec<String>> { check_outliers.then_some(extra).filter(|e| !e.is_empty()) };
+    let outliers_opt = |extra: Vec<String>| -> Option<Vec<String>> {
+        check_outliers.then_some(extra).filter(|e| !e.is_empty())
+    };
 
     match sized.len() {
         0 => {
@@ -550,8 +551,8 @@ mod tests {
             end: 34419451,
             created: None,
         };
-        // default CLI params (eps=0.2, length_weight=1.0), min_samples = support = 2
-        let s = split_dbscan(&insertions, &repeat, true, 2, 0.2, 1.0);
+        // default CLI params (eps=0.4, length_weight=0.3), min_samples = support = 2
+        let s = split_dbscan(&insertions, &repeat, true, 2, 0.4, 0.3);
         let hap2 = s.hap2.expect("expected a heterozygous call");
         // exactly the reference + expansion alleles -> two clusters
         assert_eq!(s.n_clusters, Some(2));

From 64b4eb21c545d21c5dd543ce3e5743e9d371e553 Mon Sep 17 00:00:00 2001
From: wdecoster <decosterwouter@gmail.com>
Date: Thu, 4 Jun 2026 11:54:29 +0200
Subject: [PATCH 6/8] docs: document --phasing-strategy / DBSCAN options and
 NCLUSTERS

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>
---
 README.md | 10 +++++++++-
 1 file changed, 9 insertions(+), 1 deletion(-)

diff --git a/README.md b/README.md
index a5d6c05..c11ff99 100644
--- a/README.md
+++ b/README.md
@@ -51,9 +51,13 @@ 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 <STRATEGY>  How to split unphased reads into haplotypes: 'ward' or 'dbscan' (only with --unphased) [default: ward]
+        --dbscan-eps <EPS>             DBSCAN neighbourhood radius in normalized [0,1] feature space (only with --phasing-strategy dbscan) [default: 0.4]
+        --dbscan-length-weight <W>     Weight of the length axis vs k-mer composition for DBSCAN (only with --phasing-strategy dbscan) [default: 0.3]
         --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,6 +69,9 @@ 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. Two strategies are available via `--phasing-strategy`:
+  - `ward` (default): hierarchical (Ward) clustering on a length-weighted Levenshtein distance. Robust for the common case where alleles differ mainly in length.
+  - `dbscan` (experimental): DBSCAN on length-invariant k-mer composition feature vectors. Better at keeping a length-variable expansion together as one allele (where the length-based distance tends to fragment it), at the cost of needing the reference and expanded alleles to differ in composition. 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). Tune with `--dbscan-eps` (larger = fewer, looser clusters) and `--dbscan-length-weight` (lower = composition matters more than length).
 - 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
@@ -81,6 +88,7 @@ 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">
 ##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">

From f46b60acc1ece9665d531a42176c924943466269 Mon Sep 17 00:00:00 2001
From: wdecoster <decosterwouter@gmail.com>
Date: Sun, 14 Jun 2026 10:44:41 +0200
Subject: [PATCH 7/8] style: rustfmt the seeded single-read consensus selection

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>
---
 src/consensus.rs | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/consensus.rs b/src/consensus.rs
index ea71e13..27b1acd 100644
--- a/src/consensus.rs
+++ b/src/consensus.rs
@@ -56,7 +56,10 @@ pub fn consensus(
         Consensus { seq: None, support: num_reads, std_dev, score: -1 }
     } 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 StdRng::seed_from_u64(DOWNSAMPLE_SEED)).unwrap();
+        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 }
     } else {
         // if there are more than <consensus_reads> reads, downsample before taking the consensus

From ac83a9711b0b4e32db28ce0233e9ffc727aa60d8 Mon Sep 17 00:00:00 2001
From: wdecoster <decosterwouter@gmail.com>
Date: Tue, 23 Jun 2026 11:20:19 +0200
Subject: [PATCH 8/8] feat: expansion/imprecision QC flags, MRL field, and
 --phasing both

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>
---
 Cargo.lock              |   2 +-
 Cargo.toml              |   2 +-
 README.md               |  38 +++++--
 src/consensus.rs        | 122 +++++++++++++++++-----
 src/dbscan.rs           |   2 +-
 src/genotype.rs         | 222 +++++++++++++++++++++++++++++++++++++---
 src/main.rs             |  18 ++--
 src/phase_insertions.rs |  23 +++++
 src/vcf.rs              |  64 ++++++++++--
 9 files changed, 422 insertions(+), 71 deletions(-)

diff --git a/Cargo.lock b/Cargo.lock
index d1950fd..9f67b74 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -4,7 +4,7 @@ version = 4
 
 [[package]]
 name = "STRdust"
-version = "0.18.0"
+version = "0.19.0"
 dependencies = [
  "bio",
  "clap",
diff --git a/Cargo.toml b/Cargo.toml
index 96bf9ab..3f8ecfa 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "STRdust"
-version = "0.18.0"
+version = "0.19.0"
 edition = "2024"
 
 
diff --git a/README.md b/README.md
index c11ff99..a38fda6 100644
--- a/README.md
+++ b/README.md
@@ -55,9 +55,7 @@ OPTIONS:
         --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 <STRATEGY>  How to split unphased reads into haplotypes: 'ward' or 'dbscan' (only with --unphased) [default: ward]
-        --dbscan-eps <EPS>             DBSCAN neighbourhood radius in normalized [0,1] feature space (only with --phasing-strategy dbscan) [default: 0.4]
-        --dbscan-length-weight <W>     Weight of the length axis vs k-mer composition for DBSCAN (only with --phasing-strategy dbscan) [default: 0.3]
+        --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
@@ -69,17 +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. Two strategies are available via `--phasing-strategy`:
-  - `ward` (default): hierarchical (Ward) clustering on a length-weighted Levenshtein distance. Robust for the common case where alleles differ mainly in length.
-  - `dbscan` (experimental): DBSCAN on length-invariant k-mer composition feature vectors. Better at keeping a length-variable expansion together as one allele (where the length-based distance tends to fragment it), at the cost of needing the reference and expanded alleles to differ in composition. 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). Tune with `--dbscan-eps` (larger = fewer, looser clusters) and `--dbscan-length-weight` (lower = composition matters more than length).
+- 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
@@ -90,17 +95,36 @@ Example output:
 ##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
diff --git a/src/consensus.rs b/src/consensus.rs
index 27b1acd..5db957d 100644
--- a/src/consensus.rs
+++ b/src/consensus.rs
@@ -12,17 +12,38 @@ use std::fmt;
 /// 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,
+        }
     }
 }
 
@@ -46,21 +67,32 @@ 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 StdRng::seed_from_u64(DOWNSAMPLE_SEED))
             .unwrap();
-        Consensus { seq: Some(seq.to_string()), support: num_reads, std_dev, score: 0 }
+        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
@@ -96,45 +128,53 @@ 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>>()
     }
 }
 
@@ -142,6 +182,36 @@ fn remove_outliers<'a>(
 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
diff --git a/src/dbscan.rs b/src/dbscan.rs
index d92b5e5..f0fce62 100644
--- a/src/dbscan.rs
+++ b/src/dbscan.rs
@@ -1,4 +1,4 @@
-//! DBSCAN clustering for the optional `--phasing-strategy dbscan` path.
+//! DBSCAN clustering for the optional `--phasing dbscan` path.
 //!
 //! Adapted from the `trout` cohort-outlier tool, which used a custom DBSCAN
 //! over a precomputed CSR adjacency for speed (one pairwise distance per pair,
diff --git a/src/genotype.rs b/src/genotype.rs
index c93e794..977a70b 100644
--- a/src/genotype.rs
+++ b/src/genotype.rs
@@ -10,6 +10,121 @@ use std::sync::OnceLock;
 // Compile regex once and reuse across all CS tag parsing
 static CS_REGEX: OnceLock<Regex> = OnceLock::new();
 
+/// A read counts towards EXPANSION_OUTLIER when it is at least this many times
+/// longer than the longer called allele. These capture large expansions that fell
+/// to noise/outliers or were removed as length outliers, and so were missed by both
+/// alleles. Hardcoded for now (default 2 reads, 2x longer); can be promoted to CLI
+/// arguments later if cohort experience warrants tuning.
+const EXPANSION_OUTLIER_MIN_READS: usize = 2;
+const EXPANSION_OUTLIER_RATIO: f64 = 2.0;
+
+/// In `--phasing both` QC mode, the reported (Ward) longer-allele length and the
+/// DBSCAN longer-allele length are flagged DISCORDANT_LENGTH when their ratio
+/// exceeds this. Hardcoded at 2x and deliberately fires liberally: we prioritise
+/// sensitivity, so over-flagging for review is acceptable.
+const DISCORDANT_LENGTH_RATIO: f64 = 2.0;
+
+/// Decide whether the EXPANSION_OUTLIER flag applies: at least
+/// [`EXPANSION_OUTLIER_MIN_READS`] insertions exceed [`EXPANSION_OUTLIER_RATIO`]
+/// times the longer called allele (the longest consensus length).
+///
+/// Only meaningful for `--unphased` data, where all spanning reads are pooled
+/// before clustering; it is therefore only ever called from the unphased path.
+fn expansion_outlier_flagged(
+    insertions: &[String],
+    consenses: &[crate::consensus::Consensus],
+) -> bool {
+    let longer_allele = longer_allele_len(consenses);
+    if longer_allele == 0 {
+        return false;
+    }
+    let threshold = EXPANSION_OUTLIER_RATIO * longer_allele as f64;
+    let n_big = insertions
+        .iter()
+        .filter(|s| s.len() as f64 > threshold)
+        .count();
+    n_big >= EXPANSION_OUTLIER_MIN_READS
+}
+
+/// Length of the longest called allele's consensus, in bases (0 if none called).
+fn longer_allele_len(consenses: &[crate::consensus::Consensus]) -> usize {
+    consenses
+        .iter()
+        .filter_map(|c| c.seq.as_ref().map(|s| s.len()))
+        .max()
+        .unwrap_or(0)
+}
+
+/// QC comparison for `--phasing both`: run DBSCAN alongside the reported Ward call,
+/// build its consensus alleles, and compare the longer allele of each method.
+/// Returns `(discordant, dbscan_rb)` where `discordant` is true when the longer
+/// allele lengths differ by more than [`DISCORDANT_LENGTH_RATIO`], and `dbscan_rb`
+/// is the DBSCAN allele lengths relative to the reference (for the DBSCAN_RB INFO field).
+fn dbscan_qc_comparison(
+    insertions: &[String],
+    repeat: &crate::repeats::RepeatInterval,
+    args: &Cli,
+    ward_consenses: &[crate::consensus::Consensus],
+) -> (bool, String) {
+    let dbscan = crate::phase_insertions::split_dbscan(
+        insertions,
+        repeat,
+        false,
+        args.support,
+        crate::phase_insertions::DBSCAN_EPS,
+        crate::phase_insertions::DBSCAN_LENGTH_WEIGHT,
+    );
+    let mut dbscan_consenses: Vec<crate::consensus::Consensus> = vec![];
+    match dbscan.hap2 {
+        Some(hap2) => {
+            dbscan_consenses.push(crate::consensus::consensus(
+                &dbscan.hap1,
+                args.support,
+                args.consensus_reads,
+                repeat,
+            ));
+            dbscan_consenses.push(crate::consensus::consensus(
+                &hap2,
+                args.support,
+                args.consensus_reads,
+                repeat,
+            ));
+        }
+        None => {
+            let c = crate::consensus::consensus(
+                &dbscan.hap1,
+                args.support,
+                args.consensus_reads,
+                repeat,
+            );
+            dbscan_consenses.push(c.clone());
+            dbscan_consenses.push(c);
+        }
+    }
+
+    let ref_len = (repeat.end - repeat.start) as i32;
+    let dbscan_rb = dbscan_consenses
+        .iter()
+        .map(|c| match &c.seq {
+            Some(s) => (s.len() as i32 - ref_len).to_string(),
+            None => ".".to_string(),
+        })
+        .collect::<Vec<String>>()
+        .join(",");
+
+    // Compare on full consensus lengths (always non-negative), so the 2x ratio is well-defined.
+    let discordant =
+        length_discordant(longer_allele_len(ward_consenses), longer_allele_len(&dbscan_consenses));
+    (discordant, dbscan_rb)
+}
+
+/// True when two allele lengths differ by more than [`DISCORDANT_LENGTH_RATIO`].
+/// A zero length on one side is discordant whenever the other side is non-zero.
+fn length_discordant(a: usize, b: usize) -> bool {
+    let (lo, hi) = (a.min(b), a.max(b));
+    hi as f64 > DISCORDANT_LENGTH_RATIO * lo as f64
+}
+
 /// Result of checking reads for quick reference detection
 enum ReadCheckResult {
     /// No reads found - return missing genotype immediately
@@ -211,6 +326,9 @@ pub fn genotype_with_extracted_reads(
     };
     // number of clusters found by the DBSCAN phasing strategy (None for Ward / phased input)
     let mut n_clusters: Option<usize> = None;
+    // DBSCAN allele lengths for the DBSCAN_RB INFO field, only set in `--phasing both` QC mode
+    // when the DBSCAN call is substantially length-discordant with the reported Ward call.
+    let mut dbscan_rb: Option<String> = None;
 
     // alignments can be extracted in an unphased manner, if the chromosome is --haploid or the --unphased is set
     let unphased = (args.haploid.is_some()
@@ -286,8 +404,16 @@ pub fn genotype_with_extracted_reads(
                 repeat,
                 args.find_outliers,
                 args.support,
-                args.dbscan_eps,
-                args.dbscan_length_weight,
+                crate::phase_insertions::DBSCAN_EPS,
+                crate::phase_insertions::DBSCAN_LENGTH_WEIGHT,
+            ),
+            // 'both' reports the Ward call; DBSCAN is run separately below for QC discordance.
+            crate::PhasingStrategy::Both => crate::phase_insertions::split(
+                &insertions,
+                repeat,
+                args.find_outliers,
+                args.min_haplotype_fraction,
+                args.support,
             ),
         };
         match phased.hap2 {
@@ -334,6 +460,17 @@ pub fn genotype_with_extracted_reads(
         }
         n_clusters = phased.n_clusters;
         outliers = phased.outliers;
+        if expansion_outlier_flagged(&insertions, &consenses) {
+            flags.push("EXPANSION_OUTLIER".to_string());
+        }
+        if matches!(args.phasing_strategy, crate::PhasingStrategy::Both) {
+            let (discordant, dbscan_rb_str) =
+                dbscan_qc_comparison(&insertions, repeat, args, &consenses);
+            if discordant {
+                flags.push("DISCORDANT_LENGTH".to_string());
+                dbscan_rb = Some(dbscan_rb_str);
+            }
+        }
     } else {
         // Phased
         let seq1 = &reads.phase1;
@@ -385,6 +522,7 @@ pub fn genotype_with_extracted_reads(
         all_insertions,
         outliers,
         n_clusters,
+        dbscan_rb,
         repeat,
         reads.ps,
         flags,
@@ -498,6 +636,9 @@ fn genotype_repeat(
     };
     // number of clusters found by the DBSCAN phasing strategy (None for Ward / phased input)
     let mut n_clusters: Option<usize> = None;
+    // DBSCAN allele lengths for the DBSCAN_RB INFO field, only set in `--phasing both` QC mode
+    // when the DBSCAN call is substantially length-discordant with the reported Ward call.
+    let mut dbscan_rb: Option<String> = None;
 
     // The rest of the function has three mutually exclusive options from here.
     // Either the reads are from a haploid chromosome, unphased or phased by a tool like WhatsHap/hiphase/...
@@ -581,8 +722,16 @@ fn genotype_repeat(
                 repeat,
                 args.find_outliers,
                 args.support,
-                args.dbscan_eps,
-                args.dbscan_length_weight,
+                crate::phase_insertions::DBSCAN_EPS,
+                crate::phase_insertions::DBSCAN_LENGTH_WEIGHT,
+            ),
+            // 'both' reports the Ward call; DBSCAN is run separately below for QC discordance.
+            crate::PhasingStrategy::Both => crate::phase_insertions::split(
+                &insertions,
+                repeat,
+                args.find_outliers,
+                args.min_haplotype_fraction,
+                args.support,
             ),
         };
         match phased.hap2 {
@@ -634,6 +783,17 @@ fn genotype_repeat(
         {
             outliers_vec.push(outliers_found.join(","));
         }
+        if expansion_outlier_flagged(&insertions, &consenses) {
+            flags.push("EXPANSION_OUTLIER".to_string());
+        }
+        if matches!(args.phasing_strategy, crate::PhasingStrategy::Both) {
+            let (discordant, dbscan_rb_str) =
+                dbscan_qc_comparison(&insertions, repeat, args, &consenses);
+            if discordant {
+                flags.push("DISCORDANT_LENGTH".to_string());
+                dbscan_rb = Some(dbscan_rb_str);
+            }
+        }
     } else {
         // input alignments are already phased
         for (phase, seq) in [(1, &reads.phase1), (2, &reads.phase2)] {
@@ -665,6 +825,7 @@ fn genotype_repeat(
         all_insertions,
         outliers,
         n_clusters,
+        dbscan_rb,
         repeat,
         reads.ps,
         flags,
@@ -776,6 +937,49 @@ fn parse_cs(
 #[cfg(test)]
 mod tests {
     use super::*;
+    use crate::consensus::Consensus;
+
+    fn consensus_of_len(len: usize) -> Consensus {
+        Consensus { seq: Some("A".repeat(len)), ..Default::default() }
+    }
+
+    #[test]
+    fn test_expansion_outlier_flagged() {
+        // two called alleles ~30 bp; two reads >60 bp (2x the longer allele) -> flagged
+        let consenses = vec![consensus_of_len(30), consensus_of_len(30)];
+        let insertions = vec![
+            "A".repeat(28),
+            "A".repeat(30),
+            "A".repeat(5483),
+            "A".repeat(3000),
+        ];
+        assert!(expansion_outlier_flagged(&insertions, &consenses));
+    }
+
+    #[test]
+    fn test_expansion_outlier_single_read_not_flagged() {
+        // only one read exceeds 2x the longer allele -> below the 2-read threshold (artifact case)
+        let consenses = vec![consensus_of_len(30), consensus_of_len(30)];
+        let insertions = vec!["A".repeat(28), "A".repeat(30), "A".repeat(474)];
+        assert!(!expansion_outlier_flagged(&insertions, &consenses));
+    }
+
+    #[test]
+    fn test_length_discordant() {
+        assert!(length_discordant(100, 250)); // 2.5x -> discordant
+        assert!(length_discordant(0, 50)); // one allele uncalled, other long -> discordant
+        assert!(!length_discordant(100, 150)); // 1.5x -> concordant
+        assert!(!length_discordant(0, 0)); // both uncalled -> not discordant
+        assert!(!length_discordant(100, 200)); // exactly 2x -> not over the threshold
+    }
+
+    #[test]
+    fn test_expansion_outlier_not_flagged_when_alleles_long() {
+        // the expansion is already captured as an allele; no read is 2x longer -> not flagged
+        let consenses = vec![consensus_of_len(30), consensus_of_len(3000)];
+        let insertions = vec!["A".repeat(30), "A".repeat(2900), "A".repeat(3100)];
+        assert!(!expansion_outlier_flagged(&insertions, &consenses));
+    }
 
     #[test]
     fn test_parse_cs() {
@@ -835,8 +1039,6 @@ mod tests {
             find_outliers: false,
             min_haplotype_fraction: 0.1,
             phasing_strategy: crate::PhasingStrategy::Ward,
-            dbscan_eps: 0.2,
-            dbscan_length_weight: 1.0,
             threads: 1,
             sample: None,
             haploid: None,
@@ -873,8 +1075,6 @@ mod tests {
             find_outliers: false,
             min_haplotype_fraction: 0.1,
             phasing_strategy: crate::PhasingStrategy::Ward,
-            dbscan_eps: 0.2,
-            dbscan_length_weight: 1.0,
             threads: 1,
             sample: None,
             haploid: Some(String::from("chr7")),
@@ -905,8 +1105,6 @@ mod tests {
             find_outliers: false,
             min_haplotype_fraction: 0.1,
             phasing_strategy: crate::PhasingStrategy::Ward,
-            dbscan_eps: 0.2,
-            dbscan_length_weight: 1.0,
             threads: 1,
             sample: None,
             haploid: None,
@@ -943,8 +1141,6 @@ mod tests {
             find_outliers: false,
             min_haplotype_fraction: 0.1,
             phasing_strategy: crate::PhasingStrategy::Ward,
-            dbscan_eps: 0.2,
-            dbscan_length_weight: 1.0,
             threads: 1,
             sample: None,
             haploid: None,
@@ -983,8 +1179,6 @@ mod tests {
             find_outliers: false,
             min_haplotype_fraction: 0.1,
             phasing_strategy: crate::PhasingStrategy::Ward,
-            dbscan_eps: 0.2,
-            dbscan_length_weight: 1.0,
             threads: 1,
             sample: None,
             haploid: None,
diff --git a/src/main.rs b/src/main.rs
index 67e2298..9614dfc 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -24,6 +24,9 @@ pub enum PhasingStrategy {
     Ward,
     /// k-mer composition feature vectors + DBSCAN.
     Dbscan,
+    /// QC mode: report the Ward call but additionally run DBSCAN and flag substantial
+    /// length discordance (DISCORDANT_LENGTH / DBSCAN_RB) for review.
+    Both,
 }
 
 // The arguments end up in the Cli struct
@@ -85,20 +88,11 @@ pub struct Cli {
     /// Strategy for splitting unphased reads into haplotypes (only with --unphased).
     /// 'ward': length-weighted Levenshtein + hierarchical clustering (default).
     /// 'dbscan': k-mer composition features + DBSCAN (experimental, robust to length-variable expansions).
-    #[arg(long, value_enum, default_value_t = PhasingStrategy::Ward)]
+    /// 'both': QC mode that reports the Ward call but also runs DBSCAN and flags substantial
+    /// length discordance (DISCORDANT_LENGTH / DBSCAN_RB) for review.
+    #[arg(long = "phasing", value_name = "STRATEGY", value_enum, default_value_t = PhasingStrategy::Ward)]
     phasing_strategy: PhasingStrategy,
 
-    /// DBSCAN neighbourhood radius in normalized [0,1] feature space (only with --phasing-strategy dbscan).
-    /// Smaller values split more readily. Tune together with --dbscan-length-weight.
-    #[arg(long, default_value_t = 0.4)]
-    dbscan_eps: f64,
-
-    /// Weight of the (normalized) length axis relative to k-mer composition for DBSCAN
-    /// (only with --phasing-strategy dbscan). Lower values let length-variable expansions of the
-    /// same motif cluster together; higher values better separate same-motif alleles that differ only in length.
-    #[arg(long, default_value_t = 0.3)]
-    dbscan_length_weight: f64,
-
     /// comma-separated list of haploid (sex) chromosomes
     #[arg(long)]
     haploid: Option<String>,
diff --git a/src/phase_insertions.rs b/src/phase_insertions.rs
index 937200a..3eaf253 100644
--- a/src/phase_insertions.rs
+++ b/src/phase_insertions.rs
@@ -3,6 +3,29 @@ use kodama::{Method, linkage};
 use log::{Level, debug, error, log_enabled};
 use std::{cmp::max, collections::HashMap};
 
+/// DBSCAN neighbourhood radius in the normalized feature space (squared-euclidean over
+/// `[length_weight * normalized_length, canonical_kmer_freqs...]`).
+///
+/// Tuning intuition (hardcoded for now; the `split_dbscan` signature still accepts a value
+/// so it can be swept in tests or promoted to a CLI argument on request):
+///
+/// - smaller eps -> tighter, more clusters: more splitting and more reads dropped to noise, which tends to *undercall* length-variable expansions (the long tail fragments off);
+/// - larger eps -> fewer, looser clusters: distinct alleles can be *merged* into one.
+///
+/// Adjust in small steps (±0.1) and watch `NCLUSTERS`.
+pub const DBSCAN_EPS: f64 = 0.4;
+
+/// Weight of the (normalized) length axis relative to the length-invariant k-mer composition.
+///
+/// Tuning intuition:
+///
+/// - lower (→0): composition dominates, so same-motif reads cluster together regardless of length — best for recovering length-variable / mosaic expansions as a single allele;
+/// - higher (→1): length matters more, better separating two same-motif alleles that differ only in length (behaving more like the length-dominated Ward strategy).
+///
+/// The 0.3 default favours keeping length-variable expansions together, which is the reason
+/// to reach for DBSCAN in the first place.
+pub const DBSCAN_LENGTH_WEIGHT: f64 = 0.3;
+
 pub struct SplitSequences {
     pub hap1: Vec<String>,
     pub hap2: Option<Vec<String>>,
diff --git a/src/vcf.rs b/src/vcf.rs
index cf0747e..a920410 100644
--- a/src/vcf.rs
+++ b/src/vcf.rs
@@ -10,29 +10,36 @@ use std::io::Read;
 pub struct Allele {
     pub length: String, // length of the consensus sequence minus the length of the repeat sequence
     pub full_length: String, // length of the consensus sequence
+    pub median_length: String, // median read length of the cluster, relative to the reference
     pub support: String, // number of reads supporting the allele
     pub std_dev: String, // standard deviation of the repeat length
     pub score: String,  // consensus score in the poa graph
+    pub imprecise: bool, // cluster read-length spread exceeds the imprecision threshold
     pub seq: String,    // consensus sequence
 }
 
 impl Allele {
     pub fn from_consensus(consensus: Consensus, start: u32, end: u32) -> Allele {
+        let ref_len = (end - start) as i32;
         match consensus.seq {
             Some(seq) => Allele {
-                length: (seq.len() as i32 - ((end - start) as i32)).to_string(),
+                length: (seq.len() as i32 - ref_len).to_string(),
                 full_length: seq.len().to_string(),
+                median_length: (consensus.median_length as i32 - ref_len).to_string(),
                 support: consensus.support.to_string(),
                 std_dev: consensus.std_dev.to_string(),
                 score: consensus.score.to_string(),
+                imprecise: consensus.imprecise,
                 seq,
             },
             None => Allele {
                 length: ".".to_string(),
                 full_length: ".".to_string(),
+                median_length: ".".to_string(),
                 support: consensus.support.to_string(),
                 std_dev: ".".to_string(),
                 score: ".".to_string(),
+                imprecise: false,
                 seq: ".".to_string(),
             },
         }
@@ -47,12 +54,14 @@ pub struct VCFRecord {
     pub alt_seq: Option<String>,
     pub length: (String, String),
     pub full_length: (String, String),
+    pub median_length: (String, String),
     pub support: (String, String),
     pub std_dev: (String, String),
     pub score: (String, String),
     pub somatic_info_field: String,
     pub outliers: String,
     pub n_clusters: String,
+    pub dbscan_rb: String,
     pub time_taken: String,
     pub ps: Option<u32>, // phase set identifier
     pub flags: String,
@@ -66,9 +75,10 @@ impl VCFRecord {
         all_insertions: Option<Vec<String>>,
         outlier_insertions: Option<Vec<String>>,
         n_clusters: Option<usize>,
+        dbscan_rb: Option<String>,
         repeat: &crate::repeats::RepeatInterval,
         ps: Option<u32>,
-        flag: Vec<String>,
+        mut flag: Vec<String>,
         args: &Cli,
     ) -> VCFRecord {
         // since I use .pop() to format the two consensus sequences, the order is reversed
@@ -133,6 +143,17 @@ impl VCFRecord {
             None => "".to_string(),
         };
 
+        let dbscan_rb = match dbscan_rb {
+            Some(rb) => format!(";DBSCAN_RB={rb}"),
+            None => "".to_string(),
+        };
+
+        // Flag the locus when either called allele's cluster has a wide read-length spread:
+        // a single consensus length is then not representative (e.g. continuous/long-tailed loci).
+        if (allele1.imprecise || allele2.imprecise) && !flag.iter().any(|f| f == "IMPRECISE_LENGTH")
+        {
+            flag.push("IMPRECISE_LENGTH".to_string());
+        }
         let flags = if flag.is_empty() {
             "".to_string()
         } else {
@@ -155,12 +176,14 @@ impl VCFRecord {
             alt_seq: Some(alts),
             length: (allele1.length, allele2.length),
             full_length: (allele1.full_length, allele2.full_length),
+            median_length: (allele1.median_length, allele2.median_length),
             support: (allele1.support, allele2.support),
             std_dev: (allele1.std_dev, allele2.std_dev),
             score: (allele1.score, allele2.score),
             somatic_info_field,
             outliers,
             n_clusters,
+            dbscan_rb,
             time_taken,
             ps,
             flags,
@@ -196,12 +219,14 @@ impl VCFRecord {
             alt_seq: None,
             length: ("0".to_string(), "0".to_string()),
             full_length: (repeat_ref_seq.len().to_string(), repeat_ref_seq.len().to_string()),
+            median_length: ("0".to_string(), "0".to_string()),
             support: (".".to_string(), ".".to_string()),
             std_dev: (".".to_string(), ".".to_string()),
             score: (".".to_string(), ".".to_string()),
             somatic_info_field: "".to_string(),
             outliers: "".to_string(),
             n_clusters: "".to_string(),
+            dbscan_rb: "".to_string(),
             time_taken,
             ps: None,
             flags: flags_str,
@@ -232,12 +257,14 @@ impl VCFRecord {
             alt_seq: Some(".".to_string()),
             length: (".".to_string(), ".".to_string()),
             full_length: (".".to_string(), ".".to_string()),
+            median_length: (".".to_string(), ".".to_string()),
             support: (support, ".".to_string()),
             std_dev: (".".to_string(), ".".to_string()),
             score: (".".to_string(), ".".to_string()),
             somatic_info_field: "".to_string(),
             outliers: "".to_string(),
             n_clusters: "".to_string(),
+            dbscan_rb: "".to_string(),
             time_taken,
             ps: None,
             flags: "".to_string(),
@@ -277,12 +304,17 @@ impl VCFRecord {
             alt_seq: Some(seq.to_string()),
             length: ((seq.len() as u32 - (repeat.end - repeat.start)).to_string(), ".".to_string()),
             full_length: ((seq.len() as u32).to_string(), ".".to_string()),
+            median_length: (
+                (seq.len() as u32 - (repeat.end - repeat.start)).to_string(),
+                ".".to_string(),
+            ),
             support: ("1".to_string(), ".".to_string()),
             std_dev: (".".to_string(), ".".to_string()),
             score: (".".to_string(), ".".to_string()),
             somatic_info_field,
             outliers: "".to_string(),
             n_clusters: "".to_string(),
+            dbscan_rb: "".to_string(),
             time_taken,
             ps,
             flags: if flag.is_empty() {
@@ -432,12 +464,12 @@ impl fmt::Display for VCFRecord {
         match &self.alt_seq {
             Some(alts) => {
                 let (FORMAT, ps) = match self.ps {
-                    Some(ps) => ("GT:RB:FRB:SUP:SC:PS", format!(":{}", ps)),
-                    None => ("GT:RB:FRB:SUP:SC", "".to_string()),
+                    Some(ps) => ("GT:RB:FRB:MRL:SUP:SC:PS", format!(":{}", ps)),
+                    None => ("GT:RB:FRB:MRL:SUP:SC", "".to_string()),
                 };
                 write!(
                     f,
-                    "{chrom}\t{start}\t.\t{ref}\t{alt}\t.\t.\t{flags}END={end};STDEV={sd1},{sd2}{somatic}{outliers}{n_clusters}{time_taken}\t{FORMAT}\t{genotype1}|{genotype2}:{l1},{l2}:{fl1},{fl2}:{sup1},{sup2}:{score1},{score2}{ps}",
+                    "{chrom}\t{start}\t.\t{ref}\t{alt}\t.\t.\t{flags}END={end};STDEV={sd1},{sd2}{somatic}{outliers}{n_clusters}{dbscan_rb}{time_taken}\t{FORMAT}\t{genotype1}|{genotype2}:{l1},{l2}:{fl1},{fl2}:{ml1},{ml2}:{sup1},{sup2}:{score1},{score2}{ps}",
                     chrom = self.chrom,
                     start = self.start,
                     flags = self.flags,
@@ -448,11 +480,14 @@ impl fmt::Display for VCFRecord {
                     l2 = self.length.1,
                     fl1 = self.full_length.0,
                     fl2 = self.full_length.1,
+                    ml1 = self.median_length.0,
+                    ml2 = self.median_length.1,
                     sd1 = self.std_dev.0,
                     sd2 = self.std_dev.1,
                     somatic = self.somatic_info_field,
                     outliers = self.outliers,
                     n_clusters = self.n_clusters,
+                    dbscan_rb = self.dbscan_rb,
                     time_taken = self.time_taken,
                     genotype1 = self.allele.0,
                     genotype2 = self.allele.1,
@@ -548,6 +583,18 @@ pub fn write_vcf_header(args: &Cli) {
     println!(
         r#"##INFO=<ID=CLUSTERFAILURE,Number=0,Type=Flag,Description="If unphased input failed to cluster in two haplotype">"#
     );
+    println!(
+        r#"##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)">"#
+    );
+    println!(
+        r#"##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">"#
+    );
+    println!(
+        r#"##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)">"#
+    );
+    println!(
+        r#"##INFO=<ID=DBSCAN_RB,Number=2,Type=Integer,Description="With --phasing both: DBSCAN allele lengths relative to reference, reported when DISCORDANT_LENGTH is set">"#
+    );
     println!(
         r#"##INFO=<ID=QUICKREF,Number=0,Type=Flag,Description="Locus identified as homozygous reference via fast CIGAR check, alignment skipped">"#
     );
@@ -563,6 +610,9 @@ pub fn write_vcf_header(args: &Cli) {
     println!(
         r#"##FORMAT=<ID=FRB,Number=2,Type=Integer,Description="Full repeat length of the two alleles in bases">"#
     );
+    println!(
+        r#"##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">"#
+    );
     println!(r#"##FORMAT=<ID=PS,Number=1,Type=Integer,Description="Phase set identifier">"#);
     println!(r#"##FORMAT=<ID=SUP,Number=2,Type=Integer,Description="Read support per allele">"#);
     println!(r#"##FORMAT=<ID=SC,Number=2,Type=Integer,Description="Consensus score per allele">"#);
@@ -597,8 +647,6 @@ fn test_write_vcf_header_from_bam() {
         find_outliers: false,
         min_haplotype_fraction: 0.1,
         phasing_strategy: crate::PhasingStrategy::Ward,
-        dbscan_eps: 0.2,
-        dbscan_length_weight: 1.0,
         threads: 1,
         sample: None,
         haploid: Some(String::from("chr7")),
@@ -627,8 +675,6 @@ fn test_write_vcf_header_from_name() {
         find_outliers: false,
         min_haplotype_fraction: 0.1,
         phasing_strategy: crate::PhasingStrategy::Ward,
-        dbscan_eps: 0.2,
-        dbscan_length_weight: 1.0,
         threads: 1,
         sample: Some("test_sample".to_string()),
         haploid: Some(String::from("chr7")),