feat(dna): reimplement Qualimap bamqc

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>

Author: BenjaminDEMAILLE

AGENTS.md

@@ -82,6 +82,8 @@ src/
     hs_metrics.rs     — Picard CollectHsMetrics reimplementation (targeted mode)
     insert_size.rs    — Picard CollectInsertSizeMetrics reimplementation
     intervals.rs      — BED interval parsing and merging for targeted mode
+    qualimap.rs       — Qualimap bamqc accumulation (windows, coverage, composition)
+    qualimap_output.rs — genome_results.txt, the raw data tables and the HTML report
     wgs_metrics.rs    — Picard CollectWgsMetrics reimplementation
     mosdepth/
       mod.rs          — Per-contig summarisation feeding the mosdepth outputs

CHANGELOG.md

@@ -10,8 +10,9 @@
   alignment with one worker per contig, plus Picard-compatible
   CollectWgsMetrics, CollectInsertSizeMetrics and CollectGcBiasMetrics.
   Passing `--targets` switches on targeted mode and Picard-compatible
-  CollectHsMetrics. Validated for exact parity against mosdepth 0.3.14,
-  samtools 1.24 and Picard 3.4.0.
+  CollectHsMetrics. Qualimap-compatible `bamqc` output rounds it out, with
+  `genome_results.txt`, the raw data tables and an HTML summary. Validated
+  against mosdepth 0.3.14, samtools 1.24, Picard 3.4.0 and Qualimap 2.3.
 
 ### Changed
 

docs/astro.config.mjs

@@ -78,6 +78,7 @@ export default defineConfig({
             { label: "Overview", slug: "dna/overview" },
             { label: "mosdepth", slug: "dna/mosdepth" },
             { label: "Picard metrics", slug: "dna/picard" },
+            { label: "Qualimap bamqc", slug: "dna/qualimap" },
           ],
         },
         {

docs/src/content/docs/dna/qualimap.mdx

@@ -0,0 +1,69 @@
+---
+title: Qualimap bamqc
+description: The bamqc outputs RustQC produces, how Qualimap's coverage differs from every other tool here, and which figures are not reproduced.
+---
+
+import { Aside } from "@astrojs/starlight/components";
+
+<Aside type="note" title="About Qualimap">
+Qualimap examines alignment data to detect biases in sequencing and mapping,
+and reports them as tables and plots.
+
+[Documentation](http://qualimap.conesalab.org/) | Okonechnikov K, Conesa A, García-Alcalde F. *Bioinformatics*. 2016;32(2):292-294.
+</Aside>
+
+RustQC writes `qualimap/genome_results.txt`, the
+`raw_data_qualimapReport/` tables and an HTML summary.
+
+<Aside type="caution" title="Qualimap's coverage is not the others' coverage">
+Qualimap counts every primary mapped record with no duplicate, mapping quality
+or base quality filtering, counts deletions as covered, and does **not**
+correct mate overlaps. On the project's test data it reports a mean coverage of
+16.77 where mosdepth reports 6.20 and `CollectWgsMetrics` reports 3.53. All
+three are right about different questions. Do not compare them to each other.
+</Aside>
+
+## Two figures that surprise people
+
+**Mean mapping quality reads about 2.4, not about 60.** It is the mean of the
+per-window means, and a window with no reads contributes zero. On a targeted or
+low-coverage library, most windows are empty, so the figure is closer to the
+fraction of the genome covered than to the quality of the alignments. The
+per-position histogram, which only counts covered positions, is the one to read
+for that.
+
+**Base composition is reported in reference orientation.** Reverse-strand reads
+are reverse-complemented before counting, so the A and T columns are not the
+counts of A and T in the sequencer's output.
+
+## What matches Qualimap and what does not
+
+`genome_results.txt` matches on every line but four, and three of the raw
+tables match byte for byte. The residuals, each with its cause:
+
+| Figure | Difference |
+| --- | --- |
+| `mean mapping quality` | fourth decimal; 393 of 397 windows match exactly |
+| `std coverageData` | fourth decimal, same cause |
+| `homopolymer indels` | differs outright, see below |
+| coverage histogram and what derives from it | about five reference positions of 40001 sit one deeper |
+| `genome_fraction_coverage` | last two digits of the double, Qualimap accumulates per window |
+| `insert_size_histogram` | one extra row: Qualimap trims the largest insert from the plotted table while still counting it in the statistics |
+
+Qualimap classifies an indel as a homopolymer indel against a reference context
+RustQC does not reconstruct. It reports two polyC indels on the test data, and
+no rule derived from the read alone produces them, because the deleted bases
+are not in the read. RustQC uses a run-of-four rule instead, so this one figure
+will differ.
+
+Qualimap's GC content distribution and duplication rate histogram are not
+written at all. The first is computed over a 679-read subsample whose selection
+rule is not documented; the second uses a definition that does not match a
+read-start-position count. Emitting tables under those names with different
+numbers would be worse than leaving them out.
+
+## The HTML report
+
+RustQC writes its own summary page rather than a copy of Qualimap's, which
+ships a bundle of images, CSS and JavaScript. It carries the same numbers as
+`genome_results.txt`; the raw tables remain the machine-readable source.

src/citations.rs

@@ -148,6 +148,9 @@ pub fn write_dna_citations(
     if config.samtools.enabled {
         write_citation(&mut w, &SAMTOOLS_DNA)?;
     }
+    if config.qualimap.enabled {
+        write_citation(&mut w, &QUALIMAP)?;
+    }
     if config.preseq.enabled {
         write_citation(&mut w, &PRESEQ)?;
     }

src/config.rs

@@ -988,6 +988,10 @@ pub struct DnaConfig {
     #[serde(default)]
     pub hs_metrics: HsMetricsConfig,
 
+    /// Qualimap bamqc configuration.
+    #[serde(default)]
+    pub qualimap: BamqcConfig,
+
     /// preseq lc_extrap library complexity extrapolation configuration.
     ///
     /// Reuses the same type as the `rna` pipeline; the implementation is shared.
@@ -1032,6 +1036,38 @@ impl Default for WgsMetricsConfig {
     }
 }
 
+/// Configuration for the Qualimap-compatible bamqc report.
+///
+/// Named apart from the `rna` pipeline's [`QualimapConfig`], which configures
+/// a different Qualimap analysis entirely: gene body coverage rather than
+/// bamqc.
+///
+/// Example:
+/// ```yaml
+/// qualimap:
+///   enabled: true
+///   num_windows: 400
+/// ```
+#[derive(Debug, Deserialize)]
+#[serde(default)]
+pub struct BamqcConfig {
+    /// Whether to produce the bamqc outputs. Defaults to true.
+    pub enabled: bool,
+    /// Target number of windows the reference is split into. The realised
+    /// count is usually a little lower, because the window width is rounded up
+    /// first.
+    pub num_windows: usize,
+}
+
+impl Default for BamqcConfig {
+    fn default() -> Self {
+        Self {
+            enabled: true,
+            num_windows: 400,
+        }
+    }
+}
+
 /// Configuration for the Picard-compatible GC bias metrics.
 ///
 /// Requires a reference FASTA: the analysis bins reference windows by GC.

src/dna/mod.rs

@@ -10,4 +10,6 @@ pub mod hs_metrics;
 pub mod insert_size;
 pub mod intervals;
 pub mod mosdepth;
+pub mod qualimap;
+pub mod qualimap_output;
 pub mod wgs_metrics;