Sequences::region_map_into_granges(), GRanges::map_into_array*() …

…added, etc.

 - Renamed `region_apply` in `region_map` to be more consistent with Rust's naming.
 - `retain_seqnames`, etc for `TsvRecordIterator`.
 - New `Sequences::region_map_into_granges()`, which applies a function to
   regions of a sequence, and puts the results in a new `GRanges`.
 - New `Sequences` test cases corresponding to `granges_test_case_01()`,
   and lots of new `Sequence`-related tests.
 - Added `map_into_array1()` and `map_into_array2()` for the
   --feature=ndarray, with tests.
 - `take_ranges()` and `take_both()` added.
 - New `GRanges::midpoints()` method, with tests.
 - New `GRanges::data_indices()` method, to build a `GenomeMap` out of
   the data indices.
 - New `GRanges::data_by_seqnames()` which builds a `GenomeMap` of
   `Vec<U>` of the data, by sequence name. `GRanges::data_refs_by_seqnames()`
   also added, for references.
 - `GRanges::into_vecranegs()` added, with `From` method
   for `COITreesEmpty` to `VecRangesEmpty` added since it's needed.
 - Cleaned lifetime out of `IndexedDataContainer`, brought down
   to assoc. type level. Cleans up things a lot!
 - `into_array1()` and `into_array2()` and tests.
 - Remove section on readme on slow benchmarks -- this is now
   GH issue #4.

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "granges"
-version = "0.1.0"
+version = "0.2.0"
 edition = "2021"
 license = "MIT"
 authors = ["Vince Buffalo <[email protected]>"]
@@ -31,9 +31,8 @@ csv = "1.3.0"
 serde = { version = "1.0.197", features = ["derive"] }
 
 [features]
-# cli = [ "clap" ] # TODO make feature
 dev-commands = [ ]
-test_bigranges = [] # NOTE: not used yet, for tests on large files
+bench-big = []
 polars = ["dep:polars"]
 ndarray = ["dep:ndarray", "dep:ndarray-npy"]
 big-position = []

README.md

