Make synced_bcf_reader work with 64 bit positions.

MAX_CSI_COOR is now about 12 Tbp, which is the limit for
a CSI index with min_shift 14.

Author: daviesrob

Diff for: bcf_sr_sort.c

@@ -288,7 +288,7 @@ void debug_vbuf(sr_sort_t *srt)
         for (i=0; i<srt->sr->nreaders; i++)
         {
             vcf_buf_t *buf = &srt->vcf_buf[i];
-            fprintf(stderr,"\t%d", buf->rec[j] ? buf->rec[j]->pos+1 : 0);
+            fprintf(stderr,"\t%"PRIhts_pos, buf->rec[j] ? buf->rec[j]->pos+1 : 0);
         }
         fprintf(stderr,"\n");
     }
@@ -330,7 +330,7 @@ int bcf_sr_sort_add_active(sr_sort_t *srt, int idx)
     srt->active[srt->nactive - 1] = idx;
     return 0; // FIXME: check for errs in this function
 }
-static int bcf_sr_sort_set(bcf_srs_t *readers, sr_sort_t *srt, const char *chr, int min_pos)
+static int bcf_sr_sort_set(bcf_srs_t *readers, sr_sort_t *srt, const char *chr, hts_pos_t min_pos)
 {
     if ( !srt->grp_str2int )
     {
@@ -556,7 +556,7 @@ static int bcf_sr_sort_set(bcf_srs_t *readers, sr_sort_t *srt, const char *chr,
     return 0;  // FIXME: check for errs in this function
 }
 
-int bcf_sr_sort_next(bcf_srs_t *readers, sr_sort_t *srt, const char *chr, int min_pos)
+int bcf_sr_sort_next(bcf_srs_t *readers, sr_sort_t *srt, const char *chr, hts_pos_t min_pos)
 {
     int i,j;
     assert( srt->nactive>0 );

Diff for: bcf_sr_sort.h

@@ -90,15 +90,16 @@ typedef struct
     int moff, noff, *off, mcharp;
     char **charp;
     const char *chr;
-    int pos, nsr, msr;
+    hts_pos_t pos;
+    int nsr, msr;
     int pair;
     int nactive, mactive, *active;  // list of readers with lines at the current pos
 }
 sr_sort_t;
 
 sr_sort_t *bcf_sr_sort_init(sr_sort_t *srt);
 void bcf_sr_sort_reset(sr_sort_t *srt);
-int bcf_sr_sort_next(bcf_srs_t *readers, sr_sort_t *srt, const char *chr, int pos);
+int bcf_sr_sort_next(bcf_srs_t *readers, sr_sort_t *srt, const char *chr, hts_pos_t pos);
 int bcf_sr_sort_set_active(sr_sort_t *srt, int i);
 int bcf_sr_sort_add_active(sr_sort_t *srt, int i);
 void bcf_sr_sort_destroy(sr_sort_t *srt);

Diff for: htslib/synced_bcf_reader.h

@@ -125,8 +125,9 @@ typedef struct _bcf_sr_regions_t
     char **seq_names;       // sequence names
     int nseqs;              // number of sequences (chromosomes) in the file
     int iseq;               // current position: chr name, index to snames
-    int start, end;         // current position: start, end of the region (0-based)
-    int prev_seq, prev_start;
+    hts_pos_t start, end;   // current position: start, end of the region (0-based)
+    int prev_seq;
+    hts_pos_t prev_start;
 }
 bcf_sr_regions_t;
 
@@ -241,7 +242,7 @@ int bcf_sr_next_line(bcf_srs_t *readers);
  *  @seq:  sequence name; NULL to seek to start
  *  @pos:  0-based coordinate
  */
-int bcf_sr_seek(bcf_srs_t *readers, const char *seq, int pos);
+int bcf_sr_seek(bcf_srs_t *readers, const char *seq, hts_pos_t pos);
 
 /**
  * bcf_sr_set_samples() - sets active samples
@@ -336,7 +337,7 @@ int bcf_sr_regions_next(bcf_sr_regions_t *reg);
  *  regions and more regions exist; -2 if not in the regions and there are no more
  *  regions left.
  */
-int bcf_sr_regions_overlap(bcf_sr_regions_t *reg, const char *seq, int start, int end);
+int bcf_sr_regions_overlap(bcf_sr_regions_t *reg, const char *seq, hts_pos_t start, hts_pos_t end);
 
 /*
  *  bcf_sr_regions_flush() - calls repeatedly regs->missed_reg_handler() until

Diff for: synced_bcf_reader.c

@@ -39,11 +39,15 @@ DEALINGS IN THE SOFTWARE.  */
 #include "htslib/thread_pool.h"
 #include "bcf_sr_sort.h"
 
-#define MAX_CSI_COOR 0x7fffffff     // maximum indexable coordinate of .csi
+// Maximum indexable coordinate of .csi, for default min_shift of 14.
+// This comes out to about 17 Tbp.  Limiting factor is the bin number,
+// which is a uint32_t in CSI.  The highest number of levels compatible
+// with this is 10 (needs 31 bits).
+#define MAX_CSI_COOR ((1LL << (14 + 30)) - 1)
 
 typedef struct
 {
-    uint32_t start, end;
+    hts_pos_t start, end;
 }
 region1_t;
 
@@ -61,7 +65,7 @@ typedef struct
 }
 aux_t;
 
