Commit 67861457 authored by zhanxw's avatar zhanxw
Browse files

add boltCov model

parent 5a7f17eb
......@@ -740,6 +740,28 @@ class BoltLMM::BoltLMMImpl {
double GetEffect() { return effect_; }
double GetPvalue() { return pvalue_; }
void GetCovXX(const std::vector<double>& g1, const std::vector<double>& g2,
double* out) {
assert(g1.size() == g2.size());
assert(g1.size() == N_);
Eigen::MatrixXf g(g1.size() + C_, 2); //
// copy in data
for (size_t i = 0; i != N_; ++i) {
g(i, 0) = g1[i];
g(i, 1) = g2[i];
}
// project
pl.projectCovariate(&g);
// calculate
(*out) = (double)projDot(g.col(0), g.col(1))(0);
// scale
(*out) *= xVx_xx_ratio_;
}
private:
int EstimateHeritability() {
TIMER(__PRETTY_FUNCTION__);
......@@ -820,8 +842,14 @@ class BoltLMM::BoltLMMImpl {
for (i = 2; i < 7; ++i) {
logDelta[i] = (logDelta[i - 2] * f[i - 1] - logDelta[i - 1] * f[i - 2]) /
(f[i - 1] - f[i - 2]);
if (logDelta[i] > 10) {
logDelta[i] = 10;
if (!std::isfinite(logDelta[i])) {
--i;
break;
}
// changed the upper threshold from 10 to 5
// as it does not seem to be stable, and causes f(10) = nan
if (logDelta[i] > 5) {
logDelta[i] = 5;
}
if (logDelta[i] < -10) {
logDelta[i] = -10;
......@@ -971,6 +999,13 @@ class BoltLMM::BoltLMMImpl {
for (int i = 0; i < maxIter; ++i) {
computeHx(delta, p, &ap);
alpha = rsold.array() / projDot(p, ap).array();
#ifdef DEBUG
fprintf(stderr, "alpha:");
for (int ii = 0; ii != NUM_COL; ++ii) {
fprintf(stderr, "\t%f", alpha(ii));
}
fprintf(stderr, "\n");
#endif
for (int ii = 0; ii != NUM_COL; ++ii) {
if (!std::isfinite(alpha(ii)) || fabs(alpha(ii)) < Tol) {
alpha(ii) = 0.0;
......@@ -983,7 +1018,27 @@ class BoltLMM::BoltLMMImpl {
if ((rsnew.array() < Tol).all()) {
break;
}
if ((rsnew - rsold).array().abs().maxCoeff() < Tol) {
break;
}
ratio = rsnew.array() / rsold.array();
#ifdef DEBUG
fprintf(stderr, "ratio:");
for (int ii = 0; ii != NUM_COL; ++ii) {
fprintf(stderr, "\t%f", ratio(ii));
}
fprintf(stderr, "\n");
fprintf(stderr, "rsnew:");
for (int ii = 0; ii != NUM_COL; ++ii) {
fprintf(stderr, "\t%f", rsnew(ii));
}
fprintf(stderr, "\n");
fprintf(stderr, "rsold:");
for (int ii = 0; ii != NUM_COL; ++ii) {
fprintf(stderr, "\t%f", rsold(ii));
}
fprintf(stderr, "\n");
#endif
for (int ii = 0; ii != NUM_COL; ++ii) {
if (rsnew(ii) < Tol || !std::isfinite(ratio(ii))) {
ratio(ii) = 0.0;
......@@ -1072,6 +1127,7 @@ class BoltLMM::BoltLMMImpl {
// X: [X; Z'X] [ (N+C) x M ]
//
// X_y: [X' -X'Z] [ M by (MCtrial+1) ]
ret->resize(N_ + C_, y.cols());
ret->setZero();
Eigen::MatrixXf X_y =
Eigen::MatrixXf::Zero(M2_, y.cols()); // X' * y: [ M * (MCtrial + 1) ]
......@@ -1235,6 +1291,19 @@ class BoltLMM::BoltLMMImpl {
prospectiveStat(i) / uncalibratedRetrospectiveStat(i));
}
fprintf(stderr, "infStatCalibration_ = %f\n", infStatCalibration_);
// calculate empirically X'HX / X'X
xVx_xx_ratio_ = x_V_inv_x.sum() / x_x.sum();
if (!std::isfinite(xVx_xx_ratio_)) {
xVx_xx_ratio_ = 1.0;
}
fprintf(stderr, "\ni\tx_v_inv_x\tx_x\tratio\n");
for (int i = 0; i < nSnp; ++i) {
fprintf(stderr, "%d\t%f\t%f\t%f\n", i, x_V_inv_x(i), x_x(i),
x_x(i) == 0 ? 0.0 : x_V_inv_x(i) / x_x(i));
}
fprintf(stderr, "ratio = %f\n", x_V_inv_x.sum() / x_x.sum());
return 0;
}
......@@ -1269,6 +1338,7 @@ class BoltLMM::BoltLMMImpl {
Eigen::MatrixXf H_inv_y_;
double H_inv_y_norm2_;
double infStatCalibration_;
double xVx_xx_ratio_;
Eigen::MatrixXf g_test_;
Eigen::MatrixXf alpha_;
......@@ -1294,6 +1364,10 @@ double BoltLMM::GetU() { return impl_->GetU(); }
double BoltLMM::GetV() { return impl_->GetV(); }
double BoltLMM::GetEffect() { return impl_->GetEffect(); }
double BoltLMM::GetPvalue() { return impl_->GetPvalue(); }
void BoltLMM::GetCovXX(const std::vector<double>& g1,
const std::vector<double>& g2, double* out) {
impl_->GetCovXX(g1, g2, out);
}
//////////////////////////////////////////////////
// BoltLMM::BoltLMMImpl class
......
......@@ -2,6 +2,7 @@
#define _BOLTLMM_H_
#include <string>
#include <vector>
class Matrix;
......@@ -29,6 +30,10 @@ class BoltLMM {
double GetEffect();
double GetPvalue();
// calculate covariances
void GetCovXX(const std::vector<double>& g1, const std::vector<double>& g2,
double* out);
private:
// don't copy
BoltLMM(const BoltLMM&);
......
......@@ -335,7 +335,7 @@ int DataConsolidator::prepareBoltModel(
error(
"%s:%d Fatal error, PLINK file [ %s ] does not include sample [ %s "
"]!\n",
__FILE__, __LINE__, (prefix + "fam").c_str(), sampleName[i].c_str());
__FILE__, __LINE__, (prefix + ".fam").c_str(), sampleName[i].c_str());
}
}
......
......@@ -30,7 +30,7 @@
Logger* logger = NULL;
const char* VERSION = "20170609";
const char* VERSION = "20170613";
void banner(FILE* fp) {
const char* string =
......
......@@ -3837,8 +3837,12 @@ class MetaScoreTest : public ModelFitter {
MetaBase* model;
MetaBase* modelAuto;
MetaBase* modelX;
protected:
// allow inherited-class to change
bool useBolt;
private:
double af; // overall af (unadjust or adjusted by family structure)
bool fitOK;
......@@ -3878,6 +3882,16 @@ class MetaRecessiveTest : public MetaScoreTest {
Matrix geno;
};
class MetaScoreBoltTest : public MetaScoreTest {
public:
MetaScoreBoltTest() : MetaScoreTest() { this->modelName = "MetaScoreBolt"; }
virtual int setParameter(const ModelParser& parser) {
MetaScoreTest::setParameter(parser);
this->useBolt = true;
return 0;
}
};
class MetaCovTest : public ModelFitter {
public:
typedef std::vector<double> Genotype;
......@@ -3896,6 +3910,7 @@ class MetaCovTest : public ModelFitter {
: model(NULL),
modelAuto(NULL),
modelX(NULL),
useBolt(false),
fitOK(false),
useFamilyModel(false),
isHemiRegion(false) {
......@@ -4002,11 +4017,12 @@ class MetaCovTest : public ModelFitter {
for (int i = 0; i < nSample; ++i) {
loci.geno[i] = genotype[i][0];
}
model->transformGenotype(&loci.geno, dc);
model->calculateXZ(loci.geno, &loci.covXZ);
model->calculateZZ(&this->covZZ);
CholeskyInverseMatrix(this->covZZ, &this->covZZInv);
if (!useBolt) {
model->transformGenotype(&loci.geno, dc);
model->calculateXZ(loci.geno, &loci.covXZ);
model->calculateZZ(&this->covZZ);
CholeskyInverseMatrix(this->covZZ, &this->covZZInv);
}
fitOK = true;
return 0;
} // fitWithGivenGenotype
......@@ -4091,9 +4107,14 @@ class MetaCovTest : public ModelFitter {
double covX1X2;
for (int idx = 0; iter != lociQueue.end(); ++iter, ++idx) {
position[idx] = iter->pos.pos;
model->calculateXX(lociQueue.front().geno, iter->geno, &covX1X2);
this->covXX[idx] = computeScaledXX(covX1X2, lociQueue.front().covXZ,
iter->covXZ, this->covZZInv);
if (!useBolt) {
model->calculateXX(lociQueue.front().geno, iter->geno, &covX1X2);
this->covXX[idx] = computeScaledXX(covX1X2, lociQueue.front().covXZ,
iter->covXZ, this->covZZInv);
} else {
model->calculateXX(lociQueue.front().geno, iter->geno,
&this->covXX[idx]);
}
}
result.updateValue("CHROM", lociQueue.front().pos.chrom);
......@@ -4461,9 +4482,43 @@ class MetaCovTest : public ModelFitter {
LogisticRegression logistic;
Vector weight;
}; // end class MetaCovUnrelatedBinary
class MetaCovFamQtlBolt : public MetaCovBase {
public:
MetaCovFamQtlBolt() {
fprintf(stderr, "MetaCovFamQtlBolt model started\n");
}
int FitNullModel(Matrix& genotype, DataConsolidator* dc) {
const std::string& fn = dc->getBoltGenotypeFilePrefix();
// fit null model
bool fitOK = 0 == bolt_.FitNullModel(fn, &dc->getPhenotype());
if (!fitOK) return -1;
needToFitNullModel = false;
return 0;
}
int transformGenotype(Genotype* out, DataConsolidator* dc) { return 0; }
int calculateXX(const Genotype& x1, const Genotype& x2, double* covXX) {
bolt_.GetCovXX(x1, x2, covXX);
return 0;
}
int calculateXZ(const Genotype& x, std::vector<double>* covXZ) { return 0; }
int calculateZZ(Matrix* covZZ) { return 0; }
private:
BoltLMM bolt_;
}; // class MetaCovFamQtlBolt
MetaCovBase* createModel(bool familyModel, bool binaryOutcome) {
MetaCovBase* ret = NULL;
if (this->useBolt) {
if (binaryOutcome) {
fprintf(stderr, "BoltLMM does not support binary outcomes! Exit...\n");
exit(1);
}
ret = new MetaCovFamQtlBolt;
return ret;
}
if (familyModel && !binaryOutcome) {
ret = new MetaCovFamQtl;
}
......@@ -4484,6 +4539,11 @@ class MetaCovTest : public ModelFitter {
MetaCovBase* modelAuto;
MetaCovBase* modelX;
protected:
// allow inherited-class to change
bool useBolt;
private:
std::deque<Loci> queue;
int numVariant;
int nSample;
......@@ -4530,6 +4590,18 @@ class MetaRecessiveCovTest : public MetaCovTest {
Matrix geno;
};
class MetaCovBoltTest : public MetaCovTest {
public:
MetaCovBoltTest(int windowSize) : MetaCovTest(windowSize) {
this->modelName = "MetaCovBolt";
}
virtual int setParameter(const ModelParser& parser) {
MetaCovTest::setParameter(parser);
MetaCovTest::useBolt = true;
return 0;
}
};
#if 0
class MetaSkewTest: public ModelFitter{
private:
......
......@@ -214,6 +214,15 @@ int ModelManager::create(const std::string& modelType,
"under additive model",
toStringWithComma(windowSize).c_str());
model.push_back(new MetaCovTest(windowSize));
} else if (modelName == "bolt") {
model.push_back(new MetaScoreBoltTest());
} else if (modelName == "boltcov") {
parser.assign("windowSize", &windowSize, 1000000);
logger->info(
"Meta analysis uses window size %s to produce covariance statistics "
"under additive model",
toStringWithComma(windowSize).c_str());
model.push_back(new MetaCovBoltTest(windowSize));
}
#if 0
else if (modelName == "skew") {
......
......@@ -45,6 +45,7 @@ class ModelParser {
static const std::string PARAM_DELIM;
private:
// Model name
std::string name;
// Store parsed parameters
// use this->value() and this->hasTag() to access this.
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment