edge_rs/limma/marray.rs
1//! The fit object the limma pipeline passes along.
2//!
3//! limma's `MArrayLM` is an untyped tagged list whose field set depends on
4//! which functions have run: `lmFit` sets one group, `contrasts.fit` replaces
5//! three, `eBayes` adds eleven more. The field-to-dimension map is in
6//! `R/subsetting.R:115-121`.
7//!
8//! Here it is one struct whose later fields are `Option`, populated in the same
9//! order. Each stage consumes and returns it, so nothing is cloned, and a fit
10//! that has not been through `eBayes` cannot be passed where a moderated t is
11//! needed: the `Option` is the check.
12
13use crate::limma::lm_fit::LmFitResult;
14use crate::prelude::*;
15
16//////////////
17// MArrayLm //
18//////////////
19
20/// A linear model fit, growing as the pipeline runs over it.
21///
22/// Everything before [`MArrayLm::contrasts`] comes from `lm_fit`. After
23/// `contrasts_fit` the coefficient axis is the contrast axis, and
24/// [`MArrayLm::n_coef`] counts contrasts rather than design columns.
25#[derive(Clone, Debug)]
26pub struct MArrayLm {
27 /// Number of genes.
28 pub n_genes: usize,
29 /// Number of samples in the design.
30 pub n_samples: usize,
31 /// Width of the coefficient axis: design columns, or contrasts once
32 /// `contrasts_fit` has run.
33 pub n_coef: usize,
34 /// Coefficients, row-major `n_genes * n_coef`. `NaN` where not estimable.
35 pub coefficients: Vec<f64>,
36 /// Unscaled standard deviations, row-major `n_genes * n_coef`.
37 pub stdev_unscaled: Vec<f64>,
38 /// Residual standard deviation per gene.
39 pub sigma: Vec<f64>,
40 /// Residual degrees of freedom per gene.
41 pub df_residual: Vec<f64>,
42 /// Unscaled covariance of the estimable coefficients, row-major
43 /// `rank * rank`, or `n_coef * n_coef` after `contrasts_fit`.
44 pub cov_coefficients: Vec<f64>,
45 /// Design column indices, accepted ones first. Length is the original
46 /// `n_coef`, which after `contrasts_fit` is no longer this struct's.
47 pub pivot: Vec<usize>,
48 /// Rank of the design.
49 pub rank: usize,
50 /// Design matrix, row-major `n_samples * pivot.len()`.
51 pub design: Vec<f64>,
52 /// Row means of the response, one per gene. limma's `Amean`.
53 ///
54 /// `None` when the caller had none to give, which makes
55 /// `EBayesTrend::Amean` an error rather than a silent no-trend fit.
56 pub amean: Option<Vec<f64>>,
57 /// Contrast matrix, column-major `original n_coef * n_coef`, once
58 /// `contrasts_fit` has run.
59 pub contrasts: Option<Vec<f64>>,
60 /// Prior degrees of freedom. Length one, or `n_genes` for a robust fit.
61 pub df_prior: Option<Vec<f64>>,
62 /// Prior variance. Length one, or `n_genes` for a trended fit.
63 pub s2_prior: Option<Vec<f64>>,
64 /// Prior variance of the coefficients, one per coefficient. The
65 /// B-statistic's `v0`.
66 pub var_prior: Option<Vec<f64>>,
67 /// Assumed proportion of differentially expressed genes.
68 pub proportion: Option<f64>,
69 /// Posterior variance per gene.
70 pub s2_post: Option<Vec<f64>>,
71 /// Moderated t, row-major `n_genes * n_coef`.
72 pub t: Option<Vec<f64>>,
73 /// Total degrees of freedom per gene, residual plus prior, capped.
74 pub df_total: Option<Vec<f64>>,
75 /// Two-sided p-values, row-major `n_genes * n_coef`.
76 pub p_value: Option<Vec<f64>>,
77 /// Log-odds of differential expression, row-major `n_genes * n_coef`.
78 pub lods: Option<Vec<f64>>,
79 /// Moderated F per gene. `None` when the design is not full rank.
80 pub f_stat: Option<Vec<f64>>,
81 /// P-value of [`MArrayLm::f_stat`].
82 pub f_p_value: Option<Vec<f64>>,
83}
84
85impl MArrayLm {
86 /// Wraps an `lm_fit` result, with the response's row means alongside.
87 ///
88 /// ### Params
89 ///
90 /// * `fit` - What `lm_fit` returned
91 /// * `design` - Row-major design, `n_samples * n_coef`
92 /// * `n_coef` - Number of design columns
93 /// * `n_samples` - Number of samples
94 /// * `amean` - Row means of the response, one per gene, or `None`.
95 /// `voom` returns them as [`crate::limma::voom::VoomResult::amean`].
96 ///
97 /// ### Returns
98 ///
99 /// The fit object with every empirical Bayes field unset, or
100 /// [`EdgeErrors::LengthMismatch`] if `design` or `amean` disagrees with the
101 /// fit's shape.
102 pub fn from_lm_fit(
103 fit: LmFitResult,
104 design: &[f64],
105 n_coef: usize,
106 n_samples: usize,
107 amean: Option<Vec<f64>>,
108 ) -> Result<Self, EdgeErrors> {
109 if n_coef == 0 {
110 return Err(EdgeErrors::MustBePositive("n_coef".to_string()));
111 }
112 let n_genes = fit.sigma.len();
113 if fit.coefficients.len() != n_genes * n_coef {
114 return Err(EdgeErrors::LengthMismatch {
115 name: "coefficients",
116 expected: n_genes * n_coef,
117 got: fit.coefficients.len(),
118 });
119 }
120 if design.len() != n_samples * n_coef {
121 return Err(EdgeErrors::LengthMismatch {
122 name: "design",
123 expected: n_samples * n_coef,
124 got: design.len(),
125 });
126 }
127 if let Some(a) = amean.as_ref()
128 && a.len() != n_genes
129 {
130 return Err(EdgeErrors::LengthMismatch {
131 name: "amean",
132 expected: n_genes,
133 got: a.len(),
134 });
135 }
136
137 Ok(Self {
138 n_genes,
139 n_samples,
140 n_coef,
141 coefficients: fit.coefficients,
142 stdev_unscaled: fit.stdev_unscaled,
143 sigma: fit.sigma,
144 df_residual: fit.df_residual,
145 cov_coefficients: fit.cov_coefficients,
146 pivot: fit.pivot,
147 rank: fit.rank,
148 design: design.to_vec(),
149 amean,
150 contrasts: None,
151 df_prior: None,
152 s2_prior: None,
153 var_prior: None,
154 proportion: None,
155 s2_post: None,
156 t: None,
157 df_total: None,
158 p_value: None,
159 lods: None,
160 f_stat: None,
161 f_p_value: None,
162 })
163 }
164
165 /// Clears everything `eBayes` and `treat` set.
166 ///
167 /// `contrasts.fit` drops the test statistics because they belong to the old
168 /// coefficient axis (`R/contrasts.R:21-25`).
169 pub(crate) fn clear_tests(&mut self) {
170 self.df_prior = None;
171 self.s2_prior = None;
172 self.var_prior = None;
173 self.proportion = None;
174 self.s2_post = None;
175 self.t = None;
176 self.df_total = None;
177 self.p_value = None;
178 self.lods = None;
179 self.f_stat = None;
180 self.f_p_value = None;
181 }
182}
183
184///////////
185// Tests //
186///////////
187
188#[cfg(test)]
189mod tests {
190 use super::*;
191
192 fn fit(n_genes: usize, n_coef: usize) -> LmFitResult {
193 LmFitResult {
194 coefficients: vec![0.0; n_genes * n_coef],
195 stdev_unscaled: vec![1.0; n_genes * n_coef],
196 sigma: vec![1.0; n_genes],
197 df_residual: vec![3.0; n_genes],
198 fitted: vec![0.0; n_genes * 5],
199 cov_coefficients: vec![1.0; n_coef * n_coef],
200 pivot: (0..n_coef).collect(),
201 rank: n_coef,
202 }
203 }
204
205 #[test]
206 fn test_from_lm_fit_carries_the_shape() {
207 let m = MArrayLm::from_lm_fit(fit(4, 2), &[1.0; 10], 2, 5, None).unwrap();
208 assert_eq!(m.n_genes, 4);
209 assert_eq!(m.n_coef, 2);
210 assert!(m.t.is_none());
211 }
212
213 #[test]
214 fn test_from_lm_fit_rejects_a_bad_amean() {
215 let e = MArrayLm::from_lm_fit(fit(4, 2), &[1.0; 10], 2, 5, Some(vec![0.0; 3]));
216 assert!(e.is_err());
217 }
218
219 #[test]
220 fn test_from_lm_fit_rejects_a_bad_design() {
221 assert!(MArrayLm::from_lm_fit(fit(4, 2), &[1.0; 8], 2, 5, None).is_err());
222 }
223
224 #[test]
225 fn test_clear_tests_empties_every_ebayes_field() {
226 let mut m = MArrayLm::from_lm_fit(fit(4, 2), &[1.0; 10], 2, 5, None).unwrap();
227 m.t = Some(vec![1.0; 8]);
228 m.f_stat = Some(vec![1.0; 4]);
229 m.clear_tests();
230 assert!(m.t.is_none());
231 assert!(m.f_stat.is_none());
232 }
233}