Skip to content

Commit

Permalink
Merge pull request #333 from DEploid-dev/style
Browse files Browse the repository at this point in the history
check style
  • Loading branch information
shajoezhu authored Dec 28, 2020
2 parents f37d8e7 + 6eeaf11 commit 4ed652e
Show file tree
Hide file tree
Showing 6 changed files with 235 additions and 198 deletions.
2 changes: 2 additions & 0 deletions .ci/checkedList
Original file line number Diff line number Diff line change
Expand Up @@ -16,3 +16,5 @@ tests/unittest/test_vcfReader.cpp
src/chooseK.hpp
src/chooseK.cpp
tests/unittest/test_dEploidIO.cpp
src/panel.hpp
src/panel.cpp
6 changes: 2 additions & 4 deletions .ci/mycppfileList
Original file line number Diff line number Diff line change
@@ -1,12 +1,8 @@
src/dEploidIO.cpp
src/dEploidIO.hpp
src/mcmc.cpp
src/panel.cpp
src/panel.hpp
src/updateHap.cpp
src/updateHap.hpp
src/utility.cpp
src/utility.hpp
src/export/dEploidIOExport.cpp
src/export/dEploidIOExportPosteriorProb.cpp
src/export/writeMcmcRelated.cpp
Expand All @@ -18,3 +14,5 @@ tests/unittest/test_txtReader.cpp
tests/unittest/test_updatePairHap.cpp
tests/unittest/test_updateSingleHap.cpp
tests/unittest/test_utilities.cpp
src/utility.hpp
src/utility.cpp
192 changes: 107 additions & 85 deletions src/panel.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,7 @@ Panel::Panel(vector < double > pRec,
}


