aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/gemma.cpp24
-rw-r--r--src/gemma.h6
-rw-r--r--src/lmm.cpp21
3 files changed, 30 insertions, 21 deletions
diff --git a/src/gemma.cpp b/src/gemma.cpp
index 1938e56..5a6d4d8 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -2550,8 +2550,8 @@ void GEMMA::BatchRun(PARAM &cPar) {
// LMM or mvLMM or Eigen-Decomposition
if (cPar.a_mode == M_LMM1 || cPar.a_mode == M_LMM2 || cPar.a_mode == M_LMM3 ||
- cPar.a_mode == M_LMM4 || cPar.a_mode == M_LMM5 ||
- cPar.a_mode == M_EIGEN) { // Fit LMM or mvLMM or eigen
+ 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
gsl_matrix *Y = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_ph);
enforce_msg(Y, "allocate Y"); // just to be sure
gsl_matrix *W = gsl_matrix_safe_alloc(Y->size1, cPar.n_cvt);
@@ -2565,6 +2565,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
gsl_vector *eval = gsl_vector_calloc(Y->size1);
gsl_vector *env = gsl_vector_safe_alloc(Y->size1);
gsl_vector *weight = gsl_vector_safe_alloc(Y->size1);
+ debug_msg("Started on LMM");
assert_issue(is_issue(26), UtY->data[0] == 0.0);
// set covariates matrix W and phenotype matrix Y
@@ -2578,6 +2579,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
if (!(cPar.file_kin).empty()) {
ReadFile_kin(cPar.file_kin, cPar.indicator_idv, cPar.mapID2num,
cPar.k_mode, cPar.error, G);
+ debug_msg("Read K/GRM file");
if (cPar.error == true) {
cout << "error! fail to read kinship/relatedness file. " << endl;
return;
@@ -2609,7 +2611,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
}
}
- // eigen-decomposition and calculate trace_G
+ // eigen-decomposition and calculate trace_G - main track
cout << "Start Eigen-Decomposition..." << endl;
time_start = clock();
@@ -2663,7 +2665,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
if (cPar.a_mode == M_EIGEN) {
cPar.WriteMatrix(U, "eigenU");
cPar.WriteVector(eval, "eigenD");
- } else if (!cPar.file_gene.empty()) {
+ } else if (!cPar.file_gene.empty()) { // Run with gene file
// calculate UtW and Uty
CalcUtX(U, W, UtW);
CalcUtX(U, Y, UtY);
@@ -2682,6 +2684,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
cLmm.WriteFiles();
cLmm.CopyToParam(cPar);
} else {
+ debug_msg("Main LMM track");
// calculate UtW and Uty
CalcUtX(U, W, UtW);
CalcUtX(U, Y, UtY);
@@ -2736,6 +2739,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
CalcPve(eval, UtW, &UtY_col.vector, cPar.l_remle_null, cPar.trace_G,
cPar.pve_null, cPar.pve_se_null);
+ debug_msg("main print summary");
cPar.PrintSummary();
// calculate and output residuals
@@ -2771,15 +2775,16 @@ void GEMMA::BatchRun(PARAM &cPar) {
gsl_vector_safe_free(u_hat);
gsl_vector_safe_free(e_hat);
gsl_vector_safe_free(y_hat);
- }
+ } // output residuals
}
// Fit LMM or mvLMM (w. LOCO)
- if (cPar.a_mode == 1 || cPar.a_mode == 2 || cPar.a_mode == 3 ||
- cPar.a_mode == 4) {
+ if (cPar.a_mode == M_LMM1 || cPar.a_mode == M_LMM2 || cPar.a_mode == M_LMM3 ||
+ cPar.a_mode == M_LMM4 || cPar.a_mode == M_LMM9) {
if (cPar.n_ph == 1) {
+ debug_msg("fit LMM (one phenotype)");
LMM cLmm;
- cLmm.CopyFromParam(cPar);
+ cLmm.CopyFromParam(cPar); // set parameters
// if (is_check_mode()) disable_segfpe(); // disable fast NaN checking for now
@@ -2811,8 +2816,9 @@ void GEMMA::BatchRun(PARAM &cPar) {
cLmm.WriteFiles();
cLmm.CopyToParam(cPar);
} else {
+ debug_msg("fit mvLMM (multiple phenotypes)");
MVLMM cMvlmm;
- cMvlmm.CopyFromParam(cPar);
+ cMvlmm.CopyFromParam(cPar); // set parameters
// if (is_check_mode()) disable_segfpe(); // disable fast NaN checking
diff --git a/src/gemma.h b/src/gemma.h
index 67a7f58..b34a958 100644
--- a/src/gemma.h
+++ b/src/gemma.h
@@ -41,8 +41,10 @@ using namespace std;
// gw: 72
enum M_MODE { M_LMM1=1, M_LMM2=2, M_LMM3=3, M_LMM4=4, M_LMM5=5,
- M_BSLMM5=15,
- M_KIN=21, M_KIN2=22, M_EIGEN=31 };
+ M_LMM9=9, // GeneNetwork mode
+ M_BSLMM5=15,
+ M_KIN=21, M_KIN2=22, M_EIGEN=31
+};
class GEMMA {
diff --git a/src/lmm.cpp b/src/lmm.cpp
index e10c0e9..9dd4aa6 100644
--- a/src/lmm.cpp
+++ b/src/lmm.cpp
@@ -42,6 +42,7 @@
#include "gsl/gsl_vector.h"
#include "gzstream.h"
+#include "gemma.h"
#include "gemma_io.h"
#include "fastblas.h"
#include "lapack.h"
@@ -109,27 +110,27 @@ void LMM::WriteFiles() {
}
auto common_header = [&] () {
- if (a_mode != 2) {
+ if (a_mode != M_LMM2) {
outfile << "beta" << "\t";
outfile << "se" << "\t";
}
- if (a_mode != 3)
- outfile << "logl_H1" << "\t"; // we may make this an option
+ if (a_mode != M_LMM3)
+ outfile << "logl_H1" << "\t";
switch(a_mode) {
- case 1:
+ case M_LMM1:
outfile << "l_remle" << "\t"
<< "p_wald" << endl;
break;
- case 2:
+ case M_LMM2:
outfile << "l_mle" << "\t"
<< "p_lrt" << endl;
break;
- case 3:
+ case M_LMM3:
outfile << "p_score" << endl;
break;
- case 4:
+ case M_LL4:
outfile << "l_remle" << "\t"
<< "l_mle" << "\t"
<< "p_wald" << "\t"
@@ -1848,17 +1849,17 @@ void LMM::AnalyzePlink(const gsl_matrix *U, const gsl_vector *eval,
FUNC_PARAM param1 = {false, ni_test, n_cvt, eval, Uab, ab, 0};
// 3 is before 1, for beta.
- if (a_mode == 3 || a_mode == 4) {
+ if (a_mode == M_LMM3 || a_mode == M_LMM4) {
CalcRLScore(l_mle_null, param1, beta, se, p_score);
}
- if (a_mode == 1 || a_mode == 4) {
+ if (a_mode == M_LMM1 || a_mode == M_LMM4) {
CalcLambda('R', param1, l_min, l_max, n_region, lambda_remle,
logl_H1);
CalcRLWald(lambda_remle, param1, beta, se, p_wald);
}
- if (a_mode == 2 || a_mode == 4) {
+ if (a_mode == M_LMM2 || a_mode == M_LMM4) {
CalcLambda('L', param1, l_min, l_max, n_region, lambda_mle, logl_H1);
p_lrt = gsl_cdf_chisq_Q(2.0 * (logl_H1 - logl_mle_H0), 1);
}
generated by cgit v1.2.3 (git 2.44.0) at 2025-01-18 22:50:45 +0000