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

use crate::{
    aggregation::compute_aggregation,
    enrich::enrichment_testing,
    io::{GeneFrame, SgrnaFrame, SimpleFrame},
    model::model_mean_variance,
    norm::normalize_counts,
    utils::{config::Configuration, logging::Logger},
};
use anyhow::Result;

/// Performs the `MAGeCK` Differential Expression and Gene Aggregation Algorithm
pub fn mageck(
    frame: &SimpleFrame,
    labels_controls: &[String],
    labels_treatments: &[String],
    config: &Configuration,
    logger: &Logger,
) -> Result<()> {
    let labels = [labels_controls, labels_treatments].concat();
    let count_matrix = frame.data_matrix(&labels)?;
    let sgrna_names = frame.get_sgrna_names();
    let gene_names = frame.get_gene_names();
    let n_controls = labels_controls.len();

    logger.start_mageck();
    logger.num_sgrnas(sgrna_names);
    logger.num_genes(gene_names);
    logger.norm_method(config.normalization());
    logger.aggregation_method(config.aggregation());
    logger.correction(config.correction());

    // Normalize
    let normed_matrix = normalize_counts(&count_matrix, config.normalization(), logger);

    // Mean-Variance Modeling
    let adj_var = model_mean_variance(&normed_matrix, n_controls, config.model_choice(), logger);

    // sgRNA Ranking (Enrichment)
    let sgrna_results =
        enrichment_testing(&normed_matrix, &adj_var, n_controls, config.correction());

    // Gene Ranking (Aggregation)
    let aggregation_results = compute_aggregation(
        config.aggregation(),
        &sgrna_results,
        gene_names,
        logger,
        config.correction(),
    );

    // Build sgRNA DataFrame
    let sgrna_frame = SgrnaFrame::new(sgrna_names, gene_names, &adj_var, &sgrna_results);
    sgrna_frame.write(config.prefix())?;

    // Build Gene DataFrame
    let gene_frame = GeneFrame::new(&aggregation_results);
    gene_frame.write(config.prefix())?;

    Ok(())
}