Merge branch 'develop'

Author: Rob Patro

include/RapMapUtils.hpp

@@ -36,6 +36,7 @@
 #include "chobo/small_vector.hpp"
 #include "RapMapConfig.hpp"
 #include "nonstd/optional.hpp"
+#include "FastxParser.hpp"
 
 #ifdef RAPMAP_SALMON_SUPPORT
 #include "LibraryFormat.hpp"
@@ -821,6 +822,10 @@ namespace rapmap {
 
         std::string reverseComplement(std::string& seq);
 
+
+        uint32_t writeUnalignedSingleToStream(fastx_parser::ReadSeq& r, fmt::MemoryWriter& sstream);
+        uint32_t writeUnalignedPairToStream(fastx_parser::ReadPair& r, fmt::MemoryWriter& sstream);
+
         template <typename ReadPairT, typename IndexT>
         uint32_t writeAlignmentsToStream(
                 ReadPairT& r,
@@ -1262,7 +1267,7 @@ namespace rapmap {
     template <typename Archive>
     void save(Archive& archive, const my_mer& mer);
 
-    template <typename Archive>
+    teplate <typename Archive>
     void load(Archive& archive, my_mer& mer);
     */
     }

src/FastxParser.cpp

@@ -158,7 +158,7 @@ int parseReads(
     while (!seqContainerQueue_.try_dequeue(*cCont, local)) {
       fastx_parser::thread_utils::backoffOrYield(curMaxDelay);
       // Think of a way to do this that wouldn't be loud (or would allow a user-definable logging mechanism)
-      // std::cerr << "couldn't dequeue read chunk\n";
+      //std::cerr << "couldn't dequeue read chunk\n";
     }
     size_t numObtained{local->size()};
     // open the file and init the parser
@@ -196,7 +196,10 @@ int parseReads(
     if (ksv == -3) {
       --numParsing;
       return -3;
-    }
+    } else if (ksv < -1) {
+      --numParsing;
+      return ksv;
+     }
 
     // If we hit the end of the file and have any reads in our local buffer
     // then dump them here.
@@ -207,7 +210,13 @@ int parseReads(
         fastx_parser::thread_utils::backoffOrYield(curMaxDelay);
       }
       numWaiting = 0;
+    } else if (numObtained > 0){
+      curMaxDelay = MIN_BACKOFF_ITERS;
+      while (!seqContainerQueue_.try_enqueue(std::move(local))) {
+        fastx_parser::thread_utils::backoffOrYield(curMaxDelay);
+      }
     }
+
     // destroy the parser and close the file
     kseq_destroy(seq);
     gzclose(fp);
@@ -300,7 +309,13 @@ int parseReadPair(
         fastx_parser::thread_utils::backoffOrYield(curMaxDelay);
       }
       numWaiting = 0;
+    } else if (numObtained > 0){
+      curMaxDelay = MIN_BACKOFF_ITERS;
+      while (!seqContainerQueue_.try_enqueue(std::move(local))) {
+        fastx_parser::thread_utils::backoffOrYield(curMaxDelay);
+      }
     }
+
     // destroy the parser and close the file
     kseq_destroy(seq);
     gzclose(fp);

src/RapMapSAMapper.cpp

@@ -148,6 +148,7 @@ struct MappingOpts {
   bool noOrphans{false};
   bool noDovetail{false};
   bool recoverOrphans{false};
+  bool writeUnmapped{false};
 };
 
 using AlnCacheMap = selective_alignment::utils::AlnCacheMap;
@@ -315,15 +316,13 @@ void processReadsSingleSA(single_parser * parser,
               }
             }
 
