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