Adding ability to provide bed file for sambamba slice

.test_suite.sh

@@ -26,6 +26,12 @@ testSliceMultipleRegions() {
     assertEquals 0 $?
 }
 
+testSliceMultipleRegionsBed() {
+    ./build/sambamba slice -L test/chr2_test_region.bed test/issue_204.bam |\
+      ./build/sambamba view -c /dev/stdin > /dev/null
+    assertEquals 0 $?
+}
+
 testView() {
     assertEquals "1806" `./build/sambamba view -c ex1_header.sorted.bam chr2`
     assertEquals "1464" `./build/sambamba view -c ex1_header.sorted.bam chr1`

sambamba/slice.d

@@ -35,6 +35,7 @@ import std.parallelism;
 import std.conv;
 import std.stdio;
 import std.exception;
+import sambamba.utils.common.bed;
 
 import sambamba.utils.common.overwrite;
 
@@ -167,6 +168,7 @@ void fetchRegions(BamReader bam, Region[] regions, ref Stream stream)
         // write R1
         foreach (read; r1.data)
             writer.writeRecord(read);
+        // Flush the BamWriter before our copyAsIs call.
         writer.flush();
 
         copyAsIs(bam, stream, s2_start_offset, s2_end_offset);
@@ -271,6 +273,8 @@ void printUsage()
     stderr.writeln();
     stderr.writeln("OPTIONS: -o, --output-filename=OUTPUT_FILENAME");
     stderr.writeln("            output BAM filename");
+    stderr.writeln("         -L, --regions=FILENAME");
+    stderr.writeln("            output only reads overlapping one of regions from the BED file");
 }
 
 version(standalone) {
@@ -283,17 +287,24 @@ int slice_main(string[] args) {
 
     // at least two arguments must be presented
     string output_filename = null;
+    string bed_filename = null;
 
     try {
         getopt(args,
                std.getopt.config.caseSensitive,
-               "output-filename|o", &output_filename);
+               "output-filename|o", &output_filename,
+               "regions|L",         &bed_filename);
 
-        if (args.length < 3) {
+        if ( (bed_filename.length == 0 && args.length < 3) ||
+             (args.length < 2) )  {
             printUsage();
             return 0;
         }
 
+        if (bed_filename.length > 0 && args.length > 2) {
+            throw new Exception("specifying both region and BED filename is disallowed");
+        }
+
         protectFromOverwrite(args[1], output_filename);
 
         import std.parallelism;
@@ -321,7 +332,13 @@ int slice_main(string[] args) {
         if (args[2] == "*") {
             fetchUnmapped(bam, stream);
         } else {
-            auto regions = map!parseRegion(args[2 .. $]).array;
+            Region[] regions;
+            if (bed_filename !is null) {
+                auto bam_regions = parseBed(bed_filename, bam);
+                regions = cast(Region[])bam_regions.map!(r => Region(bam.reference(r.ref_id).name, r.start, r.end)).array;
+            } else {
+                regions = cast(Region[])map!parseRegion(args[2 .. $]).array;
+            }
             fetchRegions(bam, regions, stream);
         }
 

test/chr2_test_region.bed

@@ -0,0 +1,3 @@
+2	166860000	166870000
+2	166875000	166890000
+3	100000000	200000000