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 __BSLMMDAP_H__ 20 #define __BSLMMDAP_H__ 21 22 #include "param.h" 23 #include <gsl/gsl_randist.h> 24 // #include <gsl/gsl_rng.h> 25 #include <map> 26 #include <vector> 27 28 using namespace std; 29 30 class BSLMMDAP { 31 32 public: 33 // IO-related parameters. 34 int a_mode; 35 size_t d_pace; 36 37 string file_bfile; 38 string file_geno; 39 string file_out; 40 string path_out; 41 42 // LMM related parameters 43 double pve_null; 44 double pheno_mean; 45 46 // BSLMM MCMC related parameters 47 // long int randseed; 48 double trace_G; 49 50 HYPBSLMM cHyp_initial; 51 52 // Summary statistics 53 size_t ni_total, ns_total; // Number of total individuals and SNPs. 54 size_t ni_test, ns_test; // Number of individuals and SNPs 55 // used for analysis. 56 57 double h_min, h_max, rho_min, rho_max; 58 size_t h_ngrid, rho_ngrid; 59 60 double time_UtZ; 61 double time_Omega; // Time spent on optimization iterations. 62 double time_Proposal; // Time spent on constructing the 63 // proposal distribution for gamma 64 // (i.e., lmm or lm analysis). 65 66 // Indicator for individuals (phenotypes): 0 missing, 1 67 // available for analysis. 68 vector<int> indicator_idv; 69 70 // Sequence indicator for SNPs: 0 ignored because of (a) maf, 71 // (b) miss, (c) non-poly; 1 available for analysis. 72 vector<int> indicator_snp; 73 74 vector<SNPINFO> snpInfo; // Record SNP information. 75 76 // Main functions. 77 void CopyFromParam(PARAM &cPar); 78 void CopyToParam(PARAM &cPar); 79 80 void WriteResult(const gsl_matrix *Hyper, const gsl_matrix *BF); 81 void WriteResult(const vector<string> &vec_rs, const gsl_matrix *Hyper, 82 const gsl_vector *pip, const gsl_vector *coef); 83 double CalcMarginal(const gsl_vector *Uty, const gsl_vector *K_eval, 84 const double sigma_b2, const double tau); 85 double CalcMarginal(const gsl_matrix *UtXgamma, const gsl_vector *Uty, 86 const gsl_vector *K_eval, const double sigma_a2, 87 const double sigma_b2, const double tau); 88 double CalcPrior(class HYPBSLMM &cHyp); 89 90 void DAP_CalcBF(const gsl_matrix *U, const gsl_matrix *UtX, 91 const gsl_vector *Uty, const gsl_vector *K_eval, 92 const gsl_vector *y); 93 void 94 DAP_EstimateHyper(const size_t kc, const size_t kd, 95 const vector<string> &vec_rs, const vector<double> &vec_sa2, 96 const vector<double> &vec_sb2, const vector<double> &wab, 97 const vector<vector<vector<double>>> &BF, gsl_matrix *Ac, 98 gsl_matrix_int *Ad, gsl_vector_int *dlevel); 99 }; 100 101 void ReadFile_hyb(const string &file_hyp, vector<double> &vec_sa2, 102 vector<double> &vec_sb2, vector<double> &vec_wab); 103 void ReadFile_bf(const string &file_bf, vector<string> &vec_rs, 104 vector<vector<vector<double>>> &BF); 105 void ReadFile_cat(const string &file_cat, const vector<string> &vec_rs, 106 gsl_matrix *Ac, gsl_matrix_int *Ad, gsl_vector_int *dlevel, 107 size_t &kc, size_t &kd); 108 109 #endif 110