Implement Cigar::edit_distance() to count edits correctly in ksw_extend

ksahlin · Feb 17, 2023 · 37f8322 · 37f8322
1 parent d88ba30
commit 37f8322
Show file tree

Hide file tree

Showing 3 changed files with 20 additions and 3 deletions.
diff --git a/src/cigar.hpp b/src/cigar.hpp
@@ -33,6 +33,19 @@ class Cigar {
         }
     }
 
+    /* This works only if I, D, X, = are the only operations used */
+    int edit_distance() const {
+        auto dist = 0;
+        for (auto op_len : m_ops) {
+            auto op = op_len & 0xf;
+            auto len = op_len >> 4;
+            if (op == CIGAR_INS || op == CIGAR_DEL || op == CIGAR_X) {
+                dist += len;
+            }
+        }
+        return dist;
+    }
+
     /* 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
@@ -90,9 +90,9 @@ aln_info ksw_extend(const std::string& query, const std::string& ref, int8_t mat
     global_aln(query, ref, 5, mat, gap_open, gap_extend, w, zdrop, flag, &ez);
 
     aln_info info;
-    auto edits = std::count_if(ez.cigar, ez.cigar + ez.n_cigar, [](char c) { return c == 'X' || c == 'I' || c == 'D'; });
-    info.cigar = Cigar(ez.cigar, ez.n_cigar).to_eqx(query, ref).to_string();
-    info.edit_distance = edits;
+    auto cigar = Cigar(ez.cigar, ez.n_cigar).to_eqx(query, ref);
+    info.cigar = cigar.to_string();
+    info.edit_distance = cigar.edit_distance();
     info.ref_start = 0;
     info.ref_end = ez.max_t + 1;
     info.query_end = ez.max_q + 1;

diff --git a/tests/test_cigar.cpp b/tests/test_cigar.cpp
@@ -67,3 +67,7 @@ TEST_CASE("concatenate Cigar") {
     c += Cigar{"2M1X"};
     CHECK(c.to_string() == "5M1X");
 }
+
+TEST_CASE("edit distance") {
+    CHECK(Cigar("3=1X4D5I7=").edit_distance() == 10);
+}