diff options
Diffstat (limited to 'src/lmm.cpp')
| -rw-r--r-- | src/lmm.cpp | 50 |
1 files changed, 18 insertions, 32 deletions
diff --git a/src/lmm.cpp b/src/lmm.cpp index 129cbd6..2b730f3 100644 --- a/src/lmm.cpp +++ b/src/lmm.cpp @@ -2012,23 +2012,17 @@ void LMM::mdb_analyze(std::function< SnpNameValues2(size_t) >& fetch_snp, // continue; auto tup = fetch_snp(t); // use the callback - auto success = get<0>(tup); - if (!success) + auto state = get<0>(tup); + if (state == SKIP) continue; - // typedef tuple< bool,MarkerInfo,vector<double> > SnpNameValues2; - // auto marker = get<1>(tup); - // auto chr = get<2>(tup); - // auto mpos = get<3>(tup); + if (state == LAST) + break; // marker loop because of LOCO + auto markerinfo = get<1>(tup); auto gs = get<2>(tup); markers.push_back(markerinfo); - // check whether SNP is included in gwasnps (used by LOCO) - /* - if (process_gwasnps && gwasnps.count(snp) == 0) - continue; - */ // 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 @@ -2067,7 +2061,6 @@ void LMM::mdb_analyze(std::function< SnpNameValues2(size_t) >& fetch_snp, gsl_vector_safe_memcpy(&Xlarge_col.vector, x); c++; // count markers going in - if (c % msize == 0) { batch_compute(msize,markers); markers.clear(); @@ -2212,10 +2205,6 @@ void LMM::mdb_calc_gwa(const gsl_matrix *U, const gsl_vector *eval, mdb_stat(rtxn, geno_mdb, &stat); auto num_markers = stat.ms_entries; - // fetch_snp is a callback function for every SNP row - // returns typedef std::tuple<bool,string,std::vector<double> > SnpNameValues; - - // size_t prev_line = 0; auto mdb_fetch = MDB_FIRST; auto cursor = lmdb::cursor::open(rtxn, geno_mdb); @@ -2227,21 +2216,16 @@ void LMM::mdb_calc_gwa(const gsl_matrix *U, const gsl_vector *eval, string_view key,value; - /* - while (prev_line <= num) { - cursor.get(key, value, MDB_NEXT); - prev_line++; - } - */ - auto success = cursor.get(key, value, mdb_fetch); + auto mdb_success = cursor.get(key, value, mdb_fetch); mdb_fetch = MDB_NEXT; // uint8_t chr; vector<double> gs; MarkerInfo markerinfo; - if (success) { + if (mdb_success) { size_t size = 0; + // ---- Depending on the format we get different buffers - currently float and byte are supported: if (format == "Gb") { size_t num_bytes = value.size() / sizeof(uint8_t); assert(num_bytes == ni_total); @@ -2270,13 +2254,6 @@ void LMM::mdb_calc_gwa(const gsl_matrix *U, const gsl_vector *eval, // "S>L>L>" const uint8_t* data = reinterpret_cast<const uint8_t*>(key.data()); auto chr = static_cast<uint8_t>(data[1]); - - // printf("%#02x %#02x\n", chr, loco_chr); - - if (is_loco && loco_chr != chr) { - return make_tuple(false, markerinfo, gs); - } - // Extract big-endian uint32 // uint32_t rest = static_cast<uint32_t>(data[2]); uint32_t pos = (data[2] << 24) | (data[3] << 16) | @@ -2285,6 +2262,15 @@ void LMM::mdb_calc_gwa(const gsl_matrix *U, const gsl_vector *eval, uint32_t num = (data[6] << 24) | (data[7] << 16) | (data[8] << 8) | data[9]; + // printf("%#02x %#02x\n", chr, loco_chr); + + if (is_loco && loco_chr != chr) { + if (chr > loco_chr) + return make_tuple(LAST, MarkerInfo { .name="", .chr=chr, .pos=pos } , gs); + else + return make_tuple(SKIP, MarkerInfo { .name="", .chr=chr, .pos=pos } , gs); + } + string_view value2; marker_mdb.get(rtxn,key,value2); auto marker = string(value2); @@ -2299,7 +2285,7 @@ void LMM::mdb_calc_gwa(const gsl_matrix *U, const gsl_vector *eval, // cout << "!!!!" << size << marker << ": af" << maf << " " << gs[0] << "," << gs[1] << "," << gs[2] << "," << gs[3] << endl; } - return make_tuple(success, markerinfo, gs); + return make_tuple(COMPUTE, markerinfo, gs); }; LMM::mdb_analyze(fetch_snp,U,eval,UtW,Uty,W,y,num_markers); |
