[feature] minimum bit-score

Author: h-2

src/search_algo.hpp

@@ -2431,17 +2431,28 @@ iterateMatchesFullSimd(TLocalHolder & lH)
     {
         TBlastMatch & bm = *it;
 
-        computeEValueThreadSafe(bm,
-                                qIsTranslated(TGlobalHolder::blastProgram)
-                                    ? lH.gH.untransQrySeqLengths[bm._n_qId]
-                                    : length(lH.gH.qrySeqs[_untrueQryId(bm, lH)]),
-                                context(lH.gH.outfile));
+        if (lH.options.minBitScore >= 0)
+        {
+            seqan::computeBitScore(bm, seqan::context(lH.gH.outfileBlastTab));
 
-        if (bm.eValue > lH.options.eCutOff)
+            if (bm.bitScore < lH.options.minBitScore)
+            {
+                ++lH.stats.hitsFailedExtendBitScoreTest;
+                it = blastMatches.erase(it);
+                continue;
+            }
+        }
+
+        if (lH.options.maxEValue >= 0)
         {
-            ++lH.stats.hitsFailedExtendEValueTest;
-            it = blastMatches.erase(it);
-            continue;
+            computeEValueThreadSafe(bm, bm.qLength, seqan::context(lH.gH.outfileBlastTab));
+
+            if (bm.eValue > lH.options.maxEValue)
+            {
+                ++lH.stats.hitsFailedExtendEValueTest;
+                it = blastMatches.erase(it);
+                continue;
+            }
         }
 
         ++it;
@@ -2484,9 +2495,12 @@ iterateMatchesFullSimd(TLocalHolder & lH)
             continue;
         }
 
-        computeBitScore(bm, context(lH.gH.outfile));
+        // not computed previously
+        if (lH.options.minBitScore < 0)
+            seqan::computeBitScore(bm, seqan::context(lH.gH.outfileBlastTab));
 
-        // evalue computed previously
+        if (lH.options.maxEValue < 0)
+            computeEValueThreadSafe(bm, bm.qLength, seqan::context(lH.gH.outfileBlastTab));
 
         ++it;
     }

src/search_datastructures.hpp

@@ -111,6 +111,7 @@ struct StatsHolder
 
 // post-extension
     uint64_t hitsFailedExtendPercentIdentTest;
+    uint64_t hitsFailedExtendBitScoreTest;
     uint64_t hitsFailedExtendEValueTest;
     uint64_t hitsAbundant;
     uint64_t hitsDuplicate;
@@ -150,6 +151,7 @@ struct StatsHolder
         hitsPutativeAbundant = 0;
 
         hitsFailedExtendPercentIdentTest = 0;
+        hitsFailedExtendBitScoreTest = 0;
         hitsFailedExtendEValueTest = 0;
         hitsAbundant = 0;
         hitsDuplicate = 0;
@@ -183,6 +185,7 @@ struct StatsHolder
         hitsPutativeAbundant += rhs.hitsPutativeAbundant;
 
         hitsFailedExtendPercentIdentTest += rhs.hitsFailedExtendPercentIdentTest;
+        hitsFailedExtendBitScoreTest += rhs.hitsFailedExtendBitScoreTest;
         hitsFailedExtendEValueTest += rhs.hitsFailedExtendEValueTest;
         hitsAbundant += rhs.hitsAbundant;
         hitsDuplicate += rhs.hitsDuplicate;
@@ -253,12 +256,15 @@ void printStats(StatsHolder const & stats, LambdaOptions const & options)
             std::cout << "\n - failed pre-extend test   " << R
                       << stats.hitsFailedPreExtendTest  << RR
                       << (rem -= stats.hitsFailedPreExtendTest);
-        std::cout << "\n - failed %-identity test   " << R
-                  << stats.hitsFailedExtendPercentIdentTest << RR
-                  << (rem -= stats.hitsFailedExtendPercentIdentTest);
         std::cout << "\n - failed e-value test      " << R
                   << stats.hitsFailedExtendEValueTest << RR
                   << (rem -= stats.hitsFailedExtendEValueTest);
+        std::cout << "\n - failed bitScore test     " << R
+                  << stats.hitsFailedExtendBitScoreTest << RR
+                  << (rem -= stats.hitsFailedExtendBitScoreTest);
+        std::cout << "\n - failed %-identity test   " << R
+                  << stats.hitsFailedExtendPercentIdentTest << RR
+                  << (rem -= stats.hitsFailedExtendPercentIdentTest);
         std::cout << "\n - duplicates               " << R
                   << stats.hitsDuplicate              << RR
                   << (rem -= stats.hitsDuplicate);

src/search_options.hpp

@@ -95,8 +95,8 @@ struct LambdaOptions : public SharedOptions
     int             misMatch        = 0; // only for manual
 
     int             xDropOff    = 0;
-    int             band        = -1;
-    double          eCutOff     = 0;
+    int32_t         minBitScore = -1;
+    double          maxEValue   = 1e-04;
     int             idCutOff    = 0;
     unsigned long   maxMatches  = 500;
 
@@ -903,8 +903,9 @@ printOptions(LambdaOptions const & options)
               << "  db index type:            " << _indexEnumToName(options.dbIndexType) << "\n"
               << " OUTPUT (file)\n"
               << "  output file:              " << options.output << "\n"
+              << "  maximum e-value:          " << options.maxEValue << "\n"
+              << "  minimum bit-score:        " << options.minBitScore << "\n"
               << "  minimum % identity:       " << options.idCutOff << "\n"
-              << "  maximum e-value:          " << options.eCutOff << "\n"
               << "  max #matches per query:   " << options.maxMatches << "\n"
               << "  include subj names in sam:" << options.samWithRefHeader << "\n"
               << "  include seq in sam/bam:   " << options.samBamSeq << "\n"