Release 1.4.1: summary

Makefile

@@ -24,6 +24,7 @@
 CC       = gcc
 AR       = ar
 CPPFLAGS =
+#CFLAGS   = -g -Wall -O2 -pedantic -std=c99 -D_XOPEN_SOURCE=600
 CFLAGS   = -g -Wall -O2
 LDFLAGS  =
 LIBS     =
@@ -102,7 +103,7 @@ config.h:
 include config.mk
 
 
-PACKAGE_VERSION = 1.4
+PACKAGE_VERSION = 1.4.1
 
 # If building from a Git repository, replace $(PACKAGE_VERSION) with the Git
 # description of the working tree: either a release tag with the same value
@@ -142,7 +143,7 @@ libbam.a:$(LOBJS)
 	$(AR) -csru $@ $(LOBJS)
 
 samtools: $(AOBJS) libbam.a libst.a $(HTSLIB)
-	$(CC) -pthread $(ALL_LDFLAGS) -o $@ $(AOBJS) libbam.a libst.a $(HTSLIB_LIB) $(CURSES_LIB) -lm $(ALL_LIBS)
+	$(CC) $(ALL_LDFLAGS) -o $@ $(AOBJS) libbam.a libst.a $(HTSLIB_LIB) $(CURSES_LIB) -lm $(ALL_LIBS) -lpthread
 
 # For building samtools and its test suite only: NOT to be installed.
 libst.a: $(LIBST_OBJS)
@@ -223,28 +224,28 @@ check test: samtools $(BGZIP) $(TEST_PROGRAMS)
 
 
 test/merge/test_bam_translate: test/merge/test_bam_translate.o test/test.o libst.a $(HTSLIB)
-	$(CC) -pthread $(ALL_LDFLAGS) -o $@ test/merge/test_bam_translate.o test/test.o libst.a $(HTSLIB_LIB) $(ALL_LIBS)
+	$(CC) $(ALL_LDFLAGS) -o $@ test/merge/test_bam_translate.o test/test.o libst.a $(HTSLIB_LIB) $(ALL_LIBS) -lpthread
 
 test/merge/test_rtrans_build: test/merge/test_rtrans_build.o libst.a $(HTSLIB)
-	$(CC) -pthread $(ALL_LDFLAGS) -o $@ test/merge/test_rtrans_build.o libst.a $(HTSLIB_LIB) $(ALL_LIBS)
+	$(CC) $(ALL_LDFLAGS) -o $@ test/merge/test_rtrans_build.o libst.a $(HTSLIB_LIB) $(ALL_LIBS) -lpthread
 
 test/merge/test_trans_tbl_init: test/merge/test_trans_tbl_init.o libst.a $(HTSLIB)
-	$(CC) -pthread $(ALL_LDFLAGS) -o $@ test/merge/test_trans_tbl_init.o libst.a $(HTSLIB_LIB) $(ALL_LIBS)
+	$(CC) $(ALL_LDFLAGS) -o $@ test/merge/test_trans_tbl_init.o libst.a $(HTSLIB_LIB) $(ALL_LIBS) -lpthread
 
 test/split/test_count_rg: test/split/test_count_rg.o test/test.o libst.a $(HTSLIB)
-	$(CC) -pthread $(ALL_LDFLAGS) -o $@ test/split/test_count_rg.o test/test.o libst.a $(HTSLIB_LIB) $(ALL_LIBS)
+	$(CC) $(ALL_LDFLAGS) -o $@ test/split/test_count_rg.o test/test.o libst.a $(HTSLIB_LIB) $(ALL_LIBS) -lpthread
 
 test/split/test_expand_format_string: test/split/test_expand_format_string.o test/test.o libst.a $(HTSLIB)
-	$(CC) -pthread $(ALL_LDFLAGS) -o $@ test/split/test_expand_format_string.o test/test.o libst.a $(HTSLIB_LIB) $(ALL_LIBS)
+	$(CC) $(ALL_LDFLAGS) -o $@ test/split/test_expand_format_string.o test/test.o libst.a $(HTSLIB_LIB) $(ALL_LIBS) -lpthread
 
 test/split/test_filter_header_rg: test/split/test_filter_header_rg.o test/test.o libst.a $(HTSLIB)
