Skip to content

Commit

Permalink
make src/dEploidIO-workflow.cpp
Browse files Browse the repository at this point in the history
  • Loading branch information
@shajoezhu
shajoezhu committed Jan 1, 2021
1 parent 6323f7c commit 4409c03
Show file tree
Hide file tree
Showing 7 changed files with 300 additions and 251 deletions.
1 change: 1 addition & 0 deletions .ci/checkedList
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,3 +22,4 @@ tests/unittest/test_runner.cpp
tests/unittest/test_utilities.cpp
src/export/dEploidIOExportPosteriorProb.cpp
src/export/writeMcmcRelated.cpp
src/dEploidIO-workflow.cpp
1 change: 1 addition & 0 deletions Makefile.am
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@ common_src = src/random/fastfunc.cpp \
src/panel.cpp \
src/utility.cpp \
src/dEploidIO.cpp \
src/dEploidIO-workflow.cpp \
src/updateHap.cpp \
src/txtReader.cpp \
src/vcfReader.cpp \
Expand Down
252 changes: 3 additions & 249 deletions src/dEploid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,6 @@
*/

#include <iostream> // std::cout
#include "mcmc.hpp"
#include "dEploidIO.hpp"

int main(int argc, char *argv[]) {
Expand All @@ -51,256 +50,11 @@ int main(int argc, char *argv[]) {
} else if (dEploidIO.doIbdViterbiPainting()) {
dEploidIO.paintIBDviterbi();
} else if (dEploidIO.useLasso()) { // DEploid-Lasso
dEploidIO.dEploidLasso();
MersenneTwister lassoRg(dEploidIO.randomSeed());
DEploidIO tmpIO(dEploidIO);
vector < vector <double> > hap;
for (size_t chromi = 0;
chromi < dEploidIO.indexOfChromStarts_.size();
chromi++ ) {
tmpIO.position_.clear();
tmpIO.position_.push_back(dEploidIO.position_.at(chromi));
tmpIO.indexOfChromStarts_.clear();
tmpIO.indexOfChromStarts_.push_back(0);
string job = string("DEploid-Lasso learning chromosome ");
job.append(dEploidIO.chrom_[chromi]).append(" haplotypes");
McmcSample * lassoMcmcSample = new McmcSample();
McmcMachinery lassoMcmcMachinery(
&dEploidIO.lassoPlafs.at(chromi),
&dEploidIO.lassoRefCount.at(chromi),
&dEploidIO.lassoAltCount.at(chromi),
dEploidIO.lassoPanels.at(chromi),
&tmpIO,
job,
"lasso",
lassoMcmcSample,
&lassoRg,
false);
lassoMcmcMachinery.runMcmcChain(true, // show progress
false); // use IBD
for (size_t snpi = 0;
snpi < lassoMcmcSample->hap.size(); snpi++) {
hap.push_back(vector <double> (
lassoMcmcSample->hap[snpi].begin(),
lassoMcmcSample->hap[snpi].end()));
}
delete lassoMcmcSample;
}
dEploidIO.writeHap(hap, "lasso");
dEploidIO.workflow_lasso();
} else if (dEploidIO.useBestPractice()) { // best practice
MersenneTwister rg(dEploidIO.randomSeed());

// #############################################################
// ################# DEploid-LASSO to learn K ##################
// #############################################################

vector < vector <double> > chooseKhap;

DEploidIO toLearnK(dEploidIO);
toLearnK.dEploidLassoTrimfirst();
DEploidIO toLearnKtmp(dEploidIO);
for (size_t chromi = 0;
chromi < toLearnK.indexOfChromStarts_.size();
chromi++ ) {
bool tryAgain = true;
while (tryAgain) {
toLearnKtmp.position_.clear();
toLearnKtmp.position_.push_back(
toLearnK.position_.at(chromi));
toLearnKtmp.indexOfChromStarts_.clear();
toLearnKtmp.indexOfChromStarts_.push_back(0);

McmcSample * lassoMcmcSample = new McmcSample();
string job = string("DEploid-Lasso learning K ");
job.append(toLearnK.chrom_[chromi]).append(" leanrning K");
McmcMachinery lassoMcmcMachinery(
&toLearnK.lassoPlafs.at(chromi),
&toLearnK.lassoRefCount.at(chromi),
&toLearnK.lassoAltCount.at(chromi),
toLearnK.lassoPanels.at(chromi),
&toLearnKtmp,
job,
job,
lassoMcmcSample,
&rg,
false);
lassoMcmcMachinery.runMcmcChain(true, // show progress
false, // use IBD
true, // notInR
false); // averageP
toLearnKtmp.initialProp = toLearnKtmp.finalProp;
toLearnKtmp.setInitialPropWasGiven(true);
if (toLearnKtmp.acceptRatio() != 0) {
dEploidIO.chooseK.appendProportions(toLearnKtmp.finalProp);
tryAgain = false;
} else {
toLearnKtmp.setInitialPropWasGiven(false);
tryAgain = true;
}

if (!tryAgain) {
for (size_t snpi = 0;
snpi < lassoMcmcSample->hap.size(); snpi++) {
chooseKhap.push_back(vector <double> (
lassoMcmcSample->hap[snpi].begin(),
lassoMcmcSample->hap[snpi].end()));
}
}
delete lassoMcmcSample;
}
}
// Gathering haplotype information
toLearnK.writeHap(chooseKhap, "chooseK");

// Gathering proportion information
dEploidIO.initialProp.clear();
vector <double> initialP = dEploidIO.chooseK.chosenP();
for (auto const& value : initialP) {
if (value > 0.01) {
dEploidIO.initialProp.push_back(value);
}
}
dEploidIO.setKstrain(dEploidIO.initialProp.size());
dEploidIO.setInitialPropWasGiven(true);

if (dEploidIO.inferBestPracticeP()) {
if (dEploidIO.initialProp.size() > 1) {
// #############################################################
// ###################### DEploid-IBD ########################
// #############################################################
McmcSample * ibdMcmcSample = new McmcSample();
DEploidIO tmpIO2(dEploidIO);
tmpIO2.ibdTrimming();
McmcMachinery ibdMcmcMachinery(&tmpIO2.plaf_,
&tmpIO2.refCount_,
&tmpIO2.altCount_,
tmpIO2.panel,
&tmpIO2,
string("DEploid-IBD learning proportion"),
string("ibd"),
ibdMcmcSample,
&rg,
true);
ibdMcmcMachinery.runMcmcChain(true, // show progress
true); // use IBD
// if (dEploidIO.useIbdOnly()) {
// tmpIO.paintIBD();
// dEploidIO.finalProp = tmpIO.initialProp;
// }

tmpIO2.paintIBD();

vector <double> initialP;
for (auto const& value : tmpIO2.finalProp) {
if (value > 0.01) {
initialP.push_back(value);
}
}
cout << endl;
dEploidIO.initialProp = initialP;
dEploidIO.finalProp = initialP;
dEploidIO.setKstrain(initialP.size());
dEploidIO.setInitialPropWasGiven(true);
dEploidIO.setDoUpdateProp(false);
delete ibdMcmcSample;
} else {
dEploidIO.finalProp.clear();
dEploidIO.finalProp.push_back(1.0);
dEploidIO.setKstrain(1);
}
}

if (dEploidIO.inferBestPracticeHap()) {
// #################################################################
// ###################### DEploid-LASSO ############################
// #################################################################
// *** Frist split the reference panel etc
dEploidIO.dEploidLasso();

DEploidIO dEploidLassoIO(dEploidIO);
dEploidLassoIO.initialProp = dEploidIO.initialProp;
dEploidLassoIO.setDoUpdateProp(false);
dEploidLassoIO.setInitialPropWasGiven(true);
dEploidLassoIO.setKstrain(dEploidIO.kStrain());
vector < vector <double> > hap;
for (size_t chromi = 0;
chromi < dEploidIO.indexOfChromStarts_.size();
chromi++ ) {
dEploidLassoIO.position_.clear();
dEploidLassoIO.position_.push_back(
dEploidIO.position_.at(chromi));
dEploidLassoIO.indexOfChromStarts_.clear();
dEploidLassoIO.indexOfChromStarts_.push_back(0);

McmcSample * lassoMcmcSample = new McmcSample();
string job = string("DEploid-Lasso learning chromosome ");
job.append(dEploidIO.chrom_[chromi]).append(" haplotypes");
McmcMachinery lassoMcmcMachinery(
&dEploidIO.lassoPlafs.at(chromi),
&dEploidIO.lassoRefCount.at(chromi),
&dEploidIO.lassoAltCount.at(chromi),
dEploidIO.lassoPanels.at(chromi),
&dEploidLassoIO,
job,
job,
lassoMcmcSample,
&rg,
false);
lassoMcmcMachinery.runMcmcChain(true, // show progress
false); // use IBD
// *** Append haplotypes to the end
for (size_t snpi = 0;
snpi < lassoMcmcSample->hap.size(); snpi++) {
hap.push_back(vector <double> (
lassoMcmcSample->hap[snpi].begin(),
lassoMcmcSample->hap[snpi].end()));
}
delete lassoMcmcSample;
}
dEploidIO.writeHap(hap, "final");
dEploidIO.writeVcf(hap, dEploidLassoIO.initialProp, "final");
}
if (dEploidIO.inferBestPracticeP() &
(dEploidIO.initialProp.size() > 1)) {
dEploidIO.paintIBD();
}

dEploidIO.workflow_best();
} else { // classic version, and DEploid-IBD
if (dEploidIO.useIBD()) { // ibd
McmcSample * ibdMcmcSample = new McmcSample();
MersenneTwister ibdRg(dEploidIO.randomSeed());
McmcMachinery ibdMcmcMachinery(&dEploidIO.plaf_,
&dEploidIO.refCount_,
&dEploidIO.altCount_,
dEploidIO.panel,
&dEploidIO,
"DEploid-IBD",
"ibd",
ibdMcmcSample,
&ibdRg,
true);
ibdMcmcMachinery.runMcmcChain(true, // show progress
true); // use IBD
delete ibdMcmcSample;
}

McmcSample * mcmcSample = new McmcSample();
MersenneTwister rg(dEploidIO.randomSeed());
McmcMachinery mcmcMachinery(&dEploidIO.plaf_,
&dEploidIO.refCount_,
&dEploidIO.altCount_,
dEploidIO.panel,
&dEploidIO,
"DEploid classic version",
"classic", // brief
mcmcSample,
&rg,
false); // use IBD
mcmcMachinery.runMcmcChain(true, // show progress
false); // use IBD
dEploidIO.paintIBD();
dEploidIO.writeHap(mcmcSample->hap, "final");
delete mcmcSample;
dEploidIO.workflow_ibd();
}
// Finishing, write log
dEploidIO.wrapUp();
Expand Down
Loading

0 comments on commit 4409c03

Please sign in to comment.