Unpack key to get chr,pos,marker

1 files changed, 25 insertions, 7 deletions

diff --git a/src/lmm.cpp b/src/lmm.cpp
index 761eacf..da3e100 100644
--- a/src/lmm.cpp
+++ b/src/lmm.cpp
@@ -2095,7 +2095,6 @@ void LMM::mdb_calc_gwa(const gsl_matrix *U, const gsl_vector *eval,
 
   // const auto num_snps = indicator_snp.size();
   // enforce_msg(num_snps > 0,"Zero SNPs to process - data corrupt?");
-
   auto env = lmdb::env::create();
   env.set_mapsize(1UL * 1024UL * 1024UL * 1024UL * 1024UL); /* 10 GiB */
   env.set_max_dbs(10);
@@ -2112,9 +2111,10 @@ void LMM::mdb_calc_gwa(const gsl_matrix *U, const gsl_vector *eval,
   auto rtxn = lmdb::txn::begin(env, nullptr, MDB_RDONLY);
   auto geno_mdb = lmdb::dbi::open(rtxn, "geno");
 
+  auto marker_mdb = lmdb::dbi::open(rtxn, "marker");
+
   MDB_stat stat;
   mdb_stat(rtxn, geno_mdb, &stat);
-  // cout << "Number of records: " << stat.ms_entries << endl;
   auto num_markers = stat.ms_entries;
 
   // fetch_snp is a callback function for every SNP row
@@ -2141,15 +2141,34 @@ void LMM::mdb_calc_gwa(const gsl_matrix *U, const gsl_vector *eval,
     auto success = cursor.get(key, value, mdb_fetch);
     mdb_fetch = MDB_NEXT;
 
-    string snp;
+    string marker;
     vector<double> gs;
 
     if (success) {
       size_t num_floats = value.size() / sizeof(float);
       assert(num_floats == ni_total);
       const float* gsbuf = reinterpret_cast<const float*>(value.data());
-      auto snp = string(key);
-      // std::vector<double> gs(gsbuf,gsbuf+sizeof(float)*value.size());
+
+      // unpack chr+pos from key
+      if (key.size() != 5) {
+        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]);
+
+      // Extract big-endian uint32 manually
+      auto 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]));
+
+      string_view value2;
+      marker_mdb.get(rtxn,key,value2);
+      marker = string(value2);
+      // cout << static_cast<int>(chr) << ":" << pos << " " << marker << endl ;
+
       auto size = num_floats;
       gs.reserve(size);
       for (size_t i = 0; i < size; ++i) {
@@ -2157,9 +2176,8 @@ 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, snp, gs);
+    return make_tuple(success, marker, gs);
   };
-
   LMM::mdb_analyze(fetch_snp,U,eval,UtW,Uty,W,y,gwasnps,num_markers);
 }