adding new method to write the full consensus sequence

cpp/Vihgo.cpp

@@ -31,7 +31,7 @@ Vihgo::Vihgo() {
 
 Vihgo::Vihgo(string iBF, string gF, string rtF, string intF, string protF,
              string refF, string outF, string cResF, string csF, string iRCB,
-             int dT, int nT, int lL, double rT, double pvT, bool countStop) {
+             int dT, int nT, int lL, double rT, double pvT, bool countStop, bool fullCS) {
   this->inputBamFile = iBF;
   this->geneFile = gF;
   this->integraseFile = intF;
@@ -48,6 +48,7 @@ Vihgo::Vihgo(string iBF, string gF, string rtF, string intF, string protF,
   this->ratioThreshold = rT;
   this->pvalueThreshold = pvT;
   this->countStop = countStop;
+  this->fullCS = fullCS;
 
   this->codonTable.insert(pair<string, string>("TTT", "F"));
   this->codonTable.insert(pair<string, string>("TTC", "F"));
@@ -130,7 +131,14 @@ int Vihgo::mainLoop() {
   this->displayData("control");
   this->displayData("sample");
   this->writeCodonFile("sample");
-  this->writeConsensus();
+
+  if (this->fullCS){
+    this->writeFullConsensus();
+  }
+  else{
+    this->writeConsensus();
+  }
+
   this->variantCalling();
 
   return 0;
@@ -1727,6 +1735,70 @@ void Vihgo::writeConsensus() {
   csStream.close();
 }
 
+// Method to compute and write the full consensus sequence
+void Vihgo::writeFullConsensus() {
+  cerr << "Writing full consensus sequence ... " << endl;
+  ofstream csStream;
+  csStream.open((this->getConsensusFile()).c_str(), ios::out);
+
+  map<string, map<int, map<string, int>>>::iterator myIterA;
+  map<int, map<string, int>>::iterator myIterB;
+  map<string, int>::iterator myIterC;
+
+  map<string, map<int, string>>::iterator myIterE;
+  map<int, string>::iterator myIterF;
+
+  for (myIterE = this->protCodon.begin(); myIterE != this->protCodon.end();
+       myIterE++) {
+    string prot = myIterE->first;
+    string chrom = getChrFromProt(prot);
+
+    if (!(this->sampleData.count(chrom) > 0)) {
+      continue;
+    }
+    csStream << ">" << prot << endl;
+
+    for (myIterF = (myIterE->second).begin();
+         myIterF != (myIterE->second).end(); myIterF++) {
+      int aaPosition = myIterF->first;
+      int gPosition;
+      map<int, int> tmpMap = this->refProtCorresp[prot];
+      map<int, int>::iterator myIterG;
+      for (myIterG = tmpMap.begin(); myIterG != tmpMap.end(); myIterG++) {
+        if (myIterG->second == aaPosition) {
+          gPosition = myIterG->first;
+          break;
+        }
+      }
+
+      // cerr << "ppos : " << aaPosition << " ; genomic position is : " <<
+      // gPosition << endl;
+
+      // determine most frequent codon
+      string codon = "NNN";
+      int count = 0;
+
+      if (this->sampleData[chrom].count(gPosition) > 0) {
+        map<string, int> tmpMapCodon;
+        tmpMapCodon = this->sampleData[chrom][gPosition];
+        map<string, int>::iterator myIterH;
+        for (myIterH = tmpMapCodon.begin(); myIterH != tmpMapCodon.end();
+             myIterH++) {
+          if (myIterH->second > count) {
+            codon = myIterH->first;
+            count = myIterH->second;
+          }
+        }
+      }
+      csStream << codon;
+    }
+    csStream << endl;
+  }
+  csStream.close();
+}
+
+
+
 // method to get the chromosome in which the prot is included
 string Vihgo::getChrFromProt(string incProt) {
   return this->genes[incProt]["chrom"];