feat: BED-interval rRNA quantification (split_bam) and featureCounts -M/-O

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`,

docs/src/content/docs/rna/featurecounts.mdx

@@ -115,6 +115,47 @@ The biotype is extracted from the GTF attribute specified by `biotype_attribute`
 - `<sample>.biotype_counts_mqc.tsv` -- Biotype counts formatted as a MultiQC bargraph data file, suitable for visualizing the distribution of reads across biotypes.
 - `<sample>.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
+```
+
+<Aside type="caution" title="Biotype %rRNA under-reports on stock GRCh38/GENCODE">
+Biotype counting can only see rRNA that is annotated in the GTF and captured by
+the counting rules, and on a stock GRCh38 / GENCODE setup neither holds for the
+bulk of rRNA reads:
+
+- The 45S rDNA repeat (18S / 5.8S / 28S) is **not annotated** in GENCODE GRCh38 —
+  it lives on the unassembled acrocentric short arms. Uniquely-mapped rRNA reads
+  therefore land in `Unassigned_NoFeatures`.
+- rDNA is a high-copy repeat, so most rRNA reads multi-map and are dropped as
+  `Unassigned_MultiMapping` under the default rules.
+
+The result is a biotype %rRNA close to zero even when true residual rRNA is
+several percent, so it cannot be used to judge depletion. Use `-M` (and
+optionally `-O`) to stop discarding those reads, and prefer the interval-based
+[`split_bam`](../rseqc/#split_bam-bed-interval-classification) route with an
+rRNA BED file for a robust rRNA fraction.
+</Aside>
+
 ## 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.

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?".
+
+<Aside type="note" title="Multi-mappers are included on purpose">
+Unlike the featureCounts pass, split_bam does **not** filter secondary or
+supplementary alignments. That is what lets it see the rDNA multi-mappers that
+featureCounts' default single-hit rule discards.
+</Aside>
+
+<Aside type="caution" title="No split BAM files are written">
+RSeQC's `split_bam.py` writes three BAM files (`.in.bam`, `.ex.bam`,
+`.junk.bam`). RustQC reports the counts only — writing three extra BAM copies
+would defeat the point of the single-pass design. If you need the reads
+themselves, use `samtools view -L`.
+</Aside>
+
+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

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>`
+
+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 <NAME>`
 
 Override the sample name used for output filenames. By default, the sample name

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

src/cli.rs

@@ -155,6 +155,33 @@ pub struct RnaArgs {
     )]
     pub biotype_attribute: Option<String>,
 
+    /// 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<String>,
+
+    /// 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,

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<String>,
+}
+
 /// Configuration for samtools idxstats-compatible output.
 ///
 /// When enabled, produces a file matching `samtools idxstats` output format,