feat(rna): RSeQC geneBody_coverage and read_GC - #161
Open
BenjaminDEMAILLE wants to merge 22 commits into
Open
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 the two RSeQC tools RustQC was missing, bringing the count from eight to ten. Both are validated against RSeQC 5.0.4's own output. read_GC reproduces the GC distribution exactly: all 724 distinct values on the fixture. Its semantics came from qcmodule.SAM.readGC rather than from its output, and two of them are easy to get wrong. The denominator is the full sequence length, so an N pads it without ever counting as GC and pulls the percentage down; dividing by the called bases is the more defensible measurement and the wrong one to report under this name. And the percentage is formatted with %4.2f, which rounds ties to the even digit: a read with 9 of 32 bases G or C is exactly 28.125 percent and belongs in the 28.12 bin, which on this fixture is fourteen reads. geneBody_coverage reproduces all 100 percentile points exactly. Its algorithm came from the source too, and three of its behaviours are defects being reproduced deliberately rather than design decisions: The interpolation rounds with Python's round, which is half-to-even. Rounding away from zero moves a position by one base wherever the interpolation lands halfway, which changes the 23rd percentile by a whole read here. The concatenated mRNA bases are never sorted. A transcript with overlapping exons therefore feeds a list that dips backwards into an interpolation that assumes sorted input, and the resulting positions repeat. Coverage is kept in a dictionary keyed by position, so those repeats collapse and the transcript's curve is aggregated from fewer than 100 values, quietly compressing it towards the 5' end. One gene of eight in the fixture has 100 points across 88 distinct positions. Each of these is spelled out where it is implemented, so anyone deciding to diverge from RSeQC does so knowingly rather than by accident. The read_GC fixture is the DNA alignment rather than the RNA one: the RNA test BAM carries synthetic all-A sequences, so every read there is 0 percent GC and the histogram collapses to a single row that would prove nothing. 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 #127. Stacked on #156, the tip of the DNA stack.
Adds the two RSeQC tools RustQC was missing, taking the count from eight to ten. Both are validated against RSeQC 5.0.4's own output on committed fixtures.
read_GC.pygeneBody_coverage.pySemantics taken from the source, not the output
An earlier pass in this series established that inferring a tool's rules from its output is slower and less reliable than reading it. Both tools bear that out.
read_GC's denominator is the full sequence length. AnNpads it without ever counting as GC, pulling the percentage down. Dividing by the called bases is the more defensible measurement and the wrong one to report under this name.Its percentages round half to even, because RSeQC formats with
%4.2f. A read with 9 of 32 bases G or C is exactly 28.125 percent and belongs in the28.12bin, not28.13. That is fourteen reads on this fixture.Three RSeQC defects reproduced on purpose
geneBody_coveragehas three behaviours that are bugs rather than design decisions. They are reproduced because the point of the module is to agree with the tool people already have, and each is spelled out where it is implemented so that anyone choosing to diverge does so knowingly.The interpolation rounds with Python's
round, which is half-to-even. Rounding away from zero moves an interpolated position by one base wherever the interpolation lands exactly halfway, and that alone changes the 23rd percentile by a whole read here.The concatenated mRNA bases are never sorted. A transcript whose exons overlap therefore feeds a list that dips backwards into an interpolation that assumes sorted input, and the positions it returns repeat.
Coverage is kept in a dictionary keyed by position, so those repeats collapse. The transcript's curve is then aggregated from fewer than 100 values, added at indices 0 upwards, quietly compressing it towards the 5' end. One gene of eight in the fixture has 100 percentile points across 88 distinct positions, and it was the single differing value that led to all three findings.
A fourth quirk is documented but harmless: because the interpolation runs between neighbouring mRNA bases, a percentile point can land in an intron, and coverage is then read at a base the transcript does not contain.
Fixtures
read_GCis validated against the DNA alignment rather than the RNA one. The RNA test BAM carries synthetic all-A sequences, so every read there is 0 percent GC and the histogram collapses to a single row that would prove nothing.geneBody_coverageuses the RNA BAM with a BED12 gene model derived from the existing test GTF, 4 kB.Tests
376 green, with 13 unit tests covering the percentile interpolation, the strand reversal, the repeat collapse and the rounding, plus two parity tests.
🤖 Generated with Claude Code