Hi Wei and team,
I've been reimplementing SAIGE in Rust and noticed what appears to be a bug in the fast score test for common variants.
When --is_fastTest=TRUE (the default), all variants with AF > 0.5 show massively inflated variance estimates, which deflates effect sizes and p-values.
For a variant with AF = 0.9515:
- is_fastTest=TRUE: var=717.7, BETA=-0.018, p=0.64
- is_fastTest=FALSE: var=20.1, BETA=-0.626, p=0.005
That's a 35× difference in variance..
I think the root cause is this: In scoreTestFast_noadjCov, when AF > 0.5, genotypes are flipped to the minor allele (g_flip = 2 - g), but the allele frequency used for mean-centering isn't updated? The code computes (g_flip - 2*AF)^2 instead of (g_flip - 2*(1-AF))^2. Since variance is symmetric under allele flipping (2pq = 2qp), I believe this is unintentional.
I verified this by computing both formulas across 2,000 variants.. the "buggy" formula matches SAIGE's output exactly.
Workaround
As noted above, setting --is_fastTest=FALSE produces results that correlate r > 0.99 with my implementation. I tested this on SAIGE v1.3.0 and v1.5.0; same behaviuor in both.
Happy to provide more details, but i think the fix could simply be (scoreTestFast_noadjCov) to use MAF for centring:
Before fix:
double var2_a = dot(m_mu21,pow(g1,2));
double var2_b = dot(m_mu21, 2*2*t_altFreq*g1);
double var2_c = arma::accu(m_mu2)*pow(2*t_altFreq, 2);
After fix:
double maf = (t_altFreq > 0.5) ? (1.0 - t_altFreq) : t_altFreq;
double var2_a = dot(m_mu21,pow(g1,2));
double var2_b = dot(m_mu21, 2*2*maf*g1);
double var2_c = arma::accu(m_mu2)*pow(2*maf, 2);
Thanks,
Fred
Hi Wei and team,
I've been reimplementing SAIGE in Rust and noticed what appears to be a bug in the fast score test for common variants.
When
--is_fastTest=TRUE(the default), all variants with AF > 0.5 show massively inflated variance estimates, which deflates effect sizes and p-values.For a variant with AF = 0.9515:
That's a 35× difference in variance..
I think the root cause is this: In
scoreTestFast_noadjCov, when AF > 0.5, genotypes are flipped to the minor allele (g_flip = 2 - g), but the allele frequency used for mean-centering isn't updated? The code computes(g_flip - 2*AF)^2instead of(g_flip - 2*(1-AF))^2. Since variance is symmetric under allele flipping (2pq = 2qp), I believe this is unintentional.I verified this by computing both formulas across 2,000 variants.. the "buggy" formula matches SAIGE's output exactly.
Workaround
As noted above, setting --is_fastTest=FALSE produces results that correlate r > 0.99 with my implementation. I tested this on SAIGE v1.3.0 and v1.5.0; same behaviuor in both.
Happy to provide more details, but i think the fix could simply be (scoreTestFast_noadjCov) to use MAF for centring:
Before fix:
After fix:
Thanks,
Fred