Skip to content

Commit 91481ab

Browse files
weizhouUMICHweizhouUMICH
committed
save Sigma_iXXSigma_iX when sparse GRM was used for fitting the null model and no variance ratio was used for step 2
1 parent 8e78d18 commit 91481ab

File tree

6 files changed

+71
-16
lines changed

6 files changed

+71
-16
lines changed

‎R/SAIGE_Test_main.R‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -305,6 +305,7 @@ SPAGMMATtest = function(bgenFile = "",
305305
t_XXVX_inv=obj.model$obj.noK$XXVX_inv,
306306
t_XV=obj.model$obj.noK$XV,
307307
t_XVX_inv_XV=obj.model$obj.noK$XVX_inv_XV,
308+
t_Sigma_iXXSigma_iX=obj.model$Sigma_iXXSigma_iX,
308309
t_X=obj.model$X,
309310
t_S_a=obj.model$obj.noK$S_a,
310311
t_res=obj.model$residuals,

‎R/SAIGE_fitGLMM_fast.R‎

Lines changed: 37 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1144,8 +1144,22 @@ fitNULLGLMM = function(plinkFile = "",
11441144
modglmm$offset = covoffset
11451145
}
11461146
}
1147-
1148-
modglmm$useSparseGRMforVarRatio = useSparseGRMforVarRatio
1147+
1148+
if(skipVarianceRatioEstimation & useSparseGRMtoFitNULL){
1149+
family = fit0$family
1150+
eta = modglmm$linear.predictors
1151+
mu = modglmm$fitted.values
1152+
mu.eta = family$mu.eta(eta)
1153+
sqrtW = mu.eta/sqrt(family$variance(mu))
1154+
W = sqrtW^2
1155+
tauVecNew = modglmm$theta
1156+
Sigma_iX = getSigma_X(W, tauVecNew, modglmm$X, maxiterPCG, tolPCG)
1157+
Sigma_iXXSigma_iX = Sigma_iX%*%(solve(t(modglmm$X)%*%Sigma_iX))
1158+
modglmm$Sigma_iXXSigma_iX = Sigma_iXXSigma_iX
1159+
}
1160+
1161+
1162+
modglmm$useSparseGRMtoFitNULL = useSparseGRMtoFitNULL
11491163
save(modglmm, file = modelOut)
11501164
tau = modglmm$theta
11511165

@@ -1267,7 +1281,23 @@ fitNULLGLMM = function(plinkFile = "",
12671281
}
12681282

