.appveyor.yml

@@ -5,10 +5,6 @@ version: 'vers.{build}'
 
 # branches to build
 branches:
-    # Whitelist
-    only:
-      - develop
-
     # Blacklist
     except:
       - gh-pages

Makefile

@@ -175,8 +175,9 @@ cram_sam_header_h = cram/sam_header.h cram/string_alloc.h cram/pooled_alloc.h $(
 cram_samtools_h = cram/cram_samtools.h $(htslib_sam_h) $(cram_sam_header_h)
 cram_structs_h = cram/cram_structs.h $(htslib_thread_pool_h) cram/string_alloc.h $(htslib_khash_h)
 cram_open_trace_file_h = cram/open_trace_file.h cram/mFILE.h
-hfile_internal_h = hfile_internal.h $(htslib_hfile_h)
-hts_internal_h = hts_internal.h $(htslib_hts_h)
+hfile_internal_h = hfile_internal.h $(htslib_hfile_h) $(textutils_internal_h)
+hts_internal_h = hts_internal.h $(htslib_hts_h) $(textutils_internal_h)
+textutils_internal_h = textutils_internal.h $(htslib_kstring_h)
 thread_pool_internal_h = thread_pool_internal.h $(htslib_thread_pool_h)
 
 
@@ -195,6 +196,7 @@ config.h:
 	echo '/* Default config.h generated by Makefile */' > $@
 	echo '#define HAVE_LIBBZ2 1' >> $@
 	echo '#define HAVE_LIBLZMA 1' >> $@
+	echo '#define HAVE_LZMA_H 1' >> $@
 	echo '#define HAVE_FSEEKO 1' >> $@
 	echo '#define HAVE_DRAND48 1' >> $@
 
@@ -293,8 +295,9 @@ hfile.o hfile.pico: hfile.c config.h $(htslib_hfile_h) $(hfile_internal_h) $(hts
 hfile_gcs.o hfile_gcs.pico: hfile_gcs.c config.h $(htslib_hts_h) $(htslib_kstring_h) $(hfile_internal_h)
 hfile_libcurl.o hfile_libcurl.pico: hfile_libcurl.c config.h $(hfile_internal_h) $(htslib_hts_h) $(htslib_kstring_h)
 hfile_net.o hfile_net.pico: hfile_net.c config.h $(hfile_internal_h) $(htslib_knetfile_h)
-hfile_s3.o hfile_s3.pico: hfile_s3.c config.h $(hts_internal_h) $(hfile_internal_h) $(htslib_hts_h) $(htslib_kstring_h)
+hfile_s3.o hfile_s3.pico: hfile_s3.c config.h $(hfile_internal_h) $(htslib_hts_h) $(htslib_kstring_h)
 hts.o hts.pico: hts.c config.h $(htslib_hts_h) $(htslib_bgzf_h) $(cram_h) $(hfile_internal_h) $(htslib_hfile_h) version.h $(hts_internal_h) $(htslib_khash_h) $(htslib_kseq_h) $(htslib_ksort_h)
+hts_os.o hts_os.pico: hts_os.c config.h os/rand.c
 vcf.o vcf.pico: vcf.c config.h $(htslib_vcf_h) $(htslib_bgzf_h) $(htslib_tbx_h) $(htslib_hfile_h) $(hts_internal_h) $(htslib_khash_str2int_h) $(htslib_kstring_h) $(htslib_khash_h) $(htslib_kseq_h) $(htslib_hts_endian_h)
 sam.o sam.pico: sam.c config.h $(htslib_sam_h) $(htslib_bgzf_h) $(cram_h) $(hts_internal_h) $(htslib_hfile_h) $(htslib_khash_h) $(htslib_kseq_h) $(htslib_kstring_h) $(htslib_hts_endian_h)
 tbx.o tbx.pico: tbx.c config.h $(htslib_tbx_h) $(htslib_bgzf_h) $(hts_internal_h) $(htslib_khash_h)
@@ -317,7 +320,7 @@ cram/cram_decode.o cram/cram_decode.pico: cram/cram_decode.c config.h $(cram_h)
 cram/cram_encode.o cram/cram_encode.pico: cram/cram_encode.c config.h $(cram_h) $(cram_os_h) $(htslib_hts_h) $(htslib_hts_endian_h)
 cram/cram_external.o cram/cram_external.pico: cram/cram_external.c config.h $(htslib_hfile_h) $(cram_h)
 cram/cram_index.o cram/cram_index.pico: cram/cram_index.c config.h $(htslib_bgzf_h) $(htslib_hfile_h) $(hts_internal_h) $(cram_h) $(cram_os_h)
-cram/cram_io.o cram/cram_io.pico: cram/cram_io.c config.h $(cram_h) $(cram_os_h) $(htslib_hts_h) $(cram_open_trace_file_h) cram/rANS_static.h $(htslib_hfile_h) $(htslib_bgzf_h) $(htslib_faidx_h) $(hts_internal_h)
+cram/cram_io.o cram/cram_io.pico: cram/cram_io.c config.h os/lzma_stub.h $(cram_h) $(cram_os_h) $(htslib_hts_h) $(cram_open_trace_file_h) cram/rANS_static.h $(htslib_hfile_h) $(htslib_bgzf_h) $(htslib_faidx_h) $(hts_internal_h)
 cram/cram_samtools.o cram/cram_samtools.pico: cram/cram_samtools.c config.h $(cram_h) $(htslib_sam_h)
 cram/cram_stats.o cram/cram_stats.pico: cram/cram_stats.c config.h $(cram_h) $(cram_os_h)
 cram/files.o cram/files.pico: cram/files.c config.h $(cram_misc_h)

NEWS

+Noteworthy changes in release 1.7 (26th January 2018)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+* BAM: HTSlib now supports BAMs which include CIGARs with more than 
+  65535 operations as per HTS-Specs 18th November (dab57f4 and 2f915a8).
+
+* BCF/VCF:
+  - Removed the need for long double in pileup calculations.
+  - Sped up the synced reader in some situations.
+  - Bug fixing: removed memory leak in bcf_copy.
+
+* CRAM:
+  - Added support for HTS_IDX_START in cram iterators.
+  - Easier to build when lzma header files are absent.
+  - Bug fixing: a region query with REQUIRED_FIELDS option to
+    disable sequence retrieval now gives correct results.
+  - Bug fixing: stop queries to regions starting after the last
+    read on a chromosome from incorrectly reporting errors
+    (#651, #653; reported by Imran Haque and @egafni via pysam).
+
+* Multi-region iterator: The new structure takes a list of regions and
+  iterates over all, deduplicating reads in the process, and producing a 
+  full list of file offset intervals. This is usually much faster than 
+  repeatedly using the old single-region iterator on a series of regions.
+
+* Curl improvements:
+  - Add Bearer token support via HTS_AUTH_LOCATION env (#600).
+  - Use CURL_CA_BUNDLE environment variable to override the CA (#622; 
+    thanks to Garret Kelly & David Alexander).
+  - Speed up (removal of excessive waiting) for both http(s) and ftp.
+  - Avoid repeatedly reconnecting by removal of unnecessary seeks.
+  - Bug fixing: double free when libcurl_open fails.
+
+* BGZF block caching, if enabled, now performs far better (#629; reported
+  by Ram Yalamanchili).
+
+* Added an hFILE layer for in-memory I/O buffers (#590; thanks to Thomas 
+  Hickman).
+
+* Tidied up the drand48 support (intended for systems that do not
+  provide this function).
+
 Noteworthy changes in release 1.6 (28th September 2017)
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 

bcf_sr_sort.c

@@ -327,6 +327,20 @@ char *grp_create_key(sr_sort_t *srt)
     }
     return ret;
 }
+int bcf_sr_sort_set_active(sr_sort_t *srt, int idx)
+{
+    hts_expand(int,idx+1,srt->mactive,srt->active);
+    srt->nactive = 1;
+    srt->active[srt->nactive - 1] = idx;
+    return 0;
+}
+int bcf_sr_sort_add_active(sr_sort_t *srt, int idx)
+{
+    hts_expand(int,idx+1,srt->mactive,srt->active);
+    srt->nactive++;
+    srt->active[srt->nactive - 1] = idx;
+    return 0;
+}
 static void bcf_sr_sort_set(bcf_srs_t *readers, sr_sort_t *srt, const char *chr, int min_pos)
 {
     if ( !srt->grp_str2int )
@@ -357,11 +371,11 @@ static void bcf_sr_sort_set(bcf_srs_t *readers, sr_sort_t *srt, const char *chr,
     memset(&grp,0,sizeof(grp_t));
 
     // group VCFs into groups, each with a unique combination of variants in the duplicate lines
-    int ireader,ivar,irec,igrp,ivset;
-    for (ireader=0; ireader<readers->nreaders; ireader++)
+    int ireader,ivar,irec,igrp,ivset,iact;
+    for (ireader=0; ireader<readers->nreaders; ireader++) srt->vcf_buf[ireader].nrec = 0;
+    for (iact=0; iact<srt->nactive; iact++)
     {
-        srt->vcf_buf[ireader].nrec = 0;
-
+        ireader = srt->active[iact];
         bcf_sr_t *reader = &readers->readers[ireader];
         int rid   = bcf_hdr_name2id(reader->header, chr);
         grp.nvar  = 0;
@@ -546,23 +560,8 @@ static void bcf_sr_sort_set(bcf_srs_t *readers, sr_sort_t *srt, const char *chr,
 int bcf_sr_sort_next(bcf_srs_t *readers, sr_sort_t *srt, const char *chr, int min_pos)
 {
     int i,j;
-    if ( readers->nreaders == 1 )
-    {
-        srt->nsr = 1;
-        if ( !srt->msr )
-        {
-            srt->vcf_buf = (vcf_buf_t*) calloc(1,sizeof(vcf_buf_t));   // first time here
-            srt->msr = 1;
-        }
-        bcf_sr_t *reader = &readers->readers[0];
-        assert( reader->buffer[1]->pos==min_pos );
-        bcf1_t *tmp = reader->buffer[0];
-        for (j=1; j<=reader->nbuffer; j++) reader->buffer[j-1] = reader->buffer[j];
-        reader->buffer[ reader->nbuffer ] = tmp;
-        reader->nbuffer--;
-        readers->has_line[0] = 1;
-        return 1;
-    }
+    assert( srt->nactive>0 );
+
     if ( srt->nsr != readers->nreaders )
     {
         srt->sr = readers;
@@ -575,6 +574,19 @@ int bcf_sr_sort_next(bcf_srs_t *readers, sr_sort_t *srt, const char *chr, int mi
         srt->nsr = readers->nreaders;
         srt->chr = NULL;
     }
+    if ( srt->nactive == 1 )
+    {
+        if ( readers->nreaders>1 )
+            memset(readers->has_line, 0, readers->nreaders*sizeof(*readers->has_line));
+        bcf_sr_t *reader = &readers->readers[srt->active[0]];
+        assert( reader->buffer[1]->pos==min_pos );
+        bcf1_t *tmp = reader->buffer[0];
+        for (j=1; j<=reader->nbuffer; j++) reader->buffer[j-1] = reader->buffer[j];
+        reader->buffer[ reader->nbuffer ] = tmp;
+        reader->nbuffer--;
+        readers->has_line[srt->active[0]] = 1;
+        return 1;
+    }
     if ( !srt->chr || srt->pos!=min_pos || strcmp(srt->chr,chr) ) bcf_sr_sort_set(readers, srt, chr, min_pos);
 
     if ( !srt->vcf_buf[0].nrec ) return 0;
@@ -629,6 +641,7 @@ sr_sort_t *bcf_sr_sort_init(sr_sort_t *srt)
 }
 void bcf_sr_sort_destroy(sr_sort_t *srt)
 {
+    free(srt->active);
     if ( srt->var_str2int ) khash_str2int_destroy_free(srt->var_str2int);
     if ( srt->grp_str2int ) khash_str2int_destroy_free(srt->grp_str2int);
     int i;

bcf_sr_sort.h

@@ -92,11 +92,14 @@ typedef struct
     const char *chr;
     int pos, 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);
 int bcf_sr_sort_next(bcf_srs_t *readers, sr_sort_t *srt, const char *chr, int 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);
 void bcf_sr_sort_remove_reader(bcf_srs_t *readers, sr_sort_t *srt, int i);
 

bgzf.c

@@ -71,10 +71,16 @@ typedef struct {
     uint8_t *block;
     int64_t end_offset;
 } cache_t;
+
 #include "htslib/khash.h"
 KHASH_MAP_INIT_INT64(cache, cache_t)
 #endif
 
+struct bgzf_cache_t {
+    khash_t(cache) *h;
+    khint_t last_pos;
+};
+
 #ifdef BGZF_MT
 
 typedef struct bgzf_job {
@@ -215,7 +221,16 @@ static BGZF *bgzf_read_init(hFILE *hfpr)
     fp->is_compressed = (n==18 && magic[0]==0x1f && magic[1]==0x8b);
     fp->is_gzip = ( !fp->is_compressed || ((magic[3]&4) && memcmp(&magic[12], "BC\2\0",4)==0) ) ? 0 : 1;
 #ifdef BGZF_CACHE
-    fp->cache = kh_init(cache);
+    if (!(fp->cache = malloc(sizeof(*fp->cache)))) {
+        free(fp);
+        return NULL;
+    }
+    if (!(fp->cache->h = kh_init(cache))) {
+        free(fp->cache);
+        free(fp);
+        return NULL;
+    }
+    fp->cache->last_pos = 0;
 #endif
     return fp;
 }
@@ -524,11 +539,12 @@ static int check_header(const uint8_t *header)
 static void free_cache(BGZF *fp)
 {
     khint_t k;
-    khash_t(cache) *h = (khash_t(cache)*)fp->cache;
     if (fp->is_write) return;
+    khash_t(cache) *h = fp->cache->h;
     for (k = kh_begin(h); k < kh_end(h); ++k)
         if (kh_exist(h, k)) free(kh_val(h, k).block);
     kh_destroy(cache, h);
+    free(fp->cache);
 }
 
 static int load_block_from_cache(BGZF *fp, int64_t block_address)
@@ -536,15 +552,14 @@ static int load_block_from_cache(BGZF *fp, int64_t block_address)
     khint_t k;
     cache_t *p;
 
-    khash_t(cache) *h = (khash_t(cache)*)fp->cache;
+    khash_t(cache) *h = fp->cache->h;
     k = kh_get(cache, h, block_address);
     if (k == kh_end(h)) return 0;
     p = &kh_val(h, k);
     if (fp->block_length != 0) fp->block_offset = 0;
     fp->block_address = block_address;
     fp->block_length = p->size;
-    // FIXME: why BGZF_MAX_BLOCK_SIZE and not p->size?
-    memcpy(fp->uncompressed_block, p->block, BGZF_MAX_BLOCK_SIZE);
+    memcpy(fp->uncompressed_block, p->block, p->size);
     if ( hseek(fp->fp, p->end_offset, SEEK_SET) < 0 )
     {
         // todo: move the error up
@@ -557,29 +572,48 @@ static int load_block_from_cache(BGZF *fp, int64_t block_address)
 static void cache_block(BGZF *fp, int size)
 {
     int ret;
-    khint_t k;
+    khint_t k, k_orig;
+    uint8_t *block = NULL;
     cache_t *p;
     //fprintf(stderr, "Cache block at %llx\n", (int)fp->block_address);
-    khash_t(cache) *h = (khash_t(cache)*)fp->cache;
+    khash_t(cache) *h = fp->cache->h;
     if (BGZF_MAX_BLOCK_SIZE >= fp->cache_size) return;
+    if (fp->block_length < 0 || fp->block_length > BGZF_MAX_BLOCK_SIZE) return;
     if ((kh_size(h) + 1) * BGZF_MAX_BLOCK_SIZE > (uint32_t)fp->cache_size) {
-        /* A better way would be to remove the oldest block in the
-         * cache, but here we remove a random one for simplicity. This
-         * should not have a big impact on performance. */
-        for (k = kh_begin(h); k < kh_end(h); ++k)
-            if (kh_exist(h, k)) break;
-        if (k < kh_end(h)) {
-            free(kh_val(h, k).block);
+        /* Remove uniformly from any position in the hash by a simple
+         * round-robin approach.  An alternative strategy would be to
+         * remove the least recently accessed block, but the round-robin
+         * removal is simpler and is not expected to have a big impact
+         * on performance */
+        if (fp->cache->last_pos >= kh_end(h)) fp->cache->last_pos = kh_begin(h);
+        k_orig = k = fp->cache->last_pos;
+        if (++k >= kh_end(h)) k = kh_begin(h);
+        while (k != k_orig) {
+            if (kh_exist(h, k))
+                break;
+            if (++k == kh_end(h))
+                k = kh_begin(h);
+        }
+        fp->cache->last_pos = k;
+
+        if (k != k_orig) {
+            block = kh_val(h, k).block;
             kh_del(cache, h, k);
         }
+    } else {
+        block = (uint8_t*)malloc(BGZF_MAX_BLOCK_SIZE);
     }
+    if (!block) return;
     k = kh_put(cache, h, fp->block_address, &ret);
-    if (ret == 0) return; // if this happens, a bug!
+    if (ret <= 0) { // kh_put failed, or in there already (shouldn't happen)
+        free(block);
+        return;
+    }
     p = &kh_val(h, k);
     p->size = fp->block_length;
     p->end_offset = fp->block_address + size;
-    p->block = (uint8_t*)malloc(BGZF_MAX_BLOCK_SIZE);
-    memcpy(kh_val(h, k).block, fp->uncompressed_block, BGZF_MAX_BLOCK_SIZE);
+    p->block = block;
+    memcpy(p->block, fp->uncompressed_block, p->size);
 }
 #else
 static void free_cache(BGZF *fp) {}
@@ -1489,7 +1523,7 @@ int bgzf_close(BGZF* fp)
 
 void bgzf_set_cache_size(BGZF *fp, int cache_size)
 {
-    if (fp) fp->cache_size = cache_size;
+    if (fp && fp->cache) fp->cache_size = cache_size;
 }
 
 int bgzf_check_EOF(BGZF *fp) {

bgzip.c

 /* bgzip.c -- Block compression/decompression utility.
 
    Copyright (C) 2008, 2009 Broad Institute / Massachusetts Institute of Technology
-   Copyright (C) 2010, 2013-2017 Genome Research Ltd.
+   Copyright (C) 2010, 2013-2018 Genome Research Ltd.
 
    Permission is hereby granted, free of charge, to any person obtaining a copy
    of this software and associated documentation files (the "Software"), to deal
@@ -130,7 +130,7 @@ int main(int argc, char **argv)
         case 1:
             printf(
 "bgzip (htslib) %s\n"
-"Copyright (C) 2017 Genome Research Ltd.\n", hts_version());
+"Copyright (C) 2018 Genome Research Ltd.\n", hts_version());
             return EXIT_SUCCESS;
         case 'h':
         case '?': return bgzip_main_usage();

configure.ac

@@ -29,6 +29,8 @@ AC_PREREQ(2.63)  dnl This version introduced 4-argument AC_CHECK_HEADER
 AC_CONFIG_SRCDIR(hts.c)
 AC_CONFIG_HEADERS(config.h)
 
+m4_include([m4/hts_prog_cc_warnings.m4])
+
 dnl Copyright notice to be copied into the generated configure script
 AC_COPYRIGHT([Portions copyright (C) 2016 Genome Research Ltd.
 
@@ -49,6 +51,12 @@ AH_TOP([/* If you use configure, this file provides @%:@defines reflecting your
 AC_PROG_CC
 AC_PROG_RANLIB
 
+dnl Turn on compiler warnings, if possible
+HTS_PROG_CC_WARNINGS
+dnl Flags to treat warnings as errors.  These need to be applied to CFLAGS
+dnl later as they can interfere with some of the tests (notably AC_SEARCH_LIBS)
+HTS_PROG_CC_WERROR(hts_late_cflags)
+
 dnl Avoid chicken-and-egg problem where pkg-config supplies the
 dnl PKG_PROG_PKG_CONFIG macro, but we want to use it to check
 dnl for pkg-config...
@@ -229,9 +237,9 @@ fi
 
 if test "$enable_lzma" != no; then
   lzma_devel=ok
-  AC_CHECK_HEADER([lzma.h], [], [lzma_devel=missing], [;])
+  AC_CHECK_HEADERS([lzma.h], [], [lzma_devel=header-missing], [;])
   AC_CHECK_LIB([lzma], [lzma_easy_buffer_encode], [], [lzma_devel=missing])
-  if test $lzma_devel != ok; then
+  if test $lzma_devel = missing; then
     AC_MSG_ERROR([liblzma development files not found
 
 The CRAM format may use LZMA2 compression, which is implemented in HTSlib
@@ -354,6 +362,10 @@ AC_CHECK_DECL(PTHREAD_MUTEX_RECURSIVE, [], [AC_DEFINE([_XOPEN_SOURCE],[600], [Ne
 if test "$s3" = enabled ; then
    AC_DEFINE([ENABLE_S3], 1, [Define if HTSlib should enable S3 support.])
 fi
+
+dnl Apply value from HTS_PROG_CC_WERROR (if set)
+AS_IF([test "x$hts_late_cflags" != x],[CFLAGS="$CFLAGS $hts_late_cflags"])
+
 AC_SUBST([s3])
 AC_SUBST([CRYPTO_LIBS])
 

cram/cram_decode.c

@@ -1798,6 +1798,9 @@ static int cram_decode_seq(cram_fd *fd, cram_container *c, cram_slice *s,
 		}
 		ref_pos += cr->len - seq_pos + 1;
 	    }
+	} else if (cr->ref_id >= 0) {
+	    // So alignment end can be computed even when not decoding sequence
+	    ref_pos += cr->len - seq_pos + 1;
 	}
 
 	if (ncigar+1 >= cigar_alloc) {
@@ -2938,7 +2941,8 @@ static cram_slice *cram_next_slice(cram_fd *fd, cram_container **cp) {
 	if (fd->range.refid != -2) {
 	    while (c->ref_seq_id != -2 &&
 		   (c->ref_seq_id < fd->range.refid ||
-		    c->ref_seq_start + c->ref_seq_span-1 < fd->range.start)) {
+                    (fd->range.refid >= 0 && c->ref_seq_id == fd->range.refid
+                     && c->ref_seq_start + c->ref_seq_span-1 < fd->range.start))) {
 		if (0 != cram_seek(fd, c->length, SEEK_CUR))
 		    return NULL;
 		cram_free_container(fd->ctr);
@@ -2948,9 +2952,11 @@ static cram_slice *cram_next_slice(cram_fd *fd, cram_container **cp) {
 		} while (c->length == 0);
 	    }
 
-	    if (c->ref_seq_id != -2 && c->ref_seq_id != fd->range.refid)
+	    if (c->ref_seq_id != -2 && c->ref_seq_id != fd->range.refid) {
+                fd->eof = 1;
 		return NULL;
             }
+	}
 
 	if (!(c->comp_hdr_block = cram_read_block(fd)))
 	    return NULL;
@@ -3225,7 +3231,7 @@ cram_record *cram_get_seq(cram_fd *fd) {
 /*
  * Read the next cram record and convert it to a bam_seq_t struct.
  *
- * Returns 0 on success
+ * Returns >= 0 success (number of bytes written to *bam)
  *        -1 on EOF or failure (check fd->err)
  */
 int cram_get_bam_seq(cram_fd *fd, bam_seq_t **bam) {

cram/cram_index.c

@@ -345,7 +345,10 @@ void cram_index_free(cram_fd *fd) {
  * "from" as the last slice we checked to find the next one. Otherwise
  * set "from" to be NULL to find the first one.
  *
- * Returns the cram_index pointer on sucess
+ * Refid can also be any of the special HTS_IDX_ values.
+ * For backwards compatibility, refid -1 is equivalent to HTS_IDX_NOCOOR.
+ *
+ * Returns the cram_index pointer on success
  *         NULL on failure
  */
 cram_index *cram_index_query(cram_fd *fd, int refid, int pos, 
@@ -353,8 +356,33 @@ cram_index *cram_index_query(cram_fd *fd, int refid, int pos,
     int i, j, k;
     cram_index *e;
 
-    if (refid+1 < 0 || refid+1 >= fd->index_sz)
+    switch(refid) {
+    case HTS_IDX_NONE:
+    case HTS_IDX_REST:
+	// fail, or already there, dealt with elsewhere.
+	return NULL;
+
+    case HTS_IDX_NOCOOR:
+	refid = -1;
+	break;
+
+    case HTS_IDX_START: {
+	int64_t min_idx = INT64_MAX;
+	for (i = 0, j = -1; i < fd->index_sz; i++) {
+	    if (fd->index[i].e && fd->index[i].e[0].offset < min_idx) {
+		min_idx = fd->index[i].e[0].offset;
+		j = i;
+	    }
+	}
+	if (j < 0)
 	    return NULL;
+	return fd->index[j].e;
+    }
+
+    default:
+	if (refid < HTS_IDX_NONE || refid+1 >= fd->index_sz)
+	    return NULL;
+    }
 
     if (!from)
 	from = &fd->index[refid+1];
@@ -406,6 +434,24 @@ cram_index *cram_index_query(cram_fd *fd, int refid, int pos,
     return e;
 }
 
+// Return the index entry for last slice on a specific reference.
+cram_index *cram_index_last(cram_fd *fd, int refid, cram_index *from) {
+    int slice;
+
+    if (refid+1 < 0 || refid+1 >= fd->index_sz)
+        return NULL;
+
+    if (!from)
+        from = &fd->index[refid+1];
+
+    // Ref with nothing aligned against it.
+    if (!from->e)
+        return NULL;
+
+    slice = fd->index[refid+1].nslice - 1;
+
+    return &from->e[slice];
+}
 
 /*
  * Skips to a container overlapping the start coordinate listed in
@@ -424,6 +470,9 @@ cram_index *cram_index_query(cram_fd *fd, int refid, int pos,
 int cram_seek_to_refpos(cram_fd *fd, cram_range *r) {
     cram_index *e;
 
+    if (r->refid == HTS_IDX_NONE)
+	return -2;
+
     // Ideally use an index, so see if we have one.
     if ((e = cram_index_query(fd, r->refid, r->start, NULL))) {
 	if (0 != cram_seek(fd, e->offset, SEEK_SET))
@@ -434,6 +483,9 @@ int cram_seek_to_refpos(cram_fd *fd, cram_range *r) {
 	return -2;
     }
 
+    if (r->refid == HTS_IDX_START || r->refid == HTS_IDX_REST)
+	fd->range.refid = -2; // special case in cram_next_slice
+
     if (fd->ctr) {
 	cram_free_container(fd->ctr);
 	fd->ctr = NULL;

cram/cram_index.h

@@ -52,6 +52,7 @@ void cram_index_free(cram_fd *fd);
  *         NULL on failure
  */
 cram_index *cram_index_query(cram_fd *fd, int refid, int pos, cram_index *frm);
+cram_index *cram_index_last(cram_fd *fd, int refid, cram_index *from);
 
 /*
  * Skips to a container overlapping the start coordinate listed in

cram/cram_io.c

@@ -57,7 +57,11 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #include <bzlib.h>
 #endif
 #ifdef HAVE_LIBLZMA
+#ifdef HAVE_LZMA_H
 #include <lzma.h>
+#else
+#include "os/lzma_stub.h"
+#endif
 #endif
 #include <sys/types.h>
 #include <sys/stat.h>
@@ -4129,8 +4133,9 @@ int cram_seek(cram_fd *fd, off_t offset, int whence) {
 
     fd->ooc = 0;
 
-    if (hseek(fd->fp, offset, whence) >= 0)
+    if (hseek(fd->fp, offset, whence) >= 0) {
         return 0;
+    }
 
     if (!(whence == SEEK_CUR && offset >= 0))
 	return -1;

cram/cram_io.h

@@ -842,6 +842,8 @@ int cram_close(cram_fd *fd);
  */
 int cram_seek(cram_fd *fd, off_t offset, int whence);
 
+int64_t cram_tell(cram_fd *fd);
+
 /*
  * Flushes a CRAM file.
  * Useful for when writing to stdout without wishing to close the stream.

cram/cram_samtools.c

@@ -124,7 +124,7 @@ int bam_construct_seq(bam_seq_t **bp, size_t extra_len,
     else
 	memset(cp, '\xff', len);
 
-    return 0;
+    return bam_len;
 }
 
 bam_hdr_t *cram_header_to_bam(SAM_hdr *h) {

cram/mFILE.h

@@ -81,7 +81,6 @@ void mfrecreate(mFILE *mf, char *data, int size);
 void *mfsteal(mFILE *mf, size_t *size_out);
 char *mfgets(char *s, int size, mFILE *mf);
 int mfflush(mFILE *mf);
-int mfprintf(mFILE *mf, char *fmt, ...);
 mFILE *mstdin(void);
 mFILE *mstdout(void);
 mFILE *mstderr(void);

debian/changelog

+htslib (1.7-1) unstable; urgency=medium
+
+  * Team upload.
+
+  * New upstream version.
+  * Standards-Version: 4.1.3
+
+  * win/* files no longer in upstream source tree, affecting -test package
+
+ -- Steffen Moeller <moeller@debian.org>  Mon, 12 Feb 2018 18:59:19 +0100
+
 htslib (1.6-4) unstable; urgency=medium
 
   * Fix autopkgtest: Ship new files win/rand.c and win/rand.h

debian/control

@@ -11,7 +11,7 @@ Build-Depends: debhelper (>= 10),
                liblzma-dev,
                libssl-dev,
                zlib1g-dev
-Standards-Version: 4.1.2
+Standards-Version: 4.1.3
 Vcs-Browser: https://anonscm.debian.org/cgit/debian-med/htslib.git
 Vcs-Git: https://anonscm.debian.org/git/debian-med/htslib.git -b debian/unstable
 Homepage: https://github.com/samtools/htslib

debian/copyright

 Format: http://www.debian.org/doc/packaging-manuals/copyright-format/1.0/
-Source: https://github.com/samtools/htslib/archive/1.3.1.tar.gz
+Source: https://github.com/samtools/htslib/releases
 
 Files: *
 Copyright: © 2012-2016 Genome Research Ltd.

debian/htslib-test.install

@@ -5,5 +5,3 @@ Makefile /usr/share/htslib-test
 cram/*.c /usr/share/htslib-test/cram
 cram/*.h /usr/share/htslib-test/cram
 test/* /usr/share/htslib-test/test
-win/*.c usr/share/htslib-test/win
-win/*.h usr/share/htslib-test/win

debian/patches/series

 define_PATH_MAX.patch
 htslib-add-cram_to_bam.patch
 fPIC.patch
-877670.patch
+#877670.patch	# applied by upstream