/*
Genome-wide Efficient Mixed Model Association (GEMMA)
Copyright © 2011-2017, Xiang Zhou
Copyright © 2017, Peter Carbonetto
Copyright © 2017, Pjotr Prins
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <iostream>
#include <iomanip>
#include <string>
#include <algorithm>
#include <cmath>
#include <cstring>
#include <fstream>
#include <sys/stat.h>
#include "gsl/gsl_blas.h"
#include "gsl/gsl_linalg.h"
#include "gsl/gsl_matrix.h"
#include "gsl/gsl_matrix.h"
#include "gsl/gsl_randist.h"
#include "gsl/gsl_vector.h"
#include "eigenlib.h"
#include "gemma.h"
#include "gemma_io.h"
#include "mathfunc.h"
#include "param.h"
#include "fastblas.h"
using namespace std;
// ---- Helper functions which do not use the PARAM scope
// Calculate the SNP sets as used in LOCO from mapchr which was filled
// from the annotation file. Returns ksnps (used for K) and gwasnps (used for
// GWA). Both are complimentary with each other and subsets of setSnps.
void LOCO_set_Snps(set<string> &ksnps, set<string> &gwasnps,
const map<string, string> mapchr, const string loco) {
enforce_msg(ksnps.size() == 0, "make sure knsps is not initialized twice");
enforce_msg(gwasnps.size() == 0,
"make sure gwasnps is not initialized twice");
for (auto &kv : mapchr) {
auto snp = kv.first;
auto chr = kv.second;
if (chr != loco) {
ksnps.insert(snp);
} else {
gwasnps.insert(snp);
}
}
}
// Trim #individuals to size which is used to write tests that run faster
//
// Note it actually trims the number of functional individuals
// (indicator_idv[x] == 1). This should match indicator_cvt etc. If
// this gives problems with certain sets we can simply trim to size.
void trim_individuals(vector<int> &idvs, size_t ni_max) {
if (ni_max) {
size_t count = 0;
for (auto ind = idvs.begin(); ind != idvs.end(); ++ind) {
if (*ind)
count++;
if (count >= ni_max)
break;
}
if (count != idvs.size()) {
if (is_debug_mode())
cout << "**** TEST MODE: trim individuals from " << idvs.size()
<< " to " << count << endl;
idvs.resize(count);
}
}
}
// ---- PARAM class implementation
PARAM::PARAM(void)
: a_mode(0), k_mode(1), d_pace(DEFAULT_PACE),
file_out("result"), path_out("./output/"), miss_level(0.05),
maf_level(0.01), hwe_level(0), r2_level(0.9999), l_min(1e-5), l_max(1e5),
n_region(10), p_nr(0.001), em_prec(0.0001), nr_prec(0.0001),
em_iter(10000), nr_iter(100), crt(0), pheno_mean(0), noconstrain(false),
h_min(-1), h_max(-1), h_scale(-1), rho_min(0.0), rho_max(1.0),
rho_scale(-1), logp_min(0.0), logp_max(0.0), logp_scale(-1), h_ngrid(10),
rho_ngrid(10), s_min(0), s_max(300), w_step(100000), s_step(1000000),
r_pace(10), w_pace(1000), n_accept(0), n_mh(10), geo_mean(2000.0),
randseed(-1), window_cm(0), window_bp(0), window_ns(0), n_block(200),
error(false), ni_subsample(0), n_cvt(1), n_cat(0), n_vc(1),
time_total(0.0), time_G(0.0), time_eigen(0.0), time_UtX(0.0),
time_UtZ(0.0), time_opt(0.0), time_Omega(0.0) {}
PARAM::~PARAM() {
if (gsl_r)
gsl_rng_free(gsl_r);
}
// Read files: obtain ns_total, ng_total, ns_test, ni_test.
void PARAM::ReadFiles(void) {
string file_str;
// Read cat file.
if (!file_mcat.empty()) {
if (ReadFile_mcat(file_mcat, mapRS2cat, n_vc) == false) {
error = true;
}
} else if (!file_cat.empty()) {
if (ReadFile_cat(file_cat, mapRS2cat, n_vc) == false) {
error = true;
}
}
// Read snp weight files.
if (!file_wcat.empty()) {
if (ReadFile_wsnp(file_wcat, n_vc, mapRS2wcat) == false) {
error = true;
}
}
if (!file_wsnp.empty()) {
if (ReadFile_wsnp(file_wsnp, mapRS2wsnp) == false) {
error = true;
}
}
// Count number of kinship files.
if (!file_mk.empty()) {
if (CountFileLines(file_mk, n_vc) == false) {
error = true;
}
}
// Read SNP set.
if (!file_snps.empty()) {
if (ReadFile_snps(file_snps, setSnps) == false) {
error = true;
}
} else {
setSnps.clear();
}
// Read KSNP set.
if (!file_ksnps.empty()) {
if (ReadFile_snps(file_ksnps, setKSnps) == false) {
error = true;
}
} else {
setKSnps.clear();
}
// For prediction.
if (!file_epm.empty()) {
if (ReadFile_est(file_epm, est_column, mapRS2est) == false) {
error = true;
}
if (!file_bfile.empty()) {
file_str = file_bfile + ".bim";
if (ReadFile_bim(file_str, snpInfo) == false) {
error = true;
}
file_str = file_bfile + ".fam";
if (ReadFile_fam(file_str, indicator_pheno, pheno, mapID2num, p_column) ==
false) {
error = true;
}
}
if (!file_geno.empty()) {
if (ReadFile_pheno(file_pheno, indicator_pheno, pheno, p_column) ==
false) {
error = true;
}
if (CountFileLines(file_geno, ns_total) == false) {
error = true;
}
}
if (!file_ebv.empty()) {
if (ReadFile_column(file_ebv, indicator_bv, vec_bv, 1) == false) {
error = true;
}
}
if (!file_log.empty()) {
if (ReadFile_log(file_log, pheno_mean) == false) {
error = true;
}
}
// Convert indicator_pheno to indicator_idv.
int k = 1;
for (size_t i = 0; i < indicator_pheno.size(); i++) {
k = 1;
for (size_t j = 0; j < indicator_pheno[i].size(); j++) {
if (indicator_pheno[i][j] == 0) {
k = 0;
}
}
indicator_idv.push_back(k);
}
ns_test = 0;
return;
}
// Read covariates before the genotype files.
if (!file_cvt.empty()) {
if (ReadFile_cvt(file_cvt, indicator_cvt, cvt, n_cvt) == false) {
error = true;
}
if ((indicator_cvt).size() == 0) {
n_cvt = 1;
}
} else {
n_cvt = 1;
}
trim_individuals(indicator_cvt, ni_max);
if (!file_gxe.empty()) {
if (ReadFile_column(file_gxe, indicator_gxe, gxe, 1) == false) {
error = true;
}
}
if (!file_weight.empty()) {
if (ReadFile_column(file_weight, indicator_weight, weight, 1) == false) {
error = true;
}
}
trim_individuals(indicator_idv, ni_max);
// Read genotype and phenotype file for PLINK format.
if (!file_bfile.empty()) {
file_str = file_bfile + ".bim";
snpInfo.clear();
if (ReadFile_bim(file_str, snpInfo) == false) {
error = true;
}
// If both fam file and pheno files are used, use
// phenotypes inside the pheno file.
if (!file_pheno.empty()) {
// Phenotype file before genotype file.
if (ReadFile_pheno(file_pheno, indicator_pheno, pheno, p_column) ==
false) {
error = true;
}
} else {
file_str = file_bfile + ".fam";
if (ReadFile_fam(file_str, indicator_pheno, pheno, mapID2num, p_column) ==
false) {
error = true;
}
}
// Post-process covariates and phenotypes, obtain
// ni_test, save all useful covariates.
ProcessCvtPhen();
// Obtain covariate matrix.
auto W1 = gsl_matrix_safe_alloc(ni_test, n_cvt);
CopyCvt(W1);
file_str = file_bfile + ".bed";
if (ReadFile_bed(file_str, setSnps, W1, indicator_idv, indicator_snp,
snpInfo, maf_level, miss_level, hwe_level, r2_level,
ns_test) == false) {
error = true;
}
gsl_matrix_free(W1);
ns_total = indicator_snp.size();
}
// Read genotype and phenotype file for BIMBAM format.
if (!file_geno.empty()) {
// Annotation file before genotype file.
if (!file_anno.empty()) {
if (ReadFile_anno(file_anno, mapRS2chr, mapRS2bp, mapRS2cM) == false) {
error = true;
}
}
// Phenotype file before genotype file.
if (ReadFile_pheno(file_pheno, indicator_pheno, pheno, p_column) == false) {
error = true;
}
// Post-process covariates and phenotypes, obtain
// ni_test, save all useful covariates.
ProcessCvtPhen();
// Obtain covariate matrix.
auto W2 = gsl_matrix_safe_alloc(ni_test, n_cvt);
CopyCvt(W2);
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) {
error = true;
return;
}
gsl_matrix_free(W2);
ns_total = indicator_snp.size();
}
// Read genotype file for multiple PLINK files.
if (!file_mbfile.empty()) {
igzstream infile(file_mbfile.c_str(), igzstream::in);
enforce_msg(infile,"fail to open mbfile file");
string file_name;
size_t t = 0, ns_test_tmp = 0;
gsl_matrix *W3 = NULL;
while (!safeGetline(infile, file_name).eof()) {
file_str = file_name + ".bim";
if (ReadFile_bim(file_str, snpInfo) == false) {
error = true;
}
if (t == 0) {
// If both fam file and pheno files are used, use
// phenotypes inside the pheno file.
if (!file_pheno.empty()) {
// Phenotype file before genotype file.
if (ReadFile_pheno(file_pheno, indicator_pheno, pheno, p_column) ==
false) {
error = true;
}
} else {
file_str = file_name + ".fam";
if (ReadFile_fam(file_str, indicator_pheno, pheno, mapID2num,
p_column) == false) {
error = true;
}
}
// Post-process covariates and phenotypes, obtain
// ni_test, save all useful covariates.
ProcessCvtPhen();
// Obtain covariate matrix.
W3 = gsl_matrix_safe_alloc(ni_test, n_cvt);
CopyCvt(W3);
}
file_str = file_name + ".bed";
if (ReadFile_bed(file_str, setSnps, W3, indicator_idv, indicator_snp,
snpInfo, maf_level, miss_level, hwe_level, r2_level,
ns_test_tmp) == false) {
error = true;
}
mindicator_snp.push_back(indicator_snp);
msnpInfo.push_back(snpInfo);
ns_test += ns_test_tmp;
ns_total += indicator_snp.size();
t++;
}
if (W3) gsl_matrix_free(W3);
infile.close();
infile.clear();
}
// Read genotype and phenotype file for multiple BIMBAM files.
if (!file_mgeno.empty()) {
// Annotation file before genotype file.
if (!file_anno.empty()) {
if (ReadFile_anno(file_anno, mapRS2chr, mapRS2bp, mapRS2cM) == false) {
error = true;
}
}
// Phenotype file before genotype file.
if (ReadFile_pheno(file_pheno, indicator_pheno, pheno, p_column) == false) {
error = true;
}
// Post-process covariates and phenotypes, obtain ni_test,
// save all useful covariates.
ProcessCvtPhen();
// Obtain covariate matrix.
gsl_matrix *W4 = gsl_matrix_safe_alloc(ni_test, n_cvt);
CopyCvt(W4);
igzstream infile(file_mgeno.c_str(), igzstream::in);
if (!infile) {
cout << "error! fail to open mgeno file: " << file_mgeno << endl;
error = true;
return;
}
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,
maf_level, miss_level, hwe_level, r2_level, mapRS2chr,
mapRS2bp, mapRS2cM, snpInfo, ns_test_tmp) == false) {
error = true;
}
mindicator_snp.push_back(indicator_snp);
msnpInfo.push_back(snpInfo);
ns_test += ns_test_tmp;
ns_total += indicator_snp.size();
}
gsl_matrix_free(W4);
infile.close();
infile.clear();
}
if (!file_gene.empty()) {
if (ReadFile_pheno(file_pheno, indicator_pheno, pheno, p_column) == false) {
error = true;
}
// Convert indicator_pheno to indicator_idv.
int k = 1;
for (size_t i = 0; i < indicator_pheno.size(); i++) {
k = 1;
for (size_t j = 0; j < indicator_pheno[i].size(); j++) {
if (indicator_pheno[i][j] == 0) {
k = 0;
}
}
indicator_idv.push_back(k);
}
// Post-process covariates and phenotypes, obtain
// ni_test, save all useful covariates.
ProcessCvtPhen();
// Obtain covariate matrix.
// gsl_matrix *W5 = gsl_matrix_alloc(ni_test, n_cvt);
// CopyCvt(W5);
if (ReadFile_gene(file_gene, vec_read, snpInfo, ng_total) == false) {
error = true;
}
}
// Read is after gene file.
if (!file_read.empty()) {
if (ReadFile_column(file_read, indicator_read, vec_read, 1) == false) {
error = true;
}
ni_test = 0;
for (vector<int>::size_type i = 0; i < (indicator_idv).size(); ++i) {
indicator_idv[i] *= indicator_read[i];
ni_test += indicator_idv[i];
}
enforce_msg(ni_test,"number of analyzed individuals equals 0.");
}
// For ridge prediction, read phenotype only.
if (file_geno.empty() && file_gene.empty() && !file_pheno.empty()) {
if (ReadFile_pheno(file_pheno, indicator_pheno, pheno, p_column) == false) {
error = true;
}
// Post-process covariates and phenotypes, obtain
// ni_test, save all useful covariates.
ProcessCvtPhen();
}
// Compute setKSnps when -loco is passed in
if (!loco.empty()) {
LOCO_set_Snps(setKSnps, setGWASnps, mapRS2chr, loco);
}
return;
}
void PARAM::CheckParam(void) {
struct stat fileInfo;
string str;
// Check parameters.
if (k_mode != 1 && k_mode != 2) {
cout << "error! unknown kinship/relatedness input mode: " << k_mode << endl;
error = true;
}
if (a_mode != 1 && a_mode != 2 && a_mode != 3 && a_mode != 4 && a_mode != 5 &&
a_mode != M_LMM9 &&
a_mode != 11 && a_mode != 12 && a_mode != 13 && a_mode != 14 &&
a_mode != 15 && a_mode != 21 && a_mode != 22 && a_mode != 25 &&
a_mode != 26 && a_mode != 27 && a_mode != 28 && a_mode != 31 &&
a_mode != 41 && a_mode != 42 && a_mode != 43 && a_mode != 51 &&
a_mode != 52 && a_mode != 53 && a_mode != 54 && a_mode != 61 &&
a_mode != 62 && a_mode != 63 && a_mode != 66 && a_mode != 67 &&
a_mode != 71) {
cout << "error! unknown analysis mode: " << a_mode
<< ". make sure -gk or -eigen or -lmm or -bslmm -predict or "
<< "-calccov is specified correctly." << endl;
error = true;
}
if (miss_level > 1) {
cout << "error! missing level needs to be between 0 and 1. "
<< "current value = " << miss_level << endl;
error = true;
}
if (maf_level > 0.5) {
cout << "error! maf level needs to be between 0 and 0.5. "
<< "current value = " << maf_level << endl;
error = true;
}
if (hwe_level > 1) {
cout << "error! hwe level needs to be between 0 and 1. "
<< "current value = " << hwe_level << endl;
error = true;
}
if (r2_level > 1) {
cout << "error! r2 level needs to be between 0 and 1. "
<< "current value = " << r2_level << endl;
error = true;
}
if (l_max < l_min) {
cout << "error! maximum lambda value must be larger than the "
<< "minimal value. current values = " << l_max << " and " << l_min
<< endl;
error = true;
}
if (h_max < h_min) {
cout << "error! maximum h value must be larger than the minimal "
<< "value. current values = " << h_max << " and " << h_min << endl;
error = true;
}
if (s_max < s_min) {
cout << "error! maximum s value must be larger than the minimal "
<< "value. current values = " << s_max << " and " << s_min << endl;
error = true;
}
if (rho_max < rho_min) {
cout << "error! maximum rho value must be larger than the"
<< "minimal value. current values = " << rho_max << " and " << rho_min
<< endl;
error = true;
}
if (logp_max < logp_min) {
cout << "error! maximum logp value must be larger than the "
<< "minimal value. current values = " << logp_max / log(10) << " and "
<< logp_min / log(10) << endl;
error = true;
}
if (h_max > 1) {
cout << "error! h values must be bewtween 0 and 1. current "
<< "values = " << h_max << " and " << h_min << endl;
error = true;
}
if (rho_max > 1) {
cout << "error! rho values must be between 0 and 1. current "
<< "values = " << rho_max << " and " << rho_min << endl;
error = true;
}
if (logp_max > 0) {
cout << "error! maximum logp value must be smaller than 0. "
<< "current values = " << logp_max / log(10) << " and "
<< logp_min / log(10) << endl;
error = true;
}
if (l_max < l_min) {
cout << "error! maximum lambda value must be larger than the "
<< "minimal value. current values = " << l_max << " and " << l_min
<< endl;
error = true;
}
if (h_scale > 1.0) {
cout << "error! hscale value must be between 0 and 1. "
<< "current value = " << h_scale << endl;
error = true;
}
if (rho_scale > 1.0) {
cout << "error! rscale value must be between 0 and 1. "
<< "current value = " << rho_scale << endl;
error = true;
}
if (logp_scale > 1.0) {
cout << "error! pscale value must be between 0 and 1. "
<< "current value = " << logp_scale << endl;
error = true;
}
if (rho_max == 1 && rho_min == 1 && a_mode == 12) {
cout << "error! ridge regression does not support a rho "
<< "parameter. current values = " << rho_max << " and " << rho_min
<< endl;
error = true;
}
if (window_cm < 0) {
cout << "error! windowcm values must be non-negative. "
<< "current values = " << window_cm << endl;
error = true;
}
if (window_cm == 0 && window_bp == 0 && window_ns == 0) {
window_bp = 1000000;
}
// Check p_column, and (no need to) sort p_column into
// ascending order.
if (p_column.size() == 0) {
p_column.push_back(1);
} else {
for (size_t i = 0; i < p_column.size(); i++) {
for (size_t j = 0; j < i; j++) {
if (p_column[i] == p_column[j]) {
cout << "error! identical phenotype "
<< "columns: " << p_column[i] << endl;
error = true;
}
}
}
}
n_ph = p_column.size();
// Only LMM option (and one prediction option) can deal with
// multiple phenotypes and no gene expression files.
if (n_ph > 1 && a_mode != 1 && a_mode != 2 && a_mode != 3 && a_mode != 4 &&
a_mode != 9 && a_mode != 43) {
cout << "error! the current analysis mode " << a_mode
<< " can not deal with multiple phenotypes." << endl;
error = true;
}
if (n_ph > 1 && !file_gene.empty()) {
cout << "error! multiple phenotype analysis option not "
<< "allowed with gene expression files. " << endl;
error = true;
}
if (p_nr > 1) {
cout << "error! pnr value must be between 0 and 1. current value = " << p_nr
<< endl;
error = true;
}
// check est_column
if (est_column.size() == 0) {
if (file_ebv.empty()) {
est_column.push_back(2);
est_column.push_back(5);
est_column.push_back(6);
est_column.push_back(7);
} else {
est_column.push_back(2);
est_column.push_back(0);
est_column.push_back(6);
est_column.push_back(7);
}
}
if (est_column.size() != 4) {
cout << "error! -en not followed by four numbers. current number = "
<< est_column.size() << endl;
error = true;
}
if (est_column[0] == 0) {
cout << "error! -en rs column can not be zero. current number = "
<< est_column.size() << endl;
error = true;
}
// Check if files are compatible with each other, and if files exist.
if (!file_bfile.empty()) {
str = file_bfile + ".bim";
if (stat(str.c_str(), &fileInfo) == -1) {
cout << "error! fail to open .bim file: " << str << endl;
error = true;
}
str = file_bfile + ".bed";
if (stat(str.c_str(), &fileInfo) == -1) {
cout << "error! fail to open .bed file: " << str << endl;
error = true;
}
str = file_bfile + ".fam";
if (stat(str.c_str(), &fileInfo) == -1) {
cout << "error! fail to open .fam file: " << str << endl;
error = true;
}
}
if ((!file_geno.empty() || !file_gene.empty())) {
str = file_pheno;
if (stat(str.c_str(), &fileInfo) == -1) {
cout << "error! fail to open phenotype file: " << str << endl;
error = true;
}
}
str = file_geno;
if (!str.empty() && stat(str.c_str(), &fileInfo) == -1) {
cout << "error! fail to open mean genotype file: " << str << endl;
error = true;
}
str = file_gene;
if (!str.empty() && stat(str.c_str(), &fileInfo) == -1) {
cout << "error! fail to open gene expression file: " << str << endl;
error = true;
}
str = file_cat;
if (!str.empty() && stat(str.c_str(), &fileInfo) == -1) {
cout << "error! fail to open category file: " << str << endl;
error = true;
}
str = file_mcat;
if (!str.empty() && stat(str.c_str(), &fileInfo) == -1) {
cout << "error! fail to open mcategory file: " << str << endl;
error = true;
}
str = file_beta;
if (!str.empty() && stat(str.c_str(), &fileInfo) == -1) {
cout << "error! fail to open beta file: " << str << endl;
error = true;
}
str = file_cor;
if (!str.empty() && stat(str.c_str(), &fileInfo) == -1) {
cout << "error! fail to open correlation file: " << str << endl;
error = true;
}
if (!file_study.empty()) {
str = file_study + ".Vq.txt";
if (stat(str.c_str(), &fileInfo) == -1) {
cout << "error! fail to open .Vq.txt file: " << str << endl;
error = true;
}
str = file_study + ".q.txt";
if (stat(str.c_str(), &fileInfo) == -1) {
cout << "error! fail to open .q.txt file: " << str << endl;
error = true;
}
str = file_study + ".size.txt";
if (stat(str.c_str(), &fileInfo) == -1) {
cout << "error! fail to open .size.txt file: " << str << endl;
error = true;
}
}
if (!file_ref.empty()) {
str = file_ref + ".S.txt";
if (stat(str.c_str(), &fileInfo) == -1) {
cout << "error! fail to open .S.txt file: " << str << endl;
error = true;
}
str = file_ref + ".size.txt";
if (stat(str.c_str(), &fileInfo) == -1) {
cout << "error! fail to open .size.txt file: " << str << endl;
error = true;
}
}
str = file_mstudy;
if (!str.empty() && stat(str.c_str(), &fileInfo) == -1) {
cout << "error! fail to open mstudy file: " << str << endl;
error = true;
}
str = file_mref;
if (!str.empty() && stat(str.c_str(), &fileInfo) == -1) {
cout << "error! fail to open mref file: " << str << endl;
error = true;
}
str = file_mgeno;
if (!str.empty() && stat(str.c_str(), &fileInfo) == -1) {
cout << "error! fail to open mgeno file: " << str << endl;
error = true;
}
str = file_mbfile;
if (!str.empty() && stat(str.c_str(), &fileInfo) == -1) {
cout << "error! fail to open mbfile file: " << str << endl;
error = true;
}
size_t flag = 0;
if (!file_bfile.empty()) {
flag++;
}
if (!file_geno.empty()) {
flag++;
}
if (!file_gene.empty()) {
flag++;
}
// Always set up random environment.
gsl_rng_env_setup(); // sets gsl_rng_default_seed
const gsl_rng_type *T = gsl_rng_default; // pick up environment GSL_RNG_SEED
if (randseed >= 0)
gsl_rng_default_seed = randseed; // CLI option used
else if (gsl_rng_default_seed == 0) { // by default we will randomize the seed
time_t rawtime;
time(&rawtime);
tm *ptm = gmtime(&rawtime);
gsl_rng_default_seed =
(unsigned)(ptm->tm_hour % 24 * 3600 + ptm->tm_min * 60 + ptm->tm_sec);
}
gsl_r = gsl_rng_alloc(T);
if (is_debug_mode()) {
printf ("GSL random generator type: %s; ", gsl_rng_name (gsl_r));
printf ("seed = %lu (option %li); ", gsl_rng_default_seed, randseed);
printf ("first value = %lu\n", gsl_rng_get (gsl_r));
}
if (flag != 1 && a_mode != 15 && a_mode != 27 && a_mode != 28 &&
a_mode != 43 && a_mode != 5 && a_mode != 61 && a_mode != 62 &&
a_mode != 63 && a_mode != 66 && a_mode != 67) {
cout << "error! either plink binary files, or bimbam mean"
<< "genotype files, or gene expression files are required." << endl;
error = true;
}
if (file_pheno.empty() && (a_mode == 43 || a_mode == 5)) {
cout << "error! phenotype file is required." << endl;
error = true;
}
if (a_mode == 61 || a_mode == 62) {
if (!file_beta.empty()) {
if (file_mbfile.empty() && file_bfile.empty() && file_mgeno.empty() &&
file_geno.empty() && file_mref.empty() && file_ref.empty()) {
cout << "error! missing genotype file or ref/mref file." << endl;
error = true;
}
} else if (!file_pheno.empty()) {
if (file_kin.empty() && (file_ku.empty() || file_kd.empty()) &&
file_mk.empty()) {
cout << "error! missing relatedness file. " << endl;
error = true;
}
} else if ((file_mstudy.empty() && file_study.empty()) ||
(file_mref.empty() && file_ref.empty())) {
cout << "error! either beta file, or phenotype files or "
<< "study/ref mstudy/mref files are required." << endl;
error = true;
}
}
if (a_mode == 63) {
if (file_kin.empty() && (file_ku.empty() || file_kd.empty()) &&
file_mk.empty()) {
cout << "error! missing relatedness file. " << endl;
error = true;
}
if (file_pheno.empty()) {
cout << "error! missing phenotype file." << endl;
error = true;
}
}
if (a_mode == 66 || a_mode == 67) {
if (file_beta.empty() || (file_mbfile.empty() && file_bfile.empty() &&
file_mgeno.empty() && file_geno.empty())) {
cout << "error! missing beta file or genotype file." << endl;
error = true;
}
}
if (!file_epm.empty() && file_bfile.empty() && file_geno.empty()) {
cout << "error! estimated parameter file also requires genotype "
<< "file." << endl;
error = true;
}
if (!file_ebv.empty() && file_kin.empty()) {
cout << "error! estimated breeding value file also requires "
<< "relatedness file." << endl;
error = true;
}
if (!file_log.empty() && pheno_mean != 0) {
cout << "error! either log file or mu value can be provide." << endl;
error = true;
}
enforce_fexists(file_snps, "open file");
enforce_fexists(file_ksnps, "open file");
enforce_fexists(file_gwasnps, "open file");
enforce_fexists(file_anno, "open file");
if (!loco.empty()) {
enforce_msg((a_mode >= 1 && a_mode <= 4) || a_mode == 9 || a_mode == 21 || a_mode == 22,
"LOCO only works with LMM and K");
// enforce_msg(file_bfile.empty(), "LOCO does not work with PLink (yet)");
enforce_msg(file_gxe.empty(), "LOCO does not support GXE (yet)");
enforce_msg(!file_anno.empty(),
"LOCO requires annotation file (-a switch)");
enforce_msg(file_ksnps.empty(), "LOCO does not allow -ksnps switch");
enforce_msg(file_gwasnps.empty(), "LOCO does not allow -gwasnps switch");
}
enforce_fexists(file_kin, "open file");
enforce_fexists(file_mk, "open file");
enforce_fexists(file_cvt, "open file");
enforce_fexists(file_gxe, "open file");
enforce_fexists(file_log, "open file");
enforce_fexists(file_weight, "open file");
enforce_fexists(file_epm, "open file");
enforce_fexists(file_ebv, "open file");
enforce_fexists(file_read, "open file");
// Check if files are compatible with analysis mode.
if (k_mode == 2 && !file_geno.empty()) {
cout << "error! use \"-km 1\" when using bimbam mean genotype "
<< "file. " << endl;
error = true;
}
if ((a_mode == 1 || a_mode == 2 || a_mode == 3 || a_mode == 4 ||
a_mode == 5 || a_mode == 9 || a_mode == 31) &&
(file_kin.empty() && (file_ku.empty() || file_kd.empty()))) {
cout << "error! missing relatedness file. " << endl;
error = true;
}
if ((a_mode == 43) && file_kin.empty()) {
cout << "error! missing relatedness file. -predict option requires "
<< "-k option to provide a relatedness file." << endl;
error = true;
}
if ((a_mode == 11 || a_mode == 12 || a_mode == 13 || a_mode == 14 ||
a_mode == 16) &&
!file_cvt.empty()) {
cout << "error! -bslmm option does not support covariates files." << endl;
error = true;
}
if (a_mode == 41 || a_mode == 42) {
if (!file_cvt.empty()) {
cout << "error! -predict option does not support "
<< "covariates files." << endl;
error = true;
}
if (file_epm.empty()) {
cout << "error! -predict option requires estimated "
<< "parameter files." << endl;
error = true;
}
}
if (file_beta.empty() && (a_mode == 27 || a_mode == 28)) {
cout << "error! beta effects file is required." << endl;
error = true;
}
return;
}
void PARAM::CheckData(void) {
if ((a_mode == 66 || a_mode == 67) && (v_pve.size() != n_vc)) {
cout << "error! the number of pve estimates does not equal to "
<< "the number of categories in the cat file:" << v_pve.size() << " "
<< n_vc << endl;
error = true;
}
if ((indicator_cvt).size() != 0 &&
(indicator_cvt).size() != (indicator_idv).size()) {
error = true;
cout << "error! number of rows in the covariates file do not "
<< "match the number of individuals. " << indicator_cvt.size() << endl;
return;
}
if ((indicator_gxe).size() != 0 &&
(indicator_gxe).size() != (indicator_idv).size()) {
error = true;
cout << "error! number of rows in the gxe file do not match the number "
<< "of individuals. " << endl;
return;
}
if ((indicator_weight).size() != 0 &&
(indicator_weight).size() != (indicator_idv).size()) {
error = true;
cout << "error! number of rows in the weight file do not match "
<< "the number of individuals. " << endl;
return;
}
if ((indicator_read).size() != 0 &&
(indicator_read).size() != (indicator_idv).size()) {
error = true;
cout << "error! number of rows in the total read file do not "
<< "match the number of individuals. " << endl;
return;
}
// Calculate ni_total and ni_test, and set indicator_idv to 0
// whenever indicator_cvt=0, and calculate np_obs and np_miss.
ni_total = (indicator_idv).size();
ni_test = 0;
for (vector<int>::size_type i = 0; i < (indicator_idv).size(); ++i) {
if (indicator_idv[i] == 0) {
continue;
}
ni_test++;
}
ni_cvt = 0;
for (size_t i = 0; i < indicator_cvt.size(); i++) {
if (indicator_cvt[i] == 0) {
continue;
}
ni_cvt++;
}
np_obs = 0;
np_miss = 0;
for (size_t i = 0; i < indicator_pheno.size(); i++) {
if (indicator_cvt.size() != 0) {
if (indicator_cvt[i] == 0) {
continue;
}
}
if (indicator_gxe.size() != 0) {
if (indicator_gxe[i] == 0) {
continue;
}
}
if (indicator_weight.size() != 0) {
if (indicator_weight[i] == 0) {
continue;
}
}
for (size_t j = 0; j < indicator_pheno[i].size(); j++) {
if (indicator_pheno[i][j] == 0) {
np_miss++;
} else {
np_obs++;
}
}
}
enforce_msg(ni_test,"number of analyzed individuals equals 0.");
if (ni_test == 0 && file_cor.empty() && file_mstudy.empty() &&
file_study.empty() && file_beta.empty() && file_bf.empty()) {
error = true;
cout << "error! number of analyzed individuals equals 0. " << endl;
return;
}
if (a_mode == 43) {
if (ni_cvt == ni_test) {
error = true;
cout << "error! no individual has missing "
<< "phenotypes." << endl;
return;
}
if ((np_obs + np_miss) != (ni_cvt * n_ph)) {
error = true;
cout << "error! number of phenotypes do not match the "
<< "summation of missing and observed phenotypes." << endl;
return;
}
}
// Output some information.
if (file_cor.empty() && file_mstudy.empty() && file_study.empty() &&
a_mode != 15 && a_mode != 27 && a_mode != 28) {
cout << "## number of total individuals = " << ni_total << endl;
if (a_mode == 43) {
cout << "## number of analyzed individuals = " << ni_cvt << endl;
cout << "## number of individuals with full phenotypes = " << ni_test
<< endl;
} else {
cout << "## number of analyzed individuals = " << ni_test << endl;
}
cout << "## number of covariates = " << n_cvt << endl;
cout << "## number of phenotypes = " << n_ph << endl;
if (a_mode == 43) {
cout << "## number of observed data = " << np_obs << endl;
cout << "## number of missing data = " << np_miss << endl;
}
if (!file_gene.empty()) {
cout << "## number of total genes = " << ng_total << endl;
} else if (file_epm.empty() && a_mode != 43 && a_mode != 5) {
if (!loco.empty())
cout << "## leave one chromosome out (LOCO) = " << setw(8) << loco << endl;
cout << "## number of total SNPs/var = " << setw(8) << ns_total << endl;
if (setSnps.size())
cout << "## number of considered SNPS = " << setw(8) << setSnps.size() << endl;
if (setKSnps.size())
cout << "## number of SNPS for K = " << setw(8) << setKSnps.size() << endl;
if (setGWASnps.size())
cout << "## number of SNPS for GWAS = " << setw(8) << setGWASnps.size() << endl;
cout << "## number of analyzed SNPs = " << setw(8) << ns_test << endl;
} else {
}
}
// Set d_pace to 1000 for gene expression.
if (!file_gene.empty() && d_pace == DEFAULT_PACE) {
d_pace = 1000;
}
// For case-control studies, count # cases and # controls.
int flag_cc = 0;
if (a_mode == 13) {
ni_case = 0;
ni_control = 0;
for (size_t i = 0; i < indicator_idv.size(); i++) {
if (indicator_idv[i] == 0) {
continue;
}
if (pheno[i][0] == 0) {
ni_control++;
} else if (pheno[i][0] == 1) {
ni_case++;
} else {
flag_cc = 1;
}
}
cout << "## number of cases = " << ni_case << endl;
cout << "## number of controls = " << ni_control << endl;
}
if (flag_cc == 1) {
cout << "Unexpected non-binary phenotypes for "
<< "case/control analysis. Use default (BSLMM) analysis instead."
<< endl;
a_mode = 11;
}
// Set parameters for BSLMM and check for predict.
if (a_mode == 11 || a_mode == 12 || a_mode == 13 || a_mode == 14) {
if (a_mode == 11) {
n_mh = 1;
}
if (logp_min == 0) {
logp_min = -1.0 * log((double)ns_test);
}
if (h_scale == -1) {
h_scale = min(1.0, 10.0 / sqrt((double)ni_test));
}
if (rho_scale == -1) {
rho_scale = min(1.0, 10.0 / sqrt((double)ni_test));
}
if (logp_scale == -1) {
logp_scale = min(1.0, 5.0 / sqrt((double)ni_test));
}
if (h_min == -1) {
h_min = 0.0;
}
if (h_max == -1) {
h_max = 1.0;
}
if (s_max > ns_test) {
s_max = ns_test;
cout << "s_max is re-set to the number of analyzed SNPs." << endl;
}
if (s_max < s_min) {
cout << "error! maximum s value must be larger than the "
<< "minimal value. current values = " << s_max << " and " << s_min
<< endl;
error = true;
}
} else if (a_mode == 41 || a_mode == 42) {
if (indicator_bv.size() != 0) {
if (indicator_idv.size() != indicator_bv.size()) {
cout << "error! number of rows in the "
<< "phenotype file does not match that in the "
<< "estimated breeding value file: " << indicator_idv.size()
<< "\t" << indicator_bv.size() << endl;
error = true;
} else {
size_t flag_bv = 0;
for (size_t i = 0; i < (indicator_bv).size(); ++i) {
if (indicator_idv[i] != indicator_bv[i]) {
flag_bv++;
}
}
if (flag_bv != 0) {
cout << "error! individuals with missing value in the "
<< "phenotype file does not match that in the "
<< "estimated breeding value file: " << flag_bv << endl;
error = true;
}
}
}
}
if (a_mode == 62 && !file_beta.empty() && mapRS2wcat.size() == 0) {
cout << "vc analysis with beta files requires -wcat file." << endl;
error = true;
}
if (a_mode == 67 && mapRS2wcat.size() == 0) {
cout << "ci analysis with beta files requires -wcat file." << endl;
error = true;
}
// File_mk needs to contain more than one line.
if (n_vc == 1 && !file_mk.empty()) {
cout << "error! -mk file should contain more than one line." << endl;
error = true;
}
return;
}
void PARAM::PrintSummary() {
if (n_ph == 1) {
cout << "pve estimate =" << pve_null << endl;
cout << "se(pve) =" << pve_se_null << endl;
} else {
}
return;
}
void PARAM::ReadGenotypes(gsl_matrix *UtX, gsl_matrix *K, const bool calc_K) {
string file_str;
if (!file_bfile.empty()) {
file_str = file_bfile + ".bed";
if (ReadFile_bed(file_str, indicator_idv, indicator_snp, UtX, K, calc_K) ==
false) {
error = true;
}
} else {
if (ReadFile_geno(file_geno, indicator_idv, indicator_snp, UtX, K,
calc_K) == false) {
error = true;
}
}
return;
}
void PARAM::ReadGenotypes(vector<vector<unsigned char>> &Xt, gsl_matrix *K,
const bool calc_K) {
string file_str;
if (!file_bfile.empty()) {
file_str = file_bfile + ".bed";
if (ReadFile_bed(file_str, indicator_idv, indicator_snp, Xt, K, calc_K,
ni_test, ns_test) == false) {
error = true;
}
} else {
if (ReadFile_geno(file_geno, indicator_idv, indicator_snp, Xt, K, calc_K,
ni_test, ns_test) == false) {
error = true;
}
}
return;
}
void PARAM::CalcKin(gsl_matrix *matrix_kin) {
string file_str;
gsl_matrix_set_zero(matrix_kin);
if (!file_bfile.empty()) {
file_str = file_bfile + ".bed";
// enforce_msg(loco.empty(), "FIXME: LOCO nyi");
if (PlinkKin(file_str, indicator_snp, a_mode - 20, d_pace, matrix_kin) ==
false) {
error = true;
}
} else {
file_str = file_geno;
if (BimbamKin(file_str, setKSnps, indicator_snp, a_mode - 20, d_pace,
matrix_kin, ni_max == 0) == false) {
error = true;
}
}
return;
}
// From an existing n by nd A and K matrices, compute the d by d S
// matrix (which is not necessary symmetric).
void compAKtoS(const gsl_matrix *A, const gsl_matrix *K, const size_t n_cvt,
gsl_matrix *S) {
size_t n_vc = S->size1, ni_test = A->size1;
double di, dj, tr_AK, sum_A, sum_K, s_A, s_K, sum_AK, tr_A, tr_K, d;
for (size_t i = 0; i < n_vc; i++) {
for (size_t j = 0; j < n_vc; j++) {
tr_AK = 0;
sum_A = 0;
sum_K = 0;
sum_AK = 0;
tr_A = 0;
tr_K = 0;
for (size_t l = 0; l < ni_test; l++) {
s_A = 0;
s_K = 0;
for (size_t k = 0; k < ni_test; k++) {
di = gsl_matrix_get(A, l, k + ni_test * i);
dj = gsl_matrix_get(K, l, k + ni_test * j);
s_A += di;
s_K += dj;
tr_AK += di * dj;
sum_A += di;
sum_K += dj;
if (l == k) {
tr_A += di;
tr_K += dj;
}
}
sum_AK += s_A * s_K;
}
sum_A /= (double)ni_test;
sum_K /= (double)ni_test;
sum_AK /= (double)ni_test;
tr_A -= sum_A;
tr_K -= sum_K;
d = tr_AK - 2 * sum_AK + sum_A * sum_K;
if (tr_A == 0 || tr_K == 0) {
d = 0;
} else {
assert((tr_A * tr_K) - 1 != 0);
assert(ni_test - n_cvt != 0);
d = d / (tr_A * tr_K) - 1 / (double)(ni_test - n_cvt);
}
gsl_matrix_set(S, i, j, d);
}
}
return;
}
/*
// Copied from lmm.cpp; is used in the following function compKtoV
// map a number 1..(n_cvt+2) to an index between 0 and [(n_cvt+2)*2+(n_cvt+2)]/2-1
// or 1..cols to 0..(cols*2+cols)/2-1.
For a 3x3 matrix the following index gets returned to CalcPab:
CalcPab
* 1,1:0
* 1,2:1
* 1,3:2
* 2,2:5
* 2,3:6
* 3,3:9
which is really the iteration moving forward along the diagonal and
items to the right of it.
*/
size_t GetabIndex(const size_t a, const size_t b, const size_t n_cvt) {
auto cols = n_cvt + 2;
enforce_msg(a<=cols && b<=cols,"GetabIndex problem");
size_t a1 = a, b1 = b;
if (b <= a) {
a1 = b;
b1 = a;
}
size_t index = (2 * cols - a1 + 2) * (a1 - 1) / 2 + b1 - a1;
// cout << "* GetabIndx " << a1 << "," << b1 << "," << cols << ":" << index << endl;
// FIXME: should add a contract for index range
return index;
// return ( b < a ? ((2 * cols - b + 2) * (b - 1) / 2 + a - b ): ((2 * cols - a + 2) * (a - 1) / 2 + b - a) );
}
// From an existing n by nd (centered) G matrix, compute the d+1 by
// d*(d-1)/2*(d+1) Q matrix where inside i'th d+1 by d+1 matrix, each
// element is tr(KiKlKjKm)-r*tr(KmKiKl)-r*tr(KlKjKm)+r^2*tr(KlKm),
// where r=n/(n-1)
void compKtoV(const gsl_matrix *G, gsl_matrix *V) {
size_t n_vc = G->size2 / G->size1, ni_test = G->size1;
gsl_matrix *KiKj =
gsl_matrix_alloc(ni_test, (n_vc * (n_vc + 1)) / 2 * ni_test);
gsl_vector *trKiKj = gsl_vector_alloc(n_vc * (n_vc + 1) / 2);
gsl_vector *trKi = gsl_vector_alloc(n_vc);
assert(ni_test != 1);
double d, tr, r = (double)ni_test / (double)(ni_test - 1);
size_t t, t_il, t_jm, t_lm, t_im, t_jl, t_ij;
// Compute KiKj for all pairs of i and j (not including the identity
// matrix).
t = 0;
for (size_t i = 0; i < n_vc; i++) {
gsl_matrix_const_view Ki =
gsl_matrix_const_submatrix(G, 0, i * ni_test, ni_test, ni_test);
for (size_t j = i; j < n_vc; j++) {
gsl_matrix_const_view Kj =
gsl_matrix_const_submatrix(G, 0, j * ni_test, ni_test, ni_test);
gsl_matrix_view KiKj_sub =
gsl_matrix_submatrix(KiKj, 0, t * ni_test, ni_test, ni_test);
fast_dgemm("N", "N", 1.0, &Ki.matrix, &Kj.matrix, 0.0,
&KiKj_sub.matrix);
t++;
}
}
// Compute trKi, trKiKj.
t = 0;
for (size_t i = 0; i < n_vc; i++) {
for (size_t j = i; j < n_vc; j++) {
tr = 0;
for (size_t k = 0; k < ni_test; k++) {
tr += gsl_matrix_get(KiKj, k, t * ni_test + k);
}
gsl_vector_set(trKiKj, t, tr);
t++;
}
tr = 0;
for (size_t k = 0; k < ni_test; k++) {
tr += gsl_matrix_get(G, k, i * ni_test + k);
}
gsl_vector_set(trKi, i, tr);
}
// Compute V.
for (size_t i = 0; i < n_vc; i++) {
for (size_t j = i; j < n_vc; j++) {
t_ij = GetabIndex(i + 1, j + 1, n_vc - 2);
for (size_t l = 0; l < n_vc + 1; l++) {
for (size_t m = 0; m < n_vc + 1; m++) {
if (l != n_vc && m != n_vc) {
t_il = GetabIndex(i + 1, l + 1, n_vc - 2);
t_jm = GetabIndex(j + 1, m + 1, n_vc - 2);
t_lm = GetabIndex(l + 1, m + 1, n_vc - 2);
tr = 0;
for (size_t k = 0; k < ni_test; k++) {
gsl_vector_const_view KiKl_row =
gsl_matrix_const_subrow(KiKj, k, t_il * ni_test, ni_test);
gsl_vector_const_view KiKl_col =
gsl_matrix_const_column(KiKj, t_il * ni_test + k);
gsl_vector_const_view KjKm_row =
gsl_matrix_const_subrow(KiKj, k, t_jm * ni_test, ni_test);
gsl_vector_const_view KjKm_col =
gsl_matrix_const_column(KiKj, t_jm * ni_test + k);
gsl_vector_const_view Kl_row =
gsl_matrix_const_subrow(G, k, l * ni_test, ni_test);
gsl_vector_const_view Km_row =
gsl_matrix_const_subrow(G, k, m * ni_test, ni_test);
if (i <= l && j <= m) {
gsl_blas_ddot(&KiKl_row.vector, &KjKm_col.vector, &d);
tr += d;
gsl_blas_ddot(&Km_row.vector, &KiKl_col.vector, &d);
tr -= r * d;
gsl_blas_ddot(&Kl_row.vector, &KjKm_col.vector, &d);
tr -= r * d;
} else if (i <= l && j > m) {
gsl_blas_ddot(&KiKl_row.vector, &KjKm_row.vector, &d);
tr += d;
gsl_blas_ddot(&Km_row.vector, &KiKl_col.vector, &d);
tr -= r * d;
gsl_blas_ddot(&Kl_row.vector, &KjKm_row.vector, &d);
tr -= r * d;
} else if (i > l && j <= m) {
gsl_blas_ddot(&KiKl_col.vector, &KjKm_col.vector, &d);
tr += d;
gsl_blas_ddot(&Km_row.vector, &KiKl_row.vector, &d);
tr -= r * d;
gsl_blas_ddot(&Kl_row.vector, &KjKm_col.vector, &d);
tr -= r * d;
} else {
gsl_blas_ddot(&KiKl_col.vector, &KjKm_row.vector, &d);
tr += d;
gsl_blas_ddot(&Km_row.vector, &KiKl_row.vector, &d);
tr -= r * d;
gsl_blas_ddot(&Kl_row.vector, &KjKm_row.vector, &d);
tr -= r * d;
}
}
tr += r * r * gsl_vector_get(trKiKj, t_lm);
} else if (l != n_vc && m == n_vc) {
t_il = GetabIndex(i + 1, l + 1, n_vc - 2);
t_jl = GetabIndex(j + 1, l + 1, n_vc - 2);
tr = 0;
for (size_t k = 0; k < ni_test; k++) {
gsl_vector_const_view KiKl_row =
gsl_matrix_const_subrow(KiKj, k, t_il * ni_test, ni_test);
gsl_vector_const_view KiKl_col =
gsl_matrix_const_column(KiKj, t_il * ni_test + k);
gsl_vector_const_view Kj_row =
gsl_matrix_const_subrow(G, k, j * ni_test, ni_test);
if (i <= l) {
gsl_blas_ddot(&KiKl_row.vector, &Kj_row.vector, &d);
tr += d;
} else {
gsl_blas_ddot(&KiKl_col.vector, &Kj_row.vector, &d);
tr += d;
}
}
tr += -r * gsl_vector_get(trKiKj, t_il) -
r * gsl_vector_get(trKiKj, t_jl) +
r * r * gsl_vector_get(trKi, l);
} else if (l == n_vc && m != n_vc) {
t_jm = GetabIndex(j + 1, m + 1, n_vc - 2);
t_im = GetabIndex(i + 1, m + 1, n_vc - 2);
tr = 0;
for (size_t k = 0; k < ni_test; k++) {
gsl_vector_const_view KjKm_row =
gsl_matrix_const_subrow(KiKj, k, t_jm * ni_test, ni_test);
gsl_vector_const_view KjKm_col =
gsl_matrix_const_column(KiKj, t_jm * ni_test + k);
gsl_vector_const_view Ki_row =
gsl_matrix_const_subrow(G, k, i * ni_test, ni_test);
if (j <= m) {
gsl_blas_ddot(&KjKm_row.vector, &Ki_row.vector, &d);
tr += d;
} else {
gsl_blas_ddot(&KjKm_col.vector, &Ki_row.vector, &d);
tr += d;
}
}
tr += -r * gsl_vector_get(trKiKj, t_im) -
r * gsl_vector_get(trKiKj, t_jm) +
r * r * gsl_vector_get(trKi, m);
} else {
tr = gsl_vector_get(trKiKj, t_ij) - r * gsl_vector_get(trKi, i) -
r * gsl_vector_get(trKi, j) + r * r * (double)(ni_test - 1);
}
gsl_matrix_set(V, l, t_ij * (n_vc + 1) + m, tr);
}
}
}
}
assert(ni_test != 0);
gsl_matrix_scale(V, 1.0 / pow((double)ni_test, 2));
gsl_matrix_free(KiKj);
gsl_vector_free(trKiKj);
gsl_vector_free(trKi);
return;
}
// Perform Jacknife sampling for variance of S.
void JackknifeAKtoS(const gsl_matrix *W, const gsl_matrix *A,
const gsl_matrix *K, gsl_matrix *S, gsl_matrix *Svar) {
size_t n_vc = Svar->size1, ni_test = A->size1, n_cvt = W->size2;
vector<vector<vector<double>>> trAK, sumAK;
vector<vector<double>> sumA, sumK, trA, trK, sA, sK;
vector<double> vec_tmp;
double di, dj, d, m, v;
// Initialize and set all elements to zero.
for (size_t i = 0; i < ni_test; i++) {
vec_tmp.push_back(0);
}
for (size_t i = 0; i < n_vc; i++) {
sumA.push_back(vec_tmp);
sumK.push_back(vec_tmp);
trA.push_back(vec_tmp);
trK.push_back(vec_tmp);
sA.push_back(vec_tmp);
sK.push_back(vec_tmp);
}
for (size_t i = 0; i < n_vc; i++) {
trAK.push_back(sumK);
sumAK.push_back(sumK);
}
// Run jackknife.
for (size_t i = 0; i < n_vc; i++) {
for (size_t l = 0; l < ni_test; l++) {
for (size_t k = 0; k < ni_test; k++) {
di = gsl_matrix_get(A, l, k + ni_test * i);
dj = gsl_matrix_get(K, l, k + ni_test * i);
for (size_t t = 0; t < ni_test; t++) {
if (t == l || t == k) {
continue;
}
sumA[i][t] += di;
sumK[i][t] += dj;
if (l == k) {
trA[i][t] += di;
trK[i][t] += dj;
}
}
sA[i][l] += di;
sK[i][l] += dj;
}
}
for (size_t t = 0; t < ni_test; t++) {
assert(ni_test != 1);
sumA[i][t] /= (double)(ni_test - 1);
sumK[i][t] /= (double)(ni_test - 1);
}
}
for (size_t i = 0; i < n_vc; i++) {
for (size_t j = 0; j < n_vc; j++) {
for (size_t l = 0; l < ni_test; l++) {
for (size_t k = 0; k < ni_test; k++) {
di = gsl_matrix_get(A, l, k + ni_test * i);
dj = gsl_matrix_get(K, l, k + ni_test * j);
d = di * dj;
for (size_t t = 0; t < ni_test; t++) {
if (t == l || t == k) {
continue;
}
trAK[i][j][t] += d;
}
}
for (size_t t = 0; t < ni_test; t++) {
if (t == l) {
continue;
}
di = gsl_matrix_get(A, l, t + ni_test * i);
dj = gsl_matrix_get(K, l, t + ni_test * j);
sumAK[i][j][t] += (sA[i][l] - di) * (sK[j][l] - dj);
}
}
for (size_t t = 0; t < ni_test; t++) {
assert(ni_test != 1);
sumAK[i][j][t] /= (double)(ni_test - 1);
}
m = 0;
v = 0;
for (size_t t = 0; t < ni_test; t++) {
d = trAK[i][j][t] - 2 * sumAK[i][j][t] + sumA[i][t] * sumK[j][t];
if ((trA[i][t] - sumA[i][t]) == 0 || (trK[j][t] - sumK[j][t]) == 0) {
d = 0;
} else {
d /= (trA[i][t] - sumA[i][t]) * (trK[j][t] - sumK[j][t]);
d -= 1 / (double)(ni_test - n_cvt - 1);
}
m += d;
v += d * d;
}
m /= (double)ni_test;
v /= (double)ni_test;
v -= m * m;
v *= (double)(ni_test - 1);
gsl_matrix_set(Svar, i, j, v);
if (n_cvt == 1) {
d = gsl_matrix_get(S, i, j);
d = (double)ni_test * d - (double)(ni_test - 1) * m;
gsl_matrix_set(S, i, j, d);
}
}
}
return;
}
// Compute the d by d S matrix with its d by d variance matrix of
// Svar, and the d+1 by d(d+1) matrix of Q for V(q).
void PARAM::CalcS(const map<string, double> &mapRS2wA,
const map<string, double> &mapRS2wK, const gsl_matrix *W,
gsl_matrix *A, gsl_matrix *K, gsl_matrix *S, gsl_matrix *Svar,
gsl_vector *ns) {
string file_str;
gsl_matrix_set_zero(S);
gsl_matrix_set_zero(Svar);
gsl_vector_set_zero(ns);
// Compute the kinship matrix G for multiple categories; these
// matrices are not centered, for convienence of Jacknife sampling.
if (!file_bfile.empty()) {
file_str = file_bfile + ".bed";
if (mapRS2wA.size() == 0) {
if (PlinkKin(file_str, d_pace, indicator_idv, indicator_snp, mapRS2wK,
mapRS2cat, snpInfo, W, K, ns) == false) {
error = true;
}
} else {
if (PlinkKin(file_str, d_pace, indicator_idv, indicator_snp, mapRS2wA,
mapRS2cat, snpInfo, W, A, ns) == false) {
error = true;
}
}
} else if (!file_geno.empty()) {
file_str = file_geno;
if (mapRS2wA.size() == 0) {
if (BimbamKinUncentered(file_str, setKSnps, d_pace, indicator_idv,
indicator_snp, mapRS2wK, mapRS2cat, snpInfo, W, K,
ns) == false) {
error = true;
}
} else {
if (BimbamKinUncentered(file_str, setKSnps, d_pace, indicator_idv,
indicator_snp, mapRS2wA, mapRS2cat, snpInfo, W, A,
ns) == false) {
error = true;
}
}
} else if (!file_mbfile.empty()) {
if (mapRS2wA.size() == 0) {
if (MFILEKin(1, file_mbfile, setKSnps, d_pace, indicator_idv,
mindicator_snp, mapRS2wK, mapRS2cat, msnpInfo, W, K,
ns) == false) {
error = true;
}
} else {
if (MFILEKin(1, file_mbfile, setKSnps, d_pace, indicator_idv,
mindicator_snp, mapRS2wA, mapRS2cat, msnpInfo, W, A,
ns) == false) {
error = true;
}
}
} else if (!file_mgeno.empty()) {
if (mapRS2wA.size() == 0) {
if (MFILEKin(0, file_mgeno, setKSnps, d_pace, indicator_idv,
mindicator_snp, mapRS2wK, mapRS2cat, msnpInfo, W, K,
ns) == false) {
error = true;
}
} else {
if (MFILEKin(0, file_mgeno, setKSnps, d_pace, indicator_idv,
mindicator_snp, mapRS2wA, mapRS2cat, msnpInfo, W, A,
ns) == false) {
error = true;
}
}
}
if (mapRS2wA.size() == 0) {
gsl_matrix_memcpy(A, K);
}
// Center and scale every kinship matrix inside G.
for (size_t i = 0; i < n_vc; i++) {
gsl_matrix_view Ksub =
gsl_matrix_submatrix(K, 0, i * ni_test, ni_test, ni_test);
CenterMatrix(&Ksub.matrix);
// Scale the matrix G such that the mean diagonal = 1.
ScaleMatrix(&Ksub.matrix);
gsl_matrix_view Asub =
gsl_matrix_submatrix(A, 0, i * ni_test, ni_test, ni_test);
CenterMatrix(&Asub.matrix);
ScaleMatrix(&Asub.matrix);
}
// Cased on G, compute S.
compAKtoS(A, K, W->size2, S);
// Compute Svar and update S with Jacknife.
JackknifeAKtoS(W, A, K, S, Svar);
return;
}
void PARAM::WriteVector(const gsl_vector *q, const gsl_vector *s,
const size_t n_total, const string suffix) {
string file_str;
file_str = path_out + "/" + file_out;
file_str += ".";
file_str += suffix;
file_str += ".txt";
ofstream outfile(file_str.c_str(), ofstream::out);
if (!outfile) {
cout << "error writing file: " << file_str.c_str() << endl;
return;
}
outfile.precision(10);
for (size_t i = 0; i < q->size; ++i) {
outfile << gsl_vector_get(q, i) << endl;
}
for (size_t i = 0; i < s->size; ++i) {
outfile << gsl_vector_get(s, i) << endl;
}
outfile << n_total << endl;
outfile.close();
outfile.clear();
return;
}
void PARAM::WriteVar(const string suffix) {
string file_str, rs;
file_str = path_out + "/" + file_out;
file_str += ".";
file_str += suffix;
file_str += ".txt.gz";
ogzstream outfile(file_str.c_str(), ogzstream::out);
if (!outfile) {
cout << "error writing file: " << file_str.c_str() << endl;
return;
}
outfile.precision(10);
if (mindicator_snp.size() != 0) {
for (size_t t = 0; t < mindicator_snp.size(); t++) {
indicator_snp = mindicator_snp[t];
for (size_t i = 0; i < indicator_snp.size(); i++) {
if (indicator_snp[i] == 0) {
continue;
}
rs = snpInfo[i].rs_number;
outfile << rs << endl;
}
}
} else {
for (size_t i = 0; i < indicator_snp.size(); i++) {
if (indicator_snp[i] == 0) {
continue;
}
rs = snpInfo[i].rs_number;
outfile << rs << endl;
}
}
outfile.close();
outfile.clear();
return;
}
void PARAM::WriteMatrix(const gsl_matrix *matrix_U, const string suffix) {
string file_str;
file_str = path_out + "/" + file_out;
file_str += ".";
file_str += suffix;
file_str += ".txt";
ofstream outfile(file_str.c_str(), ofstream::out);
if (!outfile) {
cout << "error writing file: " << file_str.c_str() << endl;
return;
}
outfile.precision(10);
for (size_t i = 0; i < matrix_U->size1; ++i) {
for (size_t j = 0; j < matrix_U->size2; ++j) {
outfile << tab(j) << gsl_matrix_get(matrix_U, i, j);
}
outfile << endl;
}
outfile.close();
outfile.clear();
return;
}
void PARAM::WriteVector(const gsl_vector *vector_D, const string suffix) {
string file_str;
file_str = path_out + "/" + file_out;
file_str += ".";
file_str += suffix;
file_str += ".txt";
ofstream outfile(file_str.c_str(), ofstream::out);
if (!outfile) {
cout << "error writing file: " << file_str.c_str() << endl;
return;
}
outfile.precision(10);
for (size_t i = 0; i < vector_D->size; ++i) {
outfile << gsl_vector_get(vector_D, i) << endl;
}
outfile.close();
outfile.clear();
return;
}
void PARAM::CheckCvt() {
if (indicator_cvt.size() == 0) {
return;
}
size_t ci_test = 0;
gsl_matrix *W = gsl_matrix_alloc(ni_test, n_cvt);
for (vector<int>::size_type i = 0; i < indicator_idv.size(); ++i) {
if (indicator_idv[i] == 0 || indicator_cvt[i] == 0) {
continue;
}
for (size_t j = 0; j < n_cvt; ++j) {
gsl_matrix_set(W, ci_test, j, (cvt)[i][j]);
}
ci_test++;
}
size_t flag_ipt = 0;
double v_min, v_max;
set<size_t> set_remove;
// Check if any columns is an intercept.
for (size_t i = 0; i < W->size2; i++) {
gsl_vector_view w_col = gsl_matrix_column(W, i);
gsl_vector_minmax(&w_col.vector, &v_min, &v_max);
if (v_min == v_max) {
flag_ipt = 1;
set_remove.insert(i);
}
}
// Add an intercept term if needed.
if (n_cvt == set_remove.size()) {
indicator_cvt.clear();
n_cvt = 1;
} else if (flag_ipt == 0) {
info_msg("no intercept term is found in the cvt file: a column of 1s is added");
for (vector<int>::size_type i = 0; i < indicator_idv.size(); ++i) {
if (indicator_idv[i] == 0 || indicator_cvt[i] == 0) {
continue;
}
cvt[i].push_back(1.0);
}
n_cvt++;
} else {
}
gsl_matrix_free(W);
return;
}
// Post-process phenotypes and covariates.
void PARAM::ProcessCvtPhen() {
// Convert indicator_pheno to indicator_idv.
int k = 1;
indicator_idv.clear();
for (size_t i = 0; i < indicator_pheno.size(); i++) {
k = 1;
for (size_t j = 0; j < indicator_pheno[i].size(); j++) {
if (indicator_pheno[i][j] == 0) {
k = 0;
}
}
indicator_idv.push_back(k);
}
// Remove individuals with missing covariates.
if ((indicator_cvt).size() != 0) {
for (vector<int>::size_type i = 0; i < (indicator_idv).size(); ++i) {
indicator_idv[i] *= indicator_cvt[i];
}
}
// Remove individuals with missing gxe variables.
if ((indicator_gxe).size() != 0) {
for (vector<int>::size_type i = 0; i < (indicator_idv).size(); ++i) {
indicator_idv[i] *= indicator_gxe[i];
}
}
// Remove individuals with missing residual weights.
if ((indicator_weight).size() != 0) {
for (vector<int>::size_type i = 0; i < (indicator_idv).size(); ++i) {
indicator_idv[i] *= indicator_weight[i];
}
}
// Obtain ni_test.
ni_test = 0;
for (vector<int>::size_type i = 0; i < (indicator_idv).size(); ++i) {
if (indicator_idv[i] == 0) {
continue;
}
ni_test++;
}
// If subsample number is set, perform a random sub-sampling
// to determine the subsampled ids.
if (ni_subsample != 0) {
if (ni_test < ni_subsample) {
cout << "error! number of subsamples is less than number of"
<< "analyzed individuals. " << endl;
} else {
// From ni_test, sub-sample ni_subsample.
vector<size_t> a, b;
for (size_t i = 0; i < ni_subsample; i++) {
a.push_back(0);
}
for (size_t i = 0; i < ni_test; i++) {
b.push_back(i);
}
gsl_ran_choose(gsl_r, static_cast<void *>(&a[0]), ni_subsample,
static_cast<void *>(&b[0]), ni_test, sizeof(size_t));
// Re-set indicator_idv and ni_test.
int j = 0;
for (vector<int>::size_type i = 0; i < (indicator_idv).size(); ++i) {
if (indicator_idv[i] == 0) {
continue;
}
if (find(a.begin(), a.end(), j) == a.end()) {
indicator_idv[i] = 0;
}
j++;
}
ni_test = ni_subsample;
}
}
// Check ni_test.
if (a_mode != M_BSLMM5)
enforce_msg(ni_test,"number of analyzed individuals equals 0.");
if (ni_test == 0 && a_mode != M_BSLMM5) {
error = true;
cout << "error! number of analyzed individuals equals 0. " << endl;
}
// Check covariates to see if they are correlated with each
// other, and to see if the intercept term is included.
// After getting ni_test.
// Add or remove covariates.
if (indicator_cvt.size() != 0) {
CheckCvt();
} else {
vector<double> cvt_row;
cvt_row.push_back(1);
for (vector<int>::size_type i = 0; i < (indicator_idv).size(); ++i) {
indicator_cvt.push_back(1);
cvt.push_back(cvt_row);
}
}
return;
}
void PARAM::CopyCvt(gsl_matrix *W) {
size_t ci_test = 0;
for (vector<int>::size_type i = 0; i < indicator_idv.size(); ++i) {
if (indicator_idv[i] == 0 || indicator_cvt[i] == 0) {
continue;
}
for (size_t j = 0; j < n_cvt; ++j) {
gsl_matrix_set(W, ci_test, j, (cvt)[i][j]);
}
ci_test++;
}
return;
}
void PARAM::CopyGxe(gsl_vector *env) {
size_t ci_test = 0;
for (vector<int>::size_type i = 0; i < indicator_idv.size(); ++i) {
if (indicator_idv[i] == 0 || indicator_gxe[i] == 0) {
continue;
}
gsl_vector_set(env, ci_test, gxe[i]);
ci_test++;
}
return;
}
void PARAM::CopyWeight(gsl_vector *w) {
size_t ci_test = 0;
for (vector<int>::size_type i = 0; i < indicator_idv.size(); ++i) {
if (indicator_idv[i] == 0 || indicator_weight[i] == 0) {
continue;
}
gsl_vector_set(w, ci_test, weight[i]);
ci_test++;
}
return;
}
// If flag=0, then use indicator_idv to load W and Y;
// else, use indicator_cvt to load them.
void PARAM::CopyCvtPhen(gsl_matrix *W, gsl_vector *y, size_t flag) {
size_t ci_test = 0;
for (vector<int>::size_type i = 0; i < indicator_idv.size(); ++i) {
if (flag == 0) {
if (indicator_idv[i] == 0) {
continue;
}
} else {
if (indicator_cvt[i] == 0) {
continue;
}
}
gsl_vector_set(y, ci_test, (pheno)[i][0]);
for (size_t j = 0; j < n_cvt; ++j) {
gsl_matrix_set(W, ci_test, j, (cvt)[i][j]);
}
ci_test++;
}
return;
}
// If flag=0, then use indicator_idv to load W and Y;
// else, use indicator_cvt to load them.
void PARAM::CopyCvtPhen(gsl_matrix *W, gsl_matrix *Y, size_t flag) {
size_t ci_test = 0;
for (vector<int>::size_type i = 0; i < indicator_idv.size(); ++i) {
if (flag == 0) {
if (indicator_idv[i] == 0) {
continue;
}
} else {
if (indicator_cvt[i] == 0) {
continue;
}
}
for (size_t j = 0; j < n_ph; ++j) {
gsl_matrix_set(Y, ci_test, j, (pheno)[i][j]);
}
for (size_t j = 0; j < n_cvt; ++j) {
gsl_matrix_set(W, ci_test, j, (cvt)[i][j]);
}
ci_test++;
}
return;
}
void PARAM::CopyRead(gsl_vector *log_N) {
size_t ci_test = 0;
for (vector<int>::size_type i = 0; i < indicator_idv.size(); ++i) {
if (indicator_idv[i] == 0) {
continue;
}
gsl_vector_set(log_N, ci_test, log(vec_read[i]));
ci_test++;
}
return;
}
void PARAM::ObtainWeight(const set<string> &setSnps_beta,
map<string, double> &mapRS2wK) {
mapRS2wK.clear();
vector<double> wsum, wcount;
for (size_t i = 0; i < n_vc; i++) {
wsum.push_back(0.0);
wcount.push_back(0.0);
}
string rs;
if (msnpInfo.size() == 0) {
for (size_t i = 0; i < snpInfo.size(); i++) {
if (indicator_snp[i] == 0) {
continue;
}
rs = snpInfo[i].rs_number;
if ((setSnps_beta.size() == 0 || setSnps_beta.count(rs) != 0) &&
(mapRS2wsnp.size() == 0 || mapRS2wsnp.count(rs) != 0) &&
(mapRS2wcat.size() == 0 || mapRS2wcat.count(rs) != 0) &&
(mapRS2cat.size() == 0 || mapRS2cat.count(rs) != 0)) {
if (mapRS2wsnp.size() != 0) {
mapRS2wK[rs] = mapRS2wsnp[rs];
if (mapRS2cat.size() == 0) {
wsum[0] += mapRS2wsnp[rs];
} else {
wsum[mapRS2cat[rs]] += mapRS2wsnp[rs];
}
wcount[0]++;
} else {
mapRS2wK[rs] = 1;
}
}
}
} else {
for (size_t t = 0; t < msnpInfo.size(); t++) {
snpInfo = msnpInfo[t];
indicator_snp = mindicator_snp[t];
for (size_t i = 0; i < snpInfo.size(); i++) {
if (indicator_snp[i] == 0) {
continue;
}
rs = snpInfo[i].rs_number;
if ((setSnps_beta.size() == 0 || setSnps_beta.count(rs) != 0) &&
(mapRS2wsnp.size() == 0 || mapRS2wsnp.count(rs) != 0) &&
(mapRS2wcat.size() == 0 || mapRS2wcat.count(rs) != 0) &&
(mapRS2cat.size() == 0 || mapRS2cat.count(rs) != 0)) {
if (mapRS2wsnp.size() != 0) {
mapRS2wK[rs] = mapRS2wsnp[rs];
if (mapRS2cat.size() == 0) {
wsum[0] += mapRS2wsnp[rs];
} else {
wsum[mapRS2cat[rs]] += mapRS2wsnp[rs];
}
wcount[0]++;
} else {
mapRS2wK[rs] = 1;
}
}
}
}
}
if (mapRS2wsnp.size() != 0) {
for (size_t i = 0; i < n_vc; i++) {
wsum[i] /= wcount[i];
}
for (map<string, double>::iterator it = mapRS2wK.begin();
it != mapRS2wK.end(); ++it) {
if (mapRS2cat.size() == 0) {
it->second /= wsum[0];
} else {
it->second /= wsum[mapRS2cat[it->first]];
}
}
}
return;
}
// If pve_flag=0 then do not change pve; pve_flag==1, then change pve
// to 0 if pve < 0 and pve to 1 if pve > 1.
void PARAM::UpdateWeight(const size_t pve_flag,
const map<string, double> &mapRS2wK,
const size_t ni_test, const gsl_vector *ns,
map<string, double> &mapRS2wA) {
double d;
vector<double> wsum, wcount;
for (size_t i = 0; i < n_vc; i++) {
wsum.push_back(0.0);
wcount.push_back(0.0);
}
for (map<string, double>::const_iterator it = mapRS2wK.begin();
it != mapRS2wK.end(); ++it) {
d = 1;
for (size_t i = 0; i < n_vc; i++) {
if (v_pve[i] >= 1 && pve_flag == 1) {
d += (double)ni_test / gsl_vector_get(ns, i) * mapRS2wcat[it->first][i];
} else if (v_pve[i] <= 0 && pve_flag == 1) {
d += 0;
} else {
d += (double)ni_test / gsl_vector_get(ns, i) *
mapRS2wcat[it->first][i] * v_pve[i];
}
}
mapRS2wA[it->first] = 1 / (d * d);
if (mapRS2cat.size() == 0) {
wsum[0] += mapRS2wA[it->first];
wcount[0]++;
} else {
wsum[mapRS2cat[it->first]] += mapRS2wA[it->first];
wcount[mapRS2cat[it->first]]++;
}
}
for (size_t i = 0; i < n_vc; i++) {
wsum[i] /= wcount[i];
}
for (map<string, double>::iterator it = mapRS2wA.begin();
it != mapRS2wA.end(); ++it) {
if (mapRS2cat.size() == 0) {
it->second /= wsum[0];
} else {
it->second /= wsum[mapRS2cat[it->first]];
}
}
return;
}
// This function updates indicator_snp, and save z-scores and other
// values into vectors.
void PARAM::UpdateSNPnZ(const map<string, double> &mapRS2wA,
const map<string, string> &mapRS2A1,
const map<string, double> &mapRS2z, gsl_vector *w,
gsl_vector *z, vector<size_t> &vec_cat) {
gsl_vector_set_zero(w);
gsl_vector_set_zero(z);
vec_cat.clear();
string rs, a1;
size_t c = 0;
if (msnpInfo.size() == 0) {
for (size_t i = 0; i < snpInfo.size(); i++) {
if (indicator_snp[i] == 0) {
continue;
}
rs = snpInfo[i].rs_number;
a1 = snpInfo[i].a_minor;
if (mapRS2wA.count(rs) != 0) {
if (a1 == mapRS2A1.at(rs)) {
gsl_vector_set(z, c, mapRS2z.at(rs));
} else {
gsl_vector_set(z, c, -1 * mapRS2z.at(rs));
}
vec_cat.push_back(mapRS2cat.at(rs));
gsl_vector_set(w, c, mapRS2wA.at(rs));
c++;
} else {
indicator_snp[i] = 0;
}
}
} else {
for (size_t t = 0; t < msnpInfo.size(); t++) {
snpInfo = msnpInfo[t];
for (size_t i = 0; i < snpInfo.size(); i++) {
if (mindicator_snp[t][i] == 0) {
continue;
}
rs = snpInfo[i].rs_number;
a1 = snpInfo[i].a_minor;
if (mapRS2wA.count(rs) != 0) {
if (a1 == mapRS2A1.at(rs)) {
gsl_vector_set(z, c, mapRS2z.at(rs));
} else {
gsl_vector_set(z, c, -1 * mapRS2z.at(rs));
}
vec_cat.push_back(mapRS2cat.at(rs));
gsl_vector_set(w, c, mapRS2wA.at(rs));
c++;
} else {
mindicator_snp[t][i] = 0;
}
}
}
}
return;
}
// This function updates indicator_snp, and save z-scores and other
// values into vectors.
void PARAM::UpdateSNP(const map<string, double> &mapRS2wA) {
string rs;
if (msnpInfo.size() == 0) {
for (size_t i = 0; i < snpInfo.size(); i++) {
if (indicator_snp[i] == 0) {
continue;
}
rs = snpInfo[i].rs_number;
if (mapRS2wA.count(rs) == 0) {
indicator_snp[i] = 0;
}
}
} else {
for (size_t t = 0; t < msnpInfo.size(); t++) {
snpInfo = msnpInfo[t];
for (size_t i = 0; i < mindicator_snp[t].size(); i++) {
if (mindicator_snp[t][i] == 0) {
continue;
}
rs = snpInfo[i].rs_number;
if (mapRS2wA.count(rs) == 0) {
mindicator_snp[t][i] = 0;
}
}
}
}
return;
}
