Ensuring preseq mode and orig-dups are exclusive; improving verbose output

cmake/tests.cmake

@@ -33,3 +33,4 @@ add_test(NAME "Kmers" COMMAND bash test_scripts/kmers.sh)
 add_test(NAME "Groups" COMMAND bash test_scripts/groups.sh)
 add_test(NAME "SAM input" COMMAND bash test_scripts/sam.sh)
 add_test(NAME "Orig dups" COMMAND bash test_scripts/orig_dups.sh)
+add_test(NAME "Preseq output" COMMAND bash test_scripts/preseq.sh)

data/test_data/md5sum.txt

@@ -12,3 +12,4 @@ c2d51938c568b3c224efe9134baddb03  fastq_bgzf_threads_out/fastq_bgzip_1/fastqc_da
 be637ed4a4ebf2f1a64edf30ef1f8bb8  fastq_plain_out/fastq_1/fastqc_data.txt
 bc8fd6e40a0ca55cb00634bc306a896c  groups_out/bam_1/fastqc_data.txt
 a1952ae366bd7c7207f40db833b1d16b  sam_out/sam_1/fastqc_data.txt
+3fe18e2ee85e3912ede8c57559bf9f88  preseq_out/bam_1/preseq_hist.txt

data/test_scripts/preseq.sh

@@ -0,0 +1,19 @@
+#!/usr/bin/env bash
+# SPDX-License-Identifier: MIT
+
+prog=./falco
+infile=test_data/bam_1.bam
+outdir=preseq_out
+if [[ -e "${infile}" ]]; then
+    ${prog} --preseq -o ${outdir} ${infile}
+    x=$(md5sum --ignore-missing -c test_data/md5sum.txt | \
+            grep "${outdir}" | \
+            grep -c "OK$")
+    if [[ "${x}" != "1" ]]; then
+        exit 1;
+    fi
+    rm -r ${outdir}
+else
+    echo "${infile} not found; skipping remaining tests";
+    exit 77;
+fi

src/falco.cpp

@@ -7,8 +7,12 @@ static constexpr auto description =
 
 EXAMPLES:
 
-Use all defaults making output files in 'results/SRX081761_1':
+Use all defaults making output files in a dir named 'results/SRX081761_1':
 $ falco -o results SRX081761_1.fastq
+Output files will be:
+results/SRX081761_1/summary.txt
+results/SRX081761_1/fastqc_data.txt
+results/SRX081761_1/fastqc_results.html
 
 Use 8 cores and analyze all files in 'project' with the fq suffix, and results
 for each input file in their own subdirectory of 'results':
@@ -23,6 +27,11 @@ Skip adapter content and tile-specific analysis:
 Use configuration settings from 'my_limits.txt':
 $ falco --config my_limits.txt -o results project_sample2.fq.gz
 
+Generate a file with duplication info for preseq analysis:
+$ falco --preseq -o results SRX081761_1.fq.gz
+Output files include:
+results/SRX081761_1/preseq_hist.txt
+
 Default configuration files can be found here:
 {}
 Use these as templates. Copy and modify them to customize your analysis.
