Covbuildmerge

src/cli.rs

@@ -84,6 +84,31 @@
     }
 }
 
+#[doc(hidden)]
+pub fn valid_min_kmer(s: &str) -> Result<ValidMinKmer, String> {
+    match s {
+        s if s.eq(&String::from("auto")) => Ok(ValidMinKmer::Auto),
+        s => {
+            // Throw error if it cannot be parsed to u16
+            let x: u16 = s.parse().expect("Invalid minimum kmer count");
+            if x.ge(&1) {
+                log::info!("Using provided minimum kmer count of {}", x);
+                Ok(ValidMinKmer::Val(x))
+            } else {
+                Err("Minimum kmer count must be >= 1".to_string())
+            }
+        }
+    }
+}
+/// Possible user input for minimum kmer count threshold
+#[derive(Copy, Clone, PartialEq, Eq, PartialOrd, Ord)]
+pub enum ValidMinKmer {
+    /// Attempt to calculate using cov command
+    Auto,
+    /// User provided u16 value for threshold
+    Val(u16),
+}
+
 /// Possible output file types
 #[derive(Copy, Clone, PartialEq, Eq, PartialOrd, Ord, ValueEnum)]
 pub enum FileType {
@@ -167,8 +192,8 @@
         single_strand: bool,
 
         /// Minimum k-mer count (with reads)
-        #[arg(long, default_value_t = DEFAULT_MINCOUNT)]
-        min_count: u16,
+        #[arg(long, value_parser = valid_min_kmer)]
+        min_count: Option<ValidMinKmer>,
 
         /// Minimum k-mer quality (with reads)
         #[arg(long, default_value_t = DEFAULT_MINQUAL)]

src/coverage.rs

@@ -392,7 +392,7 @@ mod tests {
         };
 
         let cutoff = test_obj.fit_histogram();
-        assert_eq!(cutoff.is_ok(), true);
+        assert!(cutoff.is_ok());
         assert_eq!(cutoff.unwrap(), 9);
 
         // The other template

src/io_utils.rs

@@ -14,11 +14,13 @@
 use crate::merge_ska_array::MergeSkaArray;
 use crate::merge_ska_dict::{build_and_merge, InputFastx};
 use crate::ska_dict::bit_encoding::UInt;
+use crate::CoverageHistogram;
 
 use crate::cli::{
     DEFAULT_KMER, DEFAULT_MINCOUNT, DEFAULT_MINQUAL, DEFAULT_PROPORTION_READS, DEFAULT_QUALFILTER,
     DEFAULT_STRAND,
 };
+use crate::ValidMinKmer;
 
 /// Given a list of input files, parses them into triples of name, filename and
 /// [`None`] to be used with [SkaDict](`crate::ska_dict::SkaDict::new()`).
@@ -140,12 +142,69 @@
 
 /// Checks if any input files are fastq
 pub fn any_fastq(files: &[InputFastx]) -> bool {
-    let mut fastq = false;
-    for file in files {
-        if file.2.is_some() {
-            fastq = true;
-            break;
+    files.iter().any(|file| file.2.is_some())
+}
+
+/// Counts number of fastq files
+pub fn count_fastq(files: &[InputFastx]) -> usize {
+    files.iter().filter(|file| file.2.is_some()).count()
+}
+
+/// Collects filepaths to the first 2 fastq files into a tuple
+pub fn get_2_fastq_path(files: &[InputFastx]) -> (String, String) {
+    let out: Vec<String> = files
+        .iter()
+        .filter(|file| file.2.is_some())
+        .take(2)
+        .map(|x| x.1.clone())
+        .collect();
+
+    if out.len().gt(&1) {
+        (out[0].clone(), out[1].clone())
+    } else {
+        panic!("Trying to get 2 fastq files from a vector with <2 elements");
+    }
+}
+
+/// Calculates minimum kmer cutoff depending on user provided argument
+pub fn kmer_min_cutoff<IntT: for<'a> UInt<'a>>(
+    v: &Option<ValidMinKmer>,
+    files: &[InputFastx],
+    k: &usize,
+    rc: bool,
+    verbose: bool,
+) -> u16 {
+    // Minimum kmer cutoff logic
+    if v.is_none() {
+        log::info!(
+            "Using user-provided minimum kmer value of {}",
+            DEFAULT_MINCOUNT
+        );
+        DEFAULT_MINCOUNT
+    } else {
+        match v.unwrap() {
+            // User-provided value (already checked by cli parser)
+            ValidMinKmer::Val(x) => {
+                log::info!("Using user-provided minimum kmer value of {}", x);
+                x
+            }
+            // auto-calculate & there are enough fastq files
+            ValidMinKmer::Auto if count_fastq(files).ge(&2) => {
+                let (fastq_fwd, fastq_rev) = get_2_fastq_path(files);
+                let mut cov =
+                    CoverageHistogram::<IntT>::new(&fastq_fwd, &fastq_rev, *k, rc, verbose);
+                let out = cov.fit_histogram().expect("Couldn't fit coverage model") as u16;
+                cov.plot_hist();
+                log::info!("Using inferred minimum kmer value of {}", out);
+                out
+            }
+            // Not enough fastq files, use default and warn user
+            ValidMinKmer::Auto => {
+                log::info!(
+                    "Not enough fastq files to fit mixture model, using default kmer count of 5"
+                );
+                DEFAULT_MINCOUNT
+            }
         }
     }
-    fastq
 }

src/lib.rs

@@ -503,17 +503,18 @@ pub fn main() {
 
             // Read input
             let input_files = get_input_list(file_list, seq_files);
-            let quality = QualOpts {
-                min_count: *min_count,
+            let rc = !*single_strand;
+            // Build, merge
+            // check for 64 or 128 bit representation
+            let mut quality = QualOpts {
+                min_count: 0,
                 min_qual: *min_qual,
                 qual_filter: *qual_filter,
             };
-
-            // Build, merge
-            let rc = !*single_strand;
-
-            if *k <= 31 {
+            if k.le(&31) {
                 log::info!("k={}: using 64-bit representation", *k);
+                quality.min_count =
+                    kmer_min_cutoff::<u64>(min_count, &input_files, k, rc, args.verbose);
                 let merged_dict = build_and_merge::<u64>(
                     &input_files,
                     *k,
@@ -527,6 +528,8 @@ pub fn main() {
                 save_skf(&merged_dict, output);
             } else {
                 log::info!("k={}: using 128-bit representation", *k);
+                quality.min_count =
+                    kmer_min_cutoff::<u128>(min_count, &input_files, k, rc, args.verbose);
                 let merged_dict = build_and_merge::<u128>(
                     &input_files,
                     *k,

src/merge_ska_array.rs

@@ -625,7 +625,7 @@ pub struct KmerIter<'a, IntT> {
     index: usize,
 }
 
-impl<'a, IntT> Iterator for KmerIter<'a, IntT>
+impl<IntT> Iterator for KmerIter<'_, IntT>
 where
     IntT: for<'b> UInt<'b>,
 {

src/ska_dict/bit_encoding.rs

@@ -124,16 +124,16 @@
     fn hash_val(self, hash_fn: &RandomState) -> u64;
 }
 
-impl<'a> UInt<'a> for u64 {
+impl UInt<'_> for u64 {
     #[inline(always)]
     fn rev_comp(mut self, k_size: usize) -> Self {
         // This part reverses the bases by shuffling them using an on/off pattern
         // of bits
         self = (self >> 2 & 0x3333333333333333) | (self & 0x3333333333333333) << 2;
         self = (self >> 4 & 0x0F0F0F0F0F0F0F0F) | (self & 0x0F0F0F0F0F0F0F0F) << 4;
         self = (self >> 8 & 0x00FF00FF00FF00FF) | (self & 0x00FF00FF00FF00FF) << 8;
         self = (self >> 16 & 0x0000FFFF0000FFFF) | (self & 0x0000FFFF0000FFFF) << 16;
         self = (self >> 32 & 0x00000000FFFFFFFF) | (self & 0x00000000FFFFFFFF) << 32;
         // This reverse complements
         self ^= 0xAAAAAAAAAAAAAAAA;
 
@@ -182,22 +182,22 @@
     }
 }
 
-impl<'a> UInt<'a> for u128 {
+impl UInt<'_> for u128 {
     #[inline(always)]
     fn rev_comp(mut self, k_size: usize) -> u128 {
         // This part reverses the bases by shuffling them using an on/off pattern
         // of bits
         self = (self >> 2 & 0x33333333333333333333333333333333)
             | (self & 0x33333333333333333333333333333333) << 2;
         self = (self >> 4 & 0x0F0F0F0F0F0F0F0F0F0F0F0F0F0F0F0F)
             | (self & 0x0F0F0F0F0F0F0F0F0F0F0F0F0F0F0F0F) << 4;
         self = (self >> 8 & 0x00FF00FF00FF00FF00FF00FF00FF00FF)
             | (self & 0x00FF00FF00FF00FF00FF00FF00FF00FF) << 8;
         self = (self >> 16 & 0x0000FFFF0000FFFF0000FFFF0000FFFF)
             | (self & 0x0000FFFF0000FFFF0000FFFF0000FFFF) << 16;
         self = (self >> 32 & 0x00000000FFFFFFFF00000000FFFFFFFF)
             | (self & 0x00000000FFFFFFFF00000000FFFFFFFF) << 32;
         self = (self >> 64 & 0x0000000000000000FFFFFFFFFFFFFFFF)
             | (self & 0x0000000000000000FFFFFFFFFFFFFFFF) << 64;
         // This reverse complements
         self ^= 0xAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA;

src/ska_ref.rs

@@ -99,7 +99,6 @@ where
     ambig_mask: bool,
 
     /// Input sequence
-
     /// Chromosome names
     chrom_names: Vec<String>,
     /// Sequence, indexed by chromosome, then position
@@ -108,7 +107,6 @@ where
     repeat_coors: Vec<usize>,
 
     /// Mapping information
-
     /// Positions of mapped bases as (chrom, pos)
     mapped_pos: Vec<(usize, usize)>,
     /// Array of mapped bases, rows loci, columns samples

src/ska_ref/idx_check.rs

@@ -47,7 +47,7 @@ pub struct IdxCheckIter<'a> {
     idx: usize,
 }
 
-impl<'a> Iterator for IdxCheckIter<'a> {
+impl Iterator for IdxCheckIter<'_> {
     type Item = (usize, usize);
 
     fn next(&mut self) -> Option<(usize, usize)> {

tests/align.rs

@@ -14,7 +14,7 @@ use pretty_assertions::assert_eq;
 fn build_cli() {
     let sandbox = TestSetup::setup();
     // Create an rfile in the tmp dir
-    let rfile_name = sandbox.create_rfile(&"test", FxType::Fasta);
+    let rfile_name = sandbox.create_rfile("test", FxType::Fasta);
 
     Command::new(cargo_bin("ska"))
         .current_dir(sandbox.get_wd())
@@ -23,7 +23,7 @@ fn build_cli() {
         .arg(rfile_name)
         .arg("-o")
         .arg("basic_build_opts")
-        .args(&["-v", "--threads", "2", "-k", "31"])
+        .args(["-v", "--threads", "2", "-k", "31"])
         .assert()
         .success();
 
@@ -68,7 +68,7 @@ fn align_cli() {
         .arg(sandbox.file_string("test_2.fa", TestDir::Input))
         .arg("-o")
         .arg("basic.aln")
-        .args(&[
+        .args([
             "-v",
             "--threads",
             "2",
@@ -359,7 +359,7 @@ fn parallel_align() {
         .arg(rfile_name)
         .arg("-o")
         .arg("serial_build")
-        .args(&["-v", "--threads", "1", "-k", "15"])
+        .args(["-v", "--threads", "1", "-k", "15"])
         .assert()
         .success();
 
@@ -379,7 +379,7 @@ fn parallel_align() {
         .arg(rfile_name)
         .arg("-o")
         .arg("parallel_build")
-        .args(&["-v", "--threads", "4", "-k", "15"])
+        .args(["-v", "--threads", "4", "-k", "15"])
         .assert()
         .success();
 

tests/common/mod.rs

@@ -13,11 +13,11 @@ use hashbrown::HashSet;
 use pretty_assertions::assert_eq;
 
 // Creates correct path for input/output files
-static FILE_IN: &'static str = "tests/test_files_in";
-static FILE_TEST: &'static str = "tests/test_results_correct";
-static SYM_IN: &'static str = "input";
-static SYM_TEST: &'static str = "correct";
-static RFILE_NAME: &'static str = "file_list.txt";
+static FILE_IN: &str = "tests/test_files_in";
+static FILE_TEST: &str = "tests/test_results_correct";
+static SYM_IN: &str = "input";
+static SYM_TEST: &str = "correct";
+static RFILE_NAME: &str = "file_list.txt";
 
 #[derive(Debug, PartialEq, Copy, Clone)]
 pub enum TestDir {
@@ -185,7 +185,7 @@ pub fn var_hash(aln_string: &[u8]) -> HashSet<Vec<char>> {
         variant_pairs.insert(var_vec);
     }
 
-    return variant_pairs;
+    variant_pairs
 }
 
 // Helper for comparing mapped alignments with different sample names

tests/distance.rs

@@ -16,15 +16,15 @@ fn basic_dists() {
         .arg("distance")
         .arg(sandbox.file_string("merge.skf", TestDir::Input))
         .arg("-v")
-        .args(&["--threads", "2"])
+        .args(["--threads", "2"])
         .assert()
         .stdout_eq_path(sandbox.file_string("merge.dist.stdout", TestDir::Correct));
 
     Command::new(cargo_bin("ska"))
         .current_dir(sandbox.get_wd())
         .arg("distance")
         .arg(sandbox.file_string("merge.skf", TestDir::Input))
-        .args(&["-o", "dists.txt"])
+        .args(["-o", "dists.txt"])
         .assert()
         .success();
     sandbox.file_check("dists.txt", "merge.dist.stdout");
@@ -35,7 +35,7 @@ fn basic_dists() {
         .arg("distance")
         .arg(sandbox.file_string("merge_k41.skf", TestDir::Input))
         .arg("-v")
-        .args(&["--threads", "2"])
+        .args(["--threads", "2"])
         .assert()
         .stdout_eq_path(sandbox.file_string("merge_k41.dist.stdout", TestDir::Correct));
 }
@@ -52,7 +52,7 @@ fn dist_filter() {
         .arg(sandbox.file_string("merge_k9.skf", TestDir::Input))
         .arg("--allow-ambiguous")
         .arg("-v")
-        .args(&["--threads", "2"])
+        .args(["--threads", "2"])
         .assert()
         .stdout_eq_path(sandbox.file_string("merge_k9.dist.stdout", TestDir::Correct));
 
@@ -70,7 +70,7 @@ fn dist_filter() {
         .current_dir(sandbox.get_wd())
         .arg("distance")
         .arg(sandbox.file_string("merge_k9.skf", TestDir::Input))
-        .args(&["--min-freq", "1"])
+        .args(["--min-freq", "1"])
         .arg("-v")
         .assert()
         .stdout_eq_path(sandbox.file_string("merge_k9_min_freq.dist.stdout", TestDir::Correct));
@@ -103,9 +103,9 @@ fn multisample_dists() {
         .arg("distance")
         .arg("multidist.skf")
         .arg("-v")
-        .args(&["--min-freq", "0"])
+        .args(["--min-freq", "0"])
         .arg("--allow-ambiguous")
-        .args(&["--threads", "2"])
+        .args(["--threads", "2"])
         .assert()
         .stdout_eq_path(sandbox.file_string("multidist.stdout", TestDir::Correct));
 

tests/fasta_input.rs

@@ -185,7 +185,7 @@ fn repeats() {
         .arg("dup_ss.skf")
         .arg("--filter")
         .arg("no-const")
-        .args(&["--min-freq", "1"])
+        .args(["--min-freq", "1"])
         .assert()
         .success();
 
@@ -241,7 +241,7 @@ fn palindromes() {
         .current_dir(sandbox.get_wd())
         .arg("align")
         .arg("otto.skf")
-        .args(&["--filter", "no-filter"])
+        .args(["--filter", "no-filter"])
         .assert()
         .stdout_eq_path(sandbox.file_string("palindrome.stdout", TestDir::Correct));
 
@@ -285,7 +285,7 @@ fn palindromes() {
         .current_dir(sandbox.get_wd())
         .arg("align")
         .arg("ottootto.skf")
-        .args(&["--filter", "no-filter"])
+        .args(["--filter", "no-filter"])
         .assert()
         .stdout_eq_path(sandbox.file_string("palindrome_reps.stdout", TestDir::Correct));
 }