From bfe06c298c1b9e8be12df91e4f1c9a8e5c459562 Mon Sep 17 00:00:00 2001 From: Peter Carbonetto Date: Fri, 2 Jun 2017 00:07:12 -0500 Subject: Removed FORCE_FLOAT from lmm.h,cpp. --- src/lapack.cpp | 32 ++--- src/lm.cpp | 416 +++++++++++++++++++++++++++++---------------------------- src/lm.h | 54 ++++---- 3 files changed, 257 insertions(+), 245 deletions(-) (limited to 'src') diff --git a/src/lapack.cpp b/src/lapack.cpp index 2bbdf62..01d2039 100644 --- a/src/lapack.cpp +++ b/src/lapack.cpp @@ -62,14 +62,14 @@ void lapack_float_cholesky_decomp (gsl_matrix_float *A) { char UPLO='L'; if (N!=(int)A->size2) { - cout << "Matrix needs to be symmetric and same dimension in" << + cout << "Matrix needs to be symmetric and same dimension in " << "lapack_cholesky_decomp." << endl; return; } spotrf_(&UPLO, &N, A->data, &LDA, &INFO); if (INFO!=0) { - cout << "Cholesky decomposition unsuccessful in" << + cout << "Cholesky decomposition unsuccessful in " << "lapack_cholesky_decomp." << endl; return; } @@ -83,14 +83,14 @@ void lapack_cholesky_decomp (gsl_matrix *A) { char UPLO='L'; if (N!=(int)A->size2) { - cout << "Matrix needs to be symmetric and same dimension in" << + cout << "Matrix needs to be symmetric and same dimension in " << "lapack_cholesky_decomp." << endl; return; } dpotrf_(&UPLO, &N, A->data, &LDA, &INFO); if (INFO!=0) { - cout << "Cholesky decomposition unsuccessful in" << + cout << "Cholesky decomposition unsuccessful in " << "lapack_cholesky_decomp."<size2 || N!=LDB) { - cout << "Matrix needs to be symmetric and same dimension in" << + cout << "Matrix needs to be symmetric and same dimension in " <size2 || N!=LDB) { - cout << "Matrix needs to be symmetric and same dimension in" << + cout << "Matrix needs to be symmetric and same dimension in " << "lapack_cholesky_solve." << endl; return; } @@ -236,7 +236,7 @@ void lapack_float_eigen_symmv (gsl_matrix_float *A, gsl_vector_float *eval, char JOBZ='V', UPLO='L'; if (N!=(int)A->size2 || N!=(int)eval->size) { - cout << "Matrix needs to be symmetric and same" << + cout << "Matrix needs to be symmetric and same " << "dimension in lapack_eigen_symmv."<data, &LDA, eval->data, WORK, &LWORK, &INFO); if (INFO!=0) { - cout << "Eigen decomposition unsuccessful in" << + cout << "Eigen decomposition unsuccessful in " << "lapack_eigen_symmv."<size2 || N!=(int)eval->size) { - cout << "Matrix needs to be symmetric and same" << + cout << "Matrix needs to be symmetric and same " << "dimension in lapack_float_eigen_symmv." << endl; return; } @@ -282,7 +282,7 @@ void lapack_float_eigen_symmv (gsl_matrix_float *A, gsl_vector_float *eval, evec->data, &LDZ, ISUPPZ, WORK_temp, &LWORK, IWORK_temp, &LIWORK, &INFO); if (INFO!=0) { - cout << "Work space estimate unsuccessful in" << + cout << "Work space estimate unsuccessful in " << "lapack_float_eigen_symmv." << endl; return; } @@ -295,7 +295,7 @@ void lapack_float_eigen_symmv (gsl_matrix_float *A, gsl_vector_float *eval, &VU, &IL, &IU, &ABSTOL, &M, eval->data, evec->data, &LDZ, ISUPPZ, WORK, &LWORK, IWORK, &LIWORK, &INFO); if (INFO!=0) { - cout << "Eigen decomposition unsuccessful in" << + cout << "Eigen decomposition unsuccessful in " << "lapack_float_eigen_symmv." << endl; return; } @@ -321,7 +321,7 @@ void lapack_eigen_symmv (gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec, char JOBZ='V', UPLO='L'; if (N!=(int)A->size2 || N!=(int)eval->size) { - cout << "Matrix needs to be symmetric and same" << + cout << "Matrix needs to be symmetric and same " << "dimension in lapack_eigen_symmv." << endl; return; } @@ -331,7 +331,7 @@ void lapack_eigen_symmv (gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec, dsyev_(&JOBZ, &UPLO, &N, A->data, &LDA, eval->data, WORK, &LWORK, &INFO); if (INFO!=0) { - cout<<"Eigen decomposition unsuccessful in" << + cout<<"Eigen decomposition unsuccessful in " << "lapack_eigen_symmv." << endl; return; } @@ -352,7 +352,7 @@ void lapack_eigen_symmv (gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec, int IL=0, IU=0, M; if (N!=(int)A->size2 || N!=(int)eval->size) { - cout << "Matrix needs to be symmetric and same" << + cout << "Matrix needs to be symmetric and same " << "dimension in lapack_eigen_symmv." << endl; return; } @@ -367,7 +367,7 @@ void lapack_eigen_symmv (gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec, &LDZ, ISUPPZ, WORK_temp, &LWORK, IWORK_temp, &LIWORK, &INFO); if (INFO!=0) { - cout << "Work space estimate unsuccessful in" << + cout << "Work space estimate unsuccessful in " << "lapack_eigen_symmv." << endl; return; } @@ -380,7 +380,7 @@ void lapack_eigen_symmv (gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec, &IL, &IU, &ABSTOL, &M, eval->data, evec->data, &LDZ, ISUPPZ, WORK, &LWORK, IWORK, &LIWORK, &INFO); if (INFO!=0) { - cout << "Eigen decomposition unsuccessful in" << + cout << "Eigen decomposition unsuccessful in " << "lapack_eigen_symmv." << endl; return; } diff --git a/src/lm.cpp b/src/lm.cpp index f8cb974..d3ad7f3 100644 --- a/src/lm.cpp +++ b/src/lm.cpp @@ -1,6 +1,6 @@ /* Genome-wide Efficient Mixed Model Association (GEMMA) - Copyright (C) 2011 Xiang Zhou + Copyright (C) 2011-2017 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 @@ -14,9 +14,7 @@ You should have received a copy of the GNU General Public License along with this program. If not, see . - */ - - +*/ #include #include @@ -27,6 +25,7 @@ #include #include #include +#include #include #include @@ -35,7 +34,6 @@ #include "gsl/gsl_linalg.h" #include "gsl/gsl_blas.h" - #include "gsl/gsl_cdf.h" #include "gsl/gsl_roots.h" #include "gsl/gsl_min.h" @@ -44,22 +42,11 @@ #include "eigenlib.h" #include "gzstream.h" #include "lapack.h" - -#ifdef FORCE_FLOAT -#include "lm_float.h" -#else #include "lm.h" -#endif - using namespace std; - - - - -void LM::CopyFromParam (PARAM &cPar) -{ +void LM::CopyFromParam (PARAM &cPar) { a_mode=cPar.a_mode; d_pace=cPar.d_pace; @@ -89,26 +76,22 @@ void LM::CopyFromParam (PARAM &cPar) return; } - -void LM::CopyToParam (PARAM &cPar) -{ +void LM::CopyToParam (PARAM &cPar) { cPar.time_opt=time_opt; - cPar.ng_test=ng_test; - return; } - - -void LM::WriteFiles () -{ +void LM::WriteFiles () { string file_str; file_str=path_out+"/"+file_out; file_str+=".assoc.txt"; ofstream outfile (file_str.c_str(), ofstream::out); - if (!outfile) {cout<<"error writing file: "<::size_type t=0; tsize; double d; @@ -200,9 +205,8 @@ void CalcvPv(const gsl_matrix *WtWi, const gsl_vector *Wty, const gsl_vector *Wt return; } - -void CalcvPv(const gsl_matrix *WtWi, const gsl_vector *Wty, const gsl_vector *y, double &yPwy) -{ +void CalcvPv(const gsl_matrix *WtWi, const gsl_vector *Wty, + const gsl_vector *y, double &yPwy) { size_t c_size=Wty->size; double d; @@ -219,11 +223,11 @@ void CalcvPv(const gsl_matrix *WtWi, const gsl_vector *Wty, const gsl_vector *y, return; } - - -//calculate p values and beta/se in a linear model -void LmCalcP (const size_t test_mode, const double yPwy, const double xPwy, const double xPwx, const double df, const size_t n_size, double &beta, double &se, double &p_wald, double &p_lrt, double &p_score) -{ +// Calculate p-values and beta/se in a linear model. +void LmCalcP (const size_t test_mode, const double yPwy, + const double xPwy, const double xPwx, const double df, + const size_t n_size, double &beta, double &se, + double &p_wald, double &p_lrt, double &p_score) { double yPxy=yPwy-xPwy*xPwy/xPwx; double se_wald, se_score; @@ -240,13 +244,12 @@ void LmCalcP (const size_t test_mode, const double yPwy, const double xPwy, cons return; } - - - -void LM::AnalyzeGene (const gsl_matrix *W, const gsl_vector *x) -{ +void LM::AnalyzeGene (const gsl_matrix *W, const gsl_vector *x) { ifstream infile (file_gene.c_str(), ifstream::in); - if (!infile) {cout<<"error reading gene expression file:"<size1-(double)W->size2-1.0; @@ -278,12 +281,14 @@ void LM::AnalyzeGene (const gsl_matrix *W, const gsl_vector *x) gsl_blas_dgemv (CblasTrans, 1.0, W, x, 0.0, Wtx); CalcvPv(WtWi, Wtx, x, xPwx); - //header + // Header. getline(infile, line); for (size_t t=0; tsize1, beta, se, p_wald, p_lrt, p_score); + LmCalcP (a_mode-50, yPwy, xPwy, xPwx, df, W->size1, + beta, se, p_wald, p_lrt, p_score); time_opt+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0); - //store summary data + // Store summary data. SUMSTAT SNPs={beta, se, 0.0, 0.0, p_wald, p_lrt, p_score}; sumStat.push_back(SNPs); } @@ -327,17 +333,14 @@ void LM::AnalyzeGene (const gsl_matrix *W, const gsl_vector *x) return; } - - - // WJA added -#include -void LM::Analyzebgen (const gsl_matrix *W, const gsl_vector *y) -{ +void LM::Analyzebgen (const gsl_matrix *W, const gsl_vector *y) { string file_bgen=file_oxford+".bgen"; ifstream infile (file_bgen.c_str(), ios::binary); - if (!infile) {cout<<"error reading bgen file:"<size1-(double)W->size2-1.0; @@ -369,7 +372,7 @@ void LM::Analyzebgen (const gsl_matrix *W, const gsl_vector *y) gsl_blas_dgemv (CblasTrans, 1.0, W, y, 0.0, Wty); CalcvPv(WtWi, Wty, y, yPwy); - // read in header + // Read in header. uint32_t bgen_snp_block_offset; uint32_t bgen_header_length; uint32_t bgen_nsamples; @@ -387,11 +390,11 @@ void LM::Analyzebgen (const gsl_matrix *W, const gsl_vector *y) infile.read(reinterpret_cast(&bgen_flags),4); bgen_snp_block_offset-=4; bool CompressedSNPBlocks=bgen_flags&0x1; -// bool LongIds=bgen_flags&0x4; infile.ignore(bgen_snp_block_offset); - double bgen_geno_prob_AA, bgen_geno_prob_AB, bgen_geno_prob_BB, bgen_geno_prob_non_miss; + double bgen_geno_prob_AA, bgen_geno_prob_AB; + double bgen_geno_prob_BB, bgen_geno_prob_non_miss; uint32_t bgen_N; uint16_t bgen_LS; @@ -407,17 +410,16 @@ void LM::Analyzebgen (const gsl_matrix *W, const gsl_vector *y) string id; string rs; string chr; - std::cout<<"Warning: WJA hard coded SNP missingness threshold of 10%"<1) {break;} - if (t%d_pace==0 || t==(ns_total-1)) {ProgressBar ("Reading SNPs ", t, ns_total-1);} - // read SNP header + std::cout << "Warning: WJA hard coded SNP missingness " << + "threshold of 10%" << std::endl; + + // Start reading genotypes and analyze. + for (size_t t=0; t(&bgen_LB),4); bgen_B_allele.resize(bgen_LB); infile.read(&bgen_B_allele[0], bgen_LB); - - - uint16_t unzipped_data[3*bgen_N]; if (indicator_snp[t]==0) { if(CompressedSNPBlocks) - infile.read(reinterpret_cast(&bgen_P),4); + infile.read(reinterpret_cast(&bgen_P),4); else - bgen_P=6*bgen_N; + bgen_P=6*bgen_N; infile.ignore(static_cast(bgen_P)); continue; } - - if(CompressedSNPBlocks) - { - - + if(CompressedSNPBlocks) { infile.read(reinterpret_cast(&bgen_P),4); uint8_t zipped_data[bgen_P]; unzipped_data_size=6*bgen_N; - infile.read(reinterpret_cast(zipped_data),bgen_P); + infile.read(reinterpret_cast(zipped_data), + bgen_P); - int result=uncompress(reinterpret_cast(unzipped_data), reinterpret_cast(&unzipped_data_size), reinterpret_cast(zipped_data), static_cast (bgen_P)); + int result= + uncompress(reinterpret_cast(unzipped_data), + reinterpret_cast(&unzipped_data_size), + reinterpret_cast(zipped_data), + static_cast (bgen_P)); assert(result == Z_OK); } @@ -485,7 +484,8 @@ void LM::Analyzebgen (const gsl_matrix *W, const gsl_vector *y) { bgen_P=6*bgen_N; - infile.read(reinterpret_cast(unzipped_data),bgen_P); + infile.read(reinterpret_cast(unzipped_data), + bgen_P); } x_mean=0.0; c_phen=0; n_miss=0; @@ -494,23 +494,32 @@ void LM::Analyzebgen (const gsl_matrix *W, const gsl_vector *y) if (indicator_idv[i]==0) {continue;} - bgen_geno_prob_AA=static_cast(unzipped_data[i*3])/32768.0; - bgen_geno_prob_AB=static_cast(unzipped_data[i*3+1])/32768.0; - bgen_geno_prob_BB=static_cast(unzipped_data[i*3+2])/32768.0; + bgen_geno_prob_AA= + static_cast(unzipped_data[i*3])/32768.0; + bgen_geno_prob_AB= + static_cast(unzipped_data[i*3+1])/32768.0; + bgen_geno_prob_BB= + static_cast(unzipped_data[i*3+2])/32768.0; + // WJA - bgen_geno_prob_non_miss=bgen_geno_prob_AA+bgen_geno_prob_AB+bgen_geno_prob_BB; - if (bgen_geno_prob_non_miss<0.9) {gsl_vector_set(x_miss, c_phen, 0.0); n_miss++;} - else { - - bgen_geno_prob_AA/=bgen_geno_prob_non_miss; - bgen_geno_prob_AB/=bgen_geno_prob_non_miss; - bgen_geno_prob_BB/=bgen_geno_prob_non_miss; - - geno=2.0*bgen_geno_prob_BB+bgen_geno_prob_AB; + bgen_geno_prob_non_miss= + bgen_geno_prob_AA + + bgen_geno_prob_AB + + bgen_geno_prob_BB; + if (bgen_geno_prob_non_miss<0.9) { + gsl_vector_set(x_miss, c_phen, 0.0); + n_miss++; + } + else { + bgen_geno_prob_AA/=bgen_geno_prob_non_miss; + bgen_geno_prob_AB/=bgen_geno_prob_non_miss; + bgen_geno_prob_BB/=bgen_geno_prob_non_miss; + + geno=2.0*bgen_geno_prob_BB+bgen_geno_prob_AB; - gsl_vector_set(x, c_phen, geno); - gsl_vector_set(x_miss, c_phen, 1.0); - x_mean+=geno; + gsl_vector_set(x, c_phen, geno); + gsl_vector_set(x_miss, c_phen, 1.0); + x_mean+=geno; } c_phen++; } @@ -518,24 +527,23 @@ void LM::Analyzebgen (const gsl_matrix *W, const gsl_vector *y) x_mean/=static_cast(ni_test-n_miss); for (size_t i=0; i1) { - //gsl_vector_set(x, i, 2-geno); - //} } - - //calculate statistics + // Calculate statistics. time_start=clock(); gsl_blas_dgemv(CblasTrans, 1.0, W, x, 0.0, Wtx); CalcvPv(WtWi, Wty, Wtx, y, x, xPwy, xPwx); - LmCalcP (a_mode-50, yPwy, xPwy, xPwx, df, W->size1, beta, se, p_wald, p_lrt, p_score); + LmCalcP (a_mode-50, yPwy, xPwy, xPwx, df, W->size1, + beta, se, p_wald, p_lrt, p_score); time_opt+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0); - //store summary data + // Store summary data. SUMSTAT SNPs={beta, se, 0.0, 0.0, p_wald, p_lrt, p_score}; sumStat.push_back(SNPs); } @@ -556,13 +564,12 @@ void LM::Analyzebgen (const gsl_matrix *W, const gsl_vector *y) return; } - - -void LM::AnalyzeBimbam (const gsl_matrix *W, const gsl_vector *y) -{ +void LM::AnalyzeBimbam (const gsl_matrix *W, const gsl_vector *y) { igzstream infile (file_geno.c_str(), igzstream::in); - // ifstream infile (file_geno.c_str(), ifstream::in); - if (!infile) {cout<<"error reading genotype file:"<size1-(double)W->size2-1.0; @@ -594,11 +601,12 @@ void LM::AnalyzeBimbam (const gsl_matrix *W, const gsl_vector *y) gsl_blas_dgemv (CblasTrans, 1.0, W, y, 0.0, Wty); CalcvPv(WtWi, Wty, y, yPwy); - //start reading genotypes and analyze + // Start reading genotypes and analyze. for (size_t t=0; t1) {break;} getline(infile, line); - if (t%d_pace==0 || t==(ns_total-1)) {ProgressBar ("Reading SNPs ", t, ns_total-1);} + if (t%d_pace==0 || t==(ns_total-1)) { + ProgressBar ("Reading SNPs ", t, ns_total-1); + } if (indicator_snp[t]==0) {continue;} ch_ptr=strtok ((char *)line.c_str(), " , \t"); @@ -611,7 +619,10 @@ void LM::AnalyzeBimbam (const gsl_matrix *W, const gsl_vector *y) ch_ptr=strtok (NULL, " , \t"); if (indicator_idv[i]==0) {continue;} - if (strcmp(ch_ptr, "NA")==0) {gsl_vector_set(x_miss, c_phen, 0.0); n_miss++;} + if (strcmp(ch_ptr, "NA")==0) { + gsl_vector_set(x_miss, c_phen, 0.0); + n_miss++; + } else { geno=atof(ch_ptr); @@ -625,23 +636,23 @@ void LM::AnalyzeBimbam (const gsl_matrix *W, const gsl_vector *y) x_mean/=(double)(ni_test-n_miss); for (size_t i=0; i1) { - //gsl_vector_set(x, i, 2-geno); - //} } - //calculate statistics + // Calculate statistics. time_start=clock(); gsl_blas_dgemv(CblasTrans, 1.0, W, x, 0.0, Wtx); CalcvPv(WtWi, Wty, Wtx, y, x, xPwy, xPwx); - LmCalcP (a_mode-50, yPwy, xPwy, xPwx, df, W->size1, beta, se, p_wald, p_lrt, p_score); + LmCalcP (a_mode-50, yPwy, xPwy, xPwx, df, W->size1, + beta, se, p_wald, p_lrt, p_score); time_opt+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0); - //store summary data + // Store summary data. SUMSTAT SNPs={beta, se, 0.0, 0.0, p_wald, p_lrt, p_score}; sumStat.push_back(SNPs); } @@ -662,17 +673,13 @@ void LM::AnalyzeBimbam (const gsl_matrix *W, const gsl_vector *y) return; } - - - - - - -void LM::AnalyzePlink (const gsl_matrix *W, const gsl_vector *y) -{ +void LM::AnalyzePlink (const gsl_matrix *W, const gsl_vector *y) { string file_bed=file_bfile+".bed"; ifstream infile (file_bed.c_str(), ios::binary); - if (!infile) {cout<<"error reading bed file:"<size1-(double)W->size2-1.0; @@ -703,42 +710,62 @@ void LM::AnalyzePlink (const gsl_matrix *W, const gsl_vector *y) gsl_blas_dgemv (CblasTrans, 1.0, W, y, 0.0, Wty); CalcvPv(WtWi, Wty, y, yPwy); - //calculate n_bit and c, the number of bit for each snp + // Calculate n_bit and c, the number of bit for each SNP. if (ni_total%4==0) {n_bit=ni_total/4;} - else {n_bit=ni_total/4+1; } + else {n_bit=ni_total/4+1;} - //print the first three majic numbers + // Print the first three magic numbers. for (int i=0; i<3; ++i) { infile.read(ch,1); b=ch[0]; } for (vector::size_type t=0; t1) { - //gsl_vector_set(x, i, 2-geno); - //} + if (geno==-9) { + gsl_vector_set(x, i, x_mean); + geno=x_mean; + } } - //calculate statistics + // Calculate statistics. time_start=clock(); gsl_blas_dgemv (CblasTrans, 1.0, W, x, 0.0, Wtx); CalcvPv(WtWi, Wty, Wtx, y, x, xPwy, xPwx); - LmCalcP (a_mode-50, yPwy, xPwy, xPwx, df, W->size1, beta, se, p_wald, p_lrt, p_score); + LmCalcP (a_mode-50, yPwy, xPwy, xPwx, df, W->size1, + beta, se, p_wald, p_lrt, p_score); //store summary data SUMSTAT SNPs={beta, se, 0.0, 0.0, p_wald, p_lrt, p_score}; @@ -781,25 +809,9 @@ void LM::AnalyzePlink (const gsl_matrix *W, const gsl_vector *y) return; } - - - - - - - - - - - - - - - - -//make sure that both y and X are centered already -void MatrixCalcLmLR (const gsl_matrix *X, const gsl_vector *y, vector > &pos_loglr) -{ +// Make sure that both y and X are centered already. +void MatrixCalcLmLR (const gsl_matrix *X, const gsl_vector *y, + vector > &pos_loglr) { double yty, xty, xtx, log_lr; gsl_blas_ddot(y, y, &yty); diff --git a/src/lm.h b/src/lm.h index 656dd52..fac84e1 100644 --- a/src/lm.h +++ b/src/lm.h @@ -1,6 +1,6 @@ /* Genome-wide Efficient Mixed Model Association (GEMMA) - Copyright (C) 2011 Xiang Zhou + Copyright (C) 2011-2017 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 @@ -13,33 +13,25 @@ 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 . - */ + along with this program. If not, see . +*/ #ifndef __LM_H__ #define __LM_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 LM { public: - // IO related parameters - int a_mode; //analysis mode, 50+1/2/3/4 for Frequentist tests - size_t d_pace; //display pace + // IO-related parameters. + int a_mode; // Analysis mode: 50+1/2/3/4 for Frequentist tests. + size_t d_pace; // Display pace. string file_bfile; string file_geno; @@ -49,31 +41,39 @@ public: string file_gene; - // Summary statistics - size_t ni_total, ni_test; //number of individuals - size_t ns_total, ns_test; //number of snps - size_t ng_total, ng_test; //number of genes + // Summary statistics. + size_t ni_total, ni_test; // Number of individuals. + size_t ns_total, ns_test; // Number of SNPs. + size_t ng_total, ng_test; // Number of genes. size_t n_cvt; - double time_opt; //time spent + double time_opt; // Time spent. - vector indicator_idv; //indicator for individuals (phenotypes), 0 missing, 1 available for analysis - vector indicator_snp; //sequence indicator for SNPs: 0 ignored because of (a) maf, (b) miss, (c) non-poly; 1 available for analysis + // Indicator for individuals (phenotypes): 0 missing, 1 + // available for analysis. + vector indicator_idv; - vector snpInfo; //record SNP information + // Sequence indicator for SNPs: 0 ignored because of (a) maf, + // (b) miss, (c) non-poly; 1 available for analysis. + vector indicator_snp; - // Not included in PARAM - vector sumStat; //Output SNPSummary Data + vector snpInfo; // Record SNP information. - // Main functions + // Not included in PARAM. + vector sumStat; // Output SNPSummary Data. + + // Main functions. void CopyFromParam (PARAM &cPar); void CopyToParam (PARAM &cPar); void AnalyzeGene (const gsl_matrix *W, const gsl_vector *x); void AnalyzePlink (const gsl_matrix *W, const gsl_vector *y); void AnalyzeBimbam (const gsl_matrix *W, const gsl_vector *y); - // WJA added + // WJA added. void Analyzebgen (const gsl_matrix *W, const gsl_vector *y); void WriteFiles (); }; -void MatrixCalcLmLR (const gsl_matrix *X, const gsl_vector *y, vector > &pos_loglr); + +void MatrixCalcLmLR (const gsl_matrix *X, const gsl_vector *y, + vector > &pos_loglr); + #endif -- cgit v1.2.3