rna-gc-bias-metrics

Measure GC-content bias in RNA-seq coverage.

Given a transcriptome FASTA and a BAM of reads aligned to it, rna-gc-bias-metrics computes per-bin GC content and coverage depth, then summarizes how normalized coverage varies with GC fraction. This reveals whether a library is over- or under-covering GC-rich (or GC-poor) regions — a common artifact of library preparation and sequencing chemistry.

How it works

For each transcript the pipeline:

1. Reads the BAM (via polars-bio) and drops unmapped, secondary, and supplementary alignments, as well as paired reads that are not properly paired (mate-unmapped or discordant mates). Single-end reads and proper pairs are kept.
2. Expands each read's CIGAR into the reference positions it covers, and deduplicates the overlap between properly paired mates so shared bases are counted once.
3. Collapses coverage into a per-range bedgraph, then redistributes it into fixed-length bins (default 100 bp), splitting partial-bin coverage proportionally.
4. Computes each bin's GC fraction from the transcript sequence and normalizes its depth against the transcript's mean bin depth (so 1.0 is the transcript average).
5. Aggregates mean normalized depth and bin counts per rounded GC fraction across the whole input.

Requirements

• Python ≥ 3.12
• numpy, polars, polars-bio

Dependencies are managed with either uv (PyPI) or pixi (conda + PyPI). Pick whichever you already use.

Installation

git clone git@github.com:nebiolabs/rna-gc-bias-metrics.git
cd rna-gc-bias-metrics

# with uv
uv sync

# or with pixi
pixi install

Input requirements

• Transcriptome FASTA, either line-wrapped (as samtools faidx, GENCODE and Ensembl emit by default) or with single-line sequences. Both parse identically.
• A FASTA index named <name>.fa.fai alongside it (created with samtools faidx).
• A BAM aligned to that transcriptome (each RNAME is a transcript ID). BAM reference names must match the FASTA names, or the run aborts with an error.

Every reference name — in the FASTA, the .fa.fai and the BAM — is taken as the header up to the first whitespace, so descriptive headers need no cleanup:

>ENST00000227525.8 cdna chromosome:GRCh38:12:6534517:6538371:1 gene:ENSG00000111640.15

is matched against a BAM RNAME of ENST00000227525.8. This mirrors the aligners (the SAM spec forbids whitespace in RNAME, so only the leading word survives alignment) and samtools faidx, which truncates the same way.

To index a FASTA:

samtools faidx transcripts.fa    # creates transcripts.fa.fai

Usage

The tool runs as a module. Under uv, prefix commands with uv run; under pixi, use pixi run.

uv run python -m rna_gc_bias_metrics.calculate_gc_coverage \
    transcripts.fa \
    reads.bam \
    --report_bin_count_for_full_transcriptome \
    -o gc_bias_profile.tsv

Arguments

Output

A tab-separated table, one row per rounded GC fraction. Every transcript with any coverage contributes all of its bins, so bins no read reached count as zero depth rather than being dropped — that is what makes a systematically uncovered GC range visible as a depleted signal instead of an absent one.

Example (run against the bundled test fixtures — two of the four transcripts have coverage, contributing 30 bins between them):

gc_fraction  mean_normalized_depth  transcriptome_bin_count  transcriptome_bin_count_all_transcripts
0.24         0.0                    1                        1
0.25                                                         1
0.28         0.0                    1                        1
0.33         0.0                    1                        2
0.36                                                         1
0.43         0.0                    1                        1
0.46                                                         1
0.47         10.888888888888891     1                        1
0.48         0.0                    1                        2
0.49                                                         2
0.5          0.0                    1                        1
0.51         0.0                    1                        1
0.52                                                         1
0.53         0.0                    3                        5
0.55                                                         1
0.56         0.0                    3                        4
0.57                                                         1
0.58         0.0                    1                        5
0.59         5.5555555555555545     2                        3
0.6          0.0                    1                        2
0.61         0.0                    1                        3
0.62         0.0                    2                        3
0.64                                                         2
0.65         1.208888888888889      3                        3
0.66                                                         2
0.67         0.0                    2                        3
0.68         0.5333333333333337     1                        1
0.7          0.37333333333333474    1                        1
0.71         0.0                    1                        2
0.77         3.466666666666665      1                        1

(GC fractions occurring only in transcripts with no coverage at all have empty depth/count columns, but still appear when the full-transcriptome background is reported.)

Python API

The pipeline steps are importable for use in scripts or notebooks:

from rna_gc_bias_metrics.calculate_gc_coverage import (
    load_sequences,
    calculate_gc_coverage,
    calculate_gc_pct_coverage,
)

sequences, faidx = load_sequences("transcripts.fa", "transcripts.fa.fai")
bin_cov_with_gc = calculate_gc_coverage("reads.bam", sequences, faidx, fixed_length_bin_bp=100)

# mean normalized depth per rounded GC fraction
per_gc_coverage = calculate_gc_pct_coverage(bin_cov_with_gc).collect()

calculate_gc_coverage returns a per-bin table (rname, bin_start, depth_fractional, gc_frac, depth_normalized, gc_frac_rounded) as a Polars LazyFrame; call .collect() to materialize it. It holds every bin of every transcript that has coverage somewhere, with depth_fractional 0 for bins no read reached. load_bam, expand_cigar, bam_to_bedgraph, get_binned_coverage, and get_bin_gc expose the individual stages — note that get_binned_coverage alone emits only bins that received coverage; the zero bins are filled in by get_bin_gc.

Notes and limitations

• Coverage is deduplicated across the overlap of properly paired mates so shared bases are counted once. Paired reads that are not properly paired (mate-unmapped or discordant mates) are dropped. Read-through ("dovetail") mate tails past the mate's end are dropped.

Development

uv run pytest test/ -v      # with uv
pixi run test               # with pixi

Tests should pass under both workflows. See AGENTS.md for project layout and details on keeping the two dependency ecosystems in sync.

License

Distributed under the GNU Affero General Public License v3.0. See LICENSE.txt.