diff --git a/src/dEploidIO.cpp b/src/dEploidIO.cpp index c4309c9e..f57d32ab 100644 --- a/src/dEploidIO.cpp +++ b/src/dEploidIO.cpp @@ -1032,6 +1032,7 @@ void DEploidIO::dEploidLassoFullPanel() { this->computeObsWsaf(); Panel tmpPanel(*panel); + tmpPanel.computeRecombProbs(this->averageCentimorganDistance(), this->parameterG(), true, 0.0000001, this->forbidCopyFromSame()); tmpPanel.findAndKeepMarkersGivenIndex(this->vcfReaderPtr_->legitVqslodAt); DEploidLASSO dummy(tmpPanel.content_, this->obsWsaf_, 250); @@ -1041,6 +1042,8 @@ void DEploidIO::dEploidLassoFullPanel() { for (size_t i = 0; i < min(dummy.choiceIdx.size(), maxNumPanel); i++) { newHeader.push_back(panel->header_[dummy.choiceIdx[i]]); } +newHeader.push_back("3d7"); +//cout < > newPanel; for (size_t i = 0; i < dummy.reducedPanel.size(); i++) { @@ -1048,6 +1051,7 @@ void DEploidIO::dEploidLassoFullPanel() { for (size_t j = 0; j < min(dummy.choiceIdx.size(), maxNumPanel); j++) { tmpRow.push_back(dummy.reducedPanel[i][j]); } +tmpRow.push_back(static_cast(0)); newPanel.push_back(tmpRow); } diff --git a/src/mcmc.cpp b/src/mcmc.cpp index 64447134..142b9a85 100644 --- a/src/mcmc.cpp +++ b/src/mcmc.cpp @@ -105,7 +105,10 @@ void McmcMachinery::initializeMcmcChain(bool useIBD) { dout << "###########################################"<< endl; dout << "# Initialization #"<< endl; dout << "###########################################"<< endl; - cout << "plaf.size() " << plaf_ptr_->size() <<" " << "Panel size = " << this->panel_->truePanelSize() << endl; + cout << "plaf.size() " << plaf_ptr_->size() <<" "; + if (panel_ != NULL) { + cout << "Panel size = " << this->panel_->truePanelSize() << endl; + } this->initializeTitre(); this->currentLogPriorTitre_ = this->calcLogPriorTitre(this->currentTitre_); this->initializeHap(); @@ -205,7 +208,7 @@ void McmcMachinery::initializellk() { } -void McmcMachinery::initializeProp( ) { +void McmcMachinery::initializeProp() { assert( this->currentProp_.size() == (size_t)0 ); this->currentProp_ = ( this->dEploidIO_ -> initialPropWasGiven()) ? this->dEploidIO_ ->initialProp: @@ -237,6 +240,7 @@ void McmcMachinery::initializeTitre() { if ( this->dEploidIO_->doUpdateProp() ) { for ( size_t k = 0; k < this->kStrain_; k++) { this->currentTitre_[k] = this->initialTitreNormalVariable(); + //this->currentTitre_[k] = 10.0; } } assert( currentTitre_.size() == this->kStrain_); @@ -281,6 +285,46 @@ void McmcMachinery::runMcmcChain( bool showProgress, bool useIBD, bool notInR ) cout << value << " "; } cout << endl; + + //if (useIBD != true){ + //#ifndef RBUILD + //clog << "Pretraining ..." << endl; + //#endif + + //for ( this->currentMcmcIteration_ = 0 ; currentMcmcIteration_ < this->mcmcThresh_ ; currentMcmcIteration_++) { + //if ( this->currentMcmcIteration_ > 0 && this->currentMcmcIteration_%30 == 0 && showProgress ) { + //#ifndef RBUILD + //clog << "\r" << " MCMC step" << setw(4) << int(currentMcmcIteration_ * 100 / this->mcmcThresh_) << "% completed ("<mcmcJob<<")"<sampleMcmcEvent(useIBD); + //for (auto const& vv : this->currentProp_){ + //cout << vv << " "; + //} + //cout << endl; + + //} + + //for (auto const& value : this->currentProp_) { + //if ((value > 0.01) & (value < 0.2)) { + //cout << "Switch sigma " << endl; + //this->currentTitre_.clear(); + //this->initializeTitre(); + //this->currentProp_.clear(); + //this->initializeProp( ); + //this->initializeExpectedWsaf(); // This requires currentHap_ and currentProp_ + + //for (auto const& vv : this->currentProp_){ + //cout << vv << " "; + //} + //cout << endl; + //this->SD_LOG_TITRE = 10.0; + //break; + //} + //} + + //} + for ( this->currentMcmcIteration_ = 0 ; currentMcmcIteration_ < this->maxIteration_ ; currentMcmcIteration_++) { dout << endl; dout << "MCMC iteration: " << this->currentMcmcIteration_ << endl; @@ -290,6 +334,10 @@ void McmcMachinery::runMcmcChain( bool showProgress, bool useIBD, bool notInR ) #endif } this->sampleMcmcEvent(useIBD); + //for (auto const& vv : this->currentProp_){ + //cout << vv << " "; + //} + //cout << endl; } #ifndef RBUILD @@ -409,13 +457,25 @@ void McmcMachinery::sampleMcmcEvent( bool useIBD ) { ibdSampleMcmcEventStep(); assert(doutProp()); } else { - this->eventInt_ = this->mcmcEventRg_->sampleInt(3); + this->eventInt_ = this->mcmcEventRg_->sampleInt(4); if ( (this->eventInt_ == 0) && (this->dEploidIO_->doUpdateProp() == true) ) { this->updateProportion(); } else if ( (this->eventInt_ == 1) && (this->dEploidIO_->doUpdateSingle() == true) ) { - this->updateSingleHap(); + this->updateSingleHap(this->panel_); } else if ( (this->eventInt_ == 2) && (this->dEploidIO_->doUpdatePair() == true) ) { - this->updatePairHaps(); + this->updatePairHaps(this->panel_); + } else if ( (this->eventInt_ == 3) && (this->dEploidIO_->doUpdateSingle() == true) ) { + this->updateSingleHap(NULL); + this->updateSingleHap(NULL); + this->updateSingleHap(NULL); + this->updateSingleHap(NULL); + //this->updatePairHaps(NULL); + //this->updatePairHaps(this->panel_); + //} else if ( (this->eventInt_ == 4) && (this->dEploidIO_->doUpdatePair() == true) ) { + ////this->updatePairHaps(NULL); + //this->updateSingleHap(NULL); + //this->updatePairHaps(this->panel_); + //this->updateProportion(); } } @@ -600,7 +660,7 @@ vector McmcMachinery::calcTmpTitre() { } -void McmcMachinery::updateSingleHap() { +void McmcMachinery::updateSingleHap(Panel *useThisPanel) { dout << " Update Single Hap "<findUpdatingStrainSingle(); @@ -618,7 +678,7 @@ void McmcMachinery::updateSingleHap() { this->currentExpectedWsaf_, this->currentProp_, this->currentHap_, this->hapRg_, start, length, - this->panel_, this->dEploidIO_->missCopyProb_, this->dEploidIO_->scalingFactor(), + useThisPanel, this->dEploidIO_->missCopyProb_, this->dEploidIO_->scalingFactor(), this->strainIndex_); if ( this->dEploidIO_->doAllowInbreeding() == true ) { @@ -645,7 +705,7 @@ void McmcMachinery::updateSingleHap() { } -void McmcMachinery::updatePairHaps() { +void McmcMachinery::updatePairHaps(Panel *useThisPanel) { if ( this->kStrain() == 1 ) { return; } @@ -664,7 +724,7 @@ void McmcMachinery::updatePairHaps() { this->currentExpectedWsaf_, this->currentProp_, this->currentHap_, this->hapRg_, start, length, - this->panel_, this->dEploidIO_->missCopyProb_, this->dEploidIO_->scalingFactor(), + useThisPanel, this->dEploidIO_->missCopyProb_, this->dEploidIO_->scalingFactor(), this->dEploidIO_->forbidCopyFromSame(), this->strainIndex1_, this->strainIndex2_); diff --git a/src/mcmc.hpp b/src/mcmc.hpp index b854a6e8..dd5228a7 100644 --- a/src/mcmc.hpp +++ b/src/mcmc.hpp @@ -204,12 +204,12 @@ class McmcMachinery { vector calcTmpTitre(); double deltaLLKs ( vector &newLLKs ); - void updateSingleHap(); - void findUpdatingStrainSingle( ); + void updateSingleHap(Panel *useThisPanel); + void findUpdatingStrainSingle(); - void updatePairHaps(); + void updatePairHaps(Panel *useThisPanel); //vector sampleNoReplace(MersenneTwister* rg, vector & proportion, size_t nSample ); - void findUpdatingStrainPair( ); + void findUpdatingStrainPair(); /* Debug */ bool doutProp(); diff --git a/src/panel.cpp b/src/panel.cpp index 3f6493ef..1c2dbb45 100644 --- a/src/panel.cpp +++ b/src/panel.cpp @@ -126,6 +126,13 @@ void Panel::checkForExceptions( size_t nLoci, string panelFileName ) { void Panel::computeRecombProbs( double averageCentimorganDistance, double G, bool useConstRecomb, double constRecombProb, bool forbidCopyFromSame ) { + pRec_.clear(); + pRecEachHap_.clear(); + pNoRec_.clear(); + pRecRec_.clear(); + pRecNoRec_.clear(); + pNoRecNoRec_.clear(); + assert(pRec_.size() == 0 ); assert(pRecEachHap_.size() == 0 ); assert(pNoRec_.size() == 0 ); diff --git a/utilities/dEploidTools.r b/utilities/dEploidTools.r index d413250e..0f2cb9af 100644 --- a/utilities/dEploidTools.r +++ b/utilities/dEploidTools.r @@ -400,6 +400,8 @@ fun.interpretDEploid.LassoIBD <- function (coverage, PLAF, dEploidPrefix, prefix includeindex = which(!excludeLogic) obsWSAF = obsWSAF[includeindex] PLAF = PLAF[includeindex] + includeLogic.lassoK = which( paste(coverage$CHROM[includeindex], coverage$POS[includeindex]) %in% paste(hap.lassoK$CHROM, hap.lassoK$POS) ) + includeLogic.ibd = which( paste(coverage$CHROM[includeindex], coverage$POS[includeindex]) %in% paste(hap.ibd$CHROM, hap.ibd$POS) ) # ref = ref[includeindex] # alt = alt[includeindex] }