From b775f49057c7e96ecba692d9eafa38dea6791a2d Mon Sep 17 00:00:00 2001
From: Jacob Morrison <Jacob.Morrison@vai.org>
Date: Thu, 22 Feb 2024 15:59:15 -0500
Subject: [PATCH] Add extended CIGAR operations to BISCUIT

Only Match (M) was allowed before. Now Equal (=) and Different (X) are
able to work as well. While in theory, Equal should mean the query base
is the same as the reference base, some aligners use Equal to cover
bases where cytosine conversion occurred (i.e., C>T). Therefore, Match,
Equal, and Different are all treated the same during processing.
---
 src/bisc_utils.c | 6 ++++--
 src/bsconv.c     | 2 ++
 src/bsstrand.c   | 2 ++
 src/cinread.c    | 2 ++
 src/epiread.c    | 2 ++
 src/pileup.c     | 2 ++
 src/tview.c      | 2 ++
 7 files changed, 16 insertions(+), 2 deletions(-)

diff --git a/src/bisc_utils.c b/src/bisc_utils.c
index d27f5b2..6779270 100644
--- a/src/bisc_utils.c
+++ b/src/bisc_utils.c
@@ -86,6 +86,8 @@ uint32_t cnt_retention(refcache_t *rs, bam1_t *b, uint8_t bsstrand) {
         oplen = bam_cigar_oplen(bam_get_cigar(b)[i]);
         switch(op) {
             case BAM_CMATCH:
+            case BAM_CEQUAL:
+            case BAM_CDIFF:
                 for (j=0; j<oplen; ++j) {
                     rb = refcache_getbase_upcase(rs, rpos+j);
                     qb = bscall(b, qpos+j);
@@ -159,9 +161,7 @@ uint32_t get_mate_length(char *m_cigar) {
 }
 
 uint8_t infer_bsstrand(refcache_t *rs, bam1_t *b, uint32_t min_base_qual) {
-
     /* infer bsstrand from nC2T and nG2A on high quality bases */
-
     bam1_core_t *c = &b->core;
     uint32_t i, rpos = c->pos+1, qpos = 0;
     uint32_t op, oplen;
@@ -172,6 +172,8 @@ uint8_t infer_bsstrand(refcache_t *rs, bam1_t *b, uint32_t min_base_qual) {
         oplen = bam_cigar_oplen(bam_get_cigar(b)[i]);
         switch(op) {
             case BAM_CMATCH:
+            case BAM_CEQUAL:
+            case BAM_CDIFF:
                 for (j=0; j<oplen; ++j) {
                     rb = refcache_getbase_upcase(rs, rpos+j);
                     qb = bscall(b, qpos+j);
diff --git a/src/bsconv.c b/src/bsconv.c
index 860ae0f..a0ee771 100644
--- a/src/bsconv.c
+++ b/src/bsconv.c
@@ -65,6 +65,8 @@ int bsconv_func(bam1_t *b, samFile *out, bam_hdr_t *hdr, void *data) {
         uint32_t oplen = bam_cigar_oplen(bam_get_cigar(b)[i]);
         switch(op) {
             case BAM_CMATCH:
+            case BAM_CEQUAL:
+            case BAM_CDIFF:
                 for(j=0; j<oplen; ++j) {
                     rb = refcache_getbase_upcase(d->rs, rpos+j);
 
diff --git a/src/bsstrand.c b/src/bsstrand.c
index c617f51..f829e2e 100755
--- a/src/bsstrand.c
+++ b/src/bsstrand.c
@@ -82,6 +82,8 @@ int bsstrand_func(bam1_t *b, samFile *out, bam_hdr_t *header, void *data) {
         uint32_t oplen = bam_cigar_oplen(bam_get_cigar(b)[i]);
         switch(op) {
             case BAM_CMATCH:
+            case BAM_CEQUAL:
+            case BAM_CDIFF:
                 for(j=0; j<oplen; ++j) {
                     rbase = toupper(refcache_getbase(d->rs, rpos+j));
                     qbase = bscall(b, qpos+j);
diff --git a/src/cinread.c b/src/cinread.c
index f2ffa32..b1d0a20 100644
--- a/src/cinread.c
+++ b/src/cinread.c
@@ -72,6 +72,8 @@ int cinread_func(bam1_t *b, samFile *out, bam_hdr_t *hdr, void *data) {
         uint32_t oplen = bam_cigar_oplen(bam_get_cigar(b)[i]);
         switch(op) {
             case BAM_CMATCH:
+            case BAM_CEQUAL:
+            case BAM_CDIFF:
                 for(j=0; j<oplen; ++j) {
                     rb = refcache_getbase_upcase(d->rs, rpos+j);
 
diff --git a/src/epiread.c b/src/epiread.c
index f95c800..3a354ae 100644
--- a/src/epiread.c
+++ b/src/epiread.c
@@ -675,6 +675,8 @@ static void *process_func(void *data) {
                 uint32_t oplen = bam_cigar_oplen(bam_get_cigar(b)[i]);
                 switch(op) {
                     case BAM_CMATCH:
+                    case BAM_CEQUAL:
+                    case BAM_CDIFF:
                         for (j=0; j<oplen; ++j) {
                             // Query positions
                             qj = qpos + j;
diff --git a/src/pileup.c b/src/pileup.c
index aef6035..7fa8da0 100644
--- a/src/pileup.c
+++ b/src/pileup.c
@@ -752,6 +752,8 @@ static void *process_func(void *_result) {
                     char qb, rb;
                     switch(op) {
                         case BAM_CMATCH:
+                        case BAM_CEQUAL:
+                        case BAM_CDIFF:
                             for (j=0; j<oplen; ++j) {
 
                                 if (rpos+j<w.beg || rpos+j>=w.end) continue; /* include begin but not end */
diff --git a/src/tview.c b/src/tview.c
index 10c0018..c451b29 100644
--- a/src/tview.c
+++ b/src/tview.c
@@ -335,6 +335,8 @@ static void draw_read1(rnode_t *nd, btview_t *tv, int readattr, int bss) {
 
     switch(op) {
     case BAM_CMATCH:
+    case BAM_CEQUAL:
+    case BAM_CDIFF:
       for (j=0; j<oplen; ++j) {
         if (rpos + j < (unsigned) tv->left_pos) continue;
         qb = toupper(bscall(b, qpos + j));