Merge pull request #255 from waveygang/inversion-patch

Inversion patching
waveygang · Jul 1, 2024 · e3dba1b · e3dba1b
2 parents a8e012f + c3ffa66
commit e3dba1b
Show file tree

Hide file tree

Showing 6 changed files with 323 additions and 78 deletions.
diff --git a/src/common/wflign/src/wflign.cpp b/src/common/wflign/src/wflign.cpp
@@ -91,6 +91,7 @@ WFlign::WFlign(
     this->erode_k = erode_k;
     this->chain_gap = chain_gap;
     this->max_patching_score = max_patching_score;
+    this->min_inversion_length = 23;
     // Query
     this->query_name = nullptr;
     this->query = nullptr;
@@ -539,6 +540,7 @@ void WFlign::wflign_affine_wavefront(
                 erode_k,
                 chain_gap,
                 max_patching_score,
+                min_inversion_length,
                 MIN_WF_LENGTH,
                 wf_max_dist_threshold
 #ifdef WFA_PNG_TSV_TIMING
@@ -1007,6 +1009,7 @@ void WFlign::wflign_affine_wavefront(
                         erode_k,
                         chain_gap,
                         max_patching_score,
+                        min_inversion_length,
                         MIN_WF_LENGTH,
                         wf_max_dist_threshold
 #ifdef WFA_PNG_TSV_TIMING

diff --git a/src/common/wflign/src/wflign.hpp b/src/common/wflign/src/wflign.hpp
@@ -64,6 +64,7 @@ namespace wflign {
             int erode_k;
             int64_t chain_gap;
             int max_patching_score;
+            uint64_t min_inversion_length;
             // Query
             const std::string* query_name;
             char* query;

diff --git a/src/common/wflign/src/wflign_alignment.cpp b/src/common/wflign/src/wflign_alignment.cpp
@@ -11,18 +11,82 @@
 /*
  * Wflign Alignment
  */
-alignment_t::alignment_t() {
-    j = 0;
-    i = 0;
-    query_length = 0;
-    target_length = 0;
-    ok = false;
-    keep = false;
-    edit_cigar = {NULL, 0, 0};
+
+// Default constructor
+alignment_t::alignment_t()
+    : j(0), i(0), query_length(0), target_length(0), ok(false), keep(false) {
+    edit_cigar = {nullptr, 0, 0};
 }
+
+// Destructor
 alignment_t::~alignment_t() {
-    free(edit_cigar.cigar_ops);
+        free(edit_cigar.cigar_ops);
+    }
+
+// Copy constructor
+alignment_t::alignment_t(const alignment_t& other)
+    : j(other.j), i(other.i), query_length(other.query_length),
+      target_length(other.target_length), ok(other.ok), keep(other.keep) {
+    if (other.edit_cigar.cigar_ops) {
+        edit_cigar.cigar_ops = (char*)malloc((other.edit_cigar.end_offset - other.edit_cigar.begin_offset) * sizeof(char));
+        memcpy(edit_cigar.cigar_ops, other.edit_cigar.cigar_ops + other.edit_cigar.begin_offset, 
+               (other.edit_cigar.end_offset - other.edit_cigar.begin_offset) * sizeof(char));
+    } else {
+        edit_cigar.cigar_ops = nullptr;
+    }
+    edit_cigar.begin_offset = 0;
+    edit_cigar.end_offset = other.edit_cigar.end_offset - other.edit_cigar.begin_offset;
+}
+
+// Move constructor
+alignment_t::alignment_t(alignment_t&& other) noexcept
+    : j(other.j), i(other.i), query_length(other.query_length),
+      target_length(other.target_length), ok(other.ok), keep(other.keep),
+      edit_cigar(other.edit_cigar) {
+    other.edit_cigar = {nullptr, 0, 0};
+}
+
+// Copy assignment operator
+alignment_t& alignment_t::operator=(const alignment_t& other) {
+    if (this != &other) {
+        j = other.j;
+        i = other.i;
+        query_length = other.query_length;
+        target_length = other.target_length;
+        ok = other.ok;
+        keep = other.keep;
+
+        free(edit_cigar.cigar_ops);
+        if (other.edit_cigar.cigar_ops) {
+            edit_cigar.cigar_ops = (char*)malloc((other.edit_cigar.end_offset - other.edit_cigar.begin_offset) * sizeof(char));
+            memcpy(edit_cigar.cigar_ops, other.edit_cigar.cigar_ops + other.edit_cigar.begin_offset, 
+                   (other.edit_cigar.end_offset - other.edit_cigar.begin_offset) * sizeof(char));
+        } else {
+            edit_cigar.cigar_ops = nullptr;
+        }
+        edit_cigar.begin_offset = 0;
+        edit_cigar.end_offset = other.edit_cigar.end_offset - other.edit_cigar.begin_offset;
+    }
+    return *this;
 }
+
+// Move assignment operator
+alignment_t& alignment_t::operator=(alignment_t&& other) noexcept {
+    if (this != &other) {
+        j = other.j;
+        i = other.i;
+        query_length = other.query_length;
+        target_length = other.target_length;
+        ok = other.ok;
+        keep = other.keep;
+
+        free(edit_cigar.cigar_ops);
+        edit_cigar = other.edit_cigar;
+        other.edit_cigar = {nullptr, 0, 0};
+    }
+    return *this;
+}
+
 //void alignment_t::display(void) {
 //    std::cerr << j << " " << i << " " << query_length << " "
 //              << target_length << " " << ok << std::endl;
@@ -581,6 +645,51 @@ void wflign_edit_cigar_copy(
     cigar_dst->end_offset = cigar_length;
     memcpy(cigar_dst->cigar_ops,cigar_ops,cigar_length);
 }
+
+int calculate_alignment_score(const wflign_cigar_t& cigar, const wflign_penalties_t& penalties) {
+    int score = 0;
+    char prev_op = '\0';
+    int gap_length = 0;
+
+    auto process_gap = [&](char op, int length) {
+        switch (op) {
+            case 'M':
+                // Match is free (best case)
+                break;
+            case 'X':
+                score += length * penalties.mismatch;
+                break;
+            case 'I':
+            case 'D':
+                score += penalties.gap_opening1 + penalties.gap_extension1;
+                if (length > 1) {
+                    score += std::min(
+                        penalties.gap_extension1 * (length - 1),
+                        penalties.gap_opening2 + penalties.gap_extension2 * (length - 1)
+                    );
+                }
+                break;
+        }
+    };
+
+    for (int i = cigar.begin_offset; i <= cigar.end_offset; ++i) {
+        char op = (i < cigar.end_offset) ? cigar.cigar_ops[i] : '\0';  // '\0' to process last gap
+
+        if (op != prev_op || i == cigar.end_offset) {
+            if (gap_length > 0) {
+                process_gap(prev_op, gap_length);
+            }
+            gap_length = 1;
+        } else {
+            ++gap_length;
+        }
+
+        prev_op = op;
+    }
+
+    return score;
+}
+
 /*
 // No more necessary
 bool hack_cigar(wfa::cigar_t &cigar, const char *query, const char *target,

diff --git a/src/common/wflign/src/wflign_alignment.hpp b/src/common/wflign/src/wflign_alignment.hpp
@@ -33,23 +33,29 @@ class alignment_t {
 public:
     int j;
     int i;
-    uint16_t query_length;
-    uint16_t target_length;
-    bool ok = false;
-    bool keep = false;
+    int query_length;
+    int target_length;
+    bool ok;
+    bool keep;
     wflign_cigar_t edit_cigar;
-    // Setup
+    // Default constructor
     alignment_t();
+    // Destructor
     ~alignment_t();
-    // Utils
-//    bool validate(
-//            const char* query,
-//            const char* target);
+    // Copy constructor
+    alignment_t(const alignment_t& other);
+    // Move constructor
+    alignment_t(alignment_t&& other) noexcept;
+    // Copy assignment operator
+    alignment_t& operator=(const alignment_t& other);
+    // Move assignment operator
+    alignment_t& operator=(alignment_t&& other) noexcept;
+    // Trim functions
     void trim_front(int query_trim);
     void trim_back(int query_trim);
-//    // Display
-//    void display(void);
 };
+
+
 /*
  * Wflign Trace-Pos: Links a position in a traceback matrix to its edit
  */
@@ -133,4 +139,6 @@ void wflign_edit_cigar_copy(
         wfa::WFAligner& wf_aligner,
         wflign_cigar_t* const cigar_dst);
 
-#endif /* WFLIGN_ALIGNMENT_HPP_ */
+int calculate_alignment_score(const wflign_cigar_t& cigar, const wflign_penalties_t& penalties);
+
+#endif /* WFLIGN_ALIGNMENT_HPP_ */