Skip to content

Commit

Permalink
update a couple of things
Browse files Browse the repository at this point in the history
  • Loading branch information
@shajoezhu
shajoezhu committed Nov 16, 2018
1 parent f18b075 commit 4d73536
Show file tree
Hide file tree
Showing 5 changed files with 86 additions and 13 deletions.
4 changes: 4 additions & 0 deletions src/dEploidIO.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);

Expand All @@ -1041,13 +1042,16 @@ 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 <<newHeader.size()<<endl;

vector < vector <double> > newPanel;
for (size_t i = 0; i < dummy.reducedPanel.size(); i++) {
vector <double> tmpRow;
for (size_t j = 0; j < min(dummy.choiceIdx.size(), maxNumPanel); j++) {
tmpRow.push_back(dummy.reducedPanel[i][j]);
}
tmpRow.push_back(static_cast<int>(0));
newPanel.push_back(tmpRow);
}

Expand Down
78 changes: 69 additions & 9 deletions src/mcmc.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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_);
Expand Down Expand Up @@ -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 ("<<this->mcmcJob<<")"<<flush;
//#endif
//}
//this->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;
Expand All @@ -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
Expand Down Expand Up @@ -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();
}
}

Expand Down Expand Up @@ -600,7 +660,7 @@ vector <double> McmcMachinery::calcTmpTitre() {
}


void McmcMachinery::updateSingleHap() {
void McmcMachinery::updateSingleHap(Panel *useThisPanel) {
dout << " Update Single Hap "<<endl;
this->findUpdatingStrainSingle();

Expand All @@ -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 ) {
Expand All @@ -645,7 +705,7 @@ void McmcMachinery::updateSingleHap() {
}


void McmcMachinery::updatePairHaps() {
void McmcMachinery::updatePairHaps(Panel *useThisPanel) {
if ( this->kStrain() == 1 ) {
return;
}
Expand All @@ -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_);
Expand Down
8 changes: 4 additions & 4 deletions src/mcmc.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -204,12 +204,12 @@ class McmcMachinery {
vector <double> calcTmpTitre();
double deltaLLKs ( vector <double> &newLLKs );

void updateSingleHap();
void findUpdatingStrainSingle( );
void updateSingleHap(Panel *useThisPanel);
void findUpdatingStrainSingle();

void updatePairHaps();
void updatePairHaps(Panel *useThisPanel);
//vector <size_t> sampleNoReplace(MersenneTwister* rg, vector <double> & proportion, size_t nSample );
void findUpdatingStrainPair( );
void findUpdatingStrainPair();

/* Debug */
bool doutProp();
Expand Down
7 changes: 7 additions & 0 deletions src/panel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 );
Expand Down
2 changes: 2 additions & 0 deletions utilities/dEploidTools.r
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]
}
Expand Down

0 comments on commit 4d73536

Please sign in to comment.