From 011fa103c7af400670b23c23bea2e8c8c1391798 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Fri, 17 Jul 2009 13:46:56 +0000 Subject: [PATCH] * samtools-0.1.5-10 (r407) * implemented bam_aux_del() to remove a tag * fillmd: generate the NM tag * fillmd: cmd interface improvement --- bam.h | 2 +- bam_aux.c | 31 +++++++++++++++++++------ bam_md.c | 68 +++++++++++++++++++++++++++++++++++++------------------ bamtk.c | 4 ++-- sam.h | 2 +- 5 files changed, 74 insertions(+), 33 deletions(-) diff --git a/bam.h b/bam.h index 64867f1..451b65a 100644 --- a/bam.h +++ b/bam.h @@ -604,8 +604,8 @@ extern "C" { char bam_aux2A(const uint8_t *s); char *bam_aux2Z(const uint8_t *s); + int bam_aux_del(bam1_t *b, uint8_t *s); void bam_aux_append(bam1_t *b, const char tag[2], char type, int len, uint8_t *data); - uint8_t *bam_aux_get_core(bam1_t *b, const char tag[2]); // an alias of bam_aux_get() /*! diff --git a/bam_aux.c b/bam_aux.c index 7482500..894ae62 100644 --- a/bam_aux.c +++ b/bam_aux.c @@ -25,24 +25,41 @@ uint8_t *bam_aux_get_core(bam1_t *b, const char tag[2]) return bam_aux_get(b, tag); } +#define __skip_tag(s) do { \ + int 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)) ++(s); ++(s); } \ + } while (0) + uint8_t *bam_aux_get(const 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]; + int 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) ++s; ++s; } + __skip_tag(s); } return 0; } +// s MUST BE returned by bam_aux_get() +int bam_aux_del(bam1_t *b, uint8_t *s) +{ + uint8_t *p, *aux; + aux = bam1_aux(b); + p = s - 2; + __skip_tag(s); + memmove(p, s, b->l_aux - (s - aux)); + b->data_len -= s - p; + b->l_aux -= s - p; + return 0; +} void bam_init_header_hash(bam_header_t *header) { diff --git a/bam_md.c b/bam_md.c index a20f9b3..45bb1db 100644 --- a/bam_md.c +++ b/bam_md.c @@ -3,7 +3,7 @@ #include #include #include "faidx.h" -#include "bam.h" +#include "sam.h" #include "kstring.h" void bam_fillmd1(bam1_t *b, char *ref, int is_equal) @@ -13,11 +13,13 @@ void bam_fillmd1(bam1_t *b, char *ref, int is_equal) bam1_core_t *c = &b->core; int i, x, y, u = 0; kstring_t *str; - uint8_t *old_md; + uint8_t *old_md, *old_nm; + int32_t old_nm_i = -1, nm = 0; old_md = bam_aux_get(b, "MD"); + old_nm = bam_aux_get(b, "NM"); if (c->flag & BAM_FUNMAP) return; - if (old_md && !is_equal) return; // no need to add MD + if (old_nm) old_nm_i = bam_aux2i(old_nm); str = (kstring_t*)calloc(1, sizeof(kstring_t)); for (i = y = 0, x = c->pos; i < c->n_cigar; ++i) { int j, l = cigar[i]>>4, op = cigar[i]&0xf; @@ -26,13 +28,13 @@ void bam_fillmd1(bam1_t *b, char *ref, int is_equal) int z = y + j; int c1 = bam1_seqi(seq, z), c2 = bam_nt16_table[(int)ref[x+j]]; if (ref[x+j] == 0) break; // out of boundary - if ((c1 == c2 && c1 != 15 && c2 != 15) || c1 == 0) { + if ((c1 == c2 && c1 != 15 && c2 != 15) || c1 == 0) { // a match if (is_equal) seq[z/2] &= (z&1)? 0xf0 : 0x0f; ++u; } else { ksprintf(str, "%d", u); kputc(ref[x+j], str); - u = 0; + u = 0; ++nm; } } if (j < l) break; @@ -46,14 +48,21 @@ void bam_fillmd1(bam1_t *b, char *ref, int is_equal) } u = 0; if (j < l) break; - x += l; + x += l; nm += l; } else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) { y += l; + if (op == BAM_CINS) nm += l; } else if (op == BAM_CREF_SKIP) { x += l; } } ksprintf(str, "%d", u); + if (!old_nm) bam_aux_append(b, "NM", 'i', 4, (uint8_t*)&nm); + else if (nm != old_nm_i) { + fprintf(stderr, "[bam_fillmd1] different NM for read '%s': %d != %d\n", bam1_qname(b), old_nm_i, nm); + bam_aux_del(b, old_nm); + bam_aux_append(b, "NM", 'i', 4, (uint8_t*)&nm); + } if (!old_md) bam_aux_append(b, "MD", 'Z', str->l + 1, (uint8_t*)str->s); else { int is_diff = 0; @@ -71,47 +80,62 @@ void bam_fillmd1(bam1_t *b, char *ref, int is_equal) int bam_fillmd(int argc, char *argv[]) { - int c, is_equal = 0, tid = -2, ret, len; - bamFile fp, fpout = 0; - bam_header_t *header; + int c, is_equal = 0, tid = -2, ret, len, is_bam_out, is_sam_in, is_uncompressed; + samfile_t *fp, *fpout = 0; faidx_t *fai; - char *ref = 0; + char *ref = 0, mode_w[8], mode_r[8]; bam1_t *b; - while ((c = getopt(argc, argv, "e")) >= 0) { + is_bam_out = is_sam_in = is_uncompressed = 0; + mode_w[0] = mode_r[0] = 0; + strcpy(mode_r, "r"); strcpy(mode_w, "w"); + while ((c = getopt(argc, argv, "eubS")) >= 0) { switch (c) { case 'e': is_equal = 1; break; + case 'b': is_bam_out = 1; break; + case 'u': is_uncompressed = is_bam_out = 1; break; + case 'S': is_sam_in = 1; break; default: fprintf(stderr, "[bam_fillmd] unrecognized option '-%c'\n", c); return 1; } } + if (!is_sam_in) strcat(mode_r, "b"); + if (is_bam_out) strcat(mode_w, "b"); + else strcat(mode_w, "h"); + if (is_uncompressed) strcat(mode_w, "u"); if (optind + 1 >= argc) { - fprintf(stderr, "Usage: bam fillmd [-e] \n"); + fprintf(stderr, "\n"); + fprintf(stderr, "Usage: samtools fillmd [-eubS] \n\n"); + fprintf(stderr, "Options: -e change identical bases to '='\n"); + fprintf(stderr, " -u uncompressed BAM output (for piping)\n"); + fprintf(stderr, " -b compressed BAM output\n"); + fprintf(stderr, " -S the input is SAM with header\n\n"); + return 1; + } + fp = samopen(argv[optind], mode_r, 0); + if (fp == 0) return 1; + if (is_sam_in && (fp->header == 0 || fp->header->n_targets == 0)) { + fprintf(stderr, "[bam_fillmd] input SAM does not have header. Abort!\n"); return 1; } - fp = strcmp(argv[optind], "-")? bam_open(argv[optind], "r") : bam_dopen(fileno(stdin), "r"); - assert(fp); - header = bam_header_read(fp); - fpout = bam_dopen(fileno(stdout), "w"); - bam_header_write(fpout, header); + fpout = samopen("-", mode_w, fp->header); fai = fai_load(argv[optind+1]); b = bam_init1(); - while ((ret = bam_read1(fp, b)) >= 0) { + while ((ret = samread(fp, b)) >= 0) { if (b->core.tid >= 0) { if (tid != b->core.tid) { free(ref); - ref = fai_fetch(fai, header->target_name[b->core.tid], &len); + ref = fai_fetch(fai, fp->header->target_name[b->core.tid], &len); tid = b->core.tid; } bam_fillmd1(b, ref, is_equal); } - bam_write1(fpout, b); + samwrite(fpout, b); } bam_destroy1(b); free(ref); fai_destroy(fai); - bam_header_destroy(header); - bam_close(fp); bam_close(fpout); + samclose(fp); samclose(fpout); return 0; } diff --git a/bamtk.c b/bamtk.c index fb000a6..f162bc3 100644 --- a/bamtk.c +++ b/bamtk.c @@ -4,7 +4,7 @@ #include "bam.h" #ifndef PACKAGE_VERSION -#define PACKAGE_VERSION "0.1.5-9 (r405)" +#define PACKAGE_VERSION "0.1.5-10 (r407)" #endif int bam_taf2baf(int argc, char *argv[]); @@ -85,7 +85,7 @@ static int usage() fprintf(stderr, " rmdup remove PCR duplicates\n"); fprintf(stderr, " glfview print GLFv3 file\n"); fprintf(stderr, " flagstat simple stats\n"); - fprintf(stderr, " fillmd fill the MD tag and change identical base to =\n"); + fprintf(stderr, " fillmd fill the MD/NM tag and generate '=' bases\n"); fprintf(stderr, "\n"); return 1; } diff --git a/sam.h b/sam.h index 970cf2d..86b7576 100644 --- a/sam.h +++ b/sam.h @@ -48,7 +48,7 @@ extern "C" { "wbu" exclusively. @param aux auxiliary data; if mode[0]=='w', aux points to - bam_header_t; if strcmp(mode, "rb")==0 and @SQ header lines in SAM + bam_header_t; if strcmp(mode, "rb")!=0 and @SQ header lines in SAM are absent, aux points the file name of the list of the reference; aux is not used otherwise. -- 2.39.2