Add mode to automatically determine a count cutoff

Cargo.toml

@@ -1,10 +1,13 @@
 [package]
 name = "ska"
-version = "0.2.5"
+version = "0.3.0"
 authors = [
     "John Lees <[email protected]>",
     "Simon Harris <[email protected]>",
     "Johanna von Wachsmann <[email protected]>",
+    "Tommi Maklin <[email protected]>",
+    "Joel Hellewell <[email protected]>",
+    "Timothy Russell <[email protected]>"
 ]
 edition = "2021"
 description = "Split k-mer analysis"
@@ -44,6 +47,10 @@ hashbrown = "0.12"
 ahash = "0.8.2"
 ndarray = { version = "0.15.6", features = ["serde"] }
 num-traits = "0.2.15"
+# coverage model
+libm = "0.2.7"
+argmin = { version = "0.8.1", features = ["slog-logger"] }
+argmin-math = "0.3.0"
 
 [dev-dependencies]
 # testing

scripts/plot_cov.py

@@ -0,0 +1,79 @@
+#!/usr/bin/env python
+# vim: set fileencoding=<utf-8> :
+
+import matplotlib as mpl
+mpl.use('Agg')
+mpl.rcParams.update({'font.size': 18})
+import matplotlib.pyplot as plt
+from math import log
+
+import argparse
+
+def norm_series(series: list()):
+    norm = max(series)
+    series = [item/norm for item in series]
+    return series
+
+def main():
+    parser = argparse.ArgumentParser(
+                    prog='plot_cov',
+                    description='Plot the results of `ska cov`',
+                    epilog='Requires matplotlib')
+    parser.add_argument('histfile', help="Input table (stdout from `ska cov`)")
+    parser.add_argument('--output', help="Output prefix", default="coverage_histogram")
+    args = parser.parse_args()
+
+    cutoff = 0
+    kmers = list()
+    density = list()
+    idx_xseries = list()
+    with open(args.histfile, 'r') as hfile:
+        hfile.readline()
+        for line in hfile:
+            (idx, count, ll, comp) = line.rstrip().split("\t")
+            kmers.append(int(count))
+            density.append(float(ll))
+            idx_xseries.append(int(idx))
+            if comp == "Coverage" and cutoff == 0:
+                cutoff = int(idx)
+
+    k_norm = norm_series(kmers)
+    #d_norm = norm_series(density)
+
+    fig, (ax1, ax2) = plt.subplots(2)
+    fig.suptitle('Coverage histogram fit')
+    fig.set_dpi(160)
+    fig.set_facecolor('w')
+    fig.set_edgecolor('k')
+    fig.set_figwidth(11)
+    fig.set_figheight(11)
+    plt.tight_layout()
+
+    ax1.set_xlabel('K-mer count')
+    ax1.set_ylabel('Frequency')
+    ax1.set_ylim(0, k_norm[1])
+    ax1.plot(idx_xseries, k_norm, color='black', linewidth=2,
+             label='K-mer count frequency')
+    ax1.plot(idx_xseries, density, color='red', linewidth=2, linestyle='--',
+             label='Mixture model fit')
+    ax1.plot([cutoff, cutoff], [0, 1], color='darkgray', linewidth=1,
+                    linestyle='-.', label=f'Count cutoff ({cutoff})')
+    ax1.legend(loc='upper right')
+
+    ax2.set_yscale("log")
+    ax2.set_xlabel('K-mer count')
+    ax2.set_ylabel('log(Frequency)')
+    ax2.set_ylim(min(k_norm), k_norm[1])
+    ax2.plot(idx_xseries, k_norm, color='black', linewidth=2,
+             label='K-mer count frequency')
+    ax2.plot(idx_xseries, density, color='red', linewidth=2, linestyle='--',
+             label='Mixture model fit')
+    ax2.plot([cutoff, cutoff], [0, 1], color='darkgray', linewidth=1,
+                    linestyle='-.', label=f'Count cutoff ({cutoff})')
+
+    plt.savefig(args.output + ".png",
+                bbox_inches='tight')
+    plt.close()
+
+if __name__ == '__main__':
+    main()

src/cli.rs

@@ -249,6 +249,22 @@ pub enum Commands {
         #[arg(long, default_value_t = false)]
         full_info: bool,
     },
+    /// Estimate a coverage cutoff using a k-mer count profile (FASTQ only)
+    Cov {
+        /// FASTQ file (or .fastq.gz) with forward reads
+        fastq_fwd: String,
+
+        /// FASTQ file (or .fastq.gz) with reverse reads
+        fastq_rev: String,
+
+        /// K-mer size
+        #[arg(short, value_parser = valid_kmer, default_value_t = DEFAULT_KMER)]
+        k: usize,
+
+        /// Ignore reverse complement (all reads are oriented along same strand)
+        #[arg(long, default_value_t = DEFAULT_STRAND)]
+        single_strand: bool,
+    },
 }
 
 /// Function to parse command line args into [`Args`] struct