adding functionality to extract consensus sequence for blocks

Author: sumit-walia

src/consensus.cpp

@@ -0,0 +1,34 @@
+#include "panmanUtils.hpp"
+
+void panmanUtils::Tree::printConsensus(std::ostream& fout){
+    std::mutex printMutex;
+
+    tbb::parallel_for_each(blocks.begin(), blocks.end(), [&](const auto& eachBlock) {
+    // for(size_t i = 0; i < blocks.size(); i++) {
+
+        std::string consensus;
+
+        int32_t primaryBlockId = ((int32_t)eachBlock.primaryBlockId);
+        int32_t secondaryBlockId = ((int32_t)eachBlock.secondaryBlockId);
+
+        for(size_t j = 0; j < eachBlock.consensusSeq.size(); j++) {
+            bool endFlag = false;
+            for(size_t k = 0; k < 8; k++) {
+                const int nucCode = (((eachBlock.consensusSeq[j]) >> (4*(7 - k))) & 15);
+
+                if(nucCode == 0) {
+                    endFlag = true;
+                    break;
+                }
+                const char nucleotide = panmanUtils::getNucleotideFromCode(nucCode);
+                consensus.push_back(nucleotide);
+            }
+            if(endFlag) {
+                break;
+            }
+        }
+        std::lock_guard<std::mutex> guard(printMutex);
+        fout << ">" << primaryBlockId << "\n" << consensus << "\n";
+    });
+
+}

src/panman.cpp

@@ -20,7 +20,9 @@
 #include <set>
 #include <boost/iostreams/filter/lzma.hpp>
 #include <typeinfo>
-
+#include <capnp/serialize.h>
+#include <kj/io.h>
+#include <limits>
 
 #include "chaining.cpp"
 #include "rotation.cpp"
@@ -35,6 +37,7 @@
 #include "reroot.cpp"
 #include "aaTrans.cpp"
 #include "panman2usher.cpp"
