Skip to content

feat(reads): raw FASTQ quality control - #160

Open
BenjaminDEMAILLE wants to merge 25 commits into
seqeralabs:mainfrom
BenjaminDEMAILLE:feat/reads-fastq
Open

feat(reads): raw FASTQ quality control#160
BenjaminDEMAILLE wants to merge 25 commits into
seqeralabs:mainfrom
BenjaminDEMAILLE:feat/reads-fastq

Conversation

@BenjaminDEMAILLE

@BenjaminDEMAILLE BenjaminDEMAILLE commented Aug 29, 2026

Copy link
Copy Markdown

Part of #163.

Stacked on #159.

Adds rustqc reads: quality control on raw FASTQ, before alignment. It sits apart from rna, dna and protein because the question is different. Those ask what the alignment says; this asks whether the reads are usable at all, whatever the library was for.

Writes a seqkit stats-compatible table and a FastQC-compatible fastqc_data.txt, so existing parsers keep working. Records stream through the accumulator and are never all held at once.

Parity

Full parity with seqkit 2.13.0 on all nineteen columns, and with FastQC 0.12.1 on four modules: per base sequence quality, per sequence quality scores, per base sequence content, and per base N content.

Four conventions that would have passed a less careful test

N50_num counts distinct lengths, not sequences. Three reads of 10, 10 and 3 give 1, not 2, because the two tens are one length. A file whose lengths are all distinct hides this completely, which is exactly how it slipped past the protein fixture in #158 where every length happened to be unique. That implementation is corrected here too, and the discriminating case now has its own test.

AvgQual averages error probabilities, not Phred scores. On this fixture the arithmetic mean is 34.00 while seqkit reports 25.60, because averaging in probability space is dominated by the worst bases. Reporting the arithmetic mean would flatter every run.

FastQC's binned rows average their positions rather than pooling their bases, for the quantiles as well as the mean. That is why a bin's tenth percentile can read 35.2 when every individual position's is a whole number. The two definitions agree until reads start running out, so the difference first appears at the 40 to 44 bin, in the seventh decimal for the mean and outright for the percentiles.

FastQC excludes N from the per-base composition denominator, so the four percentages sum to 100 even where the instrument called nothing. Including them shifts every percentage by the N rate, which on this fixture moves the fourth decimal.

One module written but not asserted

Per sequence GC content. FastQC does not report the plain distribution of per-read GC percentage: it spreads each read's contribution across neighbouring bins so a coarse discrete distribution plots smoothly, which makes its counts fractional and puts reads in bins none actually occupies. RustQC reports the plain rounded distribution, which is a different and defensible figure, so asserting equality would be asserting the wrong thing. It is called out in the code, the test and the report.

Module verdicts are written as pass uniformly rather than reproducing FastQC's pass/warn/fail thresholds, which are judgement rather than data. That is the one place the file is not a drop-in replacement.

Fixture

test.umi_1.fastq.gz from nf-core/test-datasets, 302 kB, 7890 reads of 41 to 151 bases with a few N calls, which is enough to exercise the ragged-length and N-content paths that a uniform file would not.

Tests

416 green, up from 397.

🤖 Generated with Claude Code

BenjaminDEMAILLE and others added 25 commits August 28, 2026 18:26
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>
rustqc protein sequence reads protein FASTA and reports length statistics,
amino acid composition and defects. The statistics table reproduces
seqkit 2.13.0's `stats -a -T` byte for byte on both fixtures.

Two of seqkit's conventions had to be recovered from its behaviour, and
neither is what a statistics library gives by default. Its quartiles are
Tukey's halves, so Q1 is the median of the lower half rather than an
interpolated value: on the yeast fixture that is 157 where linear
interpolation gives 165. And those halves round half-to-even, so a median of
235.5 is reported as 236 while 376.5 and 516.5 are reported as 376 and 516.
Ordinary rounding gets the first right and the other two wrong.

