Skip to main content

bixverse_methods/methods/
diffcor.rs

1//! Implementation of differential correlation method based on the Fisher method
2
3use faer::Mat;
4use rayon::prelude::*;
5
6use crate::prelude::*;
7use bixverse_core::core::math::stats::z_scores_to_pval;
8
9////////////////
10// Structures //
11////////////////
12
13/// Structure for DiffCor results
14#[derive(Clone, Debug)]
15pub struct DiffCorRes<T> {
16    /// Correlation coefficients of a
17    pub r_a: Vec<T>,
18    /// Correlation coefficients of b
19    pub r_b: Vec<T>,
20    /// Z-scores of the differential correlation
21    pub z_score: Vec<T>,
22    /// Calculated p-values from the Z-scores
23    pub p_vals: Vec<T>,
24}
25
26//////////
27// Main //
28//////////
29
30/// Calculate differential correlations
31///
32/// The function will panic if the two correlation matrices are not symmetric
33/// and do not have the same dimensions.
34///
35/// ### Params
36///
37/// * `mat_a` - The first correlation matrix.
38/// * `mat_b` - The second correlation matrix.
39/// * `no_sample_a` - Number of samples that were present to calculate mat_a.
40/// * `no_sample_b` - Number of samples that were present to calculate mat_b.
41/// * `spearman` - Was Spearman correlation used.
42///
43/// ### Returns
44///
45/// The resulting differential correlation results as a structure.
46pub fn calculate_diff_correlation<T>(
47    mat_a: &Mat<T>,
48    mat_b: &Mat<T>,
49    no_sample_a: usize,
50    no_sample_b: usize,
51    spearman: bool,
52) -> DiffCorRes<T>
53where
54    T: BixverseFloat,
55{
56    assert_symmetric_mat!(mat_a);
57    assert_symmetric_mat!(mat_b);
58    // assert_same_dims!(mat_a, mat_b);
59
60    let mut cors_a: Vec<T> = Vec::new();
61    let mut cors_b: Vec<T> = Vec::new();
62
63    let upper_triangle_indices = upper_triangle_indices(mat_a.ncols(), 1);
64
65    for (&r, &c) in upper_triangle_indices
66        .0
67        .iter()
68        .zip(upper_triangle_indices.1.iter())
69    {
70        cors_a.push(*mat_a.get(r, c));
71        cors_b.push(*mat_b.get(r, c));
72    }
73
74    // Maybe save the original correlations... Note to myself.
75    let original_cor_a = cors_a.to_vec();
76    let original_cor_b = cors_b.to_vec();
77
78    cors_a.par_iter_mut().for_each(|x| *x = x.atanh());
79    cors_b.par_iter_mut().for_each(|x| *x = x.atanh());
80
81    // Constant will depend on if Spearman or Pearson
82    let constant = if spearman {
83        T::from_f32(1.06).unwrap()
84    } else {
85        T::from_f32(1.0).unwrap()
86    };
87    let denominator: T = ((constant
88        / (T::from_usize(no_sample_a).unwrap() - T::from_f32(3.0).unwrap()))
89        + (constant / (T::from_usize(no_sample_b).unwrap() - T::from_f32(3.0).unwrap())))
90    .sqrt();
91
92    let z_scores: Vec<T> = cors_a
93        .par_iter()
94        .zip(cors_b.par_iter())
95        .map(|(a, b)| (*a - *b) / denominator)
96        .collect();
97
98    let p_values = z_scores_to_pval(&z_scores, "twosided");
99
100    DiffCorRes {
101        r_a: original_cor_a,
102        r_b: original_cor_b,
103        z_score: z_scores,
104        p_vals: p_values,
105    }
106}