Tracking: DNA QC subcommand PR series (#152, #153, #154, #155, #156)

[BenjaminDEMAILLE]

Adds a rustqc dna subcommand: single-pass DNA (whole-genome and targeted) quality control, needing no gene annotation. Five stacked PRs, each independently buildable and tested.

Addresses #128, which asks for Qualimap bamqc mode alongside the existing rnaseq mode. PR #156 delivers it.

Review order

The branches are stacked, so each diff assumes the previous one. Reviewing in order keeps them small; each will shrink further as its predecessor merges.

PR	What it does	Tests
#152	Move assay-agnostic analyses into src/common/. Pure refactor, no behaviour change.	233
#153	rustqc dna: depth engine, mosdepth outputs, samtools trio, preseq	285
#154	Picard CollectWgsMetrics and CollectInsertSizeMetrics	306
#155	--targets mode, Picard CollectHsMetrics and CollectGcBiasMetrics	330
#156	Qualimap bamqc, plus the DNA documentation pages	353

No new dependencies; Cargo.toml and Cargo.lock are untouched.

Parity

Every output is validated against the real upstream tools on a committed fixture: a real public human chr22 slice from nf-core/test-datasets, duplicate-marked locally. tests/create_dna_test_data.sh regenerates inputs and reference outputs, pinning mosdepth 0.3.14, samtools 1.24, Picard 3.4.0 and Qualimap 2.3, and refusing to run against other versions. Fixtures total 496 kB.

Tool	Result
mosdepth	all six outputs identical, including the 1094-line global distribution
samtools flagstat, idxstats	identical
samtools stats	all 1889 data lines identical; RustQC writes its own # header by design
Picard CollectInsertSizeMetrics	identical, metrics row and all 170 histogram lines
Picard CollectWgsMetrics	identical on every column and all 251 histogram lines but two
Picard CollectGcBiasMetrics	both files identical
Picard CollectHsMetrics	all 58 computable columns identical
Qualimap genome_results.txt	identical on every line but four
Qualimap raw tables	three identical byte for byte, the rest within a tight tolerance

Known gaps

These are stated in the code, the docs and the PR bodies, not hidden behind a ?.

Picard's Monte Carlo simulation. HET_SNP_SENSITIVITY and HET_SNP_Q come from TheoreticalSensitivity, whose draws would have to be reproduced bit for bit. HS_PENALTY_10X through HS_PENALTY_100X and FOLD_80_BASE_PENALTY derive from it. All are written the way Picard writes a value it cannot compute.

CollectHsMetrics per-target GC. AT_DROPOUT and GC_DROPOUT need a per-target GC binning that is not implemented. The genome-wide GC machinery from #155 is most of what it needs.

Qualimap residuals. mean mapping quality and std coverageData differ in the fourth decimal; 393 of the 397 windows match exactly and the four that do not differ by at most 0.053, consistent with Qualimap accumulating per window. About five reference positions of 40001 sit one deeper here, which shows in 7 of roughly 590 coverage histogram bins and carries into the fractions derived from it. genome_fraction_coverage differs only in the last two digits of the double.

Qualimap homopolymer indels. Classified here by a run-of-four rule. Qualimap classifies against a reference context this does not reconstruct, and reports two polyC indels that no read-derived rule produces, because the deleted bases are not in the read. Reconstructing the local reference from the MD tag is the obvious way in.

Two Qualimap tables are not written. Its GC content distribution is computed over a 679-read subsample whose selection rule is undocumented; its duplication rate histogram uses a definition that is not a read-start-position count. Emitting tables under those names with different numbers seemed worse than leaving them out.

Things worth knowing before reading the code

Several upstream rules are counter-intuitive, and each is commented where it is implemented:

• Mate-overlap correction is usually the largest single effect on a coverage number. On this fixture it is the difference between 469875 and 247878 covered bases, exactly the gap between mosdepth --fast-mode and its default.
• mosdepth's distribution files always emit depths 0 through 300 whether or not any base sits there, only observed depths above that, and never the maximum when it lies above that range.
• Picard bins reverse-strand reads for GC by alignment_end - window_size, which is not the read's 5' end.
• CollectHsMetrics clips overlapping mates before applying the base quality floor; CollectWgsMetrics does the opposite. Their coverage figures are therefore not comparable to each other, and RustQC keeps a separate accumulator for each.
• Qualimap counts coverage with no filtering at all and no overlap correction, so it reports 16.77 where mosdepth reports 6.20 and CollectWgsMetrics reports 3.53. All three are right about different questions.
• Qualimap counts base composition in reference orientation while taking the clipped span that selects positions in sequencing orientation. Reproducing that inconsistency is what makes the composition figures agree.

The first four of these did not yield to inference from the tools' outputs and came from reading the upstream Java instead.

Follow-ups not in this series

• src/rna keeps re-export shims for everything moved in refactor: extract assay-agnostic analyses into src/common/ #152, so no published path breaks on 0.2.x. They can go at 1.0.
• -Q/--mapq defaults to 0 for dna and 30 for rna. Deliberate: 0 is mosdepth's default and matching it is what makes the depth figures agree. Worth a second opinion.

🤖 Generated with Claude Code