bam_md.c

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