bam_color.c

   1 #include <ctype.h>
   2 #include "bam.h"
   3
   4 /*!
   5  @abstract     Get the color encoding the previous and current base
   6  @param b      pointer to an alignment
   7  @param i      The i-th position, 0-based
   8  @return       color
   9
  10  @discussion   Returns 0 no color information is found.
  11  */
  12 char bam_aux_getCSi(bam1_t *b, int i)
  13 {
  14         uint8_t *c = bam_aux_get(b, "CS");
  15         char *cs = NULL;
  16
  17         // return the base if the tag was not found
  18         if(0 == c) return 0;
  19
  20         cs = bam_aux2Z(c);
  21         // adjust for strandedness and leading adaptor
  22         if(bam1_strand(b)) {
  23             i = strlen(cs) - 1 - i;
  24             // adjust for leading hard clip
  25             uint32_t cigar = bam1_cigar(b)[0];
  26             if((cigar & BAM_CIGAR_MASK) == BAM_CHARD_CLIP) {
  27                 i -= cigar >> BAM_CIGAR_SHIFT;
  28             }
  29         } else { i++; }
  30         return cs[i];
  31 }
  32
  33 /*!
  34  @abstract     Get the color quality of the color encoding the previous and current base
  35  @param b      pointer to an alignment
  36  @param i      The i-th position, 0-based
  37  @return       color quality
  38
  39  @discussion   Returns 0 no color information is found.
  40  */
  41 char bam_aux_getCQi(bam1_t *b, int i)
  42 {
  43         uint8_t *c = bam_aux_get(b, "CQ");
  44         char *cq = NULL;
  45
  46         // return the base if the tag was not found
  47         if(0 == c) return 0;
  48
  49         cq = bam_aux2Z(c);
  50         // adjust for strandedness
  51         if(bam1_strand(b)) {
  52             i = strlen(cq) - 1 - i;
  53             // adjust for leading hard clip
  54             uint32_t cigar = bam1_cigar(b)[0];
  55             if((cigar & BAM_CIGAR_MASK) == BAM_CHARD_CLIP) {
  56                 i -= (cigar >> BAM_CIGAR_SHIFT);
  57             }
  58         }
  59         return cq[i];
  60 }
  61
  62 char bam_aux_nt2int(char a)
  63 {
  64         switch(toupper(a)) {
  65                 case 'A':
  66                         return 0;
  67                         break;
  68                 case 'C':
  69                         return 1;
  70                         break;
  71                 case 'G':
  72                         return 2;
  73                         break;
  74                 case 'T':
  75                         return 3;
  76                         break;
  77                 default:
  78                         return 4;
  79                         break;
  80         }
  81 }
  82
  83 char bam_aux_ntnt2cs(char a, char b)
  84 {
  85         a = bam_aux_nt2int(a);
  86         b = bam_aux_nt2int(b);
  87         if(4 == a || 4 == b) return '4';
  88         return "0123"[(int)(a ^ b)];
  89 }
  90
  91 /*!
  92  @abstract     Get the color error profile at the give position
  93  @param b      pointer to an alignment
  94  @return       the original color if the color was an error, '-' (dash) otherwise
  95
  96  @discussion   Returns 0 no color information is found.
  97  */
  98 char bam_aux_getCEi(bam1_t *b, int i)
  99 {
 100         int cs_i;
 101         uint8_t *c = bam_aux_get(b, "CS");
 102         char *cs = NULL;
 103         char prev_b, cur_b;
 104         char cur_color, cor_color;
 105
 106         // return the base if the tag was not found
 107         if(0 == c) return 0;
 108
 109         cs = bam_aux2Z(c);
 110
 111         // adjust for strandedness and leading adaptor
 112         if(bam1_strand(b)) { //reverse strand
 113                 cs_i = strlen(cs) - 1 - i;
 114                 // adjust for leading hard clip
 115                 uint32_t cigar = bam1_cigar(b)[0];
 116                 if((cigar & BAM_CIGAR_MASK) == BAM_CHARD_CLIP) {
 117                         cs_i -= cigar >> BAM_CIGAR_SHIFT;
 118                 }
 119                 // get current color
 120                 cur_color = cs[cs_i];
 121                 // get previous base.  Note: must rc adaptor
 122                 prev_b = (cs_i == 1) ? "TGCAN"[(int)bam_aux_nt2int(cs[0])] : bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i+1)];
 123                 // get current base
 124                 cur_b = bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i)];
 125         }
 126         else {
 127                 cs_i=i+1;
 128                 // get current color
 129                 cur_color = cs[cs_i];
 130                 // get previous base
 131                 prev_b = (0 == i) ? cs[0] : bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i-1)];
 132                 // get current base
 133                 cur_b = bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i)];
 134         }
 135
 136         // corrected color
 137         cor_color = bam_aux_ntnt2cs(prev_b, cur_b);
 138
 139         if(cur_color == cor_color) {
 140                 return '-';
 141         }
 142         else {
 143                 return cur_color;
 144         }
 145 }