add GXE test

1 files changed, 920 insertions, 202 deletions

diff --git a/src/param.cpp b/src/param.cpp
index 7a89ff8..c4b234a 100644
--- a/src/param.cpp
+++ b/src/param.cpp
@@ -24,6 +24,15 @@
 #include <cmath>
 #include <algorithm>
 
+#include "gsl/gsl_randist.h"
+#include "gsl/gsl_matrix.h"
+#include "gsl/gsl_vector.h"
+#include "gsl/gsl_matrix.h"
+#include "gsl/gsl_linalg.h"
+#include "gsl/gsl_blas.h"
+
+#include "eigenlib.h"
+#include "mathfunc.h"
 
 #ifdef FORCE_FLOAT
 #include "param_float.h"
@@ -39,12 +48,12 @@ using namespace std;
 
 
 
-PARAM::PARAM(void):	
+PARAM::PARAM(void):
 mode_silence (false), a_mode (0), k_mode(1), d_pace (100000),
 file_out("result"), path_out("./output/"),
 miss_level(0.05), maf_level(0.01), hwe_level(0), r2_level(0.9999),
 l_min(1e-5), l_max(1e5), n_region(10),p_nr(0.001),em_prec(0.0001),nr_prec(0.0001),em_iter(10000),nr_iter(100),crt(0),
-pheno_mean(0),
+pheno_mean(0), noconstrain (false),
 h_min(-1), h_max(-1),	h_scale(-1),
 rho_min(0.0), rho_max(1.0),	rho_scale(-1),
 logp_min(0.0), logp_max(0.0), logp_scale(-1),
@@ -55,53 +64,64 @@ n_accept(0),
 n_mh(10),
 geo_mean(2000.0),
 randseed(-1),
+window_cm(0), window_bp(0), window_ns(0),
 error(false),
-  n_cvt(1), n_vc(1),
+ni_subsample(0), n_cvt(1), n_vc(1),
 time_total(0.0), time_G(0.0), time_eigen(0.0), time_UtX(0.0), time_UtZ(0.0), time_opt(0.0), time_Omega(0.0)
 {}
 
 
 //read files
 //obtain ns_total, ng_total, ns_test, ni_test
-void PARAM::ReadFiles (void) 
+void PARAM::ReadFiles (void)
 {
 	string file_str;
-	if (!file_mk.empty()) {				
+
+
+	if (!file_cat.empty()) {
+	  if (ReadFile_cat (file_cat, mapRS2cat, n_vc)==false) {error=true;}
+	}
+
+	if (!file_var.empty()) {
+	  if (ReadFile_var (file_var, mapRS2var)==false) {error=true;}
+	}
+
+	if (!file_mk.empty()) {
 	  if (CountFileLines (file_mk, n_vc)==false) {error=true;}
 	}
-	
+
 	if (!file_snps.empty()) {
 		if (ReadFile_snps (file_snps, setSnps)==false) {error=true;}
 	} else {
 		setSnps.clear();
 	}
-	
+
 	//for prediction
 	if (!file_epm.empty()) {
 		if (ReadFile_est (file_epm, est_column, mapRS2est)==false) {error=true;}
-		
+
 		if (!file_bfile.empty()) {
 			file_str=file_bfile+".bim";
-			if (ReadFile_bim (file_str, snpInfo)==false) {error=true;}		
-			
+			if (ReadFile_bim (file_str, snpInfo)==false) {error=true;}
+
 			file_str=file_bfile+".fam";
-			if (ReadFile_fam (file_str, indicator_pheno, pheno, mapID2num, p_column)==false) {error=true;}			
+			if (ReadFile_fam (file_str, indicator_pheno, pheno, mapID2num, p_column)==false) {error=true;}
 		}
-		
-		if (!file_geno.empty()) {			
-			if (ReadFile_pheno (file_pheno, indicator_pheno, pheno, p_column)==false) {error=true;}		
-			
-			if (CountFileLines (file_geno, ns_total)==false) {error=true;}	
+
+		if (!file_geno.empty()) {
+			if (ReadFile_pheno (file_pheno, indicator_pheno, pheno, p_column)==false) {error=true;}
+
+			if (CountFileLines (file_geno, ns_total)==false) {error=true;}
 		}
-		
+
 		if (!file_ebv.empty() ) {
 			if (ReadFile_column (file_ebv, indicator_bv, vec_bv, 1)==false) {error=true;}
 		}
-		
+
 		if (!file_log.empty() ) {
 			if (ReadFile_log (file_log, pheno_mean)==false) {error=true;}
 		}
-		
+
 		//convert indicator_pheno to indicator_idv
 		int k=1;
 		for (size_t i=0; i<indicator_pheno.size(); i++) {
@@ -111,46 +131,80 @@ void PARAM::ReadFiles (void)
 			}
 			indicator_idv.push_back(k);
 		}
-		
+
 		ns_test=0;
-		
+
 		return;
 	}
-	
+
 	//read covariates before the genotype files
 	if (!file_cvt.empty() ) {
 		if (ReadFile_cvt (file_cvt, indicator_cvt, cvt, n_cvt)==false) {error=true;}
 
 		if ((indicator_cvt).size()==0) {
 			n_cvt=1;
-		} 		
+		}
 	} else {
 		n_cvt=1;
 	}
 
+	if (!file_gxe.empty() ) {
+	  if (ReadFile_column (file_gxe, indicator_gxe, gxe, 1)==false) {error=true;}
+	}
+	if (!file_weight.empty() ) {
+	  if (ReadFile_column (file_weight, indicator_weight, weight, 1)==false) {error=true;}
+	}
+
+
+	// WJA added
+	//read genotype and phenotype file for bgen format
+	if (!file_oxford.empty()) {
+		file_str=file_oxford+".sample";
+		if (ReadFile_sample(file_str, indicator_pheno, pheno, p_column,indicator_cvt, cvt, n_cvt)==false) {error=true;}
+		if ((indicator_cvt).size()==0) {
+			n_cvt=1;
+		}
+	//	n_cvt=1;
+
+		//post-process covariates and phenotypes, obtain ni_test, save all useful covariates
+		ProcessCvtPhen();
+
+
+		//obtain covariate matrix
+		gsl_matrix *W=gsl_matrix_alloc (ni_test, n_cvt);
+		CopyCvt (W);
+
+		file_str=file_oxford+".bgen";
+		if (ReadFile_bgen (file_str, setSnps, W, indicator_idv, indicator_snp, snpInfo, maf_level, miss_level, hwe_level, r2_level, ns_test)==false) {error=true;}
+		gsl_matrix_free(W);
+
+		ns_total=indicator_snp.size();
+	}
+
+
 	//read genotype and phenotype file for plink format
 	if (!file_bfile.empty()) {
 		file_str=file_bfile+".bim";
-		if (ReadFile_bim (file_str, snpInfo)==false) {error=true;}		
-		
+		if (ReadFile_bim (file_str, snpInfo)==false) {error=true;}
+
 		file_str=file_bfile+".fam";
 		if (ReadFile_fam (file_str, indicator_pheno, pheno, mapID2num, p_column)==false) {error=true;}
-		
+
 		//post-process covariates and phenotypes, obtain ni_test, save all useful covariates
 		ProcessCvtPhen();
-		
+
 		//obtain covariate matrix
 		gsl_matrix *W=gsl_matrix_alloc (ni_test, n_cvt);
 		CopyCvt (W);
-		
+
 		file_str=file_bfile+".bed";
 		if (ReadFile_bed (file_str, setSnps, W, indicator_idv, indicator_snp, snpInfo, maf_level, miss_level, hwe_level, r2_level, ns_test)==false) {error=true;}
-		
+
 		gsl_matrix_free(W);
-		
+
 		ns_total=indicator_snp.size();
 	}
-	
+
 	//read genotype and phenotype file for bimbam format
 	if (!file_geno.empty()) {
 		//annotation file before genotype file
@@ -163,7 +217,7 @@ void PARAM::ReadFiles (void)
 
 		//post-process covariates and phenotypes, obtain ni_test, save all useful covariates
 		ProcessCvtPhen();
-		
+
 		//obtain covariate matrix
 		gsl_matrix *W=gsl_matrix_alloc (ni_test, n_cvt);
 		CopyCvt (W);
@@ -171,13 +225,13 @@ void PARAM::ReadFiles (void)
 		if (ReadFile_geno (file_geno, setSnps, W, indicator_idv, indicator_snp, maf_level, miss_level, hwe_level, r2_level, mapRS2chr, mapRS2bp, mapRS2cM, snpInfo, ns_test)==false) {error=true;}
 
 		gsl_matrix_free(W);
-		
+
 		ns_total=indicator_snp.size();
 	}
-	
+
 	if (!file_gene.empty()) {
 		if (ReadFile_pheno (file_pheno, indicator_pheno, pheno, p_column)==false) {error=true;}
-		
+
 		//convert indicator_pheno to indicator_idv
 		int k=1;
 		for (size_t i=0; i<indicator_pheno.size(); i++) {
@@ -187,32 +241,39 @@ void PARAM::ReadFiles (void)
 			}
 			indicator_idv.push_back(k);
 		}
-		
-		if (ReadFile_gene (file_gene, vec_read, snpInfo, ng_total)==false) {error=true;}	
+
+        //post-process covariates and phenotypes, obtain ni_test, save all useful covariates
+        ProcessCvtPhen();
+
+        //obtain covariate matrix
+        gsl_matrix *W=gsl_matrix_alloc (ni_test, n_cvt);
+        CopyCvt (W);
+
+		if (ReadFile_gene (file_gene, vec_read, snpInfo, ng_total)==false) {error=true;}
 	}