-static int _regions_add(bcf_sr_regions_t *reg, const char *chr, int start, int end);
+static int _regions_add(bcf_sr_regions_t *reg, const char *chr, hts_pos_t start, hts_pos_t end);
 static bcf_sr_regions_t *_regions_init_string(const char *str);
 static int _regions_match_alleles(bcf_sr_regions_t *reg, int als_idx, bcf1_t *rec);
 
@@ -420,11 +424,11 @@ static inline int has_filter(bcf_sr_t *reader, bcf1_t *line)
     return 0;
 }
 
-static int _reader_seek(bcf_sr_t *reader, const char *seq, int start, int end)
+static int _reader_seek(bcf_sr_t *reader, const char *seq, hts_pos_t start, hts_pos_t end)
 {
     if ( end>=MAX_CSI_COOR )
     {
-        hts_log_error("The coordinate is out of csi index limit: %d", end+1);
+        hts_log_error("The coordinate is out of csi index limit: %"PRIhts_pos, end+1);
         exit(1);
     }
     if ( reader->itr )
@@ -446,7 +450,7 @@ static int _reader_seek(bcf_sr_t *reader, const char *seq, int start, int end)
         reader->itr = bcf_itr_queryi(reader->bcf_idx,tid,start,end+1);
     }
     if (!reader->itr) {
-        hts_log_error("Could not seek: %s:%d-%d", seq, start + 1, end + 1);
+        hts_log_error("Could not seek: %s:%"PRIhts_pos"-%"PRIhts_pos, seq, start + 1, end + 1);
         assert(0);
     }
     return 0;
@@ -581,7 +585,8 @@ static void _reader_shift_buffer(bcf_sr_t *reader)
 
 static int next_line(bcf_srs_t *files)
 {
-    int i, min_pos = INT_MAX;
+    int i;
+    hts_pos_t min_pos = HTS_POS_MAX;
     const char *chr = NULL;
 
     // Loop until next suitable line is found or all readers have finished
@@ -606,7 +611,7 @@ static int next_line(bcf_srs_t *files)
             else if ( min_pos==files->readers[i].buffer[1]->pos )
                 bcf_sr_sort_add_active(&BCF_SR_AUX(files)->sort, i);
         }
