diff --git a/src/dEploidIO.cpp b/src/dEploidIO.cpp index 79f978e2..aad28b4c 100644 --- a/src/dEploidIO.cpp +++ b/src/dEploidIO.cpp @@ -94,6 +94,7 @@ void DEploidIO::core() { void DEploidIO::init() { + this->setDoPrintLassoPanel(false); this->setIsCopied(false); this->setDoExportRecombProb(false); this->setrandomSeedWasGiven(false); @@ -450,6 +451,8 @@ void DEploidIO::parse () { } else if ( *argv_i == "-lasso" ) { this->setUseLasso(true); this->setDoUpdateProp(false); + } else if ( *argv_i == "-writePanel" ) { + this->setDoPrintLassoPanel(true); } else if ( *argv_i == "-computeLLK" ) { this->setDoComputeLLK( true ); } else if ( *argv_i == "-ibdPainting" ) { @@ -592,7 +595,7 @@ void DEploidIO::printHelp(std::ostream& out) { out << "./dEploid -vcf data/testData/PG0390-C.test.vcf -exclude data/testData/labStrains.test.exclude.txt -plaf data/testData/labStrains.test.PLAF.txt -o PG0390-CPanelExclude -panel data/testData/labStrains.test.panel.txt -painting PG0390-CPanelExclude.hap" << endl; out << "./dEploid -vcf data/testData/PG0390-C.test.vcf -plaf data/testData/labStrains.test.PLAF.txt -o PG0390-CNopanel -noPanel -k 2 -ibd -nSample 250 -rate 8 -burn 0.67" < wsaf = lassoComputeObsWsaf(start, length); vector < vector > tmpPanel = lassoSubsetPanel(start, length); DEploidLASSO dummy(tmpPanel, wsaf, 250); - + vector newHeader; + for (size_t i = 0; i < dummy.choiceIdx.size(); i++) { + newHeader.push_back(panel->header_[dummy.choiceIdx[i]]); + } Panel * tmp = new Panel(vecFromTo(this->panel->pRec_, start, end), vecFromTo(this->panel->pRecEachHap_, start, end), vecFromTo(this->panel->pNoRec_, start, end), @@ -890,7 +896,9 @@ void DEploidIO::dEploidLasso() { lassoPlafs.push_back(vecFromTo(plaf_, start, end)); lassoRefCount.push_back(vecFromTo(refCount_, start, end)); lassoAltCount.push_back(vecFromTo(altCount_, start, end)); - this->writePanel(tmp, chromi); + if (this->doPrintLassoPanel()) { + this->writePanel(tmp, chromi, newHeader); + } } for ( auto const& value: this->initialProp ) { this->finalProp.push_back(value); diff --git a/src/dEploidIO.hpp b/src/dEploidIO.hpp index d047f9bc..ce5592f7 100644 --- a/src/dEploidIO.hpp +++ b/src/dEploidIO.hpp @@ -102,6 +102,7 @@ class DEploidIO{ vector < vector > lassoAltCount; void writeHap (vector < vector > &hap, bool useIBD = false); + bool doPrintLassoPanel_; private: void core(); @@ -279,6 +280,9 @@ class DEploidIO{ } // Getters and Setters + void setDoPrintLassoPanel ( const bool setTo ) { this->doPrintLassoPanel_ = setTo; } + bool doPrintLassoPanel() const { return this->doPrintLassoPanel_; } + void setIsCopied ( const bool setTo ) { this->isCopied_ = setTo; } bool isCopied() const { return this->isCopied_; } @@ -401,7 +405,7 @@ class DEploidIO{ // Lasso related vector lassoComputeObsWsaf(size_t segmentStartIndex, size_t nLoci); vector < vector > lassoSubsetPanel(size_t segmentStartIndex, size_t nLoci); - void writePanel(Panel *panel, size_t chromi); + void writePanel(Panel *panel, size_t chromi, vector hdr); }; #endif diff --git a/src/export/writeMcmcRelated.cpp b/src/export/writeMcmcRelated.cpp index 2c863bde..5146fa9b 100644 --- a/src/export/writeMcmcRelated.cpp +++ b/src/export/writeMcmcRelated.cpp @@ -114,29 +114,28 @@ void DEploidIO::writeHap(vector < vector > &hap, bool useIBD){ } -void DEploidIO::writePanel(Panel *panel, size_t chromI){ +void DEploidIO::writePanel(Panel *panel, size_t chromI, vector hdr){ string strExport = this->prefix_ + ".panel." + to_string(chromI); + remove(strExport.c_str()); + ofstreamExportTmp.open( strExport.c_str(), ios::out | ios::app | ios::binary ); - // TODO need to redo header // HEADER ofstreamExportTmp << "CHROM" << "\t" << "POS" << "\t";; for ( size_t ii = 0; ii < panel->truePanelSize_; ii++){ - ofstreamExportTmp << "h" << (ii+1) ; + ofstreamExportTmp << hdr[ii]; ofstreamExportTmp << ((ii < (panel->truePanelSize_-1)) ? "\t" : "\n") ; } size_t siteIndex = 0; - //for ( size_t chromI = 0; chromI < chrom_.size(); chromI++ ){ - for ( size_t posI = 0; posI < position_[chromI].size(); posI++){ - ofstreamExportTmp << chrom_[chromI] << "\t" << (int)position_[chromI][posI] << "\t"; - for ( size_t ii = 0; ii < panel->content_[siteIndex].size(); ii++){ - ofstreamExportTmp << panel->content_[siteIndex][ii]; - ofstreamExportTmp << ((ii < (panel->content_[siteIndex].size()-1)) ? "\t" : "\n") ; - } - siteIndex++; + for ( size_t posI = 0; posI < position_[chromI].size(); posI++){ + ofstreamExportTmp << chrom_[chromI] << "\t" << (int)position_[chromI][posI] << "\t"; + for ( size_t ii = 0; ii < panel->content_[siteIndex].size(); ii++){ + ofstreamExportTmp << panel->content_[siteIndex][ii]; + ofstreamExportTmp << ((ii < (panel->content_[siteIndex].size()-1)) ? "\t" : "\n") ; } - //} + siteIndex++; + } assert ( siteIndex == panel->content_.size()); ofstreamExportTmp.close(); diff --git a/src/txtReader.cpp b/src/txtReader.cpp index 4046267c..08ccea31 100644 --- a/src/txtReader.cpp +++ b/src/txtReader.cpp @@ -41,6 +41,7 @@ void TxtReader::readFromFileBase(const char inchar[]) { if (in_file.good()) { // skip the first line, which is the header getline(in_file, tmp_line); + this->extractHeader(tmp_line); getline(in_file, tmp_line); while (tmp_line.size() > 0) { size_t field_start = 0; @@ -91,6 +92,29 @@ void TxtReader::readFromFileBase(const char inchar[]) { } +void TxtReader::extractHeader(const string &line) { + this->header_.clear(); + size_t field_start = 0; + size_t field_end = 0; + size_t field_index = 0; + while (field_end < line.size()) { + field_end = min( + min(line.find(',', field_start), + line.find('\t', field_start)), + line.find('\n', field_start)); + + string tmp_str = line.substr(field_start, + field_end - field_start); + if (field_index > 1) { + this->header_.push_back(tmp_str); + } + + field_start = field_end+1; + field_index++; + } +} + + void TxtReader::extractChrom(const string & tmp_str) { if (tmpChromInex_ >= 0) { if (tmp_str != this->chrom_.back()) { diff --git a/src/txtReader.hpp b/src/txtReader.hpp index 06b7e99b..4fd2c732 100644 --- a/src/txtReader.hpp +++ b/src/txtReader.hpp @@ -53,6 +53,7 @@ class TxtReader : public VariantIndex { // info_ only refers to the first column of the content vector info_; + vector header_; size_t nInfoLines_; int tmpChromInex_; @@ -61,6 +62,7 @@ class TxtReader : public VariantIndex { // Methods void extractChrom(const string & tmp_str); void extractPOS(const string & tmp_str); + void extractHeader(const string &line); void reshapeContentToInfo(); string fileName; diff --git a/tests/unittest/test_dEploidIO.cpp b/tests/unittest/test_dEploidIO.cpp index 4d612bbd..560f6145 100644 --- a/tests/unittest/test_dEploidIO.cpp +++ b/tests/unittest/test_dEploidIO.cpp @@ -837,15 +837,11 @@ class TestIO : public CppUnit::TestCase { "-lasso" }; CPPUNIT_ASSERT_NO_THROW(DEploidIO(10, argv)); DEploidIO dEploidIOlasso(10, argv); - // CPPUNIT_ASSERT_EQUAL(dEploidIOlasso.obsWsaf_.size(), (size_t)594); - // CPPUNIT_ASSERT_DOUBLES_EQUAL(dEploidIOlasso.obsWsaf_[0], - // 0, epsilon3); - // CPPUNIT_ASSERT_DOUBLES_EQUAL(dEploidIOlasso.obsWsaf_[4], - // 0.22772277, epsilon2); - // CPPUNIT_ASSERT_DOUBLES_EQUAL(dEploidIOlasso.obsWsaf_[11], - // 0.21111111, epsilon2); - // CPPUNIT_ASSERT_DOUBLES_EQUAL(dEploidIOlasso.obsWsaf_[589], - // 0.18324607, epsilon2); + CPPUNIT_ASSERT_NO_THROW(dEploidIOlasso.dEploidLasso()); + CPPUNIT_ASSERT_EQUAL(dEploidIOlasso.lassoPanels.size(), (size_t)14); + CPPUNIT_ASSERT_EQUAL(dEploidIOlasso.lassoPlafs.size(), (size_t)14); + CPPUNIT_ASSERT_EQUAL(dEploidIOlasso.lassoRefCount.size(), (size_t)14); + CPPUNIT_ASSERT_EQUAL(dEploidIOlasso.lassoAltCount.size(), (size_t)14); } };