Wiring in mdb GWA

1 files changed, 37 insertions, 29 deletions

diff --git a/src/gemma.cpp b/src/gemma.cpp
index 1cced77..fb2227e 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -2591,20 +2591,24 @@ void GEMMA::BatchRun(PARAM &cPar) {
       cPar.a_mode == M_LMM4 || cPar.a_mode == M_LMM5 || cPar.a_mode == M_LMM9 ||
       cPar.a_mode == M_EIGEN) { // Fit LMM or mvLMM or eigen
     write(cPar.a_mode, "Mode");
-    gsl_matrix *Y = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_ph); // Y is ind x phenotypes
+    auto ni_test = cPar.ni_test;
+    auto n_ph = cPar.n_ph;
+    enforce(ni_test > 0);
+    enforce(n_ph > 0);
+    gsl_matrix *Y = gsl_matrix_safe_alloc(ni_test, n_ph); // Y is ind x phenotypes
     enforce_msg(Y, "allocate Y"); // just to be sure
-    gsl_matrix *W = gsl_matrix_safe_alloc(Y->size1, cPar.n_cvt);    // W is ind x covariates
+    gsl_matrix *W = gsl_matrix_safe_alloc(ni_test, cPar.n_cvt);    // W is ind x covariates
     gsl_matrix *B = gsl_matrix_safe_alloc(Y->size2, W->size2);      // B is a d by c (effect size)
                                                           // matrix
     gsl_matrix *se_B = gsl_matrix_safe_alloc(Y->size2, W->size2);   // se
-    enforce(Y->size1 > 0);
-    gsl_matrix *K = gsl_matrix_safe_alloc(Y->size1, Y->size1);      // Kinship aka as G/GRM in the equation ind x ind
-    gsl_matrix *U = gsl_matrix_safe_alloc(Y->size1, Y->size1);
-    gsl_matrix *UtW = gsl_matrix_calloc(Y->size1, W->size2);        // Transformed ind x covariates
-    gsl_matrix *UtY = gsl_matrix_calloc(Y->size1, Y->size2);        // Transforemd ind x phenotypes
-    gsl_vector *eval = gsl_vector_calloc(Y->size1);                 // eigen values
-    gsl_vector *env = gsl_vector_safe_alloc(Y->size1);              // E environment
-    gsl_vector *weight = gsl_vector_safe_alloc(Y->size1);           // weights
+    enforce(ni_test > 0);
+    gsl_matrix *K = gsl_matrix_safe_alloc(ni_test, ni_test);      // Kinship aka as G/GRM in the equation ind x ind
+    gsl_matrix *U = gsl_matrix_safe_alloc(ni_test, ni_test);
+    gsl_matrix *UtW = gsl_matrix_calloc(ni_test, W->size2);        // Transformed ind x covariates
+    gsl_matrix *UtY = gsl_matrix_calloc(ni_test, Y->size2);        // Transforemd ind x phenotypes
+    gsl_vector *eval = gsl_vector_calloc(ni_test);                 // eigen values
+    gsl_vector *env = gsl_vector_safe_alloc(ni_test);              // E environment
+    gsl_vector *weight = gsl_vector_safe_alloc(ni_test);           // weights
     debug_msg("Started on LMM");
     assert_issue(is_issue(26), UtY->data[0] == 0.0);
 
@@ -2838,28 +2842,32 @@ void GEMMA::BatchRun(PARAM &cPar) {
           gsl_vector_view UtY_col = gsl_matrix_column(UtY, 0); // and its transformed
           write(&Y_col.vector, "Y_col");
 
-          if (!cPar.file_bfile.empty()) {
-            // PLINK analysis
-            if (cPar.file_gxe.empty()) {
-              cLmm.AnalyzePlink(U, eval, UtW, &UtY_col.vector, W,
-                                &Y_col.vector, cPar.setGWASnps);
+          if (cPar.is_mdb)
+            cLmm.mdb_calc_gwa(U, eval, UtW, &UtY_col.vector, W,
+                              &Y_col.vector, cPar.setGWASnps);
+          else
+            if (!cPar.file_bfile.empty()) {
+              // PLINK analysis
+              if (cPar.file_gxe.empty()) {
+                cLmm.AnalyzePlink(U, eval, UtW, &UtY_col.vector, W,
+                                  &Y_col.vector, cPar.setGWASnps);
+              }
+              else {
+                cLmm.AnalyzePlinkGXE(U, eval, UtW, &UtY_col.vector, W,
+                                     &Y_col.vector, env);
+              }
             }
             else {
-              cLmm.AnalyzePlinkGXE(U, eval, UtW, &UtY_col.vector, W,
-                                   &Y_col.vector, env);
+              // BIMBAM analysis
+
+              if (cPar.file_gxe.empty()) {
+                cLmm.AnalyzeBimbam(U, eval, UtW, &UtY_col.vector, W,
+                                   &Y_col.vector, cPar.setGWASnps);
+              } else {
+                cLmm.AnalyzeBimbamGXE(U, eval, UtW, &UtY_col.vector, W,
+                                      &Y_col.vector, env);
+              }
             }
-          }
-          else {
-            // BIMBAM analysis
-
-            if (cPar.file_gxe.empty()) {
-              cLmm.AnalyzeBimbam(U, eval, UtW, &UtY_col.vector, W,
-                                 &Y_col.vector, cPar.setGWASnps);
-            } else {
-              cLmm.AnalyzeBimbamGXE(U, eval, UtW, &UtY_col.vector, W,
-                                    &Y_col.vector, env);
-            }
-          }
           cLmm.WriteFiles(); // write output files
           cLmm.CopyToParam(cPar);
         } else {