Merge pull request #212 from h-2/bits2

Bits2

Author: sarahet

src/search_algo.hpp

@@ -1736,7 +1736,7 @@ computeBlastMatch(typename TBlastRecord::TBlastMatch  & bm,
     computeBitScore(bm, context(lH.gH.outfile));
 
     computeEValueThreadSafe(bm, record.qLength, context(lH.gH.outfile));
-    if (bm.eValue > lH.options.eCutOff)
+    if (bm.eValue > lH.options.maxEValue)
         return EVALUE;
 
     _setFrames(bm, m, lH);
@@ -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.outfile));
+
+            if (bm.bitScore < lH.options.minBitScore)
+            {
+                ++lH.stats.hitsFailedExtendBitScoreTest;
+                it = blastMatches.erase(it);
+                continue;
+            }
+        }
 
-        if (bm.eValue > lH.options.eCutOff)
+        if (lH.options.maxEValue < 100)
         {
-            ++lH.stats.hitsFailedExtendEValueTest;
-            it = blastMatches.erase(it);
-            continue;
+            computeEValueThreadSafe(bm, bm.qLength, seqan::context(lH.gH.outfile));
+
+            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.outfile));
 
-        // evalue computed previously
+        if (lH.options.maxEValue == 100)
+            computeEValueThreadSafe(bm, bm.qLength, seqan::context(lH.gH.outfile));
 
         ++it;
     }
@@ -2597,9 +2611,16 @@ iterateMatchesFullSerial(TLocalHolder & lH)
                                                            -band,
                                                            +band);
 
-        computeEValueThreadSafe(bm, record.qLength, context(lH.gH.outfile));
+        computeBitScore(bm, context(lH.gH.outfile));
+        if (bm.bitScore < lH.options.minBitScore)
+        {
+            ++lH.stats.hitsFailedExtendBitScoreTest;
+            record.matches.pop_back();
+            continue;
+        }
 
-        if (bm.eValue > lH.options.eCutOff)
+        computeEValueThreadSafe(bm, record.qLength, context(lH.gH.outfile));
+        if (bm.eValue > lH.options.maxEValue)
         {
             ++lH.stats.hitsFailedExtendEValueTest;
             record.matches.pop_back();
@@ -2624,7 +2645,6 @@ iterateMatchesFullSerial(TLocalHolder & lH)
             continue;
         }
 
-        computeBitScore(bm, context(lH.gH.outfile));
 
         if (lH.options.hasSTaxIds)
             bm.sTaxIds = lH.gH.sTaxIds[bm._n_sId];

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

@@ -96,7 +96,8 @@ struct LambdaOptions : public SharedOptions
 
     int             xDropOff    = 0;
     int             band        = -1;
-    double          eCutOff     = 0;
+    double          minBitScore = 0;
+    double          maxEValue   = 1e-04;
     int             idCutOff    = 0;
     unsigned long   maxMatches  = 500;
 
@@ -247,6 +248,14 @@ parseCommandLine(LambdaOptions & options, int argc, char const ** argv)
     setMinValue(parser, "e-value", "0");
     setMaxValue(parser, "e-value", "100");
 
+    addOption(parser, ArgParseOption("", "bit-score",
+        "Output only matches that score above this threshold.",
+        ArgParseArgument::DOUBLE));
+    setDefaultValue(parser, "bit-score", "0");
+    setMinValue(parser, "bit-score", "0");
+    setMaxValue(parser, "bit-score", "1000");
+
+
     addOption(parser, ArgParseOption("n", "num-matches",
         "Print at most this number of matches per query.",
         ArgParseArgument::INTEGER));
@@ -771,7 +780,8 @@ parseCommandLine(LambdaOptions & options, int argc, char const ** argv)
 
     getOptionValue(options.seedDeltaIncreasesLength, parser, "seed-delta-increases-length");
 
-    getOptionValue(options.eCutOff, parser, "e-value");
+    getOptionValue(options.maxEValue, parser, "e-value");
+    getOptionValue(options.minBitScore, parser, "bit-score");
     getOptionValue(options.idCutOff, parser, "percent-identity");
 
     getOptionValue(options.xDropOff, parser, "x-drop");
@@ -903,8 +913,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"