bam_md.c

   1 #include <unistd.h>
   2 #include "faidx.h"
   3 #include "bam.h"
   4 #include "kstring.h"
   5
   6 void bam_fillmd1(bam1_t *b, char *ref, int is_equal)
   7 {
   8         uint8_t *seq = bam1_seq(b);
   9         uint32_t *cigar = bam1_cigar(b);
  10         bam1_core_t *c = &b->core;
  11         int i, x, y, u = 0, has_md = 0;
  12         kstring_t *str;
  13
  14         if (bam_aux_get(b, "MD")) has_md = 1;
  15         if (has_md == 1 && !is_equal) return; // no need to add MD
  16         str = (kstring_t*)calloc(1, sizeof(kstring_t));
  17         if (c->flag & BAM_FUNMAP) return;
  18         for (i = y = 0, x = c->pos; i < c->n_cigar; ++i) {
  19                 int j, l = cigar[i]>>4, op = cigar[i]&0xf;
  20                 if (op == BAM_CMATCH) {
  21                         for (j = 0; j < l; ++j) {
  22                                 int z = y + j;
  23                                 int c1 = bam1_seqi(seq, z), c2 = bam_nt16_table[(int)ref[x+j]];
  24                                 if (ref[x+j] == 0) break; // out of boundary
  25                                 if ((c1 == c2 && c1 != 15 && c2 != 15) || c1 == 0) {
  26                                         if (is_equal) seq[z/2] &= (z&1)? 0xf0 : 0x0f;
  27                                         ++u;
  28                                 } else {
  29                                         ksprintf(str, "%d", u);
  30                                         kputc(ref[x+j], str);
  31                                         u = 0;
  32                                 }
  33                         }
  34                         if (j < l) break;
  35                         x += l; y += l;
  36                 } else if (op == BAM_CDEL) {
  37                         ksprintf(str, "%d", u);
  38                         kputc('^', str);
  39                         for (j = 0; j < l; ++j) {
  40                                 if (ref[x+j] == 0) break;
  41                                 kputc(ref[x+j], str);
  42                         }
  43                         if (j < l) break;
  44                         x += l;
  45                 } else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) {
  46                         y += l;
  47                 } else if (op == BAM_CREF_SKIP) {
  48                         x += l;
  49                 }
  50         }
  51         ksprintf(str, "%d", u);
  52         if (!has_md) bam_aux_append(b, "MD", 'Z', str->l + 1, (uint8_t*)str->s);
  53         free(str->s); free(str);
  54 }
  55
  56 int bam_fillmd(int argc, char *argv[])
  57 {
  58         int c, is_equal = 0, tid = -1, ret, len;
  59         bamFile fp, fpout = 0;
  60         bam_header_t *header;
  61         faidx_t *fai;
  62         char *ref = 0;
  63         bam1_t *b;
  64
  65         while ((c = getopt(argc, argv, "e")) >= 0) {
  66                 switch (c) {
  67                 case 'e': is_equal = 1; break;
  68                 default: fprintf(stderr, "[bam_fillmd] unrecognized option '-%c'\n", c); return 1;
  69                 }
  70         }
  71         if (optind + 1 >= argc) {
  72                 fprintf(stderr, "Usage: bam fillmd [-e] <aln.bam> <ref.fasta>\n");
  73                 return 1;
  74         }
  75         fp = strcmp(argv[optind], "-")? bam_open(argv[optind], "r") : bam_dopen(fileno(stdin), "r");
  76         assert(fp);
  77         header = bam_header_read(fp);
  78         fpout = bam_dopen(fileno(stdout), "w");
  79         bam_header_write(fpout, header);
  80         fai = fai_load(argv[optind+1]);
  81
  82         b = bam_init1();
  83         while ((ret = bam_read1(fp, b)) >= 0) {
  84                 if (tid != b->core.tid) {
  85                         free(ref);
  86                         ref = fai_fetch(fai, header->target_name[b->core.tid], &len);
  87                         tid = b->core.tid;
  88                 }
  89                 bam_fillmd1(b, ref, is_equal);
  90                 bam_write1(fpout, b);
  91         }
  92         bam_destroy1(b);
  93
  94         fai_destroy(fai);
  95         bam_header_destroy(header);
  96         bam_close(fp); bam_close(fpout);
  97         return 0;
  98 }