Fixed bug when using biwfa and wf-adaptive heuristics (off-grid indel…

…2indel breakpoint)
smarco · Jun 29, 2023 · 2eae6a5 · 2eae6a5
1 parent 5860976
commit 2eae6a5
Show file tree

Hide file tree

Showing 3 changed files with 26 additions and 8 deletions.
diff --git a/alignment/cigar.c b/alignment/cigar.c
@@ -68,7 +68,10 @@ cigar_t* cigar_new(
   cigar->begin_offset = 0;
   cigar->end_offset = 0;
   cigar->score = INT32_MIN;
+  cigar->end_v = -1;
+  cigar->end_h = -1;
   // CIGAR
+  cigar->cigar_length = 0;
   cigar->cigar_buffer = calloc(max_operations,sizeof(uint32_t));
   // Return
   return cigar;
@@ -371,7 +374,7 @@ void cigar_copy(
          cigar_src->end_offset-cigar_src->begin_offset);
 }
 void cigar_discover_mismatches(
-    char* const pattern,
+    const char* const pattern,
     const int pattern_length,
     const char* const text,
     const int text_length,
@@ -758,9 +761,15 @@ int cigar_sprint_SAM_CIGAR(
   // Print CIGAR-operations
   int i, cursor = 0;
   for (i=0;i<cigar_length;++i) {
-    cursor += sprintf(buffer+cursor,"%d%c",
-        cigar_buffer[i]>>4,
-        "MIDN---=X"[cigar_buffer[i]&0xf]);
+    const int op_idx = cigar_buffer[i] & 0xf;
+    if (op_idx <= 8) {
+      cursor += sprintf(buffer+cursor,"%d%c",
+          cigar_buffer[i]>>4,
+          "MIDN---=X"[cigar_buffer[i]&0xf]);
+    } else {
+      cursor += sprintf(buffer+cursor,"%d%c",
+          cigar_buffer[i]>>4,'?');
+    }
   }
   // Return
   buffer[cursor] = '\0';

diff --git a/alignment/cigar.h b/alignment/cigar.h
@@ -127,9 +127,9 @@ void cigar_copy(
     cigar_t* const cigar_src);
 
 void cigar_discover_mismatches(
-    char* const pattern,
+    const char* const pattern,
     const int pattern_length,
-    char* const text,
+    const char* const text,
     const int text_length,
     cigar_t* const cigar);
 

diff --git a/wavefront/wavefront_bialign.c b/wavefront/wavefront_bialign.c
@@ -220,16 +220,25 @@ void wavefront_bialign_breakpoint_indel2indel(
     const int dh_0 = WAVEFRONT_H(k_0,doffset_0);
     const int dh_1 = WAVEFRONT_H(k_1,doffset_1);
     // Check breakpoint d2d
-    if (dh_0 + dh_1 >= text_length && score_0 + score_1 - gap_open < breakpoint->score &&
-        dh_0 <= text_length && dh_1 <= text_length) {
+    if (dh_0 + dh_1 >= text_length && score_0 + score_1 - gap_open < breakpoint->score) {
       if (breakpoint_forward) {
+        // Check out-of-bounds coordinates
+        const int v = WAVEFRONT_V(k_0,dh_0);
+        const int h = WAVEFRONT_H(k_0,dh_0);
+        if (v > pattern_length || h > text_length) continue;
+        // Set breakpoint
         breakpoint->score_forward = score_0;
         breakpoint->score_reverse = score_1;
         breakpoint->k_forward = k_0;
         breakpoint->k_reverse = k_1;
         breakpoint->offset_forward = dh_0;
         breakpoint->offset_reverse = dh_1;
       } else {
+        // Check out-of-bounds coordinates
+        const int v = WAVEFRONT_V(k_1,dh_1);
+        const int h = WAVEFRONT_H(k_1,dh_1);
+        if (v > pattern_length || h > text_length) continue;
+        // Set breakpoint
         breakpoint->score_forward = score_1;
         breakpoint->score_reverse = score_0;
         breakpoint->k_forward = k_1;