-	$(CC) -pthread $(ALL_LDFLAGS) -o $@ test/split/test_filter_header_rg.o test/test.o libst.a $(HTSLIB_LIB) $(ALL_LIBS)
+	$(CC) $(ALL_LDFLAGS) -o $@ test/split/test_filter_header_rg.o test/test.o libst.a $(HTSLIB_LIB) $(ALL_LIBS) -lpthread
 
 test/split/test_parse_args: test/split/test_parse_args.o test/test.o libst.a $(HTSLIB)
-	$(CC) -pthread $(ALL_LDFLAGS) -o $@ test/split/test_parse_args.o test/test.o libst.a $(HTSLIB_LIB) $(ALL_LIBS)
+	$(CC) $(ALL_LDFLAGS) -o $@ test/split/test_parse_args.o test/test.o libst.a $(HTSLIB_LIB) $(ALL_LIBS) -lpthread
 
 test/vcf-miniview: test/vcf-miniview.o $(HTSLIB)
-	$(CC) -pthread $(ALL_LDFLAGS) -o $@ test/vcf-miniview.o $(HTSLIB_LIB) $(ALL_LIBS)
+	$(CC) $(ALL_LDFLAGS) -o $@ test/vcf-miniview.o $(HTSLIB_LIB) $(ALL_LIBS) -lpthread
 
 test_test_h = test/test.h $(htslib_sam_h)
 

NEWS

