diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/bslmm.cpp | 30 | ||||
| -rw-r--r-- | src/io.cpp | 26 | ||||
| -rw-r--r-- | src/param.cpp | 1124 |
3 files changed, 699 insertions, 481 deletions
diff --git a/src/bslmm.cpp b/src/bslmm.cpp index 92762e2..563b743 100644 --- a/src/bslmm.cpp +++ b/src/bslmm.cpp @@ -1730,10 +1730,14 @@ void BSLMM::MCMC (const gsl_matrix *X, const gsl_vector *y) { gsl_vector_view Xtznew_sub=gsl_vector_subvector(Xtz_new, 0, rank_new.size()); gsl_vector_view betanew_sub=gsl_vector_subvector(beta_new, 0, rank_new.size()); - gsl_matrix_memcpy(&Xold_sub.matrix, &Xnew_sub.matrix); - gsl_matrix_memcpy(&XtXold_sub.matrix, &XtXnew_sub.matrix); - gsl_vector_memcpy(&Xtzold_sub.vector, &Xtznew_sub.vector); - gsl_vector_memcpy(&betaold_sub.vector, &betanew_sub.vector); + gsl_matrix_memcpy(&Xold_sub.matrix, + &Xnew_sub.matrix); + gsl_matrix_memcpy(&XtXold_sub.matrix, + &XtXnew_sub.matrix); + gsl_vector_memcpy(&Xtzold_sub.vector, + &Xtznew_sub.vector); + gsl_vector_memcpy(&betaold_sub.vector, + &betanew_sub.vector); } } else { cHyp_new=cHyp_old; @@ -1777,12 +1781,18 @@ void BSLMM::MCMC (const gsl_matrix *X, const gsl_vector *y) { } } - gsl_matrix_set(Result_hyp,w_col,0,cHyp_old.h); - gsl_matrix_set(Result_hyp,w_col,1,cHyp_old.pve); - gsl_matrix_set(Result_hyp,w_col,2,cHyp_old.rho); - gsl_matrix_set(Result_hyp,w_col,3,cHyp_old.pge); - gsl_matrix_set(Result_hyp,w_col,4,cHyp_old.logp); - gsl_matrix_set(Result_hyp,w_col,5,cHyp_old.n_gamma); + gsl_matrix_set(Result_hyp,w_col,0, + cHyp_old.h); + gsl_matrix_set(Result_hyp,w_col,1, + cHyp_old.pve); + gsl_matrix_set(Result_hyp,w_col,2, + cHyp_old.rho); + gsl_matrix_set(Result_hyp,w_col,3, + cHyp_old.pge); + gsl_matrix_set(Result_hyp,w_col,4, + cHyp_old.logp); + gsl_matrix_set(Result_hyp,w_col,5, + cHyp_old.n_gamma); for (size_t i=0; i<cHyp_old.n_gamma; ++i) { pos=mapRank2pos[rank_old[i]]+1; diff --git a/src/io.cpp b/src/io.cpp index e8f00ed..40f1c16 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -4089,30 +4089,8 @@ void ReadFile_mref (const string &file_mref, gsl_matrix *S_mat, gsl_matrix *Svar if (i!=j) {gsl_matrix_set(Svar_mat, j, i, d);} } } - /* - //final: update V - for (size_t i=0; i<n_vc; i++) { - d1=gsl_vector_get(s_vec, i); - if (d1==0) {continue;} - for (size_t j=i; j<n_vc; j++) { - d2=gsl_vector_get(s_vec, j); - if (d2==0) {continue;} - t_ij=GetabIndex (i+1, j+1, n_vc-2); - for (size_t l=0; l<n_vc+1; l++) { - if (l==n_vc) {d3=1;} else {d3=gsl_vector_get(s_vec, l);} - if (d3==0) {continue;} - for (size_t m=0; m<n_vc+1; m++) { - if (m==n_vc) {d4=1;} else {d4=gsl_vector_get(s_vec, m);} - if (d4==0) {continue;} - - d=gsl_matrix_get (V_mat, l, t_ij*(n_vc+1)+m)/(d1*d2*d3*d4); - gsl_matrix_set (V_mat, l, t_ij*(n_vc+1)+m, d); - } - } - } - } - */ - //free matrices + + // Free matrices. gsl_matrix_free(S_sub); gsl_matrix_free(Svar_sub); gsl_vector_free(s); diff --git a/src/param.cpp b/src/param.cpp index a7be178..a56f5eb 100644 --- a/src/param.cpp +++ b/src/param.cpp @@ -175,7 +175,6 @@ void PARAM::ReadFiles (void) { } } - // WJA added. // Read genotype and phenotype file for bgen format. if (!file_oxford.empty()) { @@ -281,7 +280,6 @@ void PARAM::ReadFiles (void) { ns_total=indicator_snp.size(); } - // Read genotype file for multiple PLINK files. if (!file_mbfile.empty()) { igzstream infile (file_mbfile.c_str(), igzstream::in); @@ -515,37 +513,94 @@ void PARAM::CheckParam (void) { 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; + 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 (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_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 (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 (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) { + 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 + // Check p_column, and (no need to) sort p_column into + // ascending order. if (p_column.size()==0) { p_column.push_back(1); } else { for (size_t i=0; i<p_column.size(); i++) { for (size_t j=0; j<i; j++) { - if (p_column[i]==p_column[j]) {cout<<"error! identical phenotype columns: "<<p_column[i]<<endl; error=true;} + if (p_column[i]==p_column[j]) { + cout<<"error! identical phenotype "<< + "columns: "<<p_column[i]<<endl; + error= + true;} } } } @@ -554,15 +609,22 @@ void PARAM::CheckParam (void) { // 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) { - cout<<"error! the current analysis mode "<<a_mode<<" can not deal with multiple phenotypes."<<endl; error=true; + if (n_ph>1 && a_mode!=1 && a_mode!=2 && a_mode!=3 && a_mode!=4 && + a_mode!=43) { + cout<<"error! the current analysis mode "<<a_mode<< + " can not deal with multiple phenotypes."<<endl; + error=true; } if (n_ph>1 && !file_gene.empty() ) { - cout<<"error! multiple phenotype analysis option not allowed with gene expression files. "<<endl; error=true; + 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; + cout<<"error! pnr value must be between 0 and 1. current value = "<< + p_nr<<endl; + error=true; } //check est_column @@ -580,86 +642,162 @@ void PARAM::CheckParam (void) { } } - 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;} + 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 + // Check if files are compatible with each other, and if files exist. if (!file_bfile.empty()) { str=file_bfile+".bim"; - if (stat(str.c_str(),&fileInfo)==-1) {cout<<"error! fail to open .bim file: "<<str<<endl; error=true;} + if (stat(str.c_str(),&fileInfo)==-1) { + cout<<"error! fail to open .bim file: "<<str<<endl; + error=true; + } str=file_bfile+".bed"; - if (stat(str.c_str(),&fileInfo)==-1) {cout<<"error! fail to open .bed file: "<<str<<endl; error=true;} + 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_oxford+".bgen"; - if (stat(str.c_str(),&fileInfo)==-1) {cout<<"error! fail to open .bgen file: "<<str<<endl; error=true;} + if (stat(str.c_str(),&fileInfo)==-1) { + cout<<"error! fail to open .bgen file: "<<str<<endl; + error=true; + } str=file_oxford+".sample"; - if (stat(str.c_str(),&fileInfo)==-1) {cout<<"error! fail to open .sample file: "<<str<<endl; error=true;} + 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;} + 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;} + 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;} + 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;} + if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) { + cout<<"error! fail to open category file: "<<str<<endl; + error=true; + } str=file_mcat; - if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open mcategory file: "<<str<<endl; error=true;} + if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) { + cout<<"error! fail to open mcategory 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;} + 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;} + if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) { + cout<<"error! fail to open correlation file: "<<str<<endl; + error=true; + } if (!file_study.empty()) { str=file_study+".Vq.txt"; - if (stat(str.c_str(),&fileInfo)==-1) {cout<<"error! fail to open .Vq.txt file: "<<str<<endl; error=true;} + if (stat(str.c_str(),&fileInfo)==-1) { + cout<<"error! fail to open .Vq.txt file: "<<str<<endl; + error=true; + } str=file_study+".q.txt"; - if (stat(str.c_str(),&fileInfo)==-1) {cout<<"error! fail to open .q.txt file: "<<str<<endl; error=true;} + if (stat(str.c_str(),&fileInfo)==-1) { + cout<<"error! fail to open .q.txt file: "<<str<<endl; + error=true; + } str=file_study+".size.txt"; - if (stat(str.c_str(),&fileInfo)==-1) {cout<<"error! fail to open .size.txt file: "<<str<<endl; error=true;} + if (stat(str.c_str(),&fileInfo)==-1) { + cout<<"error! fail to open .size.txt file: "<<str<<endl; + error=true; + } } if (!file_ref.empty()) { str=file_ref+".S.txt"; - if (stat(str.c_str(),&fileInfo)==-1) {cout<<"error! fail to open .S.txt file: "<<str<<endl; error=true;} + if (stat(str.c_str(),&fileInfo)==-1) { + cout<<"error! fail to open .S.txt file: "<<str<<endl; + error=true; + } str=file_ref+".size.txt"; - if (stat(str.c_str(),&fileInfo)==-1) {cout<<"error! fail to open .size.txt file: "<<str<<endl; error=true;} + if (stat(str.c_str(),&fileInfo)==-1) { + cout<<"error! fail to open .size.txt file: "<<str<<endl; + error=true; + } } str=file_mstudy; - if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open mstudy file: "<<str<<endl; error=true;} + if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) { + cout<<"error! fail to open mstudy file: "<<str<<endl; + error=true; + } str=file_mref; - if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open mref file: "<<str<<endl; error=true;} + if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) { + cout<<"error! fail to open mref file: "<<str<<endl; + error=true; + } str=file_mgeno; - if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open mgeno file: "<<str<<endl; error=true;} + if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) { + cout<<"error! fail to open mgeno file: "<<str<<endl; + error=true; + } str=file_mbfile; - if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open mbfile file: "<<str<<endl; error=true;} + if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) { + cout<<"error! fail to open mbfile file: "<<str<<endl; + error=true; + } size_t flag=0; if (!file_bfile.empty()) {flag++;} if (!file_geno.empty()) {flag++;} if (!file_gene.empty()) {flag++;} - // WJA added + + // WJA added. if (!file_oxford.empty()) {flag++;} - if (flag!=1 && a_mode!=15 && a_mode!=27 && a_mode!=28 && a_mode!=43 && a_mode!=5 && a_mode!=61 && a_mode!=62 && a_mode!=63 && a_mode!=66 && a_mode!=67) { - cout<<"error! either plink binary files, or bimbam mean genotype files, or gene expression files are required."<<endl; error=true; + if (flag!=1 && a_mode!=15 && a_mode!=27 && a_mode!=28 && + a_mode!=43 && a_mode!=5 && a_mode!=61 && a_mode!=62 && + a_mode!=63 && a_mode!=66 && a_mode!=67) { + 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) ) { @@ -668,28 +806,30 @@ void PARAM::CheckParam (void) { if (a_mode==61 || a_mode==62) { if (!file_beta.empty()) { - if ( file_mbfile.empty() && file_bfile.empty() && file_mgeno.empty() && file_geno.empty() && file_mref.empty() && file_ref.empty() ) { - cout<<"error! missing genotype file or ref/mref file."<<endl; error=true; + if ( file_mbfile.empty() && file_bfile.empty() && + file_mgeno.empty() && file_geno.empty() && + file_mref.empty() && file_ref.empty() ) { + cout<<"error! missing genotype file or ref/mref file."<<endl; + error=true; } } else if (!file_pheno.empty()) { - if (file_kin.empty() && (file_ku.empty()||file_kd.empty()) && file_mk.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_mstudy.empty() && file_study.empty()) || (file_mref.empty() && file_ref.empty() ) ) { - cout<<"error! either beta file, or phenotype files or study/ref mstudy/mref files are required."<<endl; error=true; + } else if ( (file_mstudy.empty() && file_study.empty()) || + (file_mref.empty() && file_ref.empty() ) ) { + cout<<"error! either beta file, or phenotype files or "<< + "study/ref mstudy/mref files are required."<<endl; + error=true; } } if (a_mode==63) { - if (file_kin.empty() && (file_ku.empty()||file_kd.empty()) && file_mk.empty() ) { - cout<<"error! missing relatedness file. "<<endl; error=true; + if (file_kin.empty() && (file_ku.empty()||file_kd.empty()) && + file_mk.empty() ) { + cout<<"error! missing relatedness file. "<<endl; error=true; } if ( file_pheno.empty() ) { cout<<"error! missing phenotype file."<<endl; error=true; @@ -697,111 +837,197 @@ void PARAM::CheckParam (void) { } if (a_mode==66 || a_mode==67) { - if (file_beta.empty() || ( file_mbfile.empty() && file_bfile.empty() && file_mgeno.empty() && file_geno.empty()) ) { - cout<<"error! missing beta file or genotype file."<<endl; error=true; + if (file_beta.empty() || + (file_mbfile.empty() && file_bfile.empty() && + file_mgeno.empty() && file_geno.empty()) ) { + cout<<"error! missing beta file or genotype file."<<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_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;} + 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;} + 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;} + 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;} + if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) { + cout<<"error! fail to open annotation file: "<<str<<endl; + error=true; + } str=file_kin; - if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open relatedness matrix file: "<<str<<endl; error=true;} + if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) { + cout<<"error! fail to open relatedness matrix file: "<<str<<endl; + error=true; + } str=file_mk; - if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open relatedness matrix file: "<<str<<endl; error=true;} + 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;} + 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;} + 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;} + 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;} + 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;} + 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;} + 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;} + // 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==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==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 || a_mode==14 || a_mode==16) && !file_cvt.empty() ) {cout<<"error! -bslmm option does not support covariates files."<<endl; error=true;} + if ((a_mode==11 || a_mode==12 || a_mode==13 || a_mode==14 || + a_mode==16) && !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; + cout<<"error! beta effects file is required."<<endl; + error=true; } return; } - - - - void PARAM::CheckData (void) { - if(file_oxford.empty()) // WJA NOTE: I added this condition so that covariates can be added through sample, probably not exactly what is wanted + + // WJA NOTE: I added this condition so that covariates can be added + // through sample, probably not exactly what is wanted. + if(file_oxford.empty()) { - if ((file_cvt).empty() || (indicator_cvt).size()==0) { - n_cvt=1; - } + if ((file_cvt).empty() || (indicator_cvt).size()==0) { + n_cvt=1; + } } if ( (a_mode==66 || a_mode==67) && (v_pve.size()!=n_vc)) { - cout<<"error! the number of pve estimates does not equal to the number of categories in the cat file:"<<v_pve.size()<<" "<<n_vc<<endl; error=true; - } - - 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; - return; - } + cout<<"error! the number of pve estimates does not equal to "<< + "the number of categories in the cat file:"<<v_pve.size()<<" "<< + n_vc<<endl; + error=true; + } - //calculate ni_total and ni_test, and set indicator_idv to 0 whenever indicator_cvt=0 - //and calculate np_obs and np_miss + 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; + return; + } + + // 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; @@ -839,24 +1065,8 @@ void PARAM::CheckData (void) { } } - /* - if ((indicator_cvt).size()!=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; - 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 && file_cor.empty() && file_mstudy.empty() && file_study.empty() && file_beta.empty() && file_bf.empty() ) { + if (ni_test==0 && file_cor.empty() && file_mstudy.empty() && + file_study.empty() && file_beta.empty() && file_bf.empty() ) { error=true; cout<<"error! number of analyzed individuals equals 0. "<<endl; return; @@ -865,23 +1075,27 @@ void PARAM::CheckData (void) { 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 - if (file_cor.empty() && file_mstudy.empty() && file_study.empty() && a_mode!=15 && a_mode!=27 && a_mode!=28) { + // Output some information. + if (file_cor.empty() && file_mstudy.empty() && file_study.empty() && + a_mode!=15 && a_mode!=27 && a_mode!=28) { 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; + cout<<"## number of individuals with full phenotypes = "<< + ni_test<<endl; } else { cout<<"## number of analyzed individuals = "<<ni_test<<endl; } @@ -899,12 +1113,12 @@ void PARAM::CheckData (void) { } else {} } - //set d_pace to 1000 for gene expression + // 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 + // For case-control studies, count # cases and # controls. int flag_cc=0; if (a_mode==13) { ni_case=0; @@ -920,53 +1134,88 @@ void PARAM::CheckData (void) { 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;} + 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 + // Set parameters for BSLMM and check for predict. if (a_mode==11 || a_mode==12 || a_mode==13 || a_mode==14) { - if (a_mode==11) {n_mh=1;} - if (logp_min==0) {logp_min=-1.0*log((double)ns_test);} + 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_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;} + 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) { 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; - error=true; - } else { - size_t flag_bv=0; - for (size_t i=0; i<(indicator_bv).size(); ++i) { - if (indicator_idv[i]!=indicator_bv[i]) {flag_bv++;} - } - if (flag_bv!=0) { - cout<<"error! individuals with missing value in the phenotype file does not match that in the estimated breeding value file: "<<flag_bv<<endl; - error=true; - } - } + 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; + error=true; + } else { + size_t flag_bv=0; + for (size_t i=0; i<(indicator_bv).size(); ++i) { + if (indicator_idv[i]!=indicator_bv[i]) {flag_bv++;} + } + if (flag_bv!=0) { + cout<<"error! individuals with missing value in the "<< + "phenotype file does not match that in the "<< + "estimated breeding value file: "<<flag_bv<<endl; + error=true; + } + } } } - if (a_mode==62 && !file_beta.empty() && mapRS2wcat.size()==0) {cout<<"vc analysis with beta files requires -wcat file."<<endl; error=true;} - if (a_mode==67 && mapRS2wcat.size()==0) {cout<<"ci analysis with beta files requires -wcat file."<<endl; error=true;} + if (a_mode==62 && !file_beta.empty() && mapRS2wcat.size()==0) { + cout<<"vc analysis with beta files requires -wcat file."<<endl; + error=true; + } + if (a_mode==67 && mapRS2wcat.size()==0) { + cout<<"ci analysis with beta files requires -wcat file."<<endl; + error=true; + } - //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;} + // 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; @@ -976,39 +1225,46 @@ void PARAM::PrintSummary () return; } - - 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;} + if (ReadFile_bed (file_str, indicator_idv, indicator_snp, + UtX, K, calc_K)==false) { + error=true; + } } else { - if (ReadFile_geno (file_geno, indicator_idv, indicator_snp, UtX, K, calc_K)==false) {error=true;} + 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) { +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;} + 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;} + 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; @@ -1016,24 +1272,33 @@ void PARAM::CalcKin (gsl_matrix *matrix_kin) { 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;} + 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;} + 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;} + if (BimbamKin (file_str, indicator_snp, a_mode-20, d_pace, + matrix_kin)==false) { + error=true; + } } return; } - - -//from an existing n by nd A and K matrices, compute the d by d S matrix (which is not necessary symmetric) -void compAKtoS (const gsl_matrix *A, const gsl_matrix *K, const size_t n_cvt, gsl_matrix *S) { +// From an existing n by nd A and K matrices, compute the d by d S +// matrix (which is not necessary symmetric). +void compAKtoS (const gsl_matrix *A, const gsl_matrix *K, const size_t n_cvt, + gsl_matrix *S) { size_t n_vc=S->size1, ni_test=A->size1; double di, dj, tr_AK, sum_A, sum_K, s_A, s_K, sum_AK, tr_A, tr_K, d; @@ -1070,17 +1335,16 @@ void compAKtoS (const gsl_matrix *A, const gsl_matrix *K, const size_t n_cvt, gs } } - //eigenlib_invert(Si); - //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 compKtoV -//map a number 1-(n_cvt+2) to an index between 0 and [(n_c+2)^2+(n_c+2)]/2-1 +// Copied from lmm.cpp; is used in the following function compKtoV +// 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;} + 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;} @@ -1091,8 +1355,10 @@ size_t GetabIndex (const size_t a, const size_t b, const size_t n_cvt) { return index; } -//from an existing n by nd (centered) G matrix, compute the d+1 by d*(d-1)/2*(d+1) Q matrix -//where inside i'th d+1 by d+1 matrix, each element is tr(KiKlKjKm)-r*tr(KmKiKl)-r*tr(KlKjKm)+r^2*tr(KlKm), where r=n/(n-1) +// From an existing n by nd (centered) G matrix, compute the d+1 by +// d*(d-1)/2*(d+1) Q matrix where inside i'th d+1 by d+1 matrix, each +// element is tr(KiKlKjKm)-r*tr(KmKiKl)-r*tr(KlKjKm)+r^2*tr(KlKm), +// where r=n/(n-1) void compKtoV (const gsl_matrix *G, gsl_matrix *V) { size_t n_vc=G->size2/G->size1, ni_test=G->size1; @@ -1103,27 +1369,24 @@ void compKtoV (const gsl_matrix *G, gsl_matrix *V) { double d, tr, r=(double)ni_test/(double)(ni_test-1); size_t t, t_il, t_jm, t_lm, t_im, t_jl, t_ij; - //compute KiKj for all pairs of i and j (not including the identity matrix) + // Compute KiKj for all pairs of i and j (not 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); + 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); + 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 + // Compute trKi, trKiKj. t=0; for (size_t i=0; i<n_vc; i++) { for (size_t j=i; j<n_vc; j++) { @@ -1143,7 +1406,7 @@ void compKtoV (const gsl_matrix *G, gsl_matrix *V) { gsl_vector_set (trKi, i, tr); } - //compute V + // Compute V. for (size_t i=0; i<n_vc; i++) { for (size_t j=i; j<n_vc; j++) { t_ij=GetabIndex (i+1, j+1, n_vc-2); @@ -1153,16 +1416,21 @@ void compKtoV (const gsl_matrix *G, gsl_matrix *V) { t_il=GetabIndex (i+1, l+1, n_vc-2); t_jm=GetabIndex (j+1, m+1, n_vc-2); t_lm=GetabIndex (l+1, m+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 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 KjKm_row=gsl_matrix_const_subrow (KiKj, k, t_jm*ni_test, ni_test); - gsl_vector_const_view KjKm_col=gsl_matrix_const_column (KiKj, t_jm*ni_test+k); - - gsl_vector_const_view Kl_row=gsl_matrix_const_subrow (G, k, l*ni_test, ni_test); - gsl_vector_const_view Km_row=gsl_matrix_const_subrow (G, k, m*ni_test, ni_test); + 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 KjKm_row= + gsl_matrix_const_subrow (KiKj, k, t_jm*ni_test, ni_test); + gsl_vector_const_view KjKm_col= + gsl_matrix_const_column (KiKj, t_jm*ni_test+k); + + gsl_vector_const_view Kl_row= + gsl_matrix_const_subrow (G, k, l*ni_test, ni_test); + gsl_vector_const_view Km_row= + gsl_matrix_const_subrow (G, k, m*ni_test, ni_test); if (i<=l && j<=m) { gsl_blas_ddot (&KiKl_row.vector, &KjKm_col.vector, &d); @@ -1201,9 +1469,12 @@ void compKtoV (const gsl_matrix *G, gsl_matrix *V) { t_jl=GetabIndex (j+1, l+1, n_vc-2); tr=0; for (size_t k=0; k<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 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); if (i<=l) { gsl_blas_ddot (&KiKl_row.vector, &Kj_row.vector, &d); @@ -1213,15 +1484,19 @@ void compKtoV (const gsl_matrix *G, gsl_matrix *V) { tr+=d; } } - tr+=-r*gsl_vector_get (trKiKj, t_il)-r*gsl_vector_get (trKiKj, t_jl)+r*r*gsl_vector_get (trKi, l); + tr+=-r*gsl_vector_get (trKiKj, t_il) - + r*gsl_vector_get (trKiKj, t_jl)+r*r*gsl_vector_get (trKi, l); } else if (l==n_vc && m!=n_vc) { t_jm=GetabIndex (j+1, m+1, n_vc-2); t_im=GetabIndex (i+1, m+1, n_vc-2); tr=0; for (size_t k=0; k<ni_test; k++) { - gsl_vector_const_view KjKm_row=gsl_matrix_const_subrow (KiKj, k, t_jm*ni_test, ni_test); - gsl_vector_const_view KjKm_col=gsl_matrix_const_column (KiKj, t_jm*ni_test+k); - gsl_vector_const_view Ki_row=gsl_matrix_const_subrow (G, k, i*ni_test, ni_test); + gsl_vector_const_view KjKm_row= + gsl_matrix_const_subrow (KiKj, k, t_jm*ni_test, ni_test); + gsl_vector_const_view KjKm_col= + gsl_matrix_const_column (KiKj, t_jm*ni_test+k); + gsl_vector_const_view Ki_row= + gsl_matrix_const_subrow (G, k, i*ni_test, ni_test); if (j<=m) { gsl_blas_ddot (&KjKm_row.vector, &Ki_row.vector, &d); @@ -1231,9 +1506,12 @@ void compKtoV (const gsl_matrix *G, gsl_matrix *V) { tr+=d; } } - tr+=-r*gsl_vector_get (trKiKj, t_im)-r*gsl_vector_get (trKiKj, t_jm)+r*r*gsl_vector_get (trKi, m); + tr+=-r*gsl_vector_get (trKiKj, t_im) - + r*gsl_vector_get (trKiKj, t_jm)+r*r*gsl_vector_get (trKi, m); } else { - tr=gsl_vector_get (trKiKj, t_ij)-r*gsl_vector_get (trKi, i)-r*gsl_vector_get (trKi, j)+r*r*(double)(ni_test-1); + tr=gsl_vector_get (trKiKj, t_ij) - + r*gsl_vector_get (trKi, i) - + r*gsl_vector_get (trKi, j)+r*r*(double)(ni_test-1); } gsl_matrix_set (V, l, t_ij*(n_vc+1)+m, tr); @@ -1251,10 +1529,9 @@ void compKtoV (const gsl_matrix *G, gsl_matrix *V) { return; } - - -//perform Jacknife sampling for variance of S -void JackknifeAKtoS (const gsl_matrix *W, const gsl_matrix *A, const gsl_matrix *K, gsl_matrix *S, gsl_matrix *Svar) { +// Perform Jacknife sampling for variance of S. +void JackknifeAKtoS (const gsl_matrix *W, const gsl_matrix *A, + const gsl_matrix *K, gsl_matrix *S, gsl_matrix *Svar) { size_t n_vc=Svar->size1, ni_test=A->size1, n_cvt=W->size2; vector<vector<vector<double> > > trAK, sumAK; @@ -1262,10 +1539,7 @@ void JackknifeAKtoS (const gsl_matrix *W, const gsl_matrix *A, const gsl_matrix vector<double> vec_tmp; double di, dj, d, m, v; - //gsl_matrix *Stmp=gsl_matrix_alloc (n_vc, ni_test*n_vc); - //gsl_matrix *Stmp_sub=gsl_matrix_alloc (n_vc, n_vc); - - //initialize and set all elements to zero + // Initialize and set all elements to zero. for (size_t i=0; i<ni_test; i++) { vec_tmp.push_back(0); } @@ -1284,7 +1558,7 @@ void JackknifeAKtoS (const gsl_matrix *W, const gsl_matrix *A, const gsl_matrix sumAK.push_back(sumK); } - //run jackknife + // Run jackknife. 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++) { @@ -1344,8 +1618,6 @@ void JackknifeAKtoS (const gsl_matrix *W, const gsl_matrix *A, const gsl_matrix d/=(trA[i][t]-sumA[i][t])*(trK[j][t]-sumK[j][t]); d-=1/(double)(ni_test-n_cvt-1); } - //gsl_matrix_set (Stmp, i, t*n_vc+j, d); - //gsl_matrix_set (Stmp, j, t*n_vc+i, d); m+=d; v+=d*d; } m/=(double)ni_test; @@ -1361,76 +1633,73 @@ void JackknifeAKtoS (const gsl_matrix *W, const gsl_matrix *A, const gsl_matrix } } - /* - for (size_t t=0; t<ni_test; t++) { - gsl_matrix_view Stmp_view=gsl_matrix_submatrix(Stmp, 0, t*n_vc, n_vc, n_vc); - gsl_matrix_memcpy (Stmp_sub, &Stmp_view.matrix); - eigenlib_invert(Stmp_sub); - gsl_matrix_memcpy (&Stmp_view.matrix, Stmp_sub); - } - - 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<ni_test; t++) { - d=gsl_matrix_get (Stmp, i, t*n_vc+j); - 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 (Si, i, j); - d=(double)ni_test*d-(double)(ni_test-1)*m; - gsl_matrix_set (Si, i, j, d); - if (i!=j) {gsl_matrix_set (Svar, j, i, v); gsl_matrix_set (Si, j, i, d);} - } - } - - gsl_matrix_free (Stmp); - */ 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 (const map<string, double> &mapRS2wA, const map<string, double> &mapRS2wK, const gsl_matrix *W, gsl_matrix *A, gsl_matrix *K, gsl_matrix *S, gsl_matrix *Svar, gsl_vector *ns) { +// 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 (const map<string, double> &mapRS2wA, + const map<string, double> &mapRS2wK, + const gsl_matrix *W, gsl_matrix *A, + gsl_matrix *K, gsl_matrix *S, + gsl_matrix *Svar, gsl_vector *ns) { string file_str; gsl_matrix_set_zero (S); gsl_matrix_set_zero (Svar); gsl_vector_set_zero (ns); - //compute the kinship matrix G for multiple categories; these matrices are not centered, for convienence of Jacknife sampling + // Compute the kinship matrix G for multiple categories; these + // matrices are not centered, for convienence of Jacknife sampling. if (!file_bfile.empty() ) { file_str=file_bfile+".bed"; if (mapRS2wA.size()==0) { - if (PlinkKin (file_str, d_pace, indicator_idv, indicator_snp, mapRS2wK, mapRS2cat, snpInfo, W, K, ns)==false) {error=true;} + if (PlinkKin (file_str, d_pace, indicator_idv, indicator_snp, mapRS2wK, + mapRS2cat, snpInfo, W, K, ns)==false) { + error=true; + } } else { - if (PlinkKin (file_str, d_pace, indicator_idv, indicator_snp, mapRS2wA, mapRS2cat, snpInfo, W, A, ns)==false) {error=true;} + if (PlinkKin (file_str, d_pace, indicator_idv, indicator_snp, mapRS2wA, + mapRS2cat, snpInfo, W, A, ns)==false) { + error=true; + } } } else if (!file_geno.empty()) { file_str=file_geno; if (mapRS2wA.size()==0) { - if (BimbamKin (file_str, d_pace, indicator_idv, indicator_snp, mapRS2wK, mapRS2cat, snpInfo, W, K, ns)==false) {error=true;} + if (BimbamKin (file_str, d_pace, indicator_idv, indicator_snp, + mapRS2wK, mapRS2cat, snpInfo, W, K, ns)==false) { + error=true; + } } else { - if (BimbamKin (file_str, d_pace, indicator_idv, indicator_snp, mapRS2wA, mapRS2cat, snpInfo, W, A, ns)==false) {error=true;} + if (BimbamKin (file_str, d_pace, indicator_idv, indicator_snp, + mapRS2wA, mapRS2cat, snpInfo, W, A, ns)==false) { + error=true; + } } } else if (!file_mbfile.empty() ){ if (mapRS2wA.size()==0) { - if (MFILEKin (1, file_mbfile, d_pace, indicator_idv, mindicator_snp, mapRS2wK, mapRS2cat, msnpInfo, W, K, ns)==false) {error=true;} + if (MFILEKin (1, file_mbfile, d_pace, indicator_idv, mindicator_snp, + mapRS2wK, mapRS2cat, msnpInfo, W, K, ns)==false) { + error=true; + } } else { - if (MFILEKin (1, file_mbfile, d_pace, indicator_idv, mindicator_snp, mapRS2wA, mapRS2cat, msnpInfo, W, A, ns)==false) {error=true;} + if (MFILEKin (1, file_mbfile, d_pace, indicator_idv, mindicator_snp, + mapRS2wA, mapRS2cat, msnpInfo, W, A, ns)==false) { + error=true; + } } } else if (!file_mgeno.empty()) { if (mapRS2wA.size()==0) { - if (MFILEKin (0, file_mgeno, d_pace, indicator_idv, mindicator_snp, mapRS2wK, mapRS2cat, msnpInfo, W, K, ns)==false) {error=true;} + if (MFILEKin (0, file_mgeno, d_pace, indicator_idv, mindicator_snp, + mapRS2wK, mapRS2cat, msnpInfo, W, K, ns)==false) { + error=true; + } } else { - if (MFILEKin (0, file_mgeno, d_pace, indicator_idv, mindicator_snp, mapRS2wA, mapRS2cat, msnpInfo, W, A, ns)==false) {error=true;} + if (MFILEKin (0, file_mgeno, d_pace, indicator_idv, mindicator_snp, + mapRS2wA, mapRS2cat, msnpInfo, W, A, ns)==false) { + error=true; + } } } @@ -1438,33 +1707,28 @@ void PARAM::CalcS (const map<string, double> &mapRS2wA, const map<string, double gsl_matrix_memcpy (A, K); } - //center and scale every kinship matrix inside G + // Center and scale every kinship matrix inside G. for (size_t i=0; i<n_vc; i++) { - gsl_matrix_view Ksub=gsl_matrix_submatrix(K, 0, i*ni_test, ni_test, ni_test); + gsl_matrix_view Ksub=gsl_matrix_submatrix(K,0,i*ni_test,ni_test,ni_test); CenterMatrix(&Ksub.matrix); ScaleMatrix(&Ksub.matrix); - gsl_matrix_view Asub=gsl_matrix_submatrix(A, 0, i*ni_test, ni_test, ni_test); + gsl_matrix_view Asub=gsl_matrix_submatrix(A,0,i*ni_test,ni_test,ni_test); CenterMatrix(&Asub.matrix); ScaleMatrix(&Asub.matrix); } - //based on G, compute S + // Cased on G, compute S. compAKtoS (A, K, W->size2, S); - //compute Svar and update S with Jacknife + // Compute Svar and update S with Jacknife. JackknifeAKtoS (W, A, K, S, Svar); - //based on G, compute a matrix Q that can be used to calculate the variance of q - //compKtoV (G, V); - return; } - - -void PARAM::WriteVector (const gsl_vector *q, const gsl_vector *s, const size_t n_total, const string suffix) -{ +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+="."; @@ -1472,7 +1736,10 @@ void PARAM::WriteVector (const gsl_vector *q, const gsl_vector *s, const size_t file_str+=".txt"; ofstream outfile (file_str.c_str(), ofstream::out); - if (!outfile) {cout<<"error writing file: "<<file_str.c_str()<<endl; return;} + if (!outfile) { + cout<<"error writing file: "<<file_str.c_str()<<endl; + return; + } outfile.precision(10); @@ -1491,10 +1758,7 @@ void PARAM::WriteVector (const gsl_vector *q, const gsl_vector *s, const size_t return; } - - -void PARAM::WriteVar (const string suffix) -{ +void PARAM::WriteVar (const string suffix) { string file_str, rs; file_str=path_out+"/"+file_out; file_str+="."; @@ -1502,7 +1766,10 @@ void PARAM::WriteVar (const string 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;} + if (!outfile) { + cout<<"error writing file: "<<file_str.c_str()<<endl; + return; + } outfile.precision(10); @@ -1528,11 +1795,7 @@ void PARAM::WriteVar (const string suffix) 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+="."; @@ -1540,7 +1803,10 @@ void PARAM::WriteMatrix (const gsl_matrix *matrix_U, const string suffix) file_str+=".txt"; ofstream outfile (file_str.c_str(), ofstream::out); - if (!outfile) {cout<<"error writing file: "<<file_str.c_str()<<endl; return;} + if (!outfile) { + cout<<"error writing file: "<<file_str.c_str()<<endl; + return; + } outfile.precision(10); @@ -1556,9 +1822,7 @@ void PARAM::WriteMatrix (const gsl_matrix *matrix_U, const string suffix) 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+="."; @@ -1566,7 +1830,10 @@ void PARAM::WriteVector (const gsl_vector *vector_D, const string suffix) file_str+=".txt"; ofstream outfile (file_str.c_str(), ofstream::out); - if (!outfile) {cout<<"error writing file: "<<file_str.c_str()<<endl; return;} + if (!outfile) { + cout<<"error writing file: "<<file_str.c_str()<<endl; + return; + } outfile.precision(10); @@ -1579,9 +1846,7 @@ void PARAM::WriteVector (const gsl_vector *vector_D, const string suffix) return; } - -void PARAM::CheckCvt () -{ +void PARAM::CheckCvt () { if (indicator_cvt.size()==0) {return;} size_t ci_test=0; @@ -1600,22 +1865,25 @@ void PARAM::CheckCvt () double v_min, v_max; set<size_t> set_remove; - //check if any columns is an intercept + // 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 + // Add an intecept term if needed. if (n_cvt==set_remove.size()) { indicator_cvt.clear(); n_cvt=1; } else if (flag_ipt==0) { - cout<<"no intecept term is found in the cvt file. a column of 1s is added."<<endl; + cout<<"no intecept term is found in the cvt file. "<< + "a column of 1s is added."<<endl; for (vector<int>::size_type i=0; i<indicator_idv.size(); ++i) { - if (indicator_idv[i]==0 || indicator_cvt[i]==0) {continue;} - cvt[i].push_back(1.0); + if (indicator_idv[i]==0 || indicator_cvt[i]==0) { + continue; + } + cvt[i].push_back(1.0); } n_cvt++; @@ -1626,11 +1894,10 @@ void PARAM::CheckCvt () return; } - -//post-process phentoypes, covariates -void PARAM::ProcessCvtPhen () -{ - //convert indicator_pheno to indicator_idv +// Post-process phentoypes and covariates. +void PARAM::ProcessCvtPhen () { + + // Convert indicator_pheno to indicator_idv. int k=1; indicator_idv.clear(); for (size_t i=0; i<indicator_pheno.size(); i++) { @@ -1641,42 +1908,49 @@ void PARAM::ProcessCvtPhen () indicator_idv.push_back(k); } - //remove individuals with missing covariates + // Remove individuals with missing covariates. if ((indicator_cvt).size()!=0) { - for (vector<int>::size_type i=0; i<(indicator_idv).size(); ++i) { + for (vector<int>::size_type i=0; + i<(indicator_idv).size(); + ++i) { indicator_idv[i]*=indicator_cvt[i]; } } - //remove individuals with missing gxe variables + // Remove individuals with missing gxe variables. if ((indicator_gxe).size()!=0) { - for (vector<int>::size_type i=0; i<(indicator_idv).size(); ++i) { + for (vector<int>::size_type i=0; + i<(indicator_idv).size(); + ++i) { indicator_idv[i]*=indicator_gxe[i]; } } - //remove individuals with missing residual weights + // Remove individuals with missing residual weights. if ((indicator_weight).size()!=0) { - for (vector<int>::size_type i=0; i<(indicator_idv).size(); ++i) { + for (vector<int>::size_type i=0; + i<(indicator_idv).size(); + ++i) { indicator_idv[i]*=indicator_weight[i]; } } - //obtain ni_test + // Obtain ni_test. ni_test=0; for (vector<int>::size_type i=0; i<(indicator_idv).size(); ++i) { if (indicator_idv[i]==0) {continue;} ni_test++; } - - - //if subsample number is set, perform a random sub-sampling to determine the subsampled ids + // 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; + cout<<"error! number of subsamples is less than number of"<< + "analyzed individuals. "<<endl; } else { - //set up random environment + + // Set up random environment. gsl_rng_env_setup(); gsl_rng *gsl_r; const gsl_rng_type * gslType; @@ -1686,12 +1960,13 @@ void PARAM::ProcessCvtPhen () time (&rawtime); tm * ptm = gmtime (&rawtime); - randseed = (unsigned) (ptm->tm_hour%24*3600+ptm->tm_min*60+ptm->tm_sec); + 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 + // 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); @@ -1700,9 +1975,10 @@ void PARAM::ProcessCvtPhen () 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) ); + 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 + // 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;} @@ -1715,25 +1991,27 @@ void PARAM::ProcessCvtPhen () } } - //check ni_test + // Check ni_test. if (ni_test==0 && a_mode!=15) { 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 + // 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. if (indicator_cvt.size()!=0) { CheckCvt(); } else { vector<double> cvt_row; cvt_row.push_back(1); - for (vector<int>::size_type i=0; i<(indicator_idv).size(); ++i) { + for (vector<int>::size_type i=0; + i<(indicator_idv).size(); + ++i) { indicator_cvt.push_back(1); - cvt.push_back(cvt_row); } } @@ -1741,11 +2019,7 @@ void PARAM::ProcessCvtPhen () 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) { @@ -1759,57 +2033,7 @@ void PARAM::CopyCvt (gsl_matrix *W) return; } -/* -//if there is a column contains all 1's, then flag==1; otherwise flag=0 -void PARAM::CopyA (size_t flag, gsl_matrix *A) -{ - size_t ci_test=0; - string rs; - vector<size_t> flag_vec; - vector<double> catc; - - for (size_t j=0; j<n_cat; j++) { - flag_vec.push_back(0); - } - - for (vector<int>::size_type i=0; i<indicator_snp.size(); ++i) { - if (indicator_snp[i]==0) {continue;} - rs=snpInfo[i].rs_number; - - if (mapRS2catc.count(rs)==0) { - for (size_t j=0; j<n_cat; j++) { - gsl_matrix_set (A, ci_test, j, 0); - } - } else { - for (size_t j=0; j<n_cat; j++) { - gsl_matrix_set (A, ci_test, j, mapRS2catc[rs][j]); - } - } - - if (ci_test==0) { - for (size_t j=0; j<n_cat; j++) { - catc.push_back(mapRS2catc[rs][j]); - } - } else { - for (size_t j=0; j<n_cat; j++) { - if (catc[j]==mapRS2catc[rs][j]) {flag_vec[j]++;}; - } - } - - ci_test++; - } - - flag=0; - for (size_t j=0; j<n_cat; j++) { - if (flag_vec[j]==0) {flag++;} - } - - return; -} -*/ - -void PARAM::CopyGxe (gsl_vector *env) -{ +void PARAM::CopyGxe (gsl_vector *env) { size_t ci_test=0; for (vector<int>::size_type i=0; i<indicator_idv.size(); ++i) { @@ -1821,8 +2045,7 @@ void PARAM::CopyGxe (gsl_vector *env) return; } -void PARAM::CopyWeight (gsl_vector *w) -{ +void PARAM::CopyWeight (gsl_vector *w) { size_t ci_test=0; for (vector<int>::size_type i=0; i<indicator_idv.size(); ++i) { @@ -1834,11 +2057,9 @@ void PARAM::CopyWeight (gsl_vector *w) 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) -{ +// 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) { size_t ci_test=0; for (vector<int>::size_type i=0; i<indicator_idv.size(); ++i) { @@ -1859,10 +2080,9 @@ void PARAM::CopyCvtPhen (gsl_matrix *W, gsl_vector *y, size_t flag) 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) -{ +// 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) { size_t ci_test=0; for (vector<int>::size_type i=0; i<indicator_idv.size(); ++i) { @@ -1885,12 +2105,7 @@ void PARAM::CopyCvtPhen (gsl_matrix *W, gsl_matrix *Y, size_t flag) return; } - - - - -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) { @@ -1902,10 +2117,8 @@ void PARAM::CopyRead (gsl_vector *log_N) return; } - - -void PARAM::ObtainWeight (const set<string> &setSnps_beta, map<string, double> &mapRS2wK) -{ +void PARAM::ObtainWeight (const set<string> &setSnps_beta, + map<string, double> &mapRS2wK) { mapRS2wK.clear(); vector<double> wsum, wcount; @@ -1921,7 +2134,10 @@ void PARAM::ObtainWeight (const set<string> &setSnps_beta, map<string, double> & if (indicator_snp[i]==0) {continue;} rs=snpInfo[i].rs_number; - if ( (setSnps_beta.size()==0 || setSnps_beta.count(rs)!=0) && (mapRS2wsnp.size()==0 || mapRS2wsnp.count(rs)!=0) && (mapRS2wcat.size()==0 || mapRS2wcat.count(rs)!=0) && (mapRS2cat.size()==0 || mapRS2cat.count(rs)!=0) ) { + if ( (setSnps_beta.size()==0 || setSnps_beta.count(rs)!=0) && + (mapRS2wsnp.size()==0 || mapRS2wsnp.count(rs)!=0) && + (mapRS2wcat.size()==0 || mapRS2wcat.count(rs)!=0) && + (mapRS2cat.size()==0 || mapRS2cat.count(rs)!=0) ) { if (mapRS2wsnp.size()!=0) { mapRS2wK[rs]=mapRS2wsnp[rs]; if (mapRS2cat.size()==0) { @@ -1945,7 +2161,10 @@ void PARAM::ObtainWeight (const set<string> &setSnps_beta, map<string, double> & if (indicator_snp[i]==0) {continue;} rs=snpInfo[i].rs_number; - if ( (setSnps_beta.size()==0 || setSnps_beta.count(rs)!=0) && (mapRS2wsnp.size()==0 || mapRS2wsnp.count(rs)!=0) && (mapRS2wcat.size()==0 || mapRS2wcat.count(rs)!=0) && (mapRS2cat.size()==0 || mapRS2cat.count(rs)!=0) ) { + if ((setSnps_beta.size()==0 || setSnps_beta.count(rs)!=0) && + (mapRS2wsnp.size()==0 || mapRS2wsnp.count(rs)!=0) && + (mapRS2wcat.size()==0 || mapRS2wcat.count(rs)!=0) && + (mapRS2cat.size()==0 || mapRS2cat.count(rs)!=0) ) { if (mapRS2wsnp.size()!=0) { mapRS2wK[rs]=mapRS2wsnp[rs]; if (mapRS2cat.size()==0) { @@ -1967,7 +2186,9 @@ void PARAM::ObtainWeight (const set<string> &setSnps_beta, map<string, double> & wsum[i]/=wcount[i]; } - for (map<string, double>::iterator it=mapRS2wK.begin(); it!=mapRS2wK.end(); ++it) { + for (map<string, double>::iterator it=mapRS2wK.begin(); + it!=mapRS2wK.end(); + ++it) { if (mapRS2cat.size()==0) { it->second/=wsum[0]; } else { @@ -1978,10 +2199,12 @@ void PARAM::ObtainWeight (const set<string> &setSnps_beta, map<string, double> & return; } - -//pve_flag=0 then do not change pve; pve_flag==1, then change pve to 0 if pve < 0 and pve to 1 if pve > 1 -void PARAM::UpdateWeight (const size_t pve_flag, const map<string, double> &mapRS2wK, const size_t ni_test, const gsl_vector *ns, map<string, double> &mapRS2wA) -{ +// If pve_flag=0 then do not change pve; pve_flag==1, then change pve +// to 0 if pve < 0 and pve to 1 if pve > 1. +void PARAM::UpdateWeight (const size_t pve_flag, + const map<string, double> &mapRS2wK, + const size_t ni_test, const gsl_vector *ns, + map<string, double> &mapRS2wA) { double d; vector<double> wsum, wcount; @@ -1990,7 +2213,9 @@ void PARAM::UpdateWeight (const size_t pve_flag, const map<string, double> &mapR wcount.push_back(0.0); } - for (map<string, double>::const_iterator it=mapRS2wK.begin(); it!=mapRS2wK.end(); ++it) { + for (map<string, double>::const_iterator it=mapRS2wK.begin(); + it!=mapRS2wK.end(); + ++it) { d=1; for (size_t i=0; i<n_vc; i++) { if (v_pve[i]>=1 && pve_flag==1) { @@ -1998,7 +2223,8 @@ void PARAM::UpdateWeight (const size_t pve_flag, const map<string, double> &mapR } else if (v_pve[i]<=0 && pve_flag==1) { d+=0; } else { - d+=(double)ni_test/gsl_vector_get(ns, i)*mapRS2wcat[it->first][i]*v_pve[i]; + d+=(double)ni_test/gsl_vector_get(ns, i)* + mapRS2wcat[it->first][i]*v_pve[i]; } } mapRS2wA[it->first]=1/(d*d); @@ -2016,7 +2242,9 @@ void PARAM::UpdateWeight (const size_t pve_flag, const map<string, double> &mapR wsum[i]/=wcount[i]; } - for (map<string, double>::iterator it=mapRS2wA.begin(); it!=mapRS2wA.end(); ++it) { + for (map<string, double>::iterator it=mapRS2wA.begin(); + it!=mapRS2wA.end(); + ++it) { if (mapRS2cat.size()==0) { it->second/=wsum[0]; } else { @@ -2026,9 +2254,13 @@ void PARAM::UpdateWeight (const size_t pve_flag, const map<string, double> &mapR return; } -// this function updates indicator_snp, and save z-scores and other values into vectors -void PARAM::UpdateSNPnZ (const map<string, double> &mapRS2wA, const map<string, string> &mapRS2A1, const map<string, double> &mapRS2z, gsl_vector *w, gsl_vector *z, vector<size_t> &vec_cat) -{ +// This function updates indicator_snp, and save z-scores and other +// values into vectors. +void PARAM::UpdateSNPnZ (const map<string, double> &mapRS2wA, + const map<string, string> &mapRS2A1, + const map<string, double> &mapRS2z, + gsl_vector *w, gsl_vector *z, + vector<size_t> &vec_cat) { gsl_vector_set_zero (w); gsl_vector_set_zero (z); vec_cat.clear(); @@ -2086,11 +2318,9 @@ void PARAM::UpdateSNPnZ (const map<string, double> &mapRS2wA, const map<string, return; } - - -// this function updates indicator_snp, and save z-scores and other values into vectors -void PARAM::UpdateSNP (const map<string, double> &mapRS2wA) -{ +// This function updates indicator_snp, and save z-scores and other +// values into vectors. +void PARAM::UpdateSNP (const map<string, double> &mapRS2wA) { string rs; if (msnpInfo.size()==0) { for (size_t i=0; i<snpInfo.size(); i++) { |