Adding an option to output duplication info in a format that preseq can read directly

src/duplication_results.cpp

@@ -89,13 +89,25 @@ duplication_results::get_n_counted_reads() const -> std::uint64_t {
   return std::reduce(std::cbegin(dups_v), std::cend(dups_v));
 }
 
+[[nodiscard]] auto
+duplication_results::get_preseq_hist() const
+  -> std::vector<std::pair<std::uint32_t, std::uint32_t>> {
+  std::unordered_map<std::uint32_t, std::uint32_t> counter;
+  std::ranges::for_each(std::views::values(dups),
+                        [&](const auto v) { ++counter[v]; });
+  std::vector<std::pair<std::uint32_t, std::uint32_t>> counts;
+  std::ranges::copy(counter, std::back_inserter(counts));
+  std::ranges::sort(counts);
+  return counts;
+}
+
 [[nodiscard]] auto
 duplication_results::get_overrepresented(const std::uint64_t n_reads) const
   -> std::vector<overrep_t> {
   const auto cutoff = static_cast<double>(n_reads) * overrep_cutoff;
-  const auto gt_cutoff = [&](const auto p) { return p.second >= cutoff; };
+  const auto gte_cutoff = [&](const auto p) { return p.second >= cutoff; };
   const auto rev_p = [&](const auto p) { return std::pair{p.second, p.first}; };
-  auto overrep = dups | std::views::filter(gt_cutoff) |
+  auto overrep = dups | std::views::filter(gte_cutoff) |
                  std::views::transform(rev_p) | std::ranges::to<std::vector>();
   std::ranges::sort(overrep, std::greater{});
   std::vector<overrep_t> ret;

src/duplication_results.hpp

@@ -75,6 +75,10 @@ struct duplication_results {
     dups_map_t().swap(dups);
   }
 
+  [[nodiscard]] auto
+  get_preseq_hist() const
+    -> std::vector<std::pair<std::uint32_t, std::uint32_t>>;
+
   [[nodiscard]] auto
   get_dups_summary(const std::uint64_t n_reads) const -> dup_summary_t;
 

src/falco.cpp

@@ -92,13 +92,16 @@ write_output(
   static constexpr auto report_filename = "fastqc_data.txt";
   static constexpr auto html_filename = "fastqc_report.html";
   static constexpr auto summary_filename = "summary.txt";
+  static constexpr auto preseq_filename = "preseq_hist.txt";
   for (const auto [result, info, outdir] :
        std::views::zip(results, infos, outdirs)) {
     const auto outdir_path = std::filesystem::path{outdir};
     const auto summary = results_summary(std::move(result), mode, info);
     write_file(outdir_path / report_filename, summary.get_report());
     write_file(outdir_path / html_filename, summary.get_html());
     write_file(outdir_path / summary_filename, summary.get_summary());
+    if (mode.do_preseq())
+      write_file(outdir_path / preseq_filename, summary.get_preseq_hist());
   }
 }
 
@@ -222,6 +225,7 @@ main(int argc, char *argv[]) {
 
     int do_groups{};
     int do_bisulfite{};
+    int do_preseq{};
     int do_original_dups{};
 
     std::uint32_t n_threads{1};
@@ -304,6 +308,9 @@ main(int argc, char *argv[]) {
     app.add_flag("--bisulfite", do_bisulfite,
                  "Assume bisulfite when grading sequence content")
       ->option_text(" ");
+    app.add_flag("--preseq", do_preseq,
+                 "Write duplication info file for analysis by preseq")
+      ->option_text(" ");
     const auto orig_dups_opt =
       app.add_flag("--orig-dups", [&](const auto x) {
         do_original_dups = x;
@@ -352,6 +359,7 @@ main(int argc, char *argv[]) {
     mode.set_do_tiles(do_tiles);
     mode.set_do_groups(do_groups);
     mode.set_do_bisulfite(do_bisulfite);
+    mode.set_do_preseq(do_preseq);
     mode.set_do_original_dups(do_original_dups);
     mode.set_unassigned();
 

src/results_summary.cpp

@@ -235,3 +235,12 @@ results_summary::get_html() const -> std::string {
 results_summary::get_summary() const -> std::string {
   return grades.to_string(info.name);
 }
+
+[[nodiscard]] auto
+results_summary::get_preseq_hist() const -> std::string {
+  const auto hist = dr.get_preseq_hist();
+  std::string r;
+  for (const auto &[n_times_seen, n_reads] : hist)
+    r += std::format("{}\t{}\n", n_times_seen, n_reads);
+  return r;
+}

src/results_summary.hpp

@@ -84,6 +84,9 @@ struct results_summary {
 
   [[nodiscard]] auto
   get_summary() const -> std::string;
+
+  [[nodiscard]] auto
+  get_preseq_hist() const -> std::string;
 };
 
 #endif  // SRC_RESULTS_SUMMARY_HPP_

src/run_mode.cpp

@@ -84,6 +84,7 @@ run_mode::set_unassigned() -> void {
   // settings below are not in config file
   if (do_groups_ == 0) do_groups_ = do_groups_default;
   if (do_bisulfite_ == 0) do_bisulfite_ = do_bisulfite_default;
+  if (do_preseq_ == 0) do_preseq_ = do_preseq_default;
   if (do_original_dups_ == 0) do_original_dups_ = do_original_dups_default;
   // clang-format on
 }

src/run_mode.hpp

@@ -21,6 +21,7 @@ class run_mode {
   // ADS: these first params are not set in config file
   [[nodiscard]] auto do_groups() const -> bool { return do_groups_ == 1; }
   [[nodiscard]] auto do_bisulfite() const -> bool { return do_bisulfite_ == 1; }
+  [[nodiscard]] auto do_preseq() const -> bool { return do_preseq_ == 1; }
   [[nodiscard]] auto do_original_dups() const -> bool { return do_original_dups_ == 1; }
   //
   [[nodiscard]] auto do_adap() const -> bool { return do_adap_ == 1; }
@@ -43,6 +44,7 @@ class run_mode {
   // clang-format off
   auto set_do_groups(const int x) { if (x) do_groups_ = x; }
   auto set_do_bisulfite(const int x) { if (x) do_bisulfite_ = x; }
+  auto set_do_preseq(const int x) { if (x) do_preseq_ = x; }
   auto set_do_original_dups(const int x) { if (x) do_original_dups_ = x; }
   //
   auto set_do_adap(const int x) { if (x) do_adap_ = x; }
@@ -69,6 +71,7 @@ class run_mode {
   // first settings are not in config file
   static constexpr auto do_groups_default = -1;         // OFF
   static constexpr auto do_bisulfite_default = -1;      // OFF
+  static constexpr auto do_preseq_default = -1;         // OFF
   static constexpr auto do_original_dups_default = -1;  // OFF
 
   static constexpr auto do_adap_default = 1;  // affects processing
@@ -86,6 +89,7 @@ class run_mode {
   // not in config file
   int do_groups_{};
   int do_bisulfite_{};
+  int do_preseq_{};
   int do_original_dups_{};
 
   int do_adap_{};