feat(dna): Qualimap bamqc, plus the DNA documentation pages

AGENTS.md

@@ -62,10 +62,35 @@ src/
   config.rs         — YAML configuration loading (serde), nested tool configs
   io.rs             — Shared I/O utilities (gzip-transparent file reading)
   gtf.rs            — GTF annotation file parser (with configurable attribute extraction)
+  common/
+    mod.rs            — Re-exports the shared modules
+    bam_flags.rs      — BAM flag constants and aux-tag helpers
+    bam_stat.rs       — bam_stat.py reimplementation, result types
+    bam_stat_accum.rs — Read-level counter accumulator feeding bam_stat and samtools
+    cpp_rng.rs        — C++ RNG FFI shim for preseq bootstrap reproducibility
+    preseq.rs         — preseq lc_extrap library complexity extrapolation
+    samtools/
+      mod.rs          — Re-exports the samtools writers
+      stats.rs        — samtools stats full output (SN + all histogram sections)
+      flagstat.rs     — samtools flagstat-compatible output
+      idxstats.rs     — samtools idxstats-compatible output
+  dna/
+    mod.rs            — Re-exports the DNA submodules
+    depth.rs          — Per-contig depth accumulator (delta array, CIGAR walk,
+                        mate-overlap correction, prefix sum)
+    gc_bias.rs        — Picard CollectGcBiasMetrics reimplementation
+    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
+      output.rs       — The six mosdepth-compatible writers (bgzf for the BED outputs)
   rna/
-    mod.rs          — Re-exports all submodules (dupradar, featurecounts, rseqc, bam_flags, cpp_rng, preseq, qualimap)
-    bam_flags.rs    — BAM flag constants
-    cpp_rng.rs      — C++ RNG FFI shim for preseq bootstrap reproducibility
+    mod.rs          — Re-exports the RNA submodules (dupradar, featurecounts, rseqc, qualimap)
+                      and re-exports the shared ones from `common` for compatibility
     dupradar/
       mod.rs        — Re-exports counting, dupmatrix, fitting, plots
       counting.rs   — BAM read counting engine (largest module)
@@ -75,7 +100,6 @@ src/
     featurecounts/
       mod.rs        — Re-exports output
       output.rs     — featureCounts-format output & biotype counting
-    preseq.rs       — preseq lc_extrap library complexity extrapolation
     qualimap/
       mod.rs        — Re-exports all Qualimap modules
       accumulator.rs — Gene body coverage accumulation logic
@@ -88,17 +112,13 @@ src/
       mod.rs                — Re-exports all RSeQC modules + common helpers
       accumulators.rs       — Shared RSeQC accumulator infrastructure (read dispatch)
       common.rs             — Shared junction/intron extraction, from_genes builders
-      bam_stat.rs           — bam_stat.py reimplementation
-      flagstat.rs           — samtools flagstat-compatible output
-      idxstats.rs           — samtools idxstats-compatible output
       infer_experiment.rs   — infer_experiment.py reimplementation
       inner_distance.rs     — inner_distance.py reimplementation
       junction_annotation.rs — junction_annotation.py reimplementation
       junction_saturation.rs — junction_saturation.py reimplementation
       plots.rs              — RSeQC plot generation (duplication, junctions, etc.)
       read_distribution.rs  — read_distribution.py reimplementation
       read_duplication.rs   — read_duplication.py reimplementation
-      stats.rs              — samtools stats full output (SN + all histogram sections)
       tin.rs                — TIN (Transcript Integrity Number) analysis
 tests/
   integration_test.rs  — 12 integration tests vs R dupRadar reference output
