1 /* 2 Genome-wide Efficient Mixed Model Association (GEMMA) 3 Copyright (C) 2011-2017, Xiang Zhou 4 5 This program is free software: you can redistribute it and/or modify 6 it under the terms of the GNU General Public License as published by 7 the Free Software Foundation, either version 3 of the License, or 8 (at your option) any later version. 9 10 This program is distributed in the hope that it will be useful, 11 but WITHOUT ANY WARRANTY; without even the implied warranty of 12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 GNU General Public License for more details. 14 15 You should have received a copy of the GNU General Public License 16 along with this program. If not, see <http://www.gnu.org/licenses/>. 17 */ 18 19 #ifndef __MVLMM_H__ 20 #define __MVLMM_H__ 21 22 #include "gsl/gsl_matrix.h" 23 #include "gsl/gsl_vector.h" 24 #include "gemma_io.h" 25 #include "param.h" 26 27 using namespace std; 28 29 class MVLMM { 30 31 public: 32 // IO-related parameters. 33 int a_mode; // Analysis mode: 1/2/3/4 for Frequentist tests. 34 size_t d_pace; // Display pace. 35 36 string file_bfile; 37 string file_geno; 38 string file_oxford; 39 string file_out; 40 string path_out; 41 42 // MVLMM-related parameters. 43 double l_min; 44 double l_max; 45 size_t n_region; 46 double logl_remle_H0, logl_mle_H0; 47 vector<double> Vg_remle_null, Ve_remle_null, Vg_mle_null, Ve_mle_null; 48 vector<double> VVg_remle_null, VVe_remle_null, VVg_mle_null; 49 vector<double> VVe_mle_null; 50 vector<double> beta_remle_null, se_beta_remle_null, beta_mle_null; 51 vector<double> se_beta_mle_null; 52 double p_nr; 53 size_t em_iter, nr_iter; 54 double em_prec, nr_prec; 55 size_t crt; 56 57 // Summary statistics. 58 size_t ni_total, ni_test; // Number of individuals. 59 size_t ns_total, ns_test; // Number of SNPs. 60 size_t n_cvt; 61 size_t n_ph; 62 double time_UtX; // Time spent on optimization iterations. 63 double time_opt; // Time spent on optimization iterations. 64 65 // Indicator for individuals (phenotypes): 0 missing, 1 66 // available for analysis. 67 vector<int> indicator_idv; 68 69 // Sequence indicator for SNPs: 0 ignored because of (a) maf, 70 // (b) miss, (c) non-poly; 1 available for analysis. 71 vector<int> indicator_snp; 72 73 vector<SNPINFO> snpInfo; // Record SNP information. 74 75 // Not included in PARAM. 76 vector<MPHSUMSTAT> sumStat; // Output SNPSummary Data. 77 78 // Main functions 79 void CopyFromParam(PARAM &cPar); 80 void CopyToParam(PARAM &cPar); 81 void AnalyzeBimbam(const gsl_matrix *U, const gsl_vector *eval, 82 const gsl_matrix *UtW, const gsl_matrix *UtY); 83 void AnalyzePlink(const gsl_matrix *U, const gsl_vector *eval, 84 const gsl_matrix *UtW, const gsl_matrix *UtY); 85 void Analyzebgen(const gsl_matrix *U, const gsl_vector *eval, 86 const gsl_matrix *UtW, const gsl_matrix *UtY); 87 void AnalyzeBimbamGXE(const gsl_matrix *U, const gsl_vector *eval, 88 const gsl_matrix *UtW, const gsl_matrix *UtY, 89 const gsl_vector *env); 90 void AnalyzePlinkGXE(const gsl_matrix *U, const gsl_vector *eval, 91 const gsl_matrix *UtW, const gsl_matrix *UtY, 92 const gsl_vector *env); 93 void WriteFiles(); 94 }; 95 96 void CalcMvLmmVgVeBeta(const gsl_vector *eval, const gsl_matrix *UtW, 97 const gsl_matrix *UtY, const size_t em_iter, 98 const size_t nr_iter, const double em_prec, 99 const double nr_prec, const double l_min, 100 const double l_max, const size_t n_region, 101 gsl_matrix *V_g, gsl_matrix *V_e, gsl_matrix *B, 102 gsl_matrix *se_B); 103 104 #endif 105