Also clean *.exe (for Cygwin users using this makefile).

[samtools.git] / bam_md.c

diff --git a/bam_md.c b/bam_md.c
index 45bb1db7531202caba7a68f5084fd58fa77de710..3ca730993cbc01065eb2c76bd66df1640808e032 100644 (file)
--- a/bam_md.c
+++ b/bam_md.c
@@ -16,10 +16,6 @@ void bam_fillmd1(bam1_t *b, char *ref, int is_equal)
        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_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;
@@ -57,12 +53,18 @@ void bam_fillmd1(bam1_t *b, char *ref, int is_equal)
                }
        }
        ksprintf(str, "%d", u);
+       // update NM
+       old_nm = bam_aux_get(b, "NM");
+       if (c->flag & BAM_FUNMAP) return;
+       if (old_nm) old_nm_i = bam_aux2i(old_nm);
        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);
+               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);
        }
+       // update MD
+       old_md = bam_aux_get(b, "MD");
        if (!old_md) bam_aux_append(b, "MD", 'Z', str->l + 1, (uint8_t*)str->s);
        else {
                int is_diff = 0;
@@ -72,8 +74,11 @@ void bam_fillmd1(bam1_t *b, char *ref, int is_equal)
                                        break;
                        if (i < str->l) is_diff = 1;
                } else is_diff = 1;
-               if (is_diff)
-                       fprintf(stderr, "[bam_fillmd1] different MD for read '%s': '%s' != '%s'\n", bam1_qname(b), old_md+1, str->s);
+               if (is_diff) {
+                       fprintf(stderr, "[bam_fillmd1] different MD for read '%s': '%s' -> '%s'\n", bam1_qname(b), old_md+1, str->s);
+                       bam_aux_del(b, old_md);
+                       bam_aux_append(b, "MD", 'Z', str->l + 1, (uint8_t*)str->s);
+               }
        }
        free(str->s); free(str);
 }
@@ -127,8 +132,11 @@ int bam_fillmd(int argc, char *argv[])
                                free(ref);
                                ref = fai_fetch(fai, fp->header->target_name[b->core.tid], &len);
                                tid = b->core.tid;
+                               if (ref == 0)
+                                       fprintf(stderr, "[bam_fillmd] fail to find sequence '%s' in the reference.\n",
+                                                       fp->header->target_name[tid]);
                        }
-                       bam_fillmd1(b, ref, is_equal);
+                       if (ref) bam_fillmd1(b, ref, is_equal);
                }
                samwrite(fpout, b);
        }