Skip to content

Commit

Permalink
compute llk from data
Browse files Browse the repository at this point in the history
  • Loading branch information
@shajoezhu
shajoezhu committed Nov 30, 2017
1 parent ca6ba1e commit ef0e393
Show file tree
Hide file tree
Showing 4 changed files with 90 additions and 54 deletions.
4 changes: 3 additions & 1 deletion src/dEploid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,9 @@ int main( int argc, char *argv[] ){
return EXIT_SUCCESS;
}

if ( dEploidIO.doLsPainting() ){
if ( dEploidIO.doComputeLLK() ){
dEploidIO.computeLLKfromInitialHap();
} else if ( dEploidIO.doLsPainting() ){
dEploidIO.chromPainting();
} else if ( dEploidIO.doIbdPainting() ){
dEploidIO.paintIBD();
Expand Down
51 changes: 37 additions & 14 deletions src/dEploidIO.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -132,6 +132,7 @@ void DEploidIO::init() {
this->setUseVcf(false);
this->vcfReaderPtr_ = NULL;
this->setDoExportVcf(false);
this->setDoComputeLLK(false);

#ifdef COMPILEDATE
compileTime_ = COMPILEDATE;
Expand Down Expand Up @@ -171,7 +172,7 @@ void DEploidIO::reInit() {


void DEploidIO::finalize(){
if ( this->doIbdPainting() ){
if ( this->doIbdPainting() | this->doComputeLLK() ){
if (!initialPropWasGiven()){
throw InitialPropUngiven("");
}
Expand Down Expand Up @@ -438,6 +439,8 @@ void DEploidIO::parse (){
this->readInitialHaps();
} else if ( *argv_i == "-ibd" ){
this->setUseIBD(true);
} else if ( *argv_i == "-computeLLK" ){
this->setDoComputeLLK( true );
} else if ( *argv_i == "-ibdPainting" ){
this->setDoIbdPainting( true );
} else if ( *argv_i == "-initialP" ){
Expand Down Expand Up @@ -488,7 +491,7 @@ void DEploidIO::checkInput(){
throw FileNameMissing ( "Alt count" );}
if ( this->plafFileName_.size() == 0 ){
throw FileNameMissing ( "PLAF" );}
if ( usePanel() && this->panelFileName_.size() == 0 && !this->doIbdPainting() ){
if ( usePanel() && this->panelFileName_.size() == 0 && !this->doIbdPainting() && !this->doComputeLLK() ){
throw FileNameMissing ( "Reference panel" );}
if ( this->initialPropWasGiven() && ( abs(sumOfVec(initialProp) - 1.0) > 0.00001 )){
throw SumOfPropNotOne ( to_string(sumOfVec(initialProp)) );}
Expand Down Expand Up @@ -608,6 +611,36 @@ void DEploidIO::readInitialHaps(){
}


vector <double> DEploidIO::computeExpectedWsafFromInitialHap(){
// Make this a separate function
// calculate expected wsaf
vector <double> expectedWsaf (this->initialHap.size(), 0.0);
for ( size_t i = 0; i < this->initialHap.size(); i++ ){
assert( kStrain_ == this->initialHap[i].size() );
for ( size_t k = 0; k < this->kStrain_; k++){
expectedWsaf[i] += this->initialHap[i][k] * finalProp[k];
}
assert ( expectedWsaf[i] >= 0 );
//assert ( expectedWsaf[i] <= 1.0 );
}
return expectedWsaf;
}


void DEploidIO::computeLLKfromInitialHap(){
for ( auto const& value: this->initialProp ){
this->finalProp.push_back(value);
}

vector <double> expectedWsaf = computeExpectedWsafFromInitialHap();
if (expectedWsaf.size() != this->refCount_.size()){
throw LociNumberUnequal("Hap length differs from data!");
}
vector <double> llk = calcLLKs ( this->refCount_, this->altCount_, expectedWsaf, 0, expectedWsaf.size(), this->scalingFactor());
this->llkFromInitialHap_ = sumOfVec(llk);
}


void DEploidIO::chromPainting(){
dout << "Painting haplotypes in" << this->initialHapFileName_ <<endl;

Expand Down Expand Up @@ -638,17 +671,7 @@ void DEploidIO::chromPainting(){

//vector < vector <double>> hap = decovolutedStrainsToBeRead.content_;

// Make this a separate function
// calculate expected wsaf
vector <double> expectedWsaf (this->nLoci_, 0.0);
for ( size_t i = 0; i < this->initialHap.size(); i++ ){
assert( kStrain_ == this->initialHap[i].size() );
for ( size_t k = 0; k < this->kStrain_; k++){
expectedWsaf[i] += this->initialHap[i][k] * finalProp[k];
}
assert ( expectedWsaf[i] >= 0 );
//assert ( expectedWsaf[i] <= 1.0 );
}
vector <double> expectedWsaf = computeExpectedWsafFromInitialHap();

MersenneTwister tmpRg(this->randomSeed());

Expand Down Expand Up @@ -690,7 +713,7 @@ void DEploidIO::readPanel(){
if ( this->usePanel() == false ){
return;
}
if ( this->doIbdPainting() ){
if ( this->doIbdPainting() | this->doComputeLLK() ){
return;
}

Expand Down
7 changes: 7 additions & 0 deletions src/dEploidIO.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,6 +68,8 @@ class DEploidIO{
void chromPainting ();
bool doLsPainting() const { return this->doLsPainting_; }
bool doIbdPainting() const { return this->doIbdPainting_; }
bool doComputeLLK() const { return this->doComputeLLK_; }
void computeLLKfromInitialHap();
bool useIBD() const { return this->useIBD_;}
void paintIBD();
double ibdLLK_;
Expand All @@ -81,6 +83,7 @@ class DEploidIO{

private:
void core();
double llkFromInitialHap_;

// Read in input
string plafFileName_;
Expand Down Expand Up @@ -158,6 +161,9 @@ class DEploidIO{
void setDoExportRecombProb( const bool exportRecombProb ){ this->doExportRecombProb_ = exportRecombProb; }
bool doExportRecombProb() const { return this->doExportRecombProb_; }

bool doComputeLLK_;
void setDoComputeLLK( const bool setTo ){ this->doComputeLLK_ = setTo; }

// Parameters
double missCopyProb_;
double averageCentimorganDistance_;// = 15000.0,
Expand Down Expand Up @@ -226,6 +232,7 @@ class DEploidIO{

void set_seed(const size_t seed){ this->randomSeed_ = seed; }
void removeFilesWithSameName();
vector <double> computeExpectedWsafFromInitialHap();


template<class T>
Expand Down
82 changes: 43 additions & 39 deletions src/export/dEploidIOExport.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -123,7 +123,7 @@ void DEploidIO::writeLog ( ostream * writeTo ){
}
}
(*writeTo) << "\n";
if ( (this->doLsPainting() == false) & (this->doIbdPainting() == false) ) {
if ( (this->doLsPainting() == false) & (this->doIbdPainting() == false) & (this->doComputeLLK() == false) ) {
(*writeTo) << "MCMC diagnostic:"<< "\n";
(*writeTo) << setw(19) << " Accept_ratio: " << acceptRatio_ << "\n";
(*writeTo) << setw(19) << " Max_llks: " << maxLLKs_ << "\n";
Expand All @@ -138,48 +138,52 @@ void DEploidIO::writeLog ( ostream * writeTo ){
(*writeTo) << setw(14) << "Start at: " << startingTime_ ;
(*writeTo) << setw(14) << "End at: " << endTime_ ;
(*writeTo) << "\n";
(*writeTo) << "Output saved to:\n";
if ( this->doLsPainting() ){
for ( size_t i = 0; i < kStrain(); i++ ){
(*writeTo) << "Posterior probability of strain " << i << ": "<< strExportSingleFwdProbPrefix << i <<endl;
}
} else if (this->doIbdPainting()){
if (this->ibdProbsIntegrated.size()>1){
(*writeTo) << setw(14) << "IBD probs: " << strIbdExportProbs << "\n\n";
(*writeTo) << " IBD probabilities:\n";
for ( size_t stateI = 0; stateI < this->ibdProbsHeader.size(); stateI++ ){
(*writeTo) << setw(14) << this->ibdProbsHeader[stateI] << ": " << this->ibdProbsIntegrated[stateI] << "\n";
}
}
if ( this->doComputeLLK() ){
(*writeTo) << "Input likelihood: " << llkFromInitialHap_;
(*writeTo) << "\n";
} else {
(*writeTo) << setw(14) << "Likelihood: " << strExportLLK << "\n";
(*writeTo) << setw(14) << "Proportions: " << strExportProp << "\n";
(*writeTo) << setw(14) << "Haplotypes: " << strExportHap << "\n";
if ( doExportVcf() ) { (*writeTo) << setw(14) << "Vcf: " << strExportVcf << "\n"; }
if (this->useIBD()){
(*writeTo) << " IBD method output saved to:\n";
(*writeTo) << setw(14) << "Likelihood: " << strIbdExportLLK << "\n";
(*writeTo) << setw(14) << "Proportions: " << strIbdExportProp << "\n";
(*writeTo) << setw(14) << "Haplotypes: " << strIbdExportHap << "\n";
}
if (this->ibdProbsIntegrated.size()>1){
(*writeTo) << setw(14) << "IBD probs: " << strIbdExportProbs << "\n\n";
(*writeTo) << " IBD probabilities:\n";
for ( size_t stateI = 0; stateI < this->ibdProbsHeader.size(); stateI++ ){
(*writeTo) << setw(14) << this->ibdProbsHeader[stateI] << ": " << this->ibdProbsIntegrated[stateI] << "\n";
(*writeTo) << "Output saved to:\n";
if ( this->doLsPainting() ){
for ( size_t i = 0; i < kStrain(); i++ ){
(*writeTo) << "Posterior probability of strain " << i << ": "<< strExportSingleFwdProbPrefix << i <<endl;
}
} else if (this->doIbdPainting()){
if (this->ibdProbsIntegrated.size()>1){
(*writeTo) << setw(14) << "IBD probs: " << strIbdExportProbs << "\n\n";
(*writeTo) << " IBD probabilities:\n";
for ( size_t stateI = 0; stateI < this->ibdProbsHeader.size(); stateI++ ){
(*writeTo) << setw(14) << this->ibdProbsHeader[stateI] << ": " << this->ibdProbsIntegrated[stateI] << "\n";
}
}
} else {
(*writeTo) << setw(14) << "Likelihood: " << strExportLLK << "\n";
(*writeTo) << setw(14) << "Proportions: " << strExportProp << "\n";
(*writeTo) << setw(14) << "Haplotypes: " << strExportHap << "\n";
if ( doExportVcf() ) { (*writeTo) << setw(14) << "Vcf: " << strExportVcf << "\n"; }
if (this->useIBD()){
(*writeTo) << " IBD method output saved to:\n";
(*writeTo) << setw(14) << "Likelihood: " << strIbdExportLLK << "\n";
(*writeTo) << setw(14) << "Proportions: " << strIbdExportProp << "\n";
(*writeTo) << setw(14) << "Haplotypes: " << strIbdExportHap << "\n";
}
if (this->ibdProbsIntegrated.size()>1){
(*writeTo) << setw(14) << "IBD probs: " << strIbdExportProbs << "\n\n";
(*writeTo) << " IBD probabilities:\n";
for ( size_t stateI = 0; stateI < this->ibdProbsHeader.size(); stateI++ ){
(*writeTo) << setw(14) << this->ibdProbsHeader[stateI] << ": " << this->ibdProbsIntegrated[stateI] << "\n";
}
}
}
(*writeTo) << "\n";
(*writeTo) << " IBD best path llk: " << ibdLLK_ << "\n\n";

this->computeEffectiveKstrain(this->finalProp);
(*writeTo) << " Effective_K: " << this->effectiveKstrain_ <<"\n";
this->computeInferredKstrain(this->finalProp);
(*writeTo) << " Inferred_K: " << this->inferredKstrain_ <<"\n";
this->computeAdjustedEffectiveKstrain();
(*writeTo) << "Adjusted_effective_K: " << this->adjustedEffectiveKstrain_ <<"\n";
}
(*writeTo) << "\n";
(*writeTo) << " IBD best path llk: " << ibdLLK_ << "\n\n";

this->computeEffectiveKstrain(this->finalProp);
(*writeTo) << " Effective_K: " << this->effectiveKstrain_ <<"\n";
this->computeInferredKstrain(this->finalProp);
(*writeTo) << " Inferred_K: " << this->inferredKstrain_ <<"\n";
this->computeAdjustedEffectiveKstrain();
(*writeTo) << "Adjusted_effective_K: " << this->adjustedEffectiveKstrain_ <<"\n";

(*writeTo) << "\n";
(*writeTo) << "Proportions:\n";
for ( size_t ii = 0; ii < this->finalProp.size(); ii++){
Expand Down

0 comments on commit ef0e393

Please sign in to comment.