aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/gemma.cpp19
-rw-r--r--src/lmm.cpp29
-rw-r--r--src/mvlmm.cpp27
3 files changed, 51 insertions, 24 deletions
diff --git a/src/gemma.cpp b/src/gemma.cpp
index 55eb1ee..c85670d 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -43,7 +43,7 @@
#include "lmm_float.h" //for LMM class, and functions CalcLambda, CalcPve, CalcVgVe
#include "mvlmm_float.h" //for MVLMM class
#include "prdt_float.h" //for PRDT class
-
+#include "varcov_float.h" //for MVLMM class
#include "mathfunc_float.h" //for a few functions
#else
#include "io.h"
@@ -51,11 +51,11 @@
#include "vc.h"
#include "lm.h"
#include "bslmm.h"
-
+#include "ldr.h"
#include "lmm.h"
#include "mvlmm.h"
#include "prdt.h"
-
+#include "varcov.h"
#include "mathfunc.h"
#endif
@@ -65,7 +65,7 @@ using namespace std;
GEMMA::GEMMA(void):
-version("0.95alpha"), date("07/11/2015"), year("2011")
+version("0.95alpha"), date("07/24/2016"), year("2011")
{}
void GEMMA::PrintHeader (void)
@@ -1351,7 +1351,6 @@ void GEMMA::BatchRun (PARAM &cPar)
gsl_matrix_free (G);
}
- /*
//Compute the LDSC weights (not implemented yet)
if (cPar.a_mode==72) {
cout<<"Calculating Weights ... "<<endl;
@@ -1367,7 +1366,7 @@ void GEMMA::BatchRun (PARAM &cPar)
cVarcov.CopyToParam(cPar);
}
- */
+
//Compute the S matrix (and its variance), that is used for variance component estimation using summary statistics
if (cPar.a_mode==25 || cPar.a_mode==26) {
@@ -1475,7 +1474,7 @@ void GEMMA::BatchRun (PARAM &cPar)
gsl_vector_free (s);
}
- /*
+
//Calculate SNP covariance
if (cPar.a_mode==71) {
VARCOV cVarcov;
@@ -1489,7 +1488,7 @@ void GEMMA::BatchRun (PARAM &cPar)
cVarcov.CopyToParam(cPar);
}
- */
+
//LM
if (cPar.a_mode==51 || cPar.a_mode==52 || cPar.a_mode==53 || cPar.a_mode==54) { //Fit LM
@@ -2403,7 +2402,7 @@ void GEMMA::BatchRun (PARAM &cPar)
}
- /*
+
//LDR
if (cPar.a_mode==14) {
gsl_vector *y=gsl_vector_alloc (cPar.ni_test);
@@ -2431,7 +2430,7 @@ void GEMMA::BatchRun (PARAM &cPar)
gsl_matrix_free (W);
gsl_matrix_free (G);
}
- */
+
cPar.time_total=(clock()-time_begin)/(double(CLOCKS_PER_SEC)*60.0);
return;
diff --git a/src/lmm.cpp b/src/lmm.cpp
index af6ff8a..044f33c 100644
--- a/src/lmm.cpp
+++ b/src/lmm.cpp
@@ -1243,7 +1243,11 @@ void LMM::AnalyzeBimbam (const gsl_matrix *U, const gsl_vector *eval, const gsl_
// }
//start reading genotypes and analyze
- size_t c=0;
+ size_t c=0, t_last=0;
+ for (size_t t=0; t<indicator_snp.size(); ++t) {
+ if (indicator_snp[t]==0) {continue;}
+ t_last++;
+ }
for (size_t t=0; t<indicator_snp.size(); ++t) {
// if (t>1) {break;}
!safeGetline(infile, line).eof();
@@ -1275,7 +1279,7 @@ void LMM::AnalyzeBimbam (const gsl_matrix *U, const gsl_vector *eval, const gsl_
for (size_t i=0; i<ni_test; ++i) {
if (gsl_vector_get (x_miss, i)==0) {gsl_vector_set(x, i, x_mean);}
- geno=gsl_vector_get(x, i);
+ //geno=gsl_vector_get(x, i);
//if (x_mean>1) {
// gsl_vector_set(x, i, 2-geno);
//}
@@ -1292,7 +1296,7 @@ void LMM::AnalyzeBimbam (const gsl_matrix *U, const gsl_vector *eval, const gsl_
gsl_vector_memcpy (&Xlarge_col.vector, x);
c++;
- if (c%msize==0 || t==indicator_snp.size()-1 ) {
+ if (c%msize==0 || c==t_last) {
size_t l=0;
if (c%msize==0) {l=msize;} else {l=c%msize;}
@@ -1332,12 +1336,14 @@ void LMM::AnalyzeBimbam (const gsl_matrix *U, const gsl_vector *eval, const gsl_
p_lrt=gsl_cdf_chisq_Q (2.0*(logl_H1-logl_mle_H0), 1);
}
+ //if (p_wald<0) {cout<<t<<"\t"<<i<<"\t"<<l<<endl;}
//if (x_mean>1) {beta*=-1;}
time_opt+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
//store summary data
SUMSTAT SNPs={beta, se, lambda_remle, lambda_mle, p_wald, p_lrt, p_score};
+
sumStat.push_back(SNPs);
}
}
@@ -1412,7 +1418,11 @@ void LMM::AnalyzePlink (const gsl_matrix *U, const gsl_vector *eval, const gsl_m
b=ch[0];
}
- size_t c=0;
+ size_t c=0, t_last=0;
+ for (size_t t=0; t<snpInfo.size(); ++t) {
+ if (indicator_snp[t]==0) {continue;}
+ t_last++;
+ }
for (vector<SNPINFO>::size_type t=0; t<snpInfo.size(); ++t) {
if (t%d_pace==0 || t==snpInfo.size()-1) {ProgressBar ("Reading SNPs ", t, snpInfo.size()-1);}
if (indicator_snp[t]==0) {continue;}
@@ -1463,7 +1473,7 @@ void LMM::AnalyzePlink (const gsl_matrix *U, const gsl_vector *eval, const gsl_m
gsl_vector_memcpy (&Xlarge_col.vector, x);
c++;
- if (c%msize==0 || t==indicator_snp.size()-1 ) {
+ if (c%msize==0 || c==t_last) {
size_t l=0;
if (c%msize==0) {l=msize;} else {l=c%msize;}
@@ -1612,12 +1622,17 @@ void LMM::Analyzebgen (const gsl_matrix *U, const gsl_vector *eval, const gsl_ma
//start reading genotypes and analyze
- size_t c=0;
+ size_t c=0, t_last=0;
+ for (size_t t=0; t<indicator_snp.size(); ++t) {
+ if (indicator_snp[t]==0) {continue;}
+ t_last++;
+ }
for (size_t t=0; t<indicator_snp.size(); ++t)
{
// if (t>1) {break;}
if (t%d_pace==0 || t==(ns_total-1)) {ProgressBar ("Reading SNPs ", t, ns_total-1);}
+ if (indicator_snp[t]==0) {continue;}
// read SNP header
id.clear();
rs.clear();
@@ -1737,7 +1752,7 @@ void LMM::Analyzebgen (const gsl_matrix *U, const gsl_vector *eval, const gsl_ma
gsl_vector_memcpy (&Xlarge_col.vector, x);
c++;
- if (c%msize==0 || t==indicator_snp.size()-1 ) {
+ if (c%msize==0 || c==t_last ) {
size_t l=0;
if (c%msize==0) {l=msize;} else {l=c%msize;}
diff --git a/src/mvlmm.cpp b/src/mvlmm.cpp
index 7655b50..5deda06 100644
--- a/src/mvlmm.cpp
+++ b/src/mvlmm.cpp
@@ -3185,12 +3185,17 @@ void MVLMM::Analyzebgen (const gsl_matrix *U, const gsl_vector *eval, const gsl_
//start reading genotypes and analyze
- size_t csnp=0;
+ size_t csnp=0, t_last=0;
+ for (size_t t=0; t<indicator_snp.size(); ++t) {
+ if (indicator_snp[t]==0) {continue;}
+ t_last++;
+ }
for (size_t t=0; t<indicator_snp.size(); ++t) {
// if (t>1) {break;}
if (t%d_pace==0 || t==(ns_total-1)) {ProgressBar ("Reading SNPs ", t, ns_total-1);}
+ if (indicator_snp[t]==0) {continue;}
// read SNP header
id.clear();
rs.clear();
@@ -3293,7 +3298,7 @@ void MVLMM::Analyzebgen (const gsl_matrix *U, const gsl_vector *eval, const gsl_
for (size_t i=0; i<ni_test; ++i) {
if (gsl_vector_get (x_miss, i)==0) {gsl_vector_set(x, i, x_mean);}
- geno=gsl_vector_get(x, i);
+ //geno=gsl_vector_get(x, i);
//if (x_mean>1) {
//gsl_vector_set(x, i, 2-geno);
//}
@@ -3310,7 +3315,7 @@ void MVLMM::Analyzebgen (const gsl_matrix *U, const gsl_vector *eval, const gsl_
gsl_vector_memcpy (&Xlarge_col.vector, x);
csnp++;
- if (csnp%msize==0 || t==indicator_snp.size()-1 ) {
+ if (csnp%msize==0 || c==t_last ) {
size_t l=0;
if (csnp%msize==0) {l=msize;} else {l=csnp%msize;}
@@ -3656,7 +3661,11 @@ void MVLMM::AnalyzeBimbam (const gsl_matrix *U, const gsl_vector *eval, const gs
gsl_matrix_memcpy (B_null, B);
//start reading genotypes and analyze
- size_t csnp=0;
+ size_t csnp=0, t_last=0;
+ for (size_t t=0; t<indicator_snp.size(); ++t) {
+ if (indicator_snp[t]==0) {continue;}
+ t_last++;
+ }
for (size_t t=0; t<indicator_snp.size(); ++t) {
//if (t>=1) {break;}
!safeGetline(infile, line).eof();
@@ -3705,7 +3714,7 @@ void MVLMM::AnalyzeBimbam (const gsl_matrix *U, const gsl_vector *eval, const gs
gsl_vector_memcpy (&Xlarge_col.vector, x);
csnp++;
- if (csnp%msize==0 || t==indicator_snp.size()-1 ) {
+ if (csnp%msize==0 || c==t_last ) {
size_t l=0;
if (csnp%msize==0) {l=msize;} else {l=csnp%msize;}
@@ -4068,7 +4077,11 @@ void MVLMM::AnalyzePlink (const gsl_matrix *U, const gsl_vector *eval, const gsl
b=ch[0];
}
- size_t csnp=0;
+ size_t csnp=0, t_last=0;
+ for (size_t t=0; t<indicator_snp.size(); ++t) {
+ if (indicator_snp[t]==0) {continue;}
+ t_last++;
+ }
for (vector<SNPINFO>::size_type t=0; t<snpInfo.size(); ++t) {
if (t%d_pace==0 || t==snpInfo.size()-1) {ProgressBar ("Reading SNPs ", t, snpInfo.size()-1);}
if (indicator_snp[t]==0) {continue;}
@@ -4135,7 +4148,7 @@ void MVLMM::AnalyzePlink (const gsl_matrix *U, const gsl_vector *eval, const gsl
gsl_vector_memcpy (&Xlarge_col.vector, x);
csnp++;
- if (csnp%msize==0 || t==indicator_snp.size()-1 ) {
+ if (csnp%msize==0 || c==t_last ) {
size_t l=0;
if (csnp%msize==0) {l=msize;} else {l=csnp%msize;}
generated by cgit v1.2.3 (git 2.44.0) at 2025-01-19 02:25:38 +0000