feat: add align subcommand for single-pass DNA alignment QC

[BenjaminDEMAILLE]

Addresses #18 (phases 1 and 2 of the proposal). Not marked Closes — phase 3 and the pipeline-side work remain.

rustqc align collapses three independent passes over the same CRAM/BAM into one:

rustqc align sample.cram \
  --reference genome.fa \
  --snp-bed SNP_GRCh38_hg38_wChr.bed \
  --by 500 \
  --outdir results/

Output	Equivalent to	Consumer
{sample}.stats / .flagstat / .idxstats	samtools stats etc.	MultiQC
{sample}.mosdepth.summary.txt	mosdepth summary	MultiQC
{sample}.mosdepth.global.dist.txt	mosdepth distribution	MultiQC
{sample}.regions.bed.gz	mosdepth per-window depth	MultiQC
{sample}.ngscheckmate.vcf.gz	bcftools mpileup | bcftools call	ncm.py

The samtools outputs reuse the accumulator the rna command already validates against samtools, so this PR is mostly the depth engine and the genotyping.

mosdepth output is byte-identical

Verified against mosdepth 0.3.x (mosdepth --by 500 md tests/data/test.bam) on the test dataset — summary.txt, global.dist.txt and the decompressed regions.bed.gz all diff clean. That required matching two mosdepth behaviours that are easy to miss:

• the <chrom>_region / total_region summary rows emitted whenever --by is used
• the if cum < 8e-5: continue cutoff in write_distribution, which skips the sparse tail (chr2's max depth is 40, but 39 is the deepest row mosdepth prints)

Both are pinned in an integration test and flagged in AGENTS.md so they don't get "cleaned up" later.

Depth engine

Exact, from CIGAR-aware start/end delta events: M/=/X/D contribute, N skips leave a gap. Since input is coordinate-sorted, every base before the current read's start is final and gets folded in as the file streams past — memory scales with pile-up depth, not genome size. Excluded records match mosdepth's default --flag 1796 (unmapped, secondary, QC-fail, duplicate).

Genotyping

For a sorted file the SNP panel is a sorted array, so the per-read cost is one comparison; only reads that overlap a site get CIGAR-resolved base extraction (handling soft clips, insertions, deletions and N skips — covered by unit tests).

Genotypes come from the alternate allele fraction (<0.15 → 0/0, 0.15–0.85 → 0/1, >0.85 → 1/1), which is all ncm.py distinguishes. FORMAT is GT:AD:DP.

The SNP BED must be the 6-column NGSCheckMate layout including ref/alt alleles. A shorter BED is rejected with an explicit error rather than having alleles inferred — inferring them per sample would let two samples disagree on REF/ALT and silently corrupt the correlation matrix.

Not covered

• --by takes a fixed window size only; a target BED (mosdepth's --by <bed>) is not supported, so WES target-restricted coverage still needs mosdepth
• No per-base.bed.gz, no .csi index for regions.bed.gz
• No BAQ recalculation, so allele counts can differ from bcftools mpileup in repetitive regions

The validation the issue asks for — comparing NGSCheckMate VCFs and the resulting ncm.py correlation matrices on real data — has not been done here; the test dataset has no real variation. That is the piece to run before anyone swaps this into a pipeline, and I'd suggest treating the genotyping as unvalidated until then.

Note this shares machinery with #128 (bamqc): both need an exact streaming depth tracker. If both land, the two are worth unifying.

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

🤖 Generated with Claude Code