From d07bc541efbb237939fb1e232637cb95e374ea4e Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 13 Aug 2026 10:09:26 +0200 Subject: [PATCH] feat: add align subcommand for single-pass DNA alignment QC Sarek reads each sample's CRAM 6-7 times during QC. This collapses three of those passes into one: samtools stats, mosdepth-equivalent depth, and the bcftools mpileup genotyping step that feeds NGSCheckMate all come out of a single streaming pass. The genotyping is nearly free. The NGSCheckMate panel is ~10K positions across ~3 Gb, so for a coordinate-sorted file the check is a single comparison for almost every read; only the reads that actually overlap a site need CIGAR-resolved base extraction. bcftools mpileup instead builds complete pileup columns genome-wide and runs a genotype likelihood model. Depth is exact, from CIGAR-aware start/end delta events (M/=/X/D contribute, N skips leave a gap). Because input is sorted, depth is finalised as the file streams past, so memory scales with pile-up depth rather than genome size. Records excluded: unmapped, secondary, QC-fail, duplicate (mosdepth's default --flag 1796). Outputs: {sample}.stats / .flagstat / .idxstats (reusing the accumulator the rna command already validates against samtools), .mosdepth.summary.txt, .mosdepth.global.dist.txt, .regions.bed.gz and .ngscheckmate.vcf.gz. All three depth outputs are byte-identical to mosdepth 0.3.x on the test dataset, including its quirks: the _region summary rows emitted with --by, and the 8e-5 cumulative cutoff that skips the sparse tail of the distribution. Genotypes come from the alternate allele fraction (<0.15 hom-ref, 0.15-0.85 het, >0.85 hom-alt) since ncm.py only needs those three states. The SNP BED must carry ref/alt alleles; a shorter BED is rejected rather than guessed at, because per-sample VCFs have to share a common allele set. Scope: phases 1 and 2 of the proposal. --by takes a fixed window size only (no target BED), there is no per-base.bed.gz or .csi index, and genotyping does no BAQ recalculation. The docs say so explicitly. Refs #18 Co-Authored-By: Claude Opus 5 (1M context) --- AGENTS.md | 18 +- docs/astro.config.mjs | 4 + docs/src/content/docs/align.mdx | 112 ++++++++ src/align/depth.rs | 416 +++++++++++++++++++++++++++++ src/align/mod.rs | 17 ++ src/align/output.rs | 268 +++++++++++++++++++ src/align/snp.rs | 456 ++++++++++++++++++++++++++++++++ src/cli.rs | 124 +++++++++ src/lib.rs | 1 + src/main.rs | 149 +++++++++++ tests/integration_test.rs | 123 +++++++++ 11 files changed, 1687 insertions(+), 1 deletion(-) create mode 100644 docs/src/content/docs/align.mdx create mode 100644 src/align/depth.rs create mode 100644 src/align/mod.rs create mode 100644 src/align/output.rs create mode 100644 src/align/snp.rs diff --git a/AGENTS.md b/AGENTS.md index 108414ac..61535619 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -62,6 +62,11 @@ 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) + align/ + mod.rs — Re-exports the align accumulators and writers + depth.rs — mosdepth-equivalent per-base and per-window depth + snp.rs — NGSCheckMate SNP panel parsing and allele counting + output.rs — mosdepth + NGSCheckMate VCF writers rna/ mod.rs — Re-exports all submodules (dupradar, featurecounts, rseqc, bam_flags, cpp_rng, preseq, qualimap) bam_flags.rs — BAM flag constants @@ -111,9 +116,11 @@ Nested module structure — top-level modules (`cli`, `config`, `io`, `gtf`, `rn 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;`). -The CLI uses a single subcommand: +The CLI has two subcommands: - `rustqc rna ... --gtf [OPTIONS]` +- `rustqc align [OPTIONS]` — single-pass DNA alignment QC (samtools stats, + mosdepth-compatible depth, NGSCheckMate genotyping); no annotation required A GTF gene annotation file (`--gtf`) is required. This runs all analyses: dupRadar duplicate rate analysis, featureCounts-compatible gene counting, @@ -274,6 +281,15 @@ forwarded to `count_reads()` as the `skip_dup_check: bool` parameter). ## Notes for Agents +- `rustqc align` output is byte-compatible with mosdepth on the test data, including + two quirks that must not be "cleaned up": the `_region` rows in + `mosdepth.summary.txt` when `--by` is used, and the `8e-5` cumulative cutoff that + skips the sparse tail of `mosdepth.global.dist.txt`. Depth excludes unmapped, + secondary, QC-fail and duplicate records (mosdepth's default `--flag 1796`). +- The NGSCheckMate SNP BED must be the 6-column layout with ref/alt alleles; a + shorter BED is rejected rather than guessed at, because per-sample VCFs have to + share a common allele set to be comparable. + - A `.pre-commit-config.yaml` is provided for local git hooks (fmt, clippy, file hygiene). Use [prek](https://github.com/j178/prek) (`prek install`) or the original [pre-commit](https://pre-commit.com/) to activate them. diff --git a/docs/astro.config.mjs b/docs/astro.config.mjs index 4e1c2f42..73361857 100644 --- a/docs/astro.config.mjs +++ b/docs/astro.config.mjs @@ -72,6 +72,10 @@ export default defineConfig({ { label: "Samtools", slug: "rna/samtools" }, ], }, + { + label: "DNA", + items: [{ label: "align", slug: "align" }], + }, { label: "About", items: [ diff --git a/docs/src/content/docs/align.mdx b/docs/src/content/docs/align.mdx new file mode 100644 index 00000000..b0fbb2a7 --- /dev/null +++ b/docs/src/content/docs/align.mdx @@ -0,0 +1,112 @@ +--- +title: align +description: Single-pass alignment QC for DNA pipelines — samtools stats, mosdepth-compatible depth, and NGSCheckMate genotyping from one pass over a CRAM/BAM. +--- + +import { Aside } from "@astrojs/starlight/components"; + +`rustqc align` computes in **one streaming pass** what DNA pipelines currently +get from three independent passes over the same CRAM/BAM: + +- `samtools stats` (plus `flagstat` and `idxstats`) +- `mosdepth` per-contig and per-window depth +- the `bcftools mpileup` genotyping step that feeds NGSCheckMate + +```bash +rustqc align sample.cram \ + --reference genome.fa \ + --snp-bed SNP_GRCh38_hg38_wChr.bed \ + --by 500 \ + --outdir results/ +``` + +Input must be coordinate-sorted. + +## Why the genotyping is nearly free + +The NGSCheckMate panel is ~10,000 positions across ~3 Gb. With a +coordinate-sorted file the panel is a sorted array, so the check for the vast +majority of reads is a single comparison; only the ~0.003% of reads that +actually overlap a site need CIGAR-resolved base extraction. `bcftools mpileup` +instead builds complete pileup columns genome-wide and runs a genotype +likelihood model, which is where the hours go. + +## Options + +| Flag | Description | +| ----------------- | -------------------------------------------------------------------- | +| `-r, --reference` | Reference FASTA (required for CRAM) | +| `--snp-bed` | NGSCheckMate SNP BED; enables genotyping | +| `--by` | Depth window size in bases (default 500) | +| `-Q, --mapq` | MAPQ cutoff for stats and genotyping (default 30) | +| `--min-bq` | Minimum base quality for genotyping (default 13) | +| `-t, --threads` | htslib decompression threads | + +## Output files + +| File | Equivalent to | Consumer | +| --------------------------------- | -------------------------------- | --------- | +| `{sample}.stats` | `samtools stats` | MultiQC | +| `{sample}.flagstat` | `samtools flagstat` | MultiQC | +| `{sample}.idxstats` | `samtools idxstats` | MultiQC | +| `{sample}.mosdepth.summary.txt` | `mosdepth` summary | MultiQC | +| `{sample}.mosdepth.global.dist.txt` | `mosdepth` global distribution | MultiQC | +| `{sample}.regions.bed.gz` | `mosdepth` per-window depth | MultiQC | +| `{sample}.ngscheckmate.vcf.gz` | `bcftools mpileup \| call` | `ncm.py` | + +### SNP BED format + +The 6-column NGSCheckMate layout, as shipped with NGSCheckMate: + +``` +chr17 46549406 46549407 rs201103889 A C +chr1 152308305 152308306 rs2184953 T C +``` + +Reference and alternate alleles are **required** — without them, per-sample +VCFs could not be compared against a common set of alleles. A BED without them +is rejected with an explicit error rather than guessed at. + +### Genotypes + +`ncm.py` only needs to tell hom-ref, het and hom-alt apart, so genotypes come +from the alternate allele fraction rather than a likelihood model: + +| Alt fraction | GT | +| ------------- | ----- | +| `< 0.15` | `0/0` | +| `0.15 – 0.85` | `0/1` | +| `> 0.85` | `1/1` | +| no coverage | `./.` | + +`FORMAT` is `GT:AD:DP`. + +## Depth semantics + +Depth is exact, computed from CIGAR-aware start/end delta events: `M`/`=`/`X`/`D` +contribute, `N` skips leave a gap. Because input is coordinate-sorted, depth is +finalised as the file streams past, so memory scales with pile-up depth rather +than genome size. + +Records excluded from depth: unmapped, secondary, QC-fail and duplicate — +mosdepth's default `--flag 1796`. + +## Validation + +Against `mosdepth 0.3.x` on the test dataset, all three depth outputs are +**byte-identical**: `mosdepth.summary.txt`, `mosdepth.global.dist.txt` and the +decompressed `regions.bed.gz`. This includes mosdepth's quirks — the +`_region` summary rows emitted when `--by` is used, and the sparse tail +skipped from the distribution when the cumulative fraction is below `8e-5`. + +The `samtools stats` output comes from the same accumulator the `rna` command +uses, which is validated against samtools separately. + + diff --git a/src/align/depth.rs b/src/align/depth.rs new file mode 100644 index 00000000..3245cb69 --- /dev/null +++ b/src/align/depth.rs @@ -0,0 +1,416 @@ +//! mosdepth-equivalent per-base and per-window depth. +//! +//! Depth is computed exactly from CIGAR-aware start/end delta events. Because +//! the input is coordinate-sorted, every base before the current read's start +//! is final and is folded into the per-contig summary and the current window +//! as the file streams past, so memory stays proportional to pile-up depth +//! rather than to contig length. + +use crate::rna::bam_flags::{BAM_FDUP, BAM_FQCFAIL, BAM_FSECONDARY, BAM_FUNMAP}; +use rust_htslib::bam::{self, record::Cigar}; +use std::collections::BTreeMap; + +/// Highest depth tracked in the cumulative distribution. +/// +/// Deeper positions are folded into the top bin, matching mosdepth's bounded +/// coverage array. +pub const MAX_DIST_DEPTH: u32 = 1000; + +/// Records excluded from depth, matching mosdepth's default `--flag 1796` +/// (UNMAP, SECONDARY, QCFAIL, DUP). +const EXCLUDED_FLAGS: u16 = BAM_FUNMAP | BAM_FSECONDARY | BAM_FQCFAIL | BAM_FDUP; + +/// Per-contig depth results. +#[derive(Debug, Clone)] +pub struct ContigDepth { + /// Contig name. + pub name: String, + /// Contig length from the alignment header. + pub length: u64, + /// Σ depth over the contig (aligned bases). + pub total_bases: u64, + /// Minimum per-base depth (0 unless the contig is fully covered). + pub min_depth: u32, + /// Maximum per-base depth. + pub max_depth: u32, + /// depth (capped at [`MAX_DIST_DEPTH`]) → number of bases at that depth. + pub hist: BTreeMap, +} + +impl ContigDepth { + /// Mean depth across the whole contig, uncovered bases counted as zero. + pub fn mean(&self) -> f64 { + if self.length == 0 { + 0.0 + } else { + self.total_bases as f64 / self.length as f64 + } + } +} + +/// One fixed-size window of the reference. +#[derive(Debug, Clone)] +pub struct Window { + /// Contig name. + pub chrom: String, + /// Window start (0-based). + pub start: u64, + /// Window end (exclusive). + pub end: u64, + /// Mean depth across the window. + pub mean: f64, +} + +/// Streaming depth accumulator over a coordinate-sorted alignment file. +#[derive(Debug)] +pub struct DepthAccum { + /// Contig names and lengths, in header order. + contigs: Vec<(String, u64)>, + /// Window size in bases. + window_size: u64, + /// Finished per-contig results. + results: Vec, + /// Emitted windows, in reference order. + windows: Vec, + + // --- state for the contig being processed --- + current_tid: i32, + events: BTreeMap, + last_pos: u64, + depth: i64, + hist: BTreeMap, + total_bases: u64, + max_depth: u32, + /// Σ depth per window of the contig being processed. + window_sums: Vec, +} + +impl DepthAccum { + /// Create an accumulator for a file with the given contigs. + /// + /// # Arguments + /// * `contigs` - Reference names and lengths from the alignment header + /// * `window_size` - Window size in bases (mosdepth's `--by`) + pub fn new(contigs: Vec<(String, u64)>, window_size: u64) -> Self { + Self { + contigs, + window_size: window_size.max(1), + results: Vec::new(), + windows: Vec::new(), + current_tid: -1, + events: BTreeMap::new(), + last_pos: 0, + depth: 0, + hist: BTreeMap::new(), + total_bases: 0, + max_depth: 0, + window_sums: Vec::new(), + } + } + + /// Add one alignment record. + pub fn process_read(&mut self, record: &bam::Record) { + if record.flags() & EXCLUDED_FLAGS != 0 { + return; + } + let tid = record.tid(); + if tid < 0 { + return; + } + if tid != self.current_tid { + self.finish_contig(); + self.current_tid = tid; + } + + let start = record.pos() as u64; + self.advance_to(start); + + let mut ref_pos = start; + for op in record.cigar().iter() { + match *op { + Cigar::Match(len) | Cigar::Equal(len) | Cigar::Diff(len) | Cigar::Del(len) => { + let len = len as u64; + *self.events.entry(ref_pos).or_insert(0) += 1; + *self.events.entry(ref_pos + len).or_insert(0) -= 1; + ref_pos += len; + } + Cigar::RefSkip(len) => ref_pos += len as u64, + _ => {} + } + } + } + + /// Fold all depth up to (but excluding) `pos` into the current contig. + fn advance_to(&mut self, pos: u64) { + while let Some((&event_pos, &delta)) = self.events.iter().next() { + if event_pos >= pos { + break; + } + if event_pos > self.last_pos { + let span = event_pos - self.last_pos; + let depth = self.depth.max(0) as u32; + self.record(depth, span); + } + self.depth += delta; + self.last_pos = event_pos; + self.events.remove(&event_pos); + } + if pos > self.last_pos { + let span = pos - self.last_pos; + let depth = self.depth.max(0) as u32; + self.record(depth, span); + self.last_pos = pos; + } + } + + /// Record `span` consecutive bases at `depth`, splitting across windows. + fn record(&mut self, depth: u32, span: u64) { + if span == 0 { + return; + } + if depth > 0 { + *self.hist.entry(depth.min(MAX_DIST_DEPTH)).or_insert(0) += span; + self.total_bases += depth as u64 * span; + self.max_depth = self.max_depth.max(depth); + } + + if depth == 0 { + return; + } + + // Split the span across window boundaries + let mut pos = self.last_pos; + let end = self.last_pos + span; + while pos < end { + let idx = (pos / self.window_size) as usize; + let window_end = (idx as u64 + 1) * self.window_size; + let chunk_end = window_end.min(end); + if idx >= self.window_sums.len() { + self.window_sums.resize(idx + 1, 0); + } + self.window_sums[idx] += depth as u64 * (chunk_end - pos); + pos = chunk_end; + } + } + + /// Emit every window of the contig being processed, including empty ones. + /// + /// mosdepth tiles the whole reference when `--by ` is a fixed window + /// size, so windows with no coverage are written with a mean of 0.00. + fn flush_windows(&mut self) { + if self.current_tid < 0 || (self.current_tid as usize) >= self.contigs.len() { + self.window_sums.clear(); + return; + } + let (name, contig_len) = self.contigs[self.current_tid as usize].clone(); + let num_windows = contig_len.div_ceil(self.window_size); + for idx in 0..num_windows { + let start = idx * self.window_size; + let end = ((idx + 1) * self.window_size).min(contig_len); + let span = end - start; + if span == 0 { + continue; + } + let sum = self.window_sums.get(idx as usize).copied().unwrap_or(0); + self.windows.push(Window { + chrom: name.clone(), + start, + end, + mean: sum as f64 / span as f64, + }); + } + self.window_sums.clear(); + } + + /// Finalise the contig currently being processed. + fn finish_contig(&mut self) { + if self.current_tid >= 0 { + let last_event = self.events.keys().next_back().copied().unwrap_or(0); + self.advance_to(last_event + 1); + self.flush_windows(); + + if (self.current_tid as usize) < self.contigs.len() { + let (name, length) = self.contigs[self.current_tid as usize].clone(); + let covered: u64 = self.hist.values().sum(); + let min_depth = if covered < length { + 0 + } else { + self.hist.keys().next().copied().unwrap_or(0) + }; + self.results.push(ContigDepth { + name, + length, + total_bases: self.total_bases, + min_depth, + max_depth: self.max_depth, + hist: std::mem::take(&mut self.hist), + }); + } + } + + self.events.clear(); + self.hist.clear(); + self.depth = 0; + self.last_pos = 0; + self.total_bases = 0; + self.max_depth = 0; + self.window_sums.clear(); + } + + /// Finalise everything and return the per-contig results and windows. + /// + /// Contigs with no alignments are reported with zero coverage so the + /// summary lists every reference in the header, as mosdepth does. + pub fn finish(mut self) -> (Vec, Vec) { + self.finish_contig(); + + let seen: std::collections::HashSet<&str> = + self.results.iter().map(|r| r.name.as_str()).collect(); + let mut results = Vec::with_capacity(self.contigs.len()); + for (name, length) in &self.contigs { + if seen.contains(name.as_str()) { + continue; + } + results.push(ContigDepth { + name: name.clone(), + length: *length, + total_bases: 0, + min_depth: 0, + max_depth: 0, + hist: BTreeMap::new(), + }); + let num_windows = length.div_ceil(self.window_size); + for idx in 0..num_windows { + let start = idx * self.window_size; + let end = ((idx + 1) * self.window_size).min(*length); + if end > start { + self.windows.push(Window { + chrom: name.clone(), + start, + end, + mean: 0.0, + }); + } + } + } + // Keep header order + let mut all = self.results; + all.extend(results); + let order: std::collections::HashMap<&str, usize> = self + .contigs + .iter() + .enumerate() + .map(|(i, (name, _))| (name.as_str(), i)) + .collect(); + all.sort_by_key(|r| order.get(r.name.as_str()).copied().unwrap_or(usize::MAX)); + + (all, self.windows) + } +} + +// =================================================================== +// Tests +// =================================================================== + +#[cfg(test)] +mod tests { + use super::*; + use rust_htslib::bam::Read as BamRead; + + fn accumulate( + sam: &str, + contigs: Vec<(String, u64)>, + by: u64, + ) -> (Vec, Vec) { + use std::sync::atomic::{AtomicU64, Ordering}; + static COUNTER: AtomicU64 = AtomicU64::new(0); + let id = COUNTER.fetch_add(1, Ordering::Relaxed); + let path = std::env::temp_dir().join(format!( + "rustqc_depth_test_{:?}_{}.sam", + std::thread::current().id(), + id + )); + std::fs::write(&path, sam).unwrap(); + + let mut reader = bam::Reader::from_path(&path).unwrap(); + let mut accum = DepthAccum::new(contigs, by); + let mut record = bam::Record::new(); + while let Some(res) = reader.read(&mut record) { + res.unwrap(); + accum.process_read(&record); + } + let _ = std::fs::remove_file(&path); + accum.finish() + } + + #[test] + fn test_windows_and_contig_summary() { + // Contig of 100 bases, window size 50. + // r1 covers 1..=10, r2 covers 6..=15 (1-based) -> window 0 gets + // 20 aligned bases, window 1 gets none. + let sam = "\ +@HD\tVN:1.6\tSO:coordinate\n\ +@SQ\tSN:chr1\tLN:100\n\ +r1\t0\tchr1\t1\t60\t10M\t*\t0\t0\tACGTACGTAC\tIIIIIIIIII\n\ +r2\t0\tchr1\t6\t60\t10M\t*\t0\t0\tACGTACGTAC\tIIIIIIIIII\n"; + let (contigs, windows) = accumulate(sam, vec![("chr1".to_string(), 100)], 50); + + assert_eq!(contigs.len(), 1); + assert_eq!(contigs[0].total_bases, 20); + assert!((contigs[0].mean() - 0.2).abs() < 1e-9); + assert_eq!(contigs[0].max_depth, 2); + assert_eq!(contigs[0].min_depth, 0, "contig is not fully covered"); + assert_eq!(contigs[0].hist.get(&1), Some(&10)); + assert_eq!(contigs[0].hist.get(&2), Some(&5)); + + // mosdepth tiles the whole contig, so both windows are emitted + assert_eq!(windows.len(), 2); + assert_eq!(windows[0].start, 0); + assert_eq!(windows[0].end, 50); + assert!((windows[0].mean - 20.0 / 50.0).abs() < 1e-9); + assert_eq!(windows[1].start, 50); + assert_eq!(windows[1].mean, 0.0, "empty windows are emitted as 0"); + } + + #[test] + fn test_window_split_across_boundary() { + // 20-base read starting at 1-based 41 spans windows [0,50) and [50,100) + let sam = "\ +@HD\tVN:1.6\tSO:coordinate\n\ +@SQ\tSN:chr1\tLN:100\n\ +r1\t0\tchr1\t41\t60\t20M\t*\t0\t0\tACGTACGTACACGTACGTAC\tIIIIIIIIIIIIIIIIIIII\n"; + let (_, windows) = accumulate(sam, vec![("chr1".to_string(), 100)], 50); + + assert_eq!(windows.len(), 2); + // 10 bases in each window + assert!((windows[0].mean - 10.0 / 50.0).abs() < 1e-9); + assert!((windows[1].mean - 10.0 / 50.0).abs() < 1e-9); + } + + #[test] + fn test_excluded_flags_and_empty_contigs() { + // r1 is a duplicate, r2 secondary, r3 unmapped: none contribute. + // chr2 has no reads at all and must still be reported. + let sam = "\ +@HD\tVN:1.6\tSO:coordinate\n\ +@SQ\tSN:chr1\tLN:100\n\ +@SQ\tSN:chr2\tLN:200\n\ +r1\t1024\tchr1\t1\t60\t10M\t*\t0\t0\tACGTACGTAC\tIIIIIIIIII\n\ +r2\t256\tchr1\t1\t60\t10M\t*\t0\t0\tACGTACGTAC\tIIIIIIIIII\n\ +r3\t4\tchr1\t1\t0\t*\t*\t0\t0\tACGTACGTAC\tIIIIIIIIII\n"; + let (contigs, windows) = accumulate( + sam, + vec![("chr1".to_string(), 100), ("chr2".to_string(), 200)], + 50, + ); + + assert_eq!(contigs.len(), 2, "every header contig is reported"); + assert_eq!(contigs[0].name, "chr1"); + assert_eq!(contigs[0].total_bases, 0); + assert_eq!(contigs[1].name, "chr2"); + assert_eq!(contigs[1].total_bases, 0); + // 2 windows over chr1 (100 bp) + 4 over chr2 (200 bp), all empty + assert_eq!(windows.len(), 6); + assert!(windows.iter().all(|w| w.mean == 0.0)); + } +} diff --git a/src/align/mod.rs b/src/align/mod.rs new file mode 100644 index 00000000..4207a8aa --- /dev/null +++ b/src/align/mod.rs @@ -0,0 +1,17 @@ +//! Single-pass alignment QC for DNA pipelines (`rustqc align`). +//! +//! Replaces three independent passes over a CRAM/BAM — `samtools stats`, +//! `mosdepth`, and the `bcftools mpileup` genotyping step feeding +//! NGSCheckMate — with one streaming pass that computes all three. +//! +//! The genotyping component is close to free: for a coordinate-sorted file the +//! SNP panel is a sorted array, so most reads cost a single comparison, and +//! only the tiny fraction overlapping a site needs CIGAR-resolved base +//! extraction. + +pub mod depth; +pub mod output; +pub mod snp; + +pub use depth::{ContigDepth, DepthAccum, Window}; +pub use snp::{parse_snp_bed, AlleleCounts, SnpAccum, SnpPanel}; diff --git a/src/align/output.rs b/src/align/output.rs new file mode 100644 index 00000000..c83b8f9c --- /dev/null +++ b/src/align/output.rs @@ -0,0 +1,268 @@ +//! mosdepth- and NGSCheckMate-compatible output files. + +use std::io::Write; +use std::path::Path; + +use anyhow::{Context, Result}; +use log::debug; + +use super::depth::{ContigDepth, Window, MAX_DIST_DEPTH}; +use super::snp::{SnpAccum, SnpPanel}; + +/// Write `{prefix}.mosdepth.summary.txt`. +/// +/// Columns are mosdepth's: `chrom length bases mean min max`, one row per +/// contig plus a `total` row. Because `--by ` tiles the whole reference, +/// each contig also gets the `_region` row mosdepth emits for the +/// region set, with the same values. +pub fn write_summary(contigs: &[ContigDepth], path: &Path) -> Result<()> { + let mut out = std::fs::File::create(path) + .with_context(|| format!("Failed to create {}", path.display()))?; + writeln!(out, "chrom\tlength\tbases\tmean\tmin\tmax")?; + + let mut total_length = 0u64; + let mut total_bases = 0u64; + let mut total_max = 0u32; + for contig in contigs { + for name in [contig.name.clone(), format!("{}_region", contig.name)] { + writeln!( + out, + "{}\t{}\t{}\t{:.2}\t{}\t{}", + name, + contig.length, + contig.total_bases, + contig.mean(), + contig.min_depth, + contig.max_depth + )?; + } + total_length += contig.length; + total_bases += contig.total_bases; + total_max = total_max.max(contig.max_depth); + } + + let total_mean = if total_length == 0 { + 0.0 + } else { + total_bases as f64 / total_length as f64 + }; + for name in ["total", "total_region"] { + writeln!( + out, + "{}\t{}\t{}\t{:.2}\t{}\t{}", + name, total_length, total_bases, total_mean, 0, total_max + )?; + } + + debug!("Wrote mosdepth summary to {}", path.display()); + Ok(()) +} + +/// Write `{prefix}.mosdepth.global.dist.txt`. +/// +/// For each contig (and `total`), emits `chrom depth proportion` from the +/// deepest observed level down to 0, where `proportion` is the fraction of +/// that contig's bases covered at **at least** that depth. mosdepth's format. +pub fn write_global_dist(contigs: &[ContigDepth], path: &Path) -> Result<()> { + let mut out = std::fs::File::create(path) + .with_context(|| format!("Failed to create {}", path.display()))?; + + let mut total_hist: std::collections::BTreeMap = std::collections::BTreeMap::new(); + let mut total_length = 0u64; + + for contig in contigs { + write_dist_block(&mut out, &contig.name, &contig.hist, contig.length)?; + for (&depth, &bases) in &contig.hist { + *total_hist.entry(depth).or_insert(0) += bases; + } + total_length += contig.length; + } + write_dist_block(&mut out, "total", &total_hist, total_length)?; + + debug!("Wrote mosdepth global distribution to {}", path.display()); + Ok(()) +} + +/// Emit the cumulative distribution rows for one contig. +fn write_dist_block( + out: &mut impl Write, + name: &str, + hist: &std::collections::BTreeMap, + length: u64, +) -> Result<()> { + if length == 0 { + return Ok(()); + } + let max_depth = hist + .keys() + .next_back() + .copied() + .unwrap_or(0) + .min(MAX_DIST_DEPTH); + let mut cumulative = 0u64; + for depth in (0..=max_depth).rev() { + cumulative += hist.get(&depth).copied().unwrap_or(0); + let proportion = if depth == 0 { + 1.0 + } else { + cumulative as f64 / length as f64 + }; + // mosdepth skips the sparse tail at the top of the distribution + // (`if cum < 8e-5: continue` in its write_distribution) + if proportion < 8e-5 { + continue; + } + writeln!(out, "{name}\t{depth}\t{proportion:.2}")?; + } + Ok(()) +} + +/// Write `{prefix}.regions.bed.gz` (BGZF-compressed per-window depth). +pub fn write_regions(windows: &[Window], path: &Path) -> Result<()> { + use rust_htslib::bgzf; + + let mut writer = bgzf::Writer::from_path(path) + .with_context(|| format!("Failed to create {}", path.display()))?; + for window in windows { + let line = format!( + "{}\t{}\t{}\t{:.2}\n", + window.chrom, window.start, window.end, window.mean + ); + writer + .write_all(line.as_bytes()) + .with_context(|| format!("Failed to write {}", path.display()))?; + } + + debug!( + "Wrote {} depth windows to {}", + windows.len(), + path.display() + ); + Ok(()) +} + +/// Write the NGSCheckMate VCF (BGZF-compressed). +/// +/// Emits `GT:AD:DP` per site, with the genotype derived from the alternate +/// allele fraction. `ncm.py` only needs to tell hom-ref, het and hom-alt +/// apart, so no genotype likelihood model is involved. +pub fn write_ngscheckmate_vcf( + panel: &SnpPanel, + accum: &SnpAccum, + sample_name: &str, + contig_order: &[(String, u64)], + path: &Path, +) -> Result<()> { + use rust_htslib::bgzf; + + let mut writer = bgzf::Writer::from_path(path) + .with_context(|| format!("Failed to create {}", path.display()))?; + + let mut header = String::new(); + header.push_str("##fileformat=VCFv4.2\n"); + header.push_str("##source=rustqc align\n"); + for (name, length) in contig_order { + header.push_str(&format!("##contig=\n")); + } + header.push_str("##INFO=\n"); + header.push_str("##FORMAT=\n"); + header.push_str( + "##FORMAT=\n", + ); + header.push_str("##FORMAT=\n"); + header.push_str(&format!( + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t{sample_name}\n" + )); + writer + .write_all(header.as_bytes()) + .with_context(|| format!("Failed to write {}", path.display()))?; + + // Emit in header contig order, positions ascending, so the file is sorted + for (chrom, _) in contig_order { + let (Some(sites), Some(counts)) = (panel.by_chrom.get(chrom), accum.counts_for(chrom)) + else { + continue; + }; + for (site, count) in sites.iter().zip(counts.iter()) { + let line = format!( + "{}\t{}\t{}\t{}\t{}\t.\t.\tDP={}\tGT:AD:DP\t{}:{},{}:{}\n", + chrom, + site.pos + 1, + site.id, + site.ref_base as char, + site.alt_base as char, + count.depth(), + count.genotype(), + count.ref_count, + count.alt_count, + count.depth(), + ); + writer + .write_all(line.as_bytes()) + .with_context(|| format!("Failed to write {}", path.display()))?; + } + } + + debug!("Wrote NGSCheckMate VCF to {}", path.display()); + Ok(()) +} + +// =================================================================== +// Tests +// =================================================================== + +#[cfg(test)] +mod tests { + use super::*; + use std::collections::BTreeMap; + + #[test] + fn test_dist_block_is_cumulative() { + // 100-base contig: 10 bases at depth 1, 5 at depth 2 + let mut hist = BTreeMap::new(); + hist.insert(1u32, 10u64); + hist.insert(2u32, 5u64); + + let mut out: Vec = Vec::new(); + write_dist_block(&mut out, "chr1", &hist, 100).unwrap(); + let text = String::from_utf8(out).unwrap(); + + // depth 2 -> 5/100, depth 1 -> 15/100, depth 0 -> 1.00 + assert_eq!(text, "chr1\t2\t0.05\nchr1\t1\t0.15\nchr1\t0\t1.00\n"); + } + + #[test] + fn test_summary_has_total_row() { + let contigs = vec![ + ContigDepth { + name: "chr1".to_string(), + length: 100, + total_bases: 20, + min_depth: 0, + max_depth: 2, + hist: BTreeMap::new(), + }, + ContigDepth { + name: "chr2".to_string(), + length: 100, + total_bases: 0, + min_depth: 0, + max_depth: 0, + hist: BTreeMap::new(), + }, + ]; + let path = std::env::temp_dir().join(format!( + "rustqc_summary_test_{:?}.txt", + std::thread::current().id() + )); + write_summary(&contigs, &path).unwrap(); + let text = std::fs::read_to_string(&path).unwrap(); + let _ = std::fs::remove_file(&path); + + assert!(text.starts_with("chrom\tlength\tbases\tmean\tmin\tmax\n")); + assert!(text.contains("chr1\t100\t20\t0.20\t0\t2\n")); + assert!(text.contains("chr1_region\t100\t20\t0.20\t0\t2\n")); + assert!(text.contains("total\t200\t20\t0.10\t0\t2\n")); + assert!(text.contains("total_region\t200\t20\t0.10\t0\t2\n")); + } +} diff --git a/src/align/snp.rs b/src/align/snp.rs new file mode 100644 index 00000000..0efa8f7a --- /dev/null +++ b/src/align/snp.rs @@ -0,0 +1,456 @@ +//! Targeted genotyping at known SNP positions (NGSCheckMate). +//! +//! NGSCheckMate needs allele counts at ~10K known positions, which +//! `bcftools mpileup` currently supplies by building complete pileup columns +//! across the whole genome. Only the reads overlapping those positions matter, +//! so the work folds into an existing streaming pass for almost nothing: a +//! pointer advance per read, and CIGAR-resolved base extraction for the +//! ~0.003% of reads that actually overlap a site. + +use std::collections::HashMap; +use std::io::BufRead; + +use anyhow::{Context, Result}; +use log::debug; +use rust_htslib::bam::{self, record::Cigar}; + +use crate::rna::bam_flags::{BAM_FDUP, BAM_FQCFAIL, BAM_FSECONDARY, BAM_FUNMAP}; + +/// One SNP site to genotype. +#[derive(Debug, Clone)] +pub struct SnpSite { + /// 0-based reference position. + pub pos: u64, + /// Variant identifier (BED column 4), used as the VCF ID field. + pub id: String, + /// Reference allele (BED column 5). + pub ref_base: u8, + /// Alternate allele (BED column 6). + pub alt_base: u8, +} + +/// SNP sites grouped by chromosome, each list sorted by position. +#[derive(Debug, Default)] +pub struct SnpPanel { + /// Chromosome name → sorted sites. + pub by_chrom: HashMap>, + /// Total number of sites loaded. + pub num_sites: usize, +} + +/// Parse an NGSCheckMate SNP BED file (plain or gzip-compressed). +/// +/// The file must have the 6-column NGSCheckMate layout: +/// `chrom start end id ref alt`. The reference and alternate alleles are +/// required — without them, samples could not be compared against a common +/// set of alleles. +pub fn parse_snp_bed(path: &str) -> Result { + let reader = crate::io::open_reader(path) + .with_context(|| format!("Failed to open SNP BED file: {path}"))?; + + let mut by_chrom: HashMap> = HashMap::new(); + let mut num_sites = 0usize; + + for (lineno, line) in reader.lines().enumerate() { + let line = line.with_context(|| format!("Failed to read line from SNP BED: {path}"))?; + let trimmed = line.trim_end(); + if trimmed.is_empty() || trimmed.starts_with('#') || trimmed.starts_with("track") { + continue; + } + + let fields: Vec<&str> = trimmed.split('\t').collect(); + anyhow::ensure!( + fields.len() >= 6, + "Malformed SNP BED '{}' at line {}: expected the 6-column NGSCheckMate layout \ + (chrom, start, end, id, ref, alt), found {} column(s)", + path, + lineno + 1, + fields.len() + ); + + let end: u64 = fields[2].parse().with_context(|| { + format!( + "Malformed SNP BED '{}' at line {}: invalid end position '{}'", + path, + lineno + 1, + fields[2] + ) + })?; + anyhow::ensure!( + end > 0, + "Malformed SNP BED '{}' at line {}: end position must be greater than 0", + path, + lineno + 1 + ); + + let ref_base = single_base(fields[4], path, lineno + 1, "reference")?; + let alt_base = single_base(fields[5], path, lineno + 1, "alternate")?; + + by_chrom + .entry(fields[0].to_string()) + .or_default() + .push(SnpSite { + // BED is half-open, so the SNP base is the last one in the interval + pos: end - 1, + id: fields[3].to_string(), + ref_base, + alt_base, + }); + num_sites += 1; + } + + anyhow::ensure!(num_sites > 0, "No SNP sites found in BED file '{}'", path); + + for sites in by_chrom.values_mut() { + sites.sort_unstable_by_key(|s| s.pos); + } + + debug!("Loaded {num_sites} SNP sites from {path}"); + Ok(SnpPanel { + by_chrom, + num_sites, + }) +} + +/// Validate and upper-case a single-base allele field. +fn single_base(field: &str, path: &str, lineno: usize, which: &str) -> Result { + let bytes = field.as_bytes(); + anyhow::ensure!( + bytes.len() == 1 && matches!(bytes[0].to_ascii_uppercase(), b'A' | b'C' | b'G' | b'T'), + "Malformed SNP BED '{}' at line {}: {} allele must be a single A/C/G/T base, found '{}'", + path, + lineno, + which, + field + ); + Ok(bytes[0].to_ascii_uppercase()) +} + +// =================================================================== +// Allele counting +// =================================================================== + +/// Per-site allele counts. +#[derive(Debug, Clone, Copy, Default)] +pub struct AlleleCounts { + /// Reads supporting the reference allele. + pub ref_count: u32, + /// Reads supporting the alternate allele. + pub alt_count: u32, + /// Reads carrying any other base at the site. + pub other_count: u32, +} + +impl AlleleCounts { + /// Total reads observed at the site. + pub fn depth(&self) -> u32 { + self.ref_count + self.alt_count + self.other_count + } + + /// Alternate allele fraction over ref + alt reads. + pub fn alt_fraction(&self) -> f64 { + let informative = self.ref_count + self.alt_count; + if informative == 0 { + 0.0 + } else { + self.alt_count as f64 / informative as f64 + } + } + + /// Genotype call from the alternate allele fraction. + /// + /// NGSCheckMate only needs to tell hom-ref, het and hom-alt apart, so the + /// simple fraction thresholds used here stand in for a full genotype + /// likelihood model. + pub fn genotype(&self) -> &'static str { + if self.depth() == 0 { + return "./."; + } + let frac = self.alt_fraction(); + if frac < 0.15 { + "0/0" + } else if frac <= 0.85 { + "0/1" + } else { + "1/1" + } + } +} + +/// Accumulates allele counts at the panel's sites during a streaming pass. +#[derive(Debug)] +pub struct SnpAccum { + /// Counts per chromosome, parallel to `SnpPanel::by_chrom` entries. + counts: HashMap>, + /// Minimum base quality for a base to be counted. + min_base_quality: u8, + /// Minimum MAPQ for a read to be considered. + min_mapq: u8, +} + +impl SnpAccum { + /// Create an accumulator for the given panel. + pub fn new(panel: &SnpPanel, min_mapq: u8, min_base_quality: u8) -> Self { + let counts = panel + .by_chrom + .iter() + .map(|(chrom, sites)| (chrom.clone(), vec![AlleleCounts::default(); sites.len()])) + .collect(); + Self { + counts, + min_base_quality, + min_mapq, + } + } + + /// Add one alignment record's base calls at any overlapping SNP sites. + /// + /// Records that are unmapped, secondary, QC-failed, duplicate-flagged or + /// below the MAPQ cutoff are ignored. + pub fn process_read(&mut self, record: &bam::Record, chrom: &str, panel: &SnpPanel) { + let flags = record.flags(); + if flags & (BAM_FUNMAP | BAM_FSECONDARY | BAM_FQCFAIL | BAM_FDUP) != 0 { + return; + } + if record.mapq() < self.min_mapq { + return; + } + + let (Some(sites), Some(counts)) = (panel.by_chrom.get(chrom), self.counts.get_mut(chrom)) + else { + return; + }; + if sites.is_empty() { + return; + } + + let start = record.pos() as u64; + // Fast reject: almost every read lands past the current site pointer + let first = sites.partition_point(|s| s.pos < start); + if first >= sites.len() { + return; + } + + let seq = record.seq(); + let quals = record.qual(); + let mut ref_pos = start; + let mut read_pos = 0usize; + + for op in record.cigar().iter() { + match *op { + Cigar::Match(len) | Cigar::Equal(len) | Cigar::Diff(len) => { + let len = len as u64; + let block_end = ref_pos + len; + for (i, site) in sites.iter().enumerate().skip(first) { + if site.pos >= block_end { + break; + } + if site.pos < ref_pos { + continue; + } + let offset = (site.pos - ref_pos) as usize; + let idx = read_pos + offset; + if idx >= seq.len() { + continue; + } + if quals.get(idx).copied().unwrap_or(0) < self.min_base_quality { + continue; + } + let base = seq[idx].to_ascii_uppercase(); + let entry = &mut counts[i]; + if base == site.ref_base { + entry.ref_count += 1; + } else if base == site.alt_base { + entry.alt_count += 1; + } else { + entry.other_count += 1; + } + } + ref_pos = block_end; + read_pos += len as usize; + } + Cigar::Ins(len) | Cigar::SoftClip(len) => read_pos += len as usize, + Cigar::Del(len) | Cigar::RefSkip(len) => ref_pos += len as u64, + _ => {} + } + } + } + + /// Counts for one chromosome, in panel order. + pub fn counts_for(&self, chrom: &str) -> Option<&[AlleleCounts]> { + self.counts.get(chrom).map(|v| v.as_slice()) + } + + /// Number of sites with at least one supporting read. + pub fn covered_sites(&self) -> usize { + self.counts + .values() + .flat_map(|v| v.iter()) + .filter(|c| c.depth() > 0) + .count() + } +} + +// =================================================================== +// Tests +// =================================================================== + +#[cfg(test)] +mod tests { + use super::*; + use rust_htslib::bam::Read as BamRead; + + fn write_temp(content: &str, ext: &str) -> std::path::PathBuf { + use std::sync::atomic::{AtomicU64, Ordering}; + static COUNTER: AtomicU64 = AtomicU64::new(0); + let id = COUNTER.fetch_add(1, Ordering::Relaxed); + let path = std::env::temp_dir().join(format!( + "rustqc_snp_test_{:?}_{}.{}", + std::thread::current().id(), + id, + ext + )); + std::fs::write(&path, content).unwrap(); + path + } + + #[test] + fn test_parse_snp_bed() { + let path = write_temp( + "chr1\t99\t100\trs1\tA\tG\n\ + chr1\t199\t200\trs2\tC\tT\n\ + chr2\t9\t10\trs3\tG\tA\n", + "bed", + ); + let panel = parse_snp_bed(path.to_str().unwrap()).unwrap(); + let _ = std::fs::remove_file(&path); + + assert_eq!(panel.num_sites, 3); + let chr1 = &panel.by_chrom["chr1"]; + assert_eq!(chr1.len(), 2); + assert_eq!(chr1[0].pos, 99, "BED end - 1 is the 0-based SNP position"); + assert_eq!(chr1[0].ref_base, b'A'); + assert_eq!(chr1[0].alt_base, b'G'); + assert_eq!(chr1[0].id, "rs1"); + } + + #[test] + fn test_parse_snp_bed_requires_alleles() { + let path = write_temp("chr1\t99\t100\trs1\n", "bed"); + let err = parse_snp_bed(path.to_str().unwrap()).unwrap_err(); + let _ = std::fs::remove_file(&path); + assert!( + err.to_string().contains("6-column NGSCheckMate layout"), + "unexpected error: {err}" + ); + } + + #[test] + fn test_allele_counting_and_genotype() { + // SNP at 0-based 100 (1-based 101), ref A alt G + let bed = write_temp("chr1\t100\t101\trs1\tA\tG\n", "bed"); + let panel = parse_snp_bed(bed.to_str().unwrap()).unwrap(); + let _ = std::fs::remove_file(&bed); + + // Reads start at 1-based 96, so the SNP is the 6th base of each read. + // r1/r2 carry G (alt), r3 carries A (ref), r4 is a duplicate and is + // ignored, r5 has MAPQ 5 and is ignored. + let sam = "\ +@HD\tVN:1.6\tSO:coordinate\n\ +@SQ\tSN:chr1\tLN:20000\n\ +r1\t0\tchr1\t96\t60\t10M\t*\t0\t0\tTTTTTGTTTT\tIIIIIIIIII\n\ +r2\t0\tchr1\t96\t60\t10M\t*\t0\t0\tTTTTTGTTTT\tIIIIIIIIII\n\ +r3\t0\tchr1\t96\t60\t10M\t*\t0\t0\tTTTTTATTTT\tIIIIIIIIII\n\ +r4\t1024\tchr1\t96\t60\t10M\t*\t0\t0\tTTTTTGTTTT\tIIIIIIIIII\n\ +r5\t0\tchr1\t96\t5\t10M\t*\t0\t0\tTTTTTGTTTT\tIIIIIIIIII\n"; + let sam_path = write_temp(sam, "sam"); + + let mut reader = bam::Reader::from_path(&sam_path).unwrap(); + let mut accum = SnpAccum::new(&panel, 30, 13); + let mut record = bam::Record::new(); + while let Some(res) = reader.read(&mut record) { + res.unwrap(); + accum.process_read(&record, "chr1", &panel); + } + let _ = std::fs::remove_file(&sam_path); + + let counts = accum.counts_for("chr1").unwrap(); + assert_eq!(counts[0].ref_count, 1, "one ref-supporting read"); + assert_eq!(counts[0].alt_count, 2, "two alt-supporting reads"); + assert_eq!(counts[0].other_count, 0); + assert_eq!(counts[0].depth(), 3); + assert!((counts[0].alt_fraction() - 2.0 / 3.0).abs() < 1e-9); + assert_eq!(counts[0].genotype(), "0/1"); + assert_eq!(accum.covered_sites(), 1); + } + + #[test] + fn test_genotype_thresholds() { + let hom_ref = AlleleCounts { + ref_count: 100, + alt_count: 5, + other_count: 0, + }; + assert_eq!(hom_ref.genotype(), "0/0"); + + let het = AlleleCounts { + ref_count: 50, + alt_count: 50, + other_count: 0, + }; + assert_eq!(het.genotype(), "0/1"); + + let hom_alt = AlleleCounts { + ref_count: 2, + alt_count: 100, + other_count: 0, + }; + assert_eq!(hom_alt.genotype(), "1/1"); + + assert_eq!(AlleleCounts::default().genotype(), "./."); + } + + #[test] + fn test_soft_clip_and_deletion_offsets() { + // Read at 1-based 100 (0-based 99) with CIGAR 5S5M2D5M and sequence + // CCCCC TTTTG AAAAA: + // 5S -> read bases 0..5, no reference + // 5M -> reference 99..103, read bases 5..10 ("TTTTG") + // 2D -> reference 104..105, no read bases + // 5M -> reference 106..110, read bases 10..15 ("AAAAA") + // rs1 at 0-based 103 must pick up the 'G' (alt) despite the soft clip; + // rs2 at 0-based 104 falls inside the deletion and must count nothing. + let bed = write_temp( + "chr1\t103\t104\trs1\tA\tG\n\ + chr1\t104\t105\trs2\tA\tG\n", + "bed", + ); + let panel = parse_snp_bed(bed.to_str().unwrap()).unwrap(); + let _ = std::fs::remove_file(&bed); + + // Reference positions 99..103 are the first 5M (read bases 5..9); + // position 104 is the 5th base of that block -> read base index 9 = 'G' + let sam = "\ +@HD\tVN:1.6\tSO:coordinate\n\ +@SQ\tSN:chr1\tLN:20000\n\ +r1\t0\tchr1\t100\t60\t5S5M2D5M\t*\t0\t0\tCCCCCTTTTGAAAAA\tIIIIIIIIIIIIIII\n"; + let sam_path = write_temp(sam, "sam"); + + let mut reader = bam::Reader::from_path(&sam_path).unwrap(); + let mut accum = SnpAccum::new(&panel, 0, 0); + let mut record = bam::Record::new(); + while let Some(res) = reader.read(&mut record) { + res.unwrap(); + accum.process_read(&record, "chr1", &panel); + } + let _ = std::fs::remove_file(&sam_path); + + let counts = accum.counts_for("chr1").unwrap(); + assert_eq!( + counts[0].alt_count, 1, + "soft clip must not shift the offset" + ); + assert_eq!(counts[0].ref_count, 0); + assert_eq!(counts[1].depth(), 0, "a site inside a deletion has no base"); + } +} diff --git a/src/cli.rs b/src/cli.rs index 6e6459e8..9617b98b 100644 --- a/src/cli.rs +++ b/src/cli.rs @@ -32,6 +32,130 @@ pub enum Commands { /// analyses in one pass. Requires a GTF annotation and duplicate-marked /// (not removed) alignments. Rna(RnaArgs), + + /// Single-pass alignment QC — samtools stats, mosdepth and NGSCheckMate + /// genotyping from one pass over a CRAM/BAM. + Align(AlignArgs), +} + +/// Arguments for the `align` subcommand. +#[derive(Parser, Debug)] +#[command( + next_line_help = false, + term_width = 120, + help_template = "\ +{about-with-newline} +{usage-heading} {usage} + +{all-args}" +)] +pub struct AlignArgs { + /// Coordinate-sorted alignment file (BAM/SAM/CRAM) + #[arg(value_name = "INPUT", required = true, help_heading = "Input / Output")] + pub input: String, + + /// Reference FASTA (required for CRAM) + #[arg( + short, + long, + value_name = "FASTA", + env = "RUSTQC_REFERENCE", + help_heading = "Input / Output" + )] + pub reference: Option, + + /// NGSCheckMate SNP BED (6 columns: chrom, start, end, id, ref, alt) + #[arg( + long = "snp-bed", + value_name = "BED", + env = "RUSTQC_SNP_BED", + help_heading = "Input / Output" + )] + pub snp_bed: Option, + + /// Output directory [default: .] + #[arg( + short, + long, + default_value = ".", + hide_default_value = true, + env = "RUSTQC_OUTDIR", + help_heading = "Input / Output" + )] + pub outdir: String, + + /// Override sample name (default: derived from the input filename) + #[arg( + long, + value_name = "NAME", + env = "RUSTQC_SAMPLE_NAME", + help_heading = "Input / Output" + )] + pub sample_name: Option, + + /// Depth window size in bases [default: 500] + #[arg( + long = "by", + value_name = "N", + default_value_t = 500, + hide_default_value = true, + env = "RUSTQC_DEPTH_WINDOW", + help_heading = "General" + )] + pub by: u64, + + /// MAPQ cutoff for genotyping and stats [default: 30] + #[arg( + short = 'Q', + long = "mapq", + default_value_t = 30, + hide_default_value = true, + env = "RUSTQC_MAPQ", + help_heading = "General" + )] + pub mapq_cut: u8, + + /// Minimum base quality for genotyping [default: 13] + #[arg( + long = "min-bq", + value_name = "N", + default_value_t = 13, + hide_default_value = true, + env = "RUSTQC_MIN_BQ", + help_heading = "General" + )] + pub min_base_quality: u8, + + /// Number of htslib decompression threads [default: 1] + #[arg( + short, + long, + default_value_t = 1, + hide_default_value = true, + env = "RUSTQC_THREADS", + help_heading = "General" + )] + pub threads: usize, + + /// Suppress output except warnings/errors + #[arg( + short = 'q', + long, + conflicts_with = "verbose", + env = "RUSTQC_QUIET", + help_heading = "General" + )] + pub quiet: bool, + + /// Show additional detail + #[arg( + short = 'v', + long, + conflicts_with = "quiet", + env = "RUSTQC_VERBOSE", + help_heading = "General" + )] + pub verbose: bool, } /// Arguments for the `rna` subcommand. diff --git a/src/lib.rs b/src/lib.rs index 9a228cae..74b139ae 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -65,6 +65,7 @@ use clap::ValueEnum; use serde::Deserialize; +pub mod align; pub mod config; pub mod cpu; pub mod gtf; diff --git a/src/main.rs b/src/main.rs index 4c66c173..46b075f3 100644 --- a/src/main.rs +++ b/src/main.rs @@ -76,6 +76,8 @@ fn main() -> Result<()> { let verbosity = match &cli.command { cli::Commands::Rna(args) if args.quiet => Verbosity::Quiet, cli::Commands::Rna(args) if args.verbose => Verbosity::Verbose, + cli::Commands::Align(args) if args.quiet => Verbosity::Quiet, + cli::Commands::Align(args) if args.verbose => Verbosity::Verbose, _ => Verbosity::Normal, }; @@ -94,9 +96,156 @@ fn main() -> Result<()> { match cli.command { cli::Commands::Rna(args) => run_rna(args, &ui), + cli::Commands::Align(args) => run_align(args, &ui), } } +/// Run the `align` subcommand: samtools stats, mosdepth-equivalent depth and +/// NGSCheckMate genotyping from a single streaming pass. +/// +/// # Arguments +/// * `args` - Parsed CLI arguments +/// * `ui` - Terminal UI handle +fn run_align(args: cli::AlignArgs, ui: &Ui) -> Result<()> { + use rustqc::align; + + let start = Instant::now(); + let sample_name = args.sample_name.clone().unwrap_or_else(|| { + Path::new(&args.input) + .file_stem() + .and_then(|s| s.to_str()) + .unwrap_or("sample") + .to_string() + }); + let outdir = Path::new(&args.outdir); + std::fs::create_dir_all(outdir) + .with_context(|| format!("Failed to create output directory: {}", outdir.display()))?; + + // SNP panel for NGSCheckMate genotyping (optional) + let panel = match args.snp_bed.as_deref() { + Some(path) => { + let panel = align::parse_snp_bed(path)?; + ui.detail(&format!( + "Loaded {} SNP sites for genotyping", + format_count(panel.num_sites as u64) + )); + Some(panel) + } + None => None, + }; + + ui.blank(); + ui.step(&format!("Processing {}", args.input)); + + let mut reader = rust_htslib::bam::Reader::from_path(&args.input) + .with_context(|| format!("Failed to open alignment file: {}", args.input))?; + if let Some(ref_path) = args.reference.as_deref() { + reader + .set_reference(ref_path) + .with_context(|| format!("Failed to set reference FASTA: {}", ref_path))?; + } + if args.threads > 1 { + reader + .set_threads(args.threads.saturating_sub(1)) + .context("Failed to set htslib decompression threads")?; + } + + let header = reader.header().clone(); + let contigs: Vec<(String, u64)> = (0..header.target_count()) + .map(|tid| { + ( + String::from_utf8_lossy(header.tid2name(tid)).to_string(), + header.target_len(tid).unwrap_or(0), + ) + }) + .collect(); + ensure!( + !contigs.is_empty(), + "Alignment file has no reference sequences in its header: {}", + args.input + ); + + let mut stats_accum = rna::rseqc::accumulators::BamStatAccum::default(); + let mut depth_accum = align::DepthAccum::new(contigs.clone(), args.by); + let mut snp_accum = panel + .as_ref() + .map(|p| align::SnpAccum::new(p, args.mapq_cut, args.min_base_quality)); + + let mut record = rust_htslib::bam::Record::new(); + let mut n = 0u64; + while let Some(res) = reader.read(&mut record) { + res.context("Error reading alignment record")?; + n += 1; + + stats_accum.process_read(&record, args.mapq_cut); + depth_accum.process_read(&record); + + if let (Some(accum), Some(panel)) = (&mut snp_accum, panel.as_ref()) { + let tid = record.tid(); + if tid >= 0 && (tid as usize) < contigs.len() { + accum.process_read(&record, &contigs[tid as usize].0, panel); + } + } + } + + let stats_result = stats_accum.into_result(); + let (contig_depth, windows) = depth_accum.finish(); + + // --- samtools-compatible outputs --- + let stats_path = outdir.join(format!("{sample_name}.stats")); + rna::rseqc::stats::write_stats(&stats_result, &stats_path)?; + let flagstat_path = outdir.join(format!("{sample_name}.flagstat")); + rna::rseqc::flagstat::write_flagstat(&stats_result, &flagstat_path)?; + let bam_header_refs: Vec<(String, u64)> = contigs.clone(); + let idxstats_path = outdir.join(format!("{sample_name}.idxstats")); + rna::rseqc::idxstats::write_idxstats(&stats_result, &bam_header_refs, &idxstats_path)?; + + // --- mosdepth-compatible outputs --- + let summary_path = outdir.join(format!("{sample_name}.mosdepth.summary.txt")); + align::output::write_summary(&contig_depth, &summary_path)?; + let dist_path = outdir.join(format!("{sample_name}.mosdepth.global.dist.txt")); + align::output::write_global_dist(&contig_depth, &dist_path)?; + let regions_path = outdir.join(format!("{sample_name}.regions.bed.gz")); + align::output::write_regions(&windows, ®ions_path)?; + + ui.output_group("samtools"); + for path in [&stats_path, &flagstat_path, &idxstats_path] { + ui.output_item("samtools", &path.display().to_string()); + } + ui.output_group("mosdepth"); + for path in [&summary_path, &dist_path, ®ions_path] { + ui.output_item("mosdepth", &path.display().to_string()); + } + + // --- NGSCheckMate VCF --- + if let (Some(panel), Some(accum)) = (panel.as_ref(), snp_accum.as_ref()) { + let vcf_path = outdir.join(format!("{sample_name}.ngscheckmate.vcf.gz")); + align::output::write_ngscheckmate_vcf(panel, accum, &sample_name, &contigs, &vcf_path)?; + ui.output_group("ngscheckmate"); + ui.output_item("ngscheckmate", &vcf_path.display().to_string()); + ui.output_detail(&format!( + "{} of {} SNP sites covered", + format_count(accum.covered_sites() as u64), + format_count(panel.num_sites as u64), + )); + } + + let total_length: u64 = contigs.iter().map(|(_, len)| len).sum(); + let total_bases: u64 = contig_depth.iter().map(|c| c.total_bases).sum(); + ui.blank(); + ui.detail(&format!( + "{} records, mean depth {:.2}X, finished in {}", + format_count(n), + if total_length == 0 { + 0.0 + } else { + total_bases as f64 / total_length as f64 + }, + format_duration(start.elapsed()) + )); + Ok(()) +} + /// Reconstruct the command line for the featureCounts-compatible header comment. fn reconstruct_command_line(args: &cli::RnaArgs) -> String { let mut parts = vec![format!( diff --git a/tests/integration_test.rs b/tests/integration_test.rs index c7ea3acd..15091a08 100644 --- a/tests/integration_test.rs +++ b/tests/integration_test.rs @@ -1026,3 +1026,126 @@ fn test_dup_check_parallel_uses_global_duplicate_state() { let _ = fs::remove_dir_all(root); } + +// ============================================================================ +// align subcommand (samtools stats + mosdepth + NGSCheckMate) +// ============================================================================ + +#[test] +fn test_align_outputs_match_mosdepth_reference() { + let root = unique_test_dir("align"); + let outdir = root.display().to_string(); + + // Three SNP sites: one covered by a single read, two uncovered. + let snp_bed = root.join("snps.bed"); + fs::write( + &snp_bed, + "chr1\t1014\t1015\trs_test1\tA\tG\n\ + chr1\t2999\t3000\trs_test2\tC\tT\n\ + chr2\t999\t1000\trs_test3\tG\tA\n", + ) + .unwrap(); + + let binary = rustqc_binary(); + let output = Command::new(&binary) + .args([ + "align", + "tests/data/test.bam", + "--snp-bed", + snp_bed.to_str().unwrap(), + "--by", + "500", + "--outdir", + &outdir, + ]) + .output() + .expect("Failed to run rustqc align"); + assert!( + output.status.success(), + "rustqc align failed:\n{}", + String::from_utf8_lossy(&output.stderr) + ); + + // --- mosdepth summary --- + // Reference values from `mosdepth --by 500 md tests/data/test.bam` + // (mosdepth 0.3.x), byte-for-byte. + let summary = fs::read_to_string(format!("{outdir}/test.mosdepth.summary.txt")).unwrap(); + assert_eq!( + summary, + "chrom\tlength\tbases\tmean\tmin\tmax\n\ + chr1\t20000\t10300\t0.52\t0\t12\n\ + chr1_region\t20000\t10300\t0.52\t0\t12\n\ + chr2\t20000\t7000\t0.35\t0\t40\n\ + chr2_region\t20000\t7000\t0.35\t0\t40\n\ + total\t40000\t17300\t0.43\t0\t40\n\ + total_region\t40000\t17300\t0.43\t0\t40\n", + "summary must match mosdepth" + ); + + // --- mosdepth global distribution --- + let dist = fs::read_to_string(format!("{outdir}/test.mosdepth.global.dist.txt")).unwrap(); + // mosdepth skips the sparse tail (cumulative < 8e-5), so chr2's deepest + // reported level is 39 even though max depth is 40. + assert!( + dist.starts_with("chr1\t12\t0.00\n"), + "distribution starts at the deepest reported level" + ); + assert!( + dist.contains("\nchr2\t39\t0.00\n"), + "chr2 tail is skipped at 40" + ); + assert!( + !dist.contains("chr2\t40\t"), + "depth 40 is below mosdepth's cutoff" + ); + assert!(dist.trim_end().ends_with("total\t0\t1.00")); + + // --- samtools-compatible outputs exist --- + for file in ["test.stats", "test.flagstat", "test.idxstats"] { + assert!( + Path::new(&format!("{outdir}/{file}")).exists(), + "missing output: {file}" + ); + } + let flagstat = fs::read_to_string(format!("{outdir}/test.flagstat")).unwrap(); + assert!( + flagstat.starts_with("488 + 0 in total"), + "flagstat total should be 488, got:\n{flagstat}" + ); + + // --- regions and NGSCheckMate VCF are BGZF files that exist --- + assert!(Path::new(&format!("{outdir}/test.regions.bed.gz")).exists()); + assert!(Path::new(&format!("{outdir}/test.ngscheckmate.vcf.gz")).exists()); + + let _ = fs::remove_dir_all(&root); +} + +#[test] +fn test_align_rejects_snp_bed_without_alleles() { + let root = unique_test_dir("align-badbed"); + let outdir = root.display().to_string(); + let snp_bed = root.join("bad.bed"); + fs::write(&snp_bed, "chr1\t1014\t1015\trs_test1\n").unwrap(); + + let binary = rustqc_binary(); + let output = Command::new(&binary) + .args([ + "align", + "tests/data/test.bam", + "--snp-bed", + snp_bed.to_str().unwrap(), + "--outdir", + &outdir, + ]) + .output() + .expect("Failed to run rustqc align"); + + assert!(!output.status.success(), "should reject a 4-column SNP BED"); + let stderr = String::from_utf8_lossy(&output.stderr); + assert!( + stderr.contains("6-column NGSCheckMate layout"), + "error should name the expected layout:\n{stderr}" + ); + + let _ = fs::remove_dir_all(&root); +}