aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/io.cpp50
-rw-r--r--src/lm.cpp4
-rw-r--r--src/lmm.cpp5
-rw-r--r--src/mvlmm.cpp5
-rw-r--r--src/prdt.cpp2
-rw-r--r--src/varcov.cpp2
-rw-r--r--src/vc.cpp10
-rwxr-xr-xtest/dev_test_suite.sh7
-rwxr-xr-xtest/test_suite.sh35
9 files changed, 112 insertions, 8 deletions
diff --git a/src/io.cpp b/src/io.cpp
index be40437..6be01fd 100644
--- a/src/io.cpp
+++ b/src/io.cpp
@@ -135,6 +135,7 @@ std::istream &safeGetline(std::istream &is, std::string &t) {
// Read SNP file. A single column of SNP names.
bool ReadFile_snps(const string file_snps, set<string> &setSnps) {
+ debug_msg("entered");
setSnps.clear();
igzstream infile(file_snps.c_str(), igzstream::in);
@@ -162,6 +163,7 @@ bool ReadFile_snps(const string file_snps, set<string> &setSnps) {
// values for each row are parsed. A valid header can be, for example,
// RS POS CHR
bool ReadFile_snps_header(const string &file_snps, set<string> &setSnps) {
+ debug_msg("entered");
setSnps.clear();
igzstream infile(file_snps.c_str(), igzstream::in);
@@ -219,6 +221,7 @@ bool ReadFile_snps_header(const string &file_snps, set<string> &setSnps) {
// Read log file.
bool ReadFile_log(const string &file_log, double &pheno_mean) {
+ debug_msg("ReadFile_log");
ifstream infile(file_log.c_str(), ifstream::in);
if (!infile) {
cout << "error! fail to open log file: " << file_log << endl;
@@ -260,6 +263,7 @@ bool ReadFile_log(const string &file_log, double &pheno_mean) {
bool ReadFile_anno(const string &file_anno, map<string, string> &mapRS2chr,
map<string, long int> &mapRS2bp,
map<string, double> &mapRS2cM) {
+ debug_msg("ReadFile_anno");
mapRS2chr.clear();
mapRS2bp.clear();
@@ -319,6 +323,7 @@ bool ReadFile_anno(const string &file_anno, map<string, string> &mapRS2chr,
// Read 1 column of phenotype.
bool ReadFile_column(const string &file_pheno, vector<int> &indicator_idv,
vector<double> &pheno, const int &p_column) {
+ debug_msg("entered");
indicator_idv.clear();
pheno.clear();
@@ -361,6 +366,7 @@ bool ReadFile_pheno(const string &file_pheno,
vector<vector<int>> &indicator_pheno,
vector<vector<double>> &pheno,
const vector<size_t> &p_column) {
+ debug_msg("entered");
indicator_pheno.clear();
pheno.clear();
@@ -418,6 +424,7 @@ bool ReadFile_pheno(const string &file_pheno,
bool ReadFile_cvt(const string &file_cvt, vector<int> &indicator_cvt,
vector<vector<double>> &cvt, size_t &n_cvt) {
+ debug_msg("entered");
indicator_cvt.clear();
ifstream infile(file_cvt.c_str(), ifstream::in);
@@ -484,6 +491,7 @@ bool ReadFile_cvt(const string &file_cvt, vector<int> &indicator_cvt,
// Read .bim file.
bool ReadFile_bim(const string &file_bim, vector<SNPINFO> &snpInfo) {
+ debug_msg("entered");
snpInfo.clear();
ifstream infile(file_bim.c_str(), ifstream::in);
@@ -529,6 +537,7 @@ bool ReadFile_bim(const string &file_bim, vector<SNPINFO> &snpInfo) {
bool ReadFile_fam(const string &file_fam, vector<vector<int>> &indicator_pheno,
vector<vector<double>> &pheno, map<string, int> &mapID2num,
const vector<size_t> &p_column) {
+ debug_msg("entered");
indicator_pheno.clear();
pheno.clear();
mapID2num.clear();
@@ -612,6 +621,7 @@ bool ReadFile_geno(const string &file_geno, const set<string> &setSnps,
map<string, long int> &mapRS2bp,
map<string, double> &mapRS2cM, vector<SNPINFO> &snpInfo,
size_t &ns_test, bool debug) {
+ debug_msg("entered");
indicator_snp.clear();
snpInfo.clear();
@@ -822,6 +832,7 @@ bool ReadFile_bed(const string &file_bed, const set<string> &setSnps,
const double &maf_level, const double &miss_level,
const double &hwe_level, const double &r2_level,
size_t &ns_test) {
+ debug_msg("entered");
indicator_snp.clear();
size_t ns_total = snpInfo.size();
@@ -1010,6 +1021,7 @@ bool ReadFile_bed(const string &file_bed, const set<string> &setSnps,
// Missing values are replaced by mean.
bool Bimbam_ReadOneSNP(const size_t inc, const vector<int> &indicator_idv,
igzstream &infile, gsl_vector *geno, double &geno_mean) {
+ debug_msg("entered");
size_t ni_total = indicator_idv.size();
string line;
@@ -1060,6 +1072,7 @@ bool Bimbam_ReadOneSNP(const size_t inc, const vector<int> &indicator_idv,
// For PLINK, store SNPs as double too.
void Plink_ReadOneSNP(const int pos, const vector<int> &indicator_idv,
ifstream &infile, gsl_vector *geno, double &geno_mean) {
+ debug_msg("entered");
size_t ni_total = indicator_idv.size(), n_bit;
if (ni_total % 4 == 0) {
n_bit = ni_total / 4;
@@ -1126,6 +1139,7 @@ void Plink_ReadOneSNP(const int pos, const vector<int> &indicator_idv,
void ReadFile_kin(const string &file_kin, vector<int> &indicator_idv,
map<string, int> &mapID2num, const size_t k_mode, bool &error,
gsl_matrix *G) {
+ debug_msg("entered");
igzstream infile(file_kin.c_str(), igzstream::in);
if (!infile) {
cout << "error! fail to open kinship file: " << file_kin << endl;
@@ -1242,6 +1256,7 @@ void ReadFile_kin(const string &file_kin, vector<int> &indicator_idv,
void ReadFile_mk(const string &file_mk, vector<int> &indicator_idv,
map<string, int> &mapID2num, const size_t k_mode, bool &error,
gsl_matrix *G) {
+ debug_msg("entered");
igzstream infile(file_mk.c_str(), igzstream::in);
if (!infile) {
cout << "error! fail to open file: " << file_mk << endl;
@@ -1267,6 +1282,7 @@ void ReadFile_mk(const string &file_mk, vector<int> &indicator_idv,
}
void ReadFile_eigenU(const string &file_ku, bool &error, gsl_matrix *U) {
+ debug_msg("entered");
igzstream infile(file_ku.c_str(), igzstream::in);
if (!infile) {
cout << "error! fail to open the U file: " << file_ku << endl;
@@ -1315,6 +1331,7 @@ void ReadFile_eigenU(const string &file_ku, bool &error, gsl_matrix *U) {
}
void ReadFile_eigenD(const string &file_kd, bool &error, gsl_vector *eval) {
+ debug_msg("entered");
igzstream infile(file_kd.c_str(), igzstream::in);
if (!infile) {
cout << "error! fail to open the D file: " << file_kd << endl;
@@ -1363,6 +1380,7 @@ bool BimbamKin(const string file_geno, const set<string> ksnps,
vector<int> &indicator_snp, const int k_mode,
const int display_pace, gsl_matrix *matrix_kin,
const bool test_nind) {
+ debug_msg("entered");
igzstream infile(file_geno.c_str(), igzstream::in);
enforce_msg(infile, "error reading genotype file");
@@ -1499,6 +1517,7 @@ bool BimbamKin(const string file_geno, const set<string> ksnps,
bool PlinkKin(const string &file_bed, vector<int> &indicator_snp,
const int k_mode, const int display_pace,
gsl_matrix *matrix_kin) {
+ debug_msg("entered");
ifstream infile(file_bed.c_str(), ios::binary);
if (!infile) {
cout << "error reading bed file:" << file_bed << endl;
@@ -1641,6 +1660,7 @@ bool PlinkKin(const string &file_bed, vector<int> &indicator_snp,
bool ReadFile_geno(const string file_geno, vector<int> &indicator_idv,
vector<int> &indicator_snp, gsl_matrix *UtX, gsl_matrix *K,
const bool calc_K, bool debug) {
+ debug_msg("entered");
igzstream infile(file_geno.c_str(), igzstream::in);
if (!infile) {
cout << "error reading genotype file:" << file_geno << endl;
@@ -1745,6 +1765,7 @@ bool ReadFile_geno(const string &file_geno, vector<int> &indicator_idv,
vector<vector<unsigned char>> &Xt, gsl_matrix *K,
const bool calc_K, const size_t ni_test,
const size_t ns_test, bool debug) {
+ debug_msg("entered");
igzstream infile(file_geno.c_str(), igzstream::in);
if (!infile) {
cout << "error reading genotype file:" << file_geno << endl;
@@ -1851,6 +1872,7 @@ bool ReadFile_geno(const string &file_geno, vector<int> &indicator_idv,
bool ReadFile_bed(const string &file_bed, vector<int> &indicator_idv,
vector<int> &indicator_snp, gsl_matrix *UtX, gsl_matrix *K,
const bool calc_K) {
+ debug_msg("entered");
ifstream infile(file_bed.c_str(), ios::binary);
if (!infile) {
cout << "error reading bed file:" << file_bed << endl;
@@ -1982,6 +2004,7 @@ bool ReadFile_bed(const string &file_bed, vector<int> &indicator_idv,
vector<int> &indicator_snp, vector<vector<unsigned char>> &Xt,
gsl_matrix *K, const bool calc_K, const size_t ni_test,
const size_t ns_test) {
+ debug_msg("entered");
ifstream infile(file_bed.c_str(), ios::binary);
if (!infile) {
cout << "error reading bed file:" << file_bed << endl;
@@ -2116,6 +2139,7 @@ bool ReadFile_bed(const string &file_bed, vector<int> &indicator_idv,
bool ReadFile_est(const string &file_est, const vector<size_t> &est_column,
map<string, double> &mapRS2est) {
+ debug_msg("entered");
mapRS2est.clear();
ifstream infile(file_est.c_str(), ifstream::in);
@@ -2176,6 +2200,7 @@ bool ReadFile_est(const string &file_est, const vector<size_t> &est_column,
}
bool CountFileLines(const string &file_input, size_t &n_lines) {
+ debug_msg("entered");
igzstream infile(file_input.c_str(), igzstream::in);
if (!infile) {
cout << "error! fail to open file: " << file_input << endl;
@@ -2192,6 +2217,7 @@ bool CountFileLines(const string &file_input, size_t &n_lines) {
// Read gene expression file.
bool ReadFile_gene(const string &file_gene, vector<double> &vec_read,
vector<SNPINFO> &snpInfo, size_t &ng_total) {
+ debug_msg("entered");
vec_read.clear();
ng_total = 0;
@@ -2255,6 +2281,7 @@ bool ReadFile_sample(const string &file_sample,
vector<vector<double>> &pheno,
const vector<size_t> &p_column, vector<int> &indicator_cvt,
vector<vector<double>> &cvt, size_t &n_cvt) {
+ debug_msg("entered");
indicator_pheno.clear();
pheno.clear();
indicator_cvt.clear();
@@ -2509,6 +2536,7 @@ bool ReadFile_bgen(const string &file_bgen, const set<string> &setSnps,
const double &hwe_level, const double &r2_level,
size_t &ns_test) {
+ debug_msg("entered");
indicator_snp.clear();
ifstream infile(file_bgen.c_str(), ios::binary);
@@ -2796,6 +2824,7 @@ bool ReadFile_bgen(const string &file_bgen, const set<string> &setSnps,
// Read oxford genotype file and calculate kinship matrix.
bool bgenKin(const string &file_oxford, vector<int> &indicator_snp,
const int k_mode, const int display_pace, gsl_matrix *matrix_kin) {
+ debug_msg("entered");
string file_bgen = file_oxford;
ifstream infile(file_bgen.c_str(), ios::binary);
if (!infile) {
@@ -2854,7 +2883,7 @@ bool bgenKin(const string &file_oxford, vector<int> &indicator_snp,
for (size_t t = 0; t < indicator_snp.size(); ++t) {
if (t % display_pace == 0 || t == (indicator_snp.size() - 1)) {
- ProgressBar("Reading SNPs ", t, indicator_snp.size() - 1);
+ ProgressBar("Reading bgen SNPs ", t, indicator_snp.size() - 1);
}
id.clear();
@@ -3000,6 +3029,7 @@ bool bgenKin(const string &file_oxford, vector<int> &indicator_snp,
// Read header to determine which column contains which item.
bool ReadHeader_io(const string &line, HEADER &header) {
+ debug_msg("entered");
string rs_ptr[] = {"rs", "RS", "snp", "SNP", "snps", "SNPS",
"snpid", "SNPID", "rsid", "RSID", "MarkerName"};
set<string> rs_set(rs_ptr, rs_ptr + 11); // create a set of 11 items
@@ -3267,6 +3297,7 @@ bool ReadHeader_io(const string &line, HEADER &header) {
// it is not included in the analysis.
bool ReadFile_cat(const string &file_cat, map<string, size_t> &mapRS2cat,
size_t &n_vc) {
+ debug_msg("entered");
mapRS2cat.clear();
igzstream infile(file_cat.c_str(), igzstream::in);
@@ -3351,6 +3382,7 @@ bool ReadFile_cat(const string &file_cat, map<string, size_t> &mapRS2cat,
bool ReadFile_mcat(const string &file_mcat, map<string, size_t> &mapRS2cat,
size_t &n_vc) {
+ debug_msg("entered");
mapRS2cat.clear();
igzstream infile(file_mcat.c_str(), igzstream::in);
@@ -3389,6 +3421,8 @@ bool BimbamKinUncentered(const string &file_geno, const set<string> ksnps,
const map<string, size_t> &mapRS2cat,
const vector<SNPINFO> &snpInfo, const gsl_matrix *W,
gsl_matrix *matrix_kin, gsl_vector *vector_ns) {
+ debug_msg("entered");
+ debug_msg("BimbamKinUncentered");
igzstream infile(file_geno.c_str(), igzstream::in);
if (!infile) {
cout << "error reading genotype file:" << file_geno << endl;
@@ -3579,6 +3613,7 @@ bool PlinkKin(const string &file_bed, const int display_pace,
const map<string, size_t> &mapRS2cat,
const vector<SNPINFO> &snpInfo, const gsl_matrix *W,
gsl_matrix *matrix_kin, gsl_vector *vector_ns) {
+ debug_msg("entered");
ifstream infile(file_bed.c_str(), ios::binary);
if (!infile) {
cout << "error reading bed file:" << file_bed << endl;
@@ -3804,6 +3839,7 @@ bool MFILEKin(const size_t mfile_mode, const string &file_mfile,
const map<string, size_t> &mapRS2cat,
const vector<vector<SNPINFO>> &msnpInfo, const gsl_matrix *W,
gsl_matrix *matrix_kin, gsl_vector *vector_ns) {
+ debug_msg("entered");
size_t n_vc = vector_ns->size, ni_test = matrix_kin->size1;
gsl_matrix_set_zero(matrix_kin);
gsl_vector_set_zero(vector_ns);
@@ -3879,6 +3915,7 @@ bool MFILEKin(const size_t mfile_mode, const string &file_mfile,
// Read var file, store mapRS2wsnp.
bool ReadFile_wsnp(const string &file_wsnp, map<string, double> &mapRS2weight) {
+ debug_msg("entered");
mapRS2weight.clear();
igzstream infile(file_wsnp.c_str(), igzstream::in);
@@ -3904,6 +3941,7 @@ bool ReadFile_wsnp(const string &file_wsnp, map<string, double> &mapRS2weight) {
bool ReadFile_wsnp(const string &file_wcat, const size_t n_vc,
map<string, vector<double>> &mapRS2wvector) {
+ debug_msg("entered");
mapRS2wvector.clear();
igzstream infile(file_wcat.c_str(), igzstream::in);
@@ -3988,6 +4026,7 @@ void ReadFile_beta(const string &file_beta,
vector<size_t> &vec_ni, vector<double> &vec_weight,
vector<double> &vec_z2, size_t &ni_total, size_t &ns_total,
size_t &ns_test) {
+ debug_msg("entered");
vec_cat.clear();
vec_ni.clear();
vec_weight.clear();
@@ -4174,6 +4213,7 @@ void ReadFile_beta(const string &file_beta,
void ReadFile_beta(const string &file_beta, const map<string, double> &mapRS2wA,
map<string, string> &mapRS2A1,
map<string, double> &mapRS2z) {
+ debug_msg("entered");
mapRS2A1.clear();
mapRS2z.clear();
@@ -4335,6 +4375,7 @@ void Calcq(const size_t n_block, const vector<size_t> &vec_cat,
const vector<size_t> &vec_ni, const vector<double> &vec_weight,
const vector<double> &vec_z2, gsl_matrix *Vq, gsl_vector *q,
gsl_vector *s) {
+ debug_msg("entered");
gsl_matrix_set_zero(Vq);
gsl_vector_set_zero(q);
gsl_vector_set_zero(s);
@@ -4488,6 +4529,7 @@ void Calcq(const size_t n_block, const vector<size_t> &vec_cat,
// Read vector file.
void ReadFile_vector(const string &file_vec, gsl_vector *vec) {
+ debug_msg("entered");
igzstream infile(file_vec.c_str(), igzstream::in);
if (!infile) {
cout << "error! fail to open vector file: " << file_vec << endl;
@@ -4510,6 +4552,7 @@ void ReadFile_vector(const string &file_vec, gsl_vector *vec) {
}
void ReadFile_matrix(const string &file_mat, gsl_matrix *mat) {
+ debug_msg("entered");
igzstream infile(file_mat.c_str(), igzstream::in);
if (!infile) {
cout << "error! fail to open matrix file: " << file_mat << endl;
@@ -4536,6 +4579,7 @@ void ReadFile_matrix(const string &file_mat, gsl_matrix *mat) {
void ReadFile_matrix(const string &file_mat, gsl_matrix *mat1,
gsl_matrix *mat2) {
+ debug_msg("entered");
igzstream infile(file_mat.c_str(), igzstream::in);
if (!infile) {
cout << "error! fail to open matrix file: " << file_mat << endl;
@@ -4572,6 +4616,7 @@ void ReadFile_matrix(const string &file_mat, gsl_matrix *mat1,
// Read study file.
void ReadFile_study(const string &file_study, gsl_matrix *Vq_mat,
gsl_vector *q_vec, gsl_vector *s_vec, size_t &ni) {
+ debug_msg("entered");
string Vqfile = file_study + ".Vq.txt";
string sfile = file_study + ".size.txt";
string qfile = file_study + ".q.txt";
@@ -4597,6 +4642,7 @@ void ReadFile_study(const string &file_study, gsl_matrix *Vq_mat,
// Read reference file.
void ReadFile_ref(const string &file_ref, gsl_matrix *S_mat,
gsl_matrix *Svar_mat, gsl_vector *s_vec, size_t &ni) {
+ debug_msg("entered");
string sfile = file_ref + ".size.txt";
string Sfile = file_ref + ".S.txt";
@@ -4620,6 +4666,7 @@ void ReadFile_ref(const string &file_ref, gsl_matrix *S_mat,
// Read mstudy file.
void ReadFile_mstudy(const string &file_mstudy, gsl_matrix *Vq_mat,
gsl_vector *q_vec, gsl_vector *s_vec, size_t &ni) {
+ debug_msg("entered");
gsl_matrix_set_zero(Vq_mat);
gsl_vector_set_zero(q_vec);
gsl_vector_set_zero(s_vec);
@@ -4710,6 +4757,7 @@ void ReadFile_mstudy(const string &file_mstudy, gsl_matrix *Vq_mat,
// Read reference file.
void ReadFile_mref(const string &file_mref, gsl_matrix *S_mat,
gsl_matrix *Svar_mat, gsl_vector *s_vec, size_t &ni) {
+ debug_msg("entered");
gsl_matrix_set_zero(S_mat);
gsl_matrix_set_zero(Svar_mat);
gsl_vector_set_zero(s_vec);
diff --git a/src/lm.cpp b/src/lm.cpp
index f8fc43d..83558e4 100644
--- a/src/lm.cpp
+++ b/src/lm.cpp
@@ -290,6 +290,7 @@ void LmCalcP(const size_t test_mode, const double yPwy, const double xPwy,
}
void LM::AnalyzeGene(const gsl_matrix *W, const gsl_vector *x) {
+ debug_msg("entering");
ifstream infile(file_gene.c_str(), ifstream::in);
if (!infile) {
cout << "error reading gene expression file:" << file_gene << endl;
@@ -382,6 +383,7 @@ void LM::AnalyzeGene(const gsl_matrix *W, const gsl_vector *x) {
// WJA added
void LM::Analyzebgen(const gsl_matrix *W, const gsl_vector *y) {
+ debug_msg("entering");
string file_bgen = file_oxford + ".bgen";
ifstream infile(file_bgen.c_str(), ios::binary);
if (!infile) {
@@ -606,6 +608,7 @@ void LM::Analyzebgen(const gsl_matrix *W, const gsl_vector *y) {
}
void LM::AnalyzeBimbam(const gsl_matrix *W, const gsl_vector *y) {
+ debug_msg("entering");
igzstream infile(file_geno.c_str(), igzstream::in);
if (!infile) {
cout << "error reading genotype file:" << file_geno << endl;
@@ -720,6 +723,7 @@ void LM::AnalyzeBimbam(const gsl_matrix *W, const gsl_vector *y) {
}
void LM::AnalyzePlink(const gsl_matrix *W, const gsl_vector *y) {
+ debug_msg("entering");
string file_bed = file_bfile + ".bed";
ifstream infile(file_bed.c_str(), ios::binary);
if (!infile) {
diff --git a/src/lmm.cpp b/src/lmm.cpp
index eb76265..37f2f5b 100644
--- a/src/lmm.cpp
+++ b/src/lmm.cpp
@@ -1319,6 +1319,7 @@ void LMM::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,
const set<string> gwasnps) {
+ debug_msg("entering");
clock_t time_start = clock();
// LOCO support
@@ -1482,6 +1483,7 @@ void LMM::AnalyzeBimbam(const gsl_matrix *U, const gsl_vector *eval,
void LMM::AnalyzePlink(const gsl_matrix *U, const gsl_vector *eval,
const gsl_matrix *UtW, const gsl_vector *Uty,
const gsl_matrix *W, const gsl_vector *y) {
+ debug_msg("entering");
string file_bed = file_bfile + ".bed";
ifstream infile(file_bed.c_str(), ios::binary);
if (!infile) {
@@ -1676,6 +1678,7 @@ void LMM::AnalyzePlink(const gsl_matrix *U, const gsl_vector *eval,
void LMM::Analyzebgen(const gsl_matrix *U, const gsl_vector *eval,
const gsl_matrix *UtW, const gsl_vector *Uty,
const gsl_matrix *W, const gsl_vector *y) {
+ debug_msg("entering");
string file_bgen = file_oxford + ".bgen";
ifstream infile(file_bgen.c_str(), ios::binary);
if (!infile) {
@@ -2282,6 +2285,7 @@ void LMM::AnalyzeBimbamGXE(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_vector *env) {
+ debug_msg("entering");
igzstream infile(file_geno.c_str(), igzstream::in);
if (!infile) {
cout << "error reading genotype file:" << file_geno << endl;
@@ -2430,6 +2434,7 @@ void LMM::AnalyzePlinkGXE(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_vector *env) {
+ debug_msg("entering");
string file_bed = file_bfile + ".bed";
ifstream infile(file_bed.c_str(), ios::binary);
if (!infile) {
diff --git a/src/mvlmm.cpp b/src/mvlmm.cpp
index be9fd78..c5efb6e 100644
--- a/src/mvlmm.cpp
+++ b/src/mvlmm.cpp
@@ -2953,6 +2953,7 @@ double PCRT(const size_t mode, const size_t d_size, const double p_value,
// WJA added.
void MVLMM::Analyzebgen(const gsl_matrix *U, const gsl_vector *eval,
const gsl_matrix *UtW, const gsl_matrix *UtY) {
+ debug_msg("entering");
string file_bgen = file_oxford + ".bgen";
ifstream infile(file_bgen.c_str(), ios::binary);
if (!infile) {
@@ -3501,6 +3502,7 @@ void MVLMM::Analyzebgen(const gsl_matrix *U, const gsl_vector *eval,
void MVLMM::AnalyzeBimbam(const gsl_matrix *U, const gsl_vector *eval,
const gsl_matrix *UtW, const gsl_matrix *UtY) {
+ debug_msg("entering");
igzstream infile(file_geno.c_str(), igzstream::in);
if (!infile) {
cout << "error reading genotype file:" << file_geno << endl;
@@ -3938,6 +3940,7 @@ void MVLMM::AnalyzeBimbam(const gsl_matrix *U, const gsl_vector *eval,
void MVLMM::AnalyzePlink(const gsl_matrix *U, const gsl_vector *eval,
const gsl_matrix *UtW, const gsl_matrix *UtY) {
+ debug_msg("entering");
string file_bed = file_bfile + ".bed";
ifstream infile(file_bed.c_str(), ios::binary);
if (!infile) {
@@ -4475,6 +4478,7 @@ void CalcMvLmmVgVeBeta(const gsl_vector *eval, const gsl_matrix *UtW,
void MVLMM::AnalyzeBimbamGXE(const gsl_matrix *U, const gsl_vector *eval,
const gsl_matrix *UtW, const gsl_matrix *UtY,
const gsl_vector *env) {
+ debug_msg("entering");
igzstream infile(file_geno.c_str(), igzstream::in);
if (!infile) {
cout << "error reading genotype file:" << file_geno << endl;
@@ -4920,6 +4924,7 @@ void MVLMM::AnalyzeBimbamGXE(const gsl_matrix *U, const gsl_vector *eval,
void MVLMM::AnalyzePlinkGXE(const gsl_matrix *U, const gsl_vector *eval,
const gsl_matrix *UtW, const gsl_matrix *UtY,
const gsl_vector *env) {
+ debug_msg("entering");
string file_bed = file_bfile + ".bed";
ifstream infile(file_bed.c_str(), ios::binary);
if (!infile) {
diff --git a/src/prdt.cpp b/src/prdt.cpp
index 3e7c004..9dc84bc 100644
--- a/src/prdt.cpp
+++ b/src/prdt.cpp
@@ -206,6 +206,7 @@ void PRDT::AddBV(gsl_matrix *G, const gsl_vector *u_hat, gsl_vector *y_prdt) {
}
void PRDT::AnalyzeBimbam(gsl_vector *y_prdt) {
+ debug_msg("entering");
igzstream infile(file_geno.c_str(), igzstream::in);
if (!infile) {
cout << "error reading genotype file:" << file_geno << endl;
@@ -308,6 +309,7 @@ void PRDT::AnalyzeBimbam(gsl_vector *y_prdt) {
}
void PRDT::AnalyzePlink(gsl_vector *y_prdt) {
+ debug_msg("entering");
string file_bed = file_bfile + ".bed";
ifstream infile(file_bed.c_str(), ios::binary);
if (!infile) {
diff --git a/src/varcov.cpp b/src/varcov.cpp
index 0f87ba8..39c3523 100644
--- a/src/varcov.cpp
+++ b/src/varcov.cpp
@@ -247,6 +247,7 @@ void Calc_Cor(vector<vector<double>> &X_mat, vector<double> &cov_vec) {
// genotype file and calculate the covariance matrix for neighboring
// SNPs output values at 10000-SNP-interval.
void VARCOV::AnalyzeBimbam() {
+ debug_msg("entering");
igzstream infile(file_geno.c_str(), igzstream::in);
if (!infile) {
cout << "error reading genotype file:" << file_geno << endl;
@@ -345,6 +346,7 @@ void VARCOV::AnalyzeBimbam() {
}
void VARCOV::AnalyzePlink() {
+ debug_msg("entering");
string file_bed = file_bfile + ".bed";
ifstream infile(file_bed.c_str(), ios::binary);
if (!infile) {
diff --git a/src/vc.cpp b/src/vc.cpp
index b5f36c0..1465f16 100644
--- a/src/vc.cpp
+++ b/src/vc.cpp
@@ -407,6 +407,7 @@ int LogRL_dev12(const gsl_vector *log_sigma2, void *params, gsl_vector *dev1,
// Read header to determine which column contains which item.
bool ReadHeader_vc(const string &line, HEADER &header) {
+ debug_msg("entering");
string rs_ptr[] = {"rs", "RS", "snp", "SNP", "snps",
"SNPS", "snpid", "SNPID", "rsid", "RSID"};
set<string> rs_set(rs_ptr, rs_ptr + 10);
@@ -639,6 +640,7 @@ void ReadFile_cor(const string &file_cor, const set<string> &setSnps,
vector<double> &vec_cm, vector<double> &vec_bp,
map<string, size_t> &mapRS2in,
map<string, double> &mapRS2var) {
+ debug_msg("entering");
vec_rs.clear();
vec_n.clear();
mapRS2in.clear();
@@ -790,6 +792,7 @@ void ReadFile_beta(const bool flag_priorscale, const string &file_beta,
map<string, size_t> &mapRS2nsamp, gsl_vector *q_vec,
gsl_vector *qvar_vec, gsl_vector *s_vec, size_t &ni_total,
size_t &ns_total) {
+ debug_msg("entering");
mapRS2nsamp.clear();
igzstream infile(file_beta.c_str(), igzstream::in);
@@ -1004,6 +1007,7 @@ void ReadFile_cor(const string &file_cor, const vector<string> &vec_rs,
const double &window_ns, gsl_matrix *S_mat,
gsl_matrix *Svar_mat, gsl_vector *qvar_vec, size_t &ni_total,
size_t &ns_total, size_t &ns_test, size_t &ns_pair) {
+ debug_msg("entering");
igzstream infile(file_cor.c_str(), igzstream::in);
if (!infile) {
cout << "error! fail to open cov file: " << file_cor << endl;
@@ -2213,6 +2217,7 @@ bool BimbamXwz(const string &file_geno, const int display_pace,
vector<int> &indicator_idv, vector<int> &indicator_snp,
const vector<size_t> &vec_cat, const gsl_vector *w,
const gsl_vector *z, size_t ns_test, gsl_matrix *XWz) {
+ debug_msg("entering");
igzstream infile(file_geno.c_str(), igzstream::in);
if (!infile) {
cout << "error reading genotype file:" << file_geno << endl;
@@ -2306,6 +2311,7 @@ bool PlinkXwz(const string &file_bed, const int display_pace,
vector<int> &indicator_idv, vector<int> &indicator_snp,
const vector<size_t> &vec_cat, const gsl_vector *w,
const gsl_vector *z, size_t ns_test, gsl_matrix *XWz) {
+ debug_msg("entering");
ifstream infile(file_bed.c_str(), ios::binary);
if (!infile) {
cout << "error reading bed file:" << file_bed << endl;
@@ -2432,6 +2438,7 @@ bool MFILEXwz(const size_t mfile_mode, const string &file_mfile,
vector<vector<int>> &mindicator_snp,
const vector<size_t> &vec_cat, const gsl_vector *w,
const gsl_vector *z, gsl_matrix *XWz) {
+ debug_msg("entering");
gsl_matrix_set_zero(XWz);
igzstream infile(file_mfile.c_str(), igzstream::in);
@@ -2466,6 +2473,7 @@ bool MFILEXwz(const size_t mfile_mode, const string &file_mfile,
bool BimbamXtXwz(const string &file_geno, const int display_pace,
vector<int> &indicator_idv, vector<int> &indicator_snp,
const gsl_matrix *XWz, size_t ns_test, gsl_matrix *XtXWz) {
+ debug_msg("entering");
igzstream infile(file_geno.c_str(), igzstream::in);
if (!infile) {
cout << "error reading genotype file:" << file_geno << endl;
@@ -2556,6 +2564,7 @@ bool BimbamXtXwz(const string &file_geno, const int display_pace,
bool PlinkXtXwz(const string &file_bed, const int display_pace,
vector<int> &indicator_idv, vector<int> &indicator_snp,
const gsl_matrix *XWz, size_t ns_test, gsl_matrix *XtXWz) {
+ debug_msg("entering");
ifstream infile(file_bed.c_str(), ios::binary);
if (!infile) {
cout << "error reading bed file:" << file_bed << endl;
@@ -2679,6 +2688,7 @@ bool MFILEXtXwz(const size_t mfile_mode, const string &file_mfile,
const int display_pace, vector<int> &indicator_idv,
vector<vector<int>> &mindicator_snp, const gsl_matrix *XWz,
gsl_matrix *XtXWz) {
+ debug_msg("entering");
gsl_matrix_set_zero(XtXWz);
igzstream infile(file_mfile.c_str(), igzstream::in);
diff --git a/test/dev_test_suite.sh b/test/dev_test_suite.sh
index 9e49251..2bd432e 100755
--- a/test/dev_test_suite.sh
+++ b/test/dev_test_suite.sh
@@ -2,7 +2,6 @@
gemma=../bin/gemma
-
# Related to https://github.com/genetics-statistics/GEMMA/issues/78
testBXDStandardRelatednessMatrixKSingularError() {
outn=BXDerr
@@ -65,10 +64,12 @@ testCenteredRelatednessMatrixKLOCO1() {
testUnivariateLinearMixedModelLOCO1() {
outn=mouse_hs1940_CD8_LOCO1_lmm
rm -f output/$outn.*
- $gemma -g ../example/mouse_hs1940.geno.txt.gz -p ../example/mouse_hs1940.pheno.txt \
+ $gemma -g ../example/mouse_hs1940.geno.txt.gz \
+ -p ../example/mouse_hs1940.pheno.txt \
-n 1 \
-loco 1 \
- -a ../example/mouse_hs1940.anno.txt -k ./output/mouse_hs1940_LOCO1.cXX.txt \
+ -a ../example/mouse_hs1940.anno.txt \
+ -k ./output/mouse_hs1940_LOCO1.cXX.txt \
-snps ../example/mouse_hs1940_snps.txt -lmm \
-nind 400 \
-debug \
diff --git a/test/test_suite.sh b/test/test_suite.sh
index 44eb14c..fa66b7a 100755
--- a/test/test_suite.sh
+++ b/test/test_suite.sh
@@ -36,7 +36,7 @@ testUnivariateLinearMixedModelFullLOCO1() {
testCenteredRelatednessMatrixK() {
$gemma -g ../example/mouse_hs1940.geno.txt.gz \
-p ../example/mouse_hs1940.pheno.txt \
- -gk -o mouse_hs1940
+ -gk -o mouse_hs1940 -debug
assertEquals 0 $?
outfn=output/mouse_hs1940.cXX.txt
assertEquals "1940" `wc -l < $outfn`
@@ -52,7 +52,8 @@ testUnivariateLinearMixedModel() {
-a ../example/mouse_hs1940.anno.txt \
-k ./output/mouse_hs1940.cXX.txt \
-lmm \
- -o mouse_hs1940_CD8_lmm
+ -o mouse_hs1940_CD8_lmm \
+ -debug
assertEquals 0 $?
grep "total computation time" < output/mouse_hs1940_CD8_lmm.log.txt
assertEquals 0 $?
@@ -67,7 +68,8 @@ testMultivariateLinearMixedModel() {
-n 1 6 \
-a ../example/mouse_hs1940.anno.txt \
-k ./output/mouse_hs1940.cXX.txt \
- -lmm -o mouse_hs1940_CD8MCH_lmm
+ -lmm -o mouse_hs1940_CD8MCH_lmm \
+ -debug
assertEquals 0 $?
outfn=output/mouse_hs1940_CD8MCH_lmm.assoc.txt
@@ -81,7 +83,8 @@ testPlinkStandardRelatednessMatrixK() {
outfn=output/$testname.sXX.txt
rm -f $outfn
$gemma -bfile $datadir/HLC \
- -gk 2 -o $testname
+ -gk 2 -o $testname \
+ -debug
assertEquals 0 $?
assertEquals "427" `wc -l < $outfn`
assertEquals "-358.07" `perl -nle 'foreach $x (split(/\s+/,$_)) { $sum += sprintf("%.2f",(substr($x,,0,6))) } END { printf "%.2f",$sum }' $outfn`
@@ -97,6 +100,30 @@ testPlinkMultivariateLinearMixedModel() {
-lmm 1 \
-maf 0.1 \
-c $datadir/HLC_covariates.txt \
+ -debug \
+ -o $testname
+ assertEquals 0 $?
+ outfn=output/$testname.assoc.txt
+ assertEquals "223243" `wc -l < $outfn`
+ assertEquals "89756559859.06" `perl -nle 'foreach $x (split(/\s+/,$_)) { $sum += sprintf("%.2f",(substr($x,,0,6))) } END { printf "%.2f",$sum }' $outfn`
+}
+
+testPlinkMultivariateLinearMixedModelMultiplePhenotypes_Issue58() {
+ # n=2 is original pheno in fam file
+ # n=1 is causal1
+ # n=3..12 is causal2
+ # n=13..22 is causal3
+ # -n 1 2 3 15 is independent
+ testname=testPlinkMultivariateLinearMixedModelMultiplePhenotypes
+ datadir=../example
+ $gemma -bfile $datadir/HLC \
+ -p $datadir/HLC.simu.pheno.txt \
+ -k output/testPlinkStandardRelatednessMatrixK.sXX.txt \
+ -lmm 1 \
+ -maf 0.1 \
+ -n 1 2 3 15 \
+ -c $datadir/HLC_covariates.txt \
+ -debug \
-o $testname
assertEquals 0 $?
outfn=output/$testname.assoc.txt
generated by cgit v1.2.3 (git 2.44.0) at 2025-01-19 03:12:27 +0000