Skip to content

feat: add bamqc subcommand for generic genomic BAM/CRAM QC - #135

Open
BenjaminDEMAILLE wants to merge 1 commit into
seqeralabs:mainfrom
BenjaminDEMAILLE:feat/bamqc-mode
Open

feat: add bamqc subcommand for generic genomic BAM/CRAM QC#135
BenjaminDEMAILLE wants to merge 1 commit into
seqeralabs:mainfrom
BenjaminDEMAILLE:feat/bamqc-mode

Conversation

@BenjaminDEMAILLE

Copy link
Copy Markdown

Closes #128.

Adds rustqc bamqc — generic genomic QC in one streaming pass, no GTF required:

rustqc bamqc sample.bam --outdir results/
rustqc bamqc sample.cram --reference genome.fa --outdir results/

What it reports

Coverage depth and breadth (genome-wide and per contig), GC content, insert size, mapping quality, read lengths and base composition — the metrics side of Qualimap's bamqc mode.

Output names and layout match Qualimap, so MultiQC's Qualimap BamQC module picks them up unchanged:

genome_results.txt
raw_data_qualimapReport/coverage_histogram.txt
raw_data_qualimapReport/genome_fraction_coverage.txt
raw_data_qualimapReport/mapped_reads_gc-content_distribution.txt
raw_data_qualimapReport/mapped_reads_nucleotide_content.txt
raw_data_qualimapReport/insert_size_histogram.txt

Coverage is exact, not sampled

Each read contributes +1 over every aligned block (M/=/X/D); N skips leave a gap. Since input is coordinate-sorted, depth is finalised as the file streams past, so memory scales with pile-up depth rather than genome size — no per-base array.

Secondary alignments are excluded from coverage. Duplicates are included by default (Qualimap's default); --skip-duplicated excludes them.

Validation

Cross-checked against samtools depth on tests/data/test.bam, both modes, and pinned in integration tests:

Mode RustQC samtools
default mean 0.6038X, 12.16% ≥ 1X, 4.78% ≥ 5X samtools depth -a -J -g DUP → identical
--skip-duplicated mean 0.4325X, 11.55% ≥ 1X samtools depth -a -J → identical

(The two differ by exactly the 137 duplicate-flagged reads × 50 bp in the test file.)

Scope

This is the metrics subset, not all of Qualimap bamqc. Not included: the HTML report, per-window coverage tracks, clipping profile, homopolymer indel breakdown, PCR-bias/chimerism estimates, and coverage restricted to a target BED (--feature-file). The docs page states this explicitly. Happy to extend if you'd like any of those in the first cut — and the offer in the issue to validate against real Qualimap output on real data would be well worth taking up, since my validation here is against samtools depth rather than Qualimap itself.

Note this overlaps in machinery with #18 (single-pass CRAM QC): the depth tracker here is what a mosdepth-equivalent would build on.

  • cargo test — 204 lib + 20 integration tests pass
  • cargo fmt --check, cargo clippy -- -D warnings clean
  • New docs page docs/bamqc, AGENTS.md updated

🤖 Generated with Claude Code

RustQC's Qualimap module only covers rnaseq mode, so pipelines that also want
Qualimap's bamqc mode have to keep the Java tool around alongside RustQC.

Adds `rustqc bamqc`, a single streaming pass over a coordinate-sorted
BAM/CRAM that needs no annotation and reports the bamqc metrics: coverage
depth and breadth (genome-wide and per contig), GC content, insert size,
mapping quality, read lengths and base composition.

Coverage is exact, not sampled: each read contributes +1 over every aligned
block (M/=/X/D) while N skips leave a gap, and because the input is sorted,
depth is finalised as the file streams past. Memory scales with pile-up
depth, not genome size. Secondary alignments are excluded; duplicates are
included by default, matching Qualimap, with --skip-duplicated to exclude
them.

Output file names and layout match Qualimap so MultiQC's Qualimap BamQC
module picks the results up unchanged: genome_results.txt plus
raw_data_qualimapReport/{coverage_histogram, genome_fraction_coverage,
mapped_reads_gc-content_distribution, mapped_reads_nucleotide_content,
insert_size_histogram}.txt.

Both duplicate modes are cross-checked against samtools depth in integration
tests (mean coverage and genome fraction match to the reported precision).

This is the metrics subset, not all of Qualimap bamqc: no HTML report, no
per-window coverage track, no clipping profile, no target-BED restriction.
The docs say so explicitly.

Closes seqeralabs#128

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.

Feature request: support Qualimap bamqc mode (genomic BAM QC) in addition to rnaseq mode

1 participant