Merge Enable samtools index on multiple alignment files (PR samtools#1674)

Author: daviesrob

Makefile

@@ -173,7 +173,7 @@ bam_cat.o: bam_cat.c config.h $(htslib_bgzf_h) $(htslib_sam_h) $(htslib_cram_h)
 bam_color.o: bam_color.c config.h $(htslib_sam_h)
 bam_fastq.o: bam_fastq.c config.h $(htslib_sam_h) $(htslib_klist_h) $(htslib_kstring_h) $(htslib_bgzf_h) $(htslib_thread_pool_h) $(samtools_h) $(sam_opts_h)
 bam_import.o: bam_import.c config.h $(htslib_sam_h) $(htslib_thread_pool_h) $(samtools_h) $(sam_opts_h)
-bam_index.o: bam_index.c config.h $(htslib_hts_h) $(htslib_sam_h) $(htslib_khash_h) $(samtools_h) $(sam_opts_h)
+bam_index.o: bam_index.c config.h $(htslib_hts_h) $(htslib_sam_h) $(htslib_hfile_h) $(htslib_khash_h) $(samtools_h) $(sam_opts_h)
 bam_lpileup.o: bam_lpileup.c config.h $(bam_plbuf_h) $(bam_lpileup_h) splaysort.h
 bam_mate.o: bam_mate.c config.h $(htslib_thread_pool_h) $(sam_opts_h) $(htslib_kstring_h) $(htslib_sam_h) $(samtools_h)
 bam_md.o: bam_md.c config.h $(htslib_faidx_h) $(htslib_sam_h) $(htslib_kstring_h) $(htslib_thread_pool_h) $(sam_opts_h) $(samtools_h)

NEWS

@@ -4,6 +4,11 @@ Release a.b
  * The dict command can now read BWA's .alt file and add AH:* tags
    indicating reference sequences that represent alternate loci.
 
+ * The "samtools index" command can now accept multiple alignment filenames
+   with the new -M option, and will index each of them separately. (Specifying
+   the output index filename via out.index or the new -o option is currently
+   only applicable when there is only one alignment file to be indexed.)
+
 
 Release 1.15.1 (7th April 2022)
 -------------------------------

bam_index.c

@@ -28,10 +28,10 @@ DEALINGS IN THE SOFTWARE.  */
 
 #include <htslib/hts.h>
 #include <htslib/sam.h>
+#include <htslib/hfile.h>
 #include <htslib/khash.h>
 #include <stdlib.h>
 #include <stdio.h>
-#define __STDC_FORMAT_MACROS
 #include <inttypes.h>
 #include <unistd.h>
 #include <getopt.h>
@@ -44,63 +44,97 @@ DEALINGS IN THE SOFTWARE.  */
 static void index_usage(FILE *fp)
 {
     fprintf(fp,
-"Usage: samtools index [-bc] [-m INT] <in.bam> [out.index]\n"
+"Usage: samtools index -M [-bc] [-m INT] <in1.bam> <in2.bam>...\n"
+"   or: samtools index [-bc] [-m INT] <in.bam> [out.index]\n"
 "Options:\n"
 "  -b       Generate BAI-format index for BAM files [default]\n"
 "  -c       Generate CSI-format index for BAM files\n"
 "  -m INT   Set minimum interval size for CSI indices to 2^INT [%d]\n"
+"  -M       Interpret all filename arguments as files to be indexed\n"
+"  -o FILE  Write index to FILE [alternative to <out.index> as an argument]\n"
 "  -@ INT   Sets the number of threads [none]\n", BAM_LIDX_SHIFT);
 }
 
+// Returns 1 if the file does not exist or can be positively
+// identified as an index file.
+static int nonexistent_or_index(const char *fn)
+{
+    int ret1, ret2;
+    htsFormat fmt;
+    hFILE *fp = hopen(fn, "r");
+    if (fp == NULL) return 1;
+
+    ret1 = hts_detect_format2(fp, fn, &fmt);
+    ret2 = hclose(fp);
+    if (ret1 < 0 || ret2 < 0) return 0;
+
+    return fmt.category == index_file;
+}
+
 int bam_index(int argc, char *argv[])
 {
     int csi = 0;
     int min_shift = BAM_LIDX_SHIFT;
+    int multiple = 0;
     int n_threads = 0;
-    int c, ret;
+    int n_files, c, i, ret;
+    const char *fn_idx = NULL;
 
-    while ((c = getopt(argc, argv, "bcm:@:")) >= 0)
+    while ((c = getopt(argc, argv, "bcm:Mo:@:")) >= 0)
         switch (c) {
         case 'b': csi = 0; break;
         case 'c': csi = 1; break;
         case 'm': csi = 1; min_shift = atoi(optarg); break;
+        case 'M': multiple = 1; break;
+        case 'o': fn_idx = optarg; break;
         case '@': n_threads = atoi(optarg); break;
         default:
             index_usage(stderr);
             return 1;
         }
 
-    if (optind == argc) {
-        index_usage(stdout);
-        return 1;
-    }
+    n_files = argc - optind;
 
-    ret = sam_index_build3(argv[optind], argv[optind+1], csi? min_shift : 0, n_threads);
-    switch (ret) {
-    case 0:
+    if (n_files == 0) {
+        index_usage(stdout);
         return 0;
+    }
 
-    case -2:
-        print_error_errno("index", "failed to open \"%s\"", argv[optind]);
-        break;
+    // Handle legacy synopsis
+    if (n_files == 2 && !fn_idx && nonexistent_or_index(argv[optind+1])) {
+        n_files = 1;
+        fn_idx = argv[optind+1];
+    }
 
-    case -3:
-        print_error("index", "\"%s\" is in a format that cannot be usefully indexed", argv[optind]);
-        break;
+    if (n_files > 1 && !multiple) {
+        print_error("index", "use -M to enable indexing more than one alignment file");
+        return EXIT_FAILURE;
+    }
 
-    case -4:
-        if (argv[optind+1])
-            print_error("index", "failed to create or write index \"%s\"", argv[optind+1]);
-        else
-            print_error("index", "failed to create or write index");
-        break;
+    if (fn_idx && n_files > 1) {
+        // TODO In future we may allow %* placeholders or similar
+        print_error("index", "can't use -o with multiple input alignment files");
+        return EXIT_FAILURE;
+    }
 
-    default:
-        print_error_errno("index", "failed to create index for \"%s\"", argv[optind]);
-        break;
+    for (i = optind; i < optind + n_files; i++) {
+        ret = sam_index_build3(argv[i], fn_idx, csi? min_shift : 0, n_threads);
+        if (ret < 0) {
+            if (ret == -2)
+                print_error_errno("index", "failed to open \"%s\"", argv[i]);
+            else if (ret == -3)
+                print_error("index", "\"%s\" is in a format that cannot be usefully indexed", argv[i]);
+            else if (ret == -4 && fn_idx)
+                print_error("index", "failed to create or write index \"%s\"", fn_idx);
+            else if (ret == -4)
+                print_error("index", "failed to create or write index");
+            else
+                print_error_errno("index", "failed to create index for \"%s\"", argv[i]);
+            return EXIT_FAILURE;
+        }
     }
 
-    return EXIT_FAILURE;
+    return EXIT_SUCCESS;
 }
 
 /*

doc/samtools-index.1

@@ -43,7 +43,14 @@ samtools index \- indexes SAM/BAM/CRAM files
 .
 .SH SYNOPSIS
 .PP
-samtools index
+.B samtools index -M
+.RB [ -bc ]
+.RB [ -m
+.IR INT ]
+.I FILE FILE
+.RI [ FILE ...]
+.PP
+.B samtools index
 .RB [ -bc ]
 .RB [ -m
 .IR INT ]
@@ -52,19 +59,25 @@ samtools index
 
 .SH DESCRIPTION
 .PP
-Index a coordinate-sorted BGZIP-compressed SAM, BAM or CRAM file for fast
+Index coordinate-sorted BGZIP-compressed SAM, BAM or CRAM files for fast
 random access.
 Note for SAM this only works if the file has been BGZF compressed first.
+(The first synopsis with multiple input
+.IR FILE s
+is only available with Samtools 1.16 or later.)
 
 This index is needed when
 .I region
 arguments are used to limit
 .B samtools view
 and similar commands to particular regions of interest.
 
-If an output filename is given, the index file will be written to
-.IR out.index .
-Otherwise, for a CRAM file
+When only one alignment file is being indexed, the output index filename
+can be specified via
+.B -o
+or as shown in the second synopsis.
+
+When no output filename is specified, for a CRAM file
 .IR aln.cram ,
 index file
 .IB aln.cram .crai
@@ -101,6 +114,17 @@ as the fixed value used by the BAI format.
 .BI "-m " INT
 Create a CSI index, with a minimum interval size of 2^INT.
 .TP
+.B -M
+Interpret all filename arguments as alignment files to be indexed individually.
+(Without
+.BR -M ,
+filename arguments are interpreted solely as per the second synopsis.)
+.TP
+.BI "-o " FILE
+Write the output index to
+.IR FILE .
+(Currently may only be used when exactly one alignment file is being indexed.)
+.TP
 .BI "-@, --threads " INT
 Number of input/output compression threads to use in addition to main thread [0].
 

doc/samtools.1

@@ -233,6 +233,8 @@ samtools index
 
 Index a coordinate-sorted SAM, BAM or CRAM file for fast random access.
 Note for SAM this only works if the file has been BGZF compressed first.
+(Starting from Samtools 1.16, this command can also be given several
+alignment filenames, which are indexed individually.)
 
 This index is needed when
 .I region

test/test.pl

@@ -843,6 +843,15 @@ sub test_index
     test_cmd($opts,out=>'dat/test_input_1_b.X.expected',cmd=>"$$opts{bin}/samtools view${threads} -X $$opts{path}/dat/test_input_1_b.bam $$opts{tmp}/test_input_1_b.bam.bai ref2");
     test_cmd($opts,out=>'dat/test_input_1_ab.X.expected', ignore_pg_header => 1, cmd=>"$$opts{bin}/samtools merge${threads} -O sam - -X -cp -R ref2 $$opts{path}/dat/test_input_1_a.bam $$opts{path}/dat/test_input_1_b.bam $$opts{path}/dat/test_input_1_a.bam.bai $$opts{tmp}/test_input_1_b.bam.bai");
 
+    # Check -o option
+    cmd("$$opts{bin}/samtools index${threads} -o $$opts{tmp}/test_input_1_b_opt.bam.bai $$opts{path}/dat/test_input_1_b.bam");
+    test_cmd($opts,out=>'dat/test_input_1_b.X.expected',cmd=>"$$opts{bin}/samtools view${threads} -X $$opts{path}/dat/test_input_1_b.bam $$opts{tmp}/test_input_1_b_opt.bam.bai ref2");
+
+    # Check indexing multiple alignment files
+    test_cmd($opts,out=>'dat/empty.expected',cmd=>"$$opts{bin}/samtools index${threads} $$opts{path}/dat/test_input_1_a.bam $$opts{path}/dat/test_input_1_b.bam",want_fail=>1);
+    test_cmd($opts,out=>'dat/test_input_1_a.bam.bai.expected',cmd=>"$$opts{bin}/samtools index${threads} -M $$opts{path}/dat/test_input_1_a.bam $$opts{path}/dat/test_input_1_b.bam && cat $$opts{path}/dat/test_input_1_a.bam.bai",binary=>1);
+    test_cmd($opts,out=>'dat/test_input_1_b.X.expected',cmd=>"$$opts{bin}/samtools view${threads} -X $$opts{path}/dat/test_input_1_b.bam $$opts{path}/dat/test_input_1_b.bam.bai ref2");
+
     # Check auto-indexing
     cmd("$$opts{bin}/samtools view${threads} --write-index -o $$opts{path}/dat/auto_indexed.tmp.bam $$opts{path}/dat/mpileup.1.sam");
     test_cmd($opts,out=>"dat/auto_indexed.tmp.bam.csi", cmd=>"$$opts{bin}/samtools index${threads} -c $$opts{path}/dat/auto_indexed.tmp.bam $$opts{tmp}/auto_indexed.csi && cat $$opts{tmp}/auto_indexed.csi", binary=>1);