STRdust has no tolerance for read noise, and over-calls alternate alleles at short repeats

[wdecoster]

Found while investigating why QUICKREF almost never fires once its coordinate bug is fixed (#30). This looks like a better explanation for the benchmark result than anything found so far, and it is a bigger problem than QUICKREF.

The measurement

On 2000 catalog loci of a 30x ONT sample, with QUICKREF disabled so every locus is genotyped properly, STRdust calls 1652 of 1999 loci (83%) non-reference. The benchmark paper reports that 83.6% of loci are invariant (homozygous for the reference allele) in the R10 HPRC assemblies. Even allowing for a different sample and a 1.4 Mb catalog subset, that gap is far too large to be locus selection.

Of those non-reference calls, more than half are noise-scale:

largest |RB| per locus     count    cumulative
        0-1 bp              262       15.9%
        2 bp                263       31.8%
        3 bp                341       52.4%
        >= 4 bp             786      100.0%

The median REF length at these loci is 11 bp.

The mechanism

is_similar_to_ref (src/vcf.rs) decides whether an allele counts as reference using an edit-distance threshold of reference.len() / 20 — 5%, floor division. The comment in the code states the consequence plainly:

// For 12bp: threshold = 0, so only exact matches pass

With a median REF length of 11 bp, most loci require the consensus to match the reference exactly to be called 0/0. A POA consensus over ONT reads at a short STR essentially never does, so a one-base stutter becomes a 0/1 call.

The same noise-intolerance shows up in QUICKREF, which requires every read to differ from the reference by exactly zero. Sweeping a tolerance there does not rescue it, because a few bp is also the size of most real variation at these loci:

tolerance	QUICKREF fires	of which wrong
0	1/1999	0
3	26/1999	9
5	153/1999	99

Why this matches the published benchmark

From "A comprehensive assessment of tandem repeat genotyping methods for Nanopore long-read genomes" (Genome Biology, in press), which benchmarked v0.16.0:

• "For 0/0 genotypes, all tools except STRdust achieved >90% concordance" — and 0/0 loci are 83.6% of all loci.
• "STRdust tends to call more genotypes with an alternate allele."
• Mendelian consistency at 0/0: 99.5% against 99.8-99.99% for every other tool, "driven by a higher rate of motif-length-off calls at this class".
• "Mononucleotide repeats were the most difficult to accurately call, with STRdust below 40% consistency" — homopolymers are exactly where ONT indel noise is worst and where len()/20 is 0.
• Overall concordance 72.56% (R10), rising to 94.93% when an off-by-motif deviation is permitted.

Direction

That last figure is the most useful hint in the paper: the natural unit of tolerance is one motif, not a percentage of the interval length. STRdust currently does not know the motif at all — the catalog BED carries it in column 4 and STRdust ignores it, and src/motif.rs is an unimplemented stub.

Worth considering, roughly in order of effort:

1. Give is_similar_to_ref an absolute floor so a 1-2 bp difference at a short repeat does not force an alt call, rather than a percentage that rounds to zero.
2. Read the motif from the catalog and make both the reference-similarity test and the reported allele length motif-aware.
3. Model stutter/sequencing error properly, as Medaka Tandem (error model estimated from the reads) and LongTR (HipSTR-derived stutter model) do — the two tools that top the benchmark.

Should be measured with #23 rather than tuned against the alignment path.

[wdecoster]

The measurement here holds, but the diagnosis in this issue is wrong, and the fix it proposes would have masked the cause rather than removed it.

Benchmarked against the GIAB HG002 tandem-repeat truth set over the adotto catalog, 10,000–50,000 loci, replicated on seeds never used for tuning. False non-reference calls at loci whose truth is homozygous reference:

	false non-ref
as shipped	25.2%
--minlen 5	18.8%
--poa-gap-open 70	18.6%
after the two fixes below	7.9%

Two thirds of it was the POA consensus, not the threshold. bio::alignment::poa::Aligner::consensus() chooses its endpoint with max_by_key over a cumulative score to which every edge contributes a weight of at least 1. The score therefore rises strictly along every edge, so the argmax is always a sink: the consensus runs to the deepest point in the graph, not the best-supported one. A single read extending one base past the others creates a deeper sink and wins, with no vote taken.

Three things confirm it: the excess is 29:1 one-sided (longer, not shorter); it gets worse with read depth (14.9% at 0–4 reads, 32.6% at 15–19), which a consensus cannot do but an estimator tracking the maximum can; and MRL, the median read length of the same cluster, is 91.9% exact where the consensus is 71.2% — and is invariant across every scoring and filtering parameter tried. The reads were never wrong.

That also explains why this looked like a threshold problem: the noise-scale calls tabulated above are largely +1 and +2 artefacts of the consensus, so any tolerance would have hidden them.

Most of the rest was is_similar_to_ref, but in the opposite direction to the one proposed here. It was already too tolerant, not too strict: of 98 loci that a corrected consensus moved from variant to reference, 94 still had a called length differing from the reference — floor(len/20) with a strict < tolerated 2 edits at a median REF length of 77, so genuine 1–2 base alleles were reported 0/0. The lengths were correct in RB all along.

Adding tolerance would therefore have been sized against a length distribution distorted by the POA bug — fitting the threshold to the defect.

What shipped instead: a corrected consensus traversal (--poa-trim-fraction, +26 points of exact allele-length concordance on held-out seeds), --minlen restored to 3, and --priority {sensitive,balanced,precise} making the reference/variant decision an explicit, documented trade instead of an accidental one — the old rule sat at 52.1% recall / 83.0% precision, on no efficient frontier.

Leaving this open: the over-calling is much reduced but not gone, and the remaining is_similar_to_ref behaviour is now a deliberate setting rather than an accident.