1 /*
2 * GCTA: a tool for Genome-wide Complex Trait Analysis
3 *
4 * Implementations of functions for mixed linera model association analysis
5 *
6 * 2013 by Jian Yang <jian.yang@uq.edu.au>
7 *
8 * This file is distributed under the GNU General Public
9 * License, Version 2. Please see the file COPYING for more
10 * details
11 */
12
13 #include "gcta.h"
14
mlma(string grm_file,bool m_grm_flag,string subtract_grm_file,string phen_file,string qcovar_file,string covar_file,int mphen,int MaxIter,vector<double> reml_priors,vector<double> reml_priors_var,bool no_constrain,bool within_family,bool inbred,bool no_adj_covar)15 void gcta::mlma(string grm_file, bool m_grm_flag, string subtract_grm_file, string phen_file, string qcovar_file, string covar_file, int mphen, int MaxIter, vector<double> reml_priors, vector<double> reml_priors_var, bool no_constrain, bool within_family, bool inbred, bool no_adj_covar)
16 {
17 _within_family=within_family;
18 _reml_max_iter=MaxIter;
19 unsigned long i = 0, j = 0, k = 0;
20 bool grm_flag=(!grm_file.empty());
21 bool qcovar_flag=(!qcovar_file.empty());
22 bool covar_flag=(!covar_file.empty());
23 if (!qcovar_flag && !covar_flag) no_adj_covar=false;
24 if (m_grm_flag) grm_flag = false;
25 bool subtract_grm_flag = (!subtract_grm_file.empty());
26 if (subtract_grm_flag && m_grm_flag) throw("Error: the --mlma-subtract-grm option can not be used in combination with the --mgrm option.");
27
28 // Read data
29 int qcovar_num=0, covar_num=0;
30 vector<string> phen_ID, qcovar_ID, covar_ID, grm_id;
31 vector< vector<string> > phen_buf, qcovar, covar; // save individuals by column
32 vector<string> grm_files;
33
34 read_phen(phen_file, phen_ID, phen_buf, mphen);
35 update_id_map_kp(phen_ID, _id_map, _keep);
36 if(qcovar_flag){
37 qcovar_num=read_covar(qcovar_file, qcovar_ID, qcovar, true);
38 update_id_map_kp(qcovar_ID, _id_map, _keep);
39 }
40 if(covar_flag){
41 covar_num=read_covar(covar_file, covar_ID, covar, false);
42 update_id_map_kp(covar_ID, _id_map, _keep);
43 }
44 if(subtract_grm_flag){
45 grm_files.push_back(grm_file);
46 grm_files.push_back(subtract_grm_file);
47 for (i = 0; i < grm_files.size(); i++) {
48 read_grm(grm_files[i], grm_id, false, true, true);
49 update_id_map_kp(grm_id, _id_map, _keep);
50 }
51 }
52 else{
53 if(grm_flag){
54 grm_files.push_back(grm_file);
55 read_grm(grm_file, grm_id, true, false, true);
56 update_id_map_kp(grm_id, _id_map, _keep);
57 }
58 else if (m_grm_flag) {
59 read_grm_filenames(grm_file, grm_files, false);
60 for (i = 0; i < grm_files.size(); i++) {
61 read_grm(grm_files[i], grm_id, false, true, true);
62 update_id_map_kp(grm_id, _id_map, _keep);
63 }
64 }
65 else{
66 grm_files.push_back("NA");
67 make_grm_mkl(false, false, inbred, true, 0, true);
68 for(i=0; i<_keep.size(); i++) grm_id.push_back(_fid[_keep[i]]+":"+_pid[_keep[i]]);
69 }
70 }
71
72 vector<string> uni_id;
73 map<string, int> uni_id_map;
74 map<string, int>::iterator iter;
75 for(i=0; i<_keep.size(); i++){
76 uni_id.push_back(_fid[_keep[i]]+":"+_pid[_keep[i]]);
77 uni_id_map.insert(pair<string,int>(_fid[_keep[i]]+":"+_pid[_keep[i]], i));
78 }
79 _n=_keep.size();
80 if(_n<1) throw("Error: no individual is in common in the input files.");
81 cout<<_n<<" individuals are in common in these files."<<endl;
82
83 // construct model terms
84 _y.setZero(_n);
85 for(i=0; i<phen_ID.size(); i++){
86 iter=uni_id_map.find(phen_ID[i]);
87 if(iter==uni_id_map.end()) continue;
88 _y[iter->second]=atof(phen_buf[i][mphen-1].c_str());
89 }
90
91 _r_indx.clear();
92 vector<int> kp;
93 if (subtract_grm_flag) {
94 for(i=0; i < 2; i++) _r_indx.push_back(i);
95 _A.resize(_r_indx.size());
96
97 cout << "\nReading the primary GRM from [" << grm_files[1] << "] ..." << endl;
98 read_grm(grm_files[1], grm_id, true, false, false);
99
100 StrFunc::match(uni_id, grm_id, kp);
101 (_A[0]).resize(_n, _n);
102 MatrixXf A_N_buf(_n, _n);
103 #pragma omp parallel for private(j)
104 for (j = 0; j < _n; j++) {
105 for (k = 0; k <= j; k++) {
106 if (kp[j] >= kp[k]){
107 (_A[0])(k, j) = (_A[0])(j, k) = _grm(kp[j], kp[k]);
108 A_N_buf(k, j) = A_N_buf(j, k) = _grm_N(kp[j], kp[k]);
109 }
110 else{
111 (_A[0])(k, j) = (_A[0])(j, k) = _grm(kp[k], kp[j]);
112 A_N_buf(k, j) = A_N_buf(j, k) = _grm_N(kp[k], kp[j]);
113 }
114 }
115 }
116
117 cout << "\nReading the secondary GRM from [" << grm_files[0] << "] ..." << endl;
118 read_grm(grm_files[0], grm_id, true, false, false);
119 cout<<"\nSubtracting [" << grm_files[1] << "] from [" << grm_files[0] << "] ..." << endl;
120 StrFunc::match(uni_id, grm_id, kp);
121 #pragma omp parallel for private(j)
122 for (j = 0; j < _n; j++) {
123 for (k = 0; k <= j; k++) {
124 if (kp[j] >= kp[k]) (_A[0])(k, j) = (_A[0])(j, k) = ((_A[0])(j, k) * A_N_buf(j, k) - _grm(kp[j], kp[k]) * _grm_N(kp[j], kp[k])) / (A_N_buf(j, k) - _grm_N(kp[j], kp[k]));
125 else (_A[0])(k, j) = (_A[0])(j, k) = ((_A[0])(j, k) * A_N_buf(j, k) - _grm(kp[k], kp[j]) * _grm_N(kp[k], kp[j])) / (A_N_buf(j, k) - _grm_N(kp[k], kp[j]));
126 }
127 }
128 _grm.resize(0,0);
129 _grm_N.resize(0,0);
130 }
131 else {
132 for(i=0; i < grm_files.size() + 1; i++) _r_indx.push_back(i);
133 _A.resize(_r_indx.size());
134 if(grm_flag){
135 StrFunc::match(uni_id, grm_id, kp);
136 (_A[0]).resize(_n, _n);
137 #pragma omp parallel for private(j)
138 for(i=0; i<_n; i++){
139 for(j=0; j<=i; j++) (_A[0])(j,i)=(_A[0])(i,j)=_grm(kp[i],kp[j]);
140 }
141 _grm.resize(0,0);
142 }
143 else if(m_grm_flag){
144 cout << "There are " << grm_files.size() << " GRM file names specified in the file [" + grm_file + "]." << endl;
145 for (i = 0; i < grm_files.size(); i++) {
146 cout << "Reading the GRM from the " << i + 1 << "th file ..." << endl;
147 read_grm(grm_files[i], grm_id, true, false, true);
148 StrFunc::match(uni_id, grm_id, kp);
149 (_A[i]).resize(_n, _n);
150 #pragma omp parallel for private(j)
151 for (j = 0; j < _n; j++) {
152 for (k = 0; k <= j; k++) {
153 if (kp[j] >= kp[k]) (_A[i])(k, j) = (_A[i])(j, k) = _grm(kp[j], kp[k]);
154 else (_A[i])(k, j) = (_A[i])(j, k) = _grm(kp[k], kp[j]);
155 }
156 }
157 }
158 }
159 else{
160 StrFunc::match(uni_id, grm_id, kp);
161 (_A[0]).resize(_n, _n);
162 #pragma omp parallel for private(j)
163 for(i=0; i<_n; i++){
164 for(j=0; j<=i; j++) (_A[0])(j,i)=(_A[0])(i,j)=_grm_mkl[kp[i]*_n+kp[j]];
165 }
166 delete[] _grm_mkl;
167 }
168 }
169 _A[_r_indx.size()-1]=eigenMatrix::Identity(_n, _n);
170
171 // construct X matrix
172 vector<eigenMatrix> E_float;
173 eigenMatrix qE_float;
174 construct_X(_n, uni_id_map, qcovar_flag, qcovar_num, qcovar_ID, qcovar, covar_flag, covar_num, covar_ID, covar, E_float, qE_float);
175
176 // names of variance component
177 for (i = 0; i < grm_files.size(); i++) {
178 stringstream strstrm;
179 if (grm_files.size() == 1) strstrm << "";
180 else strstrm << i + 1;
181 _var_name.push_back("V(G" + strstrm.str() + ")");
182 _hsq_name.push_back("V(G" + strstrm.str() + ")/Vp");
183 }
184 _var_name.push_back("V(e)");
185
186 // within family
187 if(_within_family) detect_family();
188
189 // run REML algorithm
190 cout << "\nPerforming MLM association analyses" << (subtract_grm_flag?"":" (including the candidate SNP)") << " ..."<<endl;
191 unsigned long n=_keep.size(), m=_include.size();
192 reml(false, true, reml_priors, reml_priors_var, -2.0, -2.0, no_constrain, true, true);
193 _P.resize(0,0);
194 _A.clear();
195 float *y=new float[n];
196 eigenVector y_buf=_y;
197 if(!no_adj_covar) y_buf=_y.array()-(_X*_b).array(); // adjust phenotype for covariates
198 for(i=0; i<n; i++) y[i]=y_buf[i];
199
200 /* if(grm_flag || m_grm_flag){
201 cout<<endl;
202 _geno_mkl=new float[n*m];
203 make_XMat_mkl(_geno_mkl, false);
204 #pragma omp parallel for private(j, k)
205 for(i=0; i<n; i++){
206 for(j=0; j<m; j++){
207 k=i*m+j;
208 if(_geno_mkl[k]<1e5) _geno_mkl[k]-=_mu[_include[j]];
209 else _geno_mkl[k]=0.0;
210 }
211 }
212 }*/
213
214 if (_mu.empty()) calcu_mu();
215 eigenVector beta, se, pval;
216 if(no_adj_covar) mlma_calcu_stat_covar(y, _geno_mkl, n, m, beta, se, pval);
217 else mlma_calcu_stat(y, _geno_mkl, n, m, beta, se, pval);
218 delete[] y;
219 delete[] _geno_mkl;
220
221 string filename=_out+".mlma";
222 cout<<"\nSaving the results of the mixed linear model association analyses of "<<m<<" SNPs to ["+filename+"] ..."<<endl;
223 ofstream ofile(filename.c_str());
224 if(!ofile) throw("Can not open the file ["+filename+"] to write.");
225 ofile<<"Chr\tSNP\tbp\tA1\tA2\tFreq\tb\tse\tp"<<endl;
226 for(i=0; i<m; i++){
227 j=_include[i];
228 ofile<<_chr[j]<<"\t"<<_snp_name[j]<<"\t"<<_bp[j]<<"\t"<<_ref_A[j]<<"\t"<<_other_A[j]<<"\t";
229 if(pval[i]>1.5) ofile<<"NA\tNA\tNA\tNA"<<endl;
230 else ofile<<0.5*_mu[j]<<"\t"<<beta[i]<<"\t"<<se[i]<<"\t"<<pval[i]<<endl;
231 }
232 ofile.close();
233 }
234
mlma_calcu_stat(float * y,float * geno_mkl,unsigned long n,unsigned long m,eigenVector & beta,eigenVector & se,eigenVector & pval)235 void gcta::mlma_calcu_stat(float *y, float *geno_mkl, unsigned long n, unsigned long m, eigenVector &beta, eigenVector &se, eigenVector &pval)
236 {
237 int max_block_size = 10000;
238 unsigned long i=0, j=0;
239 double Xt_Vi_X=0.0, chisq=0.0;
240 float *X=new float[n];
241 float *Vi_X=new float[n];
242 float *Vi=new float[n*n];
243 #pragma omp parallel for private(j)
244 for(i=0; i<n; i++){
245 for(j=0; j<n; j++) Vi[i*n+j]=_Vi(i,j);
246 }
247 _Vi.resize(0,0);
248
249 beta.resize(m);
250 se=eigenVector::Zero(m);
251 pval=eigenVector::Constant(m,2);
252 cout<<"\nRunning association tests for "<<m<<" SNPs ..."<<endl;
253 int new_start = 0, block_size = 0, block_col = 0, k = 0, l = 0;
254 MatrixXf X_block;
255 vector<int> indx;
256 for(i = 0; i < m; i++, block_col++){
257 // get a block of SNPs
258 if(i == new_start){
259 block_col = 0;
260 new_start = i + max_block_size;
261 block_size = max_block_size;
262 if(new_start > m) block_size = m - i;
263 indx.resize(block_size);
264 for(k = i, l = 0; l < block_size; k++, l++) indx[l] = k;
265 make_XMat_subset(X_block, indx, false);
266 }
267
268 for(j = 0; j < n; j++) X[j] = X_block(j, block_col);
269 cblas_sgemv(CblasRowMajor, CblasNoTrans, n, n, 1.0, Vi, n, X, 1, 0.0, Vi_X, 1);
270 Xt_Vi_X=cblas_sdot(n, X, 1, Vi_X, 1);
271 se[i]=1.0/Xt_Vi_X;
272 beta[i]=se[i]*cblas_sdot(n, y, 1, Vi_X, 1);
273 if(se[i]>1.0e-30){
274 se[i]=sqrt(se[i]);
275 chisq=beta[i]/se[i];
276 pval[i]=StatFunc::pchisq(chisq*chisq, 1);
277 }
278 }
279 delete[] X;
280 delete[] Vi_X;
281 delete[] Vi;
282 }
283
mlma_calcu_stat_covar(float * y,float * geno_mkl,unsigned long n,unsigned long m,eigenVector & beta,eigenVector & se,eigenVector & pval)284 void gcta::mlma_calcu_stat_covar(float *y, float *geno_mkl, unsigned long n, unsigned long m, eigenVector &beta, eigenVector &se, eigenVector &pval)
285 {
286 int max_block_size = 10000;
287 unsigned long i=0, j=0, col_num=_X_c+1;
288 double chisq=0.0, d_buf=0.0;
289 float *Vi=new float[n*n];
290 float *X=new float[n*col_num];
291 float *Vi_X=new float[n*col_num];
292 float *Xt_Vi_X=new float[col_num*col_num];
293 float *Xt_Vi_y=new float[col_num];
294 float *b_vec=new float[col_num];
295 #pragma omp parallel for private(j)
296 for(i=0; i<n; i++){
297 for(j=0; j<n; j++) Vi[i*n+j]=_Vi(i,j);
298 }
299 _Vi.resize(0,0);
300 for(i=0; i<n; i++){
301 for(j=0; j<_X_c; j++) X[i*col_num+j]=_X(i,j);
302 X[i*col_num+_X_c]=0.0;
303 }
304
305 beta.resize(m);
306 se=eigenVector::Zero(m);
307 pval=eigenVector::Constant(m,2);
308 cout<<"\nRunning association tests for "<<m<<" SNPs ..."<<endl;
309 int new_start = 0, block_size = 0, block_col = 0, k = 0, l = 0;
310 MatrixXf X_block;
311 vector<int> indx;
312 for(i = 0; i < m; i++, block_col++){
313 // get a block of SNPs
314 if(i == new_start){
315 block_col = 0;
316 new_start = i + max_block_size;
317 block_size = max_block_size;
318 if(new_start > m) block_size = m - i;
319 indx.resize(block_size);
320 for(k = i, l = 0; l < block_size; k++, l++) indx[l] = k;
321 make_XMat_subset(X_block, indx, false);
322 }
323
324 for(j = 0; j < n; j++) X[j*col_num+_X_c] = X_block(j, block_col);
325 cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, n, col_num, n, 1.0, Vi, n, X, col_num, 0.0, Vi_X, col_num);
326 cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, col_num, col_num, n, 1.0, X, col_num, Vi_X, col_num, 0.0, Xt_Vi_X, col_num);
327 if(!comput_inverse_logdet_LU_mkl_array(col_num, Xt_Vi_X, d_buf)) throw("Error: Xt_Vi_X is not invertable.");
328 cblas_sgemv(CblasRowMajor, CblasTrans, n, col_num, 1.0, Vi_X, col_num, y, 1, 0.0, Xt_Vi_y, 1);
329 cblas_sgemv(CblasRowMajor, CblasNoTrans, col_num, col_num, 1.0, Xt_Vi_X, col_num, Xt_Vi_y, 1, 0.0, b_vec, 1);
330 se[i]=Xt_Vi_X[_X_c*col_num+_X_c];
331 beta[i]=b_vec[_X_c];
332 if(se[i]>1.0e-30){
333 se[i]=sqrt(se[i]);
334 chisq=beta[i]/se[i];
335 pval[i]=StatFunc::pchisq(chisq*chisq, 1);
336 }
337 }
338 delete[] Vi;
339 delete[] X;
340 delete[] Vi_X;
341 delete[] Xt_Vi_X;
342 delete[] Xt_Vi_y;
343 delete[] b_vec;
344 }
345
mlma_loco(string phen_file,string qcovar_file,string covar_file,int mphen,int MaxIter,vector<double> reml_priors,vector<double> reml_priors_var,bool no_constrain,bool inbred,bool no_adj_covar)346 void gcta::mlma_loco(string phen_file, string qcovar_file, string covar_file, int mphen, int MaxIter, vector<double> reml_priors, vector<double> reml_priors_var, bool no_constrain, bool inbred, bool no_adj_covar)
347 {
348 unsigned long i=0, j=0, k=0, c1=0, c2=0, n=0;
349 _reml_max_iter=MaxIter;
350 bool qcovar_flag=(!qcovar_file.empty());
351 bool covar_flag=(!covar_file.empty());
352 if(!qcovar_flag && !covar_flag) no_adj_covar=false;
353
354 // Read data
355 int qcovar_num=0, covar_num=0;
356 vector<string> phen_ID, qcovar_ID, covar_ID, grm_id;
357 vector< vector<string> > phen_buf, qcovar, covar; // save individuals by column
358
359 read_phen(phen_file, phen_ID, phen_buf, mphen);
360 update_id_map_kp(phen_ID, _id_map, _keep);
361 if(qcovar_flag){
362 qcovar_num=read_covar(qcovar_file, qcovar_ID, qcovar, true);
363 update_id_map_kp(qcovar_ID, _id_map, _keep);
364 }
365 if(covar_flag){
366 covar_num=read_covar(covar_file, covar_ID, covar, false);
367 update_id_map_kp(covar_ID, _id_map, _keep);
368 }
369 n=_keep.size();
370 _n=_keep.size();
371 if(_n<1) throw("Error: no individual is in common in the input files.");
372 cout<<_n<<" individuals are in common in these files."<<endl;
373
374 vector<int> chrs, vi_buf(_chr);
375 stable_sort(vi_buf.begin(), vi_buf.end());
376 vi_buf.erase(unique(vi_buf.begin(), vi_buf.end()), vi_buf.end());
377 if(vi_buf.size()<2) throw("Error: There is only one chromosome. The MLM leave-on-chromosome-out (LOCO) analysis requires at least two chromosome.");
378 for(i=0; i<vi_buf.size(); i++){
379 if(vi_buf[i]<=_autosome_num) chrs.push_back(vi_buf[i]);
380 }
381 vector<int> include_o(_include);
382 map<string, int> snp_name_map_o(_snp_name_map);
383 vector<float> m_chrs_f(chrs.size());
384 vector<float *> grm_chrs(chrs.size());
385 vector<float *> geno_chrs(chrs.size());
386 vector< vector<int> > icld_chrs(chrs.size());
387 cout<<endl;
388 if(_mu.empty()) calcu_mu();
389 cout<<"\nCalculating the genetic relationship matrix for each of the "<<chrs.size()<<" chromosomes ... "<<endl;
390 for(c1=0; c1<chrs.size(); c1++){
391 cout<<"Chr "<<chrs[c1]<<":"<<endl;
392 extract_chr(chrs[c1], chrs[c1]);
393 make_grm_mkl(false, false, inbred, true, 0, true);
394
395 m_chrs_f[c1]=(float)_include.size();
396 icld_chrs[c1]=_include;
397 _include=include_o;
398 _snp_name_map=snp_name_map_o;
399
400 geno_chrs[c1]=_geno_mkl;
401 _geno_mkl=NULL;
402 grm_chrs[c1]=_grm_mkl;
403 _grm_mkl=NULL;
404 }
405 for(i=0; i<_keep.size(); i++) grm_id.push_back(_fid[_keep[i]]+":"+_pid[_keep[i]]);
406
407 vector<string> uni_id;
408 map<string, int> uni_id_map;
409 map<string, int>::iterator iter;
410 for(i=0; i<_keep.size(); i++){
411 uni_id.push_back(_fid[_keep[i]]+":"+_pid[_keep[i]]);
412 uni_id_map.insert(pair<string,int>(_fid[_keep[i]]+":"+_pid[_keep[i]], i));
413 }
414
415 // construct model terms
416 _y.setZero(_n);
417 for(i=0; i<phen_ID.size(); i++){
418 iter=uni_id_map.find(phen_ID[i]);
419 if(iter==uni_id_map.end()) continue;
420 _y[iter->second]=atof(phen_buf[i][mphen-1].c_str());
421 }
422
423 // construct X matrix
424 vector<eigenMatrix> E_float;
425 eigenMatrix qE_float;
426 construct_X(_n, uni_id_map, qcovar_flag, qcovar_num, qcovar_ID, qcovar, covar_flag, covar_num, covar_ID, covar, E_float, qE_float);
427
428 // names of variance component
429 _var_name.push_back("V(G)");
430 _hsq_name.push_back("V(G)/Vp");
431 _var_name.push_back("V(e)");
432
433 // MLM association
434 cout<<"\nPerforming MLM association analyses (leave-one-chromosome-out) ..."<<endl;
435
436 vector<int> kp;
437 StrFunc::match(uni_id, grm_id, kp);
438 _r_indx.resize(2);
439 for(i=0; i<2; i++) _r_indx[i]=i;
440 _A.resize(_r_indx.size());
441 _A[1]=eigenMatrix::Identity(_n, _n);
442
443 eigenVector y_buf=_y;
444 float *y=new float[_n];
445 vector<eigenVector> beta(chrs.size()), se(chrs.size()), pval(chrs.size());
446 for(c1=0; c1<chrs.size(); c1++){
447 cout<<"\n-----------------------------------\n#Chr "<<chrs[c1]<<":"<<endl;
448 extract_chr(chrs[c1], chrs[c1]);
449
450 _A[0]=eigenMatrix::Zero(_n, _n);
451 double d_buf=0;
452 for(c2=0; c2<chrs.size(); c2++){
453 if(chrs[c1]==chrs[c2]) continue;
454 #pragma omp parallel for private(j)
455 for(i=0; i<_n; i++){
456 for(j=0; j<=i; j++){
457 (_A[0])(i,j)+=(grm_chrs[c2])[kp[i]*_n+kp[j]]*m_chrs_f[c2];
458 }
459 }
460 d_buf+=m_chrs_f[c2];
461 }
462
463 #pragma omp parallel for private(j)
464 for(i=0; i<_n; i++){
465 for(j=0; j<=i; j++){
466 (_A[0])(i,j)/=d_buf;
467 (_A[0])(j,i)=(_A[0])(i,j);
468 }
469 }
470
471 // run REML algorithm
472 reml(false, true, reml_priors, reml_priors_var, -2.0, -2.0, no_constrain, true, true);
473 if(!no_adj_covar) y_buf=_y.array()-(_X*_b).array(); // adjust phenotype for covariates
474 for(i=0; i<_n; i++) y[i]=y_buf[i];
475 reml_priors.clear();
476 reml_priors_var=_varcmp;
477 _P.resize(0,0);
478 _A[0].resize(0,0);
479
480 mlma_calcu_stat(y, (geno_chrs[c1]), n, _include.size(), beta[c1], se[c1], pval[c1]);
481
482 _include=include_o;
483 _snp_name_map=snp_name_map_o;
484 cout<<"-----------------------------------"<<endl;
485 }
486
487 delete[] y;
488 for(c1=0; c1<chrs.size(); c1++){
489 delete[] (grm_chrs[c1]);
490 delete[] (geno_chrs[c1]);
491 }
492
493 string filename=_out+".loco.mlma";
494 cout<<"\nSaving the results of the mixed linear model association analyses of "<<_include.size()<<" SNPs to ["+filename+"] ..."<<endl;
495 ofstream ofile(filename.c_str());
496 if(!ofile) throw("Can not open the file ["+filename+"] to write.");
497 ofile<<"Chr\tSNP\tbp\tA1\tA2\tFreq\tb\tse\tp"<<endl;
498 for(c1=0; c1<chrs.size(); c1++){
499 for(i=0; i<icld_chrs[c1].size(); i++){
500 j=icld_chrs[c1][i];
501 ofile<<_chr[j]<<"\t"<<_snp_name[j]<<"\t"<<_bp[j]<<"\t"<<_ref_A[j]<<"\t"<<_other_A[j]<<"\t";
502 if(pval[c1][i]>1.5) ofile<<"NA\tNA\tNA\tNA"<<endl;
503 else ofile<<0.5*_mu[j]<<"\t"<<beta[c1][i]<<"\t"<<se[c1][i]<<"\t"<<pval[c1][i]<<endl;
504 }
505 }
506 ofile.close();
507 }
508
grm_minus_grm(float * grm,float * sub_grm)509 void gcta::grm_minus_grm(float *grm, float *sub_grm)
510 {
511 int i=0, j=0, k=0, n=_n;
512
513 #pragma omp parallel for private(j,k)
514 for(i=0; i<n; i++){
515 for(j=0; j<=i; j++){
516 k=i*n+j;
517 sub_grm[k]=grm[k]-sub_grm[k];
518 }
519 }
520
521 }
522
523
524