Panel::Panel(const Panel &copyFrom){
Panel::Panel(const Panel &copyFrom) {
nLoci_ = copyFrom.nLoci_;

chrom_ = vector <string> (copyFrom.chrom_.begin(), copyFrom.chrom_.end());
Expand Down Expand Up @@ -99,115 +99,128 @@ Panel::Panel(const Panel &copyFrom){
inbreedingPanelSize_ = copyFrom.inbreedingPanelSize_;

this->fileName_ = copyFrom.fileName_;
for (size_t i = 0; i < copyFrom.header_.size(); i++){
for (size_t i = 0; i < copyFrom.header_.size(); i++) {
header_.push_back(copyFrom.header_[i]);
}
}



void Panel::readFromFile( const char inchar[] ) {
this->readFromFileBase( inchar );
this->setTruePanelSize( this->nInfoLines_ );
this->setInbreedingPanelSize( this->truePanelSize() );
void Panel::readFromFile(const char inchar[]) {
this->readFromFileBase(inchar);
this->setTruePanelSize(this->nInfoLines_);
this->setInbreedingPanelSize(this->truePanelSize());
}


void Panel::checkForExceptions( size_t nLoci, string panelFileName ) {
if ( this->content_.size() != nLoci ) {
throw LociNumberUnequal( panelFileName );
void Panel::checkForExceptions(size_t nLoci, string panelFileName) {
if (this->content_.size() != nLoci) {
throw LociNumberUnequal(panelFileName);
}

if ( this->pRec_.size() != nLoci ) {
throw LociNumberUnequal( panelFileName );
if (this->pRec_.size() != nLoci) {
throw LociNumberUnequal(panelFileName);
}
return;
}


void Panel::computeRecombProbs( double averageCentimorganDistance, double G, bool useConstRecomb, double constRecombProb, bool forbidCopyFromSame ) {
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 );
assert(pRecRec_.size() == 0 );
assert(pRecNoRec_.size() == 0 );
assert(pNoRecNoRec_.size() == 0 );
assert(pRec_.size() == 0);
assert(pRecEachHap_.size() == 0);
assert(pNoRec_.size() == 0);
assert(pRecRec_.size() == 0);
assert(pRecNoRec_.size() == 0);
assert(pNoRecNoRec_.size() == 0);

double averageMorganDistance = averageCentimorganDistance * 100;
double geneticDistance;
double rho;
double nPanelDouble = (double)this->truePanelSize();
double nPanelDouble = static_cast<double>(this->truePanelSize());
double nPanlelMinus1 = nPanelDouble - 1.0;
for ( size_t i = 0; i < this->position_.size(); i++) {
for ( size_t j = 1; j < this->position_[i].size(); j++) {
geneticDistance = (double)(this->position_[i][j] - this->position_[i][j-1])/averageMorganDistance ;
//rho = geneticDistance * 2 * Ne;
for (size_t i = 0; i < this->position_.size(); i++) {
for (size_t j = 1; j < this->position_[i].size(); j++) {
geneticDistance = static_cast<double>(
this->position_[i][j] - this->position_[i][j-1]) /
averageMorganDistance;
// rho = geneticDistance * 2 * Ne;
rho = geneticDistance * G;

double pRecTmp = ( useConstRecomb ) ? constRecombProb : 1.0 - exp(-rho);
this->pRec_.push_back( pRecTmp );
double pRecTmp = (useConstRecomb) ?
constRecombProb : 1.0 - exp(-rho);
this->pRec_.push_back(pRecTmp);

double pRecEachHapTmp = pRecTmp / nPanelDouble;
this->pRecEachHap_.push_back( pRecTmp / nPanelDouble );
this->pRecEachHap_.push_back(pRecTmp / nPanelDouble);

double pNoRecTmp = 1.0 - pRecTmp;
this->pNoRec_.push_back( pNoRecTmp );
this->pNoRec_.push_back(pNoRecTmp);

double secondPRecEachHapTmp = ( forbidCopyFromSame ) ? (pRecTmp / nPanlelMinus1) : pRecEachHapTmp; // allowing copy from the same
double secondPRecEachHapTmp = (forbidCopyFromSame) ?
(pRecTmp / nPanlelMinus1) :
pRecEachHapTmp; // allowing copy from the same

this->pRecRec_.push_back ( pRecEachHapTmp * secondPRecEachHapTmp );
this->pRecNoRec_.push_back ( secondPRecEachHapTmp * pNoRecTmp );
this->pNoRecNoRec_.push_back ( pNoRecTmp * pNoRecTmp );
this->pRecRec_.push_back(pRecEachHapTmp * secondPRecEachHapTmp);
this->pRecNoRec_.push_back(secondPRecEachHapTmp * pNoRecTmp);
this->pNoRecNoRec_.push_back(pNoRecTmp * pNoRecTmp);
}
this->pRec_.push_back(1.0);
this->pRecEachHap_.push_back( 1.0 / nPanelDouble );
this->pRecEachHap_.push_back(1.0 / nPanelDouble);
this->pNoRec_.push_back(0.0);
this->pRecRec_.push_back ( ( ( forbidCopyFromSame ) ? (1.0 / nPanelDouble / nPanlelMinus1) : (1.0 / nPanelDouble / nPanelDouble) ) );
this->pRecNoRec_.push_back ( 0.0 );
this->pNoRecNoRec_.push_back ( 0.0 );
this->pRecRec_.push_back(((forbidCopyFromSame) ?
(1.0 / nPanelDouble / nPanlelMinus1) :
(1.0 / nPanelDouble / nPanelDouble)));
this->pRecNoRec_.push_back(0.0);
this->pNoRecNoRec_.push_back(0.0);
}
assert(pRec_.size() == this->nLoci_ );
assert(pRecEachHap_.size() == this->nLoci_ );
assert(pNoRec_.size() == this->nLoci_ );
assert(pRecRec_.size() == this->nLoci_ );
assert(pRecNoRec_.size() == this->nLoci_ );
assert(pNoRecNoRec_.size() == this->nLoci_ );
assert(pRec_.size() == this->nLoci_);
assert(pRecEachHap_.size() == this->nLoci_);
assert(pNoRec_.size() == this->nLoci_);
assert(pRecRec_.size() == this->nLoci_);
assert(pRecNoRec_.size() == this->nLoci_);
assert(pNoRecNoRec_.size() == this->nLoci_);
}


void Panel::buildExamplePanel1() {
this->chrom_ = vector <string> ({"Pf3D7_01_v3"});
this->position_.push_back( vector <int> ({93157, 94422, 94459, 94487, 95518, 95632, 95641 }) );
this->position_.push_back(vector <int> (
{93157, 94422, 94459, 94487, 95518, 95632, 95641}));
this->indexOfChromStarts_ = vector <size_t> ({0});
this->buildExamplePanelContent();
}


void Panel::buildExamplePanel2() {
this->chrom_ = vector <string> ({"Pf3D7_01_v3", "Pf3D7_02_v3", "Pf3D7_03_v3"});
this->position_.push_back( vector <int> ({93157}) );
this->position_.push_back( vector <int> ({94422, 94459, 94487, 95518, 95632}) );
this->position_.push_back( vector <int> ({95641 }) );
this->chrom_ = vector <string> (
{"Pf3D7_01_v3", "Pf3D7_02_v3", "Pf3D7_03_v3"});
this->position_.push_back(vector <int> ({93157}));
this->position_.push_back(vector <int> (
{94422, 94459, 94487, 95518, 95632}));
this->position_.push_back(vector <int> ({95641}));
this->indexOfChromStarts_ = vector <size_t> ({0, 1, 6});
this->buildExamplePanelContent();
}


void Panel::buildExamplePanelContent() {
this->content_.push_back( vector <double> ({0,0,0,1}) );
this->content_.push_back( vector <double> ({0,0,0,1}) );
this->content_.push_back( vector <double> ({0,0,0,1}) );
this->content_.push_back( vector <double> ({0,0,0,1}) );
this->content_.push_back( vector <double> ({0,1,1,0}) );
this->content_.push_back( vector <double> ({0,0,1,0}) );
this->content_.push_back( vector <double> ({0,0,1,0}) );
this->content_.push_back(vector <double> ({0, 0, 0, 1}));
this->content_.push_back(vector <double> ({0, 0, 0, 1}));
this->content_.push_back(vector <double> ({0, 0, 0, 1}));
this->content_.push_back(vector <double> ({0, 0, 0, 1}));
this->content_.push_back(vector <double> ({0, 1, 1, 0}));
this->content_.push_back(vector <double> ({0, 0, 1, 0}));
this->content_.push_back(vector <double> ({0, 0, 1, 0}));
this->nLoci_ = this->content_.size();
this->nInfoLines_ = this->content_.back().size();
this->setTruePanelSize(this->nInfoLines_);
Expand All @@ -218,41 +231,47 @@ void Panel::buildExamplePanelContent() {
void Panel::initializeUpdatePanel(size_t inbreedingPanelSizeSetTo) {
this->setInbreedingPanelSize(inbreedingPanelSizeSetTo);

// If allows inbreeding, update reference panel by including strain haplotypes
// If allows inbreeding,
// update reference panel by including strain haplotypes
dout << "************* Allow inbreeding ************" << endl;
dout << "** Initialize inbreeding reference panel **" << endl;

if ( this->truePanelSize() == this->inbreedingPanelSize()) {
if (this->truePanelSize() == this->inbreedingPanelSize()) {
return;
}

for ( size_t siteI = 0; siteI < this->content_.size(); siteI++ ) {
for ( size_t panelStrainJ = this->truePanelSize() ; panelStrainJ < this->inbreedingPanelSize(); panelStrainJ++ ) {
this->content_[siteI].push_back( 1 );
for (size_t panelStrainJ = this->truePanelSize();
panelStrainJ < this->inbreedingPanelSize(); panelStrainJ++) {
this->content_[siteI].push_back(1);
}
assert(inbreedingPanelSizeSetTo == this->content_[siteI].size());
}
}


void Panel::updatePanelWithHaps(size_t inbreedingPanelSizeSetTo, size_t excludedStrain, vector < vector<double> > & haps) {
void Panel::updatePanelWithHaps(size_t inbreedingPanelSizeSetTo,
size_t excludedStrain, const vector < vector<double> > & haps) {
this->setInbreedingPanelSize(inbreedingPanelSizeSetTo);

// If allows inbreeding, update reference panel by including strain haplotypes
// If allows inbreeding,
// update reference panel by including strain haplotypes
dout << "*************** Allow inbreeding **************" << endl;
dout << "*** Update reference panel without strain " << excludedStrain << " ***" << endl;
dout << "*** Update reference panel without strain " <<
excludedStrain << " ***" << endl;

if ( this->truePanelSize() == this->inbreedingPanelSize()) {
if (this->truePanelSize() == this->inbreedingPanelSize()) {
return;
}

for ( size_t siteI = 0; siteI < this->content_.size(); siteI++ ) {
for (size_t siteI = 0; siteI < this->content_.size(); siteI++) {
size_t shiftAfter = this->inbreedingPanelSize();

for ( size_t panelStrainJ = this->truePanelSize() ; panelStrainJ < this->inbreedingPanelSize(); panelStrainJ++ ) {
for (size_t panelStrainJ = this->truePanelSize();
panelStrainJ < this->inbreedingPanelSize(); panelStrainJ++) {
size_t strainIndex = panelStrainJ - this->truePanelSize();

if ( strainIndex == excludedStrain ) {
if (strainIndex == excludedStrain) {
shiftAfter = panelStrainJ;
}

Expand All @@ -265,14 +284,13 @@ void Panel::updatePanelWithHaps(size_t inbreedingPanelSizeSetTo, size_t excluded
}


void Panel::trimVec(vector <double> &vec, const vector <size_t> &idx) {
void Panel::trimVec(vector <double> vec, const vector <size_t> &idx) {
vector <double> ret;
for (auto const& value : idx){
for (auto const& value : idx) {
ret.push_back(vec[value]);
}
//return ret;
vec.clear();
for (auto const& value : ret){
for (auto const& value : ret) {
vec.push_back(value);
}
}
Expand Down Expand Up @@ -309,32 +327,36 @@ void Panel::findAndKeepMarkersGivenIndexHalf(
this->trimVec(this->pNoRecNoRec_, givenIndex);
}

void IBDrecombProbs::computeRecombProbs( double averageCentimorganDistance, double G, bool useConstRecomb, double constRecombProb ) {
assert(pRec_.size() == 0 );
assert(pNoRec_.size() == 0 );
void IBDrecombProbs::computeRecombProbs(double averageCentimorganDistance,
double G, bool useConstRecomb, double constRecombProb) {
assert(pRec_.size() == 0);
assert(pNoRec_.size() == 0);
double averageMorganDistance = averageCentimorganDistance * 100;
double geneticDistance;
double rho;
for ( size_t i = 0; i < this->position_.size(); i++) {
for ( size_t j = 1; j < this->position_[i].size(); j++) {
geneticDistance = (double)(this->position_[i][j] - this->position_[i][j-1])/averageMorganDistance ;
//rho = geneticDistance * 2 * Ne;
for (size_t i = 0; i < this->position_.size(); i++) {
for (size_t j = 1; j < this->position_[i].size(); j++) {
geneticDistance =
static_cast<double>(this->position_[i][j] -
this->position_[i][j-1]) / averageMorganDistance;
// rho = geneticDistance * 2 * Ne;
rho = geneticDistance * G;
double pRecTmp = ( useConstRecomb ) ? constRecombProb : (1.0 - exp(-rho));
this->pRec_.push_back( pRecTmp );
double pRecTmp = (useConstRecomb) ?
constRecombProb : (1.0 - exp(-rho));
this->pRec_.push_back(pRecTmp);
double pNoRecTmp = 1.0 - pRecTmp;
this->pNoRec_.push_back( pNoRecTmp );
this->pNoRec_.push_back(pNoRecTmp);
}
this->pRec_.push_back(1.0);
this->pNoRec_.push_back(0.0);
}
assert(pRec_.size() == this->nLoci_ );
assert(pNoRec_.size() == this->nLoci_ );
assert(pRec_.size() == this->nLoci_);
assert(pNoRec_.size() == this->nLoci_);
}


//vector<vector<double>> outtrans(out[0].size(),
//vector<double>(out.size()));
//for (size_t i = 0; i < out.size(); ++i)
//for (size_t j = 0; j < out[0].size(); ++j)
//outtrans[j][i] = out[i][j];
// vector<vector<double>> outtrans(out[0].size(),
// vector<double>(out.size()));
// for (size_t i = 0; i < out.size(); ++i)
// for (size_t j = 0; j < out[0].size(); ++j)
// outtrans[j][i] = out[i][j];
Loading

0 comments on commit 4ed652e

Please sign in to comment.