bam_aux.c

   1 #include <ctype.h>
   2 #include "bam.h"
   3 #include "khash.h"
   4 KHASH_MAP_INIT_STR(s, int)
   5
   6 void bam_aux_append(bam1_t *b, const char tag[2], char type, int len, uint8_t *data)
   7 {
   8         int ori_len = b->data_len;
   9         b->data_len += 3 + len;
  10         b->l_aux += 3 + len;
  11         if (b->m_data < b->data_len) {
  12                 b->m_data = b->data_len;
  13                 kroundup32(b->m_data);
  14                 b->data = (uint8_t*)realloc(b->data, b->m_data);
  15         }
  16         b->data[ori_len] = tag[0]; b->data[ori_len + 1] = tag[1];
  17         b->data[ori_len + 2] = type;
  18         memcpy(b->data + ori_len + 3, data, len);
  19 }
  20
  21 uint8_t *bam_aux_get_core(bam1_t *b, const char tag[2])
  22 {
  23         return bam_aux_get(b, tag);
  24 }
  25
  26 uint8_t *bam_aux_get(bam1_t *b, const char tag[2])
  27 {
  28         uint8_t *s;
  29         int y = tag[0]<<8 | tag[1];
  30         s = bam1_aux(b);
  31         while (s < b->data + b->data_len) {
  32                 int type, x = (int)s[0]<<8 | s[1];
  33                 s += 2;
  34                 if (x == y) return s;
  35                 type = toupper(*s); ++s;
  36                 if (type == 'C') ++s;
  37                 else if (type == 'S') s += 2;
  38                 else if (type == 'I' || type == 'F') s += 4;
  39                 else if (type == 'D') s += 8;
  40                 else if (type == 'Z' || type == 'H') { while (*s) putchar(*s++); ++s; }
  41         }
  42         return 0;
  43 }
  44
  45 void bam_init_header_hash(bam_header_t *header)
  46 {
  47         if (header->hash == 0) {
  48                 int ret, i;
  49                 khiter_t iter;
  50                 khash_t(s) *h;
  51                 header->hash = h = kh_init(s);
  52                 for (i = 0; i < header->n_targets; ++i) {
  53                         iter = kh_put(s, h, header->target_name[i], &ret);
  54                         kh_value(h, iter) = i;
  55                 }
  56         }
  57 }
  58
  59 void bam_destroy_header_hash(bam_header_t *header)
  60 {
  61         if (header->hash)
  62                 kh_destroy(s, (khash_t(s)*)header->hash);
  63 }
  64
  65 int32_t bam_get_tid(const bam_header_t *header, const char *seq_name)
  66 {
  67         khint_t k;
  68         khash_t(s) *h = (khash_t(s)*)header->hash;
  69         k = kh_get(s, h, seq_name);
  70         return k == kh_end(h)? -1 : kh_value(h, k);
  71 }
  72
  73 void bam_parse_region(bam_header_t *header, const char *str, int *ref_id, int *begin, int *end)
  74 {
  75         char *s, *p;
  76         int i, l, k;
  77         khiter_t iter;
  78         khash_t(s) *h;
  79
  80         bam_init_header_hash(header);
  81         h = (khash_t(s)*)header->hash;
  82
  83         l = strlen(str);
  84         p = s = (char*)malloc(l+1);
  85         /* squeeze out "," */
  86         for (i = k = 0; i != l; ++i)
  87                 if (str[i] != ',' && !isspace(str[i])) s[k++] = str[i];
  88         s[k] = 0;
  89         for (i = 0; i != k; ++i) if (s[i] == ':') break;
  90         s[i] = 0;
  91         iter = kh_get(s, h, s); /* get the ref_id */
  92         if (iter == kh_end(h)) { // name not found
  93                 *ref_id = -1; free(s);
  94                 return;
  95         }
  96         *ref_id = kh_value(h, iter);
  97         if (i == k) { /* dump the whole sequence */
  98                 *begin = 0; *end = 1<<29; free(s);
  99                 return;
 100         }
 101         for (p = s + i + 1; i != k; ++i) if (s[i] == '-') break;
 102         *begin = atoi(p);
 103         if (i < k) {
 104                 p = s + i + 1;
 105                 *end = atoi(p);
 106         } else *end = 1<<29;
 107         if (*begin > 0) --*begin;
 108         assert(*begin <= *end);
 109         free(s);
 110 }
 111
 112 int32_t bam_aux2i(const uint8_t *s)
 113 {
 114         int type;
 115         if (s == 0) return 0;
 116         type = *s++;
 117         if (type == 'c') return (int32_t)*(int8_t*)s;
 118         else if (type == 'C') return (int32_t)*(uint8_t*)s;
 119         else if (type == 's') return (int32_t)*(int16_t*)s;
 120         else if (type == 'S') return (int32_t)*(uint16_t*)s;
 121         else if (type == 'i' || type == 'I') return *(int32_t*)s;
 122         else return 0;
 123 }
 124
 125 float bam_aux2f(const uint8_t *s)
 126 {
 127         int type;
 128         type = *s++;
 129         if (s == 0) return 0.0;
 130         if (type == 'f') return *(float*)s;
 131         else return 0.0;
 132 }
 133
 134 double bam_aux2d(const uint8_t *s)
 135 {
 136         int type;
 137         type = *s++;
 138         if (s == 0) return 0.0;
 139         if (type == 'd') return *(double*)s;
 140         else return 0.0;
 141 }
 142
 143 char bam_aux2A(const uint8_t *s)
 144 {
 145         int type;
 146         type = *s++;
 147         if (s == 0) return 0;
 148         if (type == 'A') return *(char*)s;
 149         else return 0;
 150 }
 151
 152 char *bam_aux2Z(const uint8_t *s)
 153 {
 154         int type;
 155         type = *s++;
 156         if (s == 0) return 0;
 157         if (type == 'Z' || type == 'H') return (char*)s;
 158         else return 0;
 159 }
 160
 161 char bam_aux_getCSi(bam1_t *b, int i)
 162 {
 163         uint8_t *c = bam_aux_get(b, "CS");
 164         char *cs = NULL;
 165
 166         // return the base if the tag was not found
 167         if(0 == c) return 0;
 168
 169         cs = bam_aux2Z(c);
 170         // adjust for strandedness and leading adaptor
 171         if(bam1_strand(b)) i = strlen(cs) - 1 - i;
 172         else i++;
 173         return cs[i];
 174 }
 175
 176 char bam_aux_getCQi(bam1_t *b, int i)
 177 {
 178         uint8_t *c = bam_aux_get(b, "CQ");
 179         char *cq = NULL;
 180
 181         // return the base if the tag was not found
 182         if(0 == c) return 0;
 183
 184         cq = bam_aux2Z(c);
 185         // adjust for strandedness
 186         if(bam1_strand(b)) i = strlen(cq) - 1 - i;
 187         return cq[i];
 188 }
 189
 190 char bam_aux_nt2int(char a)
 191 {
 192         switch(toupper(a)) {
 193                 case 'A':
 194                         return 0;
 195                         break;
 196                 case 'C':
 197                         return 1;
 198                         break;
 199                 case 'G':
 200                         return 2;
 201                         break;
 202                 case 'T':
 203                         return 3;
 204                         break;
 205                 default:
 206                         return 4;
 207                         break;
 208         }
 209 }
 210
 211 char bam_aux_ntnt2cs(char a, char b)
 212 {
 213         a = bam_aux_nt2int(a);
 214         b = bam_aux_nt2int(b);
 215         if(4 == a || 4 == b) return '4';
 216         return "0123"[(int)(a ^ b)];
 217 }
 218
 219 char bam_aux_getCEi(bam1_t *b, int i)
 220 {
 221         int cs_i;
 222         uint8_t *c = bam_aux_get(b, "CS");
 223         char *cs = NULL;
 224         char prev_b, cur_b;
 225         char cur_color, cor_color;
 226
 227         // return the base if the tag was not found
 228         if(0 == c) return 0;
 229
 230         cs = bam_aux2Z(c);
 231
 232         // adjust for strandedness and leading adaptor
 233         if(bam1_strand(b)) { //reverse strand
 234                 cs_i = strlen(cs) - 1 - i;
 235                 // get current color
 236                 cur_color = cs[cs_i];
 237                 // get previous base
 238                 prev_b = (0 == cs_i) ? cs[0] : bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i+1)];
 239                 // get current base
 240                 cur_b = bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i)];
 241         }
 242         else {
 243                 cs_i=i+1;
 244                 // get current color
 245                 cur_color = cs[cs_i];
 246                 // get previous base
 247                 prev_b = (0 == i) ? cs[0] : bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i-1)];
 248                 // get current base
 249                 cur_b = bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i)];
 250         }
 251
 252         // corrected color
 253         cor_color = bam_aux_ntnt2cs(prev_b, cur_b);
 254
 255         if(cur_color == cor_color) {
 256                 return '-';
 257         }
 258         else {
 259                 return cur_color;
 260         }
 261 }