diff options
Diffstat (limited to 'src/gemma.cpp')
| -rw-r--r-- | src/gemma.cpp | 34 |
1 files changed, 32 insertions, 2 deletions
diff --git a/src/gemma.cpp b/src/gemma.cpp index c82bf4a..7749bc9 100644 --- a/src/gemma.cpp +++ b/src/gemma.cpp @@ -44,11 +44,20 @@ #include "prdt.h" #include "varcov.h" #include "vc.h" +#include "debug.h" using namespace std; GEMMA::GEMMA(void) : version("0.97.1"), date("08/07/2017"), year("2017") {} +void gemma_gsl_error_handler (const char * reason, + const char * file, + int line, int gsl_errno) { + cerr << "GSL ERROR: " << reason << " in " << file + << " at line " << line << " errno " << gsl_errno <<endl; + exit(22); +} + void GEMMA::PrintHeader(void) { cout << endl; cout << "*********************************************************" << endl; @@ -247,10 +256,10 @@ void GEMMA::PrintHelp(size_t option) { << endl; cout << " to fit a multivariate linear mixed model: " << endl; cout << " ./gemma -bfile [prefix] -k [filename] -lmm [num] -n " - "[num1] [num2] -o [prefix]" + "[pheno cols...] -o [prefix]" << endl; cout << " ./gemma -g [filename] -p [filename] -a [filename] -k " - "[filename] -lmm [num] -n [num1] [num2] -o [prefix]" + "[filename] -lmm [num] -n [pheno cols...] -o [prefix]" << endl; cout << " to fit a Bayesian sparse linear mixed model: " << endl; cout << " ./gemma -bfile [prefix] -bslmm [num] -o [prefix]" << endl; @@ -536,6 +545,7 @@ void GEMMA::PrintHelp(size_t option) { if (option == 9) { cout << " MULTIVARIATE LINEAR MIXED MODEL OPTIONS" << endl; + cout << " -n [pheno cols...] - range of phenotypes" << endl; cout << " -pnr " << " specify the pvalue threshold to use the Newton-Raphson's method " "(default 0.001)" @@ -694,6 +704,7 @@ void GEMMA::PrintHelp(size_t option) { if (option == 14) { cout << " DEBUG OPTIONS" << endl; + cout << " -no-checks disable checks" << endl; cout << " -silence silent terminal display" << endl; cout << " -debug debug output" << endl; cout << " -nind [num] read up to num individuals" << endl; @@ -1578,6 +1589,8 @@ void GEMMA::Assign(int argc, char **argv, PARAM &cPar) { cPar.window_ns = atoi(str.c_str()); } else if (strcmp(argv[i], "-debug") == 0) { cPar.mode_debug = true; + } else if (strcmp(argv[i], "-no-checks") == 0) { + cPar.mode_check = false; } else { cout << "error! unrecognized option: " << argv[i] << endl; cPar.error = true; @@ -1713,6 +1726,7 @@ void GEMMA::BatchRun(PARAM &cPar) { cout << "error! fail to read kinship/relatedness file. " << endl; return; } + // FIXME: this is strange, why read twice? ReadFile_kin(cPar.file_kin, cPar.indicator_cvt, cPar.mapID2num, cPar.k_mode, cPar.error, G_full); if (cPar.error == true) { @@ -1723,6 +1737,7 @@ void GEMMA::BatchRun(PARAM &cPar) { // center matrix G CenterMatrix(G); CenterMatrix(G_full); + validate_K(G,cPar.mode_check); // eigen-decomposition and calculate trace_G cout << "Start Eigen-Decomposition..." << endl; @@ -1869,6 +1884,9 @@ void GEMMA::BatchRun(PARAM &cPar) { return; } + // Now we have the Kinship matrix test it + validate_K(G,cPar.mode_check); + if (cPar.a_mode == 21) { cPar.WriteMatrix(G, "cXX"); } else { @@ -2146,6 +2164,8 @@ void GEMMA::BatchRun(PARAM &cPar) { cPar.se_pve_total, cPar.v_sigma2, cPar.v_se_sigma2, cPar.v_enrich, cPar.v_se_enrich); + assert(!has_nan(cPar.v_se_pve)); + // if LDSC weights, then compute the weights and run the above steps again if (cPar.a_mode == 62) { // compute the weights and normalize the weights for A @@ -2175,8 +2195,10 @@ void GEMMA::BatchRun(PARAM &cPar) { cPar.ni_study, cPar.v_pve, cPar.v_se_pve, cPar.pve_total, cPar.se_pve_total, cPar.v_sigma2, cPar.v_se_sigma2, cPar.v_enrich, cPar.v_se_enrich); + assert(!has_nan(cPar.v_se_pve)); } + gsl_vector_set(s, cPar.n_vc, cPar.ni_test); cPar.WriteMatrix(S, "S"); @@ -2265,6 +2287,7 @@ void GEMMA::BatchRun(PARAM &cPar) { cPar.v_pve, cPar.v_se_pve, cPar.pve_total, cPar.se_pve_total, cPar.v_sigma2, cPar.v_se_sigma2, cPar.v_enrich, cPar.v_se_enrich); + assert(!has_nan(cPar.v_se_pve)); gsl_vector_view s_sub = gsl_vector_subvector(s, 0, cPar.n_vc); gsl_vector_memcpy(&s_sub.vector, s_ref); @@ -2325,6 +2348,7 @@ void GEMMA::BatchRun(PARAM &cPar) { // center matrix G CenterMatrix(G); + validate_K(G,cPar.mode_check); (cPar.v_traceG).clear(); double d = 0; @@ -2636,6 +2660,7 @@ return;} CalcCIss(Xz, XWz, XtXWz, S, Svar, w, z, s_vec, vec_cat, cPar.v_pve, cPar.v_se_pve, cPar.pve_total, cPar.se_pve_total, cPar.v_sigma2, cPar.v_se_sigma2, cPar.v_enrich, cPar.v_se_enrich); + assert(!has_nan(cPar.v_se_pve)); // write files // cPar.WriteMatrix (XWz, "XWz"); @@ -2692,6 +2717,7 @@ return;} // center matrix G CenterMatrix(G); + validate_K(G,cPar.mode_check); // is residual weights are provided, then if (!cPar.file_weight.empty()) { @@ -3007,6 +3033,7 @@ return;} // center matrix G CenterMatrix(G); + validate_K(G,cPar.mode_check); } else { cPar.ReadGenotypes(UtX, G, true); } @@ -3123,6 +3150,8 @@ return;} // center matrix G CenterMatrix(G); + validate_K(G,cPar.mode_check); + } else { cPar.ReadGenotypes(UtX, G, true); } @@ -3329,6 +3358,7 @@ void GEMMA::WriteLog(int argc, char **argv, PARAM &cPar) { outfile << " " << cPar.v_se_pve[i]; } outfile << endl; + assert(!has_nan(cPar.v_se_pve)); if (cPar.n_vc > 1) { outfile << "## total pve = " << cPar.pve_total << endl; |