-            if (hits.size() > 0 and !mopts->noOutput) {
-                /*
-                std::sort(hits.begin(), hits.end(),
-                            [](const QuasiAlignment& a, const QuasiAlignment& b) -> bool {
-                                return a.tid < b.tid;
-                            });
-                */
+            if (!mopts->noOutput){
+              if (hits.size() > 0) {
                 rapmap::utils::writeAlignmentsToStream(read, formatter,
                                                        hctr, hits, sstream);
+              } else {
+                rapmap::utils::writeUnalignedSingleToStream(read, sstream);
+              }
             }
 
             if (hctr.numReads > hctr.lastPrint + 1000000) {
@@ -702,9 +701,13 @@ void processReadsPairSA(paired_parser* parser,
             hctr.totHits += jointHits.size();
 
             // If we have reads to output, and we're writing output.
-            if (jointHits.size() > 0 and !mopts->noOutput and jointHits.size() <= mopts->maxNumHits) {
+            if (!mopts->noOutput) {
+              if (jointHits.size() > 0 and jointHits.size() <= mopts->maxNumHits) {
                 rapmap::utils::writeAlignmentsToStream(rpair, formatter,
                                                        hctr, jointHits, sstream);
+              } else {
+                rapmap::utils::writeUnalignedPairToStream(rpair, sstream);
+              }
             }
 
             if (hctr.numReads > hctr.lastPrint + 1000000) {
@@ -1001,10 +1004,11 @@ int rapMapSAMap(int argc, char* argv[]) {
   TCLAP::SwitchArg chain("c", "chaining", "Score the hits to find the best chain", false);
   TCLAP::SwitchArg compressedOutput("x", "compressed", "Compress the output SAM file using zlib", false);
   TCLAP::SwitchArg quiet("q", "quiet", "Disable all console output apart from warnings and errors", false);
+  TCLAP::SwitchArg writeUnmapped("u", "writeUnmapped", "include unmapped reads in the output SAM records", false);
+
   TCLAP::SwitchArg recoverOrphans("", "recoverOrphans", "perform orphan recovery to try and recover endpoints of orphaned reads", false);
   TCLAP::SwitchArg noDovetail("", "noDovetail", "discard dovetailing mappings", false);
   TCLAP::SwitchArg noOrphans("", "noOrphans", "discard orphaned mappings", false);
-
   TCLAP::SwitchArg selAln("s", "selAln", "Perform selective alignment to validate mapping hits", false);
   TCLAP::ValueArg<int16_t> gapOpenPen("", "go", "[only with selAln]: gap open penalty", false, 4, "positive integer");
   TCLAP::ValueArg<int16_t> gapExtendPen("", "ge", "[only with selAln]: gap extend penalty", false, 2, "positive integer");
@@ -1026,6 +1030,7 @@ int rapMapSAMap(int argc, char* argv[]) {
   cmd.add(fuzzy);
   cmd.add(chain);
   cmd.add(quiet);
+  cmd.add(writeUnmapped);
   cmd.add(hardFilter);
   cmd.add(recoverOrphans);
   cmd.add(noDovetail);
@@ -1099,6 +1104,7 @@ int rapMapSAMap(int argc, char* argv[]) {
     } else {
         mopts.unmatedReads = unmatedReads.getValue();
     }
+    mopts.writeUnmapped = writeUnmapped.getValue();
     mopts.numThreads = numThreads.getValue();
     mopts.maxNumHits = maxNumHits.getValue();
     mopts.outname = (outname.isSet()) ? outname.getValue() : "";

src/RapMapUtils.cpp

@@ -133,6 +133,99 @@ namespace rapmap {
             return work;
         }
 
+
+      uint32_t writeUnalignedPairToStream(fastx_parser::ReadPair& r,
+                                          fmt::MemoryWriter& sstream) {
+        constexpr uint16_t flags1 = 0x1 | 0x4 | 0x8 | 0x40;
+        constexpr uint16_t flags2 = 0x1 | 0x4 | 0x8 | 0x80;
+
+        auto processReadName = [](const std::string& name) -> fmt::StringRef {
+                                 nonstd::string_view readNameView(name);
+                                 // If the read name contains multiple space-separated parts,
+                                 // print only the first
+                                 size_t splitPos = readNameView.find(' ');
+                                 if (splitPos < readNameView.length()) {
+                                   readNameView.remove_suffix(readNameView.length() - splitPos);
+                                 } else {
+                                   splitPos = readNameView.length();
+                                 }
+
+                                 // trim /1 from the pe read
+                                 if (splitPos > 2 and readNameView[splitPos - 2] == '/') {
+                                   readNameView.remove_suffix(2);
+                                   //readName[splitPos - 2] = '\0';
+                                 }
+                                 return fmt::StringRef(readNameView.data(), readNameView.size());
+                               };
+
+        auto readNameView = processReadName(r.first.name);
+        auto mateNameView = processReadName(r.second.name);
+        std::string* readSeq1 = &(r.first.seq);
+        std::string* readSeq2 = &(r.second.seq);
+
+        sstream << readNameView << '\t' // QNAME
+                << flags1 << '\t'       // FLAGS
+                << "*\t"                // RNAME
+                << "0\t"                // POS (1-based)
+                << "255\t"              // MAPQ
+                << "*\t"                // CIGAR
+                << "*\t"                // RNEXT
+                << "*\t"                // PNEXT
+                << "0\t"                // TLEN
+                << *readSeq1 << '\t'    // SEQ
+                << "*\t"                // QUAL
+                << "NH:i:0\t"
+                << "HI:i:0\t"
+                << "AS:i:0\n";
+
+        sstream << mateNameView << '\t' // QNAME
+                << flags2 << '\t'       // FLAGS
+                << "*\t"                // RNAME
+                << "0\t"                // POS (1-based)
+                << "255\t"              // MAPQ
+                << "*\t"                // CIGAR
+                << "*\t"                // RNEXT
+                << "*\t"                // PNEXT
+                << "0\t"                // TLEN
+                << *readSeq2 << '\t'    // SEQ
+                << "*\t"                // QUAL
+                << "NH:i:0\t"
+                << "HI:i:0\t"
+                << "AS:i:0\n";
+        return 0;
+      }
+
+      uint32_t writeUnalignedSingleToStream(fastx_parser::ReadSeq& r,
+                                            fmt::MemoryWriter& sstream) {
+        constexpr uint16_t flags = 0x4;
+
+        nonstd::string_view readNameViewSV(r.name);
+        // If the read name contains multiple space-separated parts, print
+        // only the first
+        size_t splitPos = readNameViewSV.find(' ');
+        if (splitPos < readNameViewSV.length()) {
+          readNameViewSV.remove_suffix(readNameViewSV.size() - splitPos);
+        }
+        fmt::StringRef readNameView(readNameViewSV.data(), readNameViewSV.size());
+        std::string* readSeq = &(r.seq);
+
+        sstream << readNameView << '\t' // QNAME
+                << flags << '\t'    // FLAGS
+                << "*\t"            // RNAME
+                << 0 << '\t'        // POS (1-based)
+                << 255 << '\t'      // MAPQ
+                << "*\t"            // CIGAR
+                << "*\t"            // MATE NAME
+                << "0\t"            // MATE POS
+                << "0\t"            // TLEN
+                << *readSeq << '\t' // SEQ
+                << "*\t"            // QSTR
+                << "NH:i:0\t"
+                << "HI:i:0\t"
+                << "AS:i:0\n";
+        return 0;
+      }
+
         template <typename ReadT, typename IndexT>
         uint32_t writeAlignmentsToStream(
                 ReadT& r,