ggca 1.0.0

Computes efficiently the correlation (Pearson, Spearman or Kendall) and the p-value (two-sided) between all the pairs from two datasets
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
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112

use pyo3::prelude::*;

pub trait Adjustment {
    fn adjust(&mut self, p_value: f64, rank: usize) -> f64;
}

pub struct Bonferroni {
    total_number_of_elements: f64,
}

impl Bonferroni {
    fn new(total_number_of_elements: f64) -> Self {
        Bonferroni {
            total_number_of_elements,
        }
    }
}

impl Adjustment for Bonferroni {
    fn adjust(&mut self, p_value: f64, _: usize) -> f64 {
        (p_value * self.total_number_of_elements).min(1.0)
    }
}

pub struct BenjaminiHochberg {
    total_number_of_elements: f64, // To prevent casting on every adjustment
    previous_max_p_value: f64,
}

impl BenjaminiHochberg {
    fn new(total_number_of_elements: f64) -> Self {
        BenjaminiHochberg {
            total_number_of_elements,
            previous_max_p_value: f64::MAX,
        }
    }
}

impl Adjustment for BenjaminiHochberg {
    fn adjust(&mut self, p_value: f64, rank: usize) -> f64 {
        let valid_rank = self.total_number_of_elements - rank as f64;
        let q_value = p_value * (self.total_number_of_elements / valid_rank);
        let q_value = q_value.min(self.previous_max_p_value).min(1.0);
        self.previous_max_p_value = q_value;
        q_value
    }
}

pub struct BenjaminiYekutieli {
    total_number_of_elements: f64, // To prevent casting on every adjustment
    previous_max_p_value: f64,
    accumulator: f64,
}

impl BenjaminiYekutieli {
    fn new(total_number_of_elements: f64) -> Self {
        let accumulator =
            (1..total_number_of_elements as u64 + 1).fold(0.0, |acc, x| acc + (1.0 / x as f64));

        BenjaminiYekutieli {
            total_number_of_elements,
            previous_max_p_value: f64::MAX,
            accumulator,
        }
    }
}

impl Adjustment for BenjaminiYekutieli {
    fn adjust(&mut self, p_value: f64, rank: usize) -> f64 {
        let valid_rank = self.total_number_of_elements - rank as f64;
        let adjustment_factor = self.accumulator * self.total_number_of_elements / valid_rank;
        let q_value = p_value * adjustment_factor;
        let q_value = q_value.min(self.previous_max_p_value).min(1.0);
        self.previous_max_p_value = q_value;
        q_value
    }
}

#[pyclass(eq, eq_int)]
#[derive(Clone, Debug, PartialEq)]
pub enum AdjustmentMethod {
    BenjaminiHochberg = 1,
    BenjaminiYekutieli = 2,
    Bonferroni = 3,
}

impl std::fmt::Display for AdjustmentMethod {
    fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result {
        let description = match &self {
            AdjustmentMethod::BenjaminiHochberg => "Benjamini-Hochberg",
            AdjustmentMethod::BenjaminiYekutieli => "Benjamini-Yekutieli",
            AdjustmentMethod::Bonferroni => "Bonferroni",
        };

        write!(f, "{description}")
    }
}

pub fn get_adjustment_method(
    adjustment_method: &AdjustmentMethod,
    total_number_of_elements: f64,
) -> Box<dyn Adjustment> {
    match adjustment_method {
        AdjustmentMethod::Bonferroni => Box::new(Bonferroni::new(total_number_of_elements)),
        AdjustmentMethod::BenjaminiHochberg => {
            Box::new(BenjaminiHochberg::new(total_number_of_elements))
        }
        AdjustmentMethod::BenjaminiYekutieli => {
            Box::new(BenjaminiYekutieli::new(total_number_of_elements))
        }
    }
}