ksw’s built-in CIGAR-reversing capability does not seem to work

At least not the way I expected it, so we just reverse the CIGAR ourselves.
ksahlin · Feb 17, 2023 · f3c05b3 · f3c05b3
1 parent 37f8322
commit f3c05b3
Show file tree

Hide file tree

Showing 5 changed files with 17 additions and 9 deletions.
diff --git a/src/aln.cpp b/src/aln.cpp
@@ -278,8 +278,7 @@ static inline alignment get_alignment(
     std::reverse(left_ref.begin(), left_ref.end());
 
     auto left = ksw_extend(left_query, left_ref,
-        aligner.parameters.match, aligner.parameters.mismatch, aligner.parameters.gap_open - aligner.parameters.gap_extend, aligner.parameters.gap_extend,
-        true // reverse CIGAR
+        aligner.parameters.match, aligner.parameters.mismatch, aligner.parameters.gap_open - aligner.parameters.gap_extend, aligner.parameters.gap_extend
     );
 
     // build final CIGAR
@@ -288,7 +287,9 @@ static inline alignment get_alignment(
     if (left_soft_clip > 0) {
         cigar.push(CIGAR_SOFTCLIP, left_soft_clip);
     }
-    cigar += Cigar(left.cigar);
+    auto left_cigar = Cigar(left.cigar);
+    left_cigar.reverse();
+    cigar += left_cigar;
     cigar += Cigar(middle_ops);
     cigar += Cigar(right.cigar);
 

diff --git a/src/cigar.hpp b/src/cigar.hpp
@@ -4,6 +4,7 @@
 #include <cstdint>
 #include <string>
 #include <vector>
+#include <algorithm>
 
 enum CIGAR {
     CIGAR_MATCH = 0,
@@ -46,6 +47,10 @@ class Cigar {
         return dist;
     }
 
+    void reverse() {
+        std::reverse(m_ops.begin(), m_ops.end());
+    }
+
     /* Return a new Cigar that uses =/X instead of M */
     Cigar to_eqx(const std::string& query, const std::string& ref) const;
 

diff --git a/src/ksw2wrap.cpp b/src/ksw2wrap.cpp
@@ -71,14 +71,10 @@ void global_aln(
 
 }  // namespace
 
-
-aln_info ksw_extend(const std::string& query, const std::string& ref, int8_t match, int8_t mismatch, int8_t gap_open, int8_t gap_extend, bool reverse_cigar) {
+aln_info ksw_extend(const std::string& query, const std::string& ref, int8_t match, int8_t mismatch, int8_t gap_open, int8_t gap_extend) {
     int w = -1; // band width; -1 is inf
     int zdrop = -1; // -1 to disable
     int flag = KSW_EZ_EXTZ_ONLY;
-    if (reverse_cigar) {
-        flag |= KSW_EZ_REV_CIGAR;
-    }
 
     //int c, i, pair = 1;
     //char *algo = "extd", *s;

diff --git a/src/ksw2wrap.hpp b/src/ksw2wrap.hpp
@@ -5,6 +5,6 @@
 #include <cstdint>
 #include "aligner.hpp"
 
-aln_info ksw_extend(const std::string& query, const std::string& ref, int8_t match, int8_t mismatch, int8_t gap_open, int8_t gap_extend, bool reverse_cigar = false);
+aln_info ksw_extend(const std::string& query, const std::string& ref, int8_t match, int8_t mismatch, int8_t gap_open, int8_t gap_extend);
 
 #endif
diff --git a/tests/test_cigar.cpp b/tests/test_cigar.cpp
@@ -71,3 +71,9 @@ TEST_CASE("concatenate Cigar") {
 TEST_CASE("edit distance") {
     CHECK(Cigar("3=1X4D5I7=").edit_distance() == 10);
 }
+
+TEST_CASE("reversed") {
+    Cigar c{"3=1X4D5I7="};
+    c.reverse();
+    CHECK(c.to_string() == "7=5I4D1X3=");
+}