Merge pull request #36 from ekg/secondary_stuff_one_to_one

keep secondary alignment also in one-to-one mode
waveygang · Apr 20, 2021 · 2314c8d · 2314c8d
2 parents b92f076 + 9287ac4
commit 2314c8d
Show file tree

Hide file tree

Showing 2 changed files with 11 additions and 5 deletions.
diff --git a/src/map/include/computeMap.hpp b/src/map/include/computeMap.hpp
@@ -204,7 +204,10 @@ namespace skch
         //Filter over reference axis and report the mappings
         if (param.filterMode == filter::ONETOONE)
         {
-          skch::Filter::ref::filterMappings(allReadMappings, this->refSketch);
+          skch::Filter::ref::filterMappings(allReadMappings, this->refSketch,
+                                            param.secondaryToKeep
+                                           // (input->len < param.segLength ? param.shortSecondaryToKeep : param.secondaryToKeep)
+                                            );
 
           //Re-sort mappings by input order of query sequences
           //This order may be needed for any post analysis of output

diff --git a/src/map/include/filter.hpp b/src/map/include/filter.hpp
@@ -284,14 +284,17 @@ namespace skch
          * @param[in/out]   L             container with mappings
          */
         template <typename Type>
-          inline void markGood(Type &L)
+          inline void markGood(Type &L, int secondaryToKeep)
           {
             //first segment in the set order
             auto beg = L.begin();
 
+            // count how many secondary alignments we keep
+            int kept = 0;
+
             for(auto it = L.begin(); it != L.end(); it++)
             {
-              if(this->greater_score(*beg, *it) || vec[*it].discard == 0) 
+              if((this->greater_score(*beg, *it) || vec[*it].discard == 0) && ++kept > secondaryToKeep)
                 break;
 
               vec[*it].discard = 0;
@@ -325,7 +328,7 @@ namespace skch
        * @param[in]     refsketch     reference index class object, used to determine ref sequence lengths
        */
       template <typename VecIn>
-        void filterMappings(VecIn &readMappings, const skch::Sketch &refsketch)
+        void filterMappings(VecIn &readMappings, const skch::Sketch &refsketch, int secondaryToKeep)
         {
           if(readMappings.size() <= 1)
             return;
@@ -376,7 +379,7 @@ namespace skch
                                     });
 
             //mark mappings as good
-            obj.markGood(bst);
+            obj.markGood(bst, secondaryToKeep);
 
             it = it2;
           }