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

@@ -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);

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);
 

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);

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);
 

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;

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 */

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));