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 __VC_H__ 20 #define __VC_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 VC_PARAM { 30 31 public: 32 const gsl_matrix *K; 33 const gsl_matrix *W; 34 const gsl_vector *y; 35 gsl_matrix *P; 36 gsl_vector *Py; 37 gsl_matrix *KPy_mat; 38 gsl_matrix *PKPy_mat; 39 gsl_matrix *Hessian; 40 bool noconstrain; 41 }; 42 43 // Methods for variant component (VC) estimation 44 45 class VC { 46 47 public: 48 // IO-related parameters 49 size_t a_mode; 50 string file_cat; 51 string file_beta; 52 string file_cor; 53 string file_mq; 54 string file_ms; 55 56 string file_out; 57 string path_out; 58 59 set<string> setSnps; 60 61 size_t ni_total_ref, ns_total_ref, ns_pair; 62 size_t ni_total, ns_total, ns_test; 63 size_t n_vc; 64 65 double pve_total, se_pve_total; 66 vector<double> v_sigma2; 67 vector<double> v_se_sigma2; 68 vector<double> v_pve; 69 vector<double> v_se_pve; 70 vector<double> v_traceG; 71 vector<double> v_beta; 72 vector<double> v_se_beta; 73 74 size_t crt; 75 double window_cm, window_bp, window_ns; 76 77 double time_UtX; 78 double time_opt; 79 80 // Main functions. 81 void CopyFromParam(PARAM &cPar); 82 void CopyToParam(PARAM &cPar); 83 void WriteFile_qs(const gsl_vector *s_vec, const gsl_vector *q_vec, 84 const gsl_vector *qvar_vec, const gsl_matrix *S_mat, 85 const gsl_matrix *Svar_mat); 86 void CalcVChe(const gsl_matrix *K, const gsl_matrix *W, const gsl_vector *y); 87 void CalcVCreml(const bool noconstrain, const gsl_matrix *K, 88 const gsl_matrix *W, const gsl_vector *y); 89 void CalcVCacl(const gsl_matrix *K, const gsl_matrix *W, const gsl_vector *y); 90 }; 91 92 void CalcVCss(const gsl_matrix *Vq, const gsl_matrix *S_mat, 93 const gsl_matrix *Svar_mat, const gsl_vector *q_vec, 94 const gsl_vector *s_vec, const double df, vector<double> &v_pve, 95 vector<double> &v_se_pve, double &pve_total, double &se_pve_total, 96 vector<double> &v_sigma2, vector<double> &v_se_sigma2, 97 vector<double> &v_enrich, vector<double> &v_se_enrich); 98 99 bool BimbamXwz(const string &file_geno, const int display_pace, 100 vector<int> &indicator_idv, vector<int> &indicator_snp, 101 const vector<size_t> &vec_cat, const gsl_vector *w, 102 const gsl_vector *z, size_t ns_test, gsl_matrix *XWz); 103 bool PlinkXwz(const string &file_bed, const int display_pace, 104 vector<int> &indicator_idv, vector<int> &indicator_snp, 105 const vector<size_t> &vec_cat, const gsl_vector *w, 106 const gsl_vector *z, size_t ns_test, gsl_matrix *XWz); 107 bool MFILEXwz(const size_t mfile_mode, const string &file_mfile, 108 const int display_pace, vector<int> &indicator_idv, 109 vector<vector<int>> &mindicator_snp, 110 const vector<size_t> &vec_cat, const gsl_vector *w, 111 const gsl_vector *z, gsl_matrix *XWz); 112 113 bool BimbamXtXwz(const string &file_geno, const int display_pace, 114 vector<int> &indicator_idv, vector<int> &indicator_snp, 115 const gsl_matrix *XWz, size_t ns_test, gsl_matrix *XtXWz); 116 bool PlinkXtXwz(const string &file_bed, const int display_pace, 117 vector<int> &indicator_idv, vector<int> &indicator_snp, 118 const gsl_matrix *XWz, size_t ns_test, gsl_matrix *XtXWz); 119 bool MFILEXtXwz(const size_t mfile_mode, const string &file_mfile, 120 const int display_pace, vector<int> &indicator_idv, 121 vector<vector<int>> &mindicator_snp, const gsl_matrix *XWz, 122 gsl_matrix *XtXWz); 123 124 void CalcCIss(const gsl_matrix *Xz, const gsl_matrix *XWz, 125 const gsl_matrix *XtXWz, const gsl_matrix *S_mat, 126 const gsl_matrix *Svar_mat, const gsl_vector *w, 127 const gsl_vector *z, const gsl_vector *s_vec, 128 const vector<size_t> &vec_cat, const vector<double> &v_pve, 129 vector<double> &v_se_pve, double &pve_total, double &se_pve_total, 130 vector<double> &v_sigma2, vector<double> &v_se_sigma2, 131 vector<double> &v_enrich, vector<double> &v_se_enrich); 132 133 #endif 134