12691283
modglmm$offset = covoffset
1270-
modglmm$useSparseGRMforVarRatio = useSparseGRMforVarRatio
1284+
1285+
if(skipVarianceRatioEstimation & useSparseGRMtoFitNULL){
1286+
family = fit0$family
1287+
eta = modglmm$linear.predictors
1288+
mu = modglmm$fitted.values
1289+
mu.eta = family$mu.eta(eta)
1290+
sqrtW = mu.eta/sqrt(family$variance(mu))
1291+
W = sqrtW^2
1292+
tauVecNew = modglmm$theta
1293+
Sigma_iX = getSigma_X(W, tauVecNew, modglmm$X, maxiterPCG, tolPCG)
1294+
Sigma_iXXSigma_iX = Sigma_iX%*%(solve(t(modglmm$X)%*%Sigma_iX))
1295+
modglmm$Sigma_iXXSigma_iX = Sigma_iXXSigma_iX
1296+
}
1297+
1298+
1299+
modglmm$useSparseGRMtoFitNULL = useSparseGRMtoFitNULL
1300+
12711301
save(modglmm, file = modelOut)
12721302
t_end = proc.time()
12731303
print(t_end)
@@ -1774,8 +1804,8 @@ getSparseSigma = function(bedFile,
17741804
#Matrix:::writeMM(sparseGRM, sparseGRMFile)
17751805
Nval = length(W)
17761806
sparseSigma = sparseGRM * tauVecNew[2]
1777-
diag(sparseSigma) = getDiagOfSigma(W, tauVecNew)
1778-
1807+
#diag(sparseSigma) = getDiagOfSigma(W, tauVecNew)
1808+
diag(sparseSigma) = diag(sparseSigma) + (1/W)*tauVecNew[1]
17791809
#sparseSigmaFile = paste0(outputPrefix, "_relatednessCutoff_",relatednessCutoff, "_", numRandomMarkerforSparseKin, "_randomMarkersUsed.sparseSigma.mtx")
17801810
#cat("write sparse Sigma to ", sparseSigmaFile ,"\n")
17811811
#Matrix:::writeMM(sparseSigma, sparseSigmaFile)
@@ -2078,6 +2108,8 @@ extractVarianceRatio = function(obj.glmm.null,
20782108
pcginvSigma = solve(sparseSigma, g, sparse=T)
20792109
var2_a = t(g) %*% pcginvSigma
20802110
var2sparseGRM = var2_a[1,1]
2111+
x=t(G)%*%Sigma_iG/AC
2112+
cat(" x ", x, " var2 ", var2sparseGRM , "\n")
20812113
varRatio_sparseGRM_vec = c(varRatio_sparseGRM_vec, var1/var2sparseGRM)
20822114
}
20832115

‎R/readInGLMM.R‎

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -158,8 +158,11 @@ ReadModel = function(GMMATmodelFile = "", chrom="", LOCO=TRUE, is_Firth_beta=FAL
158158
if(is.null(obj.glmm.null$offset)){
159159
obj.glmm.null$offset = rep(0,nrow(obj.glmm.null$X))
160160
}
161-
162-
#}
161+
162+
163+
if(is.null(obj.glmm.null$Sigma_iXXSigma_iX)){
164+
obj.glmm.null$Sigma_iXXSigma_iX = matrix(0, nrow=1, ncol=1)
165+
}
163166
return(obj.glmm.null)
164167
}
165168

‎src/Main.cpp‎

Lines changed: 9 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -626,6 +626,7 @@ void Unified_getMarkerPval(
626626
if(t_isOnlyOutputNonZero == true)
627627
Rcpp::stop("When using SAIGE method to calculate marker-level p-values, 't_isOnlyOutputNonZero' should be false.");
628628

629+
629630
ptr_gSAIGEobj->getMarkerPval(t_GVec, t_indexForNonZero_vec, t_indexForZero_vec, t_Beta, t_seBeta, t_pval, t_pval_noSPA, t_altFreq, t_Tstat, t_gy, t_varT, t_isSPAConverge, t_gtilde, is_gtilde, is_region, t_P2Vec, t_isCondition, t_Beta_c, t_seBeta_c, t_pval_c, t_pval_noSPA_c, t_Tstat_c, t_varT_c, t_G1tilde_P_G2tilde_Vec, t_isFirth, t_isFirthConverge); //SAIGE_new.cpp
630631

631632
//t_indexForNonZero_vec.clear();
@@ -680,7 +681,9 @@ void setVCFobjInCPP(std::string t_vcfFileName,
680681
t_vcfField,
681682
false,
682683
t_SampleInModel);
683-
684+
685+
bool isEnd = ptr_gVCFobj->check_iterator_end();
686+
684687
}
685688

686689

@@ -692,6 +695,7 @@ void setSAIGEobjInCPP(arma::mat & t_XVX,
692695
arma::mat & t_XXVX_inv,
693696
arma::mat & t_XV,
694697
arma::mat & t_XVX_inv_XV,
698+
arma::mat & t_Sigma_iXXSigma_iX,
695699
arma::mat & t_X,
696700
arma::vec & t_S_a,
697701
arma::vec & t_res,
@@ -724,6 +728,7 @@ void setSAIGEobjInCPP(arma::mat & t_XVX,
724728
t_XXVX_inv,
725729
t_XV,
726730
t_XVX_inv_XV,
731+
t_Sigma_iXXSigma_iX,
727732
t_X,
728733
t_S_a,
729734
t_res,
@@ -750,7 +755,6 @@ void setSAIGEobjInCPP(arma::mat & t_XVX,
750755
t_pCutoffforFirth,
751756
t_offset);
752757
//ptr_gSAIGEobj->m_flagSparseGRM = false;
753-
754758
}
755759

756760

@@ -1885,6 +1889,7 @@ void assign_conditionMarkers_factors(
18851889
arma::vec & t_weight_cond
18861890
) // sample size
18871891
{
1892+
ptr_gSAIGEobj->set_flagSparseGRM_cur(ptr_gSAIGEobj->m_flagSparseGRM);
18881893
bool isImpute = false;
18891894
unsigned int q = t_genoIndex.size();
18901895
arma::mat P1Mat(q, t_n);
@@ -1954,7 +1959,6 @@ void assign_conditionMarkers_factors(
19541959
std::remove(end_prev);
19551960
}
19561961

1957-
19581962
bool isReadMarker = Unified_getOneMarker(t_genoType, gIndex_prev, gIndex, ref, alt, marker, pd, chr, altFreq, altCounts, missingRate, imputeInfo,
19591963
isOutputIndexForMissing, // bool t_isOutputIndexForMissing,
19601964
indexForMissing,
@@ -1981,7 +1985,7 @@ void assign_conditionMarkers_factors(
19811985
// exit(EXIT_FAILURE);
19821986
//}
19831987
}
1984-
1988+
19851989
flip = imputeGenoAndFlip(GVec, altFreq, altCounts, indexForMissing, g_impute_method, g_dosage_zerod_cutoff, g_dosage_zerod_MAC_cutoff, MAC, indexZeroVec, indexNonZeroVec);
19861990

19871991

@@ -1990,11 +1994,11 @@ void assign_conditionMarkers_factors(
19901994
indexNonZeroVec_arma = arma::conv_to<arma::uvec>::from(indexNonZeroVec);
19911995

19921996

1993-
19941997
MAF = std::min(altFreq, 1 - altFreq);
19951998
MAC = std::min(altCounts, 2*t_n-altCounts);
19961999

19972000
arma::vec gtildeVec;
2001+
19982002
Unified_getMarkerPval(
19992003
GVec,
20002004
false, // bool t_isOnlyOutputNonZero,

‎src/SAIGE_test.cpp‎

Lines changed: 16 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,7 @@ SAIGEClass::SAIGEClass(
2626
arma::mat t_XXVX_inv,
2727
arma::mat & t_XV,
2828
arma::mat & t_XVX_inv_XV,
29+
arma::mat & t_Sigma_iXXSigma_iX,
2930
arma::mat & t_X,
3031
arma::vec & t_S_a,
3132
arma::vec & t_res,
@@ -52,10 +53,18 @@ SAIGEClass::SAIGEClass(
5253
double t_pCutoffforFirth,
5354
arma::vec & t_offset){
5455

56+
5557
m_XVX = t_XVX;
5658
m_XV = t_XV;
5759
m_XXVX_inv = t_XXVX_inv;
5860
m_XVX_inv_XV = t_XVX_inv_XV;
61+
m_Sigma_iXXSigma_iX = t_Sigma_iXXSigma_iX;
62+
m_isVarPsadj = false;
63+
if(m_Sigma_iXXSigma_iX.n_cols == 1 && m_Sigma_iXXSigma_iX.n_rows == 1){
64+
m_isVarPsadj = false;
65+
}else{
66+
m_isVarPsadj = true;
67+
}
5968
m_X = t_X;
6069
m_S_a = t_S_a;
6170
m_res = t_res;
@@ -142,12 +151,16 @@ void SAIGEClass::scoreTest(arma::vec & t_GVec,
142151

143152

144153
if(!m_flagSparseGRM_cur){
145-
t_P2Vec = t_gtilde % m_mu2 *m_tauvec[0];
154+
t_P2Vec = t_gtilde % m_mu2 *m_tauvec[0];
155+
var2m = dot(t_P2Vec , t_gtilde);
146156
}else{
147157
arma::sp_mat m_SigmaMat_sp = gen_sp_SigmaMat();
148158
t_P2Vec = arma::spsolve(m_SigmaMat_sp, t_gtilde);
159+
var2m = dot(t_P2Vec , t_gtilde);
160+
if(m_isVarPsadj){
161+
var2m = var2m - t_gtilde.t() * m_Sigma_iXXSigma_iX * m_X.t() * t_P2Vec;
162+
}
149163
}
150-
var2m = dot(t_P2Vec , t_gtilde);
151164
var2 = var2m(0,0);
152165
double var1 = var2 * m_varRatioVal;
153166
double stat = S*S/var1;
@@ -349,12 +362,12 @@ void SAIGEClass::getMarkerPval(arma::vec & t_GVec,
349362
//arma::vec t_gtilde;
350363
bool isScoreFast = true;
351364

352-
353365
//if((t_altFreq >= 0.3 && t_altFreq <= 0.7) || m_flagSparseGRM || is_region){
354366
if(m_flagSparseGRM_cur){
355367
isScoreFast = false;
356368
}
357369

370+
358371
//for test
359372
//arma::vec timeoutput3 = getTime();
360373
if(!isScoreFast){

‎src/SAIGE_test.hpp‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,7 @@ class SAIGEClass
1414
arma::mat m_XVX;
1515
arma::mat m_XVX_inv_XV;
1616
arma::mat m_X;
17+
arma::mat m_Sigma_iXXSigma_iX;
1718
arma::vec m_res;
1819
arma::vec m_mu;
1920
arma::vec m_mu2;
@@ -71,7 +72,7 @@ class SAIGEClass
7172
bool m_is_Firth_beta;
7273
double m_pCutoffforFirth;
7374
arma::vec m_offset;
74-
75+
bool m_isVarPsadj;
7576
////////////////////// -------------------- functions ---------------------------------- //////////////////////
7677

7778

@@ -80,6 +81,7 @@ class SAIGEClass
8081
arma::mat t_XXVX_inv,
8182
arma::mat & t_XV,
8283
arma::mat & t_XVX_inv_XV,
84+
arma::mat & t_sigmainvX_XsignmainXtXinv,
8385
arma::mat & t_X,
8486
arma::vec & t_S_a,
8587
arma::vec & t_res,

0 commit comments

Comments
 (0)