diff options
Diffstat (limited to 'src/lmm.cpp')
| -rw-r--r-- | src/lmm.cpp | 24 |
1 files changed, 18 insertions, 6 deletions
diff --git a/src/lmm.cpp b/src/lmm.cpp index d13ce3c..66b8e37 100644 --- a/src/lmm.cpp +++ b/src/lmm.cpp @@ -96,8 +96,9 @@ void LMM::CopyToParam(PARAM &cPar) { checkpoint_nofile("lmm-copy-to-param"); cPar.time_UtX = time_UtX; cPar.time_opt = time_opt; - - cPar.ng_test = ng_test; // number of markers tested + cPar.ns_total = ns_total; + cPar.ns_test = ns_test; + cPar.ng_test = ng_test; // number of markers tested return; } @@ -2034,8 +2035,8 @@ void LMM::mdb_analyze(std::function< SnpNameValues2(size_t) >& fetch_snp, */ // drop missing idv and plug mean values for missing geno double x_total = 0.0; // sum genotype values to compute x_mean - uint vpos = 0; // position in target vector - uint n_miss = 0; + uint vpos = 0; // position in target vector + uint n_miss = 0; // count NA genotypes gsl_vector_set_zero(x_miss); for (size_t i = 0; i < ni_total; ++i) { // get the genotypes per individual and compute stats per SNP @@ -2107,6 +2108,9 @@ void LMM::mdb_analyze(std::function< SnpNameValues2(size_t) >& fetch_snp, outfile << chr << "\t"; outfile << name << "\t"; outfile << pos << "\t"; + auto n_miss = 0; + outfile << n_miss << "\t-\t-\t"; // n_miss column + // outfile << st.maf << "\t" ; if (a_mode != M_LMM2) { outfile << st.beta << "\t"; outfile << st.se << "\t"; @@ -2240,18 +2244,26 @@ void LMM::mdb_calc_gwa(const gsl_matrix *U, const gsl_vector *eval, auto marker = string(value2); // 1 rs13476251 174792257 // cout << static_cast<int>(chr) << ":" << pos2 << " line " << rest2 << ":" << marker << endl ; - markerinfo = make_tuple(marker,chr,pos,num); auto size = num_floats; gs.reserve(size); for (size_t i = 0; i < size; ++i) { gs.push_back(static_cast<double>(gsbuf[i])); } - // cout << "!!!!" << size << snp << ": " << gs[0] << "," << gs[1] << "," << gs[2] << "," << gs[3] << endl; + + // compute maf and n_miss (NAs) + size_t n_miss = 0; // count NAs: FIXME + double maf = compute_maf(ni_total, ni_test, n_miss, gs.data()); + + markerinfo = make_tuple(marker,chr,pos,num); + + // cout << "!!!!" << size << marker << ": af" << maf << " " << gs[0] << "," << gs[1] << "," << gs[2] << "," << gs[3] << endl; } return make_tuple(success, markerinfo, gs); }; LMM::mdb_analyze(fetch_snp,U,eval,UtW,Uty,W,y,gwasnps,num_markers); + + ns_total = ns_test = num_markers; // track global number of snps in original gemma (goes to cPar) } |