version 0.95alpha

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;}