kira-spliceqc 0.2.0

Deterministic, explainable splicing QC for single-cell expression data.
Documentation 
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

use std::time::Instant;

use rayon::prelude::*;
use tracing::{debug, info, warn};

use crate::expression::ExpressionMatrix;
use crate::genesets::catalog::default_catalog_path;
use crate::genesets::{GenesetCatalog, load_catalog};
use crate::input::error::InputError;
use crate::model::geneset_activity::GenesetActivityMatrix;

pub fn run_stage2(matrix: &dyn ExpressionMatrix) -> Result<GenesetActivityMatrix, InputError> {
    let catalog_path = default_catalog_path();
    let catalog = load_catalog(&catalog_path, matrix)?;
    aggregate(matrix, &catalog)
}

pub fn aggregate(
    matrix: &dyn ExpressionMatrix,
    catalog: &GenesetCatalog,
) -> Result<GenesetActivityMatrix, InputError> {
    let n_cells = matrix.n_cells();
    let n_genesets = catalog.genesets.len();
    let mut values = vec![0.0f32; n_genesets * n_cells];

    let genesets: Vec<String> = catalog.genesets.iter().map(|g| g.id.clone()).collect();
    let axes: Vec<String> = catalog.genesets.iter().map(|g| g.axis.clone()).collect();

    let libsize_scale: Vec<f32> = (0..n_cells)
        .map(|cell| 1e4_f32 / matrix.libsize(cell).max(1) as f32)
        .collect();

    let start = Instant::now();
    for (gs_idx, geneset) in catalog.genesets.iter().enumerate() {
        let row_start = gs_idx * n_cells;
        let row = &mut values[row_start..row_start + n_cells];

        if geneset.gene_ids.is_empty() {
            row.fill(f32::NAN);
            warn!(
                geneset_id = geneset.id.as_str(),
                "geneset has no resolved genes"
            );
            continue;
        }

        let panel = geneset.gene_ids.as_slice();
        let n_panel = panel.len() as f32;
        row.par_iter_mut().enumerate().for_each(|(cell, out)| {
            let scale = libsize_scale[cell];
            let sum = matrix.panel_ln1p_scaled_sum(panel, cell, scale);
            *out = sum / n_panel;
        });

        debug!(
            geneset_id = geneset.id.as_str(),
            resolved = geneset.gene_ids.len(),
            "geneset aggregated"
        );
    }

    info!(
        elapsed_ms = start.elapsed().as_millis(),
        genesets = n_genesets,
        "geneset aggregation complete"
    );

    Ok(GenesetActivityMatrix {
        genesets,
        axes,
        values,
        n_cells,
    })
}