-	
-				
+
+
 	//read is after gene file
 	if (!file_read.empty() ) {
 		if (ReadFile_column (file_read, indicator_read, vec_read, 1)==false) {error=true;}
-		
-		ni_test=0; 
+
+		ni_test=0;
 		for (vector<int>::size_type i=0; i<(indicator_idv).size(); ++i) {
 			indicator_idv[i]*=indicator_read[i];
 			ni_test+=indicator_idv[i];
 		}
-		
+
 		if (ni_test==0) {
 			error=true;
 			cout<<"error! number of analyzed individuals equals 0. "<<endl;
 			return;
 		}
 	}
-	
+
 	//for ridge prediction, read phenotype only
 	if (file_geno.empty() && file_gene.empty() && !file_pheno.empty()) {
-		if (ReadFile_pheno (file_pheno, indicator_pheno, pheno, p_column)==false) {error=true;}	
-				
+		if (ReadFile_pheno (file_pheno, indicator_pheno, pheno, p_column)==false) {error=true;}
+
 		//post-process covariates and phenotypes, obtain ni_test, save all useful covariates
 		ProcessCvtPhen();
 	}
@@ -225,37 +286,43 @@ void PARAM::ReadFiles (void)
 
 
 
-void PARAM::CheckParam (void) 
-{	
+void PARAM::CheckParam (void)
+{
 	struct stat fileInfo;
 	string str;
-	
+
 	//check parameters
 	if (k_mode!=1 && k_mode!=2) {cout<<"error! unknown kinship/relatedness input mode: "<<k_mode<<endl; error=true;}
-	if (a_mode!=1 && a_mode!=2 && a_mode!=3 && a_mode!=4 && a_mode!=5 && a_mode!=11 && a_mode!=12 && a_mode!=13 && a_mode!=21 && a_mode!=22 && a_mode!=31 && a_mode!=41 && a_mode!=42 && a_mode!=43 && a_mode!=51 && a_mode!=52 && a_mode!=53 && a_mode!=54 && a_mode!=61)   
-	{cout<<"error! unknown analysis mode: "<<a_mode<<". make sure -gk or -eigen or -lmm or -bslmm or -predict is sepcified correctly."<<endl; error=true;}
+	if (a_mode!=1 && a_mode!=2 && a_mode!=3 && a_mode!=4 && a_mode!=5 && a_mode!=11 && a_mode!=12 && a_mode!=13 && a_mode!=14 && a_mode!=21 && a_mode!=22 && a_mode!=25 && a_mode!=26 && a_mode!=27 && a_mode!=28 && a_mode!=31 && a_mode!=41 && a_mode!=42 && a_mode!=43 && a_mode!=51 && a_mode!=52 && a_mode!=53 && a_mode!=54 && a_mode!=61 && a_mode!=62 && a_mode!=71)
+	{cout<<"error! unknown analysis mode: "<<a_mode<<". make sure -gk or -eigen or -lmm or -bslmm -predict or -calccov is sepcified correctly."<<endl; error=true;}
 	if (miss_level>1) {cout<<"error! missing level needs to be between 0 and 1. current value = "<<miss_level<<endl; error=true;}
 	if (maf_level>0.5) {cout<<"error! maf level needs to be between 0 and 0.5. current value = "<<maf_level<<endl; error=true;}
 	if (hwe_level>1) {cout<<"error! hwe level needs to be between 0 and 1. current value = "<<hwe_level<<endl; error=true;}
 	if (r2_level>1) {cout<<"error! r2 level needs to be between 0 and 1. current value = "<<r2_level<<endl; error=true;}
-	
-	if (l_max<l_min) {cout<<"error! maximum lambda value must be larger than the minimal value. current values = "<<l_max<<" and "<<l_min<<endl; error=true;}	
+
+	if (l_max<l_min) {cout<<"error! maximum lambda value must be larger than the minimal value. current values = "<<l_max<<" and "<<l_min<<endl; error=true;}
 	if (h_max<h_min) {cout<<"error! maximum h value must be larger than the minimal value. current values = "<<h_max<<" and "<<h_min<<endl; error=true;}
 	if (s_max<s_min) {cout<<"error! maximum s value must be larger than the minimal value. current values = "<<s_max<<" and "<<s_min<<endl; error=true;}
 	if (rho_max<rho_min) {cout<<"error! maximum rho value must be larger than the minimal value. current values = "<<rho_max<<" and "<<rho_min<<endl; error=true;}
 	if (logp_max<logp_min) {cout<<"error! maximum logp value must be larger than the minimal value. current values = "<<logp_max/log(10)<<" and "<<logp_min/log(10)<<endl; error=true;}
-	
+
 	if (h_max>1) {cout<<"error! h values must be bewtween 0 and 1. current values = "<<h_max<<" and "<<h_min<<endl; error=true;}
 	if (rho_max>1) {cout<<"error! rho values must be between 0 and 1. current values = "<<rho_max<<" and "<<rho_min<<endl; error=true;}
 	if (logp_max>0) {cout<<"error! maximum logp value must be smaller than 0. current values = "<<logp_max/log(10)<<" and "<<logp_min/log(10)<<endl; error=true;}
 	if (l_max<l_min) {cout<<"error! maximum lambda value must be larger than the minimal value. current values = "<<l_max<<" and "<<l_min<<endl; error=true;}
-		
+
 	if (h_scale>1.0) {cout<<"error! hscale value must be between 0 and 1. current value = "<<h_scale<<endl; error=true;}
 	if (rho_scale>1.0) {cout<<"error! rscale value must be between 0 and 1. current value = "<<rho_scale<<endl; error=true;}
 	if (logp_scale>1.0) {cout<<"error! pscale value must be between 0 and 1. current value = "<<logp_scale<<endl; error=true;}
 
 	if (rho_max==1 && rho_min==1 && a_mode==12) {cout<<"error! ridge regression does not support a rho parameter. current values = "<<rho_max<<" and "<<rho_min<<endl; error=true;}
-		
+
+	if (window_cm<0) {cout<<"error! windowcm values must be non-negative. current values = "<<window_cm<<endl; error=true;}
+
+	if (window_cm==0 && window_bp==0 && window_ns==0) {
+	  window_bp=1000000;
+	}
+
 	//check p_column, and (no need to) sort p_column into ascending order
 	if (p_column.size()==0) {
 		p_column.push_back(1);
@@ -266,12 +333,12 @@ void PARAM::CheckParam (void)
 			}
 		}
 	}
-	
+
 	//sort (p_column.begin(), p_column.end() );
 	n_ph=p_column.size();
