rsomics-limma-array-weights 0.1.0

//! limma arrayWeights: per-array quality weights estimated by REML.
//!
//! Clean-room reimplementation of limma's `arrayWeights`. The method is Ritchie et al. (2006),
//! BMC Bioinformatics 7:261, doi:10.1186/1471-2105-7-261, with the REML scoring adapted from
//! Smyth (2002), J. Comput. Graph. Statist. 11:836. No limma (GPL) source was consulted; the
//! implementation follows the published method and is validated black-box against the binary.

mod matrix;
mod reml;

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

use rsomics_common::{Result, RsomicsError};

pub use matrix::{read_design, read_expr};

pub struct Options<'a> {
    pub expr: &'a Path,
    pub design: &'a Path,
    pub prior_n: f64,
    pub maxiter: usize,
    pub tol: f64,
}

pub struct Results {
    pub samples: Vec<String>,
    pub weights: Vec<f64>,
    pub n_genes: usize,
    pub iters: usize,
    pub converged: bool,
}

pub fn run(opts: &Options) -> Result<Results> {
    let expr = read_expr(opts.expr)?;
    let design = read_design(opts.design)?;
    if design.x.len() != expr.samples.len() {
        return Err(RsomicsError::InvalidInput(format!(
            "design has {} rows, expression has {} samples",
            design.x.len(),
            expr.samples.len()
        )));
    }

    let fit = reml::array_weights(&expr.y, &design.x, opts.prior_n, opts.maxiter, opts.tol)?;

    Ok(Results {
        samples: expr.samples,
        weights: fit.weights,
        n_genes: expr.y.len(),
        iters: fit.iters,
        converged: fit.converged,
    })
}

pub fn write_weights(res: &Results, out: &mut dyn Write) -> Result<()> {
    let mut buf = ryu::Buffer::new();
    writeln!(out, "sample\tweight").map_err(RsomicsError::Io)?;
    for (s, w) in res.samples.iter().zip(&res.weights) {
        writeln!(out, "{s}\t{}", buf.format(*w)).map_err(RsomicsError::Io)?;
    }
    Ok(())
}