Skip to content

Commit

Permalink
Bug fixes
Browse files Browse the repository at this point in the history
  • Loading branch information
@zhanxw
zhanxw committed Oct 7, 2014
1 parent 39a856f commit 6e1a059
Show file tree
Hide file tree
Showing 5 changed files with 49 additions and 35 deletions.
8 changes: 8 additions & 0 deletions ChangeLog
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,11 @@
2014-10-06 Xiaowei Zhan <[email protected]>

* Fix allele frequency calculation in related samples

2014-07-29 Xiaowei Zhan <[email protected]>

* Change AF = 0 to NA when the sites are all missing

2014-07-23 Xiaowei Zhan <[email protected]>

* Update U, V statistics for unrelated individual such that SE(effect_size) = 1/sqrt(V)
Expand Down
23 changes: 14 additions & 9 deletions regression/FastLMM.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,11 @@
#include "GSLMinimizer.h"
#include "gsl/gsl_cdf.h" // use gsl_cdf_chisq_Q

#if 0
// #define EIGEN_NO_DEBUG
#undef DEBUG
// #define DEBUG

#ifdef DEBUG
#include <fstream>
void dumpToFile(const Eigen::MatrixXf& mat, const char* fn) {
std::ofstream out(fn);
Expand All @@ -15,10 +19,6 @@ void dumpToFile(const Eigen::MatrixXf& mat, const char* fn) {
}
#endif

// #define EIGEN_NO_DEBUG
#undef DEBUG
// #define DEBUG

#define PI 3.1415926535897

static double goalFunction(double x, void* param);
Expand Down Expand Up @@ -94,6 +94,8 @@ class FastLMM::Impl{
this->delta = minimizer.getX();
#ifdef DEBUG
fprintf(stderr, "minimization succeed when delta = %g, sigma2 = %g\n", this->delta, this->sigma2);
dumpToFile(kinshipS.mat, "S.mat");
dumpToFile(kinshipU.mat, "U.mat");
#endif
}
}
Expand Down Expand Up @@ -290,13 +292,14 @@ class FastLMM::Impl{
static Eigen::MatrixXf g;
static Eigen::ArrayXf u1s;
static double denom;
static Eigen::ArrayXf alpha;
// static Eigen::ArrayXf alpha;
static Eigen::RowVectorXf alpha;

if (!initialized) {
Eigen::MatrixXf u1 = U.transpose().rowwise().sum(); // n by 1
u1s = u1.array() / (this->lambda.array() + delta).abs().array();
u1s = u1.array() / (this->lambda.array()).abs().array();
denom = (u1s * u1.array()).sum();
alpha = (u1.transpose() * (this->lambda.array() + delta).abs().inverse().matrix().asDiagonal() * U.transpose()).transpose();
alpha = (u1.transpose() * this->lambda.array().abs().inverse().matrix().asDiagonal() * U.transpose()).row(0);
initialized = true;
// fprintf(stderr, "initialized\n");
}
Expand All @@ -306,11 +309,13 @@ class FastLMM::Impl{
}
G_to_Eigen(Xcol, &g);

double beta = (alpha * g.array()).sum() / denom;
// double beta = (alpha * g.array()).sum() / denom;
double beta = alpha.dot(g.col(0)) / denom;

// here x is represented as 0, 1, 2, so beta(0, 0) is the mean genotype
// multiply by 0.5 to get AF
double af = 0.5 * beta;

return af;
}

Expand Down
6 changes: 3 additions & 3 deletions src/Main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@

Logger* logger = NULL;

const char* VERSION = "20140723";
const char* VERSION = "20141006";

void banner(FILE* fp) {
const char* string =
Expand All @@ -39,7 +39,7 @@ void banner(FILE* fp) {
" ... Bingshan Li, Dajiang Liu ... \n"
" ... Goncalo Abecasis ... \n"
" ... [email protected] ... \n"
" ... July 2014 ... \n"
" ... October 2014 ... \n"
" ... zhanxw.github.io/rvtests/ ... \n"
" .............................................. \n"
" \n"
Expand Down Expand Up @@ -752,7 +752,7 @@ int main(int argc, char** argv){
return -1;
}
}
logger->info("Loaded %zu sample pheontypes.", phenotype.size());
logger->info("Loaded [ %zu ] sample pheontypes.", phenotype.size());

// rearrange phenotypes
std::vector<std::string> vcfSampleNames;
Expand Down
35 changes: 18 additions & 17 deletions src/ModelFitter.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2347,7 +2347,7 @@ class MetaScoreTest: public ModelFitter{
needToFitNullModelForAuto(true),
needToFitNullModelForX(true){
this->modelName = "MetaScore";
};
}
// fitting model
virtual int fit(DataConsolidator* dc) {
Matrix& genotype = dc->getGenotype();
Expand All @@ -2373,7 +2373,7 @@ class MetaScoreTest: public ModelFitter{
if (nSample){
af = 0.5 * (het + 2. * homAlt) / (homRef + het + homAlt);
} else {
af = 0.0;
af = -1.0;
}
// fprintf(stderr, "af = %g\n", af);
// use female info to get HWE
Expand Down Expand Up @@ -2526,12 +2526,11 @@ class MetaScoreTest: public ModelFitter{
}
}
return (fitOK ? 0 : -1);
};
}

