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