mdb gwa

2 files changed, 26 insertions, 25 deletions

diff --git a/src/gemma.cpp b/src/gemma.cpp
index fb2227e..0145563 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -2591,21 +2591,20 @@ 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");
-    auto ni_test = cPar.ni_test;
-    auto n_ph = cPar.n_ph;
-    enforce(ni_test > 0);
+    auto n_ph    = cPar.n_ph;                                      // # of trait values
+    auto ni_test = cPar.ni_test;                                   // # of individuals to test
     enforce(n_ph > 0);
-    gsl_matrix *Y = gsl_matrix_safe_alloc(ni_test, n_ph); // Y is ind x phenotypes
+    enforce(ni_test > 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(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
+    gsl_matrix *B = gsl_matrix_safe_alloc(n_ph, W->size2);         // B is a d by c (effect size)
+    gsl_matrix *se_B = gsl_matrix_safe_alloc(n_ph, W->size2);      // se
     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 *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_matrix *UtY = gsl_matrix_calloc(ni_test, n_ph);            // Transformed 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
diff --git a/src/param.cpp b/src/param.cpp
index f799ef7..955da8e 100644
--- a/src/param.cpp
+++ b/src/param.cpp
@@ -149,26 +149,28 @@ void PARAM::ReadFiles(void) {
     }
   }
 
-  // Read SNP set.
-  if (!file_snps.empty()) {
-    if (ReadFile_snps(file_snps, setSnps) == false) {
-      error = true;
-    }
-  } else {
-    setSnps.clear();
+  if (file_geno.find(".mdb") != string::npos) {
+    is_mdb = true;
   }
 
-  // Read KSNP set.
-  if (!file_ksnps.empty()) {
-    if (ReadFile_snps(file_ksnps, setKSnps) == false) {
-      error = true;
+  if (!is_mdb) {
+    // Read SNP set into setSnps (without filtering)
+    if (!file_snps.empty()) {
+      if (ReadFile_snps(file_snps, setSnps) == false) {
+        error = true;
+      }
+    } else {
+      setSnps.clear();
     }
-  } else {
-    setKSnps.clear();
-  }
 
-  if (file_geno.find(".mdb") != string::npos) {
-    is_mdb = true;
+    // Also read KSNP set. Snps used by GRM.
+    if (!file_ksnps.empty()) {
+      if (ReadFile_snps(file_ksnps, setKSnps) == false) {
+        error = true;
+      }
+    } else {
+      setKSnps.clear();
+    }
   }
 
   // For prediction.