Adding loco support for mdb

4 files changed, 14 insertions, 10 deletions

diff --git a/src/gemma.cpp b/src/gemma.cpp
index 94d61a8..a39d67a 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -2843,7 +2843,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
           write(&Y_col.vector, "Y_col");
 
           if (cPar.is_mdb) {
-            cLmm.mdb_calc_gwa(U, eval, UtW, &UtY_col.vector, W, &Y_col.vector);
+            cLmm.mdb_calc_gwa(U, eval, UtW, &UtY_col.vector, W, &Y_col.vector, cPar.loco);
             cLmm.CopyToParam(cPar);
           }
           else {
diff --git a/src/lmm.cpp b/src/lmm.cpp
index 375013c..7a80b8a 100644
--- a/src/lmm.cpp
+++ b/src/lmm.cpp
@@ -2155,8 +2155,10 @@ void LMM::mdb_analyze(std::function< SnpNameValues2(size_t) >& fetch_snp,
 
 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 gsl_matrix *W, const gsl_vector *y, const string loco) {
   checkpoint("mdb-calc-gwa",file_geno);
+  bool is_loco = !loco.empty();
+  cout << loco << " CHR!!!!" << endl;
 
   // const auto num_snps = indicator_snp.size();
   // enforce_msg(num_snps > 0,"Zero SNPs to process - data corrupt?");
diff --git a/src/lmm.h b/src/lmm.h
index f10996d..f4f49aa 100644
--- a/src/lmm.h
+++ b/src/lmm.h
@@ -131,8 +131,9 @@ public:
                    const gsl_matrix *UtW, const gsl_vector *Uty,
                    const gsl_matrix *W, const gsl_vector *y, size_t num_markers);
   void 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 gsl_matrix *UtW, const gsl_vector *Uty,
+                    const gsl_matrix *W, const gsl_vector *y,
+                    const string loco);
   void AnalyzeBimbam(const gsl_matrix *U, const gsl_vector *eval,
                      const gsl_matrix *UtW, const gsl_vector *Uty,
                      const gsl_matrix *W, const gsl_vector *y,
diff --git a/src/param.cpp b/src/param.cpp
index 74b2f9b..516dfd8 100644
--- a/src/param.cpp
+++ b/src/param.cpp
@@ -149,10 +149,6 @@ void PARAM::ReadFiles(void) {
     }
   }
 
-  if (file_geno.find(".mdb") != string::npos) {
-    is_mdb = true;
-  }
-
   if (!is_mdb) {
     // Read SNP set into setSnps (without filtering)
     if (!file_snps.empty()) {
@@ -943,13 +939,18 @@ void PARAM::CheckParam(void) {
   enforce_fexists(file_gwasnps, "open file");
   enforce_fexists(file_anno, "open file");
 
+  if (file_geno.find(".mdb") != string::npos) {
+    is_mdb = true;
+  }
+
   if (!loco.empty()) {
     enforce_msg((a_mode >= 1 && a_mode <= 4) || a_mode == 9 || a_mode == 21 || a_mode == 22,
                 "LOCO only works with LMM and K");
     // enforce_msg(file_bfile.empty(), "LOCO does not work with PLink (yet)");
     enforce_msg(file_gxe.empty(), "LOCO does not support GXE (yet)");
-    enforce_msg(!file_anno.empty(),
-                "LOCO requires annotation file (-a switch)");
+    if (!is_mdb)
+      enforce_msg(!file_anno.empty(), "Without mdb LOCO requires annotation file (-a switch)");
+
     enforce_msg(file_ksnps.empty(), "LOCO does not allow -ksnps switch");
     enforce_msg(file_gwasnps.empty(), "LOCO does not allow -gwasnps switch");
   }