Tracking: protein (proteomics) QC subcommand PR series (#158, #159) and raw read QC (#160)

[BenjaminDEMAILLE]

Adds a rustqc protein subcommand for proteomics, plus a rustqc reads subcommand for raw FASTQ that came out of the same work.

Scope correction: an earlier revision of this series also put coding-region QC under protein. That was wrong. CollectRnaSeqMetrics is an RNA-seq tool and protein is for proteomics: mass spectrometry, and the protein FASTA that spectra are searched against. #162 is closed and the work moved to #164, in the rna pipeline, where it also costs no new command-line surface because rustqc rna already takes an alignment and an annotation.

Nothing in the tracker asked for these, unlike #128 which the DNA stack answered. They extend RustQC past its current "genomics QC for sequencing data" scope, so the scope question is a real one and belongs in review rather than after it.

Review order

protein sequence and protein spectra are a chain; reads branches off protein spectra.

PR	What it does	Input	Tests
#158	protein sequence: FASTA statistics, composition and defects	protein FASTA	386
#159	protein spectra: mass spectrometry run QC	mzML	398
#160	rustqc reads: raw FASTQ QC before alignment	FASTQ	416

All four are stacked on #156, the tip of the DNA series tracked in #157, because they reuse src/common from #152.

Why two modes rather than two subcommands

sequence and spectra take entirely different inputs and share no analysis code. They belong together because a protein FASTA is the search database that mass spectra are matched against; they are two halves of one workflow. A flat command would need a validation matrix explaining which flag applies when; an explicit mode makes the incompatibility structural. rustqc protein <file> is rejected rather than guessed at.

reads is separate rather than a fourth mode because the question is different in kind: the protein modes are assay-specific, raw read QC applies whatever the library was for.

Parity

Output	Reference	Result
sequence_stats.tsv	seqkit 2.13.0	all 19 columns, both fixtures
spectra_report.txt	pyteomics 5.0.1	every counter both readers can report
read_stats.tsv	seqkit 2.13.0	all 19 columns
fastqc_data.txt	FastQC 0.12.1	4 modules of 5

The mass spectrometry check is worth singling out: mzdata in Rust and pyteomics in Python agree on all 48 spectra and 305213 peaks of the fixture. Two independent readers agreeing is a stronger statement than either agreeing with itself.

The dependency question

mzdata 0.66.5 is the first dependency RustQC would carry for a single output. It sits behind a proteomics cargo feature, on by default, and --no-default-features drops the mode from the help entirely rather than offering it and failing. Both configurations build, lint and test clean. It was checked against the 1.87 MSRV before the design was committed to. Whether that trade is worth making is a maintainer's call, not an implementer's.

Upstream rules that would have passed a less careful test

Each is commented where it is implemented.

seqkit's quartiles are Tukey's halves, not interpolated: Q1 is the median of the lower half. On the yeast fixture that is 157 where linear interpolation gives 165.

Those halves round half to even. 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, which is exactly the shape of bug that passes on one file.

seqkit's N50_num counts distinct lengths, not sequences. Three reads of 10, 10 and 3 give 1, not 2. Every length in the protein fixtures happens to be distinct, which hid this completely until the same code met FASTQ, where thousands of reads share a length and the wrong definition reported 3945 against seqkit's 1. #158 carries the fix.

seqkit's AvgQual averages error probabilities, not Phred scores. On the FASTQ fixture the arithmetic mean is 34.00 while the reported figure is 25.60. Publishing the arithmetic mean would flatter every run.

FastQC's binned rows average their positions rather than pooling their bases, quantiles included, which is why a bin's tenth percentile can read 35.2 when every position's own is a whole number. The definitions agree until reads start running out.

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

Known gaps

FastQC's per sequence GC content is written but not asserted. FastQC spreads each read's contribution across neighbouring bins so a coarse discrete distribution plots smoothly, making its counts fractional and placing reads in bins none occupies. RustQC writes the plain rounded distribution, a different and defensible figure, so asserting equality would be asserting the wrong thing.

FastQC module verdicts are written as pass uniformly. Its pass, warn and fail thresholds are judgement rather than data. This is the one place the file is not a drop-in replacement.

Fixtures

Roughly 3.3 MB in total, all small, real and public: two protein FASTA files and an mzML from mzdata's test data, and a FASTQ from nf-core. tests/create_protein_test_data.sh regenerates inputs and reference outputs, pinning every tool version and refusing to run against others.

🤖 Generated with Claude Code