include panel header

DEploid-dev · Oct 10, 2018 · 68ceeb1 · 68ceeb1
1 parent d5c2657
commit 68ceeb1
Show file tree

Hide file tree

Showing 5 changed files with 53 additions and 16 deletions.
diff --git 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" <<endl;
     out << "./dEploid -vcf data/testData/PG0390-C.test.vcf -plaf data/testData/labStrains.test.PLAF.txt -o PG0390-CNopanel -ibdPainting -initialP 0.2 0.8" <<endl;
-    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-CLassoExclude -panel data/testData/labStrains.test.panel.txt -lasso -initialP 0.2 0.8" << endl;
+    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-CLassoExclude -panel data/testData/labStrains.test.panel.txt -lasso -initialP 0.2 0.8 -writePanel" << endl;
 }
 
 
@@ -878,7 +881,10 @@ void DEploidIO::dEploidLasso() {
         vector <double> wsaf = lassoComputeObsWsaf(start, length);
         vector < vector <double> > tmpPanel = lassoSubsetPanel(start, length);
         DEploidLASSO dummy(tmpPanel, wsaf, 250);
-
+        vector <string> 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
@@ -102,6 +102,7 @@ class DEploidIO{
     vector < vector <double> > lassoAltCount;
 
     void writeHap (vector < vector <double> > &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 <double> lassoComputeObsWsaf(size_t segmentStartIndex, size_t nLoci);
     vector < vector <double> > lassoSubsetPanel(size_t segmentStartIndex, size_t nLoci);
-    void writePanel(Panel *panel, size_t chromi);
+    void writePanel(Panel *panel, size_t chromi, vector <string> hdr);
 };
 
 #endif
diff --git a/src/export/writeMcmcRelated.cpp b/src/export/writeMcmcRelated.cpp
@@ -114,29 +114,28 @@ void DEploidIO::writeHap(vector < vector <double> > &hap, bool useIBD){
 }
 
 
-void DEploidIO::writePanel(Panel *panel, size_t chromI){
+void DEploidIO::writePanel(Panel *panel, size_t chromI, vector <string> 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
@@ -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
@@ -53,6 +53,7 @@ class TxtReader : public VariantIndex {
     // info_ only refers to the first column of the content
     vector <double> info_;
 
+    vector <string> 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;