Struct GcFragModel 

pub struct GcFragModel { /* private fields */ }

(condBins × numGCBins) table of fragment-GC counts (then a normalized distribution / ratio after normalize / gc_ratio).

Implementations

impl GcFragModel

pub fn new(cond_bins: usize, gc_bins: usize) -> Self

Examples found in repository

examples/conv_bench.rs (line 59)

fn main() {
    let args: Vec<String> = std::env::args().collect();
    let n: usize = args.get(1).and_then(|s| s.parse().ok()).unwrap_or(100_000);
    let passes: usize = args.get(2).and_then(|s| s.parse().ok()).unwrap_or(3);

    let mut rng = Lcg(0x9E3779B97F4A7C15);

    // Build a realistic transcript-length distribution: log-normal-ish around
    // ~1.5 kb with a long tail, clamped to [80, 15000] — close to GRCh38 cDNA.
    let bases = [b'A', b'C', b'G', b'T'];
    let mut seqs: Vec<Vec<u8>> = Vec::with_capacity(n);
    let mut prefixes: Vec<Vec<u32>> = Vec::with_capacity(n);
    for _ in 0..n {
        let u = rng.next_f64();
        // exp of a normal-ish variate -> heavy right tail
        let len = (7.3 + 0.9 * (u - 0.5) * 4.0).exp() as usize;
        let len = len.clamp(80, 15_000);
        let seq: Vec<u8> = (0..len)
            .map(|_| bases[(rng.next_u32() & 3) as usize])
            .collect();
        prefixes.push(gc_prefix(&seq));
        seqs.push(seq);
    }

    // A normalized 3×25 GC ratio model (the shape doesn't matter for timing; we
    // exercise the same binning + lookup the real run does).
    let mut obs = GcFragModel::new(DEFAULT_COND_BINS, DEFAULT_GC_BINS);
    let mut exp = GcFragModel::new(DEFAULT_COND_BINS, DEFAULT_GC_BINS);
    for ctx in 0..=100 {
        for gc in 0..=100 {
            obs.inc(gc, ctx, 1.0 + 0.3 * ((gc + ctx) as f64).sin());
            exp.inc(gc, ctx, 1.0);
        }
    }
    let gc_model = gc_ratio(&mut obs, &mut exp, 1000.0);

    // Fragment-length CDF: gaussian-ish pmf around mean 250, sd 40, over 0..1000.
    let fld_max = 1000usize;
    let mean = 250.0f64;
    let sd = 40.0f64;
    let pmf: Vec<f64> = (0..=fld_max)
        .map(|l| {
            let z = (l as f64 - mean) / sd;
            (-0.5 * z * z).exp()
        })
        .collect();
    let (cdf, fld_low, fld_high) = salmon_model::seqbias::fld_cdf_and_bounds(&pmf);

    let mut acc = 0.0f64;
    let mut best = f64::INFINITY;
    for p in 0..passes {
        let t = Instant::now();
        for (seq, prefix) in seqs.iter().zip(&prefixes) {
            let ref_len = seq.len() as f64;
            let elen = (ref_len - 200.0).max(1.0); // ensure unprocessed > 0
            let bias = BiasInputs {
                seq: None,
                gc: Some((&gc_model, salmon_model::GcView::Dense(prefix.as_slice()))),
                pos: None,
            };
            acc += corrected_effective_length_full(
                seq,
                &cdf,
                fld_low,
                fld_high,
                &bias,
                elen,
                GC_SAMP_STRIDE,
                false,
            );
        }
        let dt = t.elapsed().as_secs_f64();
        best = best.min(dt);
        eprintln!(
            "pass {p}: {:.3}s  ({:.2} µs/transcript)  acc={}",
            dt,
            dt * 1e6 / n as f64,
            black_box(acc)
        );
    }
    eprintln!(
        "BEST: {:.3}s over {} transcripts ({:.2} µs/transcript)",
        best,
        n,
        best * 1e6 / n as f64
    );
}

pub fn default_model() -> Self

salmon’s default 3 × 101 model.

pub fn dump(&self) -> &[f64]

The flattened cond_bins × gc_bins table (counts[ctx * gc_bins + gc], row-major), for dumping to the aux bias files.

pub fn inc(&mut self, gc_frac: i32, ctx_frac: i32, weight: f64)