-        if ( min_pos==INT_MAX )
+        if ( min_pos==HTS_POS_MAX )
         {
             if ( !files->regions ) break;
             continue;
@@ -622,7 +627,7 @@ static int next_line(bcf_srs_t *files)
                 for (i=0; i<files->nreaders; i++)
                     if ( files->readers[i].nbuffer && files->readers[i].buffer[1]->pos==min_pos )
                         _reader_shift_buffer(&files->readers[i]);
-                min_pos = INT_MAX;
+                min_pos = HTS_POS_MAX;
                 chr = NULL;
                 continue;
             }
@@ -672,7 +677,7 @@ static void bcf_sr_seek_start(bcf_srs_t *readers)
 }
 
 
-int bcf_sr_seek(bcf_srs_t *readers, const char *seq, int pos)
+int bcf_sr_seek(bcf_srs_t *readers, const char *seq, hts_pos_t pos)
 {
     if ( !readers->regions ) return 0;
     bcf_sr_sort_reset(&BCF_SR_AUX(readers)->sort);
@@ -767,7 +772,7 @@ int bcf_sr_set_samples(bcf_srs_t *files, const char *fname, int is_file)
 
 // Add a new region into a list sorted by start,end. On input the coordinates
 // are 1-based, stored 0-based, inclusive.
-static int _regions_add(bcf_sr_regions_t *reg, const char *chr, int start, int end)
+static int _regions_add(bcf_sr_regions_t *reg, const char *chr, hts_pos_t start, hts_pos_t end)
 {
     if ( start==-1 && end==-1 )
     {
@@ -828,7 +833,7 @@ static bcf_sr_regions_t *_regions_init_string(const char *str)
 
     kstring_t tmp = {0,0,0};
     const char *sp = str, *ep = str;
-    int from, to;
+    hts_pos_t from, to;
     while ( 1 )
     {
         while ( *ep && *ep!=',' && *ep!=':' ) ep++;
@@ -880,7 +885,7 @@ static bcf_sr_regions_t *_regions_init_string(const char *str)
 
 // ichr,ifrom,ito are 0-based;
 // returns -1 on error, 0 if the line is a comment line, 1 on success
-static int _regions_parse_line(char *line, int ichr,int ifrom,int ito, char **chr,char **chr_end,int *from,int *to)
+static int _regions_parse_line(char *line, int ichr, int ifrom, int ito, char **chr, char **chr_end, hts_pos_t *from, hts_pos_t *to)
 {
     if (ifrom < 0 || ito < 0) return -1;
     *chr_end = NULL;
@@ -970,7 +975,8 @@ bcf_sr_regions_t *bcf_sr_regions_init(const char *regions, int is_file, int ichr
         while ( hts_getline(reg->file, KS_SEP_LINE, &reg->line) > 0 )
         {
             char *chr, *chr_end;
-            int from, to, ret;
+            hts_pos_t from, to;
+            int ret;
             ret = _regions_parse_line(reg->line.s, ichr,ifrom,abs(ito), &chr,&chr_end,&from,&to);
             if ( ret < 0 )
             {
@@ -1077,7 +1083,8 @@ int bcf_sr_regions_next(bcf_sr_regions_t *reg)
 
     // reading from tabix
     char *chr, *chr_end;
-    int ichr = 0, ifrom = 1, ito = 2, is_bed = 0, from, to;
+    int ichr = 0, ifrom = 1, ito = 2, is_bed = 0;
+    hts_pos_t from, to;
     if ( reg->tbx )
     {
         ichr   = reg->tbx->conf.sc-1;
@@ -1196,7 +1203,7 @@ static int _regions_match_alleles(bcf_sr_regions_t *reg, int als_idx, bcf1_t *re
     return !(type & VCF_INDEL) ? 1 : 0;
 }
 
-int bcf_sr_regions_overlap(bcf_sr_regions_t *reg, const char *seq, int start, int end)
+int bcf_sr_regions_overlap(bcf_sr_regions_t *reg, const char *seq, hts_pos_t start, hts_pos_t end)
 {
     int iseq;
     if ( khash_str2int_get(reg->seq_hash, seq, &iseq)<0 ) return -1;    // no such sequence