Original duplicates now built-in without needing special build args

CMakeLists.txt

@@ -199,13 +199,17 @@ add_library(run_mode OBJECT src/run_mode.cpp)
 add_library(falco_config OBJECT src/falco_config.cpp)
 target_link_libraries(falco_config PUBLIC run_mode)
 
+add_library(original_duplicates OBJECT src/original_duplicates.cpp)
+target_link_libraries(original_duplicates PUBLIC HTSLIB::HTSLIB)
+
 add_executable(falco src/falco.cpp)
 target_link_libraries(falco PUBLIC
   tile_processor
   contaminants
   adapter_matcher
   adapter_set
   duplication_results
+  original_duplicates
   kmer_counter
   fastq_file
   fastq_gz_file

README.md

@@ -15,14 +15,9 @@ check large sequencing reads for common problems. Falco was rewritten for
 version 2.0 in order to facilitate incoroprating new functionality moving
 forward.
 
-Output files and formats are the same: `summary.txt`, `fastqc_report.txt` and
-`fastqc_report.html`
-
 ## Quick start
 
-You will be able to find binaries for Linux and macOS with the releases.  As of
-2026-08-03, the falco v2.0 source has been merged into the main (master) branch
-but the "release" might lag by a day or so.
+You will be able to find binaries for Linux and macOS with the releases.
 
 Example:
 ```
@@ -123,13 +118,6 @@ GitHub's macOS look
 
 ## Intended changes in Falco v2.0
 
-Whether or not changes are "correct," they can be bad for users. If you have
-many years of experience interpreting the output of FastQC or falco, that alone
-is enough to give value to those results. Issues of correctness or accuracy
-might not matter to many users. Below are intended changes. Anything else is
-likely a bug and I'm happy to fix it. I'm also happy to reconsider these
-intended changes.
-
 ### Tiles results
 
 I found that the method for tile analysis is a bit unstable, and the tile grade
@@ -163,11 +151,8 @@ the same, the numbers differ dramatically. The motivation for the change is to
 produce more useful output, and to soon build
 [preseq](https://github.com/smithlabcode/preseq) into falco.
 
-The original duplication analysis method is still implemented in falco v2.0, but
-it must be turned on at build time by supplying `-DORIGINAL_DUPS=on` to
-`cmake`. It's still there because I used a trivial implementation and helps me
-debug. Here's how the old and new analyses differ, to the best of my
-understanding. I'm happy to be corrected.
+The original duplication analysis method is still implemented in falco v2.0, and
+can be turned on with `--orig-dups`. Here's how the old and new analyses differ.
 
 #### Original method
 - 100,000 unique reads are hashed and counted (first 50nt of each read).
@@ -182,12 +167,6 @@ understanding. I'm happy to be corrected.
   will appear as though they are randomly sampled due to fluctaions in thread
   speed changing which 1M reads are hashed.
 
-Building with `-DORIGINAL_DUPS=on` enables the original duplication analysis,
-but it only gives the same results when run using one thread. The underlying
-issues with this duplication analysis were part of the motivation for preseq, so
-I'm reluctant to put effort towards improving how the original duplication
-analysis is implemented.
-
 ## Citing falco
 
 If falco was helpful for your research, you can cite us as follows:

cmake/tests.cmake

@@ -32,3 +32,4 @@ add_test(NAME "Tiles" COMMAND bash test_scripts/tiles.sh)
 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)

data/test_data/md5sum.txt

@@ -11,3 +11,4 @@ ad21a5d9eb4cac3018c61fd7afdcbaf4  bam_mult_out/bam_2/fastqc_data.txt
 4cf6ba5d87aa553a64ce5427c18ee85c  bam_tiles_out/bam_1/fastqc_data.txt
 ec81607e6c63f03e04c4c49588be816e  groups_out/bam_1/fastqc_data.txt
 a3a0045355760f0f1ec5b8ca4fc1753e  sam_out/sam_1/fastqc_data.txt
+4cf6ba5d87aa553a64ce5427c18ee85c  bam_orig_dups_out/bam_1/fastqc_data.txt

data/test_scripts/orig_dups.sh

@@ -0,0 +1,19 @@
+#!/usr/bin/env bash
+# SPDX-License-Identifier: MIT
+
+prog=./falco
+infile=test_data/bam_1.bam
+outdir=bam_orig_dups_out
+if [[ -e "${infile}" ]]; then
+    ${prog} --orig-dups -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/bam_file.cpp

@@ -2,7 +2,9 @@
 
 #include "bam_file.hpp"
 #include "bamrec.hpp"
+#include "duplication_results.hpp"
 #include "falco_utils.hpp"
+#include "falco_word.hpp"
 #include "task_queue.hpp"
 
 #include <htslib/hfile.h>
@@ -58,6 +60,52 @@ estimate_n_reads_bam(const std::string &filename)
   return {estimate, total_read_len / n_reads, filesize};
 }
 
+[[nodiscard]] auto
+init_dups(const std::string &filename,
+          const std::uint64_t n_unique) -> dups_map_t {
+  static constexpr auto complem = [](const auto x) {
+    return "TNGNNNCNNNNNNNNNNNNA"[x - 'A'];
+  };
+  static constexpr auto n_unique_multiplier = 10;
+  std::unique_ptr<htsFile, int (*)(htsFile *)> f(
+    hts_open(std::data(filename), "r"), &hts_close);
+  if (!f)
+    throw std::system_error(std::make_error_code(std::errc(errno)),
+                            "failed to open file: " + filename);
+  std::unique_ptr<sam_hdr_t, void (*)(sam_hdr_t *)> h(sam_hdr_read(f.get()),
+                                                      &sam_hdr_destroy);
+  if (!h)
+    throw std::system_error(std::make_error_code(std::errc(errno)),
+                            "failed to read header: " + filename);
+  // NOLINTNEXTLINE(cppcoreguidelines-pro-type-union-access)
+  std::unique_ptr<bam1_t, void (*)(bam1_t *)> rec(bam_init1(), &bam_destroy1);
+
+  dups_map_t dups;
+  std::vector<char> buffer;
+  int r{};
+  std::uint64_t n_reads{};
+  const auto max_n_reads = n_unique_multiplier * n_unique;
+  while (n_reads++ < max_n_reads && std::size(dups) < n_unique &&
+         (r = sam_read1(f.get(), h.get(), rec.get())) >= 0) {
+    const auto l_qseq = rec->core.l_qseq;
+    const auto seq = bam_get_seq(rec.get());
+    if (std::ssize(buffer) < l_qseq)
+      buffer.resize(l_qseq);
+    for (auto i = 0; i < l_qseq; ++i)
+      buffer[i] = seq_nt16_str[bam_seqi(seq, i)];
+    if (bam_is_rev(rec.get())) {
+      std::reverse(std::begin(buffer), std::begin(buffer) + l_qseq);
+      std::transform(std::cbegin(buffer), std::cbegin(buffer) + l_qseq,
+                     std::begin(buffer), complem);
+    }
+    ++dups[falco_word(std::data(buffer), l_qseq)];
+  }
+  if (r < -1)  // error
+    throw std::system_error(std::make_error_code(std::errc(errno)),
+                            "error reading bam record from: " + filename);
+  return dups;
+}
+
 auto
 bam_file::load_next(const std::int32_t file_id,  //
                     task_queue &tq,              //

src/bam_file.hpp

@@ -6,6 +6,7 @@
 #include "bam_header.hpp"
 #include "bgzf_block.hpp"
 #include "bgzf_reader.hpp"
+#include "duplication_results.hpp"
 
 #include <atomic>
 #include <cstdint>
@@ -119,6 +120,10 @@ class bam_file {
 estimate_n_reads_bam(const std::string &filename)
   -> std::tuple<std::uint64_t, std::uint64_t, std::int64_t>;
 
+[[nodiscard]] auto
+init_dups(const std::string &filename,
+          const std::uint64_t n_unique) -> dups_map_t;
+
 inline auto
 make_tasks(bam_file &reads_file,          //
            const std::int64_t n_threads,  //

src/duplication_results.cpp

@@ -10,12 +10,12 @@
 #include "html.hpp"
 #include "run_mode.hpp"
 
-#include "boost/boost_unordered.hpp"
-
 #define FMT_HEADER_ONLY
 #include "fmt/format.h"
 #include "fmt/ranges.h"
 
+#include "boost/boost_unordered.hpp"
+
 #include <algorithm>
 #include <array>
 #include <cstdint>
@@ -107,6 +107,21 @@ duplication_results::get_overrepresented(const std::uint64_t n_reads) const
   return ret;
 }
 
+auto
+duplication_results::initialize(const run_mode &mode, const file_info &info,
+                                const dups_init_t &dups_init) -> void {
+  read_skip =
+    info.n_reads_est < max_n_reads_total
+      ? 0
+      : static_cast<std::int32_t>(info.n_reads_est / max_n_reads_total);
+  dups = dups_init.dups_zero;
+  count_at_limit = dups_init.count_at_limit;
+  if (!mode.do_dup_analysis()) {
+    // ADS: disabling dups analysis; does nothing for original dups
+    read_idx = std::numeric_limits<std::int64_t>::max();
+  }
+}
+
 auto
 duplication_results::initialize(const run_mode &mode,
                                 const file_info &info) -> void {
@@ -117,17 +132,13 @@ duplication_results::initialize(const run_mode &mode,
   if (!mode.do_dup_analysis()) {
     // ADS: disabling dups analysis
     read_idx = std::numeric_limits<std::int64_t>::max();
-    // ADS: for ORIGINAL_DUPS this still does lots of work
   }
 }
 
 auto
 duplication_results::add_and_consume(duplication_results &rhs) -> void {
-#ifndef ORIGINAL_DUPS
   for (const auto &[k, v] : rhs.dups)
     dups[k] += v;
-#endif  // ORIGINAL_DUPS
-  // ADS: ORIGINAL_DUPS mode just uses dups from one of two duplication_results
   rhs.release();
 }
 
@@ -137,9 +148,8 @@ get_grade_overrepresented(const std::uint64_t n_reads,
   static constexpr auto label = "overrepresented";
   const auto max_n_obs =
     dr.dups.empty() ? 0LU : std::ranges::max(std::views::values(dr.dups));
-  return grader_set::get_grade(
-    label,
-    as_frac(max_n_obs, std::max(static_cast<std::uint64_t>(1), n_reads)));
+  return grader_set::get_grade(label,
+                               as_frac(max_n_obs, std::max(1LU, n_reads)));
 }
 
 [[nodiscard]] auto
@@ -174,8 +184,6 @@ make_bins(const auto &breaks, const auto &hist) {
   return binned;
 }
 
-#ifdef ORIGINAL_DUPS
-
 // ADS: for original dups, from FastQC extrapolation of dup counts.
 [[nodiscard]] auto
 get_corrected_count(const std::uint64_t count_at_limit,
@@ -209,7 +217,8 @@ get_corrected_count(const std::uint64_t count_at_limit,
 }
 
 [[nodiscard]] auto
-duplication_results::get_dups_summary() const -> dup_summary_t {
+duplication_results::get_dups_summary(const std::uint64_t n_reads) const
+  -> dup_summary_t {
   if (dups.empty())
     return {};
   const auto max_dup = std::ranges::max(std::views::values(dups));
@@ -218,22 +227,20 @@ duplication_results::get_dups_summary() const -> dup_summary_t {
     ++hist_dedup[n_copies];
   for (auto [idx, val] : std::views::enumerate(hist_dedup))
     val = static_cast<std::uint64_t>(
-      get_corrected_count(count_at_limit, read_idx, idx, val));
+      get_corrected_count(count_at_limit, n_reads, idx, val));
   auto hist_mass =
     std::views::transform(
       std::views::enumerate(hist_dedup),
       [](const auto x) { return std::get<0>(x) * std::get<1>(x); }) |
     std::ranges::to<std::vector>();
   return dup_summary_t{
     max_dup,
-    get_n_counted_reads(),  // number of counted reads
+    n_reads,                // number of counted reads
     std::move(hist_mass),   // move to avoid copying when making tuple
     std::move(hist_dedup),  // move to avoid copying when making tuple
   };
 }
 
-#else  // NOT ORIGINAL_DUPS
-
 [[nodiscard]] auto
 duplication_results::get_dups_summary() const -> dup_summary_t {
   if (dups.empty())
@@ -249,14 +256,12 @@ duplication_results::get_dups_summary() const -> dup_summary_t {
     std::ranges::to<std::vector>();
   return dup_summary_t{
     max_dup,
-    get_n_counted_reads(),  // total counted reads
+    get_n_counted_reads(),  // number of counted reads
     std::move(hist_mass),   // move to avoid copying when making tuple
     std::move(hist_dedup),  // move to avoid copying when making tuple
   };
 }
 
-#endif
-
 [[nodiscard]] auto
 get_grade_duplication(const dup_summary_t &summary) -> std::string {
   // ADS: duplication grade is based on summary stat of estimated number of