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
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94

use tracing::{debug, info};

use crate::input::error::InputError;
use crate::model::timecourse::{SplicingTrajectoryClass, TimecourseSplicingMetrics};

pub fn run_stage14(
    timepoints: Option<&[i64]>,
    sis_median: &[f32],
    entropy_median: &[f32],
    imbalance_median: &[f32],
) -> Result<Option<TimecourseSplicingMetrics>, InputError> {
    let Some(timepoints) = timepoints else {
        return Ok(None);
    };
    let metrics = compute(timepoints, sis_median, entropy_median, imbalance_median)?;
    Ok(Some(metrics))
}

pub fn compute(
    timepoints: &[i64],
    sis_median: &[f32],
    entropy_median: &[f32],
    imbalance_median: &[f32],
) -> Result<TimecourseSplicingMetrics, InputError> {
    let n = timepoints.len();
    if n == 0 {
        return Err(InputError::MissingTimecourseSignals);
    }
    if sis_median.len() != n || entropy_median.len() != n || imbalance_median.len() != n {
        return Err(InputError::LengthMismatch(
            "timecourse median length mismatch".to_string(),
        ));
    }

    let delta_sis = deltas(sis_median);
    let delta_entropy = deltas(entropy_median);
    let delta_imbalance = deltas(imbalance_median);

    let mut trajectory = SplicingTrajectoryClass::Inconclusive;
    let has_nan = sis_median
        .iter()
        .chain(entropy_median.iter())
        .chain(imbalance_median.iter())
        .any(|v| !v.is_finite())
        || delta_sis
            .iter()
            .chain(delta_entropy.iter())
            .chain(delta_imbalance.iter())
            .any(|v| !v.is_finite());

    if n >= 3 && !has_nan {
        let adaptive = delta_sis.iter().all(|d| *d > 0.0)
            && delta_entropy.iter().all(|d| *d <= 0.0)
            && delta_imbalance.iter().all(|d| *d <= 0.0);
        let degenerative = delta_sis.iter().all(|d| *d < 0.0)
            && delta_entropy.iter().all(|d| *d >= 0.0)
            && delta_imbalance.iter().all(|d| *d >= 0.0);

        trajectory = if adaptive {
            SplicingTrajectoryClass::Adaptive
        } else if degenerative {
            SplicingTrajectoryClass::Degenerative
        } else {
            SplicingTrajectoryClass::Oscillatory
        };
    }

    info!("time-course coherence analysis completed");
    debug!(
        trajectory = ?trajectory,
        delta_sis = ?delta_sis,
        delta_entropy = ?delta_entropy,
        delta_imbalance = ?delta_imbalance,
        "time-course deltas"
    );

    Ok(TimecourseSplicingMetrics {
        trajectory,
        delta_sis,
        delta_entropy,
        delta_imbalance,
    })
}

fn deltas(values: &[f32]) -> Vec<f32> {
    if values.len() < 2 {
        return Vec::new();
    }
    let mut out = Vec::with_capacity(values.len() - 1);
    for idx in 0..(values.len() - 1) {
        out.push(values[idx + 1] - values[idx]);
    }
    out
}