@@ -1,4 +1,29 @@
-Beta Release 1.4  (DD Mmmmm 2016)
+Release 1.4.1  (8th May 2017)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+* Added options to fastq to create fastq files from BC (or other)
+  tags.
+
+* Samtools view has gained a -G <flags> option to exclude on all bits
+  set.  For example to discard reads where neither end has been
+  mapped use "-G 12".
+
+* Samtools cat has a -b <fofn> option to ease concatenation of many
+  files.
+
+* Added misc/samtools_tab_completion for bash auto-completion of
+  samtools sub-commands. (#560)
+
+* Samtools tview now has J and K keys for verticale movement by 20
+  lines. (#257)
+
+* Various compilation / portability improvements.
+
+* Fixed issue with more than 65536 CIGAR operations and SAM/CRAM files.
+  (#667)
+
+
+Release 1.4  (13 March 2017)
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 Noteworthy changes in samtools:

README

@@ -9,8 +9,8 @@ Building samtools
 The typical simple case of building Samtools using the HTSlib bundled within
 this Samtools release tarball and enabling useful plugins, is done as follows:
 
-    cd .../samtools-1.4 # Within the unpacked release directory
-    ./configure --enable-plugins --enable-libcurl --with-plugin-path=$PWD/htslib-1.4
+    cd .../samtools-1.4.1 # Within the unpacked release directory
+    ./configure --enable-plugins --enable-libcurl --with-plugin-path=$PWD/htslib-1.4.1
     make all plugins-htslib
 
 You may wish to copy the resulting samtools executable into somewhere on your
@@ -22,7 +22,7 @@ installation using the HTSlib bundled within this Samtools release tarball,
 also building the various HTSlib utilities such as bgzip and enabling useful
 plugins, is done as follows:
 
-    cd .../samtools-1.4 # Within the unpacked release directory
+    cd .../samtools-1.4.1 # Within the unpacked release directory
     ./configure --enable-plugins --enable-libcurl --prefix=/path/to/location
     make all all-htslib
     make install install-htslib

bam.h

@@ -38,7 +38,7 @@ DEALINGS IN THE SOFTWARE.  */
   @copyright Genome Research Ltd.
  */
 
-#define BAM_VERSION "1.4"
+#define BAM_VERSION "1.4.1"
 
 #include <stdint.h>
 #include <stdlib.h>

bam2bcf_indel.c

@@ -439,10 +439,13 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla
     }
     free(ref2); free(query);
     { // compute indelQ
-        int *sc, tmp, *sumq;
-        sc   = alloca(n_types * sizeof(int));
-        sumq = alloca(n_types * sizeof(int));
-        memset(sumq, 0, sizeof(int) * n_types);
+        int sc_a[16], sumq_a[16];
+        int tmp, *sc = sc_a, *sumq = sumq_a;
+        if (n_types > 16) {
+            sc   = (int *)malloc(n_types * sizeof(int));
+            sumq = (int *)malloc(n_types * sizeof(int));
+        }
+        memset(sumq, 0, n_types * sizeof(int));
         for (s = K = 0; s < n; ++s) {
             for (i = 0; i < n_plp[s]; ++i, ++K) {
                 bam_pileup1_t *p = plp[s] + i;
@@ -523,6 +526,9 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla
                 //fprintf(stderr, "X pos=%d read=%d:%d name=%s call=%d type=%d seqQ=%d indelQ=%d\n", pos, s, i, bam1_qname(p->b), (p->aux>>16)&0x3f, bca->indel_types[(p->aux>>16)&0x3f], (p->aux>>8)&0xff, p->aux&0xff);
             }
         }
+
+        if (sc   != sc_a)   free(sc);
+        if (sumq != sumq_a) free(sumq);
     }
     free(score1); free(score2);
     // free

bam_cat.c

@@ -40,6 +40,7 @@ Illumina.
 #include <stdlib.h>
 #include <unistd.h>
 #include <string.h>
+#include <strings.h>
 
 #include "htslib/bgzf.h"
 #include "htslib/sam.h"
@@ -468,7 +469,7 @@ int bam_cat(int nfn, char * const *fn, const bam_hdr_t *h, const char* outbam)
         }
 
         if (in->block_offset < in->block_length) {
-            if (bgzf_write(fp, in->uncompressed_block + in->block_offset, in->block_length - in->block_offset) < 0) goto write_fail;
+            if (bgzf_write(fp, (char *)in->uncompressed_block + in->block_offset, in->block_length - in->block_offset) < 0) goto write_fail;
             if (bgzf_flush(fp) != 0) goto write_fail;
         }
 
@@ -531,10 +532,12 @@ int main_cat(int argc, char *argv[])
 {
     bam_hdr_t *h = 0;
     char *outfn = 0;
+    char **infns = NULL; // files to concatenate
+    int infns_size = 0;
     int c, ret = 0;
     samFile *in;
 
-    while ((c = getopt(argc, argv, "h:o:")) >= 0) {
+    while ((c = getopt(argc, argv, "h:o:b:")) >= 0) {
         switch (c) {
             case 'h': {
                 samFile *fph = sam_open(optarg, "r");
@@ -553,29 +556,61 @@ int main_cat(int argc, char *argv[])
                 break;
             }
             case 'o': outfn = strdup(optarg); break;
+            case 'b': {
+                // add file names in "optarg" to the list
+                // of files to concatenate
+                int nfns;
+                char **fns_read = hts_readlines(optarg, &nfns);
+                if (fns_read) {
+                    infns = realloc(infns, (infns_size + nfns) * sizeof(char*));
+                    if (infns == NULL) { ret = 1; goto end; }
+                    memcpy(infns+infns_size, fns_read, nfns * sizeof(char*));
+                    infns_size += nfns;
+                    free(fns_read);
+                } else {
+                    print_error("cat", "Invalid file list \"%s\"", optarg);
+                    ret = 1;
+                }
+                break;
+            }
         }
     }
-    if (argc - optind < 1) {
-        fprintf(stderr, "Usage: samtools cat [-h header.sam] [-o out.bam] <in1.bam> [...]\n");
+
+    // Append files specified in argv to the list.
+    int nargv_fns = argc - optind;
+    if (nargv_fns > 0) {
+        infns = realloc(infns, (infns_size + nargv_fns) * sizeof(char*));
+        if (infns == NULL) { ret = 1; goto end; }
+        memcpy(infns + infns_size, argv + optind, nargv_fns * sizeof(char*));
+    }
+
+    // Require at least one input file
+    if (infns_size + nargv_fns == 0) {
+        fprintf(stderr, "Usage: samtools cat [options] <in1.bam>  [... <inN.bam>]\n");
+        fprintf(stderr, "       samtools cat [options] <in1.cram> [... <inN.cram>]\n\n");
+        fprintf(stderr, "Concatenate BAM or CRAM files, first those in <bamlist.fofn>, then those\non the command line.\n\n");
+        fprintf(stderr, "Options: -b FILE  list of input BAM/CRAM file names, one per line\n");
+        fprintf(stderr, "         -h FILE  copy the header from FILE [default is 1st input file]\n");
+        fprintf(stderr, "         -o FILE  output BAM/CRAM\n");
         return 1;
     }
 
-    in = sam_open(argv[optind], "r");
+    in = sam_open(infns[0], "r");
     if (!in) {
-        print_error_errno("cat", "failed to open file '%s'", argv[optind]);
+        print_error_errno("cat", "failed to open file '%s'", infns[0]);
         return 1;
     }
 
     switch (hts_get_format(in)->format) {
     case bam:
         sam_close(in);
-        if (bam_cat(argc - optind, argv + optind, h, outfn? outfn : "-") < 0)
+        if (bam_cat(infns_size+nargv_fns, infns, h, outfn? outfn : "-") < 0)
             ret = 1;
         break;
 
     case cram:
         sam_close(in);
-        if (cram_cat(argc - optind, argv + optind, h, outfn? outfn : "-") < 0)
+        if (cram_cat(infns_size+nargv_fns, infns, h, outfn? outfn : "-") < 0)
             ret = 1;
         break;
 
@@ -584,7 +619,16 @@ int main_cat(int argc, char *argv[])
         fprintf(stderr, "[%s] ERROR: input is not BAM or CRAM\n", __func__);
         return 1;
     }
+
+ end:
+    if (infns_size > 0) {
+        int i;
+        for (i=0; i<infns_size; i++)
+            free(infns[i]);
+    }
+
     free(outfn);
+    free(infns);
 
     if (h)
         bam_hdr_destroy(h);

bam_plcmd.c

@@ -31,6 +31,7 @@ DEALINGS IN THE SOFTWARE.  */
 #include <unistd.h>
 #include <ctype.h>
 #include <string.h>
+#include <strings.h>
 #include <limits.h>
 #include <errno.h>
 #include <sys/stat.h>

bam_reheader.c

@@ -91,7 +91,7 @@ int bam_reheader(BGZF *in, bam_hdr_t *h, int fd,
         goto fail;
     }
     if (in->block_offset < in->block_length) {
-        if (bgzf_write(fp, in->uncompressed_block + in->block_offset, in->block_length - in->block_offset) < 0) goto write_fail;
+        if (bgzf_write(fp, (char *)in->uncompressed_block + in->block_offset, in->block_length - in->block_offset) < 0) goto write_fail;
         if (bgzf_flush(fp) < 0) goto write_fail;
     }
     while ((len = bgzf_raw_read(in, buf, BUF_SIZE)) > 0) {
@@ -246,7 +246,7 @@ int cram_reheader_inplace2(cram_fd *fd, const bam_hdr_t *h, const char *arg_list
     int32_put_blk(b, header_len);
     cram_block_append(b, sam_hdr_str(hdr), header_len);
     // Zero the remaining block
-    memset(cram_block_get_data(b)+cram_block_get_offset(b), 0,
+    memset((char *)cram_block_get_data(b)+cram_block_get_offset(b), 0,
            cram_block_get_uncomp_size(b) - cram_block_get_offset(b));
     // Make sure all sizes and byte-offsets are consistent after memset
     cram_block_set_offset(b, cram_block_get_uncomp_size(b));

bam_sort.c

@@ -1147,14 +1147,6 @@ int bam_merge_core2(int by_qname, const char *out, const char *mode,
             print_error(cmd, "couldn't read headers from \"%s\"", headers);
             goto mem_fail;
         }
-    } else  {
-        hout = bam_hdr_init();
-        if (!hout) {
-            print_error(cmd, "couldn't allocate bam header");
-            goto mem_fail;
-        }
-        hout->text = strdup("");
-        if (!hout->text) goto mem_fail;
     }
 
     g_is_by_qname = by_qname;
@@ -1227,6 +1219,16 @@ int bam_merge_core2(int by_qname, const char *out, const char *mode,
         if ((translation_tbl+i)->lost_coord_sort && !by_qname) {
             fprintf(stderr, "[bam_merge_core] Order of targets in file %s caused coordinate sort to be lost\n", fn[i]);
         }
+
+        // Potential future improvement is to share headers between CRAM files for
+        // samtools sort (where all headers are identical.
+        // Eg:
+        //
+        // if (i > 1) {
+        //     sam_hdr_free(cram_fd_get_header(fp[i]->fp.cram));
+        //     cram_fd_set_header(fp[i]->fp.cram, cram_fd_get_header(fp[0]->fp.cram));
+        //     sam_hdr_incr_ref(cram_fd_get_header(fp[0]->fp.cram));
+        // }
     }
 
     // Did we get an @HD line?
@@ -1504,6 +1506,7 @@ int bam_merge(int argc, char *argv[])
                 if (fn == NULL) { ret = 1; goto end; }
                 memcpy(fn+fn_size, fn_read, nfiles * sizeof(char*));
                 fn_size += nfiles;
+                free(fn_read);
             }
             else {
                 print_error("merge", "Invalid file list \"%s\"", optarg);
@@ -1655,18 +1658,30 @@ static void *worker(void *data)
     name = (char*)calloc(strlen(w->prefix) + 20, 1);
     if (!name) { w->error = errno; return 0; }
     sprintf(name, "%s.%.4d.bam", w->prefix, w->index);
-    if (write_buffer(name, "wbx1", w->buf_len, w->buf, w->h, 0, NULL) < 0)
-        w->error = errno;
-
-// Consider using CRAM temporary files if the final output is CRAM.
-// Typically it is comparable speed while being smaller.
-//    hts_opt opt[2] = {
-//        {"version=3.0", CRAM_OPT_VERSION, {"3.0"}, NULL},
-//        {"no_ref",      CRAM_OPT_NO_REF,  {1},     NULL}
-//    };
-//    opt[0].next = &opt[1];
-//    if (write_buffer(name, "wc1", w->buf_len, w->buf, w->h, 0, opt) < 0)
-//        w->error = errno;
+
+    uint32_t max_ncigar = 0;
+    int i;
+    for (i = 0; i < w->buf_len; i++) {
+        uint32_t nc = w->buf[i]->core.n_cigar;
+        if (max_ncigar < nc)
+            max_ncigar = nc;
+    }
+
+    if (max_ncigar > 65535) {
+        htsFormat fmt;
+        memset(&fmt, 0, sizeof(fmt));
+        if (hts_parse_format(&fmt, "cram,version=3.0,no_ref,seqs_per_slice=1000") < 0) {
+            w->error = errno;
+            free(name);
+            return 0;
+        }
+
+        if (write_buffer(name, "wcx1", w->buf_len, w->buf, w->h, 0, &fmt) < 0)
+            w->error = errno;
+    } else {
+        if (write_buffer(name, "wbx1", w->buf_len, w->buf, w->h, 0, NULL) < 0)
+            w->error = errno;
+    }
 
     free(name);
     return 0;

bam_tview_curses.c

@@ -277,8 +277,10 @@ static int curses_loop(tview_t* tv)
             case ' ': pos += tv->mcol; break;
             case KEY_UP:
             case 'j': --tv->row_shift; break;
+            case 'J': tv->row_shift -= 20; break;
             case KEY_DOWN:
             case 'k': ++tv->row_shift; break;
+            case 'K': tv->row_shift += 20; break;
             case KEY_BACKSPACE:
             case '\177': pos -= tv->mcol; break;
             case KEY_RESIZE: getmaxyx(stdscr, tv->mrow, tv->mcol); break;