-	
-		
-	
+
+
+
 	//only lmm option (and one prediction option) can deal with multiple phenotypes
 	//and no gene expression files
 	if (n_ph>1 && a_mode!=1 && a_mode!=2 && a_mode!=3 && a_mode!=4 && a_mode!=43) {
@@ -280,11 +347,11 @@ void PARAM::CheckParam (void)
 	if (n_ph>1 && !file_gene.empty() ) {
 		cout<<"error! multiple phenotype analysis option not allowed with gene expression files. "<<endl; error=true;
 	}
-	
+
 	if (p_nr>1) {
 		cout<<"error! pnr value must be between 0 and 1. current value = "<<p_nr<<endl; error=true;
 	}
-	
+
 	//check est_column
 	if (est_column.size()==0) {
 		if (file_ebv.empty()) {
@@ -299,10 +366,10 @@ void PARAM::CheckParam (void)
 			est_column.push_back(7);
 		}
 	}
-	
-	if (est_column.size()!=4) {cout<<"error! -en not followed by four numbers. current number = "<<est_column.size()<<endl; error=true;}	
+
+	if (est_column.size()!=4) {cout<<"error! -en not followed by four numbers. current number = "<<est_column.size()<<endl; error=true;}
 	if (est_column[0]==0) {cout<<"error! -en rs column can not be zero. current number = "<<est_column.size()<<endl; error=true;}
-	
+
 	//check if files are compatible with each other, and if files exist
 	if (!file_bfile.empty()) {
 		str=file_bfile+".bim";
@@ -310,44 +377,101 @@ void PARAM::CheckParam (void)
 		str=file_bfile+".bed";
 		if (stat(str.c_str(),&fileInfo)==-1) {cout<<"error! fail to open .bed file: "<<str<<endl; error=true;}
 		str=file_bfile+".fam";
-		if (stat(str.c_str(),&fileInfo)==-1) {cout<<"error! fail to open .fam file: "<<str<<endl; error=true;}			
+		if (stat(str.c_str(),&fileInfo)==-1) {cout<<"error! fail to open .fam file: "<<str<<endl; error=true;}
+	}
+
+	if (!file_oxford.empty()) {
+		str=file_bfile+".bgen";
+		if (stat(str.c_str(),&fileInfo)==-1) {cout<<"error! fail to open .bgen file: "<<str<<endl; error=true;}
+		str=file_bfile+".sample";
+		if (stat(str.c_str(),&fileInfo)==-1) {cout<<"error! fail to open .sample file: "<<str<<endl; error=true;}
 	}
-	
+
 	if ((!file_geno.empty() || !file_gene.empty()) ) {
 		str=file_pheno;
 		if (stat(str.c_str(),&fileInfo)==-1) {cout<<"error! fail to open phenotype file: "<<str<<endl; error=true;}
-	}	
-	
+	}
+
 	str=file_geno;
 	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open mean genotype file: "<<str<<endl; error=true;}
-	
+
 	str=file_gene;
 	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open gene expression file: "<<str<<endl; error=true;}
-	
+
+	str=file_cat;
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open category file: "<<str<<endl; error=true;}
+
+	str=file_var;
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open category file: "<<str<<endl; error=true;}
+
+	str=file_beta;
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open beta file: "<<str<<endl; error=true;}
+
+	str=file_cor;
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open correlation file: "<<str<<endl; error=true;}
+
+	str=file_q;
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open q file: "<<str<<endl; error=true;}
+
+	str=file_s;
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open s file: "<<str<<endl; error=true;}
+
+	str=file_v;
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open v file: "<<str<<endl; error=true;}
+
+	str=file_mq;
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open mq file: "<<str<<endl; error=true;}
+
+	str=file_ms;
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open ms file: "<<str<<endl; error=true;}
+
+	str=file_mv;
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open mv file: "<<str<<endl; error=true;}
+
 	size_t flag=0;
 	if (!file_bfile.empty()) {flag++;}
 	if (!file_geno.empty()) {flag++;}
 	if (!file_gene.empty()) {flag++;}
