gsem 0.1.3

Genomic Structural Equation Modeling from GWAS summary statistics
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
113
114
115
116
117
118
119
120
121
122
123

use std::str::FromStr;

use faer::Mat;

/// Genomic control correction mode.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum GcMode {
    Conservative,
    Standard,
    None,
}

impl GcMode {
    /// GC scaling factor for a diagonal element of the LD intercept matrix.
    fn scale(self, i_ld_diag: f64) -> f64 {
        match self {
            GcMode::Conservative => i_ld_diag,
            GcMode::Standard => i_ld_diag.sqrt(),
            GcMode::None => 1.0,
        }
    }
}

impl FromStr for GcMode {
    type Err = std::convert::Infallible;

    fn from_str(s: &str) -> Result<Self, Self::Err> {
        Ok(match s.to_lowercase().as_str() {
            "conserv" | "conservative" => GcMode::Conservative,
            "none" | "off" => GcMode::None,
            _ => GcMode::Standard,
        })
    }
}

/// Construct V_SNP matrix (k×k) for a single SNP with genomic control correction.
///
/// Port of `.get_V_SNP()` from GenomicSEM's utils.R.
pub fn build_v_snp(
    se_snp: &[f64],
    i_ld: &Mat<f64>,
    var_snp: f64,
    gc: GcMode,
    k: usize,
) -> Mat<f64> {
    let mut v_snp = Mat::zeros(k, k);

    for x in 0..k {
        for y in x..k {
            let val = if x == y {
                (se_snp[x] * var_snp * gc.scale(i_ld[(x, x)])).powi(2)
            } else {
                let base = se_snp[x] * se_snp[y] * var_snp.powi(2);
                base * i_ld[(x, y)] * gc.scale(i_ld[(x, x)]) * gc.scale(i_ld[(y, y)])
            };
            v_snp[(x, y)] = val;
            v_snp[(y, x)] = val;
        }
    }

    v_snp
}

/// Compute GC-adjusted Z-scores for a single SNP.
///
/// Port of `.get_Z_pre()` from GenomicSEM's utils.R.
pub fn gc_adjusted_z(beta: &[f64], se: &[f64], i_ld: &Mat<f64>, gc: GcMode, k: usize) -> Vec<f64> {
    (0..k)
        .map(|x| {
            let denom = se[x] * gc.scale(i_ld[(x, x)]);
            if denom.abs() > 1e-30 {
                beta[x] / denom
            } else {
                0.0
            }
        })
        .collect()
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_gc_none_diagonal() {
        let se = vec![0.1, 0.2];
        let i_ld = faer::mat![[1.05, 0.02], [0.02, 1.03]];
        let v = build_v_snp(&se, &i_ld, 0.5, GcMode::None, 2);
        // Diagonal: (se * var_snp)^2 = (0.1 * 0.5)^2 = 0.0025
        assert!((v[(0, 0)] - 0.0025).abs() < 1e-10);
    }

    #[test]
    fn test_gc_conservative_inflates() {
        let se = vec![0.1];
        let i_ld = faer::mat![[1.1]];
        let v_none = build_v_snp(&se, &i_ld, 0.5, GcMode::None, 1);
        let v_conserv = build_v_snp(&se, &i_ld, 0.5, GcMode::Conservative, 1);
        // Conservative should give larger variance
        assert!(v_conserv[(0, 0)] > v_none[(0, 0)]);
    }

    #[test]
    fn test_gc_adjusted_z() {
        let beta = vec![0.5];
        let se = vec![0.1];
        let i_ld = faer::mat![[1.0]];
        let z = gc_adjusted_z(&beta, &se, &i_ld, GcMode::None, 1);
        assert!((z[0] - 5.0).abs() < 1e-10);
    }

    #[test]
    fn test_parse_gc_mode() {
        assert_eq!(
            "conservative".parse::<GcMode>().unwrap(),
            GcMode::Conservative
        );
        assert_eq!("conserv".parse::<GcMode>().unwrap(), GcMode::Conservative);
        assert_eq!("standard".parse::<GcMode>().unwrap(), GcMode::Standard);
        assert_eq!("none".parse::<GcMode>().unwrap(), GcMode::None);
        assert_eq!("unknown".parse::<GcMode>().unwrap(), GcMode::Standard);
    }
}