@@ -107,13 +127,23 @@ tests/
   create_test_data.R   — R script to regenerate test data + references
 ```
 
-Nested module structure — top-level modules (`cli`, `config`, `io`, `gtf`, `rna`) declared
-in `main.rs`, no `lib.rs`. The `rna` module contains sub-modules for each tool group.
-Inter-module access uses `crate::` paths (e.g., `use crate::rna::dupradar::counting::GeneCounts;`).
+Nested module structure. The library crate root is `src/lib.rs`, which declares
+`common`, `config`, `cpu`, `gtf`, `io`, `rna` and `summary`; the binary
+(`src/main.rs`) additionally declares `cli`, `citations` and `ui`.
+Inter-module access uses `crate::` paths (e.g., `use crate::common::bam_stat_accum::BamStatAccum;`).
+Assay-agnostic analyses belong in `common`; put new code under `rna` only if it
+needs a gene annotation or a library strand protocol.
 
-The CLI uses a single subcommand:
+The CLI has two subcommands:
 
 - `rustqc rna <BAM>... --gtf <GTF> [OPTIONS]`
+- `rustqc dna <BAM>... [OPTIONS]`
+
+The `dna` subcommand needs no annotation. It runs depth of coverage
+(mosdepth-compatible), the samtools-compatible outputs and preseq in one pass,
+with one worker per contig. Shared flags keep their `rna` names, short forms
+and `RUSTQC_*` environment variables, with one deliberate exception:
+`-Q/--mapq` defaults to 0 for `dna`, matching mosdepth, rather than 30.
 
 A GTF gene annotation file (`--gtf`) is required. This runs all analyses:
 dupRadar duplicate rate analysis, featureCounts-compatible gene counting,

CHANGELOG.md

@@ -1,5 +1,26 @@
 # RustQC Changelog
 
+## Unreleased
+
+### Features
+
+- New `rustqc dna` subcommand for DNA (WGS) quality control: depth of coverage
+  with mosdepth-compatible outputs, samtools-compatible stats, flagstat and
+  idxstats, and preseq library complexity, all in a single pass over the
+  alignment with one worker per contig, plus Picard-compatible
+  CollectWgsMetrics, CollectInsertSizeMetrics and CollectGcBiasMetrics.
+  Passing `--targets` switches on targeted mode and Picard-compatible
+  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
+
+- Internal: assay-agnostic analyses (BAM flag helpers, read-level statistics,
+  the samtools stats/flagstat/idxstats writers, preseq) moved from `rna` to a
+  new `common` module. The old `rustqc::rna::...` paths still resolve through
+  re-exports, so this is not a breaking change for library users.
+
 ## [Version 0.2.1](https://github.com/seqeralabs/RustQC/releases/tag/v0.2.1) - 2026-04-09
 
 ### Bug fixes

docs/astro.config.mjs

@@ -72,6 +72,15 @@ export default defineConfig({
             { label: "Samtools", slug: "rna/samtools" },
           ],
         },
+        {
+          label: "DNA",
+          items: [
+            { label: "Overview", slug: "dna/overview" },
+            { label: "mosdepth", slug: "dna/mosdepth" },
+            { label: "Picard metrics", slug: "dna/picard" },
+            { label: "Qualimap bamqc", slug: "dna/qualimap" },
+          ],
+        },
         {
           label: "About",
           items: [

docs/src/content/docs/dna/mosdepth.mdx

@@ -0,0 +1,69 @@
+---
+title: mosdepth
+description: Depth of coverage outputs compatible with mosdepth, and the semantics RustQC reproduces.
+---
+
+import { Aside } from "@astrojs/starlight/components";
+
+<Aside type="note" title="About mosdepth">
+mosdepth computes depth of coverage from BAM and CRAM files quickly and
+without a per-base pileup.
+
+[Source code](https://github.com/brentp/mosdepth) | Pedersen BS, Quinlan AR. *Bioinformatics*. 2018;34(5):867-868.
+</Aside>
+
+RustQC reproduces mosdepth's outputs exactly. On the project's test alignment
+every one of the six files is identical to mosdepth 0.3.14's, including the
+1094-line global distribution and the 721-interval per-base BED.
+
+## Files
+
+| File | Written when |
+| --- | --- |
+| `{sample}.mosdepth.summary.txt` | always |
+| `{sample}.mosdepth.global.dist.txt` | always |
+| `{sample}.per-base.bed.gz` (+ `.csi`) | unless `--skip-per-base` |
+| `{sample}.regions.bed.gz` (+ `.csi`) | with `--window-size` |
+| `{sample}.mosdepth.region.dist.txt` | with `--window-size` |
+| `{sample}.thresholds.bed.gz` (+ `.csi`) | with `--window-size` and thresholds |
+
+Compressed outputs are bgzf with a CSI companion index, so `tabix` can seek
+into them just as it can into mosdepth's own.
+
+## What counts towards depth
+
+RustQC applies mosdepth's default filters:
+
+- records carrying any of `UNMAP`, `SECONDARY`, `QCFAIL` or `DUP` are skipped,
+  which is mosdepth's `-F 1796`;
+- records below `--mapq` are skipped, defaulting to 0;
+- `M`, `=` and `X` cover the reference; `D` and `N` advance without covering;
+  `I`, `S`, `H` and `P` do not advance at all;
+- **a base covered by both mates of one pair counts once.**
+
+<Aside type="tip" title="Mate overlap is not a rounding detail">
+That last rule is what mosdepth's `--fast-mode` turns off, and it is often the
+largest single effect on a coverage figure. On the project's test data,
+correcting for it takes total covered bases from 469875 down to 247878. If your
+coverage looks about twice what you expect, this is the first thing to check in
+whatever produced it.
+</Aside>
+
+## Reading the distribution files
+
+`{sample}.mosdepth.global.dist.txt` holds `chrom`, `depth` and `proportion`
+rows in descending depth order, where the proportion is the fraction of that
+contig's bases at depth **at or above** the given value, ending at depth 0 with
+`1.00`.
+
+Which depths get a row is worth knowing, because it is not simply "every depth
+seen":
+
+- depths 0 through 300 always get a row, whether or not any base sits at that
+  exact depth;
+- above 300, only depths that actually occur;
+- the maximum observed depth gets a row when it falls inside that dense range,
+  and none when it does not.
+
+The region distribution follows the same rules but is computed over windows and
+their **rounded mean** depth, not over individual bases.

docs/src/content/docs/dna/overview.mdx

@@ -0,0 +1,99 @@
+---
+title: DNA QC Overview
+description: What the rustqc dna subcommand runs, what it writes, and how it differs from the RNA pipeline.
+---
+
+import { Aside, FileTree } from "@astrojs/starlight/components";
+
+`rustqc dna` runs a DNA (whole-genome or targeted) quality control pipeline in
+a single pass over each alignment file. Unlike [`rustqc rna`](/rna/dupradar/),
+it needs no gene annotation.
+
+```bash
+rustqc dna sample.bam --reference genome.fasta --outdir results/
+```
+
+## What it runs
+
+| Upstream tool | What RustQC produces |
+| --- | --- |
+| [mosdepth](/dna/mosdepth/) | depth of coverage, per base, per window and per region |
+| [Picard `CollectWgsMetrics`](/dna/picard/) | genome-wide coverage metrics with the exclusion breakdown |
+| [Picard `CollectInsertSizeMetrics`](/dna/picard/) | insert size distribution per pair orientation |
+| [Picard `CollectGcBiasMetrics`](/dna/picard/) | coverage bias against reference GC content |
+| [Picard `CollectHsMetrics`](/dna/picard/) | targeted enrichment metrics, when `--targets` is given |
+| [Samtools](/rna/samtools/) | `stats`, `flagstat` and `idxstats` |
+| [Preseq](/rna/preseq/) | library complexity extrapolation |
+
+Every one of them is fed from the same record stream, so the alignment is read
+once no matter how many are enabled.
+
+## Requirements
+
+The input must be **duplicate-marked, not duplicate-removed**. Duplicate rate
+is a headline metric here, and several exclusion fractions are defined against
+it. RustQC refuses input with no duplicate flags at all unless you pass
+`--skip-dup-check`.
+
+A reference FASTA is needed for three things: reading CRAM, sizing
+`GENOME_TERRITORY` for `CollectWgsMetrics`, and binning reference windows for
+`CollectGcBiasMetrics`. Without one, those two analyses are skipped with a
+warning and everything else still runs.
+
+<Aside type="caution" title="One default differs from the RNA pipeline">
+`-Q/--mapq` defaults to **0** for `dna`, not to 30 as it does for `rna`. That
+is mosdepth's own default, and matching it is what makes the depth figures
+agree with mosdepth exactly. Raising it will move every coverage number.
+</Aside>
+
+## Output tree
+
+<FileTree>
+- results/
+  - mosdepth/
+    - sample.mosdepth.summary.txt
+    - sample.mosdepth.global.dist.txt
+    - sample.mosdepth.region.dist.txt
+    - sample.per-base.bed.gz
+    - sample.per-base.bed.gz.csi
+    - sample.regions.bed.gz
+    - sample.thresholds.bed.gz
+  - picard/
+    - wgs_metrics/
+    - insert_size/
+    - gc_bias/
+    - hs_metrics/
+  - samtools/
+  - preseq/
+  - rustqc_summary.json
+  - CITATIONS.md
+</FileTree>
+
+Pass `--flat-output` to write everything directly into the output directory
+instead.
+
+## Targeted mode
+
+Passing `--targets targets.bed` switches the run into targeted mode and adds
+[`CollectHsMetrics`](/dna/picard/#collecthsmetrics). `--baits` defaults to the
+same intervals; give it separately when the capture baits differ from the
+regions you want reported.
+
+Intervals are merged on load. Overlapping targets would otherwise inflate the
+reported territory and count the same base twice.
+
+## Memory
+
+The depth engine holds one array of four bytes per base for each contig being
+processed, so the largest contig sets the cost per worker: roughly 1 GB for
+GRCh38 chr1. `--max-depth-workers` bounds how many are live at once, defaulting
+to a 4 GB budget divided by the largest contig. Raise it if you have the memory
+and want more parallelism; lower it on a shared machine.
+
+## JSON summary
+
+`--json-summary` writes a machine-readable summary 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. The
+coverage thresholds are a list rather than a map so that the order you asked
+for survives.