1 files changed, 54 insertions, 0 deletions
diff --git a/src/lmm.cpp b/src/lmm.cpp
index f5b8a88..e58adc9 100644
--- a/src/lmm.cpp
+++ b/src/lmm.cpp
@@ -1867,6 +1867,60 @@ void LMM::Analyze(std::function< SnpNameValues(size_t) >& fetch_snp,
 }
 
 /*
+  This is the mirror function of below AnalyzeBimbam, but uses mdb input instead.
+ */
+
+void LMM::mdb_calc_gwa(const gsl_matrix *U, const gsl_vector *eval,
+                          const gsl_matrix *UtW, const gsl_vector *Uty,
+                          const gsl_matrix *W, const gsl_vector *y,
+                          const set<string> gwasnps) {
+  debug_msg(file_geno);
+  auto infilen = file_geno.c_str();
+  checkpoint("start-read-geno-file",file_geno);
+
+  igzstream infile(infilen, igzstream::in);
+  enforce_msg(infile, "error reading genotype file");
+  size_t prev_line = 0;
+
+  std::vector <double> gs;
+  gs.resize(ni_total);
+
+  // fetch_snp is a callback function for every SNP row
+  std::function<SnpNameValues(size_t)>  fetch_snp = [&](size_t num) {
+    string line;
+    while (prev_line <= num) {
+      // also read SNPs that were skipped
+      safeGetline(infile, line);
+      prev_line++;
+    }
+    char *ch_ptr = strtok_safe2((char *)line.c_str(), " , \t",infilen);
+    // enforce_msg(ch_ptr, "Parsing BIMBAM genofile"); // ch_ptr should not be NULL
+
+    auto snp = string(ch_ptr);
+    ch_ptr = strtok_safe2(NULL, " , \t",infilen); // skip column
+    ch_ptr = strtok_safe2(NULL, " , \t",infilen); // skip column
+
+    gs.assign (ni_total,nan("")); // wipe values
+
+    for (size_t i = 0; i < ni_total; ++i) {
+      ch_ptr = strtok_safe2(NULL, " , \t",infilen);
+      if (strcmp(ch_ptr, "NA") != 0) {
+        gs[i] = atof(ch_ptr);
+        if (is_strict_mode() && gs[i] == 0.0)
+          enforce_is_float(std::string(ch_ptr)); // only allow for NA and valid numbers
+      }
+    }
+    return std::make_tuple(snp,gs);
+  };
+
+  LMM::Analyze(fetch_snp,U,eval,UtW,Uty,W,y,gwasnps);
+
+  infile.close();
+  infile.clear();
+}
+
+
+/*
 
 Looking at the `LMM::AnalyzeBimbam` function, here are the parameters that get modified: