From 17deca2d54827a00df3ea4d98df700fc2b8ed777 Mon Sep 17 00:00:00 2001
From: xiangzhou
Date: Sat, 20 Sep 2014 10:17:34 -0400
Subject: initial upload, version 0.95alpha

---
 io.cpp | 1396 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 1396 insertions(+)
 create mode 100644 io.cpp

(limited to 'io.cpp')

diff --git a/io.cpp b/io.cpp
new file mode 100644
index 0000000..c22f668
--- /dev/null
+++ b/io.cpp
@@ -0,0 +1,1396 @@
+/*
+	Genome-wide Efficient Mixed Model Association (GEMMA)
+    Copyright (C) 2011  Xiang Zhou
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <string>
+#include <iomanip>
+#include <bitset>
+#include <vector>
+#include <map>
+#include <set>
+#include <cstring>
+#include <cmath>
+#include <stdio.h>
+#include <stdlib.h> 
+
+#include "gsl/gsl_vector.h"
+#include "gsl/gsl_matrix.h"
+#include "gsl/gsl_linalg.h"
+#include "gsl/gsl_blas.h"
+#include "gsl/gsl_cdf.h"
+
+#include "lapack.h"
+#include "gzstream.h"
+#include "mathfunc.h"
+
+#ifdef FORCE_FLOAT
+#include "io_float.h"
+#else
+#include "io.h"
+#endif
+
+
+using namespace std;
+
+
+
+//Print process bar
+void ProgressBar (string str, double p, double total)
+{
+	double progress = (100.0 * p / total); 
+	int barsize = (int) (progress / 2.0); 
+	char bar[51];
+	
+	cout<<str;
+	for (int i = 0; i <50; i++) {
+		if (i<barsize) {bar[i] = '=';}
+		else {bar[i]=' ';}
+		cout<<bar[i];
+	}
+	cout<<setprecision(2)<<fixed<<progress<<"%\r"<<flush;
+	
+	return;
+}
+
+
+//Print process bar (with acceptance ratio)
+void ProgressBar (string str, double p, double total, double ratio)
+{
+	double progress = (100.0 * p / total); 
+	int barsize = (int) (progress / 2.0); 
+	char bar[51];
+	
+	cout<<str;
+	for (int i = 0; i <50; i++) {
+		if (i<barsize) {bar[i] = '=';}
+		else {bar[i]=' ';}
+		cout<<bar[i];
+	}
+	cout<<setprecision(2)<<fixed<<progress<<"%    "<<ratio<<"\r"<<flush;
+	
+	
+	return;
+}
+
+// in case files are ended with "\r" or "\r\n"
+std::istream& safeGetline(std::istream& is, std::string& t)
+{
+    t.clear();
+
+    // The characters in the stream are read one-by-one using a std::streambuf.
+    // That is faster than reading them one-by-one using the std::istream.
+    // Code that uses streambuf this way must be guarded by a sentry object.
+    // The sentry object performs various tasks,
+    // such as thread synchronization and updating the stream state.
+
+    std::istream::sentry se(is, true);
+    std::streambuf* sb = is.rdbuf();
+
+    for(;;) {
+        int c = sb->sbumpc();
+        switch (c) {
+        case '\n':
+            return is;
+        case '\r':
+            if(sb->sgetc() == '\n')
+                sb->sbumpc();
+            return is;
+        case EOF:
+            // Also handle the case when the last line has no line ending
+            if(t.empty())
+                is.setstate(std::ios::eofbit);
+            return is;
+        default:
+            t += (char)c;
+        }
+    }
+}
+
+//Read snp file
+bool ReadFile_snps (const string &file_snps, set<string> &setSnps)
+{
+	setSnps.clear();
+
+	ifstream infile (file_snps.c_str(), ifstream::in);
+	if (!infile) {cout<<"error! fail to open snps file: "<<file_snps<<endl; return false;}
+	
+	string line;
+	char *ch_ptr;
+	
+	while (getline(infile, line)) {
+		ch_ptr=strtok ((char *)line.c_str(), " , \t");
+		setSnps.insert(ch_ptr); 
+	}
+	
+	infile.close();
+	infile.clear();	
+	
+	return true;
+}
+
+
+//Read log file
+bool ReadFile_log (const string &file_log, double &pheno_mean)
+{
+	ifstream infile (file_log.c_str(), ifstream::in);
+	if (!infile) {cout<<"error! fail to open log file: "<<file_log<<endl; return false;}
+	
+	string line;
+	char *ch_ptr;
+	size_t flag=0;
+	
+	while (getline(infile, line)) {
+		ch_ptr=strtok ((char *)line.c_str(), " , \t");
+		ch_ptr=strtok (NULL, " , \t");
+		
+		if (ch_ptr!=NULL && strcmp(ch_ptr, "estimated")==0) {
+			ch_ptr=strtok (NULL, " , \t");
+			if (ch_ptr!=NULL && strcmp(ch_ptr, "mean")==0) {
+				ch_ptr=strtok (NULL, " , \t");
+				if (ch_ptr!=NULL && strcmp(ch_ptr, "=")==0) {
+					ch_ptr=strtok (NULL, " , \t");
+					pheno_mean=atof(ch_ptr);
+					flag=1;
+				}
+			}
+		}
+		
+		if (flag==1) {break;}
+	}
+	
+	infile.close();
+	infile.clear();	
+	
+	return true;
+}
+
+
+//Read bimbam annotation file
+bool ReadFile_anno (const string &file_anno, map<string, string> &mapRS2chr, map<string, long int> &mapRS2bp, map<string, double> &mapRS2cM)
+{
+	mapRS2chr.clear();
+	mapRS2bp.clear();
+	
+	ifstream infile (file_anno.c_str(), ifstream::in);
+	if (!infile) {cout<<"error opening annotation file: "<<file_anno<<endl; return false;}
+	
+	string line;
+	char *ch_ptr;
+	
+	string rs;
+	long int b_pos;
+	string chr;
+	double cM;
+	
+	while (!safeGetline(infile, line).eof()) {
+		ch_ptr=strtok ((char *)line.c_str(), " , \t");
+		rs=ch_ptr;
+		ch_ptr=strtok (NULL, " , \t");
+		if (strcmp(ch_ptr, "NA")==0) {b_pos=-9;} else {b_pos=atol(ch_ptr);}
+		ch_ptr=strtok (NULL, " , \t");
+		if (ch_ptr==NULL || strcmp(ch_ptr, "NA")==0) {chr="-9";} else {chr=ch_ptr;}
+		ch_ptr=strtok (NULL, " , \t");
+		if (ch_ptr==NULL || strcmp(ch_ptr, "NA")==0) {cM=-9;} else {cM=atof(ch_ptr);}
+		
+		mapRS2chr[rs]=chr;
+		mapRS2bp[rs]=b_pos;
+		mapRS2cM[rs]=cM;
+	}
+	
+	infile.close();
+	infile.clear();	
+	
+	return true;
+}
+
+//read one column of phenotype
+bool ReadFile_column (const string &file_pheno, vector<int> &indicator_idv, vector<double> &pheno, const int &p_column)
+{
+	indicator_idv.clear();
+	pheno.clear();
+	
+	igzstream infile (file_pheno.c_str(), igzstream::in);
+//	ifstream infile (file_pheno.c_str(), ifstream::in);
+	if (!infile) {cout<<"error! fail to open phenotype file: "<<file_pheno<<endl; return false;}
+	
+	string line;
+	char *ch_ptr;
+	
+	string id;
+	double p;
+	while (!safeGetline(infile, line).eof()) {
+		ch_ptr=strtok ((char *)line.c_str(), " , \t");
+		for (int i=0; i<(p_column-1); ++i) {
+			ch_ptr=strtok (NULL, " , \t");	
+		}		
+		if (strcmp(ch_ptr, "NA")==0) {indicator_idv.push_back(0); pheno.push_back(-9);}		//pheno is different from pimass2
+		else {p=atof(ch_ptr); indicator_idv.push_back(1); pheno.push_back(p);}
+	}
+	
+	infile.close();
+	infile.clear();	
+	
+	return true;
+}
+
+
+
+//Read bimbam phenotype file, p_column=1, 2 ...
+bool ReadFile_pheno (const string &file_pheno, vector<vector<int> > &indicator_pheno, vector<vector<double> > &pheno, const vector<size_t> &p_column)
+{
+	indicator_pheno.clear();
+	pheno.clear();
+	
+	igzstream infile (file_pheno.c_str(), igzstream::in);
+//	ifstream infile (file_pheno.c_str(), ifstream::in);
+	if (!infile) {cout<<"error! fail to open phenotype file: "<<file_pheno<<endl; return false;}
+
+	string line;
+	char *ch_ptr;
+  
+	string id;
+	double p;
+	
+	vector<double> pheno_row;
+	vector<int> ind_pheno_row;
+	
+	size_t p_max=*max_element(p_column.begin(), p_column.end() );
+	map<size_t, size_t> mapP2c;
+	for (size_t i=0; i<p_column.size(); i++) {
+		mapP2c[p_column[i]]=i;
+		pheno_row.push_back(-9);
+		ind_pheno_row.push_back(0);
+	}	
+	
+	while (!safeGetline(infile, line).eof()) {
+		ch_ptr=strtok ((char *)line.c_str(), " , \t");
+		
+		size_t i=0;
+		while (i<p_max ) {			
+			if (mapP2c.count(i+1)!=0) {
+				if (strcmp(ch_ptr, "NA")==0) {ind_pheno_row[mapP2c[i+1]]=0; pheno_row[mapP2c[i+1]]=-9;}
+				else {p=atof(ch_ptr); ind_pheno_row[mapP2c[i+1]]=1; pheno_row[mapP2c[i+1]]=p;}
+			}
+			i++;
+			ch_ptr=strtok (NULL, " , \t");	
+		}
+		
+		indicator_pheno.push_back(ind_pheno_row);	
+		pheno.push_back(pheno_row);			
+	}
+ 
+	infile.close();
+	infile.clear();	
+	
+	return true;
+}
+
+
+bool ReadFile_cvt (const string &file_cvt, vector<int> &indicator_cvt, vector<vector<double> > &cvt, size_t &n_cvt)
+{
+	indicator_cvt.clear();
+	
+	ifstream infile (file_cvt.c_str(), ifstream::in);
+	if (!infile) {cout<<"error! fail to open covariates file: "<<file_cvt<<endl; return false;}
+	
+	string line;
+	char *ch_ptr;
+	double d;	
+	
+	int flag_na=0;	
+	
+	while (!safeGetline(infile, line).eof()) {
+		vector<double> v_d; flag_na=0;
+		ch_ptr=strtok ((char *)line.c_str(), " , \t");
+		while (ch_ptr!=NULL) {
+			if (strcmp(ch_ptr, "NA")==0) {flag_na=1; d=-9;}
+			else {d=atof(ch_ptr);}
+			
+			v_d.push_back(d);
+			ch_ptr=strtok (NULL, " , \t");	
+		}
+		if (flag_na==0) {indicator_cvt.push_back(1);} else {indicator_cvt.push_back(0);} 
+		cvt.push_back(v_d);
+	}
+	
+	if (indicator_cvt.empty()) {n_cvt=0;}
+	else {
+		flag_na=0;
+		for (vector<int>::size_type i=0; i<indicator_cvt.size(); ++i) {
+			if (indicator_cvt[i]==0) {continue;}
+			
+			if (flag_na==0) {flag_na=1; n_cvt=cvt[i].size();}
+			if (flag_na!=0 && n_cvt!=cvt[i].size()) {cout<<"error! number of covariates in row "<<i<<" do not match other rows."<<endl; return false;}
+		}
+	}
+	
+	infile.close();
+	infile.clear();	
+	
+	return true;
+}
+
+
+
+//Read .bim file
+bool ReadFile_bim (const string &file_bim, vector<SNPINFO> &snpInfo)
+{
+	snpInfo.clear();
+	
+	ifstream infile (file_bim.c_str(), ifstream::in);
+	if (!infile) {cout<<"error opening .bim file: "<<file_bim<<endl; return false;}
+	
+	string line;
+	char *ch_ptr;
+	
+	string rs;
+	long int b_pos;
+	string chr;
+	double cM;
+	string major;
+	string minor;
+	
+	while (getline(infile, line)) {
+		ch_ptr=strtok ((char *)line.c_str(), " \t");
+		chr=ch_ptr;
+		ch_ptr=strtok (NULL, " \t");
+		rs=ch_ptr;
+		ch_ptr=strtok (NULL, " \t");
+		cM=atof(ch_ptr);
+		ch_ptr=strtok (NULL, " \t");
+		b_pos=atol(ch_ptr);
+		ch_ptr=strtok (NULL, " \t");
+		minor=ch_ptr;
+		ch_ptr=strtok (NULL, " \t");
+		major=ch_ptr;
+		
+		SNPINFO sInfo={chr, rs, cM, b_pos, minor, major, -9, -9, -9};
+		snpInfo.push_back(sInfo);
+	}
+	
+	infile.close();
+	infile.clear();	
+	return true;
+}
+
+
+//Read .fam file
+bool ReadFile_fam (const string &file_fam, vector<vector<int> > &indicator_pheno, vector<vector<double> > &pheno, map<string, int> &mapID2num, const vector<size_t> &p_column)
+{
+	indicator_pheno.clear();
+	pheno.clear();
+	mapID2num.clear();	
+	
+	igzstream infile (file_fam.c_str(), igzstream::in);
+	//ifstream infile (file_fam.c_str(), ifstream::in);
+	if (!infile) {cout<<"error opening .fam file: "<<file_fam<<endl; return false;}
+
+	string line;
+	char *ch_ptr;
+
+	string id;
+	int c=0;
+	double p;
+
+	vector<double> pheno_row;
+	vector<int> ind_pheno_row;
+	
+	size_t p_max=*max_element(p_column.begin(), p_column.end() );
+	map<size_t, size_t> mapP2c;
+	for (size_t i=0; i<p_column.size(); i++) {
+		mapP2c[p_column[i]]=i;
+		pheno_row.push_back(-9);
+		ind_pheno_row.push_back(0);
+	}	
+	
+	while (!safeGetline(infile, line).eof()) {
+		ch_ptr=strtok ((char *)line.c_str(), " \t");
+		ch_ptr=strtok (NULL, " \t");
+		id=ch_ptr;
+		ch_ptr=strtok (NULL, " \t");
+		ch_ptr=strtok (NULL, " \t");
+		ch_ptr=strtok (NULL, " \t");
+		ch_ptr=strtok (NULL, " \t");
+		
+		size_t i=0;
+		while (i<p_max ) {
+			if (mapP2c.count(i+1)!=0 ) {
+				if (strcmp(ch_ptr, "NA")==0) {
+					ind_pheno_row[mapP2c[i+1]]=0; pheno_row[mapP2c[i+1]]=-9;
+				} else {
+					p=atof(ch_ptr);
+					
+					if (p==-9) {ind_pheno_row[mapP2c[i+1]]=0; pheno_row[mapP2c[i+1]]=-9;}
+					else {ind_pheno_row[mapP2c[i+1]]=1; pheno_row[mapP2c[i+1]]=p;}
+				}
+			}
+			i++;
+			ch_ptr=strtok (NULL, " , \t");	
+		}
+		
+		indicator_pheno.push_back(ind_pheno_row);
+		pheno.push_back(pheno_row);				
+		
+		mapID2num[id]=c; c++;
+	}
+ 
+	infile.close();
+	infile.clear();	
+	return true;
+}
+
+
+
+
+
+
+//Read bimbam mean genotype file, the first time, to obtain #SNPs for analysis (ns_test) and total #SNP (ns_total)
+bool ReadFile_geno (const string &file_geno, const set<string> &setSnps, const gsl_matrix *W, vector<int> &indicator_idv, vector<int> &indicator_snp, const double &maf_level, const double &miss_level, const double &hwe_level, const double &r2_level, map<string, string> &mapRS2chr, map<string, long int> &mapRS2bp, map<string, double> &mapRS2cM, vector<SNPINFO> &snpInfo, size_t &ns_test)
+{
+	indicator_snp.clear();
+	snpInfo.clear();
+	
+	igzstream infile (file_geno.c_str(), igzstream::in);
+//	ifstream infile (file_geno.c_str(), ifstream::in);
+	if (!infile) {cout<<"error reading genotype file:"<<file_geno<<endl; return false;}
+
+	gsl_vector *genotype=gsl_vector_alloc (W->size1);
+	gsl_vector *genotype_miss=gsl_vector_alloc (W->size1);
+	gsl_matrix *WtW=gsl_matrix_alloc (W->size2, W->size2);
+	gsl_matrix *WtWi=gsl_matrix_alloc (W->size2, W->size2);
+	gsl_vector *Wtx=gsl_vector_alloc (W->size2);
+	gsl_vector *WtWiWtx=gsl_vector_alloc (W->size2);
+	gsl_permutation * pmt=gsl_permutation_alloc (W->size2);
+	
+	gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, W, W, 0.0, WtW);
+	int sig;	
+	LUDecomp (WtW, pmt, &sig);
+	LUInvert (WtW, pmt, WtWi);
+	
+	double v_x, v_w;
+	int c_idv=0;
+	
+	string line;
+	char *ch_ptr;
+		
+	string rs;
+	long int b_pos;
+	string chr;
+	string major;
+	string minor;
+	double cM;
+  
+	double maf, geno, geno_old;
+	size_t n_miss;
+	size_t n_0, n_1, n_2;
+	int flag_poly;
+	
+	int ni_total=indicator_idv.size();
+	int ni_test=0;
+	for (int i=0; i<ni_total; ++i) {
+		ni_test+=indicator_idv[i];
+	}
+	ns_test=0;
+	
+	while (!safeGetline(infile, line).eof()) {		
+		ch_ptr=strtok ((char *)line.c_str(), " , \t");
+		rs=ch_ptr;
+		ch_ptr=strtok (NULL, " , \t");
+		minor=ch_ptr;
+		ch_ptr=strtok (NULL, " , \t");
+		major=ch_ptr;
+		
+		if (setSnps.size()!=0 && setSnps.count(rs)==0) {
+			SNPINFO sInfo={"-9", rs, -9, -9, minor, major, -9, -9, -9};
+			snpInfo.push_back(sInfo);
+			indicator_snp.push_back(0);
+			continue;
+		}
+				
+		if (mapRS2bp.count(rs)==0) {chr="-9"; b_pos=-9;cM=-9;}
+		else {b_pos=mapRS2bp[rs]; chr=mapRS2chr[rs]; cM=mapRS2cM[rs];}		
+				
+		maf=0; n_miss=0; flag_poly=0; geno_old=-9;
+		n_0=0; n_1=0; n_2=0;
+		c_idv=0; gsl_vector_set_zero (genotype_miss);
+		for (int i=0; i<ni_total; ++i) {
+			ch_ptr=strtok (NULL, " , \t");
+			if (indicator_idv[i]==0) {continue;}		
+
+			if (strcmp(ch_ptr, "NA")==0) {gsl_vector_set (genotype_miss, c_idv, 1); n_miss++; c_idv++; continue;}
+			
+			geno=atof(ch_ptr);
+			if (geno>=0 && geno<=0.5) {n_0++;}
+			if (geno>0.5 && geno<1.5) {n_1++;}
+			if (geno>=1.5 && geno<=2.0) {n_2++;}
+			
+			gsl_vector_set (genotype, c_idv, geno); 
+			
+//			if (geno<0) {n_miss++; continue;}
+			
+			if (flag_poly==0) {geno_old=geno; flag_poly=2;}
+			if (flag_poly==2 && geno!=geno_old) {flag_poly=1;}
+			
+			maf+=geno;
+			
+			c_idv++;
+		}
+		maf/=2.0*(double)(ni_test-n_miss);	
+		
+		SNPINFO sInfo={chr, rs, cM, b_pos, minor, major, n_miss, (double)n_miss/(double)ni_test, maf};
+		snpInfo.push_back(sInfo);
+		
+		if ( (double)n_miss/(double)ni_test > miss_level) {indicator_snp.push_back(0); continue;}
+		
+		if ( (maf<maf_level || maf> (1.0-maf_level)) && maf_level!=-1 ) {indicator_snp.push_back(0); continue;}
+		
+		if (flag_poly!=1) {indicator_snp.push_back(0); continue;}
+		
+		if (hwe_level!=0) {
+			if (CalcHWE(n_0, n_2, n_1)<hwe_level) {indicator_snp.push_back(0); continue;}
+		}
+		
+		//filter SNP if it is correlated with W
+		for (size_t i=0; i<genotype->size; ++i) {			
+			if (gsl_vector_get (genotype_miss, i)==1) {geno=maf*2.0; gsl_vector_set (genotype, i, geno);}		
+		}
+		
+		gsl_blas_dgemv (CblasTrans, 1.0, W, genotype, 0.0, Wtx);
+		gsl_blas_dgemv (CblasNoTrans, 1.0, WtWi, Wtx, 0.0, WtWiWtx);
+		gsl_blas_ddot (genotype, genotype, &v_x);
+		gsl_blas_ddot (Wtx, WtWiWtx, &v_w);
+		
+		if (v_w/v_x >= r2_level) {indicator_snp.push_back(0); continue;}
+		
+		indicator_snp.push_back(1); 
+		ns_test++;
+	}
+	
+	gsl_vector_free (genotype);
+	gsl_vector_free (genotype_miss);
+	gsl_matrix_free (WtW);
+	gsl_matrix_free (WtWi);
+	gsl_vector_free (Wtx);
+	gsl_vector_free (WtWiWtx);
+	gsl_permutation_free (pmt);
+	
+	infile.close();
+	infile.clear();	
+	
+	return true;
+}
+
+
+
+
+
+      
+//Read bed file, the first time
+bool ReadFile_bed (const string &file_bed, const set<string> &setSnps, const gsl_matrix *W, vector<int> &indicator_idv, vector<int> &indicator_snp, vector<SNPINFO> &snpInfo, const double &maf_level, const double &miss_level, const double &hwe_level, const double &r2_level, size_t &ns_test)
+{
+	indicator_snp.clear();
+	size_t ns_total=snpInfo.size();
+	
+	ifstream infile (file_bed.c_str(), ios::binary);
+	if (!infile) {cout<<"error reading bed file:"<<file_bed<<endl; return false;}
+
+	gsl_vector *genotype=gsl_vector_alloc (W->size1);
+	gsl_vector *genotype_miss=gsl_vector_alloc (W->size1);
+	gsl_matrix *WtW=gsl_matrix_alloc (W->size2, W->size2);
+	gsl_matrix *WtWi=gsl_matrix_alloc (W->size2, W->size2);
+	gsl_vector *Wtx=gsl_vector_alloc (W->size2);
+	gsl_vector *WtWiWtx=gsl_vector_alloc (W->size2);
+	gsl_permutation * pmt=gsl_permutation_alloc (W->size2);
+	
+	gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, W, W, 0.0, WtW);
+	int sig;	
+	LUDecomp (WtW, pmt, &sig);
+	LUInvert (WtW, pmt, WtWi);
+	
+	double v_x, v_w, geno;
+	size_t c_idv=0;
+	
+	char ch[1];
+	bitset<8> b;
+  	
+	size_t ni_total=indicator_idv.size();
+	size_t ni_test=0;
+	for (size_t i=0; i<ni_total; ++i) {
+		ni_test+=indicator_idv[i];
+	}
+	ns_test=0;
+	
+	//calculate n_bit and c, the number of bit for each snp
+	size_t n_bit;
+	if (ni_total%4==0) {n_bit=ni_total/4;}
+	else {n_bit=ni_total/4+1;}
+
+	//ignore the first three majic numbers
+	for (int i=0; i<3; ++i) {
+		infile.read(ch,1);
+		b=ch[0];
+	}
+	
+	double maf;
+	size_t n_miss;
+	size_t n_0, n_1, n_2, c;	
+	
+	//start reading snps and doing association test
+	for (size_t t=0; t<ns_total; ++t) {
+		infile.seekg(t*n_bit+3);		//n_bit, and 3 is the number of magic numbers
+		
+		if (setSnps.size()!=0 && setSnps.count(snpInfo[t].rs_number)==0) {
+			snpInfo[t].n_miss=-9;
+			snpInfo[t].missingness=-9;
+			snpInfo[t].maf=-9;
+			indicator_snp.push_back(0);
+			continue;
+		}
+
+		//read genotypes
+		c=0; maf=0.0; n_miss=0; n_0=0; n_1=0; n_2=0;
+		c_idv=0; gsl_vector_set_zero (genotype_miss);
+		for (size_t i=0; i<n_bit; ++i) {
+			infile.read(ch,1);
+			b=ch[0];
+			for (size_t j=0; j<4; ++j) {                //minor allele homozygous: 2.0; major: 0.0;
+				if ((i==(n_bit-1)) && c==ni_total) {break;}
+				if (indicator_idv[c]==0) {c++; continue;}
+				c++;
+				
+				if (b[2*j]==0) {
+					if (b[2*j+1]==0) {gsl_vector_set(genotype, c_idv, 2.0); maf+=2.0; n_2++;}
+					else {gsl_vector_set(genotype, c_idv, 1.0); maf+=1.0; n_1++;}
+				}
+				else {
+					if (b[2*j+1]==1) {gsl_vector_set(genotype, c_idv, 0.0); maf+=0.0; n_0++;}                                  
+					else {gsl_vector_set(genotype_miss, c_idv, 1); n_miss++; }
+				}
+				c_idv++;
+			}
+		}
+		maf/=2.0*(double)(ni_test-n_miss);
+		
+		snpInfo[t].n_miss=n_miss;
+		snpInfo[t].missingness=(double)n_miss/(double)ni_test;
+		snpInfo[t].maf=maf;
+		
+		if ( (double)n_miss/(double)ni_test > miss_level) {indicator_snp.push_back(0); continue;}
+		
+		if ( (maf<maf_level || maf> (1.0-maf_level)) && maf_level!=-1 ) {indicator_snp.push_back(0); continue;}
+		
+		if ( (n_0+n_1)==0 || (n_1+n_2)==0 || (n_2+n_0)==0) {indicator_snp.push_back(0); continue;}
+		
+		if (hwe_level!=1) {
+			if (CalcHWE(n_0, n_2, n_1)<hwe_level) {indicator_snp.push_back(0); continue;}
+		}
+			
+		
+		//filter SNP if it is correlated with W
+		for (size_t i=0; i<genotype->size; ++i) {			
+			if (gsl_vector_get (genotype_miss, i)==1) {geno=maf*2.0; gsl_vector_set (genotype, i, geno);}		
+		}
+		
+		gsl_blas_dgemv (CblasTrans, 1.0, W, genotype, 0.0, Wtx);
+		gsl_blas_dgemv (CblasNoTrans, 1.0, WtWi, Wtx, 0.0, WtWiWtx);
+		gsl_blas_ddot (genotype, genotype, &v_x);
+		gsl_blas_ddot (Wtx, WtWiWtx, &v_w);
+		
+		if (v_w/v_x > r2_level) {indicator_snp.push_back(0); continue;}
+		
+		indicator_snp.push_back(1); 
+		ns_test++;
+	}
+	
+	gsl_vector_free (genotype);
+	gsl_vector_free (genotype_miss);
+	gsl_matrix_free (WtW);
+	gsl_matrix_free (WtWi);
+	gsl_vector_free (Wtx);
+	gsl_vector_free (WtWiWtx);
+	gsl_permutation_free (pmt);
+		  
+	infile.close();
+	infile.clear();	
+	
+	return true;
+}
+
+
+
+void ReadFile_kin (const string &file_kin, vector<int> &indicator_idv, map<string, int> &mapID2num, const size_t k_mode, bool &error, gsl_matrix *G) 
+{
+	igzstream infile (file_kin.c_str(), igzstream::in);
+//	ifstream infile (file_kin.c_str(), ifstream::in);
+	if (!infile) {cout<<"error! fail to open kinship file: "<<file_kin<<endl; error=true; return;}
+	
+	size_t ni_total=indicator_idv.size();
+	
+	gsl_matrix_set_zero (G);
+	
+	string line;
+	char *ch_ptr;	
+	double d;
+	
+	if (k_mode==1) {
+		size_t i_test=0, i_total=0, j_test=0, j_total=0;
+		while (getline(infile, line)) {
+			if (i_total==ni_total) {cout<<"error! number of rows in the kinship file is larger than the number of phentypes."<<endl; error=true;}			
+			
+			if (indicator_idv[i_total]==0) {i_total++; continue;}
+			
+			j_total=0; j_test=0;
+			ch_ptr=strtok ((char *)line.c_str(), " , \t");
+			while (ch_ptr!=NULL) {
+				if (j_total==ni_total) {cout<<"error! number of columns in the kinship file is larger than the number of phentypes for row = "<<i_total<<endl; error=true;}
+				
+				d=atof(ch_ptr);
+				if (indicator_idv[j_total]==1) {gsl_matrix_set (G, i_test, j_test, d); j_test++;}				
+				j_total++;
+				
+				ch_ptr=strtok (NULL, " , \t");
+			}
+			if (j_total!=ni_total) {cout<<"error! number of columns in the kinship file do not match the number of phentypes for row = "<<i_total<<endl; error=true;}
+			i_total++; i_test++;			
+		}
+		if (i_total!=ni_total) {cout<<"error! number of rows in the kinship file do not match the number of phentypes."<<endl; error=true;}
+	}	
+	else {  
+		map<size_t, size_t> mapID2ID;
+		size_t c=0;
+		for (size_t i=0; i<indicator_idv.size(); i++) {
+			if (indicator_idv[i]==1) {mapID2ID[i]=c; c++;}
+		}
+		
+		string id1, id2;
+		double Cov_d;
+		size_t n_id1, n_id2;
+		
+		while (getline(infile, line)) {
+			ch_ptr=strtok ((char *)line.c_str(), " , \t");
+			id1=ch_ptr;
+			ch_ptr=strtok (NULL, " , \t");
+			id2=ch_ptr;
+			ch_ptr=strtok (NULL, " , \t");
+			d=atof(ch_ptr);
+			if (mapID2num.count(id1)==0 || mapID2num.count(id2)==0) {continue;}
+			if (indicator_idv[mapID2num[id1]]==0 || indicator_idv[mapID2num[id2]]==0) {continue;}
+			
+			n_id1=mapID2ID[mapID2num[id1]];
+			n_id2=mapID2ID[mapID2num[id2]];
+			
+			Cov_d=gsl_matrix_get(G, n_id1, n_id2);
+			if (Cov_d!=0 && Cov_d!=d) {cout<<"error! redundant and unequal terms in the kinship file, for id1 = "<<id1<<" and id2 = "<<id2<<endl;}
+			else {
+				gsl_matrix_set(G, n_id1, n_id2, d);
+				gsl_matrix_set(G, n_id2, n_id1, d);
+			}
+		}
+	}
+	
+	infile.close();
+	infile.clear();	
+	
+	return;
+}
+
+
+void ReadFile_mk (const string &file_mk, vector<int> &indicator_idv, map<string, int> &mapID2num, const size_t k_mode, bool &error, gsl_matrix *G) 
+{
+	igzstream infile (file_mk.c_str(), igzstream::in);
+	if (!infile) {cout<<"error! fail to open file: "<<file_mk<<endl; error=true; return;}
+
+	string file_kin, line;
+
+	size_t i=0;
+	while (getline(infile, line)) {
+	  file_kin=line.c_str();
+	  gsl_matrix_view G_sub=gsl_matrix_submatrix(G, 0, i*G->size1, G->size1, G->size1);
+	  ReadFile_kin (file_kin, indicator_idv, mapID2num, k_mode, error, &G_sub.matrix);
+	  i++;
+	}
+
+	infile.close();
+	infile.clear();	
+	return;
+}
+
+
+void ReadFile_eigenU (const string &file_ku, bool &error, gsl_matrix *U) 
+{
+	igzstream infile (file_ku.c_str(), igzstream::in);
+//	ifstream infile (file_ku.c_str(), ifstream::in);
+	if (!infile) {cout<<"error! fail to open the U file: "<<file_ku<<endl; error=true; return;}
+	
+	size_t n_row=U->size1, n_col=U->size2, i_row=0, i_col=0;
+	
+	gsl_matrix_set_zero (U);
+	
+	string line;
+	char *ch_ptr;	
+	double d;
+	
+	while (getline(infile, line)) {
+		if (i_row==n_row) {cout<<"error! number of rows in the U file is larger than expected."<<endl; error=true;}			
+				
+		i_col=0;
+		ch_ptr=strtok ((char *)line.c_str(), " , \t");
+		while (ch_ptr!=NULL) {
+			if (i_col==n_col) {cout<<"error! number of columns in the U file is larger than expected, for row = "<<i_row<<endl; error=true;}
+			
+			d=atof(ch_ptr);
+			gsl_matrix_set (U, i_row, i_col, d);			
+			i_col++;
+			
+			ch_ptr=strtok (NULL, " , \t");
+		}
+		
+		i_row++;
+	}
+		
+	infile.close();
+	infile.clear();	
+	
+	return;
+}
+
+
+
+
+void ReadFile_eigenD (const string &file_kd, bool &error, gsl_vector *eval) 
+{
+	igzstream infile (file_kd.c_str(), igzstream::in);
+//	ifstream infile (file_kd.c_str(), ifstream::in);
+	if (!infile) {cout<<"error! fail to open the D file: "<<file_kd<<endl; error=true; return;}
+	
+	size_t n_row=eval->size, i_row=0;
+	
+	gsl_vector_set_zero (eval);
+	
+	string line;
+	char *ch_ptr;	
+	double d;
+	
+	while (getline(infile, line)) {
+		if (i_row==n_row) {cout<<"error! number of rows in the D file is larger than expected."<<endl; error=true;}			
+		
+		ch_ptr=strtok ((char *)line.c_str(), " , \t");
+		d=atof(ch_ptr);
+		
+		ch_ptr=strtok (NULL, " , \t");
+		if (ch_ptr!=NULL) {cout<<"error! number of columns in the D file is larger than expected, for row = "<<i_row<<endl; error=true;}
+		
+		gsl_vector_set (eval, i_row, d);
+		
+		i_row++;
+	}
+	
+	infile.close();
+	infile.clear();	
+	
+	return;
+}
+
+
+
+//read bimbam mean genotype file and calculate kinship matrix
+bool BimbamKin (const string &file_geno, vector<int> &indicator_snp, const int k_mode, const int display_pace, gsl_matrix *matrix_kin) 
+{
+	igzstream infile (file_geno.c_str(), igzstream::in);
+	//ifstream infile (file_geno.c_str(), ifstream::in);
+	if (!infile) {cout<<"error reading genotype file:"<<file_geno<<endl; return false;}
+	
+	string line;
+	char *ch_ptr;
+	
+	size_t n_miss;
+	double d, geno_mean, geno_var;
+	
+	size_t ni_total=matrix_kin->size1;
+	gsl_vector *geno=gsl_vector_alloc (ni_total);
+	gsl_vector *geno_miss=gsl_vector_alloc (ni_total);
+
+	size_t ns_test=0;
+	for (size_t t=0; t<indicator_snp.size(); ++t) {
+		!safeGetline(infile, line).eof();
+		if (t%display_pace==0 || t==(indicator_snp.size()-1)) {ProgressBar ("Reading SNPs  ", t, indicator_snp.size()-1);}
+		if (indicator_snp[t]==0) {continue;}
+		
+		ch_ptr=strtok ((char *)line.c_str(), " , \t");
+		ch_ptr=strtok (NULL, " , \t");
+		ch_ptr=strtok (NULL, " , \t");
+		
+		geno_mean=0.0; n_miss=0; geno_var=0.0;
+		gsl_vector_set_all(geno_miss, 0);
+		for (size_t i=0; i<ni_total; ++i) {
+			ch_ptr=strtok (NULL, " , \t");
+			if (strcmp(ch_ptr, "NA")==0) {gsl_vector_set(geno_miss, i, 0); n_miss++;}
+			else {
+				d=atof(ch_ptr);
+				gsl_vector_set (geno, i, d);
+				gsl_vector_set (geno_miss, i, 1);
+				geno_mean+=d;
+				geno_var+=d*d;
+			}
+		}
+		
+		geno_mean/=(double)(ni_total-n_miss);
+		geno_var+=geno_mean*geno_mean*(double)n_miss;
+		geno_var/=(double)ni_total;
+		geno_var-=geno_mean*geno_mean;
+//		geno_var=geno_mean*(1-geno_mean*0.5);
+		
+		for (size_t i=0; i<ni_total; ++i) {
+			if (gsl_vector_get (geno_miss, i)==0) {gsl_vector_set(geno, i, geno_mean);}
+		}		
+		
+		gsl_vector_add_constant (geno, -1.0*geno_mean);
+		
+		if (geno_var!=0) {
+			if (k_mode==1) {gsl_blas_dsyr (CblasUpper, 1.0, geno, matrix_kin);}
+			else if (k_mode==2) {gsl_blas_dsyr (CblasUpper, 1.0/geno_var, geno, matrix_kin);}
+			else {cout<<"Unknown kinship mode."<<endl;}
+		}
+		
+		ns_test++;
+    }	
+	cout<<endl;
+	
+	gsl_matrix_scale (matrix_kin, 1.0/(double)ns_test);
+	
+	for (size_t i=0; i<ni_total; ++i) {
+		for (size_t j=0; j<i; ++j) {
+			d=gsl_matrix_get (matrix_kin, j, i);
+			gsl_matrix_set (matrix_kin, i, j, d);
+		}
+	}
+	
+	gsl_vector_free (geno);
+	gsl_vector_free (geno_miss);
+	
+	infile.close();
+	infile.clear();	
+	
+	return true;
+}
+
+
+
+
+
+
+
+bool PlinkKin (const string &file_bed, vector<int> &indicator_snp, const int k_mode, const int display_pace, gsl_matrix *matrix_kin) 
+{
+	ifstream infile (file_bed.c_str(), ios::binary);
+	if (!infile) {cout<<"error reading bed file:"<<file_bed<<endl; return false;}
+		
+	char ch[1];
+	bitset<8> b;
+	
+	size_t n_miss, ci_total;
+	double d, geno_mean, geno_var;
+	
+	size_t ni_total=matrix_kin->size1;
+	gsl_vector *geno=gsl_vector_alloc (ni_total);
+
+	size_t ns_test=0;
+	int n_bit;
+	
+	//calculate n_bit and c, the number of bit for each snp
+	if (ni_total%4==0) {n_bit=ni_total/4;}
+	else {n_bit=ni_total/4+1; }
+
+	//print the first three majic numbers
+	for (int i=0; i<3; ++i) {
+		infile.read(ch,1);
+		b=ch[0];
+	}	
+	
+	for (size_t t=0; t<indicator_snp.size(); ++t) {
+		if (t%display_pace==0 || t==(indicator_snp.size()-1)) {ProgressBar ("Reading SNPs  ", t, indicator_snp.size()-1);}
+		if (indicator_snp[t]==0) {continue;}
+		
+		infile.seekg(t*n_bit+3);		//n_bit, and 3 is the number of magic numbers
+		
+		//read genotypes
+		geno_mean=0.0;	n_miss=0; ci_total=0; geno_var=0.0;
+		for (int i=0; i<n_bit; ++i) {
+			infile.read(ch,1);
+			b=ch[0];
+			for (size_t j=0; j<4; ++j) {                //minor allele homozygous: 2.0; major: 0.0;
+				if ((i==(n_bit-1)) && ci_total==ni_total) {break;}
+
+				if (b[2*j]==0) {
+					if (b[2*j+1]==0) {gsl_vector_set(geno, ci_total, 2.0); geno_mean+=2.0; geno_var+=4.0; }
+					else {gsl_vector_set(geno, ci_total, 1.0); geno_mean+=1.0; geno_var+=1.0;}
+				}
+				else {
+					if (b[2*j+1]==1) {gsl_vector_set(geno, ci_total, 0.0); }      
+					else {gsl_vector_set(geno, ci_total, -9.0); n_miss++; }
+				}
+
+				ci_total++;
+			}
+		}
+				
+		geno_mean/=(double)(ni_total-n_miss);
+		geno_var+=geno_mean*geno_mean*(double)n_miss;
+		geno_var/=(double)ni_total;
+		geno_var-=geno_mean*geno_mean;
+//		geno_var=geno_mean*(1-geno_mean*0.5);
+		
+		for (size_t i=0; i<ni_total; ++i) {
+			d=gsl_vector_get(geno,i);
+			if (d==-9.0) {gsl_vector_set(geno, i, geno_mean);}
+		}		
+		
+		gsl_vector_add_constant (geno, -1.0*geno_mean);
+		
+		if (geno_var!=0) {
+			if (k_mode==1) {gsl_blas_dsyr (CblasUpper, 1.0, geno, matrix_kin);}
+			else if (k_mode==2) {gsl_blas_dsyr (CblasUpper, 1.0/geno_var, geno, matrix_kin);}
+			else {cout<<"Unknown kinship mode."<<endl;}
+		}
+		
+		ns_test++;
+    }	
+	cout<<endl;
+	
+	gsl_matrix_scale (matrix_kin, 1.0/(double)ns_test);
+	
+	for (size_t i=0; i<ni_total; ++i) {
+		for (size_t j=0; j<i; ++j) {
+			d=gsl_matrix_get (matrix_kin, j, i);
+			gsl_matrix_set (matrix_kin, i, j, d);
+		}
+	}
+	
+	gsl_vector_free (geno);
+	
+	infile.close();
+	infile.clear();	
+	
+	return true;
+}
+
+
+
+
+
+//Read bimbam mean genotype file, the second time, recode "mean" genotype and calculate K
+bool ReadFile_geno (const string &file_geno, vector<int> &indicator_idv, vector<int> &indicator_snp, gsl_matrix *UtX, gsl_matrix *K, const bool calc_K)
+{
+	igzstream infile (file_geno.c_str(), igzstream::in);
+//	ifstream infile (file_geno.c_str(), ifstream::in);
+	if (!infile) {cout<<"error reading genotype file:"<<file_geno<<endl; return false;}
+	
+	string line;
+	char *ch_ptr;
+	
+	if (calc_K==true) {gsl_matrix_set_zero (K);}
+	
+	gsl_vector *genotype=gsl_vector_alloc (UtX->size1);
+	gsl_vector *genotype_miss=gsl_vector_alloc (UtX->size1);
+	double geno, geno_mean;
+	size_t n_miss;
+	
+	int ni_total=(int)indicator_idv.size();
+	int ns_total=(int)indicator_snp.size();
+	int ni_test=UtX->size1;
+	int ns_test=UtX->size2;
+	
+	int c_idv=0, c_snp=0;
+	
+	for (int i=0; i<ns_total; ++i) {
+		!safeGetline(infile, line).eof();
+		if (indicator_snp[i]==0) {continue;}	
+		
+		ch_ptr=strtok ((char *)line.c_str(), " , \t");
+		ch_ptr=strtok (NULL, " , \t");
+		ch_ptr=strtok (NULL, " , \t");
+		
+		c_idv=0; geno_mean=0; n_miss=0;
+		gsl_vector_set_zero (genotype_miss);
+		for (int j=0; j<ni_total; ++j) {
+			ch_ptr=strtok (NULL, " , \t");
+			if (indicator_idv[j]==0) {continue;}			
+			
+			if (strcmp(ch_ptr, "NA")==0) {gsl_vector_set (genotype_miss, c_idv, 1); n_miss++;}
+			else {			
+				geno=atof(ch_ptr);
+				gsl_vector_set (genotype, c_idv, geno); 
+				geno_mean+=geno;
+			}
+			c_idv++;
+		}
+		
+		geno_mean/=(double)(ni_test-n_miss);
+		
+		for (size_t i=0; i<genotype->size; ++i) {			
+			if (gsl_vector_get (genotype_miss, i)==1) {geno=0;}
+			else {geno=gsl_vector_get (genotype, i); geno-=geno_mean;}
+			
+			gsl_vector_set (genotype, i, geno);
+			gsl_matrix_set (UtX, i, c_snp, geno);
+		}
+		
+		if (calc_K==true) {gsl_blas_dsyr (CblasUpper, 1.0, genotype, K);}
+		
+		c_snp++;
+	}	
+	
+	if (calc_K==true) {
+		gsl_matrix_scale (K, 1.0/(double)ns_test);
+		
+		for (size_t i=0; i<genotype->size; ++i) {
+			for (size_t j=0; j<i; ++j) {
+				geno=gsl_matrix_get (K, j, i);
+				gsl_matrix_set (K, i, j, geno);
+			}
+		}
+	}
+	
+	gsl_vector_free (genotype);
+	gsl_vector_free (genotype_miss);
+	
+	infile.clear();
+	infile.close();
+	
+	return true;
+}
+
+
+
+
+
+
+
+//Read bimbam mean genotype file, the second time, recode "mean" genotype and calculate K
+bool ReadFile_bed (const string &file_bed, vector<int> &indicator_idv, vector<int> &indicator_snp, gsl_matrix *UtX, gsl_matrix *K, const bool calc_K)
+{
+	ifstream infile (file_bed.c_str(), ios::binary);
+	if (!infile) {cout<<"error reading bed file:"<<file_bed<<endl; return false;}
+	
+	char ch[1];
+	bitset<8> b;
+	
+	int ni_total=(int)indicator_idv.size();
+	int ns_total=(int)indicator_snp.size();
+	int ni_test=UtX->size1;
+	int ns_test=UtX->size2;
+	int n_bit;
+	
+	if (ni_total%4==0) {n_bit=ni_total/4;}
+	else {n_bit=ni_total/4+1;}
+	
+	//print the first three majic numbers
+	for (int i=0; i<3; ++i) {
+		infile.read(ch,1);
+		b=ch[0];
+	}
+	
+	if (calc_K==true) {gsl_matrix_set_zero (K);}
+	
+	gsl_vector *genotype=gsl_vector_alloc (UtX->size1);	
+	
+	double geno, geno_mean;
+	size_t n_miss;	
+	int c_idv=0, c_snp=0, c=0;
+	
+	//start reading snps and doing association test
+	for (int t=0; t<ns_total; ++t) {
+		if (indicator_snp[t]==0) {continue;}	
+		infile.seekg(t*n_bit+3);		//n_bit, and 3 is the number of magic numbers
+		
+		//read genotypes
+		c_idv=0; geno_mean=0.0; n_miss=0; c=0;
+		for (int i=0; i<n_bit; ++i) {
+			infile.read(ch,1);
+			b=ch[0];
+			for (size_t j=0; j<4; ++j) {                //minor allele homozygous: 2.0; major: 0.0;
+				if ((i==(n_bit-1)) && c==ni_total) {break;}				
+				if (indicator_idv[c]==0) {c++; continue;}
+				c++;
+				
+				if (b[2*j]==0) {
+					if (b[2*j+1]==0) {gsl_vector_set(genotype, c_idv, 2.0); geno_mean+=2.0;}
+					else {gsl_vector_set(genotype, c_idv, 1.0); geno_mean+=1.0;}
+				}
+				else {
+					if (b[2*j+1]==1) {gsl_vector_set(genotype, c_idv, 0.0); geno_mean+=0.0;}                               
+					else {gsl_vector_set(genotype, c_idv, -9.0); n_miss++;}
+				}
+				c_idv++;
+			}
+		}
+		
+		geno_mean/=(double)(ni_test-n_miss);
+		
+		for (size_t i=0; i<genotype->size; ++i) {		
+			geno=gsl_vector_get (genotype, i);
+			if (geno==-9) {geno=0;}
+			else {geno-=geno_mean;}
+			
+			gsl_vector_set (genotype, i, geno);
+			gsl_matrix_set (UtX, i, c_snp, geno);
+		}
+		
+		if (calc_K==true) {gsl_blas_dsyr (CblasUpper, 1.0, genotype, K);}
+		
+		c_snp++;
+	}	
+	
+	if (calc_K==true) {
+		gsl_matrix_scale (K, 1.0/(double)ns_test);
+		
+		for (size_t i=0; i<genotype->size; ++i) {
+			for (size_t j=0; j<i; ++j) {
+				geno=gsl_matrix_get (K, j, i);
+				gsl_matrix_set (K, i, j, geno);
+			}
+		}
+	}
+	
+	gsl_vector_free (genotype);		  
+	infile.clear();
+	infile.close();
+	
+	return true;
+}
+
+
+
+
+
+bool ReadFile_est (const string &file_est, const vector<size_t> &est_column, map<string, double> &mapRS2est)
+{
+	mapRS2est.clear();
+	
+	ifstream infile (file_est.c_str(), ifstream::in);
+	if (!infile) {cout<<"error opening estimated parameter file: "<<file_est<<endl; return false;}
+	
+	string line;
+	char *ch_ptr;
+	
+	string rs;
+	double alpha, beta, gamma, d;
+	
+	//header
+	getline(infile, line);
+	
+	size_t n=*max_element(est_column.begin(), est_column.end());
+	
+	while (getline(infile, line)) {
+		ch_ptr=strtok ((char *)line.c_str(), " \t");		
+		
+		alpha=0.0; beta=0.0; gamma=1.0;
+		for (size_t i=0; i<n+1; ++i) {
+			if (i==est_column[0]-1) {rs=ch_ptr;}
+			if (i==est_column[1]-1) {alpha=atof(ch_ptr);}
+			if (i==est_column[2]-1) {beta=atof(ch_ptr);}
+			if (i==est_column[3]-1) {gamma=atof(ch_ptr);}
+			if (i<n) {ch_ptr=strtok (NULL, " \t");}
+		}
+		
+		d=alpha+beta*gamma;
+		
+		if (mapRS2est.count(rs)==0) {
+			mapRS2est[rs]=d;
+		}
+		else {
+			cout<<"the same SNP occurs more than once in estimated parameter file: "<<rs<<endl; return false;
+		}
+	}
+	
+	infile.clear();
+	infile.close();
+	return true;
+}
+
+
+
+bool CountFileLines (const string &file_input, size_t &n_lines)
+{
+	igzstream infile (file_input.c_str(), igzstream::in);
+	//ifstream infile (file_input.c_str(), ifstream::in);
+	if (!infile) {cout<<"error! fail to open file: "<<file_input<<endl; return false;}
+
+	n_lines=count(istreambuf_iterator<char>(infile), istreambuf_iterator<char>(), '\n');
+	infile.seekg (0, ios::beg);
+	
+	return true;
+}
+
+
+
+//Read gene expression file
+bool ReadFile_gene (const string &file_gene, vector<double> &vec_read, vector<SNPINFO> &snpInfo, size_t &ng_total)
+{
+	vec_read.clear();
+	ng_total=0;
+	
+	ifstream infile (file_gene.c_str(), ifstream::in);
+	if (!infile) {cout<<"error! fail to open gene expression file: "<<file_gene<<endl; return false;}
+	
+	string line;
+	char *ch_ptr;
+	string rs;
+	
+	size_t n_idv=0, t=0;
+	
+	//header
+	getline(infile, line);
+	
+	while (getline(infile, line)) {
+		ch_ptr=strtok ((char *)line.c_str(), " , \t");
+		rs=ch_ptr;
+		
+		ch_ptr=strtok (NULL, " , \t");	
+		
+		t=0;
+		while (ch_ptr!=NULL) {
+			if (ng_total==0) {
+				vec_read.push_back(0);
+				t++;
+				n_idv++;
+			} else {
+				vec_read[t]+=atof(ch_ptr);		
+				t++;
+			}
+			
+			ch_ptr=strtok (NULL, " , \t");	
+		}
+		
+		if (t!=n_idv) {cout<<"error! number of columns doesn't match in row: "<<ng_total<<endl; return false;}
+		
+		SNPINFO sInfo={"-9", rs, -9, -9, "-9", "-9", -9, -9, -9};
+		snpInfo.push_back(sInfo);
+		
+		ng_total++;
+	}
+	
+	infile.close();
+	infile.clear();	
+	
+	return true;
+}
+
+
-- 
cgit v1.2.3