feat(protein): add the protein subcommand and its sequence mode - #158
Open
BenjaminDEMAILLE wants to merge 23 commits into
Open
feat(protein): add the protein subcommand and its sequence mode#158BenjaminDEMAILLE wants to merge 23 commits into
BenjaminDEMAILLE wants to merge 23 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>
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>
This was referenced Aug 28, 2026
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>
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.
Part of #163.
First of a protein QC series. Stacked on #156, the tip of the DNA stack, because it reuses
src/commonfrom #152.Adds
rustqc protein sequence: protein FASTA quality control.Parity against seqkit 2.13.0
The statistics table reproduces
seqkit stats -a -Tbyte for byte on both fixtures, two small real protein FASTA files from nf-core/test-datasets totalling 16 kB.Two of seqkit's conventions had to be recovered from its behaviour, and neither is what a statistics library gives you by default:
Quartiles are Tukey's halves.
Q1is the median of the lower half of the sorted lengths, not an interpolated value. On the yeast fixture that is 157, where linear interpolation gives 165 and inverse-ECDF agrees onQ1but then disagrees onQ2.Those halves round half-to-even. A median of 235.5 is reported as 236, but 376.5 and 516.5 are reported as 376 and 516. Ordinary rounding gets the first right and the other two wrong, which is exactly the kind of near-miss that looks like a passing test on one file.
Beyond seqkit
seqkit answers "how long are these sequences". The report answers "is this proteome usable":
Ambiguity codes (
B,Z,J,X) are counted rather than flagged, since they are valid. Selenocysteine and pyrrolysine are treated as standard residues, because translated proteomes carry them.Command shape
The subcommand takes an explicit mode:
rather than inferring one from which flags were given. Two further modes are planned,
coding(BAM plus GTF) andspectra(mzML), and they take entirely different inputs; a flat command would need a validation matrix explaining which flag applies when. Onlysequenceexists so far, andrustqc protein <file>is rejected rather than silently guessing.Shared options keep the names, short flags and
RUSTQC_*environment variables they have inrnaanddna.Tests
385 green, up from 353. Six are protein parity tests; the rest are unit tests over the parsing, the statistics conventions and each defect category.
🤖 Generated with Claude Code