* no changes to samtools itself

[samtools.git] / bam_index.c

diff --git a/bam_index.c b/bam_index.c
index 6bc735eeaac4fbbb4318eabd6eb8c07c443465e5..4a41f5264a67148c75f1f9be6d7e6880e5b848a6 100644 (file)
--- a/bam_index.c
+++ b/bam_index.c
@@ -158,7 +158,7 @@ bam_index_t *bam_index_core(bamFile fp)
                        fprintf(stderr, "[bam_index_core] the alignment is not sorted. Abort!\n");
                        exit(1);
                }
-               if (b->core.bin < 4681) insert_offset2(&idx->index2[b->core.tid], b, last_off);
+               if (b->core.tid >= 0 && b->core.bin < 4681) insert_offset2(&idx->index2[b->core.tid], b, last_off);
                if (c->bin != last_bin) { // then possibly write the binning index
                        if (save_bin != 0xffffffffu) // save_bin==0xffffffffu only happens to the first record
                                insert_offset(idx->index[save_tid], save_bin, save_off, last_off);
@@ -268,7 +268,16 @@ bam_index_t *bam_index_load(const char *fn)
 
        fnidx = (char*)calloc(strlen(fn) + 5, 1);
        strcpy(fnidx, fn); strcat(fnidx, ".bai");
-       if ((fp = fopen(fnidx, "r")) == 0) {
+       fp = fopen(fnidx, "r");
+       if (fp == 0) { // try "{base}.bai"
+               char *s = strstr(fn, "bam");
+               if (s == fn + strlen(fn) - 3) {
+                       strcpy(fnidx, fn);
+                       fnidx[strlen(fn)-1] = 'i';
+                       fp = fopen(fnidx, "r");
+               }
+       }
+       if (fp == 0) {
                fprintf(stderr, "[bam_index_load] the alignment is not indexed. Please run `index' command first. Abort!\n");
                exit(1);
        }
@@ -327,7 +336,7 @@ bam_index_t *bam_index_load(const char *fn)
        return idx;
 }
 
-int bam_index_build(const char *fn)
+int bam_index_build2(const char *fn, const char *_fnidx)
 {
        char *fnidx;
        FILE *fpidx;
@@ -336,9 +345,16 @@ int bam_index_build(const char *fn)
        assert(fp = bam_open(fn, "r"));
        idx = bam_index_core(fp);
        bam_close(fp);
-       fnidx = (char*)calloc(strlen(fn) + 5, 1);
-       strcpy(fnidx, fn); strcat(fnidx, ".bai");
-       assert(fpidx = fopen(fnidx, "w"));
+       if (_fnidx == 0) {
+               fnidx = (char*)calloc(strlen(fn) + 5, 1);
+               strcpy(fnidx, fn); strcat(fnidx, ".bai");
+       } else fnidx = strdup(_fnidx);
+       fpidx = fopen(fnidx, "w");
+       if (fpidx == 0) {
+               fprintf(stderr, "[bam_index_build2] fail to create the index file.\n");
+               free(fnidx);
+               return 1;
+       }
        bam_index_save(idx, fpidx);
        bam_index_destroy(idx);
        fclose(fpidx);
@@ -346,13 +362,19 @@ int bam_index_build(const char *fn)
        return 0;
 }
 
+int bam_index_build(const char *fn)
+{
+       return bam_index_build2(fn, 0);
+}
+
 int bam_index(int argc, char *argv[])
 {
        if (argc < 2) {
-               fprintf(stderr, "Usage: samtools index <in.bam>\n");
+               fprintf(stderr, "Usage: samtools index <in.bam> [<out.index>]\n");
                return 1;
        }
-       bam_index_build(argv[1]);
+       if (argc >= 3) bam_index_build2(argv[1], argv[2]);
+       else bam_index_build(argv[1]);
        return 0;
 }
 
@@ -374,7 +396,7 @@ static inline int reg2bins(uint32_t beg, uint32_t end, uint16_t list[MAX_BIN])
 static inline int is_overlap(uint32_t beg, uint32_t end, const bam1_t *b)
 {
        uint32_t rbeg = b->core.pos;
-       uint32_t rend = bam_calend(&b->core, bam1_cigar(b));
+       uint32_t rend = b->core.n_cigar? bam_calend(&b->core, bam1_cigar(b)) : b->core.pos + 1;
        return (rend > beg && rbeg < end);
 }
 
@@ -410,16 +432,20 @@ int bam_fetch(bamFile fp, const bam_index_t *idx, int tid, int beg, int end, voi
        free(bins);
        {
                bam1_t *b;
-               int ret, n_seeks;
+               int l, ret, n_seeks;
                uint64_t curr_off;
                b = (bam1_t*)calloc(1, sizeof(bam1_t));
                ks_introsort(off, n_off, off);
-               // resolve overlaps between adjecent blocks; this may happen due to the merge in indexing
+               // resolve completely contained adjacent blocks
+               for (i = 1, l = 0; i < n_off; ++i)
+                       if (off[l].v < off[i].v)
+                               off[++l] = off[i];
+               n_off = l + 1;
+               // resolve overlaps between adjacent blocks; this may happen due to the merge in indexing
                for (i = 1; i < n_off; ++i)
                        if (off[i-1].v >= off[i].u) off[i-1].v = off[i].u;
                { // merge adjacent blocks
 #if defined(BAM_TRUE_OFFSET) || defined(BAM_VIRTUAL_OFFSET16)
-                       int l;
                        for (i = 1, l = 0; i < n_off; ++i) {
 #ifdef BAM_TRUE_OFFSET
                                if (off[l].v + BAM_MIN_CHUNK_GAP > off[i].u) off[l].v = off[i].v;