-	
-	if (flag!=1 && a_mode!=43 && a_mode!=5 && a_mode!=61) {
+	// WJA added
+	if (!file_oxford.empty()) {flag++;}
+
+	if (flag!=1 && a_mode!=27 && a_mode!=28 && a_mode!=43 && a_mode!=5 && a_mode!=61 && a_mode!=62) {
 		cout<<"error! either plink binary files, or bimbam mean genotype files, or gene expression files are required."<<endl; error=true;
 	}
-	
-	if (file_pheno.empty() && (a_mode==43 || a_mode==5 || a_mode==61) ) {
+
+	if (file_pheno.empty() && (a_mode==43 || a_mode==5) ) {
 		cout<<"error! phenotype file is required."<<endl; error=true;
 	}
-	
+
+	if (a_mode==61 || a_mode==62) {
+	  if (!file_pheno.empty()) {
+	    if (file_kin.empty() && (file_ku.empty()||file_kd.empty()) && file_mk.empty() ) {
+	      cout<<"error! missing relatedness file. "<<endl;  error=true;
+	    }
+	  } else if (!file_cor.empty()) {
+	    if (file_beta.empty() ) {
+	      cout<<"error! missing cor file."<<endl; error=true;
+	    }
+	  } else {
+	    if ( (file_mq.empty() || file_ms.empty() || file_mv.empty() ) && (file_q.empty() || file_s.empty() || file_v.empty() )  ) {
+	      cout<<"error! either phenotype/kinship files or ms/mq/mv s/q/v files are required."<<endl; error=true;
+	    }
+	  }
+	}
+
+
+
 	if (!file_epm.empty() && file_bfile.empty() && file_geno.empty() ) {cout<<"error! estimated parameter file also requires genotype file."<<endl; error=true;}
 	if (!file_ebv.empty() && file_kin.empty()) {cout<<"error! estimated breeding value file also requires relatedness file."<<endl; error=true;}
-	
+
 	if (!file_log.empty() && pheno_mean!=0) {cout<<"error! either log file or mu value can be provide."<<endl; error=true;}
-	
+
 	str=file_snps;
 	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open snps file: "<<str<<endl; error=true;}
-	
+
 	str=file_log;
 	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open log file: "<<str<<endl; error=true;}
-	
+
 	str=file_anno;
 	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open annotation file: "<<str<<endl; error=true;}
 
@@ -356,52 +480,75 @@ void PARAM::CheckParam (void)
 
 	str=file_mk;
 	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open relatedness matrix file: "<<str<<endl; error=true;}
-	
+
 	str=file_cvt;
 	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open covariates file: "<<str<<endl; error=true;}
-	
+
+	str=file_gxe;
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open environmental covariate file: "<<str<<endl; error=true;}
+
+	str=file_weight;
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open the residual weight file: "<<str<<endl; error=true;}
+
 	str=file_epm;
 	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open estimated parameter file: "<<str<<endl; error=true;}
-	
+
 	str=file_ebv;
 	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open estimated breeding value file: "<<str<<endl; error=true;}
-	
+
 	str=file_read;
 	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open total read file: "<<str<<endl; error=true;}
-		
+
 	//check if files are compatible with analysis mode
 	if (k_mode==2 && !file_geno.empty() ) {cout<<"error! use \"-km 1\" when using bimbam mean genotype file. "<<endl; error=true;}
-	
-	if ((a_mode==1 || a_mode==2 || a_mode==3 || a_mode==4 || a_mode==5 || a_mode==31) && (file_kin.empty() && (file_ku.empty()||file_kd.empty())) )  {cout<<"error! missing relatedness file. "<<endl;  error=true;}
 
-	if (a_mode==61 && (file_kin.empty() && (file_ku.empty()||file_kd.empty()) && file_mk.empty() ) )  {cout<<"error! missing relatedness file. "<<endl;  error=true;}
+	if ((a_mode==1 || a_mode==2 || a_mode==3 || a_mode==4 || a_mode==5 || a_mode==31) && (file_kin.empty() && (file_ku.empty()||file_kd.empty())) )  {cout<<"error! missing relatedness file. "<<endl;  error=true;}
 
 	if ((a_mode==43) && file_kin.empty())  {cout<<"error! missing relatedness file. -predict option requires -k option to provide a relatedness file."<<endl;  error=true;}
-	
+
 	if ((a_mode==11 || a_mode==12 || a_mode==13) && !file_cvt.empty() ) {cout<<"error! -bslmm option does not support covariates files."<<endl; error=true;}
-		
+
 	if (a_mode==41 || a_mode==42) {
-		if (!file_cvt.empty() ) {cout<<"error! -predict option does not support covariates files."<<endl; error=true;}	
-		if (file_epm.empty() ) {cout<<"error! -predict option requires estimated parameter files."<<endl; error=true;}		
+		if (!file_cvt.empty() ) {cout<<"error! -predict option does not support covariates files."<<endl; error=true;}
+		if (file_epm.empty() ) {cout<<"error! -predict option requires estimated parameter files."<<endl; error=true;}
+	}
+
+	if (file_beta.empty() && (a_mode==27 || a_mode==28) ) {
+		cout<<"error! beta effects file is required."<<endl; error=true;
 	}
 
 	return;
 }
 
 
-		
+
 
 
 void PARAM::CheckData (void) {
-	if ((file_cvt).empty() || (indicator_cvt).size()==0) {
-		n_cvt=1;
+  if(file_oxford.empty())	// WJA NOTE: I added this condition so that covariates can be added through sample, probably not exactly what is wanted
+
+	{
+ 	if ((file_cvt).empty() || (indicator_cvt).size()==0) {
+ 		n_cvt=1;
+ 	}
 	}
+
 	if ( (indicator_cvt).size()!=0 && (indicator_cvt).size()!=(indicator_idv).size()) {
 		error=true;
 		cout<<"error! number of rows in the covariates file do not match the number of individuals. "<<endl;
 		return;
 	}
-	
+	if ( (indicator_gxe).size()!=0 && (indicator_gxe).size()!=(indicator_idv).size()) {
+		error=true;
+		cout<<"error! number of rows in the gxe file do not match the number of individuals. "<<endl;
+		return;
+	}
+	if ( (indicator_weight).size()!=0 && (indicator_weight).size()!=(indicator_idv).size()) {
+		error=true;
+		cout<<"error! number of rows in the weight file do not match the number of individuals. "<<endl;
+		return;
+	}
+
 	if ( (indicator_read).size()!=0 && (indicator_read).size()!=(indicator_idv).size()) {
 		error=true;
 		cout<<"error! number of rows in the total read file do not match the number of individuals. "<<endl;
@@ -411,13 +558,13 @@ void PARAM::CheckData (void) {
 	//calculate ni_total and ni_test, and set indicator_idv to 0 whenever indicator_cvt=0
 	//and calculate np_obs and np_miss
 	ni_total=(indicator_idv).size();
-	
-	ni_test=0; 
+
+	ni_test=0;
 	for (vector<int>::size_type i=0; i<(indicator_idv).size(); ++i) {
 		if (indicator_idv[i]==0) {continue;}
 		ni_test++;
 	}
-	
+
 	ni_cvt=0;
 	for (size_t i=0; i<indicator_cvt.size(); i++) {
 		if (indicator_cvt[i]==0) {continue;}
@@ -429,8 +576,16 @@ void PARAM::CheckData (void) {
 		if (indicator_cvt.size()!=0) {
 			if (indicator_cvt[i]==0) {continue;}
 		}
-		
-		for (size_t j=0; j<indicator_pheno[i].size(); j++) {					
+
+		if (indicator_gxe.size()!=0) {
+			if (indicator_gxe[i]==0) {continue;}
+		}
+
+		if (indicator_weight.size()!=0) {
+			if (indicator_weight[i]==0) {continue;}
+		}
+
+		for (size_t j=0; j<indicator_pheno[i].size(); j++) {
 			if (indicator_pheno[i][j]==0) {
 				np_miss++;
 			} else {
@@ -441,101 +596,103 @@ void PARAM::CheckData (void) {
 
 	/*
 	if ((indicator_cvt).size()!=0) {
-		ni_test=0; 
+		ni_test=0;
 		for (vector<int>::size_type i=0; i<(indicator_idv).size(); ++i) {
 			indicator_idv[i]*=indicator_cvt[i];
 			ni_test+=indicator_idv[i];
 		}
-	}	
-	
+	}
+
 	if ((indicator_read).size()!=0) {
-		ni_test=0; 
+		ni_test=0;
 		for (vector<int>::size_type i=0; i<(indicator_idv).size(); ++i) {
 			indicator_idv[i]*=indicator_read[i];
 			ni_test+=indicator_idv[i];
 		}
 	}
 	*/
-	if (ni_test==0) {
+	if (ni_test==0 && file_cor.empty() && file_mq.empty() && file_q.empty() && file_beta.empty() ) {
 		error=true;
 		cout<<"error! number of analyzed individuals equals 0. "<<endl;
 		return;
 	}
-	
+
 	if (a_mode==43) {
 		if (ni_cvt==ni_test) {
 			error=true;
-			cout<<"error! no individual has missing phenotypes."<<endl; 
+			cout<<"error! no individual has missing phenotypes."<<endl;
 			return;
 		}
 		if ((np_obs+np_miss)!=(ni_cvt*n_ph)) {
 			error=true;
 			//cout<<ni_cvt<<"\t"<<ni_test<<"\t"<<ni_total<<"\t"<<np_obs<<"\t"<<np_miss<<"\t"<<indicator_cvt.size()<<endl;
-			cout<<"error! number of phenotypes do not match the summation of missing and observed phenotypes."<<endl; 
+			cout<<"error! number of phenotypes do not match the summation of missing and observed phenotypes."<<endl;
 			return;
 		}
 	}
 
 	//output some information
-	cout<<"## number of total individuals = "<<ni_total<<endl;
-	if (a_mode==43) {
-		cout<<"## number of analyzed individuals = "<<ni_cvt<<endl;
-		cout<<"## number of individuals with full phenotypes = "<<ni_test<<endl;
-	} else {
-		cout<<"## number of analyzed individuals = "<<ni_test<<endl;
+	if (file_cor.empty() && file_mq.empty() && file_q.empty() ) {
+	  cout<<"## number of total individuals = "<<ni_total<<endl;
+	  if (a_mode==43) {
+	    cout<<"## number of analyzed individuals = "<<ni_cvt<<endl;
+	    cout<<"## number of individuals with full phenotypes = "<<ni_test<<endl;
+	  } else {
+	    cout<<"## number of analyzed individuals = "<<ni_test<<endl;
+	  }
+	  cout<<"## number of covariates = "<<n_cvt<<endl;
+	  cout<<"## number of phenotypes = "<<n_ph<<endl;
+	  if (a_mode==43) {
+	    cout<<"## number of observed data = "<<np_obs<<endl;
+	    cout<<"## number of missing data = "<<np_miss<<endl;
+	  }
+	  if (!file_gene.empty()) {
+	    cout<<"## number of total genes = "<<ng_total<<endl;
+	  } else if (file_epm.empty() && a_mode!=43 && a_mode!=5) {
+	    cout<<"## number of total SNPs = "<<ns_total<<endl;
+	    cout<<"## number of analyzed SNPs = "<<ns_test<<endl;
+	  } else {}
 	}
-	cout<<"## number of covariates = "<<n_cvt<<endl;
-	cout<<"## number of phenotypes = "<<n_ph<<endl;
-	if (a_mode==43) {
-		cout<<"## number of observed data = "<<np_obs<<endl;
-		cout<<"## number of missing data = "<<np_miss<<endl;
-	}
-	if (!file_gene.empty()) {		
-		cout<<"## number of total genes = "<<ng_total<<endl;
-	} else if (file_epm.empty() && a_mode!=43 && a_mode!=5) {
-		cout<<"## number of total SNPs = "<<ns_total<<endl;	
-		cout<<"## number of analyzed SNPs = "<<ns_test<<endl;
-	} else {}
-	
+
 	//set d_pace to 1000 for gene expression
 	if (!file_gene.empty() && d_pace==100000) {
 		d_pace=1000;
 	}
-	
+
 	//for case-control studies, count #cases and #controls
 	int flag_cc=0;
-	if (a_mode==13) {	
+	if (a_mode==13) {
 		ni_case=0;
 		ni_control=0;
 		for (size_t i=0; i<indicator_idv.size(); i++) {
 			if (indicator_idv[i]==0) {continue;}
-		
+
 			if (pheno[i][0]==0) {ni_control++;}
 			else if (pheno[i][0]==1) {ni_case++;}
 			else {flag_cc=1;}
 		}
-		cout<<"## number of cases = "<<ni_case<<endl;	
-		cout<<"## number of controls = "<<ni_control<<endl;	
-	}	
-	
+		cout<<"## number of cases = "<<ni_case<<endl;
+		cout<<"## number of controls = "<<ni_control<<endl;
+	}
+
 	if (flag_cc==1) {cout<<"Unexpected non-binary phenotypes for case/control analysis. Use default (BSLMM) analysis instead."<<endl; a_mode=11;}
-	
+
 	//set parameters for BSLMM
 	//and check for predict
 	if (a_mode==11 || a_mode==12 || a_mode==13) {
-		if (a_mode==11) {n_mh=1;}	
+		if (a_mode==11) {n_mh=1;}
 		if (logp_min==0) {logp_min=-1.0*log((double)ns_test);}
-	
+
 		if (h_scale==-1) {h_scale=min(1.0, 10.0/sqrt((double)ni_test) );}
 		if (rho_scale==-1) {rho_scale=min(1.0, 10.0/sqrt((double)ni_test) );}
 		if (logp_scale==-1) {logp_scale=min(1.0, 5.0/sqrt((double)ni_test) );}
-		
+
 		if (h_min==-1) {h_min=0.0;}
 		if (h_max==-1) {h_max=1.0;}
-		
+
 		if (s_max>ns_test) {s_max=ns_test; cout<<"s_max is re-set to the number of analyzed SNPs."<<endl;}
 		if (s_max<s_min) {cout<<"error! maximum s value must be larger than the minimal value. current values = "<<s_max<<" and "<<s_min<<endl; error=true;}
-	} else if (a_mode==41 || a_mode==42) {		
+	} else if (a_mode==41 || a_mode==42) {
 		if (indicator_bv.size()!=0) {
 			if (indicator_idv.size()!=indicator_bv.size()) {
 				cout<<"error! number of rows in the phenotype file does not match that in the estimated breeding value file: "<<indicator_idv.size()<<"\t"<<indicator_bv.size()<<endl;
@@ -555,18 +712,18 @@ void PARAM::CheckData (void) {
 
 	//file_mk needs to contain more than one line
 	if (n_vc==1 && !file_mk.empty()) {cout<<"error! -mk file should contain more than one line."<<endl; error=true;}
-	
+
 	return;
 }
 
 
-void PARAM::PrintSummary () 
+void PARAM::PrintSummary ()
 {
 	if (n_ph==1) {
 		cout<<"pve estimate ="<<pve_null<<endl;
 		cout<<"se(pve) ="<<pve_se_null<<endl;
 	} else {
-		
+
 	}
 	return;
 }
@@ -575,7 +732,7 @@ void PARAM::PrintSummary ()
 
 void PARAM::ReadGenotypes (gsl_matrix *UtX, gsl_matrix *K, const bool calc_K) {
 	string file_str;
-	
+
 	if (!file_bfile.empty()) {
 		file_str=file_bfile+".bed";
 		if (ReadFile_bed (file_str, indicator_idv, indicator_snp, UtX, K, calc_K)==false) {error=true;}
@@ -583,91 +740,563 @@ void PARAM::ReadGenotypes (gsl_matrix *UtX, gsl_matrix *K, const bool calc_K) {
 	else {
 		if (ReadFile_geno (file_geno, indicator_idv, indicator_snp, UtX, K, calc_K)==false) {error=true;}
 	}
-	
+
 	return;
 }
-		
+
+
+void PARAM::ReadGenotypes (vector<vector<unsigned char> > &Xt, gsl_matrix *K, const bool calc_K) {
+	string file_str;
+
+	if (!file_bfile.empty()) {
+		file_str=file_bfile+".bed";
+		if (ReadFile_bed (file_str, indicator_idv, indicator_snp, Xt, K, calc_K, ni_test, ns_test)==false) {error=true;}
+	} else {
+		if (ReadFile_geno (file_geno, indicator_idv, indicator_snp, Xt, K, calc_K, ni_test, ns_test)==false) {error=true;}
+	}
+
+	return;
+}
+
 
 
 
 void PARAM::CalcKin (gsl_matrix *matrix_kin)  {
 	string file_str;
-	
+
 	gsl_matrix_set_zero (matrix_kin);
-	
-	if (!file_bfile.empty() ) {		
+
+	if (!file_bfile.empty() ) {
 		file_str=file_bfile+".bed";
 		if (PlinkKin (file_str, indicator_snp, a_mode-20, d_pace, matrix_kin)==false) {error=true;}
 	}
+	else if (!file_oxford.empty() ) {
+		file_str=file_oxford+".bgen";
+		if (bgenKin (file_str, indicator_snp, a_mode-20, d_pace, matrix_kin)==false) {error=true;}
+ 	}
 	else {
 		file_str=file_geno;
 		if (BimbamKin (file_str, indicator_snp, a_mode-20, d_pace, matrix_kin)==false) {error=true;}
 	}
-	
+
+	return;
+}
+
+
+
+//from an existing n by nd G matrix, compute the d by d S matrix
+void compKtoS (const gsl_matrix *G, gsl_matrix *S) {
+  size_t n_vc=S->size1, ni_test=G->size1;
+  double di, dj, tr_KiKj, sum_Ki, sum_Kj, s_Ki, s_Kj, s_KiKj, si, sj, d;
+
+  for (size_t i=0; i<n_vc; i++) {
+    for (size_t j=i; j<n_vc; j++) {
+      tr_KiKj=0; sum_Ki=0; sum_Kj=0; s_KiKj=0; si=0; sj=0;
+      for (size_t l=0; l<ni_test; l++) {
+	s_Ki=0; s_Kj=0;
+	for (size_t k=0; k<ni_test; k++) {
+	  di=gsl_matrix_get(G, l, k+ni_test*i);
+	  dj=gsl_matrix_get(G, l, k+ni_test*j);
+	  s_Ki+=di; s_Kj+=dj;
+
+	  tr_KiKj+=di*dj; sum_Ki+=di; sum_Kj+=dj;
+	  if (l==k) {si+=di; sj+=dj;}
+	}
+	s_KiKj+=s_Ki*s_Kj;
+      }
+
+      sum_Ki/=(double)ni_test;
+      sum_Kj/=(double)ni_test;
+      s_KiKj/=(double)ni_test;
+      si-=sum_Ki;
+      sj-=sum_Kj;
+      d=tr_KiKj-2*s_KiKj+sum_Ki*sum_Kj;
+      d=d/(si*sj)-1/(double)(ni_test-1);
+
+      gsl_matrix_set (S, i, j, d);
+      if (i!=j) {gsl_matrix_set (S, j, i, d);}
+    }
+  }
+  //cout<<tr_KiKj<<" "<<s_KiKj<<" "<<sum_Ki<<" "<<sum_Kj<<" "<<si<<" "<<sj<<" "<<d*1000000<<endl;
+  return;
+}
+
+
+
+//copied from lmm.cpp; is used in the following function compKtoQ
+//map a number 1-(n_cvt+2) to an index between 0 and [(n_c+2)^2+(n_c+2)]/2-1
+size_t GetabIndex (const size_t a, const size_t b, const size_t n_cvt) {
+	if (a>n_cvt+2 || b>n_cvt+2 || a<=0 || b<=0) {cout<<"error in GetabIndex."<<endl; return 0;}
+	size_t index;
+	size_t l, h;
+	if (b>a) {l=a; h=b;} else {l=b; h=a;}
+
+	size_t n=n_cvt+2;
+	index=(2*n-l+2)*(l-1)/2+h-l;
+
+	return index;
+}
+
+//from an existing n by nd (centered) G matrix, compute the d+1 by d*(d+1) Q matrix
+//where inside i'th d+1 by d+1 matrix, each element is tr(KiKjKiKl)-r*tr(KjKiKl)-r*tr(KlKiKj)+r^2*tr(KjKl), where r=n/(n-1)
+void compKtoQ (const gsl_matrix *G, gsl_matrix *Q) {
+  size_t n_vc=G->size2/G->size1, ni_test=G->size1;
+
+  gsl_matrix *KiKj=gsl_matrix_alloc(ni_test, n_vc*(n_vc+1)/2*ni_test);
+  gsl_vector *trKiKjKi=gsl_vector_alloc ( n_vc*n_vc );
+  gsl_vector *trKiKj=gsl_vector_alloc( n_vc*(n_vc+1)/2 );
+  gsl_vector *trKi=gsl_vector_alloc(n_vc);
+
+  double d, tr, r=(double)ni_test/(double)(ni_test-1);
+  size_t t, t_ij, t_il, t_jl, t_ii;
+
+  //compute KiKj for all pairs of i and j (including the identity matrix)
+  t=0;
+  for (size_t i=0; i<n_vc; i++) {
+    gsl_matrix_const_view Ki=gsl_matrix_const_submatrix(G, 0, i*ni_test, ni_test, ni_test);
+    for (size_t j=i; j<n_vc; j++) {
+      gsl_matrix_const_view Kj=gsl_matrix_const_submatrix(G, 0, j*ni_test, ni_test, ni_test);
+      gsl_matrix_view KiKj_sub=gsl_matrix_submatrix (KiKj, 0, t*ni_test, ni_test, ni_test);
+      eigenlib_dgemm ("N", "N", 1.0, &Ki.matrix, &Kj.matrix, 0.0, &KiKj_sub.matrix);
+      t++;
+    }
+  }
+  /*
+  for (size_t i=0; i<5; i++) {
+    for (size_t j=0; j<5; j++) {
+      cout<<gsl_matrix_get (G, i, j)<<" ";
+    }
+    cout<<endl;
+  }
+  */
+
+  //compute trKi, trKiKj
+  t=0;
+  for (size_t i=0; i<n_vc; i++) {
+    for (size_t j=i; j<n_vc; j++) {
+      tr=0;
+      for (size_t k=0; k<ni_test; k++) {
+	tr+=gsl_matrix_get (KiKj, k, t*ni_test+k);
+      }
+      gsl_vector_set (trKiKj, t, tr);
+
+      t++;
+    }
+
+    tr=0;
+    for (size_t k=0; k<ni_test; k++) {
+      tr+=gsl_matrix_get (G, k, i*ni_test+k);
+    }
+    gsl_vector_set (trKi, i, tr);
+  }
+
+  //compute trKiKjKi (it is not symmetric w.r.t. i and j)
+  for (size_t i=0; i<n_vc; i++) {
+    for (size_t j=0; j<n_vc; j++) {
+      tr=0;
+      t=GetabIndex (i+1, j+1, n_vc-2);
+      for (size_t k=0; k<ni_test; k++) {
+	gsl_vector_const_view KiKj_row=gsl_matrix_const_subrow (KiKj, k, t*ni_test, ni_test);
+	gsl_vector_const_view KiKj_col=gsl_matrix_const_column (KiKj, t*ni_test+k);
+
+	gsl_vector_const_view Ki_col=gsl_matrix_const_column (G, i*ni_test+k);
+
+	if (i<=j) {
+	  gsl_blas_ddot (&KiKj_row.vector, &Ki_col.vector, &d);
+	  tr+=d;
+	} else {
+	  gsl_blas_ddot (&KiKj_col.vector, &Ki_col.vector, &d);
+	  tr+=d;
+	}
+      }
+      gsl_vector_set (trKiKjKi, i*n_vc+j, tr);
+    }
+  }
+
+  //compute Q
+  for (size_t i=0; i<n_vc; i++) {
+    for (size_t j=0; j<n_vc+1; j++) {
+      for (size_t l=j; l<n_vc+1; l++) {
+	if (j!=n_vc && l!=n_vc) {
+	  t_ij=GetabIndex (i+1, j+1, n_vc-2);
+	  t_il=GetabIndex (i+1, l+1, n_vc-2);
+	  t_jl=GetabIndex (j+1, l+1, n_vc-2);
+
+	  //cout<<ni_test<<" "<<r<<t_ij<<" "<<t_il<<" "<<t_jl<<" "<<endl;
+	  tr=0;
+	  for (size_t k=0; k<ni_test; k++) {
+	    gsl_vector_const_view KiKj_row=gsl_matrix_const_subrow (KiKj, k, t_ij*ni_test, ni_test);
+	    gsl_vector_const_view KiKj_col=gsl_matrix_const_column (KiKj, t_ij*ni_test+k);
+	    gsl_vector_const_view KiKl_row=gsl_matrix_const_subrow (KiKj, k, t_il*ni_test, ni_test);
+	    gsl_vector_const_view KiKl_col=gsl_matrix_const_column (KiKj, t_il*ni_test+k);
+
+	    gsl_vector_const_view Kj_row=gsl_matrix_const_subrow (G, k, j*ni_test, ni_test);
+	    gsl_vector_const_view Kl_row=gsl_matrix_const_subrow (G, k, l*ni_test, ni_test);
+
+	    if (i<=j && i<=l) {
+	      gsl_blas_ddot (&KiKj_row.vector, &KiKl_col.vector, &d);
+	      tr+=d;
+	      gsl_blas_ddot (&Kj_row.vector, &KiKl_col.vector, &d);
+	      tr-=r*d;
+	      gsl_blas_ddot (&Kl_row.vector, &KiKj_col.vector, &d);
+	      tr-=r*d;
+	    } else if (i<=j && i>l) {
+	      gsl_blas_ddot (&KiKj_row.vector, &KiKl_row.vector, &d);
+	      tr+=d;
+	      gsl_blas_ddot (&Kj_row.vector, &KiKl_row.vector, &d);
+	      tr-=r*d;
+	      gsl_blas_ddot (&Kl_row.vector, &KiKj_col.vector, &d);
+	      tr-=r*d;
+	    } else if (i>j && i<=l) {
+	      gsl_blas_ddot (&KiKj_col.vector, &KiKl_col.vector, &d);
+	      tr+=d;
+	      gsl_blas_ddot (&Kj_row.vector, &KiKl_col.vector, &d);
+	      tr-=r*d;
+	      gsl_blas_ddot (&Kl_row.vector, &KiKj_row.vector, &d);
+	      tr-=r*d;
+	    } else {
+	      gsl_blas_ddot (&KiKj_col.vector, &KiKl_row.vector, &d);
+	      tr+=d;
+	      gsl_blas_ddot (&Kj_row.vector, &KiKl_row.vector, &d);
+	      tr-=r*d;
+	      gsl_blas_ddot (&Kl_row.vector, &KiKj_row.vector, &d);
+	      tr-=r*d;
+	    }
+	  }
+
+	  tr+=r*r*gsl_vector_get (trKiKj, t_jl);
+	} else if (j!=n_vc && l==n_vc) {
+	  t_ij=GetabIndex (i+1, j+1, n_vc-2);
+	  tr=gsl_vector_get (trKiKjKi, i*n_vc+j)-2*r*gsl_vector_get (trKiKj, t_ij)+r*r*gsl_vector_get (trKi, j);
+	} else if (j==n_vc && l==n_vc) {
+	  t_ii=GetabIndex (i+1, i+1, n_vc-2);
+	  tr=gsl_vector_get (trKiKj, t_ii)-2*r*gsl_vector_get (trKi, i)+r*r*(double)(ni_test-1);
+	}
+
+	gsl_matrix_set (Q, j, i*(n_vc+1)+l, tr);
+	if (l!=j) {gsl_matrix_set (Q, l, i*(n_vc+1)+j, tr);}
+      }
+    }
+  }
+
+  gsl_matrix_scale (Q, 1.0/pow((double)ni_test, 2) );
+
+  gsl_matrix_free(KiKj);
+  gsl_vector_free(trKiKjKi);
+  gsl_vector_free(trKiKj);
+  gsl_vector_free(trKi);
+
+  return;
+}
+
+
+
+//perform Jacknife sampling for variance of S
+void JacknifeGtoS (const gsl_matrix *G, gsl_matrix *S, gsl_matrix *Svar) {
+  size_t n_vc=Svar->size1, ni_test=G->size1;
+  vector<vector<vector<double> > > tr_KiKj, s_KiKj;
+  vector<vector<double> > sum_Ki, s_Ki, si;
+  vector<double> vec_tmp;
+  double di, dj, d, m, v;
+
+  //initialize and set all elements to zero
+  for (size_t i=0; i<ni_test; i++) {
+    vec_tmp.push_back(0);
+  }
+
+  for (size_t i=0; i<n_vc; i++) {
+    sum_Ki.push_back(vec_tmp);
+    s_Ki.push_back(vec_tmp);
+    si.push_back(vec_tmp);
+  }
+
+  for (size_t i=0; i<n_vc; i++) {
+    tr_KiKj.push_back(sum_Ki);
+    s_KiKj.push_back(sum_Ki);
+  }
+
+  //run jacknife
+  for (size_t i=0; i<n_vc; i++) {
+    for (size_t l=0; l<ni_test; l++) {
+      for (size_t k=0; k<ni_test; k++) {
+	di=gsl_matrix_get(G, l, k+ni_test*i);
+
+	for (size_t t=0; t<ni_test; t++) {
+	  if (t==l || t==k) {continue;}
+	  sum_Ki[i][t]+=di;
+	  if (l==k) {si[i][t]+=di;}
+	}
+	s_Ki[i][l]+=di;
+      }
+    }
+
+    for (size_t t=0; t<ni_test; t++) {
+      sum_Ki[i][t]/=(double)(ni_test-1);
+    }
+  }
+
+  for (size_t i=0; i<n_vc; i++) {
+    for (size_t j=i; j<n_vc; j++) {
+      for (size_t l=0; l<ni_test; l++) {
+	for (size_t k=0; k<ni_test; k++) {
+	  di=gsl_matrix_get(G, l, k+ni_test*i);
+	  dj=gsl_matrix_get(G, l, k+ni_test*j);
+	  d=di*dj;
+
+	  for (size_t t=0; t<ni_test; t++) {
+	    if (t==l || t==k) {continue;}
+	    tr_KiKj[i][j][t]+=d;
+          }
+	}
+
+	for (size_t t=0; t<ni_test; t++) {
+	  if (t==l) {continue;}
+	  di=gsl_matrix_get(G, l, t+ni_test*i);
+	  dj=gsl_matrix_get(G, l, t+ni_test*j);
+
+	  s_KiKj[i][j][t]+=(s_Ki[i][l]-di)*(s_Ki[j][l]-dj);
+	}
+      }
+
+      for (size_t t=0; t<ni_test; t++) {
+	s_KiKj[i][j][t]/=(double)(ni_test-1);
+      }
+
+      m=0; v=0;
+      for (size_t t=0; t<ni_test; t++) {
+	d=tr_KiKj[i][j][t]-2*s_KiKj[i][j][t]+sum_Ki[i][t]*sum_Ki[j][t];
+	d/=(si[i][t]-sum_Ki[i][t])*(si[j][t]-sum_Ki[j][t]);
+	d-=1/(double)(ni_test-2);
+
+	m+=d; v+=d*d;
+      }
+      m/=(double)ni_test;
+      v/=(double)ni_test;
+      v-=m*m;
+      v*=(double)(ni_test-1);
+
+      gsl_matrix_set (Svar, i, j, v);
+      d=gsl_matrix_get (S, i, j);
+      d=(double)ni_test*d-(double)(ni_test-1)*m;
+      gsl_matrix_set (S, i, j, d);
+      if (i!=j) {gsl_matrix_set (Svar, j, i, v); gsl_matrix_set (S, j, i, d);}
+    }
+  }
+
+  return;
+}
+
+
+
+//compute the d by d S matrix with its d by d variance matrix of Svar, and the d+1 by d(d+1) matrix of Q for V(q)
+void PARAM::CalcS (gsl_matrix *S, gsl_matrix *Svar, gsl_matrix *Q)  {
+  string file_str;
+
+  gsl_matrix_set_zero (S);
+  gsl_matrix_set_zero (Svar);
+  gsl_matrix_set_zero (Q);
+
+  //compute the kinship matrix G for multiple categories; these matrices are not centered, for convienence of Jacknife sampling
+  gsl_matrix *G=gsl_matrix_alloc (ni_test, n_vc*ni_test);
+  gsl_matrix_set_zero (G);
+
+  if (!file_bfile.empty() ) {
+    file_str=file_bfile+".bed";
+    if (PlinkKin (file_str, indicator_idv, indicator_snp, a_mode-24, d_pace, mapRS2cat, mapRS2var, snpInfo, G)==false) {error=true;}
+  } else {
+    file_str=file_geno;
+    if (BimbamKin (file_str, indicator_idv, indicator_snp, a_mode-24, d_pace, mapRS2cat, mapRS2var, snpInfo, G)==false) {error=true;}
+  }
+
+  //center and scale every kinship matrix inside G
+  double d;
+  for (size_t i=0; i<n_vc; i++) {
+    gsl_matrix_view K=gsl_matrix_submatrix(G, 0, i*ni_test, ni_test, ni_test);
+    CenterMatrix(&K.matrix);
+    d=ScaleMatrix(&K.matrix);
+  }
+
+  //based on G, compute S
+  compKtoS (G, S);
+
+  //based on G, compute a matrix Q that can be used to calculate the variance of q
+  compKtoQ (G, Q);
+
+  /*
+  //set up random environment
+  gsl_rng_env_setup();
+  gsl_rng *gsl_r;
+  const gsl_rng_type * gslType;
+  gslType = gsl_rng_default;
+  if (randseed<0) {
+    time_t rawtime;
+    time (&rawtime);
+    tm * ptm = gmtime (&rawtime);
+
+    randseed = (unsigned) (ptm->tm_hour%24*3600+ptm->tm_min*60+ptm->tm_sec);
+  }
+  gsl_r = gsl_rng_alloc(gslType);
+  gsl_rng_set(gsl_r, randseed);
+
+  //bootstrap: in each iteration, sample individuals and compute S_pmt
+  size_t n_pmt=100;
+  vector<size_t> idv_order, idv_remove;
+  for (size_t i=0; i<ni_test; i++) {
+    idv_order.push_back(i);
+  }
+  for (size_t i=0; i<n_pmt; i++) {
+    idv_remove.push_back(0);
+  }
+  gsl_ran_choose (gsl_r, static_cast<void*>(&idv_remove[0]), n_pmt, static_cast<void*>(&idv_order[0]), ni_test, sizeof(size_t));
+
+  gsl_matrix *S_pmt=gsl_matrix_alloc(n_vc, n_vc*n_pmt);
+  for (size_t i=0; i<n_pmt; i++) {
+    gsl_matrix_view S_sub=gsl_matrix_submatrix (S_pmt, 0, n_vc*i, n_vc, n_vc);
+    compKtoS (G, idv_remove[i], &S_sub.matrix);
+  }
+
+  //based on S_pmt, compute Svar
+  double m, v, d;
+  for (size_t i=0; i<n_vc; i++) {
+    for (size_t j=i; j<n_vc; j++) {
+      m=0; v=0;
+      for (size_t t=0; t<n_pmt; t++) {
+	d=gsl_matrix_get(S_pmt, i, j);
+	m+=d; v+=d*d;
+      }
+      m/=(double)n_pmt; v/=(double)n_pmt;
+      v=v-m*m;
+      gsl_matrix_set(Svar, i, j, v);
+      if (i!=j) {gsl_matrix_set(Svar, j, i, v);}
+    }
+  }
+  */
+
+  //compute Svar and update S with Jacknife
+  JacknifeGtoS (G, S, Svar);
+
+  gsl_matrix_free(G);
+  return;
+}
+
+
+
+void PARAM::WriteVector (const gsl_vector *q, const gsl_vector *s, const size_t n_total, const string suffix)
+{
+	string file_str;
+	file_str=path_out+"/"+file_out;
+	file_str+=".";
+	file_str+=suffix;
+	file_str+=".txt";
+
+	ofstream outfile (file_str.c_str(), ofstream::out);
+	if (!outfile) {cout<<"error writing file: "<<file_str.c_str()<<endl; return;}
+
+	outfile.precision(10);
+
+	for (size_t i=0; i<q->size; ++i) {
+	  outfile<<gsl_vector_get (q, i)<<endl;
+	}
+
+	for (size_t i=0; i<s->size; ++i) {
+	  outfile<<gsl_vector_get (s, i)<<endl;
+	}
+
+	outfile<<n_total<<endl;
+
+	outfile.close();
+	outfile.clear();
+	return;
+}
+
+
+
+void PARAM::WriteVar (const string suffix)
+{
+  string file_str, rs;
+	file_str=path_out+"/"+file_out;
+	file_str+=".";
+	file_str+=suffix;
+	file_str+=".txt.gz";
+
+	ogzstream outfile (file_str.c_str(), ogzstream::out);
+	if (!outfile) {cout<<"error writing file: "<<file_str.c_str()<<endl; return;}
+
+	outfile.precision(10);
+
+	for (size_t i=0; i<indicator_snp.size(); i++) {
+	  if (indicator_snp[i]==0) {continue;}
+	  rs=snpInfo[i].rs_number;
+	  if (mapRS2var.count(rs)!=0) {
+	    outfile<<rs<<"\t"<<mapRS2var.at(rs)<<endl;
+	  }
+	}
+
+	outfile.close();
+	outfile.clear();
 	return;
 }
-		
 
 
 
 
-void PARAM::WriteMatrix (const gsl_matrix *matrix_U, const string suffix) 
+void PARAM::WriteMatrix (const gsl_matrix *matrix_U, const string suffix)
 {
 	string file_str;
 	file_str=path_out+"/"+file_out;
 	file_str+=".";
 	file_str+=suffix;
-	file_str+=".txt";	
-	
+	file_str+=".txt";
+
 	ofstream outfile (file_str.c_str(), ofstream::out);
 	if (!outfile) {cout<<"error writing file: "<<file_str.c_str()<<endl; return;}
-	
+
 	outfile.precision(10);
-	
+
 	for (size_t i=0; i<matrix_U->size1; ++i) {
 		for (size_t j=0; j<matrix_U->size2; ++j) {
 			outfile<<gsl_matrix_get (matrix_U, i, j)<<"\t";
 		}
 		outfile<<endl;
 	}
-	
+
 	outfile.close();
 	outfile.clear();
 	return;
 }
 
 
-void PARAM::WriteVector (const gsl_vector *vector_D, const string suffix) 
+void PARAM::WriteVector (const gsl_vector *vector_D, const string suffix)
 {
 	string file_str;
 	file_str=path_out+"/"+file_out;
 	file_str+=".";
 	file_str+=suffix;
 	file_str+=".txt";
-	
+
 	ofstream outfile (file_str.c_str(), ofstream::out);
 	if (!outfile) {cout<<"error writing file: "<<file_str.c_str()<<endl; return;}
-	
+
 	outfile.precision(10);
-	
+
 	for (size_t i=0; i<vector_D->size; ++i) {
 		outfile<<gsl_vector_get (vector_D, i)<<endl;
 	}
-	
+
 	outfile.close();
 	outfile.clear();
 	return;
 }
 
 
-void PARAM::CheckCvt () 
+void PARAM::CheckCvt ()
 {
 	if (indicator_cvt.size()==0) {return;}
-		
+
 	size_t ci_test=0;
-	
+
 	gsl_matrix *W=gsl_matrix_alloc (ni_test, n_cvt);
-	
+
 	for (vector<int>::size_type i=0; i<indicator_idv.size(); ++i) {
 		if (indicator_idv[i]==0 || indicator_cvt[i]==0) {continue;}
 		for (size_t j=0; j<n_cvt; ++j) {
@@ -679,14 +1308,14 @@ void PARAM::CheckCvt ()
 	size_t flag_ipt=0;
 	double v_min, v_max;
 	set<size_t> set_remove;
-	
+
 	//check if any columns is an intercept
 	for (size_t i=0; i<W->size2; i++) {
 		gsl_vector_view w_col=gsl_matrix_column (W, i);
 		gsl_vector_minmax (&w_col.vector, &v_min, &v_max);
 		if (v_min==v_max) {flag_ipt=1; set_remove.insert (i);}
 	}
-	
+
 	//add an intecept term if needed
 	if (n_cvt==set_remove.size()) {
 		indicator_cvt.clear();
@@ -697,19 +1326,19 @@ void PARAM::CheckCvt ()
 			if (indicator_idv[i]==0 || indicator_cvt[i]==0) {continue;}
 			cvt[i].push_back(1.0);
 		}
-		
+
 		n_cvt++;
-	} else {}	
-	
+	} else {}
+
 	gsl_matrix_free(W);
-	
+
 	return;
 }
 
 
 //post-process phentoypes, covariates
 void PARAM::ProcessCvtPhen ()
-{	
+{
 	//convert indicator_pheno to indicator_idv
 	int k=1;
 	indicator_idv.clear();
@@ -720,27 +1349,88 @@ void PARAM::ProcessCvtPhen ()
 		}
 		indicator_idv.push_back(k);
 	}
-	
+
 	//remove individuals with missing covariates
 	if ((indicator_cvt).size()!=0) {
 		for (vector<int>::size_type i=0; i<(indicator_idv).size(); ++i) {
 			indicator_idv[i]*=indicator_cvt[i];
 		}
 	}
-	
+
+	//remove individuals with missing gxe variables
+	if ((indicator_gxe).size()!=0) {
+		for (vector<int>::size_type i=0; i<(indicator_idv).size(); ++i) {
+			indicator_idv[i]*=indicator_gxe[i];
+		}
+	}
+
+	//remove individuals with missing residual weights
+	if ((indicator_weight).size()!=0) {
+		for (vector<int>::size_type i=0; i<(indicator_idv).size(); ++i) {
+			indicator_idv[i]*=indicator_weight[i];
+		}
+	}
+
 	//obtain ni_test
-	ni_test=0; 
+	ni_test=0;
 	for (vector<int>::size_type i=0; i<(indicator_idv).size(); ++i) {
-		if (indicator_idv[i]==0) {continue;}
+	    if (indicator_idv[i]==0) {continue;}
 		ni_test++;
 	}
-	
+
+
+
+	//if subsample number is set, perform a random sub-sampling to determine the subsampled ids
+	if (ni_subsample!=0) {
+	  if (ni_test<ni_subsample) {
+	    cout<<"error! number of subsamples is less than number of analyzed individuals. "<<endl;
+	  } else {
+	    //set up random environment
+	    gsl_rng_env_setup();
+	    gsl_rng *gsl_r;
+	    const gsl_rng_type * gslType;
+	    gslType = gsl_rng_default;
+	    if (randseed<0) {
+	      time_t rawtime;
+	      time (&rawtime);
+	      tm * ptm = gmtime (&rawtime);
+
+	      randseed = (unsigned) (ptm->tm_hour%24*3600+ptm->tm_min*60+ptm->tm_sec);
+	    }
+	    gsl_r = gsl_rng_alloc(gslType);
+	    gsl_rng_set(gsl_r, randseed);
+
+	    //from ni_test, sub-sample ni_subsample
+	    vector<size_t> a, b;
+	    for (size_t i=0; i<ni_subsample; i++) {
+              a.push_back(0);
+	    }
+	    for (size_t i=0; i<ni_test; i++) {
+	      b.push_back(i);
+	    }
+
+	    gsl_ran_choose (gsl_r, static_cast<void*>(&a[0]), ni_subsample, static_cast<void*>(&b[0]), ni_test, sizeof (size_t) );
+
+	    //re-set indicator_idv and ni_test
+	    int j=0;
+	    for (vector<int>::size_type i=0; i<(indicator_idv).size(); ++i) {
+	      if (indicator_idv[i]==0) {continue;}
+	      if(find(a.begin(), a.end(), j) == a.end()) {
+		indicator_idv[i]=0;
+	      }
+	      j++;
+	    }
+	    ni_test=ni_subsample;
+	  }
+	}
+
+	//check ni_test
 	if (ni_test==0) {
 		error=true;
 		cout<<"error! number of analyzed individuals equals 0. "<<endl;
 		return;
 	}
-	
+
 	//check covariates to see if they are correlated with each other, and to see if the intercept term is included
 	//after getting ni_test
 	//add or remove covariates
@@ -749,24 +1439,24 @@ void PARAM::ProcessCvtPhen ()
 	} else {
 		vector<double> cvt_row;
 		cvt_row.push_back(1);
-		
+
 		for (vector<int>::size_type i=0; i<(indicator_idv).size(); ++i) {
 			indicator_cvt.push_back(1);
-			
+
 			cvt.push_back(cvt_row);
 		}
 	}
-	 
+
 	return;
 }
 
 
 
 
-void PARAM::CopyCvt (gsl_matrix *W) 
+void PARAM::CopyCvt (gsl_matrix *W)
 {
 	size_t ci_test=0;
-	
+
 	for (vector<int>::size_type i=0; i<indicator_idv.size(); ++i) {
 		if (indicator_idv[i]==0 || indicator_cvt[i]==0) {continue;}
 		for (size_t j=0; j<n_cvt; ++j) {
@@ -774,57 +1464,85 @@ void PARAM::CopyCvt (gsl_matrix *W)
 		}
 		ci_test++;
 	}
-	
+
+	return;
+}
+
+
+void PARAM::CopyGxe (gsl_vector *env)
+{
+	size_t ci_test=0;
+
+	for (vector<int>::size_type i=0; i<indicator_idv.size(); ++i) {
+		if (indicator_idv[i]==0 || indicator_gxe[i]==0) {continue;}
+		gsl_vector_set (env, ci_test, gxe[i]);
+		ci_test++;
+	}
+
+	return;
+}
+
+void PARAM::CopyWeight (gsl_vector *w)
+{
+	size_t ci_test=0;
+
+	for (vector<int>::size_type i=0; i<indicator_idv.size(); ++i) {
+		if (indicator_idv[i]==0 || indicator_weight[i]==0) {continue;}
+		gsl_vector_set (w, ci_test, weight[i]);
+		ci_test++;
+	}
+
 	return;
 }
 
 
 //if flag=0, then use indicator_idv to load W and Y
 //else, use indicator_cvt to load them
-void PARAM::CopyCvtPhen (gsl_matrix *W, gsl_vector *y, size_t flag) 
+void PARAM::CopyCvtPhen (gsl_matrix *W, gsl_vector *y, size_t flag)
 {
 	size_t ci_test=0;
-	
+
 	for (vector<int>::size_type i=0; i<indicator_idv.size(); ++i) {
 		if (flag==0) {
 			if (indicator_idv[i]==0) {continue;}
 		} else {
 			if (indicator_cvt[i]==0) {continue;}
 		}
-		
+
 		gsl_vector_set (y, ci_test, (pheno)[i][0]);
-		
+
 		for (size_t j=0; j<n_cvt; ++j) {
 			gsl_matrix_set (W, ci_test, j, (cvt)[i][j]);
 		}
 		ci_test++;
 	}
-	
+
 	return;
 }
 
 //if flag=0, then use indicator_idv to load W and Y
 //else, use indicator_cvt to load them
-void PARAM::CopyCvtPhen (gsl_matrix *W, gsl_matrix *Y, size_t flag) 
+void PARAM::CopyCvtPhen (gsl_matrix *W, gsl_matrix *Y, size_t flag)
 {
 	size_t ci_test=0;
-	
+
 	for (vector<int>::size_type i=0; i<indicator_idv.size(); ++i) {
 		if (flag==0) {
 			if (indicator_idv[i]==0) {continue;}
 		} else {
 			if (indicator_cvt[i]==0) {continue;}
-		}		
-		
-		for (size_t j=0; j<n_ph; ++j) {
+		}
+
+        for (size_t j=0; j<n_ph; ++j) {
 			gsl_matrix_set (Y, ci_test, j, (pheno)[i][j]);
 		}
 		for (size_t j=0; j<n_cvt; ++j) {
 			gsl_matrix_set (W, ci_test, j, (cvt)[i][j]);
 		}
+
 		ci_test++;
 	}
-	
+
 	return;
 }
 
@@ -832,18 +1550,18 @@ void PARAM::CopyCvtPhen (gsl_matrix *W, gsl_matrix *Y, size_t flag)
 
 
 
-void PARAM::CopyRead (gsl_vector *log_N) 
+void PARAM::CopyRead (gsl_vector *log_N)
 {
 	size_t ci_test=0;
-	
+
 	for (vector<int>::size_type i=0; i<indicator_idv.size(); ++i) {
 		if (indicator_idv[i]==0) {continue;}
-		gsl_vector_set (log_N, ci_test, log(vec_read[i]) );	
+		gsl_vector_set (log_N, ci_test, log(vec_read[i]) );
 		ci_test++;
 	}
-	
+
 	return;
 }
-		
-		
+
+