limma/
detectionpvalues.rs

1//! Detection p-values from negative controls (limma `detectionPValues`).
2//!
3//! Pure-Rust port of `detectionPValues.default`: for each array, the detection
4//! p-value of a probe is the right-tail rank of its intensity within the
5//! negative-control distribution. Counting controls `>=` and `>` the observed
6//! value and averaging makes ties count as one half.
7
8use ndarray::Array2;
9
10/// Detection p-values for an `nprobes × narrays` intensity matrix `x`, given a
11/// per-probe `status` and the label `negctrl` marking negative controls.
12///
13/// For probe `i` in array `j` the p-value is `(c_ge + c_gt) / (2·nneg)`, where
14/// `c_ge` / `c_gt` count negative controls with intensity `>=` / `>` the
15/// observed value. Requires at least one negative control and finite
16/// intensities.
17pub fn detection_p_values<S: AsRef<str>>(
18    x: &Array2<f64>,
19    status: &[S],
20    negctrl: &str,
21) -> Array2<f64> {
22    let (nrow, ncol) = (x.nrows(), x.ncols());
23    let isneg: Vec<i64> = status
24        .iter()
25        .map(|s| (s.as_ref() == negctrl) as i64)
26        .collect();
27    let nneg: i64 = isneg.iter().sum();
28    let denom = (2 * nneg) as f64;
29
30    let mut out = Array2::<f64>::zeros((nrow, ncol));
31    for j in 0..ncol {
32        let col: Vec<f64> = (0..nrow).map(|i| x[[i, j]]).collect();
33        // o1: value desc, negatives-first on ties; o2: negatives-last on ties.
34        // Final tie-break is ascending index, matching R's stable `order`.
35        let mut o1: Vec<usize> = (0..nrow).collect();
36        o1.sort_by(|&a, &b| {
37            col[b]
38                .partial_cmp(&col[a])
39                .unwrap()
40                .then(isneg[b].cmp(&isneg[a]))
41                .then(a.cmp(&b))
42        });
43        let mut o2: Vec<usize> = (0..nrow).collect();
44        o2.sort_by(|&a, &b| {
45            col[b]
46                .partial_cmp(&col[a])
47                .unwrap()
48                .then(isneg[a].cmp(&isneg[b]))
49                .then(a.cmp(&b))
50        });
51        // cs1[i] = #negatives >= x[i,j] (ties included),
52        // cs2[i] = #negatives >  x[i,j] (ties excluded).
53        let mut run = 0i64;
54        let mut cs = vec![0i64; nrow];
55        for &i in &o1 {
56            run += isneg[i];
57            cs[i] = run;
58        }
59        run = 0;
60        for &i in &o2 {
61            run += isneg[i];
62            cs[i] += run;
63        }
64        for i in 0..nrow {
65            out[[i, j]] = cs[i] as f64 / denom;
66        }
67    }
68    out
69}
70
71#[cfg(test)]
72mod tests {
73    use super::*;
74    use ndarray::array;
75
76    fn close(a: &Array2<f64>, b: &Array2<f64>, tol: f64) -> bool {
77        a.shape() == b.shape() && a.iter().zip(b.iter()).all(|(&x, &y)| (x - y).abs() <= tol)
78    }
79
80    #[test]
81    fn detection_with_ties_matches_r() {
82        let x = array![
83            [10.0, 12.0, 8.0],
84            [10.0, 5.0, 8.0],
85            [3.0, 20.0, 15.0],
86            [7.0, 7.0, 7.0],
87            [3.0, 4.0, 2.0],
88            [1.0, 2.0, 9.0],
89        ];
90        let status = [
91            "regular", "regular", "regular", "regular", "negative", "negative",
92        ];
93        let got = detection_p_values(&x, &status, "negative");
94        let want = array![
95            [0.0, 0.0, 0.5],
96            [0.0, 0.0, 0.5],
97            [0.25, 0.0, 0.0],
98            [0.0, 0.0, 0.5],
99            [0.5, 0.5, 1.0],
100            [1.0, 1.0, 0.5],
101        ];
102        assert!(close(&got, &want, 1e-12));
103    }
104
105    #[test]
106    fn detection_single_negative_matches_r() {
107        let x = array![[9.5, 1.0], [2.0, 8.0], [6.0, 3.0], [0.5, 0.5]];
108        let status = ["reg", "reg", "reg", "neg"];
109        let got = detection_p_values(&x, &status, "neg");
110        let want = array![[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [1.0, 1.0]];
111        assert!(close(&got, &want, 1e-12));
112    }
113}
limma/detectionpvalues.rs

limma/
detectionpvalues.rs
