diff options
| -rw-r--r-- | src/gemma.cpp | 12 | ||||
| -rw-r--r-- | src/gemma_io.cpp | 16 | ||||
| -rw-r--r-- | src/gemma_io.h | 2 | ||||
| -rw-r--r-- | src/param.cpp | 15 | ||||
| -rw-r--r-- | src/param.h | 4 |
5 files changed, 26 insertions, 23 deletions
diff --git a/src/gemma.cpp b/src/gemma.cpp index 01e80f7..38b4913 100644 --- a/src/gemma.cpp +++ b/src/gemma.cpp @@ -2901,7 +2901,7 @@ void GEMMA::BatchRun(PARAM &cPar) { // run bvsr if rho==1 if (cPar.rho_min == 1 && cPar.rho_max == 1) { // read genotypes X (not UtX) - cPar.ReadGenotypes(UtX, G, false); + cPar.ReadBIMBAMGenotypes(UtX, G, false); // perform BSLMM analysis BSLMM cBslmm; @@ -2919,7 +2919,7 @@ void GEMMA::BatchRun(PARAM &cPar) { // read relatedness matrix G if (!(cPar.file_kin).empty()) { - cPar.ReadGenotypes(UtX, G, false); + cPar.ReadBIMBAMGenotypes(UtX, G, false); // read relatedness matrix G ReadFile_kin(cPar.file_kin, cPar.indicator_idv, cPar.mapID2num, @@ -2933,7 +2933,7 @@ void GEMMA::BatchRun(PARAM &cPar) { CenterMatrix(G); validate_K(G); } else { - cPar.ReadGenotypes(UtX, G, true); + cPar.ReadBIMBAMGenotypes(UtX, G, true); } // eigen-decomposition and calculate trace_G @@ -3011,7 +3011,7 @@ void GEMMA::BatchRun(PARAM &cPar) { // run bvsr if rho==1 if (cPar.rho_min == 1 && cPar.rho_max == 1) { // read genotypes X (not UtX) - cPar.ReadGenotypes(UtX, G, false); + cPar.ReadBIMBAMGenotypes(UtX, G, false); // perform BSLMM analysis BSLMM cBslmm; @@ -3030,7 +3030,7 @@ void GEMMA::BatchRun(PARAM &cPar) { // read relatedness matrix G if (!(cPar.file_kin).empty()) { - cPar.ReadGenotypes(UtX, G, false); + cPar.ReadBIMBAMGenotypes(UtX, G, false); // read relatedness matrix G ReadFile_kin(cPar.file_kin, cPar.indicator_idv, cPar.mapID2num, @@ -3045,7 +3045,7 @@ void GEMMA::BatchRun(PARAM &cPar) { validate_K(G); } else { - cPar.ReadGenotypes(UtX, G, true); + cPar.ReadBIMBAMGenotypes(UtX, G, true); } // eigen-decomposition and calculate trace_G diff --git a/src/gemma_io.cpp b/src/gemma_io.cpp index b1bc368..6d33e00 100644 --- a/src/gemma_io.cpp +++ b/src/gemma_io.cpp @@ -688,14 +688,14 @@ bool ReadFile_fam(const string &file_fam, vector<vector<int>> &indicator_pheno, The function returns a `bool` indicating success/failure, but the real outputs are these three modified reference parameters. */ -bool ReadFile_geno(const string &file_geno, const set<string> &setSnps, - const gsl_matrix *W, const vector<int> &indicator_idv, - vector<int> &indicator_snp, const double &maf_level, - const double &miss_level, const double &hwe_level, - const double &r2_level, map<string, string> &mapRS2chr, - map<string, long int> &mapRS2bp, - map<string, double> &mapRS2cM, vector<SNPINFO> &snpInfo, - size_t &ns_test) { +bool ReadFile_bimbam_geno(const string &file_geno, const set<string> &setSnps, + const gsl_matrix *W, const vector<int> &indicator_idv, + vector<int> &indicator_snp, const double &maf_level, + const double &miss_level, const double &hwe_level, + const double &r2_level, map<string, string> &mapRS2chr, + map<string, long int> &mapRS2bp, + map<string, double> &mapRS2cM, vector<SNPINFO> &snpInfo, + size_t &ns_test) { checkpoint("read-file-geno",file_geno); debug_msg("entered"); indicator_snp.clear(); diff --git a/src/gemma_io.h b/src/gemma_io.h index 116c1a5..139f379 100644 --- a/src/gemma_io.h +++ b/src/gemma_io.h @@ -59,7 +59,7 @@ bool ReadFile_pheno(const string &file_pheno, bool ReadFile_column(const string &file_pheno, vector<int> &indicator_idv, vector<double> &pheno, const int &p_column); -bool ReadFile_geno(const string &file_geno, const set<string> &setSnps, +bool ReadFile_bimbam_geno(const string &file_geno, const set<string> &setSnps, const gsl_matrix *W, const vector<int> &indicator_idv, vector<int> &indicator_snp, const double &maf_level, const double &miss_level, const double &hwe_level, diff --git a/src/param.cpp b/src/param.cpp index 9755579..7e7ff6e 100644 --- a/src/param.cpp +++ b/src/param.cpp @@ -318,9 +318,10 @@ void PARAM::ReadFiles(void) { trim_individuals(indicator_idv, ni_max); trim_individuals(indicator_cvt, ni_max); - if (ReadFile_geno(file_geno, setSnps, W2, indicator_idv, indicator_snp, - maf_level, miss_level, hwe_level, r2_level, mapRS2chr, - mapRS2bp, mapRS2cM, snpInfo, ns_test) == false) { + // The following reads the geno file to get the SNPs + if (is_bimbam && ReadFile_bimbam_geno(file_geno, setSnps, W2, indicator_idv, indicator_snp, + maf_level, miss_level, hwe_level, r2_level, mapRS2chr, + mapRS2bp, mapRS2cM, snpInfo, ns_test) == false) { error = true; return; } @@ -424,7 +425,7 @@ void PARAM::ReadFiles(void) { string file_name; size_t ns_test_tmp; while (!safeGetline(infile, file_name).eof()) { - if (ReadFile_geno(file_name, setSnps, W4, indicator_idv, indicator_snp, + if (ReadFile_bimbam_geno(file_name, setSnps, W4, indicator_idv, indicator_snp, maf_level, miss_level, hwe_level, r2_level, mapRS2chr, mapRS2bp, mapRS2cM, snpInfo, ns_test_tmp) == false) { error = true; @@ -1262,7 +1263,7 @@ void PARAM::PrintSummary() { return; } -void PARAM::ReadGenotypes(gsl_matrix *UtX, gsl_matrix *K, const bool calc_K) { +void PARAM::ReadBIMBAMGenotypes(gsl_matrix *UtX, gsl_matrix *K, const bool calc_K) { string file_str; if (!file_bfile.empty()) { @@ -1272,8 +1273,8 @@ void PARAM::ReadGenotypes(gsl_matrix *UtX, gsl_matrix *K, const bool calc_K) { error = true; } } else { - if (ReadFile_geno(file_geno, indicator_idv, indicator_snp, UtX, K, - calc_K) == false) { + // Read BIMBAM. Not supposed to be mdb + if (is_mdb || ReadFile_geno(file_geno, indicator_idv, indicator_snp, UtX, K, calc_K) == false) { error = true; } } diff --git a/src/param.h b/src/param.h index 8d2fba3..0b8ed97 100644 --- a/src/param.h +++ b/src/param.h @@ -338,7 +338,7 @@ public: void CheckParam(); void CheckData(); void PrintSummary(); - void ReadGenotypes(gsl_matrix *UtX, gsl_matrix *K, const bool calc_K); + void ReadBIMBAMGenotypes(gsl_matrix *UtX, gsl_matrix *K, const bool calc_K); void CheckCvt(); void CopyCvt(gsl_matrix *W); void CopyA(size_t flag, gsl_matrix *A); @@ -372,4 +372,6 @@ public: size_t GetabIndex(const size_t a, const size_t b, const size_t n_cvt); +#define is_bimbam (!is_mdb) + #endif |