rsomics_compute_gc_bias/

lib.rs

//! GC-bias estimation for a BAM over a 2bit genome — deeptools computeGCBias.
//!
//! For each genomic position sampled at `stepSize`, the fragment-length window
//! contributes its G+C count to `n_gc` (expected, from the genome) and the
//! number of forward reads starting at that position to `f_gc` (observed). The
//! per-GC ratio `r_gc = (f/n) * (sum n / sum f)` is the bias estimate
//! (Benjamini & Speed 2012, doi:10.1093/nar/gks001).

#![allow(clippy::cast_possible_truncation)]
#![allow(clippy::cast_precision_loss)]
#![allow(clippy::cast_sign_loss)]

mod twobit;

use std::io::Write;
use std::path::Path;

use rsomics_bamio::raw::{self, RawRecord};
use rsomics_common::{Result, RsomicsError};

pub use twobit::TwoBit;

const FLAG_UNMAPPED: u16 = 0x4;
const FLAG_REVERSE: u16 = 0x10;

pub struct GcBiasOpts {
    pub fragment_length: u32,
    pub effective_genome_size: u64,
    /// deeptools default 5e7; sets `stepSize = max(genome_size / sample_size, 1)`.
    pub sample_size: u64,
}

/// One row per GC count 0..=fragment_length.
pub struct GcBiasTable {
    pub f_gc: Vec<u64>,
    pub n_gc: Vec<u64>,
    pub r_gc: Vec<f64>,
}

/// Per-chromosome forward-read start histogram plus the mapped-read total.
struct ReadStarts {
    starts: Vec<Vec<u32>>,
    total_mapped: u64,
}

pub fn compute(bam: &Path, genome: &Path, opts: &GcBiasOpts) -> Result<GcBiasTable> {
    let tbit = TwoBit::open(genome)?;
    let genome_size = tbit.genome_size();
    let step = (genome_size / opts.sample_size).max(1);

    let mut reader = rsomics_bamio::open_with_workers(bam, std::num::NonZero::<usize>::MIN)?;
    let header = reader.read_header().map_err(RsomicsError::Io)?;
    let refs: Vec<(String, u32)> = header
        .reference_sequences()
        .iter()
        .map(|(name, seq)| (name.to_string(), usize::from(seq.length()) as u32))
        .collect();

    let reads = collect_read_starts(reader.get_mut(), &refs)?;

    let reads_per_bp = reads.total_mapped as f64 / opts.effective_genome_size as f64;
    let max_reads = poisson_isf(
        4.0 * reads_per_bp * f64::from(opts.fragment_length),
        1.0 / opts.sample_size as f64,
    );

    let flen = opts.fragment_length;
    let mut f_gc = vec![0u64; flen as usize + 1];
    let mut n_gc = vec![0u64; flen as usize + 1];

    for (tid, (name, _len)) in refs.iter().enumerate() {
        let Some(chrom_len) = tbit.chrom_len(name) else {
            continue;
        };
        let starts = &reads.starts[tid];
        let mut i = 0u32;
        while i < chrom_len {
            if i + flen > chrom_len {
                break;
            }
            let Some(w) = tbit.window_bases(name, i, i + flen) else {
                i = i.saturating_add(step as u32);
                continue;
            };
            // deeptools skips windows with > 5% Ns
            if (w.non_n as f64) < 0.95 * f64::from(flen) {
                i = i.saturating_add(step as u32);
                continue;
            }
            let gc = w.gc as usize;
            n_gc[gc] += 1;
            let num_reads = starts.get(i as usize).copied().unwrap_or(0);
            if u64::from(num_reads) < max_reads {
                f_gc[gc] += u64::from(num_reads);
            }
            i = i.saturating_add(step as u32);
        }
    }

    let sum_n: u64 = n_gc.iter().sum();
    let sum_f: u64 = f_gc.iter().sum();
    let scaling = sum_n as f64 / sum_f as f64;
    let r_gc: Vec<f64> = (0..f_gc.len())
        .map(|x| {
            if n_gc[x] > 0 && f_gc[x] > 0 {
                (f_gc[x] as f64 / n_gc[x] as f64) * scaling
            } else {
                1.0
            }
        })
        .collect();

    Ok(GcBiasTable { f_gc, n_gc, r_gc })
}

/// Build, per reference, the count of forward (non-reverse) reads starting at
/// each position. deeptools counts a read at GC position `i` iff its `pos == i`
/// and it is forward and mapped.
fn collect_read_starts(
    reader: &mut impl std::io::Read,
    refs: &[(String, u32)],
) -> Result<ReadStarts> {
    let mut starts: Vec<Vec<u32>> = refs
        .iter()
        .map(|(_, len)| vec![0u32; *len as usize])
        .collect();
    let mut total_mapped = 0u64;
    let mut record = RawRecord::default();

    while raw::read_record(reader, &mut record)? != 0 {
        let flags = record.flags();
        if flags & FLAG_UNMAPPED != 0 {
            continue;
        }
        let tid = record.reference_sequence_id();
        if tid < 0 {
            continue;
        }
        total_mapped += 1;
        if flags & FLAG_REVERSE != 0 {
            continue;
        }
        let pos = record.alignment_start();
        if pos < 0 {
            continue;
        }
        if let Some(chrom) = starts.get_mut(tid as usize)
            && let Some(slot) = chrom.get_mut(pos as usize)
        {
            *slot += 1;
        }
    }

    Ok(ReadStarts {
        starts,
        total_mapped,
    })
}

/// scipy `poisson(mu).isf(q)`: the smallest integer k with P(X > k) <= q.
/// Computed from the exact PMF recurrence; mu here is small (a few tens).
fn poisson_isf(mu: f64, q: f64) -> u64 {
    if mu <= 0.0 {
        return 0;
    }
    let mut pmf = (-mu).exp();
    let mut cdf = pmf;
    let mut k = 0u64;
    // sf(k) = 1 - cdf(k); stop at the first k whose tail probability <= q
    while 1.0 - cdf > q {
        k += 1;
        pmf *= mu / k as f64;
        cdf += pmf;
        if pmf == 0.0 && 1.0 - cdf <= q {
            break;
        }
    }
    k
}

/// deeptools writes the table with numpy `savetxt` defaults: three
/// space-separated columns in `%.18e`, one row per GC count.
pub fn write_freq(table: &GcBiasTable, out: &mut impl Write) -> Result<()> {
    for i in 0..table.f_gc.len() {
        writeln!(
            out,
            "{} {} {}",
            fmt_e(table.f_gc[i] as f64),
            fmt_e(table.n_gc[i] as f64),
            fmt_e(table.r_gc[i]),
        )
        .map_err(RsomicsError::Io)?;
    }
    Ok(())
}

/// numpy savetxt default `%.18e`: 18 fractional digits, signed exponent of at
/// least two digits — Rust's `{:e}` drops the sign and zero-padding, so rebuild.
fn fmt_e(v: f64) -> String {
    let raw = format!("{v:.18e}");
    let (mantissa, exp) = raw.split_once('e').expect("scientific format has e");
    let exp: i32 = exp.parse().expect("exponent is an integer");
    let sign = if exp < 0 { '-' } else { '+' };
    format!("{mantissa}e{sign}{:02}", exp.abs())
}