diff --git a/AGENTS.md b/AGENTS.md index 108414ac..c29f323a 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -98,6 +98,7 @@ src/ plots.rs — RSeQC plot generation (duplication, junctions, etc.) read_distribution.rs — read_distribution.py reimplementation read_duplication.rs — read_duplication.py reimplementation + split_bam.rs — split_bam.py-style BED-interval classification (rRNA quantification) stats.rs — samtools stats full output (SN + all histogram sections) tin.rs — TIN (Transcript Integrity Number) analysis tests/ @@ -294,6 +295,15 @@ forwarded to `count_reads()` as the `skip_dup_check: bool` parameter). `infer_experiment:`, `read_duplication:`, `read_distribution:`, `junction_annotation:`, `junction_saturation:`, `inner_distance:`, `tin:`). Each has an `enabled: bool` toggle and tool-specific parameter overrides. CLI flags take precedence over config values. +- `split_bam` (`src/rna/rseqc/split_bam.rs`) classifies every alignment against a BED + file of intervals (`--rrna-bed` or `rna.split_bam.bed`) into in/ex/junk counts, the + counting side of RSeQC's `split_bam.py`. It is disabled unless a BED file is given, + and deliberately does not filter secondary/supplementary alignments so rDNA + multi-mappers are visible. No split BAM files are written. +- featureCounts `-M` / `-O` equivalents are exposed as `--count-multi-mapping` / + `--count-multi-overlapping` (and `rna.featurecounts.count_multi_mapping` / + `count_multi_overlapping`). They affect only the featureCounts gene-level and + biotype-level counts, not the dupRadar matrix. - Under `rna:`, there are also sections for `preseq:`, `qualimap:`, `flagstat:`, `idxstats:`, and `samtools_stats:`. Each has an `enabled: bool` toggle. Preseq has additional parameters: `max_extrap`, `step_size`, `n_bootstraps`, diff --git a/docs/src/content/docs/rna/featurecounts.mdx b/docs/src/content/docs/rna/featurecounts.mdx index a74a21a5..40536e56 100644 --- a/docs/src/content/docs/rna/featurecounts.mdx +++ b/docs/src/content/docs/rna/featurecounts.mdx @@ -115,6 +115,47 @@ The biotype is extracted from the GTF attribute specified by `biotype_attribute` - `.biotype_counts_mqc.tsv` -- Biotype counts formatted as a MultiQC bargraph data file, suitable for visualizing the distribution of reads across biotypes. - `.biotype_counts_rrna_mqc.tsv` -- rRNA percentage formatted as a MultiQC general statistics value, reporting the fraction of assigned reads mapping to rRNA genes. +## Counting multi-mapping and multi-overlapping reads + +By default RustQC follows featureCounts' defaults: a read with `NH` > 1 is +reported as `Unassigned_MultiMapping`, and a read overlapping more than one +feature is reported as `Unassigned_Ambiguity`. Two flags relax those rules: + +| Flag | featureCounts equivalent | Effect | +| --------------------------- | ------------------------ | ---------------------------------------------------------------- | +| `-M`, `--count-multi-mapping` | `-M` | Count every reported alignment of a multi-mapping read | +| `-O`, `--count-multi-overlapping` | `-O` | Count a read once for each feature it overlaps | + +Both apply to the gene-level and biotype-level featureCounts outputs. dupRadar's +duplicate-rate matrix is unaffected: it tracks multi-mappers separately by design. + +In the YAML config: + +```yaml +rna: + featurecounts: + count_multi_mapping: true + count_multi_overlapping: true +``` + + + ## Biotype attribute detection RustQC auto-detects whether your GTF uses `gene_biotype` (Ensembl) or `gene_type` (GENCODE). You can override this with `--biotype-attribute`. See the [CLI reference](../usage/cli-reference/) for details. diff --git a/docs/src/content/docs/rna/rseqc.mdx b/docs/src/content/docs/rna/rseqc.mdx index 0cf64346..ef0e7cfb 100644 --- a/docs/src/content/docs/rna/rseqc.mdx +++ b/docs/src/content/docs/rna/rseqc.mdx @@ -609,6 +609,65 @@ RustQC produces gene-level TIN scores (one row per gene, using the longest transcript as representative), while RSeQC's `tin.py` produces transcript-level scores. Both formats are compatible with MultiQC. +## split_bam (BED-interval classification) + +Reimplements the counting side of RSeQC's `split_bam.py`: every alignment is +classified against a BED file of genomic intervals. The usual use is **rRNA +quantification**, where interval overlap is far more robust than the GTF/biotype +route (see the [featureCounts page](../featurecounts/#counting-multi-mapping-and-multi-overlapping-reads) +for why biotype %rRNA under-reports on stock GRCh38/GENCODE builds). + +Disabled unless a BED file is supplied: + +```bash +rustqc rna sample.bam --gtf genes.gtf --rrna-bed GRCh38_rRNA.bed +``` + +| File | Description | +| ------------------------- | ------------------------------------------------------ | +| `{stem}.split_bam.tsv` | Per-category counts and percentages | + +Categories follow `split_bam.py`: + +| Category | Meaning | +| -------- | ----------------------------------------------------------- | +| `in` | Alignment (or its mate) starts inside a BED interval | +| `ex` | Mapped, QC-passing alignment that does not | +| `junk` | Unmapped or QC-failed record | + +``` +category count percent_of_total percent_of_usable +in 107 21.9262 22.1532 +ex 376 77.0492 77.8468 +junk 5 1.0246 NA +total 488 100.0000 NA +``` + +`percent_of_usable` excludes `junk` records from the denominator, so it answers +"what fraction of usable alignments fall inside the intervals?". + + + + + +In the YAML config: + +```yaml +rna: + split_bam: + enabled: true + bed: /refs/GRCh38_rRNA.bed +``` + ## Compatibility with RSeQC All output files are designed to be drop-in replacements for the corresponding diff --git a/docs/src/content/docs/usage/cli-reference.mdx b/docs/src/content/docs/usage/cli-reference.mdx index deca6b40..80191d27 100644 --- a/docs/src/content/docs/usage/cli-reference.mdx +++ b/docs/src/content/docs/usage/cli-reference.mdx @@ -134,6 +134,26 @@ If not specified, RustQC defaults to `gene_biotype` and auto-detects the attribute. If the specified attribute is not found in the GTF, a warning is printed and biotype counting is skipped. +#### `--rrna-bed ` + +BED file of genomic intervals (plain or gzip-compressed) used by the +[`split_bam`](../../rna/rseqc/#split_bam-bed-interval-classification) analysis. +Supplying this flag enables the analysis, which reports how many alignments fall +inside the intervals. The usual use is rRNA quantification, which is more +reliable than the biotype route on stock GRCh38/GENCODE builds. + +#### `-M, --count-multi-mapping` + +Count multi-mapping reads (`NH` > 1) in the featureCounts gene-level and +biotype-level outputs, one count per reported alignment. Equivalent to +`featureCounts -M`. By default such reads are reported as +`Unassigned_MultiMapping`. + +#### `-O, --count-multi-overlapping` + +Count a read once for every feature it overlaps instead of reporting it as +`Unassigned_Ambiguity`. Equivalent to `featureCounts -O`. + #### `--sample-name ` Override the sample name used for output filenames. By default, the sample name diff --git a/docs/src/content/docs/usage/configuration.md b/docs/src/content/docs/usage/configuration.md index a6c21d7d..507465b4 100644 --- a/docs/src/content/docs/usage/configuration.md +++ b/docs/src/content/docs/usage/configuration.md @@ -131,6 +131,9 @@ Run `rustqc rna --help` to see the associated environment variable for each flag | `RUSTQC_THREADS` | `--threads` | Number of threads | | `RUSTQC_MAPQ` | `--mapq` | MAPQ quality cutoff | | `RUSTQC_BIOTYPE_ATTRIBUTE` | `--biotype-attribute` | GTF biotype attribute name | +| `RUSTQC_RRNA_BED` | `--rrna-bed` | BED intervals for split_bam classification | +| `RUSTQC_COUNT_MULTI_MAPPING` | `--count-multi-mapping` | Count multi-mapping reads (featureCounts `-M`) | +| `RUSTQC_COUNT_MULTI_OVERLAPPING` | `--count-multi-overlapping` | Count multi-overlapping reads (featureCounts `-O`) | | `RUSTQC_SKIP_DUP_CHECK` | `--skip-dup-check` | Skip duplicate-marking check | | `RUSTQC_QUIET` | `--quiet` | Suppress output | | `RUSTQC_VERBOSE` | `--verbose` | Show additional detail | @@ -389,8 +392,24 @@ rna: biotype_counts_mqc: true # Biotype counts MultiQC bargraph file biotype_rrna_mqc: true # Biotype rRNA percentage MultiQC file biotype_attribute: "gene_biotype" # GTF attribute for biotype grouping + count_multi_mapping: false # featureCounts -M: count multi-mapping reads + count_multi_overlapping: false # featureCounts -O: count reads once per overlapping feature ``` +### `count_multi_mapping` / `count_multi_overlapping` + +By default a read with `NH` > 1 is reported as `Unassigned_MultiMapping`, and a +read overlapping several features as `Unassigned_Ambiguity`. These toggles are +the featureCounts `-M` and `-O` equivalents and also have CLI flags (`-M` / +`--count-multi-mapping`, `-O` / `--count-multi-overlapping`), which take +precedence over the config file. + +They matter for high-copy repeat families such as rDNA, where most reads +multi-map and are otherwise discarded — see +[split_bam](#split_bam) for the interval-based alternative. + +**Default:** `false` for both. + ### `biotype_attribute` The GTF attribute name used for biotype grouping. This controls how genes are @@ -537,6 +556,23 @@ read origin classification (exonic/intronic/intergenic), strand-specificity estimation, and splice junction motif counting. Produces Qualimap-compatible output files parseable by MultiQC. +## split_bam + +```yaml +rna: + split_bam: + enabled: true + bed: /refs/GRCh38_rRNA.bed # BED intervals (plain or .gz) +``` + +Classifies every alignment as inside (`in`) or outside (`ex`) the BED intervals, +with unmapped/QC-failed records reported as `junk` — the counting side of +RSeQC's `split_bam.py`. Typically used for rRNA quantification, which is more +robust than biotype counting on stock GRCh38/GENCODE builds. + +Disabled unless a BED file is supplied. The `--rrna-bed` CLI flag takes +precedence over the config file value and enables the tool on its own. + ## preseq ```yaml diff --git a/src/cli.rs b/src/cli.rs index 6e6459e8..37a0fac2 100644 --- a/src/cli.rs +++ b/src/cli.rs @@ -155,6 +155,33 @@ pub struct RnaArgs { )] pub biotype_attribute: Option, + /// BED file of intervals (e.g. rRNA regions) for split_bam classification + #[arg( + long = "rrna-bed", + value_name = "BED", + env = "RUSTQC_RRNA_BED", + help_heading = "General" + )] + pub rrna_bed: Option, + + /// Count multi-mapping reads (featureCounts -M) + #[arg( + short = 'M', + long = "count-multi-mapping", + env = "RUSTQC_COUNT_MULTI_MAPPING", + help_heading = "General" + )] + pub count_multi_mapping: bool, + + /// Count reads overlapping several features (featureCounts -O) + #[arg( + short = 'O', + long = "count-multi-overlapping", + env = "RUSTQC_COUNT_MULTI_OVERLAPPING", + help_heading = "General" + )] + pub count_multi_overlapping: bool, + /// Skip duplicate-marking check #[arg( long, diff --git a/src/config.rs b/src/config.rs index 4952a1fd..9ed29a12 100644 --- a/src/config.rs +++ b/src/config.rs @@ -165,6 +165,10 @@ pub struct RnaConfig { /// Qualimap RNA-Seq QC configuration. #[serde(default)] pub qualimap: QualimapConfig, + + /// split_bam BED-interval classification configuration (rRNA quantification). + #[serde(default)] + pub split_bam: SplitBamConfig, } // ============================================================================ @@ -275,6 +279,24 @@ pub struct FeatureCountsConfig { /// Defaults to `"gene_biotype"` (Ensembl convention). /// Use `"gene_type"` for GENCODE GTF files. pub biotype_attribute: String, + + /// Count multi-mapping reads (featureCounts `-M`). + /// + /// By default multi-mapping reads (`NH` > 1) are reported as + /// `Unassigned_MultiMapping` and never counted. When enabled, every + /// reported alignment of a multi-mapping read is counted, as + /// `featureCounts -M` does. Matters for high-copy repeat families such as + /// rDNA, where most reads multi-map. + /// **Default:** `false`. + pub count_multi_mapping: bool, + + /// Count reads overlapping several features (featureCounts `-O`). + /// + /// By default a read overlapping more than one gene is reported as + /// `Unassigned_Ambiguity`. When enabled, the read is counted once for + /// every feature it overlaps, as `featureCounts -O` does. + /// **Default:** `false`. + pub count_multi_overlapping: bool, } impl Default for FeatureCountsConfig { @@ -287,6 +309,8 @@ impl Default for FeatureCountsConfig { biotype_counts_mqc: true, biotype_rrna_mqc: true, biotype_attribute: "gene_biotype".to_string(), + count_multi_mapping: false, + count_multi_overlapping: false, } } } @@ -600,6 +624,33 @@ impl Default for QualimapConfig { } } +/// Configuration for split_bam BED-interval read classification. +/// +/// Reimplements the counting side of RSeQC's `split_bam.py`: every alignment +/// is classified as overlapping (`in`) or not overlapping (`ex`) a set of BED +/// intervals, with unmapped/QC-failed records reported as `junk`. The usual +/// use is rRNA quantification, which the GTF/biotype route under-reports on +/// stock GRCh38/GENCODE builds. +/// +/// Disabled unless a BED file is supplied (`--rrna-bed` or `bed:` here). +/// +/// Example: +/// ```yaml +/// split_bam: +/// enabled: true +/// bed: /refs/GRCh38_rRNA.bed +/// ``` +#[derive(Debug, Deserialize, Default)] +#[serde(default)] +pub struct SplitBamConfig { + /// Whether to run the BED-interval classification. Requires `bed`. + pub enabled: bool, + /// Path to the BED file of intervals (plain or gzip-compressed). + /// + /// The CLI `--rrna-bed` flag takes precedence over this setting. + pub bed: Option, +} + /// Configuration for samtools idxstats-compatible output. /// /// When enabled, produces a file matching `samtools idxstats` output format, diff --git a/src/main.rs b/src/main.rs index 4c66c173..3362980d 100644 --- a/src/main.rs +++ b/src/main.rs @@ -471,6 +471,32 @@ fn run_rna(args: cli::RnaArgs, ui: &Ui) -> Result<()> { None }; + // Load BED intervals for split_bam classification (rRNA quantification). + // CLI --rrna-bed takes precedence over the config file's split_bam.bed. + let bed_path: Option = args + .rrna_bed + .clone() + .or_else(|| config.split_bam.bed.clone()) + .filter(|_| args.rrna_bed.is_some() || config.split_bam.enabled); + let bed_intervals = match bed_path.as_deref() { + Some(path) => { + ui.detail(&format!("Loading BED intervals from {}...", path)); + let intervals = rna::rseqc::split_bam::parse_bed(path)?; + ui.detail(&format!( + "Loaded {} intervals for split_bam classification", + format_count(intervals.num_intervals as u64) + )); + Some(intervals) + } + None => None, + }; + + let fc_options = rna::dupradar::counting::FeatureCountsOptions { + count_multi_mapping: args.count_multi_mapping || config.featurecounts.count_multi_mapping, + count_multi_overlapping: args.count_multi_overlapping + || config.featurecounts.count_multi_overlapping, + }; + let chrom_mapping = config.alignment_to_gtf_mapping(); let chrom_prefix = config.chromosome_prefix().map(|s| s.to_owned()); @@ -570,6 +596,9 @@ fn run_rna(args: cli::RnaArgs, ui: &Ui) -> Result<()> { tin_min_coverage: config.tin.min_coverage.unwrap_or(10), gtf_path: &args.gtf, sample_name_override: effective_sample_name, + bed_intervals: bed_intervals.as_ref(), + bed_path: bed_path.as_deref(), + fc_options, }; // Step 2: Process all alignment files (in parallel when multiple) @@ -815,6 +844,12 @@ struct SharedParams<'a> { tin_min_coverage: u32, /// Path to GTF file (for Qualimap report output). gtf_path: &'a str, + /// Pre-parsed BED intervals for split_bam classification (from --rrna-bed). + bed_intervals: Option<&'a rna::rseqc::split_bam::BedIntervals>, + /// BED file path the intervals came from (recorded in the output header). + bed_path: Option<&'a str>, + /// featureCounts counting-mode options (-M / -O). + fc_options: rna::dupradar::counting::FeatureCountsOptions, } // ============================================================================ @@ -988,6 +1023,7 @@ fn process_single_bam( junction_saturation_seed: config.junction_saturation.seed.unwrap_or(42), preseq_enabled: config.preseq.enabled, preseq_max_segment_length: config.preseq.max_segment_length, + split_bam_enabled: params.bed_intervals.is_some(), }; let rseqc_annotations = RseqcAnnotations { @@ -997,6 +1033,7 @@ fn process_single_bam( exon_bitset: params.exon_bitset, transcript_tree: params.transcript_tree, tin_index: params.tin_index, + bed_intervals: params.bed_intervals, }; let any_rseqc_enabled = rseqc_config.bam_stat_enabled @@ -1007,7 +1044,8 @@ fn process_single_bam( || rseqc_config.junction_saturation_enabled || rseqc_config.inner_distance_enabled || rseqc_config.preseq_enabled - || rseqc_config.tin_enabled; + || rseqc_config.tin_enabled + || rseqc_config.split_bam_enabled; // === Build Qualimap exon index (if enabled) === let qualimap_index = if params.config.qualimap.enabled { @@ -1043,6 +1081,7 @@ fn process_single_bam( None }, qualimap_index.as_ref(), + params.fc_options, Some(&pb), )?; let count_duration = count_start.elapsed(); @@ -1871,6 +1910,26 @@ fn write_rseqc_outputs( } } + // --- split_bam (BED-interval classification, e.g. rRNA) --- + if let (Some(accum), Some(bed_path)) = (accums.split_bam, params.bed_path) { + let split_dir = if flat { + outdir.to_path_buf() + } else { + outdir.join("rseqc").join("split_bam") + }; + std::fs::create_dir_all(&split_dir)?; + let output_path = split_dir.join(format!("{}.split_bam.tsv", sample_name)); + let result = accum.into_result(); + rna::rseqc::split_bam::write_split_bam_summary(&result, bed_path, &output_path)?; + let p = output_path.display().to_string(); + ui.output_item("split_bam", &p); + ui.output_detail(&format!( + "{:.2}% of usable alignments inside intervals", + result.percent_in() + )); + written.push(("split_bam".into(), p)); + } + Ok(RseqcOutputs { written, infer_experiment: infer_experiment_result, diff --git a/src/rna/dupradar/counting.rs b/src/rna/dupradar/counting.rs index 785115c5..45d06e0b 100644 --- a/src/rna/dupradar/counting.rs +++ b/src/rna/dupradar/counting.rs @@ -301,6 +301,20 @@ pub struct CountResult { pub qualimap: Option, } +/// featureCounts counting-mode options (`-M` / `-O` equivalents). +/// +/// These affect only the featureCounts-compatible gene-level and biotype-level +/// counts. dupRadar's own duplicate-rate matrix keeps its established +/// semantics, where multi-mappers are tracked separately by design. +#[derive(Debug, Clone, Copy, Default)] +pub struct FeatureCountsOptions { + /// Count multi-mapping reads (`NH` > 1), one count per reported alignment. + pub count_multi_mapping: bool, + /// Count a read once for every feature it overlaps instead of calling it + /// ambiguous. + pub count_multi_overlapping: bool, +} + /// Metadata stored with each interval in the cache-oblivious interval tree. #[derive(Debug, Clone, Copy, Default)] struct IvMeta { @@ -677,8 +691,13 @@ impl ChromResult { /// When multiple genes are hit, the read is Unassigned_Ambiguity /// (matching default featureCounts `-g gene_id` behaviour where each /// gene is its own meta-feature). -fn classify_read_fc(is_multi: bool, gene_hits: &[GeneIdx], result: &mut ChromResult) { - if is_multi { +fn classify_read_fc( + is_multi: bool, + gene_hits: &[GeneIdx], + options: FeatureCountsOptions, + result: &mut ChromResult, +) { + if is_multi && !options.count_multi_mapping { result.fc_multimapping += 1; } else if gene_hits.is_empty() { result.fc_no_features += 1; @@ -688,6 +707,16 @@ fn classify_read_fc(is_multi: bool, gene_hits: &[GeneIdx], result: &mut ChromRes if idx < result.gene_counts.len() { result.gene_counts[idx].fc_reads += 1; } + } else if options.count_multi_overlapping { + // featureCounts -O: the read is Assigned once and counted for every + // feature it overlaps. + result.fc_assigned += 1; + for &gidx in gene_hits { + let idx = gidx as usize; + if idx < result.gene_counts.len() { + result.gene_counts[idx].fc_reads += 1; + } + } } else { // Multiple gene hits → Ambiguous (default featureCounts behaviour) result.fc_ambiguous += 1; @@ -714,10 +743,12 @@ fn classify_read_fc_biotype( gene_hits: &[GeneIdx], gene_to_biotype: &[u16], biotype_hits_buf: &mut Vec, + options: FeatureCountsOptions, result: &mut ChromResult, ) { - // Multi-mapped reads are excluded from biotype counting (same as gene-level) - if is_multi { + // Multi-mapped reads are excluded from biotype counting (same as + // gene-level) unless -M is in effect. + if is_multi && !options.count_multi_mapping { return; } if gene_hits.is_empty() { @@ -756,6 +787,16 @@ fn classify_read_fc_biotype( // meta-feature), but not tracked in named biotype counts result.fc_biotype_assigned += 1; } + } else if options.count_multi_overlapping { + // featureCounts -O at the biotype level: count the read once for each + // distinct known biotype meta-feature it overlaps. + result.fc_biotype_assigned += 1; + for &bidx in biotype_hits_buf.iter() { + let idx = bidx as usize; + if idx < result.biotype_reads.len() { + result.biotype_reads[idx] += 1; + } + } } else { // Multiple distinct meta-features → Ambiguous at biotype level result.fc_biotype_ambiguous += 1; @@ -781,6 +822,7 @@ fn process_counting_record( stranded: Strandedness, paired: bool, gene_to_biotype: &[u16], + fc_options: FeatureCountsOptions, aligned_blocks_buf: &mut Vec<(u64, u64)>, gene_hits: &mut Vec, biotype_hits_buf: &mut Vec, @@ -847,12 +889,13 @@ fn process_counting_record( } // --- Per-read featureCounts counting (independent of mate pairing) --- - classify_read_fc(is_multi, gene_hits, result); + classify_read_fc(is_multi, gene_hits, fc_options, result); classify_read_fc_biotype( is_multi, gene_hits, gene_to_biotype, biotype_hits_buf, + fc_options, result, ); @@ -981,6 +1024,7 @@ fn process_chromosome_batch( qualimap_index: Option<&crate::rna::qualimap::QualimapIndex>, gene_to_biotype: &[u16], num_biotypes: usize, + fc_options: FeatureCountsOptions, progress: Option<&ProgressBar>, ) -> Result<(ChromResult, Option)> { let mut result = ChromResult::new(num_genes, num_biotypes); @@ -1090,6 +1134,7 @@ fn process_chromosome_batch( stranded, paired, gene_to_biotype, + fc_options, &mut aligned_blocks_buf, &mut gene_hits, &mut biotype_hits_buf, @@ -1152,6 +1197,7 @@ pub fn count_reads( rseqc_config: Option<&RseqcConfig>, rseqc_annotations: Option<&RseqcAnnotations>, qualimap_index: Option<&crate::rna::qualimap::QualimapIndex>, + fc_options: FeatureCountsOptions, progress: Option<&ProgressBar>, ) -> Result { // Build gene ID interner for allocation-free lookups in the hot loop @@ -1302,6 +1348,7 @@ pub fn count_reads( qualimap_index, &gene_to_biotype, num_biotypes, + fc_options, progress, ) }) @@ -1433,6 +1480,7 @@ pub fn count_reads( stranded, paired, &gene_to_biotype, + fc_options, &mut aligned_blocks_buf, &mut gene_hits, &mut biotype_hits_buf, @@ -1783,7 +1831,12 @@ mod tests { let mut result = make_test_chrom_result(3); let gene_hits: Vec = vec![0]; // has a hit, but is_multi=true - classify_read_fc(true, &gene_hits, &mut result); + classify_read_fc( + true, + &gene_hits, + FeatureCountsOptions::default(), + &mut result, + ); assert_eq!( result.fc_multimapping, 1, @@ -1803,7 +1856,12 @@ mod tests { let mut result = make_test_chrom_result(3); let gene_hits: Vec = vec![]; - classify_read_fc(false, &gene_hits, &mut result); + classify_read_fc( + false, + &gene_hits, + FeatureCountsOptions::default(), + &mut result, + ); assert_eq!( result.fc_no_features, 1, @@ -1819,7 +1877,12 @@ mod tests { let mut result = make_test_chrom_result(3); let gene_hits: Vec = vec![1]; - classify_read_fc(false, &gene_hits, &mut result); + classify_read_fc( + false, + &gene_hits, + FeatureCountsOptions::default(), + &mut result, + ); assert_eq!( result.fc_assigned, 1, @@ -1841,7 +1904,12 @@ mod tests { let mut result = make_test_chrom_result(4); let gene_hits: Vec = vec![0, 1]; - classify_read_fc(false, &gene_hits, &mut result); + classify_read_fc( + false, + &gene_hits, + FeatureCountsOptions::default(), + &mut result, + ); assert_eq!( result.fc_ambiguous, 1, @@ -1861,7 +1929,12 @@ mod tests { let mut result = make_test_chrom_result(5); let gene_hits: Vec = vec![0, 1, 2]; - classify_read_fc(false, &gene_hits, &mut result); + classify_read_fc( + false, + &gene_hits, + FeatureCountsOptions::default(), + &mut result, + ); assert_eq!( result.fc_ambiguous, 1, @@ -1879,15 +1952,15 @@ mod tests { let mut result = make_test_chrom_result(3); // First call: single hit to gene 0 - classify_read_fc(false, &[0], &mut result); + classify_read_fc(false, &[0], FeatureCountsOptions::default(), &mut result); // Second call: single hit to gene 0 again - classify_read_fc(false, &[0], &mut result); + classify_read_fc(false, &[0], FeatureCountsOptions::default(), &mut result); // Third call: multi-mapped - classify_read_fc(true, &[0], &mut result); + classify_read_fc(true, &[0], FeatureCountsOptions::default(), &mut result); // Fourth call: no features - classify_read_fc(false, &[], &mut result); + classify_read_fc(false, &[], FeatureCountsOptions::default(), &mut result); // Fifth call: ambiguous (multiple genes) - classify_read_fc(false, &[0, 1], &mut result); + classify_read_fc(false, &[0, 1], FeatureCountsOptions::default(), &mut result); assert_eq!(result.fc_assigned, 2); assert_eq!(result.fc_multimapping, 1); @@ -1897,6 +1970,85 @@ mod tests { assert_eq!(result.gene_counts[1].fc_reads, 0); } + #[test] + fn test_fc_classify_count_multi_mapping_option() { + // -M: multi-mapping reads are counted instead of being reported as + // Unassigned_MultiMapping. + let options = FeatureCountsOptions { + count_multi_mapping: true, + count_multi_overlapping: false, + }; + let mut result = make_test_chrom_result(3); + + classify_read_fc(true, &[0], options, &mut result); + classify_read_fc(true, &[1], options, &mut result); + + assert_eq!(result.fc_multimapping, 0, "no read left unassigned"); + assert_eq!(result.fc_assigned, 2); + assert_eq!(result.gene_counts[0].fc_reads, 1); + assert_eq!(result.gene_counts[1].fc_reads, 1); + } + + #[test] + fn test_fc_classify_count_multi_overlapping_option() { + // -O: a read overlapping several genes is counted once per gene and + // reported as Assigned rather than Ambiguous. + let options = FeatureCountsOptions { + count_multi_mapping: false, + count_multi_overlapping: true, + }; + let mut result = make_test_chrom_result(3); + + classify_read_fc(false, &[0, 1, 2], options, &mut result); + + assert_eq!(result.fc_ambiguous, 0); + assert_eq!(result.fc_assigned, 1, "the read is assigned once"); + assert_eq!(result.gene_counts[0].fc_reads, 1); + assert_eq!(result.gene_counts[1].fc_reads, 1); + assert_eq!(result.gene_counts[2].fc_reads, 1); + } + + #[test] + fn test_fc_biotype_count_multi_overlapping_option() { + // -O at the biotype level: one count per distinct biotype meta-feature. + let options = FeatureCountsOptions { + count_multi_mapping: false, + count_multi_overlapping: true, + }; + let mut result = make_test_chrom_result(3); + result.biotype_reads = vec![0; 2]; + // gene 0 → biotype 0, gene 1 → biotype 1, gene 2 → biotype 0 + let gene_to_biotype: Vec = vec![0, 1, 0]; + let mut buf: Vec = Vec::new(); + + classify_read_fc_biotype( + false, + &[0, 1], + &gene_to_biotype, + &mut buf, + options, + &mut result, + ); + + assert_eq!(result.fc_biotype_ambiguous, 0); + assert_eq!(result.fc_biotype_assigned, 1); + assert_eq!(result.biotype_reads[0], 1); + assert_eq!(result.biotype_reads[1], 1); + + // Two genes of the SAME biotype remain a single meta-feature hit + classify_read_fc_biotype( + false, + &[0, 2], + &gene_to_biotype, + &mut buf, + options, + &mut result, + ); + assert_eq!(result.fc_biotype_assigned, 2); + assert_eq!(result.biotype_reads[0], 2); + assert_eq!(result.biotype_reads[1], 1); + } + // --- Chromosome partitioning tests --- #[test] diff --git a/src/rna/rseqc/accumulators.rs b/src/rna/rseqc/accumulators.rs index b91a409e..c4d45414 100644 --- a/src/rna/rseqc/accumulators.rs +++ b/src/rna/rseqc/accumulators.rs @@ -50,6 +50,8 @@ pub struct RseqcAnnotations<'a> { /// TIN index for transcript integrity number calculation. pub tin_index: Option<&'a super::tin::TinIndex>, + /// BED intervals for split_bam interval-based classification. + pub bed_intervals: Option<&'a super::split_bam::BedIntervals>, } /// Per-tool configuration parameters. @@ -97,6 +99,8 @@ pub struct RseqcConfig { pub inner_distance_enabled: bool, /// Whether TIN analysis is enabled. pub tin_enabled: bool, + /// Whether split_bam BED-interval classification is enabled. + pub split_bam_enabled: bool, /// Number of equally-spaced sampling positions per transcript for TIN. pub tin_sample_size: usize, /// Minimum number of read starts for a transcript to compute TIN. @@ -2123,6 +2127,8 @@ pub struct RseqcAccumulators { pub tin: Option, /// preseq library complexity accumulator (`None` when disabled). pub preseq: Option, + /// split_bam BED-interval classification accumulator (`None` when disabled). + pub split_bam: Option, } impl RseqcAccumulators { @@ -2138,6 +2144,7 @@ impl RseqcAccumulators { inner_dist: None, tin: None, preseq: None, + split_bam: None, } } @@ -2196,6 +2203,11 @@ impl RseqcAccumulators { } else { None }, + split_bam: if config.split_bam_enabled { + Some(super::split_bam::SplitBamAccum::default()) + } else { + None + }, } } @@ -2283,6 +2295,14 @@ impl RseqcAccumulators { if let Some(ref mut accum) = self.preseq { accum.process_read(record); } + + // split_bam: classifies every record (including secondary and + // supplementary alignments) against the BED intervals. + if let (Some(ref mut accum), Some(intervals)) = + (&mut self.split_bam, annotations.bed_intervals) + { + accum.process_read(record, chrom, intervals); + } } /// Merge another set of accumulators into this one. @@ -2314,6 +2334,9 @@ impl RseqcAccumulators { if let (Some(ref mut a), Some(b)) = (&mut self.preseq, other.preseq) { a.merge(b); } + if let (Some(ref mut a), Some(b)) = (&mut self.split_bam, other.split_bam) { + a.merge(b); + } } } diff --git a/src/rna/rseqc/mod.rs b/src/rna/rseqc/mod.rs index a4730ccb..967a2809 100644 --- a/src/rna/rseqc/mod.rs +++ b/src/rna/rseqc/mod.rs @@ -16,5 +16,6 @@ pub mod junction_annotation; pub mod junction_saturation; pub mod read_distribution; pub mod read_duplication; +pub mod split_bam; pub mod stats; pub mod tin; diff --git a/src/rna/rseqc/split_bam.rs b/src/rna/rseqc/split_bam.rs new file mode 100644 index 00000000..350eed79 --- /dev/null +++ b/src/rna/rseqc/split_bam.rs @@ -0,0 +1,446 @@ +//! BED-interval read classification, equivalent to RSeQC's `split_bam.py`. +//! +//! Classifies every alignment record against a BED file of genomic intervals +//! (typically rRNA regions) into three categories — `in` (overlapping an +//! interval), `ex` (not overlapping) and `junk` (unmapped or QC-failed) — and +//! reports the counts and percentages. +//! +//! This complements the GTF/biotype route to rRNA quantification, which +//! structurally under-reports rRNA on stock GRCh38/GENCODE builds: the 45S +//! rDNA repeat is not annotated there, and rDNA multi-mappers are dropped by +//! featureCounts' default single-hit rule. Interval overlap catches both. + +use std::collections::HashMap; +use std::io::BufRead; +use std::path::Path; + +use anyhow::{Context, Result}; +use coitrees::{BasicCOITree, Interval, IntervalTree}; +use log::debug; + +use crate::rna::bam_flags::{BAM_FQCFAIL, BAM_FUNMAP}; +use rust_htslib::bam; + +/// Per-chromosome interval tree over the BED regions. +type ChromIntervals = BasicCOITree<(), u32>; + +/// Genomic intervals loaded from a BED file, indexed per chromosome. +/// +/// `Debug` is implemented by hand because the coitrees tree type does not +/// implement it. +pub struct BedIntervals { + /// Chromosome name → interval tree of that chromosome's BED records. + trees: HashMap, + /// Total number of intervals loaded (for logging). + pub num_intervals: usize, +} + +impl std::fmt::Debug for BedIntervals { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + f.debug_struct("BedIntervals") + .field("chromosomes", &self.trees.len()) + .field("num_intervals", &self.num_intervals) + .finish() + } +} + +impl BedIntervals { + /// Return `true` when `pos` (0-based) falls inside any interval on `chrom`. + fn contains(&self, chrom: &str, pos: i64) -> bool { + match self.trees.get(chrom) { + Some(tree) => { + let p = pos as i32; + tree.query_count(p, p) > 0 + } + None => false, + } + } + + /// Return `true` when the chromosome is present in the BED file at all. + fn has_chrom(&self, chrom: &str) -> bool { + self.trees.contains_key(chrom) + } +} + +/// Parse a BED file (plain or gzip-compressed) into per-chromosome interval trees. +/// +/// Only the first three columns (chrom, start, end) are used. BED coordinates +/// are 0-based half-open, and are stored as 0-based inclusive `[start, end-1]` +/// so that point queries against BAM positions line up. +/// +/// # Arguments +/// * `path` - Path to the BED file +/// +/// # Returns +/// Indexed intervals ready for overlap queries +pub fn parse_bed(path: &str) -> Result { + let reader = + crate::io::open_reader(path).with_context(|| format!("Failed to open BED file: {path}"))?; + + let mut per_chrom: HashMap>> = HashMap::new(); + let mut num_intervals = 0usize; + + for (lineno, line) in reader.lines().enumerate() { + let line = line.with_context(|| format!("Failed to read line from BED file: {path}"))?; + let trimmed = line.trim_end(); + // Skip comments, empty lines and browser/track headers + if trimmed.is_empty() + || trimmed.starts_with('#') + || trimmed.starts_with("track") + || trimmed.starts_with("browser") + { + continue; + } + + let mut fields = trimmed.split('\t'); + let (Some(chrom), Some(start_s), Some(end_s)) = + (fields.next(), fields.next(), fields.next()) + else { + anyhow::bail!( + "Malformed BED file '{}' at line {}: expected at least 3 tab-separated columns", + path, + lineno + 1 + ); + }; + + let start: i64 = start_s.parse().with_context(|| { + format!( + "Malformed BED file '{}' at line {}: invalid start position '{}'", + path, + lineno + 1, + start_s + ) + })?; + let end: i64 = end_s.parse().with_context(|| { + format!( + "Malformed BED file '{}' at line {}: invalid end position '{}'", + path, + lineno + 1, + end_s + ) + })?; + anyhow::ensure!( + end > start, + "Malformed BED file '{}' at line {}: end ({}) must be greater than start ({})", + path, + lineno + 1, + end, + start + ); + + per_chrom + .entry(chrom.to_string()) + .or_default() + // BED is half-open; coitrees intervals are inclusive on both ends + .push(Interval::new(start as i32, (end - 1) as i32, ())); + num_intervals += 1; + } + + anyhow::ensure!( + num_intervals > 0, + "No usable intervals found in BED file '{}'", + path + ); + + let trees = per_chrom + .into_iter() + .map(|(chrom, ivs)| (chrom, ChromIntervals::new(&ivs))) + .collect(); + + debug!("Loaded {num_intervals} intervals from BED file {path}"); + Ok(BedIntervals { + trees, + num_intervals, + }) +} + +// =================================================================== +// Accumulator +// =================================================================== + +/// Per-worker accumulator classifying reads against the BED intervals. +#[derive(Debug, Default, Clone)] +pub struct SplitBamAccum { + /// Records overlapping at least one BED interval. + pub in_count: u64, + /// Mapped, QC-passing records not overlapping any BED interval. + pub ex_count: u64, + /// Unmapped or QC-failed records. + pub junk_count: u64, +} + +impl SplitBamAccum { + /// Classify a single alignment record. + /// + /// Follows RSeQC `split_bam.py`: unmapped and QC-failed records are "junk"; + /// otherwise the read's own start position, and (when the mate is mapped) + /// the mate's start position, are tested against the intervals. A hit on + /// either end puts the record in the "in" category. + /// + /// Secondary and supplementary alignments are *not* filtered out — this is + /// deliberate, and is what lets the interval route see the rDNA + /// multi-mappers that featureCounts' default single-hit rule discards. + /// + /// # Arguments + /// * `record` - The alignment record + /// * `chrom` - Chromosome name for this record, already mapped to BED naming + /// * `intervals` - The BED intervals to test against + pub fn process_read(&mut self, record: &bam::Record, chrom: &str, intervals: &BedIntervals) { + let flags = record.flags(); + if flags & (BAM_FUNMAP | BAM_FQCFAIL) != 0 { + self.junk_count += 1; + return; + } + + if !intervals.has_chrom(chrom) { + self.ex_count += 1; + return; + } + + let hit = intervals.contains(chrom, record.pos()) + || (record.mtid() == record.tid() + && record.mpos() >= 0 + && intervals.contains(chrom, record.mpos())); + + if hit { + self.in_count += 1; + } else { + self.ex_count += 1; + } + } + + /// Merge another worker's counts into this accumulator. + pub fn merge(&mut self, other: SplitBamAccum) { + self.in_count += other.in_count; + self.ex_count += other.ex_count; + self.junk_count += other.junk_count; + } + + /// Finalise into a result. + pub fn into_result(self) -> SplitBamResult { + SplitBamResult { + in_count: self.in_count, + ex_count: self.ex_count, + junk_count: self.junk_count, + } + } +} + +/// Final classification counts. +#[derive(Debug, Clone)] +pub struct SplitBamResult { + /// Records overlapping at least one BED interval. + pub in_count: u64, + /// Mapped, QC-passing records not overlapping any BED interval. + pub ex_count: u64, + /// Unmapped or QC-failed records. + pub junk_count: u64, +} + +impl SplitBamResult { + /// Total records classified. + pub fn total(&self) -> u64 { + self.in_count + self.ex_count + self.junk_count + } + + /// Percentage of mapped, QC-passing records that fall inside the intervals. + /// + /// Junk records are excluded from the denominator, so the value answers + /// "what fraction of usable alignments are rRNA?". + pub fn percent_in(&self) -> f64 { + let denom = self.in_count + self.ex_count; + if denom == 0 { + 0.0 + } else { + 100.0 * self.in_count as f64 / denom as f64 + } + } +} + +// =================================================================== +// Output +// =================================================================== + +/// Write the classification summary as a TSV file. +/// +/// # Arguments +/// * `result` - The computed classification counts +/// * `bed_path` - BED file the intervals came from (recorded in the header) +/// * `output_path` - Path to write the summary to +pub fn write_split_bam_summary( + result: &SplitBamResult, + bed_path: &str, + output_path: &Path, +) -> Result<()> { + use std::io::Write; + + let mut out = std::fs::File::create(output_path).with_context(|| { + format!( + "Failed to create split_bam summary file: {}", + output_path.display() + ) + })?; + + let total = result.total(); + let usable = result.in_count + result.ex_count; + let pct = |n: u64, d: u64| { + if d == 0 { + 0.0 + } else { + 100.0 * n as f64 / d as f64 + } + }; + + writeln!(out, "# BED intervals: {bed_path}")?; + writeln!( + out, + "# 'in' = alignment overlaps an interval, 'ex' = does not, \ + 'junk' = unmapped or QC-failed" + )?; + writeln!( + out, + "# percent_of_usable excludes junk records from the denominator" + )?; + writeln!(out, "category\tcount\tpercent_of_total\tpercent_of_usable")?; + writeln!( + out, + "in\t{}\t{:.4}\t{:.4}", + result.in_count, + pct(result.in_count, total), + pct(result.in_count, usable) + )?; + writeln!( + out, + "ex\t{}\t{:.4}\t{:.4}", + result.ex_count, + pct(result.ex_count, total), + pct(result.ex_count, usable) + )?; + writeln!( + out, + "junk\t{}\t{:.4}\tNA", + result.junk_count, + pct(result.junk_count, total) + )?; + writeln!(out, "total\t{total}\t100.0000\tNA")?; + + debug!( + "Wrote split_bam summary to {} ({} in, {} ex, {} junk)", + output_path.display(), + result.in_count, + result.ex_count, + result.junk_count + ); + Ok(()) +} + +// =================================================================== +// Tests +// =================================================================== + +#[cfg(test)] +mod tests { + use super::*; + + fn write_temp(content: &str, ext: &str) -> std::path::PathBuf { + use std::io::Write; + 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_split_bam_test_{:?}_{}.{}", + std::thread::current().id(), + id, + ext + )); + let mut f = std::fs::File::create(&path).unwrap(); + f.write_all(content.as_bytes()).unwrap(); + f.flush().unwrap(); + path + } + + #[test] + fn test_parse_bed_basic() { + let path = write_temp( + "track name=rRNA\n\ + # a comment\n\ + chr1\t100\t200\trRNA1\t0\t+\n\ + chr1\t300\t400\n\ + chr2\t50\t60\n", + "bed", + ); + let intervals = parse_bed(path.to_str().unwrap()).unwrap(); + let _ = std::fs::remove_file(&path); + + assert_eq!(intervals.num_intervals, 3); + // BED is half-open: [100, 200) covers 100..=199 + assert!(intervals.contains("chr1", 100)); + assert!(intervals.contains("chr1", 199)); + assert!(!intervals.contains("chr1", 200)); + assert!(!intervals.contains("chr1", 99)); + assert!(intervals.contains("chr2", 55)); + assert!(!intervals.has_chrom("chr3")); + } + + #[test] + fn test_parse_bed_rejects_malformed_line() { + let path = write_temp("chr1\t100\n", "bed"); + let err = parse_bed(path.to_str().unwrap()).unwrap_err(); + let _ = std::fs::remove_file(&path); + assert!( + err.to_string().contains("at least 3 tab-separated columns"), + "unexpected error: {err}" + ); + } + + #[test] + fn test_parse_bed_rejects_empty_file() { + let path = write_temp("# nothing here\n", "bed"); + let err = parse_bed(path.to_str().unwrap()).unwrap_err(); + let _ = std::fs::remove_file(&path); + assert!( + err.to_string().contains("No usable intervals"), + "unexpected error: {err}" + ); + } + + #[test] + fn test_classification_counts_multimappers_and_mates() { + use rust_htslib::bam::Read as BamRead; + + let bed = write_temp("chr1\t500\t600\n", "bed"); + let intervals = parse_bed(bed.to_str().unwrap()).unwrap(); + let _ = std::fs::remove_file(&bed); + + // r1: primary, starts inside the interval -> in + // r2: secondary (NH>1), starts inside the interval -> in (not filtered) + // r3: primary, outside -> ex + // r4: primary, outside but mate starts inside -> in + // r5: unmapped -> junk + let sam = "\ +@HD\tVN:1.6\tSO:coordinate\n\ +@SQ\tSN:chr1\tLN:20000\n\ +r1\t99\tchr1\t501\t30\t10M\t=\t900\t0\tACGTACGTAC\tIIIIIIIIII\n\ +r2\t355\tchr1\t520\t30\t10M\t=\t900\t0\tACGTACGTAC\tIIIIIIIIII\n\ +r3\t99\tchr1\t900\t30\t10M\t=\t950\t0\tACGTACGTAC\tIIIIIIIIII\n\ +r4\t147\tchr1\t950\t30\t10M\t=\t520\t0\tACGTACGTAC\tIIIIIIIIII\n\ +r5\t77\tchr1\t960\t0\t*\t*\t0\t0\tACGTACGTAC\tIIIIIIIIII\n"; + let sam_path = write_temp(sam, "sam"); + + let mut reader = bam::Reader::from_path(&sam_path).unwrap(); + let mut accum = SplitBamAccum::default(); + let mut record = bam::Record::new(); + while let Some(res) = reader.read(&mut record) { + res.unwrap(); + accum.process_read(&record, "chr1", &intervals); + } + let _ = std::fs::remove_file(&sam_path); + + let result = accum.into_result(); + assert_eq!(result.in_count, 3, "in"); + assert_eq!(result.ex_count, 1, "ex"); + assert_eq!(result.junk_count, 1, "junk"); + assert_eq!(result.total(), 5); + assert!((result.percent_in() - 75.0).abs() < 1e-9, "percent_in"); + } +}