feat(dna): Picard CollectWgsMetrics and CollectInsertSizeMetrics

[BenjaminDEMAILLE]

Part of #157.

Third of the DNA stack. Stacked on #153, which is itself stacked on #152; review those first.

Adds the two core Picard metrics to rustqc dna, both computed in the same single pass as everything else.

Parity against Picard 3.4.0

Output	Result
insert_size_metrics.txt	identical, metrics row and all 170 histogram lines
wgs_metrics.txt	identical on every column and all 251 histogram lines, except two

The two exceptions are HET_SNP_SENSITIVITY and HET_SNP_Q, which come from Picard's TheoreticalSensitivity, a Monte Carlo simulation over the base quality and depth distributions. Matching it would mean reproducing its random number generator and sampling order bit for bit, which buys nothing for quality control. Both are written as ?, the marker Picard itself uses for a value it cannot compute, and the parity test permits a difference in exactly those two columns and nowhere else. FOLD_80/90/95_BASE_PENALTY are ? on both sides: Picard could not compute them on this data either.

How the semantics were established

Neither tool's rules were written from memory. Each was derived by reproducing Picard's own numbers on the fixture until every figure matched, and the two that mattered most both turned out to be non-obvious.

Insert size inclusion. Paired, not secondary, supplementary, duplicate or unmapped, mate mapped, positive TLEN so each pair counts once. Proper-pair is deliberately not required: requiring it drops one pair and shortens the maximum from 300 to 239 on this data. MEAN and STANDARD_DEVIATION are over the histogram trimmed to DEVIATIONS median absolute deviations either side of the median with the n-1 denominator, while MIN and MAX are over the untrimmed set. WIDTH_OF_XX_PERCENT grows a window symmetrically around the median until it covers the percentile and reports 2i+1; all eleven match.

The WGS exclusion model. Unmapped, secondary and supplementary records never enter the calculation at all. Every other record's reference-consuming bases form the denominator of all the PCT_EXC_* columns, 670989 bases on the fixture. Exclusions then apply in a fixed order, each counted against that same denominator: duplicate, low mapping quality and unpaired remove a whole read; low base quality and mate overlap remove single bases; depth beyond COVERAGE_CAP is counted as excess. What survives is the high quality coverage the histogram reports, and it reconciles exactly: 670989 − 201120 − 4933 − 217866 − 105814 = 141256.

SD_COVERAGE is the sample standard deviation over every base of the territory, uncovered ones included, not over the covered ones only.

A design point that held up

CollectWgsMetrics gets its own depth accumulator rather than a correction applied to the mosdepth one, because the two tools disagree on both which reads and which bases count. Keeping them separate was the design's stated reason for accepting two accumulators, and reproducing Picard's fractions confirmed it: there is no post-hoc correction that turns one into the other.

Behaviour

CollectWgsMetrics needs --reference to count the reference's non-N bases, which is GENOME_TERRITORY. Without one it is skipped with a warning rather than reported against a wrong denominator; insert size needs no reference and is still written. There is a test for exactly that.

Fixtures are generated with the JVM locale pinned to English. A French default locale writes 3,531312 where an English one writes 3.531312, which would make the fixtures depend on the machine that produced them. Picard's four-line preamble, a command line with absolute paths and a timestamp, is stripped.

Tests

306 green, up from 285. Twenty-four of them are DNA parity tests, at both the library and the binary level.

🤖 Generated with Claude Code