diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/bslmmdap.cpp | 1015 | ||||
| -rw-r--r-- | src/bslmmdap.h | 101 | ||||
| -rw-r--r-- | src/ldr.cpp | 285 | ||||
| -rw-r--r-- | src/ldr.h | 77 | ||||
| -rw-r--r-- | src/logistic.cpp | 783 | ||||
| -rw-r--r-- | src/logistic.h | 70 | ||||
| -rw-r--r-- | src/varcov.cpp | 482 | ||||
| -rw-r--r-- | src/varcov.h | 72 |
8 files changed, 2885 insertions, 0 deletions
diff --git a/src/bslmmdap.cpp b/src/bslmmdap.cpp new file mode 100644 index 0000000..0bf0e7b --- /dev/null +++ b/src/bslmmdap.cpp @@ -0,0 +1,1015 @@ +/* + Genome-wide Efficient Mixed Model Association (GEMMA) + Copyright (C) 2011 Xiang Zhou + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include <iostream> +#include <fstream> +#include <sstream> + +#include <iomanip> +#include <cmath> +#include <iostream> +#include <stdio.h> +#include <stdlib.h> +#include <ctime> +#include <cstring> +#include <algorithm> + +#include "gsl/gsl_vector.h" +#include "gsl/gsl_matrix.h" +#include "gsl/gsl_linalg.h" +#include "gsl/gsl_blas.h" +#include "gsl/gsl_eigen.h" +#include "gsl/gsl_randist.h" +#include "gsl/gsl_cdf.h" +#include "gsl/gsl_roots.h" + + + +#include "logistic.h" +#include "lapack.h" + +#ifdef FORCE_FLOAT +#include "param_float.h" +#include "bslmmdap_float.h" +#include "lmm_float.h" //for class FUNC_PARAM and MatrixCalcLR +#include "lm_float.h" +#include "mathfunc_float.h" //for function CenterVector +#else +#include "param.h" +#include "bslmmdap.h" +#include "lmm.h" +#include "lm.h" +#include "mathfunc.h" +#endif + +using namespace std; + + + + +void BSLMMDAP::CopyFromParam (PARAM &cPar) +{ + file_out=cPar.file_out; + path_out=cPar.path_out; + + time_UtZ=0.0; + time_Omega=0.0; + + h_min=cPar.h_min; + h_max=cPar.h_max; + h_ngrid=cPar.h_ngrid; + rho_min=cPar.rho_min; + rho_max=cPar.rho_max; + rho_ngrid=cPar.rho_ngrid; + + if (h_min<=0) {h_min=0.01;} + if (h_max>=1) {h_max=0.99;} + if (rho_min<=0) {rho_min=0.01;} + if (rho_max>=1) {rho_max=0.99;} + + trace_G=cPar.trace_G; + + ni_total=cPar.ni_total; + ns_total=cPar.ns_total; + ni_test=cPar.ni_test; + ns_test=cPar.ns_test; + + indicator_idv=cPar.indicator_idv; + indicator_snp=cPar.indicator_snp; + snpInfo=cPar.snpInfo; + + return; +} + + +void BSLMMDAP::CopyToParam (PARAM &cPar) +{ + cPar.time_UtZ=time_UtZ; + cPar.time_Omega=time_Omega; + + return; +} + + + +//read hyp file +void ReadFile_hyb (const string &file_hyp, vector<double> &vec_sa2, vector<double> &vec_sb2, vector<double> &vec_wab) +{ + vec_sa2.clear(); vec_sb2.clear(); vec_wab.clear(); + + igzstream infile (file_hyp.c_str(), igzstream::in); + if (!infile) {cout<<"error! fail to open hyp file: "<<file_hyp<<endl; return;} + + string line; + char *ch_ptr; + + getline(infile, line); + + while (!safeGetline(infile, line).eof()) { + ch_ptr=strtok ((char *)line.c_str(), " , \t"); + ch_ptr=strtok (NULL, " , \t"); + + ch_ptr=strtok (NULL, " , \t"); + vec_sa2.push_back(atof(ch_ptr)); + + ch_ptr=strtok (NULL, " , \t"); + vec_sb2.push_back(atof(ch_ptr)); + + ch_ptr=strtok (NULL, " , \t"); + vec_wab.push_back(atof(ch_ptr)); + } + + infile.close(); + infile.clear(); + + return; +} + + +//read bf file +void ReadFile_bf (const string &file_bf, vector<string> &vec_rs, vector<vector<vector<double> > > &BF) +{ + BF.clear(); vec_rs.clear(); + + igzstream infile (file_bf.c_str(), igzstream::in); + if (!infile) {cout<<"error! fail to open bf file: "<<file_bf<<endl; return;} + + string line, rs, block; + vector<double> vec_bf; + vector<vector<double> > mat_bf; + char *ch_ptr; + + size_t bf_size, flag_block; + + getline(infile, line); + + size_t t=0; + while (!safeGetline(infile, line).eof()) { + flag_block=0; + + ch_ptr=strtok ((char *)line.c_str(), " , \t"); + rs=ch_ptr; + vec_rs.push_back(rs); + + ch_ptr=strtok (NULL, " , \t"); + if (t==0) { + block=ch_ptr; + } else { + if (strcmp(ch_ptr, block.c_str() )!=0) { + flag_block=1; + block=ch_ptr; + } + } + + ch_ptr=strtok (NULL, " , \t"); + while (ch_ptr!=NULL) { + vec_bf.push_back(atof(ch_ptr)); + ch_ptr=strtok (NULL, " , \t"); + } + + if (t==0) { + bf_size=vec_bf.size(); + } else { + if (bf_size!=vec_bf.size()) {cout<<"error! unequal row size in bf file."<<endl;} + } + + if (flag_block==0) { + mat_bf.push_back(vec_bf); + } else { + BF.push_back(mat_bf); + mat_bf.clear(); + } + vec_bf.clear(); + + t++; + } + + infile.close(); + infile.clear(); + + return; +} + + +//read category files +//read both continuous and discrete category file, record mapRS2catc +void ReadFile_cat (const string &file_cat, const vector<string> &vec_rs, gsl_matrix *Ac, gsl_matrix_int *Ad, gsl_vector_int *dlevel, size_t &kc, size_t &kd) +{ + igzstream infile (file_cat.c_str(), igzstream::in); + if (!infile) {cout<<"error! fail to open category file: "<<file_cat<<endl; return;} + + string line; + char *ch_ptr; + + string rs, chr, a1, a0, pos, cm; + + //read header + HEADER header; + !safeGetline(infile, line).eof(); + ReadHeader (line, header); + + //use the header to determine the number of categories + kc=header.catc_col.size(); kd=header.catd_col.size(); + + //set up storage and mapper + map<string, vector<double> > mapRS2catc; + map<string, vector<int> > mapRS2catd; + vector<double> catc; + vector<int> catd; + + //read the following lines to record mapRS2cat + while (!safeGetline(infile, line).eof()) { + ch_ptr=strtok ((char *)line.c_str(), " , \t"); + + if (header.rs_col==0) { + rs=chr+":"+pos; + } + + catc.clear(); catd.clear(); + + for (size_t i=0; i<header.coln; i++) { + if (header.rs_col!=0 && header.rs_col==i+1) { + rs=ch_ptr; + } else if (header.chr_col!=0 && header.chr_col==i+1) { + chr=ch_ptr; + } else if (header.pos_col!=0 && header.pos_col==i+1) { + pos=ch_ptr; + } else if (header.cm_col!=0 && header.cm_col==i+1) { + cm=ch_ptr; + } else if (header.a1_col!=0 && header.a1_col==i+1) { + a1=ch_ptr; + } else if (header.a0_col!=0 && header.a0_col==i+1) { + a0=ch_ptr; + } else if (header.catc_col.size()!=0 && header.catc_col.count(i+1)!=0 ) { + catc.push_back(atof(ch_ptr)); + } else if (header.catd_col.size()!=0 && header.catd_col.count(i+1)!=0 ) { + catd.push_back(atoi(ch_ptr)); + } else {} + + ch_ptr=strtok (NULL, " , \t"); + } + + if (mapRS2catc.count(rs)==0 && kc>0) {mapRS2catc[rs]=catc;} + if (mapRS2catd.count(rs)==0 && kd>0) {mapRS2catd[rs]=catd;} + } + + //load into Ad and Ac + if (kc>0) { + Ac=gsl_matrix_alloc(vec_rs.size(), kc); + for (size_t i=0; i<vec_rs.size(); i++) { + if (mapRS2catc.count(vec_rs[i])!=0) { + for (size_t j=0; j<kc; j++) { + gsl_matrix_set(Ac, i, j, mapRS2catc[vec_rs[i]][j]); + } + } else { + for (size_t j=0; j<kc; j++) { + gsl_matrix_set(Ac, i, j, 0); + } + } + } + } + + if (kd>0) { + Ad=gsl_matrix_int_alloc(vec_rs.size(), kd); + + for (size_t i=0; i<vec_rs.size(); i++) { + if (mapRS2catd.count(vec_rs[i])!=0) { + for (size_t j=0; j<kd; j++) { + gsl_matrix_int_set(Ad, i, j, mapRS2catd[vec_rs[i]][j]); + } + } else { + for (size_t j=0; j<kd; j++) { + gsl_matrix_int_set(Ad, i, j, 0); + } + } + } + + dlevel=gsl_vector_int_alloc(kd); + map<int, int> rcd; + int val; + for (size_t j=0; j<kd; j++) { + rcd.clear(); + for (size_t i=0; i<Ad->size1; i++) { + val = gsl_matrix_int_get(Ad, i, j); + rcd[val] = 1; + } + gsl_vector_int_set (dlevel, j, rcd.size()); + } + } + + infile.clear(); + infile.close(); + + return; +} + + + + + + + + +void BSLMMDAP::WriteResult (const gsl_matrix *Hyper, const gsl_matrix *BF) +{ + string file_bf, file_hyp; + file_bf=path_out+"/"+file_out; + file_bf+=".bf.txt"; + file_hyp=path_out+"/"+file_out; + file_hyp+=".hyp.txt"; + + ofstream outfile_bf, outfile_hyp; + + outfile_bf.open (file_bf.c_str(), ofstream::out); + outfile_hyp.open (file_hyp.c_str(), ofstream::out); + + if (!outfile_bf) {cout<<"error writing file: "<<file_bf<<endl; return;} + if (!outfile_hyp) {cout<<"error writing file: "<<file_hyp<<endl; return;} + + outfile_hyp<<"h"<<"\t"<<"rho"<<"\t"<<"sa2"<<"\t"<<"sb2"<<"\t"<<"weight"<<endl; + outfile_hyp<<scientific; + for (size_t i=0; i<Hyper->size1; i++) { + for (size_t j=0; j<Hyper->size2; j++) { + outfile_hyp<<setprecision(6)<<gsl_matrix_get (Hyper, i, j)<<"\t"; + } + outfile_hyp<<endl; + } + + outfile_bf<<"chr"<<"\t"<<"rs"<<"\t"<<"ps"<<"\t"<<"n_miss"; + for (size_t i=0; i<BF->size2; i++) { + outfile_bf<<"\t"<<"BF"<<i+1; + } + outfile_bf<<endl; + + size_t t=0; + for (size_t i=0; i<ns_total; ++i) { + if (indicator_snp[i]==0) {continue;} + + outfile_bf<<snpInfo[i].chr<<"\t"<<snpInfo[i].rs_number<<"\t" + <<snpInfo[i].base_position<<"\t"<<snpInfo[i].n_miss; + + outfile_bf<<scientific; + for (size_t j=0; j<BF->size2; j++) { + outfile_bf<<"\t"<<setprecision(6)<<gsl_matrix_get (BF, t, j); + } + outfile_bf<<endl; + + t++; + } + + outfile_hyp.close(); + outfile_hyp.clear(); + outfile_bf.close(); + outfile_bf.clear(); + return; +} + + + +void BSLMMDAP::WriteResult (const vector<string> &vec_rs, const gsl_matrix *Hyper, const gsl_vector *pip, const gsl_vector *coef) +{ + string file_gamma, file_hyp, file_coef; + file_gamma=path_out+"/"+file_out; + file_gamma+=".gamma.txt"; + file_hyp=path_out+"/"+file_out; + file_hyp+=".hyp.txt"; + file_coef=path_out+"/"+file_out; + file_coef+=".coef.txt"; + + ofstream outfile_gamma, outfile_hyp, outfile_coef; + + outfile_gamma.open (file_gamma.c_str(), ofstream::out); + outfile_hyp.open (file_hyp.c_str(), ofstream::out); + outfile_coef.open (file_coef.c_str(), ofstream::out); + + if (!outfile_gamma) {cout<<"error writing file: "<<file_gamma<<endl; return;} + if (!outfile_hyp) {cout<<"error writing file: "<<file_hyp<<endl; return;} + if (!outfile_coef) {cout<<"error writing file: "<<file_coef<<endl; return;} + + outfile_hyp<<"h"<<"\t"<<"rho"<<"\t"<<"sa2"<<"\t"<<"sb2"<<"\t"<<"weight"<<endl; + outfile_hyp<<scientific; + for (size_t i=0; i<Hyper->size1; i++) { + for (size_t j=0; j<Hyper->size2; j++) { + outfile_hyp<<setprecision(6)<<gsl_matrix_get (Hyper, i, j)<<"\t"; + } + outfile_hyp<<endl; + } + + + outfile_gamma<<"rs"<<"\t"<<"gamma"<<endl; + for (size_t i=0; i<vec_rs.size(); ++i) { + outfile_gamma<<vec_rs[i]<<"\t"<<scientific<<setprecision(6)<<gsl_vector_get(pip, i)<<endl; + } + + outfile_coef<<"coef"<<endl; + outfile_coef<<scientific; + for (size_t i=0; i<coef->size; i++) { + outfile_coef<<setprecision(6)<<gsl_vector_get (coef, i)<<endl; + } + + outfile_coef.close(); + outfile_coef.clear(); + outfile_hyp.close(); + outfile_hyp.clear(); + outfile_gamma.close(); + outfile_gamma.clear(); + return; +} + + + + +/* +void BSLMMDAP::SetXgamma (gsl_matrix *Xgamma, const gsl_matrix *X, vector<size_t> &rank) +{ + size_t pos; + for (size_t i=0; i<rank.size(); ++i) { + pos=mapRank2pos[rank[i]]; + gsl_vector_view Xgamma_col=gsl_matrix_column (Xgamma, i); + gsl_vector_const_view X_col=gsl_matrix_const_column (X, pos); + gsl_vector_memcpy (&Xgamma_col.vector, &X_col.vector); + } + + return; +} +*/ + +double BSLMMDAP::CalcMarginal (const gsl_vector *Uty, const gsl_vector *K_eval, const double sigma_b2, const double tau) +{ + gsl_vector *weight_Hi=gsl_vector_alloc (Uty->size); + + double logm=0.0; + double d, uy, Hi_yy=0, logdet_H=0.0; + for (size_t i=0; i<ni_test; ++i) { + d=gsl_vector_get (K_eval, i)*sigma_b2; + d=1.0/(d+1.0); + gsl_vector_set (weight_Hi, i, d); + + logdet_H-=log(d); + uy=gsl_vector_get (Uty, i); + Hi_yy+=d*uy*uy; + } + + //calculate likelihood + logm=-0.5*logdet_H-0.5*tau*Hi_yy+0.5*log(tau)*(double)ni_test; + + gsl_vector_free (weight_Hi); + + return logm; +} + + +double BSLMMDAP::CalcMarginal (const gsl_matrix *UtXgamma, const gsl_vector *Uty, const gsl_vector *K_eval, const double sigma_a2, const double sigma_b2, const double tau) +{ + clock_t time_start; + double logm=0.0; + double d, uy, P_yy=0, logdet_O=0.0, logdet_H=0.0; + + gsl_matrix *UtXgamma_eval=gsl_matrix_alloc (UtXgamma->size1, UtXgamma->size2); + gsl_matrix *Omega=gsl_matrix_alloc (UtXgamma->size2, UtXgamma->size2); + gsl_vector *XtHiy=gsl_vector_alloc (UtXgamma->size2); + gsl_vector *beta_hat=gsl_vector_alloc (UtXgamma->size2); + gsl_vector *weight_Hi=gsl_vector_alloc (UtXgamma->size1); + + gsl_matrix_memcpy (UtXgamma_eval, UtXgamma); + + logdet_H=0.0; P_yy=0.0; + for (size_t i=0; i<ni_test; ++i) { + gsl_vector_view UtXgamma_row=gsl_matrix_row (UtXgamma_eval, i); + d=gsl_vector_get (K_eval, i)*sigma_b2; + d=1.0/(d+1.0); + gsl_vector_set (weight_Hi, i, d); + + logdet_H-=log(d); + uy=gsl_vector_get (Uty, i); + P_yy+=d*uy*uy; + gsl_vector_scale (&UtXgamma_row.vector, d); + } + + //calculate Omega + gsl_matrix_set_identity (Omega); + + time_start=clock(); +#ifdef WITH_LAPACK + lapack_dgemm ((char *)"T", (char *)"N", sigma_a2, UtXgamma_eval, UtXgamma, 1.0, Omega); +#else + gsl_blas_dgemm (CblasTrans, CblasNoTrans, sigma_a2, UtXgamma_eval, UtXgamma, 1.0, Omega); +#endif + time_Omega+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0); + + //calculate beta_hat + gsl_blas_dgemv (CblasTrans, 1.0, UtXgamma_eval, Uty, 0.0, XtHiy); + + logdet_O=CholeskySolve(Omega, XtHiy, beta_hat); + + gsl_vector_scale (beta_hat, sigma_a2); + + gsl_blas_ddot (XtHiy, beta_hat, &d); + P_yy-=d; + + gsl_matrix_free (UtXgamma_eval); + gsl_matrix_free (Omega); + gsl_vector_free (XtHiy); + gsl_vector_free (beta_hat); + gsl_vector_free (weight_Hi); + + logm=-0.5*logdet_H-0.5*logdet_O-0.5*tau*P_yy+0.5*log(tau)*(double)ni_test; + + return logm; +} + + +double BSLMMDAP::CalcPrior (class HYPBSLMM &cHyp) { + double logprior=0; + logprior=((double)cHyp.n_gamma-1.0)*cHyp.logp+((double)ns_test-(double)cHyp.n_gamma)*log(1.0-exp(cHyp.logp)); + return logprior; +} + + +//where A is the ni_test by n_cat matrix of annotations +void BSLMMDAP::DAP_CalcBF (const gsl_matrix *U, const gsl_matrix *UtX, const gsl_vector *Uty, const gsl_vector *K_eval, const gsl_vector *y) { + clock_t time_start; + + //set up BF + double tau, h, rho, sigma_a2, sigma_b2, d; + size_t ns_causal=10; + size_t n_grid=h_ngrid*rho_ngrid; + vector<double> vec_sa2, vec_sb2, logm_null; + + gsl_matrix *BF=gsl_matrix_alloc(ns_test, n_grid); + gsl_matrix *Xgamma=gsl_matrix_alloc(ni_test, 1); + gsl_matrix *Hyper=gsl_matrix_alloc(n_grid, 5); + + //compute tau by using yty + gsl_blas_ddot (Uty, Uty, &tau); + tau=(double)ni_test/tau; + + //set up grid values for sigma_a2 and sigma_b2 based on an approximately even grid for h and rho, and a fixed number of causals + size_t ij=0; + for (size_t i=0; i<h_ngrid; i++) { + h=h_min+(h_max-h_min)*(double)i/((double)h_ngrid-1); + for (size_t j=0; j<rho_ngrid; j++) { + rho=rho_min+(rho_max-rho_min)*(double)j/((double)rho_ngrid-1); + + sigma_a2=h*rho/((1-h)*(double)ns_causal); + sigma_b2=h*(1.0-rho)/(trace_G*(1-h)); + + vec_sa2.push_back(sigma_a2); + vec_sb2.push_back(sigma_b2); + logm_null.push_back(CalcMarginal (Uty, K_eval, 0.0, tau)); + + gsl_matrix_set (Hyper, ij, 0, h); + gsl_matrix_set (Hyper, ij, 1, rho); + gsl_matrix_set (Hyper, ij, 2, sigma_a2); + gsl_matrix_set (Hyper, ij, 3, sigma_b2); + gsl_matrix_set (Hyper, ij, 4, 1/(double)n_grid); + ij++; + } + } + + //compute BF factors + time_start=clock(); + cout<<"Calculating BF..."<<endl; + for (size_t t=0; t<ns_test; t++) { + gsl_vector_view Xgamma_col=gsl_matrix_column (Xgamma, 0); + gsl_vector_const_view X_col=gsl_matrix_const_column (UtX, t); + gsl_vector_memcpy (&Xgamma_col.vector, &X_col.vector); + + for (size_t ij=0; ij<n_grid; ij++) { + sigma_a2=vec_sa2[ij]; + sigma_b2=vec_sb2[ij]; + + d=CalcMarginal (Xgamma, Uty, K_eval, sigma_a2, sigma_b2, tau); + d-=logm_null[ij]; + d=exp(d); + + gsl_matrix_set(BF, t, ij, d); + } + } + time_Proposal=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0); + + //save results + WriteResult (Hyper, BF); + + //free matrices and vectors + gsl_matrix_free(BF); + gsl_matrix_free(Xgamma); + gsl_matrix_free(Hyper); + return; +} + + + + + +void single_ct_regression(const gsl_matrix_int *Xd, const gsl_vector_int *dlevel, const gsl_vector *pip_vec, gsl_vector *coef, gsl_vector *prior_vec) { + + map<int,double> sum_pip; + map<int,double> sum; + + int levels = gsl_vector_int_get(dlevel,0); + + for(int i=0;i<levels;i++){ + sum_pip[i] = sum[i] = 0; + } + + for(int i=0;i<Xd->size1;i++){ + int cat = gsl_matrix_int_get(Xd,i,0); + sum_pip[cat] += gsl_vector_get(pip_vec,i); + sum[cat] += 1; + } + + for(int i=0;i<Xd->size1;i++){ + int cat = gsl_matrix_int_get(Xd,i,0); + gsl_vector_set(prior_vec,i,sum_pip[cat]/sum[cat]); + } + + //double baseline=0; + for(int i=0;i<levels;i++){ + double new_prior = sum_pip[i]/sum[i]; + //gsl_vector_set(coef, i, log(new_prior/(1-new_prior))-baseline); + //if(i==0){ + //baseline = log(new_prior/(1-new_prior)); + //} + gsl_vector_set(coef, i, log(new_prior/(1-new_prior)) ); + } + + return; +} + + + + +//where A is the ni_test by n_cat matrix of annotations +void BSLMMDAP::DAP_EstimateHyper (const size_t kc, const size_t kd, const vector<string> &vec_rs, const vector<double> &vec_sa2, const vector<double> &vec_sb2, const vector<double> &wab, const vector<vector<vector<double> > > &BF, gsl_matrix *Ac, gsl_matrix_int *Ad, gsl_vector_int *dlevel) { + clock_t time_start; + + //set up BF + double h, rho, sigma_a2, sigma_b2, d, s, logm, logm_save; + size_t t1, t2; + size_t n_grid=wab.size(), ns_test=vec_rs.size(); + + gsl_vector *prior_vec=gsl_vector_alloc(ns_test); + gsl_matrix *Hyper=gsl_matrix_alloc(n_grid, 5); + gsl_vector *pip=gsl_vector_alloc(ns_test); + gsl_vector *coef=gsl_vector_alloc(kc+kd+1); + + //perform the EM algorithm + vector<double> vec_wab, vec_wab_new; + + //initial values + for (size_t t=0; t<ns_test; t++) { + gsl_vector_set (prior_vec, t, (double)BF.size()/(double)ns_test); + } + for (size_t ij=0; ij<n_grid; ij++) { + vec_wab.push_back(wab[ij]); + vec_wab_new.push_back(wab[ij]); + } + + //EM iteration + size_t it=0; + double dif=1; + while (it<100 && dif>1e-3) { + //update E_gamma + t1=0, t2=0; + for (size_t b=0; b<BF.size(); b++) { + s=1; + for (size_t m=0; m<BF[b].size(); m++) { + d=0; + for (size_t ij=0; ij<n_grid; ij++) { + d+=vec_wab_new[ij]*BF[b][m][ij]; + } + d*=gsl_vector_get(prior_vec, t1)/(1-gsl_vector_get(prior_vec, t1)); + + gsl_vector_set(pip, t1, d); + s+=d; + t1++; + } + + for (size_t m=0; m<BF[b].size(); m++) { + d=gsl_vector_get(pip, t2)/s; + gsl_vector_set(pip, t2, d); + t2++; + } + } + + //update E_wab + s=0; + for (size_t ij=0; ij<n_grid; ij++) { + vec_wab_new[ij]=0; + + t1=0; + for (size_t b=0; b<BF.size(); b++) { + d=1; + for (size_t m=0; m<BF[b].size(); m++) { + d+=gsl_vector_get(prior_vec, t1)/(1-gsl_vector_get(prior_vec, t1))*vec_wab[ij]*BF[b][m][ij]; + t1++; + } + vec_wab_new[ij]+=log(d); + } + + s=max(s, vec_wab_new[ij]); + } + + d=0; + for (size_t ij=0; ij<n_grid; ij++) { + vec_wab_new[ij]=exp(vec_wab_new[ij]-s); + d+=vec_wab_new[ij]; + } + + for (size_t ij=0; ij<n_grid; ij++) { + vec_wab_new[ij]/=d; + // vec_wab[ij]=vec_wab_new[ij]; + } + + //update coef, and pi + if(kc==0 && kd==0){//no annotation + s=0; + for (size_t t=0; t<pip->size; t++) { + s+=gsl_vector_get(pip, t); + } + s=s/(double)pip->size; + for (size_t t=0; t<pip->size; t++) { + gsl_vector_set(prior_vec, t, s); + } + + gsl_vector_set (coef, 0, log(s/(1-s))); + } else if(kc==0 && kd!=0){//only discrete annotations + if(kd == 1){ + single_ct_regression(Ad, dlevel, pip, coef, prior_vec); + }else{ + logistic_cat_fit(coef, Ad, dlevel, pip, 0, 0); + logistic_cat_pred(coef, Ad, dlevel, prior_vec); + } + } else if (kc!=0 && kd==0) {//only continuous annotations + logistic_cont_fit(coef, Ac, pip, 0, 0); + logistic_cont_pred(coef, Ac, prior_vec); + } else if (kc!=0 && kd!=0) {//both continuous and categorical annotations + logistic_mixed_fit(coef, Ad, dlevel, Ac, pip, 0, 0); + logistic_mixed_pred(coef, Ad, dlevel, Ac, prior_vec); + } + + //compute marginal likelihood + logm=0; + + t1=0; + for (size_t b=0; b<BF.size(); b++) { + d=1; s=0; + for (size_t m=0; m<BF[b].size(); m++) { + s+=log(1-gsl_vector_get(prior_vec, t1)); + for (size_t ij=0; ij<n_grid; ij++) { + d+=gsl_vector_get(prior_vec, t1)/(1-gsl_vector_get(prior_vec, t1))*vec_wab[ij]*BF[b][m][ij]; + } + } + logm+=log(d)+s; + t1++; + } + + if (it>0) { + dif=logm-logm_save; + } + logm_save=logm; + it++; + + cout<<"iteration = "<<it<<"; marginal likelihood = "<<logm<<endl; + } + + //update h and rho that correspond to w_ab + for (size_t ij=0; ij<n_grid; ij++) { + sigma_a2=vec_sa2[ij]; + sigma_b2=vec_sb2[ij]; + + d=exp(gsl_vector_get(coef, coef->size-1))/(1+exp(gsl_vector_get(coef, coef->size-1))); + h=(d*(double)ns_test*sigma_a2+1*sigma_b2)/(1+d*(double)ns_test*sigma_a2+1*sigma_b2); + rho=d*(double)ns_test*sigma_a2/(d*(double)ns_test*sigma_a2+1*sigma_b2); + + gsl_matrix_set (Hyper, ij, 0, h); + gsl_matrix_set (Hyper, ij, 1, rho); + gsl_matrix_set (Hyper, ij, 2, sigma_a2); + gsl_matrix_set (Hyper, ij, 3, sigma_b2); + gsl_matrix_set (Hyper, ij, 4, vec_wab_new[ij]); + } + + //obtain beta and alpha parameters + + + //save results + WriteResult (vec_rs, Hyper, pip, coef); + + //free matrices and vectors + gsl_vector_free(prior_vec); + gsl_matrix_free(Hyper); + gsl_vector_free(pip); + gsl_vector_free(coef); + return; +} + +/* +//readin the estimated hyper-parameters and perform fine mapping for each region +void BSLMM::DAP_FineMapping (const gsl_matrix *U, const gsl_matrix *UtX, const gsl_matrix *A, const gsl_vector *Uty, const gsl_vector *K_eval, const gsl_vector *y, gsl_matrix *Hyper, gsl_vector *alpha, gsl_vector *pip) { + clock_t time_start; + + //two priority sets: S_1 contains all candidate causal SNPs; S_2 contains the prioritized combintion of them + //two marginal probability sets: P_1 contains marginals for S_1; P_2 contains marginals for S_2; + set<size_t> S1set, S2set; + vector<size_t> S1vec; + vector<set<size_t> > S2vec; + vector<double> P1, P2; + + //calculate P0 (null) and P1 (for every SNP) + + + + //loop through the number of combinations + for (size_t s=0; s<p; s++) { + //if (s==0), set up S_1: compute marginal of the null model, then compute P_1, then compute BF_1 and use them to select S_1; compute C_1 + + + + //if (s==1), set up S_2: compute pair-wise P_2 and use them to select S_2; compute C_2 + + //otherwise, match each combination of S_2 with each SNP from S_1, select into S_3; and replace S_2 with S_3; compute C_s + + + //stop when the stopping critieria are reached (if S_2 is empty; if t; if kappa); add the residual component R + + for (size_t t=0; t<total_step; ++t) { + if (t%d_pace==0 || t==total_step-1) {ProgressBar ("Running MCMC ", t, total_step-1, (double)n_accept/(double)(t*n_mh+1));} +// if (t>10) {break;} + + if (a_mode==13) { + SampleZ (y, z_hat, z); + mean_z=CenterVector (z); + + time_start=clock(); + gsl_blas_dgemv (CblasTrans, 1.0, U, z, 0.0, Utz); + time_UtZ+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0); + + //First proposal + if (cHyp_old.n_gamma==0 || cHyp_old.rho==0) { + logPost_old=CalcPosterior(Utz, K_eval, Utu_old, alpha_old, cHyp_old); + beta_old.clear(); + for (size_t i=0; i<cHyp_old.n_gamma; ++i) { + beta_old.push_back(0); + } + } + else { + gsl_matrix *UtXgamma=gsl_matrix_alloc (ni_test, cHyp_old.n_gamma); + gsl_vector *beta=gsl_vector_alloc (cHyp_old.n_gamma); + SetXgamma (UtXgamma, UtX, rank_old); + logPost_old=CalcPosterior(UtXgamma, Utz, K_eval, UtXb_old, Utu_old, alpha_old, beta, cHyp_old); + + beta_old.clear(); + for (size_t i=0; i<beta->size; ++i) { + beta_old.push_back(gsl_vector_get(beta, i)); + } + gsl_matrix_free (UtXgamma); + gsl_vector_free (beta); + } + } + + + delete [] p_gamma; + beta_g.clear(); + + return; +} + +*/ + + + + + + +/* +//below fits MCMC for rho=1 +void BSLMM::CalcXtX (const gsl_matrix *X, const gsl_vector *y, const size_t s_size, gsl_matrix *XtX, gsl_vector *Xty) +{ + time_t time_start=clock(); + gsl_matrix_const_view X_sub=gsl_matrix_const_submatrix(X, 0, 0, X->size1, s_size); + gsl_matrix_view XtX_sub=gsl_matrix_submatrix(XtX, 0, 0, s_size, s_size); + gsl_vector_view Xty_sub=gsl_vector_subvector(Xty, 0, s_size); + +#ifdef WITH_LAPACK + lapack_dgemm ((char *)"T", (char *)"N", 1.0, &X_sub.matrix, &X_sub.matrix, 0.0, &XtX_sub.matrix); +#else + gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, &X_sub.matrix, &X_sub.matrix, 0.0, &XtX_sub.matrix); +#endif + gsl_blas_dgemv(CblasTrans, 1.0, &X_sub.matrix, y, 0.0, &Xty_sub.vector); + + time_Omega+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0); + + return; +} + + + +double BSLMM::CalcPosterior (const double yty, class HYPBSLMM &cHyp) +{ + double logpost=0.0; + + //for quantitative traits, calculate pve and pge + //pve and pge for case/control data are calculted in CalcCC_PVEnZ + if (a_mode==11) { + cHyp.pve=0.0; + cHyp.pge=1.0; + } + + //calculate likelihood + if (a_mode==11) {logpost-=0.5*(double)ni_test*log(yty);} + else {logpost-=0.5*yty;} + + logpost+=((double)cHyp.n_gamma-1.0)*cHyp.logp+((double)ns_test-(double)cHyp.n_gamma)*log(1-exp(cHyp.logp)); + + return logpost; +} + + +double BSLMM::CalcPosterior (const gsl_matrix *Xgamma, const gsl_matrix *XtX, const gsl_vector *Xty, const double yty, const size_t s_size, gsl_vector *Xb, gsl_vector *beta, class HYPBSLMM &cHyp) +{ + double sigma_a2=cHyp.h/( (1-cHyp.h)*exp(cHyp.logp)*(double)ns_test); + double logpost=0.0; + double d, P_yy=yty, logdet_O=0.0; + + gsl_matrix_const_view Xgamma_sub=gsl_matrix_const_submatrix (Xgamma, 0, 0, Xgamma->size1, s_size); + gsl_matrix_const_view XtX_sub=gsl_matrix_const_submatrix (XtX, 0, 0, s_size, s_size); + gsl_vector_const_view Xty_sub=gsl_vector_const_subvector (Xty, 0, s_size); + + gsl_matrix *Omega=gsl_matrix_alloc (s_size, s_size); + gsl_matrix *M_temp=gsl_matrix_alloc (s_size, s_size); + gsl_vector *beta_hat=gsl_vector_alloc (s_size); + gsl_vector *Xty_temp=gsl_vector_alloc (s_size); + + gsl_vector_memcpy (Xty_temp, &Xty_sub.vector); + + //calculate Omega + gsl_matrix_memcpy (Omega, &XtX_sub.matrix); + gsl_matrix_scale (Omega, sigma_a2); + gsl_matrix_set_identity (M_temp); + gsl_matrix_add (Omega, M_temp); + + //calculate beta_hat + logdet_O=CholeskySolve(Omega, Xty_temp, beta_hat); + gsl_vector_scale (beta_hat, sigma_a2); + + gsl_blas_ddot (Xty_temp, beta_hat, &d); + P_yy-=d; + + //sample tau + double tau=1.0; + if (a_mode==11) {tau =gsl_ran_gamma (gsl_r, (double)ni_test/2.0, 2.0/P_yy); } + + //sample beta + for (size_t i=0; i<s_size; i++) + { + d=gsl_ran_gaussian(gsl_r, 1); + gsl_vector_set(beta, i, d); + } + gsl_vector_view beta_sub=gsl_vector_subvector(beta, 0, s_size); + gsl_blas_dtrsv(CblasUpper, CblasNoTrans, CblasNonUnit, Omega, &beta_sub.vector); + + //it compuates inv(L^T(Omega)) %*% beta; + gsl_vector_scale(&beta_sub.vector, sqrt(sigma_a2/tau)); + gsl_vector_add(&beta_sub.vector, beta_hat); + gsl_blas_dgemv (CblasNoTrans, 1.0, &Xgamma_sub.matrix, &beta_sub.vector, 0.0, Xb); + + //for quantitative traits, calculate pve and pge + if (a_mode==11) { + gsl_blas_ddot (Xb, Xb, &d); + cHyp.pve=d/(double)ni_test; + cHyp.pve/=cHyp.pve+1.0/tau; + cHyp.pge=1.0; + } + + logpost=-0.5*logdet_O; + if (a_mode==11) {logpost-=0.5*(double)ni_test*log(P_yy);} + else {logpost-=0.5*P_yy;} + + logpost+=((double)cHyp.n_gamma-1.0)*cHyp.logp+((double)ns_test-(double)cHyp.n_gamma)*log(1.0-exp(cHyp.logp)); + + gsl_matrix_free (Omega); + gsl_matrix_free (M_temp); + gsl_vector_free (beta_hat); + gsl_vector_free (Xty_temp); + + return logpost; +} +*/ + + diff --git a/src/bslmmdap.h b/src/bslmmdap.h new file mode 100644 index 0000000..ac78f97 --- /dev/null +++ b/src/bslmmdap.h @@ -0,0 +1,101 @@ +/* + Genome-wide Efficient Mixed Model Association (GEMMA) + Copyright (C) 2011 Xiang Zhou + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + + +#ifndef __BSLMMDAP_H__ +#define __BSLMMDAP_H__ + +#include <vector> +#include <map> +#include <gsl/gsl_rng.h> +#include <gsl/gsl_randist.h> + +#ifdef FORCE_FLOAT +#include "param_float.h" +#else +#include "param.h" +#endif + + +using namespace std; + + + + + + +class BSLMMDAP { + +public: + // IO related parameters + int a_mode; + size_t d_pace; + + string file_bfile; + string file_geno; + string file_out; + string path_out; + + // LMM related parameters + double pve_null; + double pheno_mean; + + // BSLMM MCMC related parameters + long int randseed; + double trace_G; + + HYPBSLMM cHyp_initial; + + // Summary statistics + size_t ni_total, ns_total; //number of total individuals and snps + size_t ni_test, ns_test; //number of individuals and snps used for analysis + + double h_min, h_max, rho_min, rho_max; + size_t h_ngrid, rho_ngrid; + + double time_UtZ; + double time_Omega; //time spent on optimization iterations + double time_Proposal; //time spent on constructing the proposal distribution for gamma (i.e. lmm or lm analysis) + vector<int> indicator_idv; //indicator for individuals (phenotypes), 0 missing, 1 available for analysis + vector<int> indicator_snp; //sequence indicator for SNPs: 0 ignored because of (a) maf, (b) miss, (c) non-poly; 1 available for analysis + + vector<SNPINFO> snpInfo; //record SNP information + + // Main Functions + void CopyFromParam (PARAM &cPar); + void CopyToParam (PARAM &cPar); + + void WriteResult (const gsl_matrix *Hyper, const gsl_matrix *BF); + void WriteResult (const vector<string> &vec_rs, const gsl_matrix *Hyper, const gsl_vector *pip, const gsl_vector *coef); + double CalcMarginal (const gsl_vector *Uty, const gsl_vector *K_eval, const double sigma_b2, const double tau); + double CalcMarginal (const gsl_matrix *UtXgamma, const gsl_vector *Uty, const gsl_vector *K_eval, const double sigma_a2, const double sigma_b2, const double tau); + double CalcPrior (class HYPBSLMM &cHyp); + + void DAP_CalcBF (const gsl_matrix *U, const gsl_matrix *UtX, const gsl_vector *Uty, const gsl_vector *K_eval, const gsl_vector *y); + void DAP_EstimateHyper (const size_t kc, const size_t kd, const vector<string> &vec_rs, const vector<double> &vec_sa2, const vector<double> &vec_sb2, const vector<double> &wab, const vector<vector<vector<double> > > &BF, gsl_matrix *Ac, gsl_matrix_int *Ad, gsl_vector_int *dlevel); + +}; + +void ReadFile_hyb (const string &file_hyp, vector<double> &vec_sa2, vector<double> &vec_sb2, vector<double> &vec_wab); +void ReadFile_bf (const string &file_bf, vector<string> &vec_rs, vector<vector<vector<double> > > &BF); +void ReadFile_cat (const string &file_cat, const vector<string> &vec_rs, gsl_matrix *Ac, gsl_matrix_int *Ad, gsl_vector_int *dlevel, size_t &kc, size_t &kd); + + +#endif + + diff --git a/src/ldr.cpp b/src/ldr.cpp new file mode 100644 index 0000000..e28490a --- /dev/null +++ b/src/ldr.cpp @@ -0,0 +1,285 @@ +/* + Genome-wide Efficient Mixed Model Association (GEMMA) + Copyright (C) 2011 Xiang Zhou + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include <iostream> +#include <fstream> +#include <sstream> + +#include <iomanip> +#include <cmath> +#include <iostream> +#include <stdio.h> +#include <stdlib.h> +#include <ctime> +#include <cstring> +#include <algorithm> + +#include "gsl/gsl_vector.h" +#include "gsl/gsl_matrix.h" +#include "gsl/gsl_linalg.h" +#include "gsl/gsl_blas.h" +#include "gsl/gsl_eigen.h" +#include "gsl/gsl_randist.h" +#include "gsl/gsl_cdf.h" +#include "gsl/gsl_roots.h" +#include "Eigen/Dense" + + + +#include "lapack.h" + +#ifdef FORCE_FLOAT +#include "param_float.h" +#include "ldr_float.h" +#include "lm_float.h" +#include "mathfunc_float.h" //for function CenterVector +#else +#include "param.h" +#include "ldr.h" +#include "lm.h" +#include "mathfunc.h" +#endif + +using namespace std; +using namespace Eigen; + + + +void LDR::CopyFromParam (PARAM &cPar) +{ + a_mode=cPar.a_mode; + d_pace=cPar.d_pace; + + file_bfile=cPar.file_bfile; + file_geno=cPar.file_geno; + file_out=cPar.file_out; + path_out=cPar.path_out; + + ni_total=cPar.ni_total; + ns_total=cPar.ns_total; + ni_test=cPar.ni_test; + ns_test=cPar.ns_test; + n_cvt=cPar.n_cvt; + + indicator_idv=cPar.indicator_idv; + indicator_snp=cPar.indicator_snp; + snpInfo=cPar.snpInfo; + + return; +} + + +void LDR::CopyToParam (PARAM &cPar) +{ + //cPar.pheno_mean=pheno_mean; + //cPar.randseed=randseed; + + return; +} + + +/* +void BSLMM::WriteBV (const gsl_vector *bv) +{ + string file_str; + file_str=path_out+"/"+file_out; + file_str+=".bv.txt"; + + ofstream outfile (file_str.c_str(), ofstream::out); + if (!outfile) {cout<<"error writing file: "<<file_str.c_str()<<endl; return;} + + size_t t=0; + for (size_t i=0; i<ni_total; ++i) { + if (indicator_idv[i]==0) { + outfile<<"NA"<<endl; + } + else { + outfile<<scientific<<setprecision(6)<<gsl_vector_get(bv, t)<<endl; + t++; + } + } + + outfile.clear(); + outfile.close(); + return; +} + + + + +void BSLMM::WriteParam (vector<pair<double, double> > &beta_g, const gsl_vector *alpha, const size_t w) +{ + string file_str; + file_str=path_out+"/"+file_out; + file_str+=".param.txt"; + + ofstream outfile (file_str.c_str(), ofstream::out); + if (!outfile) {cout<<"error writing file: "<<file_str.c_str()<<endl; return;} + + outfile<<"chr"<<"\t"<<"rs"<<"\t" + <<"ps"<<"\t"<<"n_miss"<<"\t"<<"alpha"<<"\t" + <<"beta"<<"\t"<<"gamma"<<endl; + + size_t t=0; + for (size_t i=0; i<ns_total; ++i) { + if (indicator_snp[i]==0) {continue;} + + outfile<<snpInfo[i].chr<<"\t"<<snpInfo[i].rs_number<<"\t" + <<snpInfo[i].base_position<<"\t"<<snpInfo[i].n_miss<<"\t"; + + outfile<<scientific<<setprecision(6)<<gsl_vector_get(alpha, t)<<"\t"; + if (beta_g[t].second!=0) { + outfile<<beta_g[t].first/beta_g[t].second<<"\t"<<beta_g[t].second/(double)w<<endl; + } + else { + outfile<<0.0<<"\t"<<0.0<<endl; + } + t++; + } + + outfile.clear(); + outfile.close(); + return; +} + + +void BSLMM::WriteParam (const gsl_vector *alpha) +{ + string file_str; + file_str=path_out+"/"+file_out; + file_str+=".param.txt"; + + ofstream outfile (file_str.c_str(), ofstream::out); + if (!outfile) {cout<<"error writing file: "<<file_str.c_str()<<endl; return;} + + outfile<<"chr"<<"\t"<<"rs"<<"\t" + <<"ps"<<"\t"<<"n_miss"<<"\t"<<"alpha"<<"\t" + <<"beta"<<"\t"<<"gamma"<<endl; + + size_t t=0; + for (size_t i=0; i<ns_total; ++i) { + if (indicator_snp[i]==0) {continue;} + + outfile<<snpInfo[i].chr<<"\t"<<snpInfo[i].rs_number<<"\t" + <<snpInfo[i].base_position<<"\t"<<snpInfo[i].n_miss<<"\t"; + outfile<<scientific<<setprecision(6)<<gsl_vector_get(alpha, t)<<"\t"; + outfile<<0.0<<"\t"<<0.0<<endl; + t++; + } + + outfile.clear(); + outfile.close(); + return; +} + + +void BSLMM::WriteResult (const int flag, const gsl_matrix *Result_hyp, const gsl_matrix *Result_gamma, const size_t w_col) +{ + string file_gamma, file_hyp; + file_gamma=path_out+"/"+file_out; + file_gamma+=".gamma.txt"; + file_hyp=path_out+"/"+file_out; + file_hyp+=".hyp.txt"; + + ofstream outfile_gamma, outfile_hyp; + + if (flag==0) { + outfile_gamma.open (file_gamma.c_str(), ofstream::out); + outfile_hyp.open (file_hyp.c_str(), ofstream::out); + if (!outfile_gamma) {cout<<"error writing file: "<<file_gamma<<endl; return;} + if (!outfile_hyp) {cout<<"error writing file: "<<file_hyp<<endl; return;} + + outfile_hyp<<"h \t pve \t rho \t pge \t pi \t n_gamma"<<endl; + + for (size_t i=0; i<s_max; ++i) { + outfile_gamma<<"s"<<i<<"\t"; + } + outfile_gamma<<endl; + } + else { + outfile_gamma.open (file_gamma.c_str(), ofstream::app); + outfile_hyp.open (file_hyp.c_str(), ofstream::app); + if (!outfile_gamma) {cout<<"error writing file: "<<file_gamma<<endl; return;} + if (!outfile_hyp) {cout<<"error writing file: "<<file_hyp<<endl; return;} + + size_t w; + if (w_col==0) {w=w_pace;} + else {w=w_col;} + + for (size_t i=0; i<w; ++i) { + outfile_hyp<<scientific; + for (size_t j=0; j<4; ++j) { + outfile_hyp<<setprecision(6)<<gsl_matrix_get (Result_hyp, i, j)<<"\t"; + } + outfile_hyp<<setprecision(6)<<exp(gsl_matrix_get (Result_hyp, i, 4))<<"\t"; + outfile_hyp<<(int)gsl_matrix_get (Result_hyp, i, 5)<<"\t"; + outfile_hyp<<endl; + } + + for (size_t i=0; i<w; ++i) { + for (size_t j=0; j<s_max; ++j) { + outfile_gamma<<(int)gsl_matrix_get (Result_gamma, i, j)<<"\t"; + } + outfile_gamma<<endl; + } + + } + + outfile_hyp.close(); + outfile_hyp.clear(); + outfile_gamma.close(); + outfile_gamma.clear(); + return; +} + +*/ + +//X is a p by n matrix +void LDR::VB (const vector<vector<unsigned char> > &Xt, const gsl_matrix *W_gsl, const gsl_vector *y_gsl) +{ + //save gsl_vector and gsl_matrix into eigen library formats + MatrixXd W(W_gsl->size1, W_gsl->size2); + VectorXd y(y_gsl->size); + VectorXd x_col(y_gsl->size); + + double d; + for (size_t i=0; i<W_gsl->size1; i++) { + d=gsl_vector_get(y_gsl, i); + y(i)=d; + for (size_t j=0; j<W_gsl->size2; j++) { + W(i,j)=gsl_matrix_get(W_gsl, i, j); + } + } + + //initial VB values by lm + cout<<indicator_snp[0]<<" "<<indicator_snp[1]<<" "<<indicator_snp[2]<<endl; + uchar_matrix_get_row (Xt, 0, x_col); + + for (size_t j=0; j<10; j++) { + cout<<x_col(j)<<endl; + } + + + //run VB iterations + + + + //save results + + return; +} diff --git a/src/ldr.h b/src/ldr.h new file mode 100644 index 0000000..d81048c --- /dev/null +++ b/src/ldr.h @@ -0,0 +1,77 @@ +/* + Genome-wide Efficient Mixed Model Association (GEMMA) + Copyright (C) 2011 Xiang Zhou + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + + +#ifndef __LDR_H__ +#define __LDR_H__ + +#include <vector> +#include <map> +#include <gsl/gsl_rng.h> +#include <gsl/gsl_randist.h> + +#ifdef FORCE_FLOAT +#include "param_float.h" +#else +#include "param.h" +#endif + + +using namespace std; + + + + + + +class LDR { + +public: + // IO related parameters + int a_mode; + size_t d_pace; + + string file_bfile; + string file_geno; + string file_out; + string path_out; + + // Summary statistics + size_t ni_total, ns_total; //number of total individuals and snps + size_t ni_test, ns_test; //number of individuals and snps used for analysis + size_t n_cvt; //number of covariates + vector<int> indicator_idv; //indicator for individuals (phenotypes), 0 missing, 1 available for analysis + vector<int> indicator_snp; //sequence indicator for SNPs: 0 ignored because of (a) maf, (b) miss, (c) non-poly; 1 available for analysis + + vector<SNPINFO> snpInfo; //record SNP information + + // Not included in PARAM + gsl_rng *gsl_r; + + // Main Functions + void CopyFromParam (PARAM &cPar); + void CopyToParam (PARAM &cPar); + + void VB(const vector<vector<unsigned char> > &Xt, const gsl_matrix *W_gsl, const gsl_vector *y_gsl); +}; + + + +#endif + + diff --git a/src/logistic.cpp b/src/logistic.cpp new file mode 100644 index 0000000..1b47946 --- /dev/null +++ b/src/logistic.cpp @@ -0,0 +1,783 @@ +#include <stdio.h> +#include <math.h> +#include <gsl/gsl_matrix.h> +#include <gsl/gsl_rng.h> +#include <gsl/gsl_multimin.h> +#include <gsl/gsl_sf.h> +#include <gsl/gsl_linalg.h> +#include "logistic.h" + +// I need to bundle all the data that goes to the function to optimze together. +typedef struct{ + gsl_matrix_int *X; + gsl_vector_int *nlev; + gsl_vector *y; + gsl_matrix *Xc; // continuous covariates Matrix Nobs x Kc (NULL if not used) + double lambdaL1; + double lambdaL2; +}fix_parm_mixed_T; + + + + + + +double fLogit_mixed(gsl_vector *beta + ,gsl_matrix_int *X + ,gsl_vector_int *nlev + ,gsl_matrix *Xc + ,gsl_vector *y + ,double lambdaL1 + ,double lambdaL2) +{ + int n = y->size; + // int k = X->size2; + int npar = beta->size; + double total = 0; + double aux = 0; + /* omp_set_num_threads(ompthr); */ + /* /\* Changed loop start at 1 instead of 0 to avoid regularization of beta_0*\/ */ + /* /\*#pragma omp parallel for reduction (+:total)*\/ */ + for(int i = 1; i < npar; ++i) + total += beta->data[i]*beta->data[i]; + total = (-total*lambdaL2/2); + /* /\*#pragma omp parallel for reduction (+:aux)*\/ */ + for(int i = 1; i < npar; ++i) + aux += (beta->data[i]>0 ? beta->data[i] : -beta->data[i]); + total = total-aux*lambdaL1; + /* #pragma omp parallel for schedule(static) shared(n,beta,X,nlev,y) reduction (+:total) */ + for(int i = 0; i < n; ++i) { + double Xbetai=beta->data[0]; + int iParm=1; + for(int k = 0; k < X->size2; ++k) { + if(gsl_matrix_int_get(X,i,k)>0) + Xbetai+=beta->data[gsl_matrix_int_get(X,i,k)-1+iParm]; + iParm+=nlev->data[k]-1; + } + for(int k = 0; k < (Xc->size2); ++k) + Xbetai+= gsl_matrix_get(Xc,i,k)*beta->data[iParm++]; + total += y->data[i]*Xbetai-gsl_sf_log_1plusx(gsl_sf_exp(Xbetai)); + } + return -total; +} + + +void logistic_mixed_pred(gsl_vector *beta // Vector of parameters length = 1 + Sum_k(C_k - 1) + ,gsl_matrix_int *X //Matrix Nobs x K + ,gsl_vector_int *nlev // Vector with number categories + ,gsl_matrix *Xc // continuous covariates Matrix Nobs x Kc (NULL if not used) + ,gsl_vector *yhat //Vector of prob. predicted by the logistic + ) +{ + for(int i = 0; i < X->size1; ++i) { + double Xbetai=beta->data[0]; + int iParm=1; + for(int k = 0; k < X->size2; ++k) { + if(gsl_matrix_int_get(X,i,k)>0) + Xbetai+=beta->data[gsl_matrix_int_get(X,i,k)-1+iParm]; + iParm+=nlev->data[k]-1; + } + // Adding the continuous + for(int k = 0; k < (Xc->size2); ++k) + Xbetai+= gsl_matrix_get(Xc,i,k)*beta->data[iParm++]; + yhat->data[i]=1/(1 + gsl_sf_exp(-Xbetai)); + } +} + + +/* The gradient of f, df = (df/dx, df/dy). */ +void +wgsl_mixed_optim_df (const gsl_vector *beta, void *params, + gsl_vector *out) +{ + fix_parm_mixed_T *p = (fix_parm_mixed_T *)params; + int n = p->y->size; + int K = p->X->size2; + int Kc = p->Xc->size2; + int npar = beta->size; + // Intitialize gradient out necessary? + for(int i = 0; i < npar; ++i) + out->data[i]= 0; + /* Changed loop start at 1 instead of 0 to avoid regularization of beta 0 */ + for(int i = 1; i < npar; ++i) + out->data[i]= p->lambdaL2*beta->data[i]; + for(int i = 1; i < npar; ++i) + out->data[i]+= p->lambdaL1*((beta->data[i]>0)-(beta->data[i]<0)); + + for(int i = 0; i < n; ++i) { + double pn=0; + double Xbetai=beta->data[0]; + int iParm=1; + for(int k = 0; k < K; ++k) { + if(gsl_matrix_int_get(p->X,i,k)>0) + Xbetai+=beta->data[gsl_matrix_int_get(p->X,i,k)-1+iParm]; + iParm+=p->nlev->data[k]-1; + } + // Adding the continuous + for(int k = 0; k < Kc; ++k) + Xbetai+= gsl_matrix_get(p->Xc,i,k)*beta->data[iParm++]; + + pn= -( p->y->data[i] - 1/(1 + gsl_sf_exp(-Xbetai)) ); + + out->data[0]+= pn; + iParm=1; + for(int k = 0; k < K; ++k) { + if(gsl_matrix_int_get(p->X,i,k)>0) + out->data[gsl_matrix_int_get(p->X,i,k)-1+iParm]+=pn; + iParm+=p->nlev->data[k]-1; + } + // Adding the continuous + for(int k = 0; k < Kc; ++k) { + out->data[iParm++] += gsl_matrix_get(p->Xc,i,k)*pn; + } + } + +} + + +/* The Hessian of f */ +void +wgsl_mixed_optim_hessian (const gsl_vector *beta, void *params, + gsl_matrix *out) +{ + fix_parm_mixed_T *p = (fix_parm_mixed_T *)params; + int n = p->y->size; + int K = p->X->size2; + int Kc = p->Xc->size2; + int npar = beta->size; + gsl_vector *gn = gsl_vector_alloc(npar); /* gn */ + // Intitialize Hessian out necessary ??? + gsl_matrix_set_zero(out); + /* Changed loop start at 1 instead of 0 to avoid regularization of beta 0*/ + for(int i = 1; i < npar; ++i) + gsl_matrix_set(out,i,i,(p->lambdaL2)); // Double check this + // L1 penalty not working yet, as not differentiable, I may need to do coordinate descent (as in glm_net) + + for(int i = 0; i < n; ++i) { + double pn=0; + double aux=0; + double Xbetai=beta->data[0]; + int iParm1=1; + for(int k = 0; k < K; ++k) { + if(gsl_matrix_int_get(p->X,i,k)>0) + Xbetai+=beta->data[gsl_matrix_int_get(p->X,i,k)-1+iParm1]; + iParm1+=p->nlev->data[k]-1; //-1? + } + // Adding the continuous + for(int k = 0; k < Kc; ++k) + Xbetai+= gsl_matrix_get(p->Xc,i,k)*beta->data[iParm1++]; + + pn= 1/(1 + gsl_sf_exp(-Xbetai)); + // Add a protection for pn very close to 0 or 1? + aux=pn*(1-pn); + + // Calculate sub-gradient vector gn + gsl_vector_set_zero(gn); + gn->data[0]= 1; + iParm1=1; + for(int k = 0; k < K; ++k) { + if(gsl_matrix_int_get(p->X,i,k)>0) + gn->data[gsl_matrix_int_get(p->X,i,k)-1+iParm1]=1; + iParm1+=p->nlev->data[k]-1; + } + // Adding the continuous + for(int k = 0; k < Kc; ++k) { + gn->data[iParm1++] = gsl_matrix_get(p->Xc,i,k); + } + + for(int k1=0;k1<npar; ++k1) + if(gn->data[k1]!=0) + for(int k2=0;k2<npar; ++k2) + if(gn->data[k2]!=0) + *gsl_matrix_ptr(out,k1,k2) += (aux * gn->data[k1] * gn->data[k2]); + } + gsl_vector_free(gn); +} + + +double wgsl_mixed_optim_f(gsl_vector *v, void *params) +{ + double mLogLik=0; + fix_parm_mixed_T *p = (fix_parm_mixed_T *)params; + mLogLik = fLogit_mixed(v,p->X,p->nlev,p->Xc,p->y,p->lambdaL1,p->lambdaL2); + return mLogLik; +} + + +/* Compute both f and df together. */ +void +wgsl_mixed_optim_fdf (gsl_vector *x, void *params, + double *f, gsl_vector *df) +{ + *f = wgsl_mixed_optim_f(x, params); + wgsl_mixed_optim_df(x, params, df); +} + + +int logistic_mixed_fit(gsl_vector *beta + ,gsl_matrix_int *X + ,gsl_vector_int *nlev + ,gsl_matrix *Xc // continuous covariates Matrix Nobs x Kc (NULL if not used) + ,gsl_vector *y + ,double lambdaL1 + ,double lambdaL2) +{ + + double mLogLik=0; + fix_parm_mixed_T p; + int npar = beta->size; + int iter=0; + double maxchange=0; + + //Intializing fix parameters + p.X=X; + p.Xc=Xc; + p.nlev=nlev; + p.y=y; + p.lambdaL1=lambdaL1; + p.lambdaL2=lambdaL2; + + //Initial fit + //#ifdef _RPR_DEBUG_ + mLogLik = wgsl_mixed_optim_f(beta,&p); + //fprintf(stderr,"#Initial -log(Lik(0))=%lf\n",mLogLik); + //#endif //_RPR_DEBUG + + gsl_matrix *myH = gsl_matrix_alloc(npar,npar); /* Hessian matrix*/ + gsl_vector *stBeta = gsl_vector_alloc(npar); /* Direction to move */ + + gsl_vector *myG = gsl_vector_alloc(npar); /* Gradient*/ + gsl_vector *tau = gsl_vector_alloc(npar); /* tau for QR*/ + + for(iter=0;iter<100;iter++){ + wgsl_mixed_optim_hessian(beta,&p,myH); //Calculate Hessian + wgsl_mixed_optim_df(beta,&p,myG); //Calculate Gradient + gsl_linalg_QR_decomp(myH,tau); //Calculate next beta + gsl_linalg_QR_solve(myH,tau,myG,stBeta); + gsl_vector_sub(beta,stBeta); + + //Monitor convergence + maxchange=0; + for(int i=0;i<npar; i++) + if(maxchange<fabs(stBeta->data[i])) + maxchange=fabs(stBeta->data[i]); + +#ifdef _RPR_DEBUG_ + //mLogLik = wgsl_mixed_optim_f(beta,&p); + //fprintf(stderr,"#iter %d, -log(Lik(0))=%lf,%lf\n",(int)iter,mLogLik,maxchange); +#endif //_RPR_DEBUG + + if(maxchange<1E-4) + break; + } + +#ifdef _RPR_DEBUG_ + //for (int i = 0; i < npar; i++) + // fprintf(stderr,"#par_%d= %lf\n",i,beta->data[i]); +#endif //_RPR_DEBUG + + //Final fit + mLogLik = wgsl_mixed_optim_f(beta,&p); + //fprintf(stderr,"#Final %d) -log(Lik(0))=%lf, maxchange %lf\n",iter,mLogLik,maxchange); + + gsl_vector_free (tau); + gsl_vector_free (stBeta); + gsl_vector_free (myG); + gsl_matrix_free (myH); + + return 0; +} + +/***************/ +/* Categorical */ +/***************/ + +// I need to bundle all the data that goes to the function to optimze together. +typedef struct{ + gsl_matrix_int *X; + gsl_vector_int *nlev; + gsl_vector *y; + double lambdaL1; + double lambdaL2; +}fix_parm_cat_T; + + +double fLogit_cat(gsl_vector *beta + ,gsl_matrix_int *X + ,gsl_vector_int *nlev + ,gsl_vector *y + ,double lambdaL1 + ,double lambdaL2) +{ + int n = y->size; + // int k = X->size2; + int npar = beta->size; + double total = 0; + double aux = 0; + /* omp_set_num_threads(ompthr); */ + /* /\* Changed loop start at 1 instead of 0 to avoid regularization of beta 0*\/ */ + /* /\*#pragma omp parallel for reduction (+:total)*\/ */ + for(int i = 1; i < npar; ++i) + total += beta->data[i]*beta->data[i]; + total = (-total*lambdaL2/2); + /* /\*#pragma omp parallel for reduction (+:aux)*\/ */ + for(int i = 1; i < npar; ++i) + aux += (beta->data[i]>0 ? beta->data[i] : -beta->data[i]); + total = total-aux*lambdaL1; + /* #pragma omp parallel for schedule(static) shared(n,beta,X,nlev,y) reduction (+:total) */ + for(int i = 0; i < n; ++i) { + double Xbetai=beta->data[0]; + int iParm=1; + for(int k = 0; k < X->size2; ++k) { + if(gsl_matrix_int_get(X,i,k)>0) + Xbetai+=beta->data[gsl_matrix_int_get(X,i,k)-1+iParm]; + iParm+=nlev->data[k]-1; + } + total += y->data[i]*Xbetai-gsl_sf_log_1plusx(gsl_sf_exp(Xbetai)); + } + return -total; +} + + +void logistic_cat_pred(gsl_vector *beta // Vector of parameters length = 1 + Sum_k(C_k - 1) + ,gsl_matrix_int *X //Matrix Nobs x K + ,gsl_vector_int *nlev // Vector with number categories + ,gsl_vector *yhat //Vector of prob. predicted by the logistic + ) +{ + for(int i = 0; i < X->size1; ++i) { + double Xbetai=beta->data[0]; + int iParm=1; + for(int k = 0; k < X->size2; ++k) { + if(gsl_matrix_int_get(X,i,k)>0) + Xbetai+=beta->data[gsl_matrix_int_get(X,i,k)-1+iParm]; + iParm+=nlev->data[k]-1; + } + yhat->data[i]=1/(1 + gsl_sf_exp(-Xbetai)); + } +} + + +/* The gradient of f, df = (df/dx, df/dy). */ +void +wgsl_cat_optim_df (const gsl_vector *beta, void *params, + gsl_vector *out) +{ + fix_parm_cat_T *p = (fix_parm_cat_T *)params; + int n = p->y->size; + int K = p->X->size2; + int npar = beta->size; + // Intitialize gradient out necessary? + for(int i = 0; i < npar; ++i) + out->data[i]= 0; + /* Changed loop start at 1 instead of 0 to avoid regularization of beta 0 */ + for(int i = 1; i < npar; ++i) + out->data[i]= p->lambdaL2*beta->data[i]; + for(int i = 1; i < npar; ++i) + out->data[i]+= p->lambdaL1*((beta->data[i]>0)-(beta->data[i]<0)); + + for(int i = 0; i < n; ++i) { + double pn=0; + double Xbetai=beta->data[0]; + int iParm=1; + for(int k = 0; k < K; ++k) { + if(gsl_matrix_int_get(p->X,i,k)>0) + Xbetai+=beta->data[gsl_matrix_int_get(p->X,i,k)-1+iParm]; + iParm+=p->nlev->data[k]-1; + } + // total += y->data[i]*Xbetai-log(1+gsl_sf_exp(Xbetai)); + pn= -( p->y->data[i] - 1/(1 + gsl_sf_exp(-Xbetai)) ); + + out->data[0]+= pn; + iParm=1; + for(int k = 0; k < K; ++k) { + if(gsl_matrix_int_get(p->X,i,k)>0) + out->data[gsl_matrix_int_get(p->X,i,k)-1+iParm]+=pn; + iParm+=p->nlev->data[k]-1; + } + } +} + + +/* The Hessian of f */ +void +wgsl_cat_optim_hessian (const gsl_vector *beta, void *params, + gsl_matrix *out) +{ + fix_parm_cat_T *p = (fix_parm_cat_T *)params; + int n = p->y->size; + int K = p->X->size2; + int npar = beta->size; + // Intitialize Hessian out necessary ??? + gsl_matrix_set_zero(out); + /* Changed loop start at 1 instead of 0 to avoid regularization of beta 0*/ + for(int i = 1; i < npar; ++i) + gsl_matrix_set(out,i,i,(p->lambdaL2)); // Double check this + // L1 penalty not working yet, as not differentiable, I may need to do coordinate descent (as in glm_net) + + + for(int i = 0; i < n; ++i) { + double pn=0; + double aux=0; + double Xbetai=beta->data[0]; + int iParm2=1; + int iParm1=1; + for(int k = 0; k < K; ++k) { + if(gsl_matrix_int_get(p->X,i,k)>0) + Xbetai+=beta->data[gsl_matrix_int_get(p->X,i,k)-1+iParm1]; + iParm1+=p->nlev->data[k]-1; //-1? + } + // total += y->data[i]*Xbetai-log(1+gsl_sf_exp(Xbetai)); + pn= 1/(1 + gsl_sf_exp(-Xbetai)); + // Add a protection for pn very close to 0 or 1? + aux=pn*(1-pn); + *gsl_matrix_ptr(out,0,0)+=aux; + iParm2=1; + for(int k2 = 0; k2 < K; ++k2) { + if(gsl_matrix_int_get(p->X,i,k2)>0) + *gsl_matrix_ptr(out,0,gsl_matrix_int_get(p->X,i,k2)-1+iParm2)+=aux; + iParm2+=p->nlev->data[k2]-1; //-1? + } + iParm1=1; + for(int k1 = 0; k1 < K; ++k1) { + if(gsl_matrix_int_get(p->X,i,k1)>0) + *gsl_matrix_ptr(out,gsl_matrix_int_get(p->X,i,k1)-1+iParm1,0)+=aux; + iParm2=1; + for(int k2 = 0; k2 < K; ++k2) { + if((gsl_matrix_int_get(p->X,i,k1)>0) && (gsl_matrix_int_get(p->X,i,k2)>0)) + *gsl_matrix_ptr(out + ,gsl_matrix_int_get(p->X,i,k1)-1+iParm1 + ,gsl_matrix_int_get(p->X,i,k2)-1+iParm2 + )+=aux; + iParm2+=p->nlev->data[k2]-1; //-1? + } + iParm1+=p->nlev->data[k1]-1; //-1? + } + } +} + + +double wgsl_cat_optim_f(gsl_vector *v, void *params) +{ + double mLogLik=0; + fix_parm_cat_T *p = (fix_parm_cat_T *)params; + mLogLik = fLogit_cat(v,p->X,p->nlev,p->y,p->lambdaL1,p->lambdaL2); + return mLogLik; +} + + +/* Compute both f and df together. */ +void +wgsl_cat_optim_fdf (gsl_vector *x, void *params, + double *f, gsl_vector *df) +{ + *f = wgsl_cat_optim_f(x, params); + wgsl_cat_optim_df(x, params, df); +} + + +int logistic_cat_fit(gsl_vector *beta + ,gsl_matrix_int *X + ,gsl_vector_int *nlev + ,gsl_vector *y + ,double lambdaL1 + ,double lambdaL2) +{ + double mLogLik=0; + fix_parm_cat_T p; + int npar = beta->size; + int iter=0; + double maxchange=0; + + //Intializing fix parameters + p.X=X; + p.nlev=nlev; + p.y=y; + p.lambdaL1=lambdaL1; + p.lambdaL2=lambdaL2; + + //Initial fit + //#ifdef _RPR_DEBUG_ + mLogLik = wgsl_cat_optim_f(beta,&p); + //fprintf(stderr,"#Initial -log(Lik(0))=%lf\n",mLogLik); + //#endif //_RPR_DEBUG + + gsl_matrix *myH = gsl_matrix_alloc(npar,npar); /* Hessian matrix*/ + gsl_vector *stBeta = gsl_vector_alloc(npar); /* Direction to move */ + + gsl_vector *myG = gsl_vector_alloc(npar); /* Gradient*/ + gsl_vector *tau = gsl_vector_alloc(npar); /* tau for QR*/ + + for(iter=0;iter<100;iter++){ + wgsl_cat_optim_hessian(beta,&p,myH); //Calculate Hessian + wgsl_cat_optim_df(beta,&p,myG); //Calculate Gradient + gsl_linalg_QR_decomp(myH,tau); //Calculate next beta + gsl_linalg_QR_solve(myH,tau,myG,stBeta); + gsl_vector_sub(beta,stBeta); + + //Monitor convergence + maxchange=0; + for(int i=0;i<npar; i++) + if(maxchange<fabs(stBeta->data[i])) + maxchange=fabs(stBeta->data[i]); + +#ifdef _RPR_DEBUG_ + mLogLik = wgsl_cat_optim_f(beta,&p); + //fprintf(stderr,"#iter %d, -log(Lik(0))=%lf,%lf\n",(int)iter,mLogLik,maxchange); +#endif //_RPR_DEBUG + + if(maxchange<1E-4) + break; + } + +#ifdef _RPR_DEBUG_ + //for (int i = 0; i < npar; i++) + // fprintf(stderr,"#par_%d= %lf\n",i,beta->data[i]); +#endif //_RPR_DEBUG + + //Final fit + mLogLik = wgsl_cat_optim_f(beta,&p); + //fprintf(stderr,"#Final %d) -log(Lik(0))=%lf, maxchange %lf\n",iter,mLogLik,maxchange); + + gsl_vector_free (tau); + gsl_vector_free (stBeta); + gsl_vector_free (myG); + gsl_matrix_free (myH); + + return 0; +} + + +/***************/ +/* Continuous */ +/***************/ + +// I need to bundle all the data that goes to the function to optimze together. +typedef struct{ + gsl_matrix *Xc; // continuous covariates Matrix Nobs x Kc + gsl_vector *y; + double lambdaL1; + double lambdaL2; +}fix_parm_cont_T; + + +double fLogit_cont(gsl_vector *beta + ,gsl_matrix *Xc + ,gsl_vector *y + ,double lambdaL1 + ,double lambdaL2) +{ + int n = y->size; + int npar = beta->size; + double total = 0; + double aux = 0; + /* omp_set_num_threads(ompthr); */ + /* /\* Changed loop start at 1 instead of 0 to avoid regularization of beta_0*\/ */ + /* /\*#pragma omp parallel for reduction (+:total)*\/ */ + for(int i = 1; i < npar; ++i) + total += beta->data[i]*beta->data[i]; + total = (-total*lambdaL2/2); + /* /\*#pragma omp parallel for reduction (+:aux)*\/ */ + for(int i = 1; i < npar; ++i) + aux += (beta->data[i]>0 ? beta->data[i] : -beta->data[i]); + total = total-aux*lambdaL1; + /* #pragma omp parallel for schedule(static) shared(n,beta,X,nlev,y) reduction (+:total) */ + for(int i = 0; i < n; ++i) { + double Xbetai=beta->data[0]; + int iParm=1; + for(int k = 0; k < (Xc->size2); ++k) + Xbetai+= gsl_matrix_get(Xc,i,k)*beta->data[iParm++]; + total += y->data[i]*Xbetai-gsl_sf_log_1plusx(gsl_sf_exp(Xbetai)); + } + return -total; +} + + +void logistic_cont_pred(gsl_vector *beta // Vector of parameters length = 1 + Sum_k(C_k - 1) + ,gsl_matrix *Xc // continuous covariates Matrix Nobs x Kc (NULL if not used) + ,gsl_vector *yhat //Vector of prob. predicted by the logistic + ) +{ + for(int i = 0; i < Xc->size1; ++i) { + double Xbetai=beta->data[0]; + int iParm=1; + for(int k = 0; k < (Xc->size2); ++k) + Xbetai+= gsl_matrix_get(Xc,i,k)*beta->data[iParm++]; + yhat->data[i]=1/(1 + gsl_sf_exp(-Xbetai)); + } +} + + +/* The gradient of f, df = (df/dx, df/dy). */ +void +wgsl_cont_optim_df (const gsl_vector *beta, void *params, + gsl_vector *out) +{ + fix_parm_cont_T *p = (fix_parm_cont_T *)params; + int n = p->y->size; + int Kc = p->Xc->size2; + int npar = beta->size; + // Intitialize gradient out necessary? + for(int i = 0; i < npar; ++i) + out->data[i]= 0; + /* Changed loop start at 1 instead of 0 to avoid regularization of beta 0 */ + for(int i = 1; i < npar; ++i) + out->data[i]= p->lambdaL2*beta->data[i]; + for(int i = 1; i < npar; ++i) + out->data[i]+= p->lambdaL1*((beta->data[i]>0)-(beta->data[i]<0)); + + for(int i = 0; i < n; ++i) { + double pn=0; + double Xbetai=beta->data[0]; + int iParm=1; + for(int k = 0; k < Kc; ++k) + Xbetai+= gsl_matrix_get(p->Xc,i,k)*beta->data[iParm++]; + + pn= -( p->y->data[i] - 1/(1 + gsl_sf_exp(-Xbetai)) ); + + out->data[0]+= pn; + iParm=1; + // Adding the continuous + for(int k = 0; k < Kc; ++k) { + out->data[iParm++] += gsl_matrix_get(p->Xc,i,k)*pn; + } + } +} + + +/* The Hessian of f */ +void +wgsl_cont_optim_hessian (const gsl_vector *beta, void *params, + gsl_matrix *out) +{ + fix_parm_cont_T *p = (fix_parm_cont_T *)params; + int n = p->y->size; + int Kc = p->Xc->size2; + int npar = beta->size; + gsl_vector *gn = gsl_vector_alloc(npar); /* gn */ + // Intitialize Hessian out necessary ??? + gsl_matrix_set_zero(out); + /* Changed loop start at 1 instead of 0 to avoid regularization of beta 0*/ + for(int i = 1; i < npar; ++i) + gsl_matrix_set(out,i,i,(p->lambdaL2)); // Double check this + // L1 penalty not working yet, as not differentiable, I may need to do coordinate descent (as in glm_net) + + for(int i = 0; i < n; ++i) { + double pn=0; + double aux=0; + double Xbetai=beta->data[0]; + int iParm1=1; + for(int k = 0; k < Kc; ++k) + Xbetai+= gsl_matrix_get(p->Xc,i,k)*beta->data[iParm1++]; + + pn= 1/(1 + gsl_sf_exp(-Xbetai)); + // Add a protection for pn very close to 0 or 1? + aux=pn*(1-pn); + + // Calculate sub-gradient vector gn + gsl_vector_set_zero(gn); + gn->data[0]= 1; + iParm1=1; + for(int k = 0; k < Kc; ++k) { + gn->data[iParm1++] = gsl_matrix_get(p->Xc,i,k); + } + + for(int k1=0;k1<npar; ++k1) + if(gn->data[k1]!=0) + for(int k2=0;k2<npar; ++k2) + if(gn->data[k2]!=0) + *gsl_matrix_ptr(out,k1,k2) += (aux * gn->data[k1] * gn->data[k2]); + } + gsl_vector_free(gn); +} + + +double wgsl_cont_optim_f(gsl_vector *v, void *params) +{ + double mLogLik=0; + fix_parm_cont_T *p = (fix_parm_cont_T *)params; + mLogLik = fLogit_cont(v,p->Xc,p->y,p->lambdaL1,p->lambdaL2); + return mLogLik; +} + + +/* Compute both f and df together. */ +void +wgsl_cont_optim_fdf (gsl_vector *x, void *params, + double *f, gsl_vector *df) +{ + *f = wgsl_cont_optim_f(x, params); + wgsl_cont_optim_df(x, params, df); +} + + +int logistic_cont_fit(gsl_vector *beta + ,gsl_matrix *Xc // continuous covariates Matrix Nobs x Kc (NULL if not used) + ,gsl_vector *y + ,double lambdaL1 + ,double lambdaL2) +{ + + double mLogLik=0; + fix_parm_cont_T p; + int npar = beta->size; + int iter=0; + double maxchange=0; + + //Intializing fix parameters + p.Xc=Xc; + p.y=y; + p.lambdaL1=lambdaL1; + p.lambdaL2=lambdaL2; + + //Initial fit + //#ifdef _RPR_DEBUG_ + mLogLik = wgsl_cont_optim_f(beta,&p); + //fprintf(stderr,"#Initial -log(Lik(0))=%lf\n",mLogLik); + //#endif //_RPR_DEBUG + + gsl_matrix *myH = gsl_matrix_alloc(npar,npar); /* Hessian matrix*/ + gsl_vector *stBeta = gsl_vector_alloc(npar); /* Direction to move */ + + gsl_vector *myG = gsl_vector_alloc(npar); /* Gradient*/ + gsl_vector *tau = gsl_vector_alloc(npar); /* tau for QR*/ + + for(iter=0;iter<100;iter++){ + wgsl_cont_optim_hessian(beta,&p,myH); //Calculate Hessian + wgsl_cont_optim_df(beta,&p,myG); //Calculate Gradient + gsl_linalg_QR_decomp(myH,tau); //Calculate next beta + gsl_linalg_QR_solve(myH,tau,myG,stBeta); + gsl_vector_sub(beta,stBeta); + + //Monitor convergence + maxchange=0; + for(int i=0;i<npar; i++) + if(maxchange<fabs(stBeta->data[i])) + maxchange=fabs(stBeta->data[i]); + +#ifdef _RPR_DEBUG_ + mLogLik = wgsl_cont_optim_f(beta,&p); + //fprintf(stderr,"#iter %d, -log(Lik(0))=%lf,%lf\n",(int)iter,mLogLik,maxchange); +#endif //_RPR_DEBUG + + if(maxchange<1E-4) + break; + } + +#ifdef _RPR_DEBUG_ + //for (int i = 0; i < npar; i++) + // fprintf(stderr,"#par_%d= %lf\n",i,beta->data[i]); +#endif //_RPR_DEBUG + + //Final fit + mLogLik = wgsl_cont_optim_f(beta,&p); + //fprintf(stderr,"#Final %d) -log(Lik(0))=%lf, maxchange %lf\n",iter,mLogLik,maxchange); + + gsl_vector_free (tau); + gsl_vector_free (stBeta); + gsl_vector_free (myG); + gsl_matrix_free (myH); + + return 0; +} + diff --git a/src/logistic.h b/src/logistic.h new file mode 100644 index 0000000..a68ee09 --- /dev/null +++ b/src/logistic.h @@ -0,0 +1,70 @@ +#ifndef LOGISTIC_H_ /* Include guard */ +#define LOGISTIC_H_ + +/* Mixed interface */ +void logistic_mixed_pred(gsl_vector *beta // Vector of parameters length = 1 + Sum_k(C_k - 1) + Kc + ,gsl_matrix_int *X //Matrix Nobs x K + ,gsl_vector_int *nlev // Vector with number categories + ,gsl_matrix *Xc // continuous covariates Matrix Nobs x Kc + ,gsl_vector *yhat //Vector of prob. predicted by the logistic + ); + +int logistic_mixed_fit(gsl_vector *beta // Vector of parameters length = 1 + Sum_k(C_k - 1) + Kc + ,gsl_matrix_int *X //Matrix Nobs x K + ,gsl_vector_int *nlev // Vector with number categories + ,gsl_matrix *Xc // continuous covariates Matrix Nobs x Kc + ,gsl_vector *y //Vector of prob. to predict + ,double lambdaL1 // Regularization L1 0.0 if not used + ,double lambdaL2); // Regularization L2 0.0 if not used + +double fLogit_mixed(gsl_vector *beta + ,gsl_matrix_int *X + ,gsl_vector_int *nlev + ,gsl_matrix *Xc // continuous covariates Matrix Nobs x Kc + ,gsl_vector *y + ,double lambdaL1 + ,double lambdaL2); + + +/* Categorical only interface */ +void logistic_cat_pred(gsl_vector *beta // Vector of parameters length = 1 + Sum_k(C_k - 1) + Kc + ,gsl_matrix_int *X //Matrix Nobs x K + ,gsl_vector_int *nlev // Vector with number categories + ,gsl_vector *yhat //Vector of prob. predicted by the logistic + ); + +int logistic_cat_fit(gsl_vector *beta // Vector of parameters length = 1 + Sum_k(C_k - 1) + Kc + ,gsl_matrix_int *X //Matrix Nobs x K + ,gsl_vector_int *nlev // Vector with number categories + ,gsl_vector *y //Vector of prob. to predict + ,double lambdaL1 // Regularization L1 0.0 if not used + ,double lambdaL2); // Regularization L2 0.0 if not used + +double fLogit_cat(gsl_vector *beta + ,gsl_matrix_int *X + ,gsl_vector_int *nlev + ,gsl_vector *y + ,double lambdaL1 + ,double lambdaL2); + + +/* Continuous only interface */ +void logistic_cont_pred(gsl_vector *beta // Vector of parameters length = 1 + Sum_k(C_k - 1) + Kc + ,gsl_matrix *Xc // continuous covariates Matrix Nobs x Kc + ,gsl_vector *yhat //Vector of prob. predicted by the logistic + ); + +int logistic_cont_fit(gsl_vector *beta // Vector of parameters length = 1 + Sum_k(C_k - 1) + Kc + ,gsl_matrix *Xc // continuous covariates Matrix Nobs x Kc + ,gsl_vector *y //Vector of prob. to predict + ,double lambdaL1 // Regularization L1 0.0 if not used + ,double lambdaL2); // Regularization L2 0.0 if not used + +double fLogit_cont(gsl_vector *beta + ,gsl_matrix *Xc // continuous covariates Matrix Nobs x Kc + ,gsl_vector *y + ,double lambdaL1 + ,double lambdaL2); + + +#endif // LOGISTIC_H_ diff --git a/src/varcov.cpp b/src/varcov.cpp new file mode 100644 index 0000000..fdc6f10 --- /dev/null +++ b/src/varcov.cpp @@ -0,0 +1,482 @@ +/* + Genome-wide Efficient Mixed Model Association (GEMMA) + Copyright (C) 2011 Xiang Zhou + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see <http://www.gnu.org/licenses/>. +*/ + +#include <iostream> +#include <fstream> +#include <sstream> +#include <string> +#include <iomanip> +#include <bitset> +#include <vector> +#include <map> +#include <set> +#include <cstring> +#include <cmath> +#include <stdio.h> +#include <stdlib.h> + +#include "gsl/gsl_vector.h" +#include "gsl/gsl_matrix.h" +#include "gsl/gsl_linalg.h" +#include "gsl/gsl_blas.h" +#include "gsl/gsl_cdf.h" + +#include "lapack.h" +#include "gzstream.h" + +#ifdef FORCE_FLOAT +#include "param_float.h" +#include "varcov_float.h" +#include "io_float.h" +#include "mathfunc_float.h" +#else +#include "param.h" +#include "varcov.h" +#include "io.h" +#include "mathfunc.h" +#endif + + +using namespace std; + + + + +void VARCOV::CopyFromParam (PARAM &cPar) +{ + d_pace=cPar.d_pace; + + file_bfile=cPar.file_bfile; + file_geno=cPar.file_geno; + file_out=cPar.file_out; + path_out=cPar.path_out; + + time_opt=0.0; + + window_cm=cPar.window_cm; + window_bp=cPar.window_bp; + window_ns=cPar.window_ns; + + indicator_idv=cPar.indicator_idv; + indicator_snp=cPar.indicator_snp; + snpInfo=cPar.snpInfo; + + return; +} + + +void VARCOV::CopyToParam (PARAM &cPar) +{ + cPar.time_opt=time_opt; + return; +} + + + +//chr rs ps n_idv allele1 allele0 af var window_size r2 (r2_11,n_11,r2_12,n_12...r2_1m,n_1m) +void VARCOV::WriteCov (const int flag, const vector<SNPINFO> &snpInfo_sub, const vector<vector<double> > &Cov_mat) +{ + string file_cov; + file_cov=path_out+"/"+file_out; + file_cov+=".cor.txt"; + + ofstream outfile; + + if (flag==0) { + outfile.open (file_cov.c_str(), ofstream::out); + if (!outfile) {cout<<"error writing file: "<<file_cov<<endl; return;} + + outfile<<"chr"<<"\t"<<"rs"<<"\t"<<"ps"<<"\t"<<"n_mis" + <<"\t"<<"n_obs"<<"\t"<<"allele1"<<"\t"<<"allele0" + <<"\t"<<"af"<<"\t"<<"window_size" + <<"\t"<<"var"<<"\t"<<"cor"<<endl; + } else { + outfile.open (file_cov.c_str(), ofstream::app); + if (!outfile) {cout<<"error writing file: "<<file_cov<<endl; return;} + + for (size_t i=0; i<Cov_mat.size(); i++) { + outfile<<snpInfo_sub[i].chr<<"\t"<<snpInfo_sub[i].rs_number<<"\t"<<snpInfo_sub[i].base_position<<"\t"<<snpInfo_sub[i].n_miss<<"\t"<<snpInfo_sub[i].n_idv<<"\t"<<snpInfo_sub[i].a_minor<<"\t"<<snpInfo_sub[i].a_major<<"\t"<<fixed<<setprecision(3)<<snpInfo_sub[i].maf<<"\t"<<Cov_mat[i].size()-1<<"\t"; + outfile<<scientific<<setprecision(6)<<Cov_mat[i][0]<<"\t"; + + if (Cov_mat[i].size()==1) { + outfile<<"NA"; + } else { + for (size_t j=1; j<Cov_mat[i].size(); j++) { + if (j==(Cov_mat[i].size()-1)) { + outfile<<Cov_mat[i][j]; + } else { + outfile<<Cov_mat[i][j]<<","; + } + } + } + + outfile<<endl; + } + } + + outfile.close(); + outfile.clear(); + return; +} + + + + +// sort SNPs first based on chromosomes then based on bp (or cm, if bp is not available) +//if chr1=chr2 and bp1=bp2, then return with the same order +bool CompareSNPinfo (const SNPINFO &snpInfo1, const SNPINFO &snpInfo2) +{ + int c_chr=snpInfo1.chr.compare(snpInfo2.chr); + long int c_bp=snpInfo1.base_position-snpInfo2.base_position; + + if(c_chr<0) { + return true; + } else if (c_chr>0) { + return false; + } else { + if (c_bp<0) { + return true; + } else if (c_bp>0) { + return false; + } else { + return true; + } + } +} + + +// do not sort SNPs (because gzip files do not support random access) +// then calculate n_nb, the number of neighbours, for each snp +void VARCOV::CalcNB (vector<SNPINFO> &snpInfo_sort) +{ + // stable_sort(snpInfo_sort.begin(), snpInfo_sort.end(), CompareSNPinfo); + + size_t t2=0, n_nb=0; + for (size_t t=0; t<indicator_snp.size(); ++t) { + if (indicator_snp[t]==0) {continue;} + + if (snpInfo_sort[t].chr=="-9" || (snpInfo_sort[t].cM==-9 && window_cm!=0) || (snpInfo_sort[t].base_position==-9 && window_bp!=0) ) { + snpInfo_sort[t].n_nb=0; continue; + } + + if (t==indicator_snp.size()-1) {snpInfo_sort[t].n_nb=0; continue;} + + t2=t+1; n_nb=0; + + while (t2<indicator_snp.size() && snpInfo_sort[t2].chr==snpInfo_sort[t].chr && indicator_snp[t2]==0) {t2++;} + + while (t2<indicator_snp.size() && snpInfo_sort[t2].chr==snpInfo_sort[t].chr && (snpInfo_sort[t2].cM-snpInfo_sort[t].cM<window_cm || window_cm==0) && (snpInfo_sort[t2].base_position-snpInfo_sort[t].base_position<window_bp || window_bp==0) && (n_nb<window_ns|| window_ns==0) ) { + t2++; n_nb++; + while (t2<indicator_snp.size() && snpInfo_sort[t2].chr==snpInfo_sort[t].chr && indicator_snp[t2]==0) {t2++;} + } + + snpInfo_sort[t].n_nb=n_nb; + } + + return; +} + + + +//vector double is centered to have mean 0 +void Calc_Cor(vector<vector<double> > &X_mat, vector<double> &cov_vec) +{ + cov_vec.clear(); + + double v1, v2, r; + vector<double> x_vec=X_mat[0]; + + lapack_ddot(x_vec, x_vec, v1); + cov_vec.push_back(v1/(double)x_vec.size() ); + + for (size_t i=1; i<X_mat.size(); i++) { + lapack_ddot(X_mat[i], x_vec, r); + lapack_ddot(X_mat[i], X_mat[i], v2); + r/=sqrt(v1*v2); + + cov_vec.push_back(r); + } + + return; +} + +/* +//somehow this can produce nan for some snps; perhaps due to missing values? +//vector int is not centered, and have 0/1/2 values only +//missing value is -9; to calculate covariance, these values are replaced by +void Calc_Cor(const vector<vector<int> > &X_mat, vector<double> &cov_vec) +{ + cov_vec.clear(); + + int d1, d2, m1, m2, n1, n2, n12; + double m1d, m2d, m12d, v1d, v2d, v; + + vector<int> x_vec=X_mat[0]; + + //calculate m2 + m2=0; n2=0; + for (size_t j=0; j<x_vec.size(); j++) { + d2=x_vec[j]; + if (d2==-9) {continue;} + m2+=d2; n2++; + } + m2d=(double)m2/(double)n2; + + for (size_t i=0; i<X_mat.size(); i++) { + //calculate m1 + m1=0; n1=0; + for (size_t j=0; j<x_vec.size(); j++) { + d1=X_mat[i][j]; + if (d1==-9) {continue;} + m1+=d1; n1++; + } + m1d=(double)m1/(double)n1; + + //calculate v1 + m1=0; m2=0; m12d=0; n12=0; + for (size_t j=0; j<x_vec.size(); j++) { + d1=X_mat[i][j]; + if (d1==-9) { + n12++; + } else if (d1!=0) { + if (d1==1) {m12d+=1;} else if (d1==2) {m12d+=4;} else {m12d+=(double)(d1*d1);} + } + } + + v1d=((double)m12d+m1d*(double)m2+m2d*(double)m1+(double)n12*m1d*m2d)/(double)x_vec.size()-m1d*m2d; + + //calculate covariance + if (i!=0) { + m1=0; m2=0; m12d=0; n12=0; + for (size_t j=0; j<x_vec.size(); j++) { + d1=X_mat[i][j]; d2=x_vec[j]; + if (d1==-9 && d2==-9) { + n12++; + } else if (d1==-9) { + m2+=d2; + } else if (d2==-9) { + m1+=d1; + } else if (d1!=0 && d2!=0) { + if (d1==1) {m12d+=d2;} else if (d1==2) {m12d+=d2+d2;} else {m12d+=(double)(d1*d2);} + } + } + + v=((double)m12d+m1d*(double)m2+m2d*(double)m1+(double)n12*m1d*m2d)/(double)x_vec.size()-m1d*m2d; + v/=sqrt(v1d*v2d); + } else { + v2d=v1d; + v=v1d/(double)x_vec.size(); + } + + cov_vec.push_back(v); + } + return; +} +*/ + + +//read the genotype file again, calculate r2 between pairs of SNPs within a window, output the file every 10K SNPs +//the output results are sorted based on chr and bp +//output format similar to assoc.txt files (remember n_miss is replaced by n_idv) + +//r2 between the current SNP and every following SNPs within the window_size (which can vary if cM was used) +//read bimbam mean genotype file and calculate the covariance matrix for neighboring SNPs +//output values at 10000-SNP-interval +void VARCOV::AnalyzeBimbam () +{ + igzstream infile (file_geno.c_str(), igzstream::in); + if (!infile) {cout<<"error reading genotype file:"<<file_geno<<endl; return;} + + //calculate the number of right-hand-side neighbours for each SNP + vector<SNPINFO> snpInfo_sub; + CalcNB(snpInfo); + + size_t ni_test=0; + for (size_t i=0; i<indicator_idv.size(); i++) { + ni_test+=indicator_idv[i]; + } + + gsl_vector *geno=gsl_vector_alloc (ni_test); + double geno_mean; + + vector<double> x_vec, cov_vec; + vector<vector<double> > X_mat, Cov_mat; + + for (size_t i=0; i<ni_test; i++) { + x_vec.push_back(0); + } + + WriteCov (0, snpInfo_sub, Cov_mat); + + size_t t2=0, inc; + int n_nb=0; + + for (size_t t=0; t<indicator_snp.size(); ++t) { + if (t%d_pace==0 || t==(indicator_snp.size()-1)) {ProgressBar ("Reading SNPs ", t, indicator_snp.size()-1);} + if (indicator_snp[t]==0) {continue;} + // if (t>=2) {break;} + + if (X_mat.size()==0) { + n_nb=snpInfo[t].n_nb+1; + } else { + n_nb=snpInfo[t].n_nb-n_nb+1; + } + + for (int i=0; i<n_nb; i++) { + if (X_mat.size()==0) {t2=t;} + + //read a line of the snp is filtered out + inc=0; + while (t2<indicator_snp.size() && indicator_snp[t2]==0) { + t2++; inc++; + } + + Bimbam_ReadOneSNP (inc, indicator_idv, infile, geno, geno_mean); + gsl_vector_add_constant (geno, -1.0*geno_mean); + + for (size_t j=0; j<geno->size; j++) { + x_vec[j]=gsl_vector_get(geno, j); + } + X_mat.push_back(x_vec); + + t2++; + } + + n_nb=snpInfo[t].n_nb; + + Calc_Cor(X_mat, cov_vec); + Cov_mat.push_back(cov_vec); + snpInfo_sub.push_back(snpInfo[t]); + + X_mat.erase(X_mat.begin()); + + //write out var/cov values + if (Cov_mat.size()==10000) { + WriteCov (1, snpInfo_sub, Cov_mat); + Cov_mat.clear(); + snpInfo_sub.clear(); + } + } + + if (Cov_mat.size()!=0) { + WriteCov (1, snpInfo_sub, Cov_mat); + Cov_mat.clear(); + snpInfo_sub.clear(); + } + + gsl_vector_free(geno); + + infile.close(); + infile.clear(); + + return; +} + + + + +void VARCOV::AnalyzePlink () +{ + string file_bed=file_bfile+".bed"; + ifstream infile (file_bed.c_str(), ios::binary); + if (!infile) {cout<<"error reading bed file:"<<file_bed<<endl; return;} + + //calculate the number of right-hand-side neighbours for each SNP + vector<SNPINFO> snpInfo_sub; + CalcNB(snpInfo); + + size_t ni_test=0; + for (size_t i=0; i<indicator_idv.size(); i++) { + ni_test+=indicator_idv[i]; + } + + gsl_vector *geno=gsl_vector_alloc (ni_test); + double geno_mean; + + vector<double> x_vec, cov_vec; + vector<vector<double> > X_mat, Cov_mat; + + for (size_t i=0; i<ni_test; i++) { + x_vec.push_back(0); + } + + WriteCov (0, snpInfo_sub, Cov_mat); + + size_t t2=0, inc; + int n_nb=0; + + for (size_t t=0; t<indicator_snp.size(); ++t) { + if (t%d_pace==0 || t==(indicator_snp.size()-1)) {ProgressBar ("Reading SNPs ", t, indicator_snp.size()-1);} + if (indicator_snp[t]==0) {continue;} + // if (t>=2) {break;} + + if (X_mat.size()==0) { + n_nb=snpInfo[t].n_nb+1; + } else { + n_nb=snpInfo[t].n_nb-n_nb+1; + } + + for (int i=0; i<n_nb; i++) { + if (X_mat.size()==0) {t2=t;} + + //read a line of the snp is filtered out + inc=0; + while (t2<indicator_snp.size() && indicator_snp[t2]==0) { + t2++; inc++; + } + + Plink_ReadOneSNP (t2, indicator_idv, infile, geno, geno_mean); + gsl_vector_add_constant (geno, -1.0*geno_mean); + + for (size_t j=0; j<geno->size; j++) { + x_vec[j]=gsl_vector_get(geno, j); + } + X_mat.push_back(x_vec); + + t2++; + } + + n_nb=snpInfo[t].n_nb; + + Calc_Cor(X_mat, cov_vec); + Cov_mat.push_back(cov_vec); + snpInfo_sub.push_back(snpInfo[t]); + + X_mat.erase(X_mat.begin()); + + //write out var/cov values + if (Cov_mat.size()==10000) { + WriteCov (1, snpInfo_sub, Cov_mat); + Cov_mat.clear(); + snpInfo_sub.clear(); + } + } + + if (Cov_mat.size()!=0) { + WriteCov (1, snpInfo_sub, Cov_mat); + Cov_mat.clear(); + snpInfo_sub.clear(); + } + + gsl_vector_free(geno); + + infile.close(); + infile.clear(); + + return; +} diff --git a/src/varcov.h b/src/varcov.h new file mode 100644 index 0000000..b380b8c --- /dev/null +++ b/src/varcov.h @@ -0,0 +1,72 @@ +/* + Genome-wide Efficient Mixed Model Association (GEMMA) + Copyright (C) 2011 Xiang Zhou + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see <http://www.gnu.org/licenses/>. +*/ + +#ifndef __VARCOV_H__ +#define __VARCOV_H__ + +#include "gsl/gsl_vector.h" +#include "gsl/gsl_matrix.h" + + +#ifdef FORCE_FLOAT +#include "param_float.h" +#include "io_float.h" +#else +#include "param.h" +#include "io.h" +#endif + +using namespace std; + + + + +class VARCOV { + +public: + // IO related parameters + string file_out; + string path_out; + string file_geno; + string file_bfile; + int d_pace; + + vector<int> indicator_idv; + vector<int> indicator_snp; + + vector<SNPINFO> snpInfo; + + double time_opt; + + // Class specific parameters + double window_cm; + size_t window_bp; + size_t window_ns; + + // Main functions + void CopyFromParam (PARAM &cPar); + void CopyToParam (PARAM &cPar); + void CalcNB (vector<SNPINFO> &snpInfo_sort); + void WriteCov (const int flag, const vector<SNPINFO> &snpInfo_sub, const vector<vector<double> > &Cov_mat); + void AnalyzeBimbam (); + void AnalyzePlink (); +}; + +#endif + + |