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) {
|
