granges filter is working, match bedtools intersect output, and is 40…

…% faster.

 - Working version of granges filter subcommand, for filtering based
   on overlaps.
 - New integration test of granges filter vs bedtools intersect.
 - New benchmarks: granges filter vs bedtools intersect. We're fast!
   e.g. random one: granges: 6.28ms, bedtools: 8.80ms - about 40% faster.
 - Renamed `to_coitrees()` to `into_coitrees()`.
 - Fixed issue with the `AsGRangesRef` not being used when it should have
   been.
 - Bunch of new `filter_overlaps()` tests.
 - Fixed off-by-one bug with `From<coitrees::Interval>` to our range types.
 - Bunch of new tests for converting between range types to check by
   off-by-one errors.
 - Reshuffled `temp_bedfile()` and `random_bedfile()` into `test_utilities.rs`.

Cargo.toml

@@ -20,6 +20,7 @@ indexmap = "2.2.3"
 ndarray = "0.15.6"
 noodles = { version = "0.63.0", features = ["core", "bed"] }
 rand = "0.8.5"
+tempfile = "3.10.0"
 thiserror = "1.0.57"
 
 [features]
@@ -38,7 +39,6 @@ path = "src/main/mod.rs"
 
 [dev-dependencies]
 criterion = { version = "0.5.1", features = ["html_reports"] }
-tempfile = "3.10.0"
 
 [[bench]]
 name = "bedtools_comparison"

benches/bedtools_comparison.rs

@@ -5,32 +5,17 @@
 //! ensure the output is the *exact* same.
 
 use criterion::{criterion_group, criterion_main, Criterion};
-use granges::{commands::granges_random_bed, test_utilities::granges_binary_path};
+use granges::test_utilities::{granges_binary_path, random_bedfile};
 use std::process::Command;
-use tempfile::NamedTempFile;
 
-const BED_LENGTH: u32 = 10_000;
-
-/// Create a random BED3 file based on the hg38 sequence lengths, and write to disk.
-pub fn random_bedfile() -> NamedTempFile {
-    let temp_file = NamedTempFile::new().expect("Failed to create temp file");
-    let random_bedfile_path = temp_file.path().to_path_buf();
-    granges_random_bed(
-        "tests_data/hg38_seqlens.tsv",
-        BED_LENGTH,
-        Some(&random_bedfile_path),
-        true,
-    )
-    .expect("could not generate random BED file");
-    temp_file
-}
+const BED_LENGTH: usize = 10_000;
 
 fn bench_range_adjustment(c: &mut Criterion) {
     // create the benchmark group
     let mut group = c.benchmark_group("slop vs adjust");
 
     // create the test data
-    let input_bedfile = random_bedfile();
+    let input_bedfile = random_bedfile(BED_LENGTH);
 
     // configure the sample size for the group
     // group.sample_size(10);
@@ -68,5 +53,51 @@ fn bench_range_adjustment(c: &mut Criterion) {
         });
     });
 }
-criterion_group!(benches, bench_range_adjustment);
+
+fn bench_filter_adjustment(c: &mut Criterion) {
+    // create the benchmark group
+    let mut group = c.benchmark_group("intersect vs filter");
+
+    // create the test data
+    let random_bedfile_left_tempfile = random_bedfile(BED_LENGTH);
+    let random_bedfile_right_tempfile = random_bedfile(BED_LENGTH);
+    let random_bedfile_left = random_bedfile_left_tempfile.path();
+    let random_bedfile_right = random_bedfile_right_tempfile.path();
+
+    // configure the sample size for the group
+    // group.sample_size(10);
+    group.bench_function("bedtools_intersect", |b| {
+        b.iter(|| {
+            let bedtools_output = Command::new("bedtools")
+                .arg("intersect")
+                .arg("-a")
+                .arg(&random_bedfile_left)
+                .arg("-b")
+                .arg(&random_bedfile_right)
+                .arg("-wa")
+                .arg("-u")
+                .output()
+                .expect("bedtools intersect failed");
+            assert!(bedtools_output.status.success());
+        });
+    });
+
+    group.bench_function("granges_filter", |b| {
+        b.iter(|| {
+            let granges_output = Command::new(granges_binary_path())
+                .arg("filter")
+                .arg("--seqlens")
+                .arg("tests_data/hg38_seqlens.tsv")
+                .arg("--left")
+                .arg(&random_bedfile_left)
+                .arg("--right")
+                .arg(&random_bedfile_right)
+                .output()
+                .expect("granges adjust failed");
+            assert!(granges_output.status.success());
+        });
+    });
+}
+
+criterion_group!(benches, bench_filter_adjustment, bench_range_adjustment);
 criterion_main!(benches);

src/commands.rs