@@ -110,23 +110,3 @@ windows 418.87 s 65.13 s 84.4515
 map_median 112.35 s 82.81 s 26.2995
 ```
 
-Note however, there is a arithmetic scaling issue. Larger range datasets
-(1,000,000) lead to map operations that are inefficient. This is almost surely 
-due to `data/operations.rs`. Here is an example benchmark:
-
-```
-command bedtools time granges time granges speedup (%)
------------- --------------- -------------- ---------------------
-map_multiple 26.89 min 688.60 s 57.3189
-map_max 21.34 min 21.54 min -0.94889 *
-adjust 21.35 min 499.54 s 60.9959
-filter 27.01 min 20.38 min 24.5583
-map_min 935.69 s 17.48 min -12.1008 * 
-flank 30.26 min 943.36 s 48.0353
-map_mean 20.13 min 931.21 s 22.913
-map_sum 17.99 min 18.44 min -2.5225 *
-windows 507.84 s 84.79 s 83.304
-map_median 16.90 min 17.81 min -5.36811 *
-```
-
-So median, sum, min, and max are slower.

benches/bedtools_comparison.rs

@@ -13,7 +13,10 @@ use std::{
  process::{Command, Stdio},
 };
 
+#[cfg(not(feature = "bench-big"))]
 const BED_LENGTH: usize = 100_000;
+#[cfg(feature = "bench-big")]
+const BED_LENGTH: usize = 1_000_000;
 
 fn bench_range_adjustment(c: &mut Criterion) {
  // create the benchmark group
@@ -229,7 +232,7 @@ fn bench_map(c: &mut Criterion) {
  let mut group = c.benchmark_group(format!("map_{}", operation).to_string());
 
  // configure the sample size for the group
- group.sample_size(10);
+ group.sample_size(100);
  group.bench_function("bedtools_map", |b| {
  b.iter(|| {
  let bedtools_output_file = File::create(&bedtools_path).unwrap();

examples/mean_score.rs

@@ -35,11 +35,9 @@ fn try_main() -> Result<(), granges::error::GRangesError> {
  .map_data(|bed5_cols| bed5_cols.score)?;
 
  // Compute overlaps and combine scores into mean.
- let results_gr = left_gr
- .left_overlaps(&right_gr)?
- .map_joins(mean_score)?;
+ let results_gr = left_gr.left_overlaps(&right_gr)?.map_joins(mean_score)?;
 
- results_gr.to_tsv(None::<String>, &BED_TSV)?;
+ results_gr.write_to_tsv(None::<String>, &BED_TSV)?;
  Ok(())
 }
 

src/commands.rs

@@ -105,19 +105,22 @@ pub fn granges_adjust(
  match ranges_iter {
  GenomicRangesParser::Bed3(iter) => {
  let gr = GRangesEmpty::from_iter(iter, &genome)?;
- gr.adjust_ranges(-both, both).to_tsv(output, &BED_TSV)?
+ gr.adjust_ranges(-both, both)
+ .write_to_tsv(output, &BED_TSV)?
  }
  GenomicRangesParser::Bed5(iter) => {
  let gr = GRanges::from_iter(iter, &genome)?;
- gr.adjust_ranges(-both, both).to_tsv(output, &BED_TSV)?
+ gr.adjust_ranges(-both, both)
+ .write_to_tsv(output, &BED_TSV)?
  }
  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
  // always creates headaches).
  let gr = GRanges::from_iter(iter.try_unwrap_data(), &genome)?;
- gr.adjust_ranges(-both, both).to_tsv(output, &BED_TSV)?
+ gr.adjust_ranges(-both, both)
+ .write_to_tsv(output, &BED_TSV)?
  }
  GenomicRangesParser::Unsupported => {
  return Err(GRangesError::UnsupportedGenomicRangesFileFormat)
@@ -181,7 +184,7 @@ pub fn granges_filter(
  let right_gr = right_gr.into_coitrees()?;
 
  let semijoin = left_gr.filter_overlaps(&right_gr)?;
- semijoin.to_tsv(output, &BED_TSV)?;
+ semijoin.write_to_tsv(output, &BED_TSV)?;
 
  Ok(CommandOutput::new((), report))
  }
@@ -203,7 +206,7 @@ pub fn granges_filter(
  let right_gr = right_gr.into_coitrees()?;
 
  let semijoin = left_gr.filter_overlaps(&right_gr)?;
- semijoin.to_tsv(output, &BED_TSV)?;
+ semijoin.write_to_tsv(output, &BED_TSV)?;
 
  Ok(CommandOutput::new((), report))
  }
@@ -223,7 +226,7 @@ pub fn granges_filter(
  let right_gr = right_gr.into_coitrees()?;
 
  let semijoin = left_gr.filter_overlaps(&right_gr)?;
- semijoin.to_tsv(output, &BED_TSV)?;
+ semijoin.write_to_tsv(output, &BED_TSV)?;
 
  Ok(CommandOutput::new((), report))
  }
@@ -246,7 +249,7 @@ pub fn granges_filter(
  let right_gr = right_gr.into_coitrees()?;
 
  let intersection = left_gr.filter_overlaps(&right_gr)?;
- intersection.to_tsv(output, &BED_TSV)?;
+ intersection.write_to_tsv(output, &BED_TSV)?;
 
  Ok(CommandOutput::new((), report))
  }
@@ -303,7 +306,8 @@ pub fn granges_flank(
  } else {
  GRangesEmpty::from_iter(iter, &genome)?
  };
- gr.flanking_ranges(left, right)?.to_tsv(output, &BED_TSV)?
+ gr.flanking_ranges(left, right)?
+ .write_to_tsv(output, &BED_TSV)?
  }
  GenomicRangesParser::Bed5(_iter) => {
  unimplemented!()
@@ -314,7 +318,8 @@ pub fn granges_flank(
  } else {
  GRanges::from_iter(iter.try_unwrap_data(), &genome)?
  };
- gr.flanking_ranges(left, right)?.to_tsv(output, &BED_TSV)?
+ gr.flanking_ranges(left, right)?
+ .write_to_tsv(output, &BED_TSV)?
  }
  GenomicRangesParser::Unsupported => {
  return Err(GRangesError::UnsupportedGenomicRangesFileFormat)
@@ -469,7 +474,7 @@ pub fn granges_map(
  .collect::<Vec<DatumType>>()
  })?;
 
- result_gr.to_tsv(output, &BED_TSV)?;
+ result_gr.write_to_tsv(output, &BED_TSV)?;
 
  Ok(CommandOutput::new((), report))
 }
@@ -483,7 +488,7 @@ pub fn granges_windows(
  output: Option<impl Into<PathBuf>>,
 ) -> Result<CommandOutput<()>, GRangesError> {
  let genome = read_seqlens(seqlens)?;
- GRangesEmpty::from_windows(&genome, width, step, chop)?.to_tsv(output, &BED_TSV)?;
+ GRangesEmpty::from_windows(&genome, width, step, chop)?.write_to_tsv(output, &BED_TSV)?;
  let report = Report::new();
  Ok(CommandOutput::new((), report))
 }
@@ -504,13 +509,13 @@ pub fn granges_random_bed(
  if sort {
  gr = gr.sort()
  }
- gr.to_tsv(output, &BED_TSV)?;
+ gr.write_to_tsv(output, &BED_TSV)?;
  } else {
  let mut gr = random_granges(&genome, num)?;
  if sort {
  gr = gr.sort();
  }
- gr.to_tsv(output, &BED_TSV)?;
+ gr.write_to_tsv(output, &BED_TSV)?;
  };
 
  let report = Report::new();

src/data/mod.rs

@@ -3,9 +3,14 @@
 
 use crate::traits::{DataContainer, IntoDatumType};
 
+#[cfg(feature = "ndarray")]
+pub mod ndarray;
 pub mod operations;
 pub mod vec;
 
+impl<U> DataContainer for Vec<U> {}
+impl DataContainer for () {}
+
 /// These are core supported data types stored in an `enum`, to
 /// unify the types that come out of standard operations like
 /// `select()`.
@@ -67,6 +72,3 @@ impl<T: IntoDatumType> From<T> for DatumType {
  item.into_data_type()
  }
 }
-
-impl<U> DataContainer for Vec<U> {}
-impl DataContainer for () {}