@@ -209,6 +218,9 @@ get_min_buffer_size(const auto max_read_length) {
 int
 main(int argc, char *argv[]) {
   try {
+    const auto start_time =
+      get_program_start_time();  // ADS: do it now or it's wrong forever
+
     static constexpr auto buffer_size_default = 256 * 1024 * 1024;
     static constexpr std::int64_t min_buf_size = 1024 * 1024;
     std::vector<std::string> infiles;
@@ -273,7 +285,7 @@ main(int argc, char *argv[]) {
       ->required()
       ->option_text(" ")
       ->check(CLI::ExistingFile);
-    app.add_option("-o,--output", outdir, "Output directory")
+    app.add_option("-o,--output", outdir, "Output directory (required)")
       ->required()
       ->option_text("DIR");
     app.add_option("-t,--threads", n_threads,
@@ -308,22 +320,24 @@ 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")
+    const auto preseq_opt =
+      app.add_flag("--preseq", do_preseq,
+                   "Make file for preseq input (excludes --orig-dups)")
       ->option_text(" ");
     const auto orig_dups_opt =
-      app.add_flag("--orig-dups", [&](const auto x) {
-        do_original_dups = x;
+      app.add_flag("--orig-dups", [&](const auto should_be_one) {
+        do_original_dups = should_be_one;
         do_dup_analysis = 1;
-      }, "Use original duplication mode (turns dups on)")
+      }, "Use original duplication mode (enables --dups)")
+      ->excludes(preseq_opt)
       ->option_text(" ");
     app.add_flag("--groups", do_groups, "Group base positions in output")
       ->option_text(" ");
     app.add_flag("--tiles,!--no-tiles", do_tiles,
                  "Toggle per-tile quality analysis (default: on)")
       ->option_text(" ");
     app.add_flag("--dups,!--no-dups", do_dup_analysis,
-                 "Toggle duplication/overrep analysis (default: on)")
+                 "Toggle sequence duplication analysis (default: on)")
       ->excludes(orig_dups_opt)
       ->option_text(" ");
     app.add_flag("--adap,!--no-adap", do_adap,
@@ -334,8 +348,6 @@ main(int argc, char *argv[]) {
       ->option_text(" ");
     // clang-format on
 
-    const auto start_time{std::chrono::high_resolution_clock::now()};
-
     if (argc < 2) {
       std::println("{}", app.help());
       return EXIT_SUCCESS;
@@ -403,38 +415,34 @@ main(int argc, char *argv[]) {
       max_read_length = get_max_read_length(buffer_size);
 
     if (verbose) {
-      std::println("threads requested: {}\n"
+      std::println("Falco v{}\n"
+                   "Run started: {}\n",
+                   VERSION, format_program_start_date_and_time());
+      std::println("Resources\n"
+                   "threads requested: {}\n"
                    "input memory buffer size: {}\n"
-                   "max analyzable read length: {}\n"
-                   "tile analysis requested: {}\n"
-                   "k-mer analysis requested: {}\n"
-                   "dups analysis requested: {}\n"
-                   "adapter analysis requested: {}\n"
-                   "use base groups in output: {}\n"
-                   "input files:",
+                   "max analyzable read length: {}\n",
                    n_threads, size_to_units(buffer_size),
-                   size_to_units(max_read_length, "bp"), mode.do_tiles(),
-                   mode.do_kmers(), mode.do_dups(), mode.do_adap(),
-                   mode.do_groups());
+                   size_to_units(max_read_length, "bp"));
+      std::println("Analyses\n{}", mode.string_verbose());
+      std::println("Input files");
       std::ranges::for_each(infos, [](const auto &info) {
         std::println("{}\t{}\t{}", info.name, info.description,
                      size_to_units(info.size, std::string{}));
       });
+      std::println();
     }
 
-    std::vector<dups_init_t> dups =
-      do_original_dups
-        ? initialize_original_duplicates(infiles, infos, n_threads)
-        : std::vector<dups_init_t>{};
+    auto dups = do_original_dups
+                  ? initialize_original_duplicates(infiles, infos, n_threads)
+                  : std::vector<dups_init_t>{};
     auto results =
       analyze(n_threads, mode, infos,
               make_reads_files(infos, infiles, buffer_size), std::move(dups));
     write_output(mode, infos, outdirs, std::move(results));
 
     if (verbose)
-      std::println(
-        "total run time: {:.6g}s",
-        duration(start_time, std::chrono::high_resolution_clock::now()));
+      std::println("Run duration: {}", get_run_duration(start_time));
   }
   catch (const std::exception &e) {
     std::println("{}", e.what());

src/falco_utils.cpp

@@ -5,6 +5,7 @@
 #include <array>
 #include <cmath>
 #include <cstdint>
+#include <ctime>  // for std::localtime
 #include <format>
 #include <span>
 #include <string>
@@ -142,3 +143,20 @@ combine_gc_content_for_lengths(const std::vector<falco::gc_content_array> &gcs)
   }
   return hist;
 }
+
+[[nodiscard]] auto
+get_program_start_time()
+  -> std::chrono::time_point<std::chrono::high_resolution_clock,
+                             std::chrono::nanoseconds> {
+  static const auto start_time = std::chrono::high_resolution_clock::now();
+  return start_time;
+}
+
+[[nodiscard]] auto
+format_program_start_date_and_time() -> std::string {
+  const auto t = get_program_start_time();
+  const auto t_c = std::chrono::system_clock::to_time_t(t);
+  std::ostringstream oss;
+  oss << std::put_time(std::localtime(&t_c), "%F %T %Z");
+  return oss.str();
+}

src/falco_utils.hpp

@@ -90,12 +90,18 @@ get_theoretical_distribution(const std::vector<double> &gc,
 [[nodiscard]] auto
 sum_deviation_from_normal(const std::vector<double> &gc) -> double;
 
-[[nodiscard]] inline constexpr auto
-duration(const auto start, const auto stop) {
-  const auto d = stop - start;
-  // ADS: 'count()' because macos has locale issues formatting times
-  return std::chrono::duration_cast<std::chrono::duration<double>>(d).count();
-};
+[[nodiscard]] auto
+get_run_duration(const auto start_time) {
+  using namespace std::literals::chrono_literals;
+  const auto d = std::chrono::high_resolution_clock::now() - start_time;
+  const auto d_ms = std::chrono::floor<std::chrono::milliseconds>(d);
+  if (d < 1min) {
+    const auto ds =
+      std::chrono::duration_cast<std::chrono::duration<double>>(d).count();
+    return std::format("0:{:05.2f}", ds);
+  }
+  return std::format("{}", std::chrono::hh_mm_ss{d_ms});
+}
 
 inline constexpr auto end_module_tag = ">>END_MODULE\n";
 
@@ -324,4 +330,12 @@ estimate_read_length_fastq_chunk(const auto &data, const auto n) {
   return total / (std::size(lines) / fastq_lines_per_read);
 }
 
+[[nodiscard]] auto
+get_program_start_time()
+  -> std::chrono::time_point<std::chrono::high_resolution_clock,
+                             std::chrono::nanoseconds>;
+
+[[nodiscard]] auto
+format_program_start_date_and_time() -> std::string;
+
 #endif  // SRC_FALCO_UTILS_HPP_

src/html.cpp

@@ -69,16 +69,12 @@ get_html_module(const std::string &label, const std::string &text,
 [[nodiscard]] auto
 falco_get_html(const file_info &info, const file_grades &grades,
                const std::string &analysis_modules) -> std::string {
-  const auto now = std::chrono::system_clock::now();
-  const auto now_c = std::chrono::system_clock::to_time_t(now);
-  std::ostringstream oss;
-  oss << std::put_time(std::localtime(&now_c), "%F %T %Z");
-  return fmt::format(falco_html_body,                           //
-                     fmt::arg("date", oss.str()),               //
-                     fmt::arg("filename", info.name),           //
-                     fmt::arg("style", style),                  //
-                     fmt::arg("summary", get_summary(grades)),  //
-                     fmt::arg("modules", analysis_modules),     //
+  return fmt::format(falco_html_body,                                         //
+                     fmt::arg("date", format_program_start_date_and_time()),  //
+                     fmt::arg("filename", info.name),                         //
+                     fmt::arg("style", style),                                //
+                     fmt::arg("summary", get_summary(grades)),                //
+                     fmt::arg("modules", analysis_modules),                   //
                      fmt::arg("version", VERSION));
 }
 

src/run_mode.cpp

@@ -2,7 +2,10 @@
 
 #include "run_mode.hpp"
 
+#include <algorithm>
+#include <format>
 #include <iterator>
+#include <ranges>
 #include <string>
 #include <unordered_map>
 #include <vector>
@@ -29,6 +32,16 @@
 // sequence_length  | do_length     |
 // tile             | do_tiles      | Yes
 // ------------------------------------------------------------
+//
+// ADS (2026-08-15 edit): additional modes have been added, but I'm not going to
+// adjust the above right now because of time. So here are the extras, none of
+// which can be set in the config file, only cli:
+//
+// do_groups (had been part of run_mode)
+// do_bisulfite (changes grading and nothing else)
+// do_original_dups (use the original duplication mode from FastQC and Falco v1)
+// do_preseq (make another output file for input to preseq)
+//
 // clang-format on
 
 // clang-format off
@@ -70,6 +83,12 @@ run_mode::assign(const std::unordered_map<std::string, bool> &modes) -> void {
 auto
 run_mode::set_unassigned() -> void {
   // clang-format off
+  // 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;
+  //
   if (do_adap_ == 0) do_adap_ = do_adap_default;
   if (do_dups_ == 0) do_dups_ = do_dups_default;
   if (do_gc_content_ == 0) do_gc_content_ = do_gc_content_default;
@@ -81,10 +100,50 @@ run_mode::set_unassigned() -> void {
   if (do_sequence_ == 0) do_sequence_ = do_sequence_default;
   if (do_length_ == 0) do_length_ = do_length_default;
   if (do_tiles_ == 0) do_tiles_ = do_tiles_default;
-  // 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
 }
+
+[[nodiscard]] auto
+run_mode::string_verbose() const -> std::string {
+  static constexpr auto spacer = 2;
+  // ADS: need a mechanism to keep this map it in sync with var names and labels
+  static const auto descriptions = std::unordered_map<std::string, std::string>{
+    {"adapter", "adapter content"},
+    {"duplication", "sequence duplication"},
+    {"gc_sequence", "GC content"},
+    {"kmer", "k-mer content"},
+    {"n_content", "N content"},
+    {"overrepresented", "overrep sequences"},
+    {"quality_base", "per-base quality"},
+    {"quality_sequence", "per-sequence quality"},
+    {"sequence", "sequence composition"},
+    {"sequence_length", "sequence length"},
+    {"tile", "per-tile quality"},
+  };
+  static const auto fmt_mode = [&](const std::string &label,
+                                   const auto the_mode) {
+    const auto itr = descriptions.find(label);
+    if (itr == std::cend(descriptions))
+      throw std::runtime_error("failed to find mode: " + label);
+    return std::pair{itr->second, std::string{the_mode ? "ON" : "OFF"}};
+  };
+  std::vector<std::pair<std::string, std::string>> lines;
+  lines.push_back(fmt_mode("adapter", do_adap()));
+  lines.push_back(fmt_mode("duplication", do_dups()));
+  lines.push_back(fmt_mode("gc_sequence", do_gc_content()));
+  lines.push_back(fmt_mode("kmer", do_kmers()));
+  lines.push_back(fmt_mode("sequence_length", do_length()));
+  lines.push_back(fmt_mode("n_content", do_n_content()));
+  lines.push_back(fmt_mode("overrepresented", do_overrep()));
+  lines.push_back(fmt_mode("quality_base", do_qual_base()));
+  lines.push_back(fmt_mode("quality_sequence", do_qual_seq()));
+  lines.push_back(fmt_mode("sequence", do_sequence()));
+  lines.push_back(fmt_mode("tile", do_tiles()));
+  const auto w = spacer + std::ranges::max(std::views::transform(
+                            std::views::values(descriptions),
+                            [](const auto &x) { return std::size(x); }));
+  std::string r;
+  for (const auto &[descr, val] : lines)
+    r += std::format("{:{}} {}\n", std::format("{}:", descr), w, val);
+  return r;
+}

src/run_mode.hpp

@@ -13,9 +13,7 @@ class run_mode {
   set_unassigned() -> void;
 
   [[nodiscard]] auto
-  string() const -> std::string {
-    return {};
-  }
+  string_verbose() const -> std::string;
 
   // clang-format off
   // ADS: these first params are not set in config file