Skip to content

feat: add geneBody_coverage and read_GC RSeQC tools - #134

Open
BenjaminDEMAILLE wants to merge 1 commit into
seqeralabs:mainfrom
BenjaminDEMAILLE:feat/genebody-coverage-read-gc
Open

feat: add geneBody_coverage and read_GC RSeQC tools#134
BenjaminDEMAILLE wants to merge 1 commit into
seqeralabs:mainfrom
BenjaminDEMAILLE:feat/genebody-coverage-read-gc

Conversation

@BenjaminDEMAILLE

Copy link
Copy Markdown

Closes #127.

Adds the two remaining RSeQC tools, both inside the existing single pass.

geneBody_coverage

src/rna/rseqc/gene_body_coverage.rs — 5' → 3' coverage in 100 percentile bins.

  • Each transcript's exonic bases are reduced to 100 percentile positions with a direct port of RSeQC's mystat.percentile_list
  • Reverse-strand transcripts are flipped so bin 1 is always the 5' end
  • Transcripts under 100 bp are skipped (RSeQC's --minimum_length default)
  • Read filters: unmapped, QC-fail, secondary, duplicate; deletions and N skips do not count as covered

Outputs {sample}.geneBodyCoverage.txt (RSeQC format), .geneBodyCoverage.r, and a native .curves.png/.svg.

One deliberate deviation, documented in the docs and module: RSeQC profiles every entry of a BED12 reference (usually a housekeeping-gene set); RustQC has no BED input here and profiles one representative transcript per gene from the GTF — the one with the most exonic bases, the same choice the TIN analysis already makes.

⚠️ Worth a maintainer's eye: the port needs Python's round-half-to-even. With f64::round (half away from zero), interpolated percentile positions land one base off and 7 of 100 bins came out ±1 against upstream. There's a unit test pinning the rounding behaviour and a note in AGENTS.md.

read_GC

src/rna/rseqc/read_gc.rs — GC distribution of mapped reads.

  • Filters match upstream ParseBAM.readGC: skip unmapped, QC-fail and below the MAPQ cutoff; secondary alignments and duplicates are counted (upstream does not filter them)
  • GC percentage formatted to two decimals, as upstream does

Outputs {sample}.GC.xls, .GC_plot.r, and a native .GC_plot.png/.svg.

Two deliberate deviations: rows are sorted ascending by GC percentage (upstream emits Python dict order, which is BAM-order dependent and not reproducible under parallel processing), and reads stored without SEQ are skipped with a warning instead of dividing by zero.

Validation against RSeQC 5.05

Ran upstream geneBody_coverage.py and read_GC.py (RSeQC 5.05, pysam) on tests/data/test.bam, using a BED12 built from tests/data/test.gtf with the same representative-transcript choice:

  • gene body coverage: all 100 bins identical
  • GC table: identical (0.00 483)

Both references are pinned in new integration tests. Also verified --threads 4 produces byte-identical output to single-threaded, and that the skip flags remove the outputs.

Other

  • Config: rna.gene_body_coverage.enabled, rna.read_gc.enabled (both default true); CLI --skip-gene-body-coverage, --skip-read-gc
  • Docs: new sections on the RSeQC page (including a note on how this differs from the Qualimap gene body coverage RustQC already produces), config/CLI reference entries, "8 RSeQC tools" → 10 where claimed
  • cargo test — 211 lib + 20 integration tests pass; cargo fmt --check, cargo clippy -- -D warnings clean

🤖 Generated with Claude Code

The two RSeQC tools missing from RustQC's coverage: gene body coverage
(5'->3' uniformity, the standard RNA degradation / 3'-bias signal) and the
read GC distribution (the standard GC-bias check). Both run inside the
existing single pass.

gene_body_coverage.rs samples 100 percentile positions per transcript with a
direct port of RSeQC's mystat.percentile_list, counts coverage at those
positions, and flips reverse-strand transcripts so bin 1 is always the 5'
end. The port includes Python's round-half-to-even; using f64::round shifts
interpolated positions by one base and changes bin counts. Transcripts under
100 bp are skipped, matching RSeQC's --minimum_length default.

read_gc.rs bins the GC percentage of each read's query sequence with
upstream's filters (skip unmapped, QC-fail and below-MAPQ; keep secondary
alignments and duplicates). Rows are sorted by GC percentage rather than
emitted in BAM order, which is not reproducible under parallel processing.
Reads stored without SEQ are skipped and reported instead of dividing by
zero.

Both write the RSeQC output files and R scripts, plus native PNG/SVG plots.
Toggles: --skip-gene-body-coverage / --skip-read-gc and the matching config
sections.

Validated against RSeQC 5.05 on tests/data/test.bam: all 100 gene body
coverage bins and the full GC table match exactly.

Closes seqeralabs#127

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: implement RSeQC's geneBody_coverage.py and read_GC.py equivalents

1 participant