feat(coverage): bigWig tracks from the alignment pass - #166
Open
BenjaminDEMAILLE wants to merge 22 commits into
Open
feat(coverage): bigWig tracks from the alignment pass#166BenjaminDEMAILLE wants to merge 22 commits into
BenjaminDEMAILLE wants to merge 22 commits into
Conversation
These three modules carry no RNA-specific logic and are needed by the forthcoming dna subcommand. src/rna re-exports them so every existing crate::rna::... path and the published 0.2.x library surface keep working. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
bam_stat is read-level and needs no annotation, and the samtools stats, flagstat and idxstats writers consume its result type, so all four move together into src/common/. src/rna/rseqc re-exports them. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
BamStatAccum gathers the read-level counters behind bam_stat and the samtools writers. Its process_read takes only a record and a MAPQ cutoff, so it is assay-agnostic and the dna pipeline will drive the same struct. The merge_vec_arrays helper moves with it, being its only consumer. rna::rseqc::accumulators re-exports the type. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Also corrects the AGENTS.md claim that the crate has no lib.rs, which has been untrue since seqeralabs#101. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
A real public human chr22 slice from nf-core/test-datasets, duplicate-marked locally with samtools, plus mosdepth 0.3.14 and samtools 1.24 reference outputs. The generation script pins both tool versions and refuses to run against others, so fixtures and tool versions cannot drift apart. 380 kB in total, well inside the 10 MB fixture budget. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Shared options keep the same long name, short flag and RUSTQC_* environment variable as their rna counterparts. The deliberate differences: no --gtf, no --stranded, and --mapq defaults to 0 rather than 30 because that is mosdepth's default. run_dna is a stub for now; the pipeline lands in the following commits. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Config gains a `dna` block alongside `rna`, with mosdepth and samtools sub-sections and a reuse of the existing PreseqConfig. Shared settings (chromosome_prefix, chromosome_mapping, sample_name, flat_output) are declared on DnaConfig itself, mirroring RnaConfig, so the two pipelines can be configured independently in one file. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
DepthAccum records aligned blocks as increments in a delta array the length of the contig, then a prefix sum turns that into per-base depth in one linear pass. Filters and CIGAR handling reproduce mosdepth 0.3.14 outside fast mode: flags 1796 excluded, MAPQ floor applied, M/=/X cover the reference, D/N advance without covering, and I/S/H/P do not advance. Mate-overlap correction, the other half of mosdepth's default behaviour, lands in the next commit. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
mosdepth counts a base once when both mates of a pair cover it, unless --fast-mode is given. On the test dataset this is the difference between 469875 and 247878 total covered bases, so it is the dominant behaviour rather than an edge case. Pending mates are held in a map keyed by read name, indexed by the position the outstanding mate was announced at so entries that can never be claimed are evicted as the coordinate-ordered scan moves past them. A test asserts the map empties. Includes an engine-level parity check against the committed mosdepth fixture: total covered bases and maximum depth both match exactly. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Six writers plus the per-contig summarisation that feeds them: summary, global and region distributions, per-base runs, per-window means and per-window threshold counts. Compressed outputs are bgzf, matching mosdepth. The distribution emission rule was reverse-engineered from the fixtures and is the non-obvious part: every depth from 0 up to min(300, max) gets a row whether or not any base sits at it, above 300 only depths that occur and lie strictly below the maximum do. So the maximum gets a row when it falls inside the dense range and none when it does not. The global distribution tops out at 866 with a maximum of 867, while the region distribution does emit its maximum of 204. The region distribution is over windows and their rounded mean depth, not over bases. Parity tests drive the library directly and compare every mosdepth output against the committed fixtures: all eight match, including the 1094-line global distribution and the 721-interval per-base BED. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
One rayon worker per contig, each holding its own depth array, feeding a DepthAccum, a BamStatAccum and a PreseqAccum from the same record stream, so the alignment is read once. A separate pass over unmapped records feeds the counters flagstat and idxstats report. Workers run longest contig first and their number is bounded by --max-depth-workers, defaulting to a 4 GB budget divided by the largest contig, because each worker costs four bytes per base. Outputs land under mosdepth/, samtools/ and preseq/, or flat with --flat-output. Input without duplicate marks is rejected unless --skip-dup-check is passed. Also fixes the samtools stats header, which hardcoded "rustqc rna" and so labelled DNA output as RNA output. End-to-end parity tests run the binary and compare against the fixtures: all six mosdepth files match byte for byte, flagstat and idxstats match exactly, and all 1889 data lines of samtools stats match. The stats header differs by design, RustQC naming itself rather than reproducing samtools' version banner, so that comparison is on data lines. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
InputSummary gains an optional dna block carrying genome length, covered bases, mean, median and maximum coverage, the percentage of the reference at or above each requested threshold, and the duplicate rate. An input carries either the RNA fields or this one, never both. Coverage thresholds are a list of objects rather than a map so the requested order survives serialisation; a map keyed by the threshold would sort "10" before "5". CITATIONS.md for a dna run cites mosdepth, samtools and preseq, and none of the RNA-only tools. The header is now shared between both writers. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…n used Closes two gaps left open in this PR. The bgzf BED outputs now get a .csi index built through htslib's tbx_index_build, as mosdepth writes and as tabix needs to seek into them. CSI rather than TBI because CSI carries no 512 Mb coordinate ceiling. Indexes are not compared byte for byte: an index is binary metadata over the compressed blocks, and two writers answering the same queries need not produce the same bytes. The test asserts instead that a region query returns the same rows through our index as through mosdepth's, going through the tabix binary because rust-htslib's tabix reader ends a fetched region with a TabixTruncatedRecord rather than stopping, and does so at different points for the two files. It skips where tabix is absent. CITATIONS.md for a dna run now cites samtools v1.24, the version its fixtures were generated with, instead of the v1.22.1 the rna pipeline was validated against. Each pipeline cites the version it was actually compared with rather than both claiming the newer one. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Matches Picard 3.4.0 byte for byte on the project fixture, metrics row and all 170 histogram lines. Every rule was measured against Picard's own output rather than recalled. The inclusion filter is paired, not secondary, supplementary, duplicate or unmapped, mate mapped, and a 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; minimum and maximum are over the untrimmed set. WIDTH_OF_XX grows a window symmetrically around the median until it covers the percentile, reporting 2i+1; all eleven widths match. The fixture does not exercise trimming, since nothing on it lies beyond ten MADs of the median, so that path has its own unit test. Fixtures are generated with the JVM locale pinned to English: a French default writes "3,531312" where an English one writes "3.531312", which would make them depend on the machine that produced them. Picard's four-line preamble is stripped, holding only a command line and a timestamp. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Matches Picard 3.4.0 on every column and every one of the 251 histogram lines, with two exceptions noted below. The exclusion model was derived by reproducing Picard's own numbers until each fraction matched, not recalled. Unmapped, secondary and supplementary records never enter the calculation. Every other record's reference-consuming bases form the denominator of all PCT_EXC_* columns, 670989 on the fixture. Exclusions then apply in order: 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. SD_COVERAGE is the sample standard deviation over every base of the territory, uncovered ones included. This needs its own depth accumulator rather than a correction applied to the mosdepth one, because the two tools do not agree on which reads or which bases count. That was the design's reason for keeping the accumulators separate and it holds up. HET_SNP_SENSITIVITY and HET_SNP_Q come from Picard's TheoreticalSensitivity, a Monte Carlo simulation whose draws would have to be reproduced bit for bit. 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. CollectWgsMetrics needs --reference to count the reference's non-N bases; without one it is skipped with a warning rather than reported against a wrong genome territory. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Matches Picard 3.4.0 byte for byte: the 101-row detail table and the summary, AT_DROPOUT, GC_DROPOUT and the GC_NC columns included. An earlier attempt to infer the rules from Picard's output failed, so these come from its source, GcBiasUtils and GcBiasMetricsCollector. Three of them would not have been guessed: 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 it. Only unmapped reads and reads with an empty sequence are skipped. Secondary and supplementary alignments count, which is the difference between 5640 and 5642 read starts here. Unmapped reads still count towards TOTAL_CLUSTERS even though they reach nothing else. Two further details came out of the fixture: 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. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
--targets switches the run into targeted mode and produces hs_metrics.txt; --baits defaults to the same intervals. BED input is parsed and merged, since overlapping targets would otherwise inflate the territory and double-count on-target bases. All 58 computable columns match Picard 3.4.0 exactly, including HS_LIBRARY_SIZE, which solves the Lander-Waterman equation by bisection the way Picard's own estimator does. The reason an earlier attempt missed by 0.8 percent is that HsMetrics does not filter the way CollectWgsMetrics does. It clips overlapping mates first, at the read level, and only then applies the base quality floor; WgsMetrics does the opposite. That is why the two report different PCT_EXC_BASEQ and PCT_EXC_OVERLAP on the same file. Two further details came from htsjdk: only the left-most mate is clipped, losing everything from its mate's start onwards, and htsjdk's MATCH_OR_MISMATCH is the M operator alone, so = and X lose their whole element rather than a partial one. Unmapped records reach no contig worker but still count towards TOTAL_READS, PF_BASES and the cluster count, so they are fed to both accumulators during the unmapped pass. Seven columns are not computed: HET_SNP_SENSITIVITY, HET_SNP_Q, the six HS_PENALTY 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. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Three pages under docs/dna: an overview of what the pipeline runs and writes, a mosdepth page, and a Picard page covering all four collectors. They document the things that actually catch people out rather than restating the flags: that --mapq defaults to 0 here and to 30 for rna, that mate-overlap correction is often a factor of two rather than a rounding detail, that the distribution files always emit depths 0 to 300 and never the maximum above that range, that reverse-strand reads are GC-binned by their far end, and that CollectHsMetrics and CollectWgsMetrics filter in opposite orders so their coverage figures are not comparable to each other. Each page also states plainly which columns are not reproduced and why. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Writes genome_results.txt, ten of the raw data tables and an HTML summary. genome_results.txt matches Qualimap 2.3 on every line but four, and the clipping profile, nucleotide content and mapping quality histogram match byte for byte. Several of Qualimap's rules are surprising and none were guessable from its output alone. The reference is split into ceil(len / ceil(len / 400)) windows, so 397 rather than 400. Coverage counts every primary mapped record with no filtering at all, counts deletions, and does not correct mate overlaps, which is why it reports 16.77 where mosdepth reports 6.20. The global mean mapping quality is the mean of the per-window means with empty windows contributing zero, hence 2.4178 rather than about 60, while the per-position histogram truncates that mean instead of rounding it. Mismatches are NM less inserted bases only. Base composition is counted in reference orientation while the clipped span selecting which positions count is taken in sequencing orientation; mixing the two is what Qualimap does and matching it is the only way the composition agrees. Four residuals are documented rather than papered over. The mean mapping quality and the coverage standard deviation differ in the fourth decimal because Qualimap accumulates them per window. About five reference positions of 40001 sit one deeper here, which carries into the coverage histogram and the fractions derived from it. The homopolymer indel classification differs outright: Qualimap reads a reference context this does not reconstruct, and reports two polyC indels that no read-derived rule produces, since the deleted bases are not in the read. Qualimap's GC content distribution and duplication rate histogram are not written. The first is computed over a 679-read subsample whose selection rule is undocumented, the second uses a definition that is not a read-start count. Emitting tables under those names with different numbers would be worse than leaving them out. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Closes seqeralabs#112. nf-core/rnaseq currently runs bedtools genomecov and then UCSC bedGraphToBigWig per strand, three file round-trips for data the alignment pass already has. This computes the coverage as another accumulator in that pass and writes the bigWig directly. Per-base coverage matches bedtools genomecov -bg -split exactly: all 738 intervals on the project fixture. The written track is then read back and checked against the same reference base by base. Its semantics deliberately differ from every other depth engine in the crate, which is why it has its own accumulator rather than reusing one. bedtools filters nothing: duplicates, secondary alignments and both mates of an overlapping pair all contribute. mosdepth excludes the first two and corrects the third; CollectWgsMetrics goes further still. And -split means an N in the CIGAR breaks a read into separate intervals rather than covering the intron, as does a deletion. Two details of the format are worth knowing. A gap between intervals is not zero coverage in a bigWig, it is undefined, and reads back as NaN; that matches what bedGraphToBigWig produces from a bedGraph that omits its zero-depth spans. And the format cannot express a track with no data at all, so a strand with no reads writes no file and says so rather than failing the run, which the fixture exercises: it has no reverse-strand reads. Positions are contig-local, so the accumulators are kept keyed by contig rather than summed across the per-chromosome workers. The growable variant exists because a worker does not know its contig length when the accumulators are built; that it is a field rather than something inferred from the array matters, since after its first read a growable accumulator is otherwise indistinguishable from a fixed one and would silently stop growing. bigtools sits behind a `bigwig` cargo feature, on by default, 19 transitive packages. Both configurations build, lint and test clean, and the MSRV holds. The tracks themselves are off by default: a bigWig is large next to the other outputs and not every run wants one. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Closes #112. Stacked on #156.
nf-core/rnaseq currently runs
bedtools genomecovand then UCSCbedGraphToBigWigper strand: three file round-trips for data the alignment pass already has. This computes the coverage as another accumulator in that pass and writes the bigWig directly.Parity
Per-base coverage matches
bedtools genomecov -bg -splitexactly, all 738 intervals on the project fixture. The written track is then read back withbigtoolsand checked against the same reference base by base, so the container is validated too, not just the numbers going into it.Why it needs its own accumulator
Its semantics differ from every other depth engine in the crate, and deliberately so:
bedtools genomecovCollectWgsMetricsDuplicates and secondary alignments contribute here. Reusing either of the existing engines would silently produce a different track from the one pipelines currently get.
-splitalso means anNin the CIGAR breaks a read into separate intervals rather than covering the intron, as does a deletion.Two things about the format
A gap between intervals is not zero coverage in a bigWig, it is undefined, and reads back as
NaN. That matches whatbedGraphToBigWigproduces from a bedGraph that omits its zero-depth spans, which is whatbedtools genomecov -bgemits. The parity test asserts on covered bases and does not pretend the gaps are zeroes.The format cannot express a track with no data at all, and
bigtoolsrejects it outright. A strand with no reads therefore writes no file and warns, rather than failing the run. The fixture exercises this: it has no reverse-strand reads, so the reverse track is correctly absent.Implementation notes
Positions are contig-local, so accumulators are kept keyed by contig rather than summed across the per-chromosome workers; summing them would merge unrelated coordinates. The growable variant exists because a worker does not know its contig length when the accumulators are built, and growable is a field rather than something inferred from the array: after its first read a growable accumulator is otherwise indistinguishable from a fixed one of that size and would silently stop growing. That was a real bug caught by the merge test.
Configuration
Off by default, since a bigWig is large next to the other outputs:
bigtoolssits behind abigwigcargo feature, on by default, 19 transitive packages. Both configurations build, lint and test clean, and the 1.87 MSRV holds.Tests
369 green with the feature, 365 without.
🤖 Generated with Claude Code