rsomics_limma_decide_tests/
lib.rs

1//! decideTests: classify each gene × contrast as up (+1) / down (-1) / notsig (0).
2//!
3//! Clean-room reimplementation of limma's decideTests. Method and behaviour
4//! follow the limma documentation and black-box observation; no GPL source was
5//! consulted. Smyth (2004), doi:10.2202/1544-6115.1027.
6
7mod adjust;
8mod table;
9
10use std::io::{BufWriter, Write};
11use std::path::Path;
12
13use rsomics_common::{Result, RsomicsError};
14
15pub use adjust::Adjust;
16pub use table::{Fit, read_fit};
17
18#[derive(Clone, Copy, Debug, PartialEq, Eq)]
19pub enum Method {
20    Separate,
21    Global,
22}
23
24pub struct Options<'a> {
25    pub fit: &'a Path,
26    pub method: Method,
27    pub adjust: Adjust,
28    pub p_value: f64,
29    pub lfc: f64,
30}
31
32pub struct Decisions {
33    pub genes: Vec<String>,
34    pub contrasts: Vec<String>,
35    /// row-major [gene][contrast], each -1 / 0 / 1
36    pub calls: Vec<Vec<i8>>,
37}
38
39pub fn run(opts: &Options) -> Result<Decisions> {
40    let fit = read_fit(opts.fit)?;
41    Ok(decide(
42        &fit,
43        opts.method,
44        opts.adjust,
45        opts.p_value,
46        opts.lfc,
47    ))
48}
49
50pub fn decide(fit: &Fit, method: Method, adjust: Adjust, p_value: f64, lfc: f64) -> Decisions {
51    let ng = fit.genes.len();
52    let nc = fit.contrasts.len();
53
54    let adj = match method {
55        Method::Separate => {
56            let mut a = vec![vec![0.0f64; nc]; ng];
57            let mut col = vec![0.0f64; ng];
58            for j in 0..nc {
59                for (c, row) in col.iter_mut().zip(&fit.p) {
60                    *c = row[j];
61                }
62                for (dst, v) in a.iter_mut().zip(adjust.apply(&col)) {
63                    dst[j] = v;
64                }
65            }
66            a
67        }
68        // limma pools p-values column-major before adjusting.
69        Method::Global => {
70            let mut pooled = Vec::with_capacity(ng * nc);
71            for j in 0..nc {
72                for row in &fit.p {
73                    pooled.push(row[j]);
74                }
75            }
76            let adjusted = adjust.apply(&pooled);
77            let mut a = vec![vec![0.0f64; nc]; ng];
78            let mut k = 0;
79            for j in 0..nc {
80                for row in a.iter_mut() {
81                    row[j] = adjusted[k];
82                    k += 1;
83                }
84            }
85            a
86        }
87    };
88
89    let mut calls = vec![vec![0i8; nc]; ng];
90    for (i, crow) in calls.iter_mut().enumerate() {
91        for (j, call) in crow.iter_mut().enumerate() {
92            let sig = adj[i][j] <= p_value && fit.logfc[i][j].abs() >= lfc;
93            if sig {
94                *call = fit.t[i][j].signum() as i8;
95            }
96        }
97    }
98
99    Decisions {
100        genes: fit.genes.clone(),
101        contrasts: fit.contrasts.clone(),
102        calls,
103    }
104}
105
106pub fn write_decisions(d: &Decisions, out: &mut dyn Write) -> Result<()> {
107    let mut w = BufWriter::with_capacity(1 << 20, out);
108    w.write_all(b"gene").map_err(RsomicsError::Io)?;
109    for c in &d.contrasts {
110        write!(w, "\t{c}").map_err(RsomicsError::Io)?;
111    }
112    w.write_all(b"\n").map_err(RsomicsError::Io)?;
113
114    let mut line = String::with_capacity(64);
115    for (gene, row) in d.genes.iter().zip(&d.calls) {
116        line.clear();
117        line.push_str(gene);
118        for &v in row {
119            line.push('\t');
120            match v {
121                1 => line.push('1'),
122                -1 => line.push_str("-1"),
123                _ => line.push('0'),
124            }
125        }
126        line.push('\n');
127        w.write_all(line.as_bytes()).map_err(RsomicsError::Io)?;
128    }
129    w.flush().map_err(RsomicsError::Io)?;
130    Ok(())
131}
rsomics_limma_decide_tests/lib.rs

rsomics_limma_decide_tests/
lib.rs
