Skip to content

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

Open
BenjaminDEMAILLE wants to merge 1 commit into
seqeralabs:mainfrom
BenjaminDEMAILLE:feat/rrna-bed-quantification
Open

feat: BED-interval rRNA quantification (split_bam) and featureCounts -M/-O#133
BenjaminDEMAILLE wants to merge 1 commit into
seqeralabs:mainfrom
BenjaminDEMAILLE:feat/rrna-bed-quantification

Conversation

@BenjaminDEMAILLE

Copy link
Copy Markdown

Closes #111.

Addresses all three asks in the issue.

1. split_bam: BED-interval classification

New module src/rna/rseqc/split_bam.rs, the counting side of RSeQC's split_bam.py. Runs inside the existing single pass, disabled unless a BED file is supplied:

rustqc rna sample.bam --gtf genes.gtf --rrna-bed GRCh38_rRNA.bed

Output rseqc/split_bam/{sample}.split_bam.tsv:

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

Classification follows split_bam.py: unmapped/QC-failed records are junk; otherwise the read's own start, and the mate's start when the mate is mapped, are tested against the intervals. Secondary and supplementary alignments are deliberately not filtered — that is exactly what lets this route see the rDNA multi-mappers that featureCounts' default single-hit rule discards.

Not implemented: the three output BAM files (.in.bam / .ex.bam / .junk.bam) that split_bam.py writes. Writing three extra BAM copies would defeat the single-pass design; the docs point at samtools view -L for that. Happy to add them behind a flag if you'd rather have full parity.

2. featureCounts -M / -O

-M, --count-multi-mapping and -O, --count-multi-overlapping (plus rna.featurecounts.count_multi_mapping / count_multi_overlapping), applied to the gene-level and biotype-level featureCounts classification. dupRadar's duplicate-rate matrix is untouched — it tracks multi-mappers separately by design.

On tests/data/test.bam, -M -O moves the 10 multi-mapping reads from Unassigned_MultiMapping into Assigned (473 → 483).

3. Docs note

A caution box on the featureCounts page spelling out why biotype %rRNA reads near zero on stock GRCh38/GENCODE (45S rDNA unannotated + multi-mapper drop), pointing to -M/-O and to split_bam. Plus a split_bam section on the RSeQC page and CLI/config reference entries.

Verification

  • samtools view -c tests/data/test.bam chr1:1001-3000 = 107, matching the in count for a chr1 1000 3000 BED
  • New unit tests: BED parsing (half-open coordinates, malformed lines, empty file), classification of a record mix including a secondary alignment and a mate-side hit, and -M / -O behaviour at gene and biotype level
  • cargo test — 207 lib + 18 integration tests pass
  • cargo fmt --check, cargo clippy -- -D warnings clean

🤖 Generated with Claude Code

Biotype counting structurally under-reports rRNA on stock GRCh38/GENCODE
builds: the 45S rDNA repeat is not annotated there (uniquely-mapped rRNA
reads land in Unassigned_NoFeatures) and rDNA is a high-copy repeat, so most
rRNA reads multi-map and are dropped as Unassigned_MultiMapping. The result
is a %rRNA near zero even when true residual rRNA is several percent.

Three changes address that:

- New split_bam module (src/rna/rseqc/split_bam.rs), the counting side of
  RSeQC's split_bam.py. Given a BED file (--rrna-bed or rna.split_bam.bed) it
  classifies every alignment as in / ex / junk and writes a TSV of counts and
  percentages. It runs inside the existing single pass, is disabled unless a
  BED file is supplied, and deliberately does not filter secondary or
  supplementary alignments, which is what lets it see the rDNA multi-mappers
  featureCounts discards. No split BAM files are written.

- featureCounts -M / -O equivalents: --count-multi-mapping and
  --count-multi-overlapping (plus config toggles). They affect the
  featureCounts gene-level and biotype-level counts only; dupRadar's matrix
  keeps its existing multi-mapper semantics.

- Docs: a caution on the featureCounts page explaining why biotype %rRNA
  under-reports and what to use instead, a split_bam section on the RSeQC
  page, and CLI/config reference entries.

Closes seqeralabs#111

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Feature request: RSeQC split_bam.py - BED-file approach to rRNA quantification

1 participant