bixverse_methods/methods/
diffcor.rs1use faer::Mat;
4use rayon::prelude::*;
5
6use crate::prelude::*;
7use bixverse_core::core::math::stats::z_scores_to_pval;
8
9#[derive(Clone, Debug)]
15pub struct DiffCorRes<T> {
16 pub r_a: Vec<T>,
18 pub r_b: Vec<T>,
20 pub z_score: Vec<T>,
22 pub p_vals: Vec<T>,
24}
25
26pub 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 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 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 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}