@@ -1,7 +1,7 @@
 use std::path::PathBuf;
 
 use crate::{
-    io::{OutputFile, parsers::GenomicRangesParser},
+    io::{parsers::GenomicRangesParser, OutputFile},
     prelude::*,
     ranges::operations::adjust_range,
     reporting::{CommandOutput, Report},
@@ -56,34 +56,33 @@ pub fn granges_adjust(
 
             if skipped_ranges > 0 {
                 report.add_issue(format!(
-                        "{} ranges were removed because their widths after adjustment were ≤ 0",
-                        skipped_ranges
-                        ))
+                    "{} ranges were removed because their widths after adjustment were ≤ 0",
+                    skipped_ranges
+                ))
             }
         }
     } else {
-        // If we do need to sort, build up a GRanges variant and adjust ranges through 
-        // the GRanges interface. Note we need to detect and build a specific iterator 
+        // If we do need to sort, build up a GRanges variant and adjust ranges through
+        // the GRanges interface. Note we need to detect and build a specific iterator
         // for the filetype.
 
         let ranges_iter = GenomicRangesFile::parsing_iterator(bedfile)?;
         match ranges_iter {
             GenomicRangesParser::Bed3(iter) => {
                 let gr = GRangesEmpty::from_iter(iter, &genome)?;
                 gr.adjust_ranges(-both, both).to_tsv(output)?
-            },
+            }
             GenomicRangesParser::Bedlike(iter) => {
                 // Note the call to try_unwrap_data() here: this is because
                 // we know that the records *do* have data. Unwrapping the Option<String>
-                // values means that writing to TSV doesn't have to deal with this (which 
+                // values means that writing to TSV doesn't have to deal with this (which
                 // always creates headaches).
                 let gr = GRanges::from_iter(iter.try_unwrap_data()?, &genome)?;
                 gr.adjust_ranges(-both, both).to_tsv(output)?
-            },
+            }
             GenomicRangesParser::Unsupported => {
                 return Err(GRangesError::UnsupportedGenomicRangesFileFormat)
-            },
-
+            }
         }
     }
     Ok(CommandOutput::new((), report))
@@ -102,47 +101,65 @@ pub fn granges_filter(
     let left_iter = GenomicRangesFile::parsing_iterator(left_path)?;
     let right_iter = GenomicRangesFile::parsing_iterator(right_path)?;
 
-    let output_stream = output.as_ref().map_or(OutputFile::new_stdout(None), |file| {
-        OutputFile::new(file, None)
-    });
-    let mut writer = output_stream.writer()?;
-
     // for reporting stuff to the user
     let mut report = Report::new();
 
     match (left_iter, right_iter) {
         (GenomicRangesParser::Bed3(left), GenomicRangesParser::Bed3(right)) => {
+            let left_gr = GRangesEmpty::from_iter(left, &genome)?;
+            let right_gr = GRangesEmpty::from_iter(right, &genome)?;
 
-            let left_gr = GRangesEmpty::from_iter(left, &genome)?; let right_gr =
-                GRangesEmpty::from_iter(right, &genome)?;
-
-            let right_gr = right_gr.to_coitrees()?;
+            let right_gr = right_gr.into_coitrees()?;
 
-            let intersection = left_gr.filter_overlaps(&right_gr)?;
+            let mut intersection = left_gr.filter_overlaps(&right_gr)?;
             intersection.to_tsv(output)?;
 
             Ok(CommandOutput::new((), report))
         }
         (GenomicRangesParser::Bed3(left), GenomicRangesParser::Bedlike(right)) => {
+            let left_gr = GRangesEmpty::from_iter(left, &genome)?;
+            let right_gr = GRanges::from_iter(right.try_unwrap_data()?, &genome)?;
+
+            let right_gr = right_gr.into_coitrees()?;
+
+            let intersection = left_gr.filter_overlaps(&right_gr)?;
+            intersection.to_tsv(output)?;
+
             Ok(CommandOutput::new((), report))
         }
         (GenomicRangesParser::Bedlike(left), GenomicRangesParser::Bed3(right)) => {
+            let left_gr = GRanges::from_iter(left.try_unwrap_data()?, &genome)?;
+            let right_gr = GRangesEmpty::from_iter(right, &genome)?;
+
+            let right_gr = right_gr.into_coitrees()?;
+
+            let intersection = left_gr.filter_overlaps(&right_gr)?;
+            intersection.to_tsv(output)?;
+
             Ok(CommandOutput::new((), report))
         }
         (GenomicRangesParser::Bedlike(left), GenomicRangesParser::Bedlike(right)) => {
+            let left_gr = GRanges::from_iter(left.try_unwrap_data()?, &genome)?;
+            let right_gr = GRanges::from_iter(right.try_unwrap_data()?, &genome)?;
+
+            let right_gr = right_gr.into_coitrees()?;
+
+            let intersection = left_gr.filter_overlaps(&right_gr)?;
+            intersection.to_tsv(output)?;
+
             Ok(CommandOutput::new((), report))
         }
-        _ => Ok(CommandOutput::new((), report)),
+        _ => Err(GRangesError::UnsupportedGenomicRangesFileFormat),
     }
 }
 
 /// Generate a random BED-like file with genomic ranges.
 pub fn granges_random_bed(
     seqlens: impl Into<PathBuf>,
-    num: u32,
+    num: usize,
     output: Option<impl Into<PathBuf>>,
     sort: bool,
-    ) -> Result<CommandOutput<()>, GRangesError> {
+) -> Result<CommandOutput<()>, GRangesError> {
     // get the genome info
     let genome = read_seqlens(seqlens)?;