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