diff --git a/Cargo.toml b/Cargo.toml index 9e5112a..ab9cf5b 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "sketchlib" -version = "0.4.1" +version = "0.4.2" authors = [ "John Lees ", "Nicholas Croucher ", diff --git a/src/hashing/aahash_iterator.rs b/src/hashing/aahash_iterator.rs index 5f87413..408a4a4 100644 --- a/src/hashing/aahash_iterator.rs +++ b/src/hashing/aahash_iterator.rs @@ -30,16 +30,15 @@ pub struct AaHashIterator { } impl RollHash for AaHashIterator { - fn set_k(&mut self, k: usize) { + fn set_k(&mut self, k: usize) -> anyhow::Result<()> { self.k = k; if let Some(new_it) = Self::new_iterator(0, &self.level, &self.seq, k) { self.fh = new_it.0; self.index = new_it.1; + Ok(()) } else { - panic!( - "K-mer larger than smallest valid sequence, which is:\n{}", - std::str::from_utf8(&self.seq).unwrap() - ); + let seq_str = std::str::from_utf8(&self.seq).unwrap_or(""); + Err(anyhow::anyhow!("K-mer size {} larger than smallest valid sequence {}", k, seq_str)) } } diff --git a/src/hashing/mod.rs b/src/hashing/mod.rs index e4cb14a..04a2359 100644 --- a/src/hashing/mod.rs +++ b/src/hashing/mod.rs @@ -114,7 +114,7 @@ fn swapbits3263(v: u64) -> u64 { /// Rolling functions supported by both ntHash and aaHash pub trait RollHash: Iterator { /// Set the k-mer size - fn set_k(&mut self, k: usize); + fn set_k(&mut self, k: usize) -> anyhow::Result<()>; /// Get the current hash fn curr_hash(&self) -> u64; /// The type of sequence being hashed diff --git a/src/hashing/nthash_iterator.rs b/src/hashing/nthash_iterator.rs index adf48ea..da0664e 100644 --- a/src/hashing/nthash_iterator.rs +++ b/src/hashing/nthash_iterator.rs @@ -46,14 +46,16 @@ pub struct NtHashIterator { } impl RollHash for NtHashIterator { - fn set_k(&mut self, k: usize) { + fn set_k(&mut self, k: usize) -> anyhow::Result<()> { if k != self.k { self.k = k; self.offset_idx = 0; // rewind: offsets must be re-traversed for each k if self.next_iterator(0).is_none() { - panic!("K-mer larger than smallest valid sequence"); + let seq_str = std::str::from_utf8(&self.seq).unwrap_or(""); + return Err(anyhow::anyhow!("K-mer size {} larger than smallest valid sequence {}", k, seq_str)) } } + Ok(()) } /// Retrieve the current hash (minimum of forward and reverse complement hashes) @@ -95,7 +97,7 @@ impl NtHashIterator { rc: bool, min_qual: u8, reads: bool, - ) -> Vec { + ) -> anyhow::Result> { let mut seq = Vec::new(); let mut offsets = Vec::new(); let mut acgt = [0, 0, 0, 0]; @@ -129,8 +131,8 @@ impl NtHashIterator { non_acgt, reads, }; - hash_it.set_k(k); - vec![hash_it] + hash_it.set_k(k)?; + Ok(vec![hash_it]) } #[cfg(target_arch = "wasm32")] @@ -448,7 +450,7 @@ impl NtHashIterator { non_acgt, reads: false, }; - hash_it.set_k(k); + hash_it.set_k(k).unwrap(); hash_it } } @@ -633,9 +635,9 @@ mod tests { // Simulate the sketch command: one iterator, multiple set_k calls in sequence. let mut it = NtHashIterator::from_seq(seq, 3, true); let hashes_k3: Vec = it.by_ref().collect(); - it.set_k(5); + it.set_k(5).unwrap(); let hashes_k5: Vec = it.by_ref().collect(); - it.set_k(7); + it.set_k(7).unwrap(); let hashes_k7: Vec = it.by_ref().collect(); assert_eq!(hashes_k3, ref_hashes(seq, 3, true), "k=3 hashes wrong"); assert_eq!( diff --git a/src/inverted.rs b/src/inverted.rs index 502e63d..0870cbb 100644 --- a/src/inverted.rs +++ b/src/inverted.rs @@ -379,7 +379,7 @@ impl Inverted { crate::io::NeedletailIterator::new(reader) }).collect::>(); - NtHashIterator::new(&mut records_readers, k, rc, min_qual, reads) + NtHashIterator::new(&mut records_readers, k, rc, min_qual, reads).unwrap() .into_iter() .map(|it| Box::new(it) as Box) .collect() @@ -428,18 +428,18 @@ impl Inverted { // Not yet written! if !multientrysamples.contains(name) { // Densifying now! - Sketch::densify_bin(&mut sketch); + Sketch::densify_bin(sketch.as_mut().unwrap()); } else { // We'll need to densify afterwards, let's save the index indexes.insert(genome_idx); } - sketch_results[genome_idx] = sketch.iter().map(|h| *h as u16).collect(); + sketch_results[genome_idx] = sketch.as_ref().unwrap().iter().map(|h| *h as u16).collect(); differentsamples.remove(name); } else { // already written! We have to merge for bin in 0..sketch_size { let saved_sketch = &mut sketch_results[genome_idx][bin as usize]; - *saved_sketch = cmp::min(*saved_sketch, sketch[bin as usize] as u16); + *saved_sketch = cmp::min(*saved_sketch, sketch.as_ref().unwrap()[bin as usize] as u16); } } } @@ -496,7 +496,7 @@ impl Inverted { }; let (signs, densified) = - Sketch::get_signs(&mut **hash_it, k, &mut read_filter, sketch_size); + Sketch::get_signs(&mut **hash_it, k, &mut read_filter, sketch_size).unwrap(); if densified { logw("The query was densified", Some("trace")); } diff --git a/src/sketch/mod.rs b/src/sketch/mod.rs index e799a9e..573c0f2 100644 --- a/src/sketch/mod.rs +++ b/src/sketch/mod.rs @@ -171,7 +171,7 @@ impl Sketch { sketch_size: u64, rc: bool, min_count: u16, - ) -> Self { + ) -> anyhow::Result { let (_sketchsize64, num_bins, usigs_size) = num_bins(sketch_size); let flattened_size_u64 = usigs_size as usize * kmer_lengths.len(); let mut usigs = aligned_sketch_vec_with_capacity(flattened_size_u64); @@ -189,7 +189,7 @@ impl Sketch { let mut densified = false; for k in kmer_lengths { log::debug!("Running sketching at k={k}"); - let (signs, k_densified) = Self::get_signs(seq_hashes, *k, &mut read_filter, num_bins); + let (signs, k_densified) = Self::get_signs(seq_hashes, *k, &mut read_filter, num_bins)?; densified |= k_densified; minhash_sum += (signs[0] as f64) / (u64::MAX as f64); @@ -208,7 +208,7 @@ impl Sketch { seq_hashes.seq_len() }; - Self { + Ok(Self { usigs, name: name.to_string(), index: None, @@ -218,7 +218,7 @@ impl Sketch { densified, acgt, non_acgt, - } + }) } /// Get the sketch bins for a sample, but do not transpose @@ -227,13 +227,13 @@ impl Sketch { kmer_size: usize, filter: &mut Option, num_bins: u64, - ) -> (Vec, bool) { + ) -> anyhow::Result<(Vec, bool)> { // Setup storage for each k let mut signs = vec![u64::MAX; num_bins as usize]; if let Some(read_filter) = filter { read_filter.clear(); } - seq_hashes.set_k(kmer_size); + seq_hashes.set_k(kmer_size)?; // Calculate bin minima across all sequence for hash in seq_hashes.iter() { @@ -241,7 +241,7 @@ impl Sketch { } // Densify let densified = Self::densify_bin(&mut signs); - (signs, densified) + Ok((signs, densified)) } /// Get the sketch bins for a sample, but do not transpose @@ -250,20 +250,20 @@ impl Sketch { kmer_size: usize, filter: &mut Option, num_bins: u64, - ) -> Vec { + ) -> anyhow::Result> { // Setup storage for each k let mut signs = vec![u64::MAX; num_bins as usize]; if let Some(read_filter) = filter { read_filter.clear(); } - seq_hashes.set_k(kmer_size); + seq_hashes.set_k(kmer_size)?; // Calculate bin minima across all sequence for hash in seq_hashes.iter() { Self::bin_sign(&mut signs, hash, num_bins, filter); } - signs + Ok(signs) } /// The name of the sample @@ -467,7 +467,7 @@ impl fmt::Display for Sketch { /// let reader = needletail::parse_fastx_file(fastx_path).unwrap(); /// let mut filtered_iters = vec![NeedletailFilterIterator::new(reader, want_ids)]; /// -/// sketch_data(&mut filtered_iters, opts) +/// sketch_data(&mut filtered_iters, opts).unwrap() /// } /// /// let fastq_path_str = "tests/test_files_in/14412_3#82.contigs_velvet.fa.gz"; @@ -490,7 +490,7 @@ pub fn sketch_data, Option>)>>( opts: SketchingOpts, #[cfg(feature = "3di")] convert_pdb: bool, #[cfg(feature = "3di")] struct_string: Option, -) -> Vec { +) -> anyhow::Result> { // Read in sequence and set up rolling hash by alphabet type let mut hash_its: Vec> = match opts.seq_type { HashType::DNA => NtHashIterator::new( @@ -499,7 +499,7 @@ pub fn sketch_data, Option>)>>( opts.add_rc, opts.min_qual, opts.is_reads, - ) + )? .into_iter() .map(|it| Box::new(it) as Box) .collect(), @@ -533,29 +533,28 @@ pub fn sketch_data, Option>)>>( } }; - hash_its - .iter_mut() - .enumerate() - .map(|(idx, hash_it)| { - let sample_name = if opts.concat_fasta { - format!("{}_{}", opts.name, idx + 1) - } else { - opts.name.to_string() - }; - if hash_it.seq_len() == 0 { - panic!("{sample_name} has no valid sequence"); - } - // Run the sketching - Sketch::new( - &mut **hash_it, - &sample_name, - &opts.k_vals, - opts.sketch_size, - opts.add_rc, - opts.min_count, - ) - }) - .collect::>() + let mut sketches: Vec = Vec::with_capacity(hash_its.len()); + for(idx, hash_it) in hash_its.iter_mut().enumerate() { + let sample_name = if opts.concat_fasta { + format!("{}_{}", opts.name, idx + 1) + } else { + opts.name.to_string() + }; + if hash_it.seq_len() == 0 { + return Err(anyhow::anyhow!("{sample_name} has no valid sequence")) + } + // Run the sketching + let sketch = Sketch::new( + &mut **hash_it, + &sample_name, + &opts.k_vals, + opts.sketch_size, + opts.add_rc, + opts.min_count, + )?; + sketches.push(sketch); + } + Ok(sketches) } #[cfg(not(target_arch = "wasm32"))] @@ -656,7 +655,7 @@ pub fn sketch_files( convert_pdb, #[cfg(feature = "3di")] di, - ) + ).unwrap() }) .for_each_with(tx, |tx, sketch| { // Emit the sketch results to the writer thread