Beyond seqkit, the report carries per-residue composition and the defects that
make a proteome unusable downstream: stop codons anywhere but the final
position, residues that are not amino acids, byte-identical duplicate
sequences and reused identifiers. Ambiguity codes are counted rather than
flagged, and selenocysteine and pyrrolysine are treated as standard.

The subcommand takes an explicit mode rather than inferring one from which
flags were given, because the three planned modes take entirely different
inputs. Only `sequence` exists so far.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
seqkit walks the distinct lengths from longest down and reports how many it
consumed. Three sequences of lengths 10, 10 and 3 give 1, not 2, because the
two tens are one length.

Both project fixtures happen to have every length distinct, which makes the
two definitions agree and hid this entirely. The discriminating case now has
its own test.

Found while extending the same statistics to FASTQ, where thousands of reads
share a length and the wrong definition reported 3945 against seqkit's 1.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
rustqc protein spectra reads mzML and reports per-level spectrum and peak
counts, total ion current, the retention time range, the fragmentation ratio
and the precursor charge distribution.

mzdata does the reading, behind a `proteomics` cargo feature that is on by
default. Building with --no-default-features drops the mode from the help
entirely rather than offering it and then failing, and drops mzdata's
transitive packages with it. Both configurations are tested.

Validated against pyteomics 5.0.1 on mzdata's own test file: two independent
readers agree on all 48 spectra, 305213 peaks, the per-level minima and
maxima, the retention time bounds and the precursor m/z range. Total ion
current is compared with a relative tolerance rather than exactly, because
the intensities are 32-bit in the file and the two readers accumulate them at
different precision.

Two reporting decisions worth noting. Total ion current is summed from the
peaks actually present rather than read from the spectrum header, so a
profile spectrum that has since been centroided reports what the file now
holds. And a run with no survey scans reports its fragmentation ratio as NA
rather than 0, since 0 would read as "no fragmentation" when the truth is
"nothing to divide by".

The fixture annotates no precursor charge states, which the run warns about
and the report shows as `without_charge`, rather than quietly reporting an
empty charge distribution.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
rustqc reads streams FASTQ and writes a seqkit-compatible statistics table
and a FastQC-compatible fastqc_data.txt. It sits apart from rna, dna and
protein because the question is different: those ask what the alignment says,
this asks whether the reads are usable at all, whatever the library was for.

Full parity with seqkit 2.13.0 on all nineteen columns, and with FastQC
0.12.1 on four modules: per base sequence quality, per sequence quality
scores, per base sequence content and per base N content.

Four upstream conventions had to be recovered, and every one of them would
have passed a less careful test:

seqkit's N50_num counts distinct lengths, not sequences. Three reads of 10,
10 and 3 give 1, not 2, because the two tens are one length. A file whose
lengths are all distinct hides this completely, which is exactly how it slipped
past the protein fixture; that implementation is corrected here too.

seqkit's AvgQual averages error probabilities, not Phred scores. On this
fixture the arithmetic mean is 34.00 while the reported figure is 25.60,
because averaging in probability space is dominated by the worst bases.
Reporting the arithmetic mean would flatter every run.

FastQC's binned rows average their positions rather than pooling their bases,
for the quantiles as well as the mean. That is why a bin's tenth percentile
can read 35.2 when every individual position's is a whole number. The two
definitions agree until reads start running out, so the difference first
appears at the 40 to 44 bin.

FastQC excludes N from the per-base composition denominator, so the four
percentages sum to 100 even where the instrument called nothing.

One module is written but not asserted against FastQC: per sequence GC
content. FastQC spreads each read's contribution across neighbouring bins so
a coarse discrete distribution plots smoothly, which makes its counts
fractional. RustQC reports the plain rounded distribution, a different and
defensible figure, so asserting equality would be asserting the wrong thing.
It is called out in the code and the report.

Module verdicts are written as `pass` uniformly rather than reproducing
FastQC's thresholds, which are judgement rather than data.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant