X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bam_aux.c;h=2fa9ac2e83a7a00e2e83a523fe687448e3c4255e;hb=41586185b43962c1ffb3f82e93066a480254aa17;hp=081f07bc367a57fe2bd139b4fd82ea73421b5d5d;hpb=f93dae0d03856955f9424e8b2aaf261304ca647e;p=samtools.git diff --git a/bam_aux.c b/bam_aux.c index 081f07b..2fa9ac2 100644 --- a/bam_aux.c +++ b/bam_aux.c @@ -1,9 +1,42 @@ #include #include "bam.h" #include "khash.h" -KHASH_MAP_INIT_INT(aux, uint8_t*) KHASH_MAP_INIT_STR(s, int) +void bam_aux_append(bam1_t *b, const char tag[2], char type, int len, uint8_t *data) +{ + int ori_len = b->data_len; + b->data_len += 3 + len; + b->l_aux += 3 + len; + if (b->m_data < b->data_len) { + b->m_data = b->data_len; + kroundup32(b->m_data); + b->data = (uint8_t*)realloc(b->data, b->m_data); + } + b->data[ori_len] = tag[0]; b->data[ori_len + 1] = tag[1]; + b->data[ori_len + 2] = type; + memcpy(b->data + ori_len + 3, data, len); +} + +uint8_t *bam_aux_get(bam1_t *b, const char tag[2]) +{ + uint8_t *s; + int y = tag[0]<<8 | tag[1]; + s = bam1_aux(b); + while (s < b->data + b->data_len) { + int type, x = (int)s[0]<<8 | s[1]; + s += 2; + if (x == y) return s; + type = toupper(*s); ++s; + if (type == 'C') ++s; + else if (type == 'S') s += 2; + else if (type == 'I' || type == 'F') s += 4; + else if (type == 'D') s += 8; + else if (type == 'Z' || type == 'H') { while (*s) putchar(*s++); ++s; } + } + return 0; +} + void bam_init_header_hash(bam_header_t *header) { if (header->hash == 0) { @@ -41,7 +74,7 @@ void bam_parse_region(bam_header_t *header, const char *str, int *ref_id, int *b bam_init_header_hash(header); h = (khash_t(s)*)header->hash; - + l = strlen(str); p = s = (char*)malloc(l+1); /* squeeze out "," */ @@ -71,90 +104,153 @@ void bam_parse_region(bam_header_t *header, const char *str, int *ref_id, int *b free(s); } -void bam_aux_init(bam1_t *b) -{ - khash_t(aux) *h; - uint8_t *s; - if (b->hash == 0) { - h = kh_init(aux); - b->hash = h; - } else { - h = (khash_t(aux)*)b->hash; - kh_clear(aux, h); - } - s = bam1_aux(b); - while (s < b->data + b->data_len) { - uint32_t x = (uint32_t)s[0]<<8 | s[1]; - int ret, type; - khint_t k; - s += 2; type = toupper(*s); ++s; - k = kh_put(aux, h, x, &ret); - kh_value(h, k) = s; - if (type == 'C') ++s; - else if (type == 'S') s += 2; - else if (type == 'I') s += 4; - else if (type == 'F') s += 4; - else if (type == 'Z') { while (*s) putchar(*s++); ++s; } - } -} -void bam_aux_destroy(bam1_t *b) -{ - khash_t(aux) *h = (khash_t(aux)*)b->hash; - kh_destroy(aux, h); - b->hash = 0; -} -static uint8_t *bam_aux_get_core(bam1_t *b, const char tag[2]) -{ - uint32_t x = (uint32_t)tag[0]<<8 | tag[1]; - khint_t k; - khash_t(aux) *h; - if (b->hash == 0) bam_aux_init(b); - h = (khash_t(aux)*)b->hash; - k = kh_get(aux, h, x); - if (k == kh_end(h)) return 0; - return kh_value(h, k); -} -int32_t bam_aux_geti(bam1_t *b, const char tag[2], int *err) +int32_t bam_aux2i(const uint8_t *s) { int type; - uint8_t *s = bam_aux_get_core(b, tag); - *err = 0; - if (s == 0) { *err = -1; return 0; } + if (s == 0) return 0; type = *s++; if (type == 'c') return (int32_t)*(int8_t*)s; else if (type == 'C') return (int32_t)*(uint8_t*)s; else if (type == 's') return (int32_t)*(int16_t*)s; else if (type == 'S') return (int32_t)*(uint16_t*)s; else if (type == 'i' || type == 'I') return *(int32_t*)s; - else { *err = -2; return 0; } + else return 0; } -float bam_aux_getf(bam1_t *b, const char tag[2], int *err) + +float bam_aux2f(const uint8_t *s) { int type; - uint8_t *s = bam_aux_get_core(b, tag); - *err = 0; type = *s++; - if (s == 0) { *err = -1; return 0; } + if (s == 0) return 0.0; if (type == 'f') return *(float*)s; - else { *err = -2; return 0; } + else return 0.0; +} + +double bam_aux2d(const uint8_t *s) +{ + int type; + type = *s++; + if (s == 0) return 0.0; + if (type == 'd') return *(double*)s; + else return 0.0; } -char bam_aux_getc(bam1_t *b, const char tag[2], int *err) + +char bam_aux2A(const uint8_t *s) { int type; - uint8_t *s = bam_aux_get_core(b, tag); - *err = 0; type = *s++; - if (s == 0) { *err = -1; return 0; } - if (type == 'c') return *(char*)s; - else { *err = -2; return 0; } + if (s == 0) return 0; + if (type == 'A') return *(char*)s; + else return 0; } -char *bam_aux_getZH(bam1_t *b, const char tag[2], int *err) + +char *bam_aux2Z(const uint8_t *s) { int type; - uint8_t *s = bam_aux_get_core(b, tag); - *err = 0; type = *s++; - if (s == 0) { *err = -1; return 0; } + if (s == 0) return 0; if (type == 'Z' || type == 'H') return (char*)s; - else { *err = -2; return 0; } + else return 0; +} + +char bam_aux_getCSi(bam1_t *b, int i) +{ + uint8_t *c = bam_aux_get(b, "CS"); + char *cs = NULL; + + // return the base if the tag was not found + if(0 == c) return 0; + + cs = bam_aux2Z(c); + // adjust for strandedness and leading adaptor + if(bam1_strand(b)) i = strlen(cs) - 1 - i; + else i++; + return cs[i]; +} + +char bam_aux_getCQi(bam1_t *b, int i) +{ + uint8_t *c = bam_aux_get(b, "CQ"); + char *cq = NULL; + + // return the base if the tag was not found + if(0 == c) return 0; + + cq = bam_aux2Z(c); + // adjust for strandedness + if(bam1_strand(b)) i = strlen(cq) - 1 - i; + return cq[i]; +} + +char bam_aux_nt2int(char a) +{ + switch(toupper(a)) { + case 'A': + return 0; + break; + case 'C': + return 1; + break; + case 'G': + return 2; + break; + case 'T': + return 3; + break; + default: + return 4; + break; + } +} + +char bam_aux_ntnt2cs(char a, char b) +{ + a = bam_aux_nt2int(a); + b = bam_aux_nt2int(b); + if(4 == a || 4 == b) return '4'; + return "0123"[(int)(a ^ b)]; +} + +char bam_aux_getCEi(bam1_t *b, int i) +{ + int cs_i; + uint8_t *c = bam_aux_get(b, "CS"); + char *cs = NULL; + char prev_b, cur_b; + char cur_color, cor_color; + + // return the base if the tag was not found + if(0 == c) return 0; + + cs = bam_aux2Z(c); + + // adjust for strandedness and leading adaptor + if(bam1_strand(b)) { //reverse strand + cs_i = strlen(cs) - 1 - i; + // get current color + cur_color = cs[cs_i]; + // get previous base + prev_b = (0 == cs_i) ? cs[0] : bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i+1)]; + // get current base + cur_b = bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i)]; + } + else { + cs_i=i+1; + // get current color + cur_color = cs[cs_i]; + // get previous base + prev_b = (0 == i) ? cs[0] : bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i-1)]; + // get current base + cur_b = bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i)]; + } + + // corrected color + cor_color = bam_aux_ntnt2cs(prev_b, cur_b); + + if(cur_color == cor_color) { + return '-'; + } + else { + return cur_color; + } }