Accumulate weight for a fragment with the given GC and context fractions.

Examples found in repository

examples/conv_bench.rs (line 63)

fn main() {
    let args: Vec<String> = std::env::args().collect();
    let n: usize = args.get(1).and_then(|s| s.parse().ok()).unwrap_or(100_000);
    let passes: usize = args.get(2).and_then(|s| s.parse().ok()).unwrap_or(3);

    let mut rng = Lcg(0x9E3779B97F4A7C15);

    // Build a realistic transcript-length distribution: log-normal-ish around
    // ~1.5 kb with a long tail, clamped to [80, 15000] — close to GRCh38 cDNA.
    let bases = [b'A', b'C', b'G', b'T'];
    let mut seqs: Vec<Vec<u8>> = Vec::with_capacity(n);
    let mut prefixes: Vec<Vec<u32>> = Vec::with_capacity(n);
    for _ in 0..n {
        let u = rng.next_f64();
        // exp of a normal-ish variate -> heavy right tail
        let len = (7.3 + 0.9 * (u - 0.5) * 4.0).exp() as usize;
        let len = len.clamp(80, 15_000);
        let seq: Vec<u8> = (0..len)
            .map(|_| bases[(rng.next_u32() & 3) as usize])
            .collect();
        prefixes.push(gc_prefix(&seq));
        seqs.push(seq);
    }

    // A normalized 3×25 GC ratio model (the shape doesn't matter for timing; we
    // exercise the same binning + lookup the real run does).
    let mut obs = GcFragModel::new(DEFAULT_COND_BINS, DEFAULT_GC_BINS);
    let mut exp = GcFragModel::new(DEFAULT_COND_BINS, DEFAULT_GC_BINS);
    for ctx in 0..=100 {
        for gc in 0..=100 {
            obs.inc(gc, ctx, 1.0 + 0.3 * ((gc + ctx) as f64).sin());
            exp.inc(gc, ctx, 1.0);
        }
    }
    let gc_model = gc_ratio(&mut obs, &mut exp, 1000.0);

    // Fragment-length CDF: gaussian-ish pmf around mean 250, sd 40, over 0..1000.
    let fld_max = 1000usize;
    let mean = 250.0f64;
    let sd = 40.0f64;
    let pmf: Vec<f64> = (0..=fld_max)
        .map(|l| {
            let z = (l as f64 - mean) / sd;
            (-0.5 * z * z).exp()
        })
        .collect();
    let (cdf, fld_low, fld_high) = salmon_model::seqbias::fld_cdf_and_bounds(&pmf);

    let mut acc = 0.0f64;
    let mut best = f64::INFINITY;
    for p in 0..passes {
        let t = Instant::now();
        for (seq, prefix) in seqs.iter().zip(&prefixes) {
            let ref_len = seq.len() as f64;
            let elen = (ref_len - 200.0).max(1.0); // ensure unprocessed > 0
            let bias = BiasInputs {
                seq: None,
                gc: Some((&gc_model, salmon_model::GcView::Dense(prefix.as_slice()))),
                pos: None,
            };
            acc += corrected_effective_length_full(
                seq,
                &cdf,
                fld_low,
                fld_high,
                &bias,
                elen,
                GC_SAMP_STRIDE,
                false,
            );
        }
        let dt = t.elapsed().as_secs_f64();
        best = best.min(dt);
        eprintln!(
            "pass {p}: {:.3}s  ({:.2} µs/transcript)  acc={}",
            dt,
            dt * 1e6 / n as f64,
            black_box(acc)
        );
    }
    eprintln!(
        "BEST: {:.3}s over {} transcripts ({:.2} µs/transcript)",
        best,
        n,
        best * 1e6 / n as f64
    );
}

pub fn get(&self, gc_frac: i32, ctx_frac: i32) -> f64

Value at a (context, GC) cell (a normalized density after normalize, or a clamped ratio for a model produced by gc_ratio).

pub fn combine_counts(&mut self, other: &GcFragModel)

Merge another (compatible) model’s counts. Both must be pre-normalization.

pub fn normalize(&mut self)

Normalize each conditioning row into a distribution over GC bins, with a pseudocount prior (salmon’s default 0.1). Idempotent.

Trait Implementations

impl Clone for GcFragModel

fn clone(&self) -> GcFragModel

Returns a duplicate of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl Debug for GcFragModel

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.