resolve ksw2 int8_t overflow with mimicStrictBT2

COMBINE-lab · Dec 6, 2018 · b189626 · b189626
1 parent 3d16f2b
commit b189626
Show file tree

Hide file tree

Showing 7 changed files with 40 additions and 37 deletions.
diff --git a/src/Alevin.cpp b/src/Alevin.cpp
@@ -926,7 +926,7 @@ void initiatePipeline(AlevinOpts<ProtocolT>& aopt,
   }
 }
 
-int salmonBarcoding(int argc, char* argv[]) {
+int salmonBarcoding(int argc, const char* argv[]) {
   namespace bfs = boost::filesystem;
   namespace po = boost::program_options;
 

diff --git a/src/BuildSalmonIndex.cpp b/src/BuildSalmonIndex.cpp
@@ -43,7 +43,7 @@
 // http://stackoverflow.com/questions/108318/whats-the-simplest-way-to-test-whether-a-number-is-a-power-of-2-in-c
 bool isPowerOfTwo(uint32_t n) { return (n > 0 and (n & (n - 1)) == 0); }
 
-int salmonIndex(int argc, char* argv[]) {
+int salmonIndex(int argc, const char* argv[]) {
 
   using std::string;
   namespace bfs = boost::filesystem;

diff --git a/src/Salmon.cpp b/src/Salmon.cpp
@@ -93,7 +93,7 @@ int dualModeMessage() {
 /**
  * Bonus!
  */
-int salmonSwim(int argc, char* argv[]) {
+int salmonSwim(int argc, const char* argv[]) {
 
   std::cout << R"(
     _____       __
@@ -144,11 +144,11 @@ Nature Methods. 2017;14(4):417-419. doi: 10.1038/nmeth.4197
 )";
 }
 
-int salmonIndex(int argc, char* argv[]);
-int salmonQuantify(int argc, char* argv[]);
-int salmonAlignmentQuantify(int argc, char* argv[]);
-int salmonBarcoding(int argc, char* argv[]);
-int salmonQuantMerge(int argc, char* argv[]);
+int salmonIndex(int argc, const char* argv[]);
+int salmonQuantify(int argc, const char* argv[]);
+int salmonAlignmentQuantify(int argc, const char* argv[]);
+int salmonBarcoding(int argc, const char* argv[]);
+int salmonQuantMerge(int argc, const char* argv[]);
 
 bool verbose = false;
 
@@ -225,7 +225,7 @@ int main(int argc, char* argv[]) {
       opts.insert(opts.begin(), "--help");
     }
 
-    std::unordered_map<string, std::function<int(int, char* [])>> cmds(
+    std::unordered_map<string, std::function<int(int, const char* [])>> cmds(
         {{"index", salmonIndex},
          {"quant", salmonQuantify},
          {"quantmerge", salmonQuantMerge},
@@ -241,10 +241,10 @@ int main(int argc, char* argv[]) {
     */
 
     int32_t subCommandArgc = opts.size() + 1;
-    std::unique_ptr<char* []> argv2(new char*[subCommandArgc]);
+    std::unique_ptr<const char* []> argv2(new const char*[subCommandArgc]);
     argv2[0] = argv[0];
     for (int32_t i = 0; i < subCommandArgc - 1; ++i) {
-      argv2[i + 1] = &*opts[i].begin();
+      argv2[i + 1] = opts[i].c_str();
     }
 
     auto cmdMain = cmds.find(cmd);
@@ -262,11 +262,11 @@ int main(int argc, char* argv[]) {
         // detect mode-specific help request
         if (strncmp(argv2[1], "--help-alignment", 16) == 0) {
           std::vector<char> helpStr{'-', '-', 'h', 'e', 'l', 'p', '\0'};
-          char* helpArgv[] = {argv[0], &helpStr[0]};
+          const char* helpArgv[] = {argv[0], &helpStr[0]};
           return salmonAlignmentQuantify(2, helpArgv);
         } else if (strncmp(argv2[1], "--help-reads", 12) == 0) {
           std::vector<char> helpStr{'-', '-', 'h', 'e', 'l', 'p', '\0'};
-          char* helpArgv[] = {argv[0], &helpStr[0]};
+          const char* helpArgv[] = {argv[0], &helpStr[0]};
           return salmonQuantify(2, helpArgv);
         }
 

diff --git a/src/SalmonQuantMerge.cpp b/src/SalmonQuantMerge.cpp
@@ -208,7 +208,7 @@ bool doMerge(QuantMergeOptions& qmOpts) {
   return true;
 }
 
-int salmonQuantMerge(int argc, char* argv[]) {
+int salmonQuantMerge(int argc, const char* argv[]) {
   using std::cerr;
   using std::vector;
   using std::string;

diff --git a/src/SalmonQuantify.cpp b/src/SalmonQuantify.cpp
@@ -801,7 +801,7 @@ inline int32_t getAlnScore(
   if (invalidStart) { rptr += -pos; rlen += pos; pos = 0; }
 
   // if we are trying to mimic Bowtie2 with RSEM params
-  if (invalidStart and mimicStrictBT2) { return 0; }
+  if (invalidStart and mimicStrictBT2) { return s; }
 
   if (pos < tlen) {
     bool doUngapped{(!invalidStart) and (chainStat == rapmap::utils::ChainStatus::UNGAPPED)};
@@ -1194,14 +1194,14 @@ void processReadsQuasi(
               // throw away dovetailed reads
               if (mimicStrictBT2) {
                 if (h.fwd == h.mateIsFwd) {
-                  s1 = 0;
-                  s2 = 0;
+                  s1 = std::numeric_limits<int32_t>::min();
+                  s2 = std::numeric_limits<int32_t>::min();
                 } else if (h.fwd and (h.pos > h.matePos)) {
-                  s1 = 0;
-                  s2 = 0;
+                  s1 = std::numeric_limits<int32_t>::min();
+                  s2 = std::numeric_limits<int32_t>::min();
                 } else if (h.mateIsFwd and (h.matePos > h.pos)) {
-                  s1 = 0;
-                  s2 = 0;
+                  s1 = std::numeric_limits<int32_t>::min();
+                  s2 = std::numeric_limits<int32_t>::min();
                 }
               }
 
@@ -1932,13 +1932,12 @@ void processReadLibrary(
    **/
   // NOTE : When we can support C++14, we can replace the entire ProcessFunctor class above with this
   // generic lambda.
-  auto processFunctor = [&](size_t i, auto& parserPtr, auto* index) {
+  auto processFunctor = [&](size_t i, auto* parserPtr, auto* index) {
     if (salmonOpts.qmFileName != "" and i == 0) {
       rapmap::utils::writeSAMHeader(*index, salmonOpts.qmLog);
     }
-    // TODO: understand why the index must be captured by value on old OSX compilers
-    auto threadFun = [&, i, index]() -> void {
-      processReadsQuasi(parserPtr.get(), readExp, rl, structureVec[i],
+    auto threadFun = [&, i, parserPtr, index]() -> void {
+      processReadsQuasi(parserPtr, readExp, rl, structureVec[i],
                         numObservedFragments, numAssignedFragments, numValidHits,
                         upperBoundHits, index, transcripts,
                         fmCalc, clusterForest, fragLengthDist, observedBiasParams[i],
@@ -1997,19 +1996,19 @@ void processReadLibrary(
         // change value before the lambda below is evaluated --- crazy!
         if (largeIndex) {
           if (perfectHashIndex) { // Perfect Hash
-            if (isPairedEnd) {processFunctor(i, pairedParserPtr, sidx->quasiIndexPerfectHash64());}
-            else if (isSingleEnd) {processFunctor(i, singleParserPtr, sidx->quasiIndexPerfectHash64());}
+            if (isPairedEnd) {processFunctor(i, pairedParserPtr.get(), sidx->quasiIndexPerfectHash64());}
+            else if (isSingleEnd) {processFunctor(i, singleParserPtr.get(), sidx->quasiIndexPerfectHash64());}
           } else { // Dense Hash
-            if (isPairedEnd) {processFunctor(i, pairedParserPtr, sidx->quasiIndex64());}
-            else if (isSingleEnd) {processFunctor(i, singleParserPtr, sidx->quasiIndex64());}
+            if (isPairedEnd) {processFunctor(i, pairedParserPtr.get(), sidx->quasiIndex64());}
+            else if (isSingleEnd) {processFunctor(i, singleParserPtr.get(), sidx->quasiIndex64());}
           }
         } else {
           if (perfectHashIndex) { // Perfect Hash
-            if (isPairedEnd) { processFunctor(i, pairedParserPtr, sidx->quasiIndexPerfectHash32()); }
-            else if (isSingleEnd) { processFunctor(i, singleParserPtr, sidx->quasiIndexPerfectHash32()); }
+            if (isPairedEnd) { processFunctor(i, pairedParserPtr.get(), sidx->quasiIndexPerfectHash32()); }
+            else if (isSingleEnd) { processFunctor(i, singleParserPtr.get(), sidx->quasiIndexPerfectHash32()); }
           } else { // Dense Hash
-            if (isPairedEnd) { processFunctor(i, pairedParserPtr, sidx->quasiIndex32()); }
-            else if (isSingleEnd) { processFunctor(i, singleParserPtr, sidx->quasiIndex32()); }
+            if (isPairedEnd) { processFunctor(i, pairedParserPtr.get(), sidx->quasiIndex32()); }
+            else if (isSingleEnd) { processFunctor(i, singleParserPtr.get(), sidx->quasiIndex32()); }
           }
         } // End spawn current thread
 
@@ -2332,7 +2331,7 @@ void quantifyLibrary(ReadExperimentT& experiment, bool greedyChain,
   jointLog->info("finished quantifyLibrary()");
 }
 
-int salmonQuantify(int argc, char* argv[]) {
+int salmonQuantify(int argc, const char* argv[]) {
   using std::cerr;
   using std::vector;
   using std::string;

diff --git a/src/SalmonQuantifyAlignments.cpp b/src/SalmonQuantifyAlignments.cpp
@@ -1373,7 +1373,7 @@ bool processSample(AlignmentLibraryT<ReadT>& alnLib, size_t requiredObservations
   return true;
 }
 
-int salmonAlignmentQuantify(int argc, char* argv[]) {
+int salmonAlignmentQuantify(int argc, const char* argv[]) {
   using std::cerr;
   using std::vector;
   using std::string;

diff --git a/src/SalmonUtils.cpp b/src/SalmonUtils.cpp
@@ -1477,8 +1477,12 @@ std::string getCurrentTimeAsString() {
       sopt.minScoreFraction = 0.8;
       sopt.matchScore = 1;
       sopt.mismatchPenalty = 0;
-      sopt.gapOpenPenalty = 50;
-      sopt.gapExtendPenalty = 50;
+      // NOTE: as a limitation of ksw2, we can't have
+      // (gapOpenPenalty + gapExtendPenalty) * 2 + matchScore < numeric_limits<int8_t>::max()
+      // these parameters below are sufficiently large penalties to
+      // prohibit gaps, while not overflowing the above condition
+      sopt.gapOpenPenalty = 25;
+      sopt.gapExtendPenalty = 25;
     }
 
     // If the consensus slack was not set explicitly, then it defaults to 1 with