diff options
| -rw-r--r-- | src/lmm.cpp | 31 | ||||
| -rw-r--r-- | src/lmm.h | 3 |
2 files changed, 18 insertions, 16 deletions
diff --git a/src/lmm.cpp b/src/lmm.cpp index da3e100..0a5d8e6 100644 --- a/src/lmm.cpp +++ b/src/lmm.cpp @@ -2014,8 +2014,10 @@ void LMM::mdb_analyze(std::function< SnpNameValues2(size_t) >& fetch_snp, auto success = get<0>(tup); if (!success) continue; - auto snp = get<1>(tup); - auto gs = get<2>(tup); + auto marker = get<1>(tup); + auto chr = get<2>(tup); + auto mpos = get<3>(tup); + auto gs = get<4>(tup); // cout << t << " SNP: " << snp << endl; // check whether SNP is included in gwasnps (used by LOCO) @@ -2025,7 +2027,7 @@ 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 pos = 0; // position in target vector + uint vpos = 0; // position in target vector uint n_miss = 0; gsl_vector_set_zero(x_miss); for (size_t i = 0; i < ni_total; ++i) { @@ -2035,15 +2037,15 @@ void LMM::mdb_analyze(std::function< SnpNameValues2(size_t) >& fetch_snp, double geno = gs[i]; if (isnan(geno)) { - gsl_vector_set(x_miss, pos, 1.0); + gsl_vector_set(x_miss, vpos, 1.0); n_miss++; } else { - gsl_vector_set(x, pos, geno); + gsl_vector_set(x, vpos, geno); x_total += geno; } - pos++; + vpos++; } - enforce(pos == ni_test); + enforce(vpos == ni_test); const double x_mean = x_total/(double)(ni_test - n_miss); @@ -2068,9 +2070,9 @@ void LMM::mdb_analyze(std::function< SnpNameValues2(size_t) >& fetch_snp, if (c % msize) batch_compute(c % msize); - ProgressBar("Reading SNPS", num_markers - 1, num_markers - 1); + ProgressBar("Reading markers", num_markers - 1, num_markers - 1); cout << endl; - cout << "Counted SNPs " << c << " sumStat " << sumStat.size() << endl; + cout << "Counted markers " << c << " sumStat " << sumStat.size() << endl; checkpoint_nofile("end-lmm-mdb-analyze"); gsl_vector_safe_free(x); @@ -2142,6 +2144,8 @@ void LMM::mdb_calc_gwa(const gsl_matrix *U, const gsl_vector *eval, mdb_fetch = MDB_NEXT; string marker; + uint8_t chr; + size_t pos; vector<double> gs; if (success) { @@ -2154,12 +2158,9 @@ void LMM::mdb_calc_gwa(const gsl_matrix *U, const gsl_vector *eval, throw std::runtime_error("Invalid key size"); } const uint8_t* data = reinterpret_cast<const uint8_t*>(key.data()); - - // Extract signed char - auto chr = static_cast<uint8_t>(data[0]); - + chr = static_cast<uint8_t>(data[0]); // Extract big-endian uint32 manually - auto pos = (static_cast<uint32_t>(data[1]) << 24) | + pos = (static_cast<uint32_t>(data[1]) << 24) | (static_cast<uint32_t>(data[2]) << 16) | (static_cast<uint32_t>(data[3]) << 8) | (static_cast<uint32_t>(data[4])); @@ -2176,7 +2177,7 @@ void LMM::mdb_calc_gwa(const gsl_matrix *U, const gsl_vector *eval, } // cout << "!!!!" << size << snp << ": " << gs[0] << "," << gs[1] << "," << gs[2] << "," << gs[3] << endl; } - return make_tuple(success, marker, gs); + return make_tuple(success, chr, pos, marker, gs); }; LMM::mdb_analyze(fetch_snp,U,eval,UtW,Uty,W,y,gwasnps,num_markers); } diff --git a/src/lmm.h b/src/lmm.h index 75ed6de..65d37b2 100644 --- a/src/lmm.h +++ b/src/lmm.h @@ -27,6 +27,7 @@ #include "param.h" #include <functional> #include <tuple> +#include <cstdint> using namespace std; @@ -45,7 +46,7 @@ public: }; typedef tuple< string,vector<double> > SnpNameValues; -typedef tuple< bool,string,vector<double> > SnpNameValues2; +typedef tuple< bool,uint8_t,size_t,string,vector<double> > SnpNameValues2; class LMM { |