diff options
| -rw-r--r-- | src/gemma_io.cpp | 27 | ||||
| -rw-r--r-- | src/gemma_io.h | 8 | ||||
| -rw-r--r-- | src/lmm.h | 5 | ||||
| -rw-r--r-- | src/param.cpp | 2 |
4 files changed, 34 insertions, 8 deletions
diff --git a/src/gemma_io.cpp b/src/gemma_io.cpp index e183e3f..3e290f4 100644 --- a/src/gemma_io.cpp +++ b/src/gemma_io.cpp @@ -1478,8 +1478,27 @@ void ReadFile_eigenD(const string &file_kd, bool &error, gsl_vector *eval) { // Read bimbam mean genotype file and calculate kinship matrix. -gsl_matrix *mdb_calc_kin(const string file_geno, bool is_loco, const int k_mode) { +gsl_matrix *mdb_calc_kin(const string file_geno, const int k_mode, const string loco) { checkpoint("mdb-kin",file_geno); + bool is_loco = !loco.empty(); + + // Convert loco string to what we use in the chrpos index + uint8_t loco_chr = 0; + if (is_loco) { + if (loco == "X") { + loco_chr = CHR_X; + } else if (loco == "Y") { + loco_chr = CHR_Y; + } else if (loco == "M") { + loco_chr = CHR_M; + } else { + try { + loco_chr = static_cast<uint8_t>(stoi(loco)); + } catch (...) { + loco_chr = 0; + } + } + } /* Create and open the LMDB environment: */ auto env = lmdb::env::create(); @@ -1552,6 +1571,12 @@ gsl_matrix *mdb_calc_kin(const string file_geno, bool is_loco, const int k_mode) if (indicator_snp[t] == 0) continue; */ + const uint8_t* data = reinterpret_cast<const uint8_t*>(key.data()); + auto chr = static_cast<uint8_t>(data[1]); + + if (is_loco && loco_chr == chr) + continue; + #if !defined NDEBUG size_t num_floats = value.size() / sizeof(float); assert(num_floats == ni_total); diff --git a/src/gemma_io.h b/src/gemma_io.h index 4e9f4c1..06cd0ee 100644 --- a/src/gemma_io.h +++ b/src/gemma_io.h @@ -26,12 +26,17 @@ #include <algorithm> #include <map> #include <vector> +#include <cstdint> #include "gzstream.h" #include "param.h" #define tab(col) ( col ? "\t" : "") +constexpr uint8_t CHR_X = 'X'; +constexpr uint8_t CHR_Y = 'Y'; +constexpr uint8_t CHR_M = 'M'; + using namespace std; void ProgressBar(string str, double p, double total, double ratio = -1.0); @@ -87,7 +92,8 @@ void ReadFile_mk(const string &file_mk, vector<int> &indicator_idv, void ReadFile_eigenU(const string &file_u, bool &error, gsl_matrix *U); void ReadFile_eigenD(const string &file_d, bool &error, gsl_vector *eval); -gsl_matrix *mdb_calc_kin(const string file_geno, const bool is_loco, const int mode); +gsl_matrix *mdb_calc_kin(const string file_geno, const int k_mode, const string loco); + 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, diff --git a/src/lmm.h b/src/lmm.h index 8a65f8d..da5ad21 100644 --- a/src/lmm.h +++ b/src/lmm.h @@ -27,7 +27,6 @@ #include "param.h" #include <functional> #include <tuple> -#include <cstdint> using namespace std; @@ -46,10 +45,6 @@ public: }; -constexpr uint8_t CHR_X = 'X'; -constexpr uint8_t CHR_Y = 'Y'; -constexpr uint8_t CHR_M = 'M'; - // typedef tuple< string, uint16_t, uint32_t, uint32_t > MarkerInfo; // name, chr, pos, line struct MarkerInfo { string name; diff --git a/src/param.cpp b/src/param.cpp index 516dfd8..96a9cd2 100644 --- a/src/param.cpp +++ b/src/param.cpp @@ -1303,7 +1303,7 @@ void PARAM::ReadBIMBAMGenotypes(gsl_matrix *UtX, gsl_matrix *K, const bool calc_ } gsl_matrix *PARAM::MdbCalcKin() { - return mdb_calc_kin(file_geno, is_loco(), a_mode - 20); + return mdb_calc_kin(file_geno, a_mode - 20, loco); } void PARAM::CalcKin(gsl_matrix *matrix_kin) { |