feat(dna): targeted mode, CollectHsMetrics and CollectGcBiasMetrics

[BenjaminDEMAILLE]

Part of #157.

Fourth of the DNA stack. Stacked on #154, itself on #153 and #152.

Adds GC bias metrics, and makes --targets real: it switches the run into targeted mode and produces Picard-compatible hs_metrics.txt.

Parity against Picard 3.4.0

Output	Result
gc_bias.detail_metrics.txt (101 rows)	identical
gc_bias.summary_metrics.txt	identical
hs_metrics.txt	all 58 computable columns identical

Seven HS metrics columns are not computed: HET_SNP_SENSITIVITY, HET_SNP_Q, the six HS_PENALTY_xX levels and FOLD_80_BASE_PENALTY all derive from Picard's Monte Carlo theoretical sensitivity, and AT_DROPOUT and GC_DROPOUT from a per-target GC binning not implemented here. Each is written the way Picard writes a value it cannot compute, -1 or ?.

These rules came from Picard's source, not from its output

A first pass tried to infer both collectors from their outputs and failed. Reading GcBiasUtils, GcBiasMetricsCollector, TargetMetricsCollector and htsjdk's SAMUtils settled it. Five of the rules would not have been guessed:

GC windows slide over positions 1 to len - window_size - 1. Both ends are clipped, so a 40001 base reference gives 39900 windows of 100 bases rather than the 39902 a naive reading produces, and the GC value truncates rather than rounds.

A read is assigned to the window at its alignment start, except on the reverse strand, where it goes to alignment_end - window_size. That is not the read's 5' end, and no offset applied to either end reproduces the distribution. Only unmapped reads and reads with an empty sequence are skipped, so secondary alignments count; unmapped reads still count towards TOTAL_CLUSTERS.

HsMetrics clips overlapping mates before applying the base quality floor, where CollectWgsMetrics does the opposite. That single ordering is why the two collectors report different PCT_EXC_BASEQ and PCT_EXC_OVERLAP on the same file (0.003982 against 0.007352, 0.330968 against 0.324694), and why reusing the WGS model put ON_TARGET_BASES at 247070 against Picard's 245122.

Only the left-most mate of an overlapping pair is clipped, losing everything from its mate's alignment start onwards. Note also that the count Picard reports is of read bases, insertions included, while what is actually removed is every base at or past that reference position; the two coincide only for a gapless read.

htsjdk's MATCH_OR_MISMATCH is the M operator alone, so = and X lose their whole CIGAR element to clipping rather than a partial one.

Two more came from the fixtures: GC_NC_x_y is a mean weighted by each bin's window count rather than a plain average over bins, and the GC tables carry two trailing blank lines where the other Picard tables carry one.

Also here

HS_LIBRARY_SIZE solves the Lander-Waterman equation by bisection exactly as Picard's estimator does, down to the forty iterations and the starting bracket, and reproduces its 3807.

BED intervals are merged on load. Overlapping targets would otherwise inflate the territory and double-count on-target bases. BED is zero-based half-open and an interval list is one-based inclusive; the conversion happens only at the edges.

Tests

330 green, up from 306. Thirty-three are DNA parity tests, at both the library and the binary level, including one asserting the targeted outputs appear only when targets are given.

🤖 Generated with Claude Code