+#include "consensus.cpp"
 #include "panmanUtils.hpp"
 
 char panmanUtils::getNucleotideFromCode(int code) {
@@ -6950,9 +6953,9 @@ panmanUtils::TreeGroup::TreeGroup(std::istream& fin, bool isOld) {
 	capnp::ReaderOptions readerOptions;
         readerOptions.traversalLimitInWords = std::numeric_limits<uint64_t>::max();
         readerOptions.nestingLimit = 1024;
-        capnp::InputStreamMessageReader messageReader(kjInputStream, readerOptions);
+	capnp::InputStreamMessageReader messageReader(kjInputStream, readerOptions);
 
-        panman::TreeGroup::Reader TG = messageReader.getRoot<panman::TreeGroup>();
+    	panman::TreeGroup::Reader TG = messageReader.getRoot<panman::TreeGroup>();
 
         int count=0;
         for (auto treeFromTG: TG.getTrees()){

src/panman.hpp

@@ -545,7 +545,7 @@ class Tree {
     void printVCFParallel(std::string reference, std::ostream& fout);
     void printVCFParallel(panmanUtils::Node* node, std::string& fileName);
     void extractAminoAcidTranslations(std::ostream& fout, int64_t start, int64_t end);
-
+    void printConsensus(std::ostream& fout);
     // Extract PanMAT representing a segment of the genome. The start and end coordinates
     // are with respect to the root sequence. The strands of the terminal blocks in all
     // sequences are assumed to be the same as their strands in the root sequence for the

src/panmanUtils.cpp

@@ -122,6 +122,7 @@ po::positional_options_description printPathsArgumentDesc;
 po::options_description indexDesc("Indexing Command Line Arguments");
 po::positional_options_description indexArgumentDesc;
 po::options_description printRootDesc("Root Printer Command Line Arguments");
+po::options_description printConsensusDesc("Consensus Printer Command Line Arguments");
 po::positional_options_description printRootPositionArgumentDesc;
 
 
@@ -158,7 +159,7 @@ void setupOptionDescriptions() {
     ("acr,q", "ACR method [fitch(default), mppa]")
     ("index",po::value< bool >(0), "Generating indexes and print sequence (passed as reference) between x:y")
     ("refFile",po::value< std::string >() ,"reference sequence file")
-    // ("printRoot", "Print root sequence")
+    ("printConsensus", "Print consensus sequences of each block in FASTA format")
     // ("printNodePaths", "Print mutations from root to each node")
     ("toUsher", "Convert a PanMAT in PanMAN to Usher-MAT")
     // ("samples",po::value< bool >() ,"Samples in the PanMAN")
@@ -260,6 +261,9 @@ void setupOptionDescriptions() {
 
     printRootDesc.add_options()
         ("output-file,o", po::value< std::string >(), "Output file name");
+
+    printConsensusDesc.add_options()
+        ("output-file,o", po::value< std::string >(), "Output file name");
 }
 
 void writePanMAN(po::variables_map &globalVm, panmanUtils::TreeGroup *TG) {
@@ -1208,6 +1212,39 @@ void printRoot(panmanUtils::TreeGroup *TG, po::variables_map &globalVm, std::ofs
 
 }
 
+void printConsensus(panmanUtils::TreeGroup *TG, po::variables_map &globalVm, std::ofstream &outputFile, std::streambuf * buf) {
+    // Print raw sequences to output file
+    if(TG == nullptr) {
+        std::cout << "No PanMAN selected" << std::endl;
+        return;
+    }
+
+    panmanUtils::TreeGroup tg = *TG;
+
+    auto fastaStart = std::chrono::high_resolution_clock::now();
+    for(int i = 0; i < tg.trees.size(); i++) {
+        panmanUtils::Tree * T = &tg.trees[i];
+        if(globalVm.count("output-file")) {
+            std::string fileName = globalVm["output-file"].as< std::string >();
+            outputFile.open("./info/" + fileName + "_" + std::to_string(i) + ".fasta");
+            buf = outputFile.rdbuf();
+        } else {
+            buf = std::cout.rdbuf();
+        }
+        std::ostream fout (buf);
+
+
+        T->printConsensus(fout);
+
+        if(globalVm.count("output-file")) outputFile.close();
+    }
+
+    auto fastaEnd = std::chrono::high_resolution_clock::now();
+    std::chrono::nanoseconds fastaTime = fastaEnd - fastaStart;
+    std::cout << "\nFASTA execution time: " << fastaTime.count() << " nanoseconds\n";
+
+}
+
 void toUsher(panmanUtils::TreeGroup *TG, po::variables_map &globalVm) {
     // Print raw sequences to output file
     if(TG == nullptr) {
@@ -1549,6 +1586,8 @@ void parseAndExecute(int argc, char* argv[]) {
     // } else if(globalVm.count("fasta-fast")){
     //     fastaFast(TG, globalVm, outputFile, buf);
     //     return;
+    } else if (globalVm.count("printConsensus")) {
+        printConsensus(TG, globalVm, outputFile, buf);
     } else {
         char** splitCommandArray;
 
@@ -1691,6 +1730,12 @@ void parseAndExecute(int argc, char* argv[]) {
                         .options(printRootDesc)
                         .run(), printRootVm);
                     printRoot(TG, printRootVm, outputFile, buf);
+                } else if(strcmp(splitCommandArray[0], "printConsensus") == 0) {
+                    po::variables_map printConsensusVm;
+                    po::store(po::command_line_parser((int)splitCommand.size(), splitCommandArray)
+                        .options(printConsensusDesc)
+                        .run(), printConsensusVm);
+                    printConsensus(TG, printConsensusVm, outputFile, buf);
                 } else if (strcmp(splitCommandArray[0], "exit") == 0 || strcmp(splitCommandArray[0], "q") == 0) {
                     return;
                 } else {