From b3b491cd9143d33bfebd4c5b26629573afcf0970 Mon Sep 17 00:00:00 2001 From: xiangzhou Date: Sat, 11 Jul 2015 12:57:37 -0400 Subject: add GXE test --- src/io.cpp | 2482 +++++++++++++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 2131 insertions(+), 351 deletions(-) (limited to 'src/io.cpp') diff --git a/src/io.cpp b/src/io.cpp index 7ed95c4..d6a70dd 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -28,7 +28,7 @@ #include #include #include -#include +#include #include "gsl/gsl_vector.h" #include "gsl/gsl_matrix.h" @@ -39,6 +39,7 @@ #include "lapack.h" #include "gzstream.h" #include "mathfunc.h" +#include "eigenlib.h" #ifdef FORCE_FLOAT #include "io_float.h" @@ -54,10 +55,10 @@ 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); + double progress = (100.0 * p / total); + int barsize = (int) (progress / 2.0); char bar[51]; - + cout< &setSnps) ifstream infile (file_snps.c_str(), ifstream::in); if (!infile) {cout<<"error! fail to open snps file: "< &mapRS2chr, map { mapRS2chr.clear(); mapRS2bp.clear(); - + ifstream infile (file_anno.c_str(), ifstream::in); if (!infile) {cout<<"error opening annotation file: "< &mapRS2chr, map 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(); - + infile.clear(); + return true; } @@ -225,28 +226,28 @@ bool ReadFile_column (const string &file_pheno, vector &indicator_idv, vect { 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: "< > &indicator_p { 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: "< pheno_row; vector ind_pheno_row; - + size_t p_max=*max_element(p_column.begin(), p_column.end() ); map mapP2c; for (size_t i=0; i > &indicator_p bool ReadFile_cvt (const string &file_cvt, vector &indicator_cvt, vector > &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: "< 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"); + ch_ptr=strtok (NULL, " , \t"); } - if (flag_na==0) {indicator_cvt.push_back(1);} else {indicator_cvt.push_back(0);} + 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::size_type i=0; i &indicator_cvt, vector &snpInfo) { snpInfo.clear(); - + ifstream infile (file_bim.c_str(), ifstream::in); if (!infile) {cout<<"error opening .bim file: "< &snpInfo) minor=ch_ptr; ch_ptr=strtok (NULL, " \t"); major=ch_ptr; - - SNPINFO sInfo={chr, rs, cM, b_pos, minor, major, -9, -9, -9}; + + SNPINFO sInfo={chr, rs, cM, b_pos, minor, major, 0, -9, -9, 0, 0, 0}; snpInfo.push_back(sInfo); } - + infile.close(); - infile.clear(); + infile.clear(); return true; } @@ -396,8 +397,8 @@ bool ReadFile_fam (const string &file_fam, vector > &indicator_pheno { indicator_pheno.clear(); pheno.clear(); - mapID2num.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: "< > &indicator_pheno vector pheno_row; vector ind_pheno_row; - + size_t p_max=*max_element(p_column.begin(), p_column.end() ); map mapP2c; for (size_t i=0; i > &indicator_pheno ch_ptr=strtok (NULL, " \t"); ch_ptr=strtok (NULL, " \t"); ch_ptr=strtok (NULL, " \t"); - + size_t i=0; while (i > &indicator_pheno 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"); + ch_ptr=strtok (NULL, " , \t"); } - + indicator_pheno.push_back(ind_pheno_row); - pheno.push_back(pheno_row); - + pheno.push_back(pheno_row); + mapID2num[id]=c; c++; } - + infile.close(); - infile.clear(); + infile.clear(); return true; } @@ -466,7 +467,7 @@ bool ReadFile_geno (const string &file_geno, const set &setSnps, const g { 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:"< &setSnps, const g 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; + //eigenlib_dgemm("T", "N", 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; - + size_t file_pos; + 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=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); - + + 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}; + 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, ni_test-n_miss, 0, file_pos}; snpInfo.push_back(sInfo); - + file_pos++; + if ( (double)n_miss/(double)ni_test > miss_level) {indicator_snp.push_back(0); continue;} - + if ( (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 && maf_level!=-1) { if (CalcHWE(n_0, n_2, n_1)size; ++i) { - if (gsl_vector_get (genotype_miss, i)==1) {geno=maf*2.0; gsl_vector_set (genotype, i, geno);} + for (size_t i=0; isize; ++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 (W->size2!=1 && v_w/v_x >= r2_level) {indicator_snp.push_back(0); continue;} - - indicator_snp.push_back(1); + + indicator_snp.push_back(1); ns_test++; } - + gsl_vector_free (genotype); gsl_vector_free (genotype_miss); gsl_matrix_free (WtW); @@ -591,10 +598,10 @@ bool ReadFile_geno (const string &file_geno, const set &setSnps, const g gsl_vector_free (Wtx); gsl_vector_free (WtWiWtx); gsl_permutation_free (pmt); - + infile.close(); - infile.clear(); - + infile.clear(); + return true; } @@ -602,13 +609,13 @@ bool ReadFile_geno (const string &file_geno, const set &setSnps, const g - + //Read bed file, the first time bool ReadFile_bed (const string &file_bed, const set &setSnps, const gsl_matrix *W, vector &indicator_idv, vector &indicator_snp, vector &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:"< &setSnps, const gsl 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; + 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 &setSnps, const gsl infile.read(ch,1); b=ch[0]; } - + double maf; size_t n_miss; - size_t n_0, n_1, n_2, c; - + size_t n_0, n_1, n_2, c; + //start reading snps and doing association test for (size_t t=0; t &setSnps, const gsl 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++;} + 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; - + snpInfo[t].n_idv=ni_test-n_miss; + snpInfo[t].n_nb=0; + snpInfo[t].file_position=t; + if ( (double)n_miss/(double)ni_test > miss_level) {indicator_snp.push_back(0); continue;} - + if ( (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 && maf_level!=-1) { if (CalcHWE(n_0, n_2, n_1)size; ++i) { - if (gsl_vector_get (genotype_miss, i)==1) {geno=maf*2.0; gsl_vector_set (genotype, i, geno);} + for (size_t i=0; isize; ++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 (W->size2!=1 && v_w/v_x > r2_level) {indicator_snp.push_back(0); continue;} - - indicator_snp.push_back(1); + + indicator_snp.push_back(1); ns_test++; } - + gsl_vector_free (genotype); gsl_vector_free (genotype_miss); gsl_matrix_free (WtW); @@ -728,63 +739,177 @@ bool ReadFile_bed (const string &file_bed, const set &setSnps, const gsl gsl_vector_free (Wtx); gsl_vector_free (WtWiWtx); gsl_permutation_free (pmt); - + infile.close(); - infile.clear(); - + infile.clear(); + return true; } -void ReadFile_kin (const string &file_kin, vector &indicator_idv, map &mapID2num, const size_t k_mode, bool &error, gsl_matrix *G) + + +//read the genotype for one SNP; remember to read empty lines +//geno stores original genotypes without centering +//missing values are replaced by mean +bool Bimbam_ReadOneSNP (const size_t inc, const vector &indicator_idv, igzstream &infile, gsl_vector *geno, double &geno_mean) +{ + size_t ni_total=indicator_idv.size(); + + // if (infile.eof()) {infile.clear();} + // infile.seekg(pos); + + string line; + char *ch_ptr; + bool flag=false; + + for (size_t i=0; i geno_miss; + + for (size_t i=0; i &indicator_idv, ifstream &infile, gsl_vector *geno, double &geno_mean) +{ + size_t ni_total=indicator_idv.size(), n_bit; + if (ni_total%4==0) {n_bit=ni_total/4;} + else {n_bit=ni_total/4+1;} + infile.seekg(pos*n_bit+3); //n_bit, and 3 is the number of magic numbers + + //read genotypes + char ch[1]; + bitset<8> b; + + geno_mean=0.0; + size_t c=0, c_idv=0; + vector geno_miss; + + for (size_t i=0; i &indicator_idv, map &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: "< mapID2ID; size_t c=0; for (size_t i=0; i &indicator_idv, map &indicator_idv, mapsize1, n_col=U->size2, i_row=0, i_col=0; - + gsl_matrix_set_zero (U); - + string line; - char *ch_ptr; + 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."<size, i_row=0; - + gsl_vector_set_zero (eval); - + string line; - char *ch_ptr; + 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."<size1; gsl_vector *geno=gsl_vector_alloc (ni_total); gsl_vector *geno_miss=gsl_vector_alloc (ni_total); @@ -934,11 +1059,11 @@ bool BimbamKin (const string &file_geno, vector &indicator_snp, const int k !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 &indicator_snp, const int k 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 &indicator_snp, const int k -bool PlinkKin (const string &file_bed, vector &indicator_snp, const int k_mode, const int display_pace, gsl_matrix *matrix_kin) +bool PlinkKin (const string &file_bed, vector &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:"< 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; } @@ -1024,14 +1154,14 @@ bool PlinkKin (const string &file_bed, vector &indicator_snp, const int k_m for (int i=0; i<3; ++i) { infile.read(ch,1); b=ch[0]; - } - + } + for (size_t t=0; t &indicator_snp, const int k_m 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); } + 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 &indicator_idv, vector< igzstream infile (file_geno.c_str(), igzstream::in); // ifstream infile (file_geno.c_str(), ifstream::in); if (!infile) {cout<<"error reading genotype file:"<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; isize; ++i) { + + for (size_t i=0; isize; ++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; isize; ++i) { for (size_t j=0; j &indicator_idv, vector< } } } - + gsl_vector_free (genotype); gsl_vector_free (genotype_miss); - + infile.clear(); infile.close(); - + return true; } +//compact version of the above function, using uchar instead of gsl_matrix +bool ReadFile_geno (const string &file_geno, vector &indicator_idv, vector &indicator_snp, vector > &Xt, gsl_matrix *K, const bool calc_K, const size_t ni_test, const size_t ns_test) +{ + igzstream infile (file_geno.c_str(), igzstream::in); + // ifstream infile (file_geno.c_str(), ifstream::in); + if (!infile) {cout<<"error reading genotype file:"< Xt_row; + for (size_t i=0; isize; ++j) { + if (gsl_vector_get (genotype_miss, j)==1) { + geno=geno_mean; + } else { + geno=gsl_vector_get (genotype, j); + } + + Xt_row[j]=Double02ToUchar(geno); + gsl_vector_set (genotype, j, (geno-geno_mean)); + } + Xt.push_back(Xt_row); + + 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; isize; ++i) { + for (size_t j=0; j &indicator_idv, vector 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; + + size_t ni_total=indicator_idv.size(); + size_t ns_total=indicator_snp.size(); + size_t ni_test=UtX->size1; + size_t 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); - + + 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; - + size_t n_miss; + size_t c_idv=0, c_snp=0, c=0; + //start reading snps and doing association test - for (int t=0; tsize; ++i) { + + for (size_t i=0; isize; ++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; isize; ++i) { + for (size_t j=0; j &indicator_idv, vector &indicator_snp, vector > &Xt, gsl_matrix *K, const bool calc_K, const size_t ni_test, const size_t ns_test) +{ + ifstream infile (file_bed.c_str(), ios::binary); + if (!infile) {cout<<"error reading bed file:"< Xt_row; + for (size_t i=0; i b; + + size_t ni_total=indicator_idv.size(); + size_t ns_total=indicator_snp.size(); + 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 (ni_test); + + double geno, geno_mean; + size_t n_miss; + size_t c_idv=0, c_snp=0, c=0; + + //start reading snps and doing association test + for (size_t t=0; tsize; ++i) { + geno=gsl_vector_get (genotype, i); + if (geno==-9) {geno=geno_mean;} + + Xt_row[i]=Double02ToUchar(geno); + + geno-=geno_mean; + + gsl_vector_set (genotype, i, geno); + } + Xt.push_back(Xt_row); + 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; isize; ++i) { for (size_t j=0; j &indicator_idv, vector &indicator_idv, vector &est_column, map &mapRS2est) { mapRS2est.clear(); - + ifstream infile (file_est.c_str(), ifstream::in); if (!infile) {cout<<"error opening estimated parameter file: "< &est_column, map if (i==est_column[3]-1) {gamma=atof(ch_ptr);} if (i &est_column, map cout<<"the same SNP occurs more than once in estimated parameter file: "<(infile), istreambuf_iterator(), '\n'); infile.seekg (0, ios::beg); - + return true; } @@ -1348,25 +1671,25 @@ bool ReadFile_gene (const string &file_gene, vector &vec_read, vector &vec_read, vector > &indicator_pheno, vector > &pheno, const vector &p_column, vector &indicator_cvt, vector > &cvt, size_t &n_cvt) +{ + indicator_pheno.clear(); + pheno.clear(); + indicator_cvt.clear(); + + igzstream infile (file_sample.c_str(), igzstream::in); + + if (!infile) {cout<<"error! fail to open sample file: "< pheno_row; + vector ind_pheno_row; + int flag_na=0; + + size_t num_cols=0; + size_t num_p_in_file=0; + size_t num_cvt_in_file=0; + +// size_t p_max=*max_element(p_column.begin(), p_column.end()); + + map mapP2c; + for (size_t i=0; i > cvt_factor_levels; + + char col_type[num_cols]; + // read header line2 + if(!safeGetline(infile, line).eof()) { + ch_ptr=strtok ((char *)line.c_str(), " "); + if(strcmp(ch_ptr, "0")!=0) {return false;} + ch_ptr=strtok(NULL, " "); + if(strcmp(ch_ptr, "0")!=0) {return false;} + ch_ptr=strtok(NULL, " "); + if(strcmp(ch_ptr, "0")!=0) {return false;} + size_t it=0; + ch_ptr=strtok (NULL, " "); + if(ch_ptr!=NULL) + while(ch_ptr!=NULL){ + col_type[it++]=ch_ptr[0]; + if(ch_ptr[0]=='D') {cvt_factor_levels.push_back(map());num_cvt_in_file++;} + if(ch_ptr[0]=='C') {num_cvt_in_file++;} + if((ch_ptr[0]=='P')||(ch_ptr[0]=='B')) {num_p_in_file++;} + ch_ptr=strtok(NULL, " "); + } + + } + + while (!safeGetline(infile, line).eof()) { + + ch_ptr=strtok ((char *)line.c_str(), " "); + + for(int it=0;it<3;it++){ch_ptr=strtok(NULL, " ");} + + + size_t i=0; + size_t p_i=0; + size_t fac_cvt_i=0; + + while (i0) + { + igzstream infile2 (file_sample.c_str(), igzstream::in); + + if (!infile2) {cout<<"error! fail to open sample file: "< v_d; flag_na=0; + ch_ptr=strtok ((char *)line.c_str(), " "); + + for(int it=0;it<3;it++){ch_ptr=strtok(NULL, " ");} + + + size_t i=0; + size_t fac_cvt_i=0; + size_t num_fac_levels; + while (i1) + { + if (strcmp(ch_ptr, "NA")==0) {flag_na=1; for(size_t it=0;it::size_type i=0; i +#include +bool ReadFile_bgen(const string &file_bgen, const set &setSnps, const gsl_matrix *W, vector &indicator_idv, vector &indicator_snp, vector &snpInfo, const double &maf_level, const double &miss_level, const double &hwe_level, const double &r2_level, size_t &ns_test) +{ + + indicator_snp.clear(); + + ifstream infile (file_bgen.c_str(), ios::binary); + if (!infile) {cout<<"error reading bgen file:"<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); + + // read in header + uint32_t bgen_snp_block_offset; + uint32_t bgen_header_length; + uint32_t bgen_nsamples; + uint32_t bgen_nsnps; + uint32_t bgen_flags; + infile.read(reinterpret_cast(&bgen_snp_block_offset),4); + infile.read(reinterpret_cast(&bgen_header_length),4); + bgen_snp_block_offset-=4; + infile.read(reinterpret_cast(&bgen_nsnps),4); + bgen_snp_block_offset-=4; + infile.read(reinterpret_cast(&bgen_nsamples),4); + bgen_snp_block_offset-=4; + infile.ignore(4+bgen_header_length-20); + bgen_snp_block_offset-=4+bgen_header_length-20; + infile.read(reinterpret_cast(&bgen_flags),4); + bgen_snp_block_offset-=4; + bool CompressedSNPBlocks=bgen_flags&0x1; + bool LongIds=bgen_flags&0x4; + + if(!LongIds) {return false;} + + infile.ignore(bgen_snp_block_offset); + + ns_test=0; + + size_t ns_total=static_cast(bgen_nsnps); + + snpInfo.clear(); + string rs; + long int b_pos; + string chr; +// double cM; + string major; + string minor; + string id; + + double v_x, v_w; + int c_idv=0; + + + double maf, geno, geno_old; + size_t n_miss; + size_t n_0, n_1, n_2; + int flag_poly; + + double bgen_geno_prob_AA, bgen_geno_prob_AB, bgen_geno_prob_BB, bgen_geno_prob_non_miss; + + + size_t ni_total=indicator_idv.size(); // total number of samples in phenotype file + size_t ni_test=0; // number of samples to use in test + + uint32_t bgen_N; + uint16_t bgen_LS; + uint16_t bgen_LR; + uint16_t bgen_LC; + uint32_t bgen_SNP_pos; + uint32_t bgen_LA; + std::string bgen_A_allele; + uint32_t bgen_LB; + std::string bgen_B_allele; + uint32_t bgen_P; + size_t unzipped_data_size; + + for (size_t i=0; i(&bgen_N),4); + infile.read(reinterpret_cast(&bgen_LS),2); + + id.resize(bgen_LS); + infile.read(&id[0], bgen_LS); + + infile.read(reinterpret_cast(&bgen_LR),2); + rs.resize(bgen_LR); + infile.read(&rs[0], bgen_LR); + + infile.read(reinterpret_cast(&bgen_LC),2); + chr.resize(bgen_LC); + infile.read(&chr[0], bgen_LC); + + infile.read(reinterpret_cast(&bgen_SNP_pos),4); + + infile.read(reinterpret_cast(&bgen_LA),4); + bgen_A_allele.resize(bgen_LA); + infile.read(&bgen_A_allele[0], bgen_LA); + + + infile.read(reinterpret_cast(&bgen_LB),4); + bgen_B_allele.resize(bgen_LB); + infile.read(&bgen_B_allele[0], bgen_LB); + + + // should we switch according to MAF? + minor=bgen_B_allele; + major=bgen_A_allele; + b_pos=static_cast(bgen_SNP_pos); + + uint16_t unzipped_data[3*bgen_N]; + + 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); + if(CompressedSNPBlocks) + infile.read(reinterpret_cast(&bgen_P),4); + else + bgen_P=6*bgen_N; + + infile.ignore(static_cast(bgen_P)); + + continue; + } + + + if(CompressedSNPBlocks) + { + infile.read(reinterpret_cast(&bgen_P),4); + uint8_t zipped_data[bgen_P]; + + unzipped_data_size=6*bgen_N; + + infile.read(reinterpret_cast(zipped_data),bgen_P); + int result=uncompress(reinterpret_cast(unzipped_data), reinterpret_cast(&unzipped_data_size), reinterpret_cast(zipped_data), static_cast (bgen_P)); + assert(result == Z_OK); + + } + else + { + bgen_P=6*bgen_N; + infile.read(reinterpret_cast(unzipped_data),bgen_P); + + } + + + 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 (size_t i=0; i(unzipped_data[i*3])/32768.0; + bgen_geno_prob_AB=static_cast(unzipped_data[i*3+1])/32768.0; + bgen_geno_prob_BB=static_cast(unzipped_data[i*3+2])/32768.0; + bgen_geno_prob_non_miss=bgen_geno_prob_AA+bgen_geno_prob_AB+bgen_geno_prob_BB; + + //CHECK 0.1 OK + if (bgen_geno_prob_non_miss<0.9) {gsl_vector_set (genotype_miss, c_idv, 1); n_miss++; c_idv++; continue;} + + + bgen_geno_prob_AA/=bgen_geno_prob_non_miss; + bgen_geno_prob_AB/=bgen_geno_prob_non_miss; + bgen_geno_prob_BB/=bgen_geno_prob_non_miss; + + geno=2.0*bgen_geno_prob_BB+bgen_geno_prob_AB; + 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); + + // CHECK WHAT THIS DOES + 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*static_cast(ni_test-n_miss); + + SNPINFO sInfo={chr, rs, -9, 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 (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 && maf_level!=-1) { + if (CalcHWE(n_0, n_2, n_1)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 (W->size2!=1 && v_w/v_x >= r2_level) {indicator_snp.push_back(0); continue;} + + indicator_snp.push_back(1); + ns_test++; + + } + + + + + return true; + +} + + +//read oxford genotype file and calculate kinship matrix +bool bgenKin (const string &file_oxford, vector &indicator_snp, const int k_mode, const int display_pace, gsl_matrix *matrix_kin) +{ + string file_bgen=file_oxford; + ifstream infile (file_bgen.c_str(), ios::binary); + if (!infile) {cout<<"error reading bgen file:"<(&bgen_snp_block_offset),4); + infile.read(reinterpret_cast(&bgen_header_length),4); + bgen_snp_block_offset-=4; + infile.read(reinterpret_cast(&bgen_nsnps),4); + bgen_snp_block_offset-=4; + infile.read(reinterpret_cast(&bgen_nsamples),4); + bgen_snp_block_offset-=4; + infile.ignore(4+bgen_header_length-20); + bgen_snp_block_offset-=4+bgen_header_length-20; + infile.read(reinterpret_cast(&bgen_flags),4); + bgen_snp_block_offset-=4; + bool CompressedSNPBlocks=bgen_flags&0x1; +// bool LongIds=bgen_flags&0x4; + + infile.ignore(bgen_snp_block_offset); + + double bgen_geno_prob_AA, bgen_geno_prob_AB, bgen_geno_prob_BB, bgen_geno_prob_non_miss; + + uint32_t bgen_N; + uint16_t bgen_LS; + uint16_t bgen_LR; + uint16_t bgen_LC; + uint32_t bgen_SNP_pos; + uint32_t bgen_LA; + std::string bgen_A_allele; + uint32_t bgen_LB; + std::string bgen_B_allele; + uint32_t bgen_P; + size_t unzipped_data_size; + string id; + string rs; + string chr; + double genotype; + + + 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(&bgen_N),4); + infile.read(reinterpret_cast(&bgen_LS),2); + + id.resize(bgen_LS); + infile.read(&id[0], bgen_LS); + + infile.read(reinterpret_cast(&bgen_LR),2); + rs.resize(bgen_LR); + infile.read(&rs[0], bgen_LR); + + infile.read(reinterpret_cast(&bgen_LC),2); + chr.resize(bgen_LC); + infile.read(&chr[0], bgen_LC); + + infile.read(reinterpret_cast(&bgen_SNP_pos),4); + + infile.read(reinterpret_cast(&bgen_LA),4); + bgen_A_allele.resize(bgen_LA); + infile.read(&bgen_A_allele[0], bgen_LA); + + + infile.read(reinterpret_cast(&bgen_LB),4); + bgen_B_allele.resize(bgen_LB); + infile.read(&bgen_B_allele[0], bgen_LB); + + + + + uint16_t unzipped_data[3*bgen_N]; + + if (indicator_snp[t]==0) { + if(CompressedSNPBlocks) + infile.read(reinterpret_cast(&bgen_P),4); + else + bgen_P=6*bgen_N; + + infile.ignore(static_cast(bgen_P)); + + continue; + } + + + + if(CompressedSNPBlocks) + { + + + infile.read(reinterpret_cast(&bgen_P),4); + uint8_t zipped_data[bgen_P]; + + unzipped_data_size=6*bgen_N; + + infile.read(reinterpret_cast(zipped_data),bgen_P); + + int result=uncompress(reinterpret_cast(unzipped_data), reinterpret_cast(&unzipped_data_size), reinterpret_cast(zipped_data), static_cast (bgen_P)); + assert(result == Z_OK); + + } + else + { + + bgen_P=6*bgen_N; + infile.read(reinterpret_cast(unzipped_data),bgen_P); + } + + + + geno_mean=0.0; n_miss=0; geno_var=0.0; + gsl_vector_set_all(geno_miss, 0); + + for (size_t i=0; i(unzipped_data[i*3])/32768.0; + bgen_geno_prob_AB=static_cast(unzipped_data[i*3+1])/32768.0; + bgen_geno_prob_BB=static_cast(unzipped_data[i*3+2])/32768.0; + // WJA + bgen_geno_prob_non_miss=bgen_geno_prob_AA+bgen_geno_prob_AB+bgen_geno_prob_BB; + if (bgen_geno_prob_non_miss<0.9) {gsl_vector_set(geno_miss, i, 0.0); n_miss++;} + else { + + bgen_geno_prob_AA/=bgen_geno_prob_non_miss; + bgen_geno_prob_AB/=bgen_geno_prob_non_miss; + bgen_geno_prob_BB/=bgen_geno_prob_non_miss; + + genotype=2.0*bgen_geno_prob_BB+bgen_geno_prob_AB; + + gsl_vector_set(geno, i, genotype); + gsl_vector_set(geno_miss, i, 1.0); + geno_mean+=genotype; + geno_var+=genotype*genotype; + } + + } + + + 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 rs_set(rs_ptr, rs_ptr+10); + string chr_ptr[]={"chr","CHR"}; + set chr_set(chr_ptr, chr_ptr+2); + string pos_ptr[]={"ps","PS","pos","POS","base_position","BASE_POSITION", "bp", "BP"}; + set pos_set(pos_ptr, pos_ptr+8); + string cm_ptr[]={"cm","CM"}; + set cm_set(cm_ptr, cm_ptr+2); + string a1_ptr[]={"a1","A1","allele1","ALLELE1"}; + set a1_set(a1_ptr, a1_ptr+4); + string a0_ptr[]={"a0","A0","allele0","ALLELE0"}; + set a0_set(a0_ptr, a0_ptr+4); + + string z_ptr[]={"z","Z","z_score","Z_SCORE","zscore","ZSCORE"}; + set z_set(z_ptr, z_ptr+6); + string beta_ptr[]={"beta","BETA","b","B"}; + set beta_set(beta_ptr, beta_ptr+4); + string sebeta_ptr[]={"se_beta","SE_BETA","se","SE"}; + set sebeta_set(sebeta_ptr, sebeta_ptr+4); + string chisq_ptr[]={"chisq","CHISQ","chisquare","CHISQUARE"}; + set chisq_set(chisq_ptr, chisq_ptr+4); + string p_ptr[]={"p","P","pvalue","PVALUE","p-value","P-VALUE"}; + set p_set(p_ptr, p_ptr+6); + + string n_ptr[]={"n","N","ntotal","NTOTAL","n_total","N_TOTAL"}; + set n_set(n_ptr, n_ptr+6); + string nmis_ptr[]={"nmis","NMIS","n_mis","N_MIS","n_miss","N_MISS"}; + set nmis_set(nmis_ptr, nmis_ptr+6); + string nobs_ptr[]={"nobs","NOBS","n_obs","N_OBS"}; + set nobs_set(nobs_ptr, nobs_ptr+4); + + string af_ptr[]={"af","AF","maf","MAF","f","F","allele_freq","ALLELE_FREQ","allele_frequency","ALLELE_FREQUENCY"}; + set af_set(af_ptr, af_ptr+10); + string var_ptr[]={"var","VAR"}; + set var_set(var_ptr, var_ptr+2); + + string ws_ptr[]={"window_size","WINDOW_SIZE","ws","WS"}; + set ws_set(ws_ptr, ws_ptr+4); + string cor_ptr[]={"cor","COR","r","R"}; + set cor_set(cor_ptr, cor_ptr+4); + + header.rs_col=0; header.chr_col=0; header.pos_col=0; header.a1_col=0; header.a0_col=0; header.z_col=0; header.beta_col=0; header.sebeta_col=0; header.chisq_col=0; header.p_col=0; header.n_col=0; header.nmis_col=0; header.nobs_col=0; header.af_col=0; header.var_col=0; header.ws_col=0; header.cor_col=0; header.coln=0; + + char *ch_ptr; + string type; + size_t n_error=0; + + ch_ptr=strtok ((char *)line.c_str(), " , \t"); + while (ch_ptr!=NULL) { + type=ch_ptr; + if (rs_set.count(type)!=0) { + if (header.rs_col==0) {header.rs_col=header.coln+1;} else {cout<<"error! more than two rs columns in the file."< &mapRS2cat, size_t &n_vc) +{ + mapRS2cat.clear(); + + igzstream infile (file_cat.c_str(), igzstream::in); + if (!infile) {cout<<"error! fail to open category file: "<0) {n_vc++;} + + infile.clear(); + infile.close(); + + return true; +} + + + + +//read bimbam mean genotype file and calculate kinship matrix; this time, the kinship matrix is not centered, and can contain multiple K matrix +bool BimbamKin (const string &file_geno, vector &indicator_idv, vector &indicator_snp, const int k_mode, const int display_pace, const map &mapRS2cat, map &mapRS2var, vector &snpInfo, 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:"<size1; + gsl_vector *geno=gsl_vector_alloc (ni_test); + gsl_vector *geno_miss=gsl_vector_alloc (ni_test); + + size_t n_vc=matrix_kin->size2/ni_test, i_vc; + string rs; + vector ns_vec; + for (size_t i=0; i &indicator_idv, vector &indicator_snp, const int k_mode, const int display_pace, const map &mapRS2cat, map &mapRS2var, vector &snpInfo, gsl_matrix *matrix_kin) +{ + ifstream infile (file_bed.c_str(), ios::binary); + if (!infile) {cout<<"error reading bed file:"< b; + + size_t n_miss, ci_total, ci_test; + double d, geno_mean, geno_var; + + size_t ni_test=matrix_kin->size1; + size_t ni_total=indicator_idv.size(); + gsl_vector *geno=gsl_vector_alloc (ni_test); + + size_t ns_test=0; + int n_bit; + + size_t n_vc=matrix_kin->size2/ni_test, i_vc; + string rs; + vector ns_vec; + for (size_t i=0; i &mapRS2cat, const map &mapRS2var, gsl_vector *q, gsl_vector *s, size_t &ni_total, size_t &ns_total, size_t &ns_test) +{ + gsl_vector_set_zero(q); + ni_total=0; ns_total=0; ns_test=0; + + igzstream infile (file_beta.c_str(), igzstream::in); + if (!infile) {cout<<"error! fail to open beta file: "< vec_q, vec_s; + for (size_t i=0; isize; i++) { + vec_q.push_back(0.0); + vec_s.push_back(0.0); + } + + //read header + HEADER header; + !safeGetline(infile, line).eof(); + ReadHeader (line, header); + + if (header.n_col==0 ) { + if (header.nobs_col==0 && header.nmis_col==0) { + cout<<"error! missing sample size in the beta file."<size; i++) { + if (vec_s[i]!=0) { + gsl_vector_set(q, i, vec_q[i]/vec_s[i]); + } + gsl_vector_set(s, i, vec_s[i]); + } + + infile.clear(); + infile.close(); + + return; +} + + + + +//read S file: S and Svar +void ReadFile_s (const string &file_s, gsl_matrix *S, gsl_matrix *Svar) +{ + igzstream infile (file_s.c_str(), igzstream::in); + if (!infile) {cout<<"error! fail to open s file: "<size1; i++) { + !safeGetline(infile, line).eof(); + ch_ptr=strtok ((char *)line.c_str(), " , \t"); + for (size_t j=0; jsize2; j++) { + d=gsl_matrix_get(S, i, j)+atof(ch_ptr); + gsl_matrix_set(S, i, j, d); + ch_ptr=strtok (NULL, " , \t"); + } + } + + for (size_t i=0; isize1; i++) { + !safeGetline(infile, line).eof(); + ch_ptr=strtok ((char *)line.c_str(), " , \t"); + for (size_t j=0; jsize2; j++) { + d=gsl_matrix_get(Svar, i, j)+atof(ch_ptr); + gsl_matrix_set(Svar, i, j, d); + ch_ptr=strtok (NULL, " , \t"); + } + } + + infile.clear(); + infile.close(); + + return; +} + + + + +void ReadFile_ms (const string &file_ms, gsl_matrix *S, gsl_matrix *Svar) +{ + gsl_matrix_set_zero(S); + gsl_matrix_set_zero(Svar); + + string file_name; + + igzstream infile (file_ms.c_str(), igzstream::in); + if (!infile) {cout<<"error! fail to open ms file: "<size1; i++) { + !safeGetline(infile, line).eof(); + ch_ptr=strtok ((char *)line.c_str(), " , \t"); + for (size_t j=0; jsize2; j++) { + d=gsl_matrix_get(V, i, j)+atof(ch_ptr); + gsl_matrix_set(V, i, j, d); + ch_ptr=strtok (NULL, " , \t"); + } + } + + infile.clear(); + infile.close(); + + return; +} + + +void ReadFile_mv (const string &file_mv, gsl_matrix *V) +{ + gsl_matrix_set_zero(V); + + string file_name; + + igzstream infile (file_mv.c_str(), igzstream::in); + if (!infile) {cout<<"error! fail to open ms file: "<size; i++) { + !safeGetline(infile, line).eof(); + ch_ptr=strtok ((char *)line.c_str(), " , \t"); + d=gsl_vector_get(q_vec, i)+atof(ch_ptr); + gsl_vector_set(q_vec, i, d); + } + + for (size_t i=0; isize; i++) { + !safeGetline(infile, line).eof(); + ch_ptr=strtok ((char *)line.c_str(), " , \t"); + d=gsl_vector_get(s_vec, i)+atof(ch_ptr); + gsl_vector_set(s_vec, i, d); + } + + !safeGetline(infile, line).eof(); + ch_ptr=strtok ((char *)line.c_str(), " , \t"); + df=atof(ch_ptr); + + infile.clear(); + infile.close(); + + return; +} + + + +void ReadFile_mq (const string &file_mq, gsl_vector *q_vec, gsl_vector *s_vec, double &df) +{ + gsl_vector_set_zero(q_vec); + gsl_vector_set_zero(s_vec); + + string file_name; + + igzstream infile (file_mq.c_str(), igzstream::in); + if (!infile) {cout<<"error! fail to open mq file: "<