Fixed a bug in knet_seek

[samtools.git] / bam_index.c

diff --git a/bam_index.c b/bam_index.c
index 329df323851f7bd61ecbbc345790f122b4303b6f..853ac705fc9bef5a653ab9aae9eb9100ec551e12 100644 (file)
--- a/bam_index.c
+++ b/bam_index.c
@@ -1,9 +1,12 @@
 #include <ctype.h>
+#include <assert.h>
 #include "bam.h"
 #include "khash.h"
 #include "ksort.h"
 #include "bam_endian.h"
+#ifdef _USE_KNETFILE
 #include "knetfile.h"
+#endif
 
 /*!
   @header
@@ -324,11 +327,10 @@ static bam_index_t *bam_index_load_core(FILE *fp)
 
 bam_index_t *bam_index_load_local(const char *_fn)
 {
-       bam_index_t *idx;
        FILE *fp;
        char *fnidx, *fn;
 
-       if (strstr(_fn, "ftp://") == _fn) {
+       if (strstr(_fn, "ftp://") == _fn || strstr(_fn, "http://") == _fn) {
                const char *p;
                int l = strlen(_fn);
                for (p = _fn + l - 1; p >= _fn; --p)
@@ -337,21 +339,24 @@ bam_index_t *bam_index_load_local(const char *_fn)
        } else fn = strdup(_fn);
        fnidx = (char*)calloc(strlen(fn) + 5, 1);
        strcpy(fnidx, fn); strcat(fnidx, ".bai");
-       fp = fopen(fnidx, "r");
+       fp = fopen(fnidx, "rb");
        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");
+                       fp = fopen(fnidx, "rb");
                }
        }
        free(fnidx); free(fn);
-       idx = bam_index_load_core(fp);
-       fclose(fp);
-       return idx;
+       if (fp) {
+               bam_index_t *idx = bam_index_load_core(fp);
+               fclose(fp);
+               return idx;
+       } else return 0;
 }
 
+#ifdef _USE_KNETFILE
 static void download_from_remote(const char *url)
 {
        const int buf_size = 1 * 1024 * 1024;
@@ -360,7 +365,7 @@ static void download_from_remote(const char *url)
        uint8_t *buf;
        knetFile *fp_remote;
        int l;
-       if (strstr(url, "ftp://") != url) return;
+       if (strstr(url, "ftp://") != url && strstr(url, "http://") != url) return;
        l = strlen(url);
        for (fn = (char*)url + l - 1; fn >= url; --fn)
                if (*fn == '/') break;
@@ -370,7 +375,7 @@ static void download_from_remote(const char *url)
                fprintf(stderr, "[download_from_remote] fail to open remote file.\n");
                return;
        }
-       if ((fp = fopen(fn, "w")) == 0) {
+       if ((fp = fopen(fn, "wb")) == 0) {
                fprintf(stderr, "[download_from_remote] fail to create file in the working directory.\n");
                knet_close(fp_remote);
                return;
@@ -382,18 +387,25 @@ static void download_from_remote(const char *url)
        fclose(fp);
        knet_close(fp_remote);
 }
+#else
+static void download_from_remote(const char *url)
+{
+       return;
+}
+#endif
 
 bam_index_t *bam_index_load(const char *fn)
 {
        bam_index_t *idx;
        idx = bam_index_load_local(fn);
-       if (idx == 0 && strstr(fn, "ftp://") == fn) {
+       if (idx == 0 && (strstr(fn, "ftp://") == fn || strstr(fn, "http://") == fn)) {
                char *fnidx = calloc(strlen(fn) + 5, 1);
                strcat(strcpy(fnidx, fn), ".bai");
                fprintf(stderr, "[bam_index_load] attempting to download the remote index file.\n");
                download_from_remote(fnidx);
                idx = bam_index_load_local(fn);
        }
+       if (idx == 0) fprintf(stderr, "[bam_index_load] fail to load BAM index.\n");
        return idx;
 }
 
@@ -403,18 +415,21 @@ int bam_index_build2(const char *fn, const char *_fnidx)
        FILE *fpidx;
        bamFile fp;
        bam_index_t *idx;
-       assert(fp = bam_open(fn, "r"));
+       if ((fp = bam_open(fn, "r")) == 0) {
+               fprintf(stderr, "[bam_index_build2] fail to open the BAM file.\n");
+               return -1;
+       }
        idx = bam_index_core(fp);
        bam_close(fp);
        if (_fnidx == 0) {
                fnidx = (char*)calloc(strlen(fn) + 5, 1);
                strcpy(fnidx, fn); strcat(fnidx, ".bai");
        } else fnidx = strdup(_fnidx);
-       fpidx = fopen(fnidx, "w");
+       fpidx = fopen(fnidx, "wb");
        if (fpidx == 0) {
                fprintf(stderr, "[bam_index_build2] fail to create the index file.\n");
                free(fnidx);
-               return 1;
+               return -1;
        }
        bam_index_save(idx, fpidx);
        bam_index_destroy(idx);
@@ -461,7 +476,8 @@ static inline int is_overlap(uint32_t beg, uint32_t end, const bam1_t *b)
        return (rend > beg && rbeg < end);
 }
 
-int bam_fetch(bamFile fp, const bam_index_t *idx, int tid, int beg, int end, void *data, bam_fetch_f func)
+// bam_fetch helper function retrieves 
+pair64_t * get_chunk_coordinates(const bam_index_t *idx, int tid, int beg, int end, int* cnt_off)
 {
        uint16_t *bins;
        int i, n_bins, n_off;
@@ -492,10 +508,8 @@ int bam_fetch(bamFile fp, const bam_index_t *idx, int tid, int beg, int end, voi
        }
        free(bins);
        {
-               bam1_t *b;
-               int l, ret, n_seeks;
-               uint64_t curr_off;
-               b = (bam1_t*)calloc(1, sizeof(bam1_t));
+               bam1_t *b = (bam1_t*)calloc(1, sizeof(bam1_t));
+               int l;
                ks_introsort(off, n_off, off);
                // resolve completely contained adjacent blocks
                for (i = 1, l = 0; i < n_off; ++i)
@@ -518,8 +532,23 @@ int bam_fetch(bamFile fp, const bam_index_t *idx, int tid, int beg, int end, voi
                        n_off = l + 1;
 #endif
                }
+               bam_destroy1(b);
+       }
+       *cnt_off = n_off;
+       return off;
+}
+
+int bam_fetch(bamFile fp, const bam_index_t *idx, int tid, int beg, int end, void *data, bam_fetch_f func)
+{
+       int n_off;
+       pair64_t *off = get_chunk_coordinates(idx, tid, beg, end, &n_off);
+       if (off == 0) return 0;
+       {
                // retrive alignments
+               uint64_t curr_off;
+               int i, ret, n_seeks;
                n_seeks = 0; i = -1; curr_off = 0;
+               bam1_t *b = (bam1_t*)calloc(1, sizeof(bam1_t));
                for (;;) {
                        if (curr_off == 0 || curr_off >= off[i].v) { // then jump to the next chunk
                                if (i == n_off - 1) break; // no more chunks