// write result header
void writeHeader(FileWriter* fp, const Result& siteInfo) {
if (g_SummaryHeader) {

g_SummaryHeader->outputHeader(fp);
}

Expand All @@ -2555,12 +2554,14 @@ class MetaScoreTest: public ModelFitter{
int informativeAC = het + 2* homAlt;

result.clearValue();
if (!isBinaryOutcome()) {
result.updateValue("AF", af);
} else {
static char afString[128];
snprintf(afString, 128, "%g:%g:%g", af, afFromCase, afFromControl);
result.updateValue("AF", afString);
if (af > 0) {
if (!isBinaryOutcome()) {
result.updateValue("AF", af);
} else {
static char afString[128];
snprintf(afString, 128, "%g:%g:%g", af, afFromCase, afFromControl);
result.updateValue("AF", afString);
}
}

result.updateValue("INFORMATIVE_ALT_AC", informativeAC);
Expand Down Expand Up @@ -2617,7 +2618,7 @@ class MetaScoreTest: public ModelFitter{
}
}
result.writeValueLine(fp);
};
}
private:
double af;
double afFromCase;
Expand Down Expand Up @@ -2700,7 +2701,7 @@ class MetaCovTest: public ModelFitter{
this->fout = NULL;
this->windowSize = windowSize;
this->needToFitNullModelForAuto = true;
this->needToFitNullModelForX = true;
this->needToFitNullModelForX = true;
result.addHeader("CHROM");
result.addHeader("START_POS");
result.addHeader("END_POS");
Expand Down Expand Up @@ -2766,15 +2767,15 @@ class MetaCovTest: public ModelFitter{
needToFitNullModelForAuto = false;
// get weight
metaCovForAuto.GetWeight(&this->weight);
}
}
} else { // hemi region
if (!dc->hasKinshipForX()) {
fitOK = false;
return -1;
}
if (this->needToFitNullModelForX || dc->isPhenotypeUpdated() || dc->isCovariateUpdated()) {
copyCovariateAndIntercept(genotype.rows, covariate, &cov);
fitOK = (0 == metaCovForX.FitNullModel(cov, phenotype, *dc->getKinshipUForX(), *dc->getKinshipSForX()));
fitOK = (0 == metaCovForX.FitNullModel(cov, phenotype, *dc->getKinshipUForX(), *dc->getKinshipSForX()));
if (!fitOK) {
return -1;
}
Expand Down Expand Up @@ -2865,8 +2866,8 @@ class MetaCovTest: public ModelFitter{
if (useFamilyModel) {
if (!isHemiRegion) {
metaCovForAuto.TransformCentered(&loci.geno,
*dc->getKinshipUForAuto(),
*dc->getKinshipSForAuto());
*dc->getKinshipUForAuto(),
*dc->getKinshipSForAuto());
} else {
metaCovForX.TransformCentered(&loci.geno,
*dc->getKinshipUForX(),
Expand Down Expand Up @@ -3091,7 +3092,7 @@ class MetaCovTest: public ModelFitter{
MetaCov metaCovForAuto;
MetaCov metaCovForX;
bool needToFitNullModelForAuto;
bool needToFitNullModelForX;
bool needToFitNullModelForX;
bool useFamilyModel;
Vector weight; // per individual weight
LogisticRegression logistic;
Expand Down
12 changes: 6 additions & 6 deletions third/Makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,11 +19,11 @@ gsl: gsl-1.16.tar.gz
tar zvxf gsl-1.16.tar.gz
-(DIR=`pwd`; cd gsl-1.16; ./configure --prefix="$${DIR}"/gsl; make -j; make install)

eigen: eigen-3.2.1.tar.bz2
eigen: eigen-3.2.2.tar.bz2
-rm -rf eigen-eigen*
tar jvxf eigen-3.2.1.tar.bz2
-mv eigen-eigen-* eigen-3.2.1
ln -s -f eigen-3.2.1 eigen
tar jvxf eigen-3.2.2.tar.bz2
-mv eigen-eigen-* eigen-3.2.2
ln -s -f eigen-3.2.2 eigen

bzip2: bzip2-1.0.6.tar.gz
tar zvxf $<
Expand All @@ -45,9 +45,9 @@ else
DOWNLOAD = curl -L $(1) -O
DOWNLOAD_RENAME = curl -L $(1) -o $(2)
endif
eigen-3.2.1.tar.bz2:
eigen-3.2.2.tar.bz2:
echo "obtain Eigen..."
$(call DOWNLOAD_RENAME,http://bitbucket.org/eigen/eigen/get/3.2.1.tar.bz2,$@)
$(call DOWNLOAD_RENAME,http://bitbucket.org/eigen/eigen/get/3.2.2.tar.bz2,$@)

gsl-1.16.tar.gz:
echo "obtain GSL"
Expand Down

0 comments on commit 6e1a059

Please sign in to comment.