Correction (real data, 2026-10-04): the reproduction below is synthetic, with all mass on the bin columns. In real trans, which is indexed by mean quality, wrong anchors land on sparsely populated columns rather than empty ones. So the model builds, in R too and identically, but is badly corrupted rather than all 1e-7: with 2,12,23,37 on NovaSeq 16S, 241 of 456 off-diagonal rates sit at the floor and 21 at the 0.25 clamp. The existing min/max warning did fire for that case. The fix in #266 adds a mass-based warning that catches it, including when both extremes sit on anchors.
Summary
binned_qual_errfun (our makeBinnedQualErrfun, src/error_models.rs) reads each transition's rate only at the anchor qualities and interpolates between adjacent anchors. When an anchor falls on a quality column with no observations, its rate is NaN, and every segment touching it is skipped. If every segment is skipped, all off-diagonal rates fall back to min_error_rate (1e-7). No warning is printed when both observed extremes happen to land on anchors.
R does not build this model: the same case reaches max(which(!is.na(pred))) on an all-NA vector and stops with a subscript error. So we silently produce a model that R refuses to build.
Reproduction
Data binned at Q2/11/25/37 (NovaSeq), with anchors that bracket it but miss the interior bins:
| anchors |
A→C rate at Q2 / Q11 / Q25 / Q37 |
2,11,25,37 (correct) |
9.8e-3 / 2.5e-3 / 1.1e-3 / 6.5e-4 |
2,12,23,37 |
1e-7 / 1e-7 / 1e-7 / 1e-7, and no warning |
The existing warning (#208) mirrors R's: it fires only when the observed minimum or maximum is not an anchor. Here Q2 and Q37 are both anchors, so it stays silent.
Why it matters
An error model with every error rate at 1e-7 treats almost every variant as real. It also feeds anything pinned to it, such as floor measurements and A/B arms. The repo's own example anchors, 2,12,23,37 in dev/concordance/run_illumina.sh and dev/run_screen_sweep.sh, match no instrument's bins and produce exactly this case on NovaSeq data. See the ITS2 re-measure tracking issue.
A partial version is also silent and R-consistent: with one empty interior anchor, the surviving segments are flat-extrapolated across the gap. That gives a usable but distorted model.
Proposed fix (to agree before implementing)
- Error when an off-diagonal transition has no finite estimate after interpolation. That matches R's effective behaviour (it stops), with an actionable message: which anchors have no observations, and that
summary --report shows the run's real bins.
- Always warn (stderr, not
--verbose) when an anchor inside the observed range has no observations, or when a populated quality column is not an anchor. R does not check this; the warning changes no output, and it catches the wrong-bins case directly.
- Regression tests for both: the all-collapse case must error, and the partial case must warn but fit unchanged.
Summary
binned_qual_errfun(ourmakeBinnedQualErrfun,src/error_models.rs) reads each transition's rate only at the anchor qualities and interpolates between adjacent anchors. When an anchor falls on a quality column with no observations, its rate isNaN, and every segment touching it is skipped. If every segment is skipped, all off-diagonal rates fall back tomin_error_rate(1e-7). No warning is printed when both observed extremes happen to land on anchors.R does not build this model: the same case reaches
max(which(!is.na(pred)))on an all-NAvector and stops with a subscript error. So we silently produce a model that R refuses to build.Reproduction
Data binned at Q2/11/25/37 (NovaSeq), with anchors that bracket it but miss the interior bins:
2,11,25,37(correct)2,12,23,37The existing warning (#208) mirrors R's: it fires only when the observed minimum or maximum is not an anchor. Here Q2 and Q37 are both anchors, so it stays silent.
Why it matters
An error model with every error rate at
1e-7treats almost every variant as real. It also feeds anything pinned to it, such as floor measurements and A/B arms. The repo's own example anchors,2,12,23,37indev/concordance/run_illumina.shanddev/run_screen_sweep.sh, match no instrument's bins and produce exactly this case on NovaSeq data. See the ITS2 re-measure tracking issue.A partial version is also silent and R-consistent: with one empty interior anchor, the surviving segments are flat-extrapolated across the gap. That gives a usable but distorted model.
Proposed fix (to agree before implementing)
summary --reportshows the run's real bins.--verbose) when an anchor inside the observed range has no observations, or when a populated quality column is not an anchor. R does not check this; the warning changes no output, and it catches the wrong-bins case directly.