crispr_screen 0.1.14

A fast and configurable differential expression analysis tool for CRISPR screens
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

use crate::enrich::EnrichmentResult;
use anyhow::Result;
use ndarray::Array1;
use std::{
    fs::File,
    io::{BufWriter, Write},
};

pub struct SgrnaFrame<'a> {
    sgrna_names: &'a [String],
    gene_names: &'a [String],
    adj_var: &'a Array1<f64>,
    base: &'a Array1<f64>,
    control: &'a Array1<f64>,
    treatment: &'a Array1<f64>,
    pvalue_low: &'a Array1<f64>,
    pvalue_high: &'a Array1<f64>,
    pvalue_twosided: &'a Array1<f64>,
    fdr: &'a Array1<f64>,
    fold_change: &'a Array1<f64>,
    log_fold_change: &'a Array1<f64>,
    size: usize,
}
impl<'a> SgrnaFrame<'a> {
    pub fn new(
        sgrna_names: &'a [String],
        gene_names: &'a [String],
        adj_var: &'a Array1<f64>,
        sgrna_results: &'a EnrichmentResult,
    ) -> Self {
        Self {
            sgrna_names,
            gene_names,
            adj_var,
            base: sgrna_results.base_means(),
            control: sgrna_results.control_means(),
            treatment: sgrna_results.treatment_means(),
            pvalue_low: sgrna_results.pvalues_low(),
            pvalue_high: sgrna_results.pvalues_high(),
            pvalue_twosided: sgrna_results.pvalues_twosided(),
            fdr: sgrna_results.fdr(),
            fold_change: sgrna_results.fold_change(),
            log_fold_change: sgrna_results.log_fold_change(),
            size: sgrna_names.len(),
        }
    }

    pub fn write(&self, prefix: &str) -> Result<()> {
        let mut writer = File::create(format!("{prefix}.sgrna_results.tab")).map(BufWriter::new)?;

        writeln!(
            writer,
            "sgrna\tgene\tbase\tcontrol\ttreatment\tadj_var\tfold_change\tlog2_fold_change\tpvalue_low\tpvalue_high\tpvalue_twosided\tfdr",
        )?;

        for idx in 0..self.size {
            writeln!(
                writer,
                "{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{:e}\t{:e}\t{:e}\t{:e}",
                self.sgrna_names[idx],
                self.gene_names[idx],
                self.base[idx],
                self.control[idx],
                self.treatment[idx],
                self.adj_var[idx],
                self.fold_change[idx],
                self.log_fold_change[idx],
                self.pvalue_low[idx],
                self.pvalue_high[idx],
                self.pvalue_twosided[idx],
                self.fdr[idx]
            )?;
        }

        Ok(())
    }
}

#[cfg(test)]
mod testing {
    use adjustp::Procedure;

    use crate::enrich::EnrichmentResult;

    #[test]
    fn test_sgrna_frame() {
        let sgrna_names = vec!["sgrna1".to_string(), "sgrna2".to_string()];
        let gene_names = vec!["gene1".to_string(), "gene2".to_string()];
        let adj_var = vec![0.1, 0.2].into_iter().collect();
        let control = vec![1.0, 2.0].into_iter().collect();
        let treatment = vec![3.0, 4.0].into_iter().collect();
        let pvalue_low = vec![0.1, 0.2].into_iter().collect();
        let pvalue_high = vec![0.3, 0.4].into_iter().collect();
        let correction = Procedure::BenjaminiHochberg;

        let enrichment =
            EnrichmentResult::new(pvalue_low, pvalue_high, control, treatment, correction);

        let sgrna_frame = super::SgrnaFrame::new(&sgrna_names, &gene_names, &adj_var, &enrichment);
        assert_eq!(sgrna_frame.size, 2);
        sgrna_frame.write("test").unwrap();
    }
}