GenomicRangeRecordUnwrappable trait added, some renaming, better CLI …

…ergo.

 - Changed --seqlens to --genome.
 - Added short arguments.
 - Added skip_missing to skip over ranges not in the genome file.
 - skip_missing required new filtering, so introduced the
   `GenomicRangeRecordUnwrappable` trait. This implements the
   `try_unwrap_data()` method (which was removed from `BedlikeIterator`.
   The previous method in `BedlikeIterator` returned `impl Iterator`
   which since it isn't a concrete type, caused lots of ergonomic issues.
 - Introduced new `UnwrappedRanges` type for the above.
 - Removed `into_variants()`.
 - Implemented `retain_seqnames()` and `exclude_seqnames()` for `UnwrappedRanges`.

benches/bedtools_comparison.rs

@@ -41,7 +41,7 @@ fn bench_range_adjustment(c: &mut Criterion) {
  b.iter(|| {
  let granges_output = Command::new(granges_binary_path())
  .arg("adjust")
- .arg("--seqlens")
+ .arg("--genome")
  .arg("tests_data/hg38_seqlens.tsv")
  .arg("--both")
  .arg("10")
@@ -86,7 +86,7 @@ fn bench_filter_adjustment(c: &mut Criterion) {
  b.iter(|| {
  let granges_output = Command::new(granges_binary_path())
  .arg("filter")
- .arg("--seqlens")
+ .arg("--genome")
  .arg("tests_data/hg38_seqlens.tsv")
  .arg("--left")
  .arg(&random_bedfile_left)

src/commands.rs

@@ -77,7 +77,7 @@ pub fn granges_adjust(
  // 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)?;
+ let gr = GRanges::from_iter(iter.try_unwrap_data(), &genome)?;
  gr.adjust_ranges(-both, both).to_tsv(output)?
  }
  GenomicRangesParser::Unsupported => {
@@ -94,9 +94,11 @@ pub fn granges_filter(
  left_path: &PathBuf,
  right_path: &PathBuf,
  output: Option<&PathBuf>,
+ skip_missing: bool,
  sort: bool,
 ) -> Result<CommandOutput<()>, GRangesError> {
  let genome = read_seqlens(seqlens)?;
+ let seqnames: Vec<String> = genome.keys().cloned().collect();
 
  let left_iter = GenomicRangesFile::parsing_iterator(left_path)?;
  let right_iter = GenomicRangesFile::parsing_iterator(right_path)?;
@@ -106,19 +108,38 @@ pub fn granges_filter(
 
  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;
+ let right_gr;
+
+ if skip_missing {
+ left_gr = GRangesEmpty::from_iter(left.retain_seqnames(&seqnames), &genome)?;
+ right_gr = GRangesEmpty::from_iter(right.retain_seqnames(&seqnames), &genome)?;
+ } else {
+ left_gr = GRangesEmpty::from_iter(left, &genome)?;
+ right_gr = GRangesEmpty::from_iter(right, &genome)?;
+ }
 
  let right_gr = right_gr.into_coitrees()?;
 
- let mut intersection = left_gr.filter_overlaps(&right_gr)?;
+ let 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 left_gr;
+ let right_gr;
+
+ if skip_missing {
+ left_gr = GRangesEmpty::from_iter(left.retain_seqnames(&seqnames), &genome)?;
+ right_gr = GRanges::from_iter(
+ right.try_unwrap_data().retain_seqnames(&seqnames),
+ &genome,
+ )?;
+ } else {
+ left_gr = GRangesEmpty::from_iter(left, &genome)?;
+ right_gr = GRanges::from_iter(right.try_unwrap_data(), &genome)?;
+ }
 
  let right_gr = right_gr.into_coitrees()?;
 
@@ -128,8 +149,17 @@ pub fn granges_filter(
  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 left_gr;
+ let right_gr;
+
+ if skip_missing {
+ left_gr =
+ GRanges::from_iter(left.try_unwrap_data().retain_seqnames(&seqnames), &genome)?;
+ right_gr = GRangesEmpty::from_iter(right.retain_seqnames(&seqnames), &genome)?;
+ } else {
+ left_gr = GRanges::from_iter(left.try_unwrap_data(), &genome)?;
+ right_gr = GRangesEmpty::from_iter(right, &genome)?;
+ }
 
  let right_gr = right_gr.into_coitrees()?;
 
@@ -139,8 +169,20 @@ pub fn granges_filter(
  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 left_gr;
+ let right_gr;
+
+ if skip_missing {
+ left_gr =
+ GRanges::from_iter(left.try_unwrap_data().retain_seqnames(&seqnames), &genome)?;
+ right_gr = GRanges::from_iter(
+ right.try_unwrap_data().retain_seqnames(&seqnames),
+ &genome,
+ )?;
+ } else {
+ left_gr = GRanges::from_iter(left.try_unwrap_data(), &genome)?;
+ right_gr = GRanges::from_iter(right.try_unwrap_data(), &genome)?;
+ }
 
  let right_gr = right_gr.into_coitrees()?;
 
@@ -163,7 +205,7 @@ pub fn granges_random_bed(
  // get the genome info
  let genome = read_seqlens(seqlens)?;
 
- let mut gr = random_granges(&genome, num.try_into().unwrap())?;
+ let mut gr = random_granges(&genome, num)?;
 
  if sort {
  gr = gr.sort();

src/error.rs

@@ -34,7 +34,7 @@ pub enum GRangesError {
  InvalidGenomicRange(Position, Position),
  #[error("Range [{0}, {1}] is invalid for sequence of length {2}")]
  InvalidGenomicRangeForSequence(Position, Position, Position),
- #[error("Sequence name '{0}' is not the ranges container")]
+ #[error("Sequence name '{0}' is not in the ranges container")]
  MissingSequence(String),
  #[error("Error encountered in genomap::GenomeMap")]
  GenomeMapError(#[from] GenomeMapError),

src/io/parsers.rs

@@ -77,7 +77,7 @@ use std::path::{Path, PathBuf};
 use crate::error::GRangesError;
 use crate::io::file::InputFile;
 use crate::ranges::{GenomicRangeEmptyRecord, GenomicRangeRecord};
-use crate::traits::GeneralRangeRecordIterator;
+use crate::traits::{GeneralRangeRecordIterator, GenomicRangeRecordUnwrappable};
 use crate::Position;
 
 /// Get the *base* extension to help infer filetype, which ignores compression-related
@@ -365,37 +365,10 @@ impl BedlikeIterator {
  self.number_columns() == 3
  }
 
- /// Try to unwrap each [`GenomicRangeRecord<Option<String>>`] iterator item into a
- /// [`GenomicRangeRecord<String>`] iterator item.
- ///
- /// This will raise errors if:
- /// 1. The detected number of columns is < 4 (i.e. there appears to be no data to unwrap).
- /// 2. During iteration, a `None` data element is encountered.
- pub fn try_unwrap_data(
- self,
- ) -> Result<impl Iterator<Item = Result<GenomicRangeRecord<String>, GRangesError>>, GRangesError>
- {
- if self.number_columns() < 4 {
- return Err(GRangesError::TooFewColumns)?;
- }
- Ok(self.iter.map(|result| {
- result.and_then(|record| {
- if let Some(data) = record.data {
- Ok(GenomicRangeRecord::new(
- record.seqname,
- record.start,
- record.end,
- data,
- ))
- } else {
- Err(GRangesError::TryUnwrapDataError)
- }
- })
- }))
- }
-
  /// Drop the data in each [`GenomicRangeRecord<Option<String>>`] iterator, converting it to a range-only
  /// [`GenomicRangeRecord<()>`] iterator item.
+ // TODO: candidate for trait? the impl Iterator isn't a concrete type and can be cumbersome
+ // in terms of ergonomics. See UnwrappedRanges for an example.
  pub fn drop_data(self) -> impl Iterator<Item = Result<GenomicRangeRecord<()>, GRangesError>> {
  self.iter.map(|result| {
  result
@@ -410,37 +383,6 @@ impl BedlikeIterator {
  .unwrap_or_else(Err) // pass through parsing errors
  })
  }
-
- /// Try to convert the iterator into one of two variants: one with data and one without.
- pub fn into_variant(self) -> Result<GenomicRangesIteratorVariant, GRangesError> {
- let number_columns = self.number_columns();
-
- if number_columns == 3 {
- let without_data_iterator = self.drop_data().map(|result| {
- result.map(|record| GenomicRangeRecord {
- seqname: record.seqname,
- start: record.start,
- end: record.end,
- data: (),
- })
- });
- Ok(GenomicRangesIteratorVariant::Empty(Box::new(
- without_data_iterator,
- )))
- } else {
- let with_data_iterator = self.try_unwrap_data()?.map(|result| {
- result.map(|record| GenomicRangeRecord {
- seqname: record.seqname,
- start: record.start,
- end: record.end,
- data: record.data,
- })
- });
- Ok(GenomicRangesIteratorVariant::WithData(Box::new(
- with_data_iterator,
- )))
- }
- }
 }
 
 impl Iterator for BedlikeIterator {
@@ -465,7 +407,7 @@ impl Iterator for BedlikeIterator {
 ///
 /// let iter = Bed3Iterator::new("tests_data/example.bed")
 /// .expect("error reading file")
-/// .exclude_seqnames(vec!["chr1".to_string()]);
+/// .exclude_seqnames(&vec!["chr1".to_string()]);
 ///
 /// let seqlens = seqlens! { "chr1" => 22, "chr2" => 10, "chr3" => 10, "chr4" => 15 };
 /// let gr = GRangesEmpty::from_iter(iter, &seqlens)
@@ -479,6 +421,7 @@ impl Iterator for BedlikeIterator {
 /// assert_eq!(iter.next().unwrap().end, 15);
 /// assert_eq!(iter.next(), None);
 /// ```
+#[derive(Debug)]
 pub struct FilteredRanges<I, R>
 where
  I: Iterator<Item = Result<R, GRangesError>>,
@@ -494,11 +437,11 @@ where
 {
  pub fn new(
  inner: I,
- retain_seqnames: Option<Vec<String>>,
- exclude_seqnames: Option<Vec<String>>,
+ retain_seqnames: Option<&Vec<String>>,
+ exclude_seqnames: Option<&Vec<String>>,
  ) -> Self {
- let retain_seqnames = retain_seqnames.map(HashSet::from_iter);
- let exclude_seqnames = exclude_seqnames.map(HashSet::from_iter);
+ let retain_seqnames = retain_seqnames.cloned().map(HashSet::from_iter);
+ let exclude_seqnames = exclude_seqnames.cloned().map(HashSet::from_iter);
  Self {
  inner,
  retain_seqnames,
@@ -575,18 +518,93 @@ where
  }
 }
 
-impl GeneralRangeRecordIterator<GenomicRangeEmptyRecord> for Bed3Iterator {
+impl GeneralRangeRecordIterator<GenomicRangeRecord<Option<String>>> for BedlikeIterator {
  fn retain_seqnames(
  self,
- seqnames: Vec<String>,
- ) -> FilteredRanges<Self, GenomicRangeEmptyRecord> {
- FilteredRanges::new(self, Some(seqnames), None)
+ seqnames: &[String],
+ ) -> FilteredRanges<Self, GenomicRangeRecord<Option<String>>> {
+ FilteredRanges::new(self, Some(&seqnames.to_vec()), None)
  }
  fn exclude_seqnames(
  self,
- seqnames: Vec<String>,
+ seqnames: &[String],
+ ) -> FilteredRanges<Self, GenomicRangeRecord<Option<String>>> {
+ FilteredRanges::new(self, None, Some(&seqnames.to_vec()))
+ }
+}
+
+impl GeneralRangeRecordIterator<GenomicRangeEmptyRecord> for Bed3Iterator {
+ fn retain_seqnames(self, seqnames: &[String]) -> FilteredRanges<Self, GenomicRangeEmptyRecord> {
+ FilteredRanges::new(self, Some(&seqnames.to_vec()), None)
+ }
+ fn exclude_seqnames(
+ self,
+ seqnames: &[String],
  ) -> FilteredRanges<Self, GenomicRangeEmptyRecord> {
- FilteredRanges::new(self, None, Some(seqnames))
+ FilteredRanges::new(self, None, Some(&seqnames.to_vec()))
+ }
+}
+
+impl<I> GeneralRangeRecordIterator<GenomicRangeRecord<String>> for UnwrappedRanges<I>
+where
+ I: Iterator<Item = Result<GenomicRangeRecord<Option<String>>, GRangesError>>,
+{
+ fn retain_seqnames(
+ self,
+ seqnames: &[String],
+ ) -> FilteredRanges<Self, GenomicRangeRecord<String>> {
+ FilteredRanges::new(self, Some(&seqnames.to_vec()), None)
+ }
+ fn exclude_seqnames(
+ self,
+ seqnames: &[String],
+ ) -> FilteredRanges<Self, GenomicRangeRecord<String>> {
+ FilteredRanges::new(self, None, Some(&seqnames.to_vec()))
+ }
+}
+
+#[derive(Debug)]
+pub struct UnwrappedRanges<I>
+where
+ I: Iterator<Item = Result<GenomicRangeRecord<Option<String>>, GRangesError>>,
+{
+ inner: I,
+}
+
+impl<I> UnwrappedRanges<I>
+where
+ I: Iterator<Item = Result<GenomicRangeRecord<Option<String>>, GRangesError>>,
+{
+ pub fn new(inner: I) -> Self {
+ Self { inner }
+ }
+}
+
+impl<I> Iterator for UnwrappedRanges<I>
+where
+ I: Iterator<Item = Result<GenomicRangeRecord<Option<String>>, GRangesError>>,
+{
+ type Item = Result<GenomicRangeRecord<String>, GRangesError>;
+
+ fn next(&mut self) -> Option<Self::Item> {
+ if let Some(result) = self.inner.next() {
+ return Some(result.and_then(|record| match record.data {
+ Some(data) => Ok(GenomicRangeRecord::new(
+ record.seqname,
+ record.start,
+ record.end,
+ data,
+ )),
+ None => Err(GRangesError::TryUnwrapDataError),
+ }));
+ }
+ None
+ }
+}
+
+impl GenomicRangeRecordUnwrappable for BedlikeIterator {
+ fn try_unwrap_data(self) -> UnwrappedRanges<Self> {
+ UnwrappedRanges::new(self)
  }
 }
 

src/lib.rs

@@ -32,8 +32,9 @@ pub mod prelude {
 
  pub use crate::ranges::vec::{VecRangesEmpty, VecRangesIndexed};
  pub use crate::traits::{
- AsGRangesRef, GeneralRangeRecordIterator, GenericRange, GenomicRangesTsvSerialize,
- IndexedDataContainer, IntoIterableRangesContainer, IterableRangeContainer, TsvSerialize,
+ AsGRangesRef, GeneralRangeRecordIterator, GenericRange, GenomicRangeRecordUnwrappable,
+ GenomicRangesTsvSerialize, IndexedDataContainer, IntoIterableRangesContainer,
+ IterableRangeContainer, TsvSerialize,
  };
 
  pub use crate::seqlens;