IUPAC symbols for orig-dups mode: non-ACGT (case insensitive) are converted to N

CMakeLists.txt

@@ -130,6 +130,8 @@ target_link_libraries(falco_file_format PRIVATE HTSLIB::HTSLIB)
 target_link_libraries(tile_processor PRIVATE HTSLIB::HTSLIB)
 target_link_libraries(quality_score PRIVATE HTSLIB::HTSLIB)
 target_link_libraries(fastq_gz_file PRIVATE HTSLIB::HTSLIB)
+target_link_libraries(original_duplicates PRIVATE HTSLIB::HTSLIB)
+
 if(USE_ISAL)
   target_link_libraries(fastq_gz_file PUBLIC ISAL::ISAL)
 endif()

src/bam_file.cpp

@@ -2,9 +2,7 @@
 
 #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/bgzf.h>  // IWYU pragma: keep
@@ -65,52 +63,6 @@ estimate_n_reads_bam(const std::string &filename)
   return {n_reads_estimate, read_len_estimate, 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;
-  const std::unique_ptr<htsFile, int (*)(htsFile *)> in(
-    hts_open(std::data(filename), "r"), &hts_close);
-  if (!in)
-    throw std::system_error(std::make_error_code(std::errc(errno)),
-                            "failed to open file: " + filename);
-  const std::unique_ptr<sam_hdr_t, void (*)(sam_hdr_t *)> h(
-    sam_hdr_read(in.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(in.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,7 +6,6 @@
 #include "bam_header.hpp"
 #include "bgzf_block.hpp"
 #include "bgzf_reader.hpp"
-#include "duplication_results.hpp"
 
 #include <atomic>
 #include <cstdint>
@@ -129,10 +128,6 @@ 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

@@ -27,6 +27,7 @@
 #include <numeric>
 #include <ranges>
 #include <span>
+#include <stdexcept>
 #include <string>
 #include <tuple>
 #include <unordered_map>
@@ -236,6 +237,8 @@ duplication_results::get_dups_summary(const std::uint64_t n_reads) const
   std::vector<std::uint64_t> hist_dedup(max_dup + 1);
   for (const auto n_copies : std::views::values(dups))
     ++hist_dedup[n_copies];
+  if (hist_dedup[0] > 0)
+    throw std::runtime_error("Non-zero count of indexed reads seen 0 times");
   for (auto [idx, val] : falco::views::enumerate(hist_dedup))
     val = static_cast<std::uint64_t>(
       get_corrected_count(count_at_limit, n_reads, idx, val));

src/fastq_bgzf_file.cpp

@@ -1,15 +1,13 @@
 // SPDX-License-Identifier: MIT; Copyright 2026 Andrew D Smith
 
 #include "fastq_bgzf_file.hpp"
-#include "duplication_results.hpp"
 #include "falco_utils.hpp"
-#include "falco_word.hpp"
 #include "fqrec.hpp"
 #include "task_queue.hpp"
 
 #include <htslib/bgzf.h>
 #include <htslib/hfile.h>
-#include <htslib/sam.h>
+#include <htslib/sam.h>  // IWYU pragma: keep
 
 #include <algorithm>
 #include <cassert>
@@ -55,38 +53,6 @@ estimate_n_reads_fastq_bgzf(const std::string &filename)
   return {static_cast<std::uint64_t>(n_reads_est), read_len_est, filesize};
 }
 
-[[nodiscard]] auto
-init_dups_fq(const std::string &filename, const std::uint64_t n_unique)
-  -> dups_map_t {
-  static constexpr auto n_unique_multiplier = 10;
-  const std::unique_ptr<BGZF, int (*)(BGZF *)> in(
-    bgzf_open(std::data(filename), "r"), &bgzf_close);
-  if (!in)
-    throw std::system_error(std::make_error_code(std::errc(errno)),
-                            "failed to open file: " + filename);
-  dups_map_t dups;
-  kstring_t s = KS_INITIALIZE;
-  int r{};
-  std::uint64_t n_lines{};
-  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 = bgzf_getline(in.get(), '\n', &s)) >= 0) {
-    if (n_lines % 4 == 1) {
-      const auto l_qseq = ks_len(&s);
-      const auto seq = ks_c_str(&s);
-      ++dups[falco_word(seq, l_qseq)];
-      ++n_reads;
-    }
-    ++n_lines;
-  }
-  if (r < -1)  // error
-    throw std::system_error(std::make_error_code(std::errc(errno)),
-                            "error reading bam record from: " + filename);
-  ks_free(&s);
-  return dups;
-}
-
 auto
 fastq_bgzf_file::load_next(const std::int32_t file_id,
                            task_queue &tq,

src/fastq_bgzf_file.hpp

@@ -5,7 +5,6 @@
 
 #include "bgzf_block.hpp"
 #include "bgzf_reader.hpp"
-#include "duplication_results.hpp"
 
 #include <atomic>
 #include <cstdint>
@@ -128,10 +127,6 @@ class fastq_bgzf_file {
 estimate_n_reads_fastq_bgzf(const std::string &filename)
   -> std::tuple<std::uint64_t, std::uint64_t, std::int64_t>;
 
-[[nodiscard]] auto
-init_dups_fq(const std::string &filename, const std::uint64_t n_unique)
-  -> dups_map_t;
-
 inline auto
 make_tasks(fastq_bgzf_file &reads_file,
            const std::int64_t n_threads,

src/original_duplicates.cpp

@@ -2,29 +2,132 @@
 
 #include "original_duplicates.hpp"
 
-#include "bam_file.hpp"
 #include "duplication_results.hpp"
 #include "falco_file_format.hpp"
-#include "fastq_bgzf_file.hpp"
+#include "falco_word.hpp"
 #include "file_info.hpp"
 
 #include "boost/boost_unordered.hpp"
 
+#include <htslib/bgzf.h>  // IWYU pragma: keep
+#include <htslib/sam.h>
+
 #include <algorithm>
 #include <cassert>
+#include <cerrno>
 #include <iterator>
+#include <memory>
 #include <mutex>
 #include <ranges>
+#include <span>
 #include <string>
+#include <system_error>
 #include <thread>
 #include <utility>  // for std::move
 #include <vector>
 
+namespace {
+[[nodiscard]] auto
+init_orig_dups_bam(const std::string &filename, const std::uint64_t n_unique)
+  -> dups_map_t {
+  // ADS: n_unique_multiplier: contributes to max reads to examine, multiplying
+  // the n_unique that is the target; impossible to know how many needed to hit
+  // the target, so try to overshoot.
+  static constexpr auto n_unique_multiplier = 10;
+  static constexpr auto complem = [](char &x) {
+    assert(x >= 'A' && x <= 'T');
+    x = "TNGNNNCNNNNNNNNNNNNA"[x - 'A'];
+  };
+  static constexpr auto to_acgtn = [](const char x) {
+    static constexpr char mask = 15;
+    static constexpr std::span seq_nt16_str_acgtn = "NACNGNNNTNNNNNNN";
+    return seq_nt16_str_acgtn[x & mask];
+  };
+  const std::unique_ptr<htsFile, int (*)(htsFile *)> in(
+    hts_open(std::data(filename), "r"), &hts_close);
+  if (!in)
+    throw std::system_error(std::make_error_code(std::errc(errno)),
+                            "failed to open file: " + filename);
+  const std::unique_ptr<sam_hdr_t, void (*)(sam_hdr_t *)> h(
+    sam_hdr_read(in.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(in.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] = to_acgtn(bam_seqi(seq, i));
+    if (bam_is_rev(rec.get())) {
+      std::reverse(std::begin(buffer), std::begin(buffer) + l_qseq);
+      std::for_each_n(std::begin(buffer), l_qseq, 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;
+}
+
+[[nodiscard]] auto
+init_orig_dups_fq(const std::string &filename, const std::uint64_t n_unique)
+  -> dups_map_t {
+  // ADS: n_unique_multiplier: contributes to max reads to examine, multiplying
+  // the n_unique that is the target; impossible to know how many needed to hit
+  // the target, so try to overshoot.
+  static constexpr auto n_unique_multiplier = 10;
+  static constexpr auto to_acgtn = [](const auto c) {
+    static constexpr auto mask = 31;
+    static constexpr std::span to_acgtn_sp = "NANCNNNGNNNNNNNNNNNNTNNNNNNNNNNN";
+    return to_acgtn_sp[c & mask];
+  };
+  const std::unique_ptr<BGZF, int (*)(BGZF *)> in(
+    bgzf_open(std::data(filename), "r"), &bgzf_close);
+  if (!in)
+    throw std::system_error(std::make_error_code(std::errc(errno)),
+                            "failed to open file: " + filename);
+  dups_map_t dups;
+  kstring_t s = KS_INITIALIZE;
+  int r{};
+  std::uint64_t n_lines{};
+  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 = bgzf_getline(in.get(), '\n', &s)) >= 0) {
+    if (n_lines % 4 == 1) {
+      const auto l_qseq = ks_len(&s);
+      std::ranges::transform(std::span(s.s, l_qseq), s.s, to_acgtn);
+      const auto seq = ks_c_str(&s);
+      ++dups[falco_word(seq, l_qseq)];
+      ++n_reads;
+    }
+    ++n_lines;
+  }
+  if (r < -1)  // error
+    throw std::system_error(std::make_error_code(std::errc(errno)),
+                            "error reading bam record from: " + filename);
+  ks_free(&s);
+  return dups;
+}
+}  // namespace
+
 [[nodiscard]] auto
-initialize_original_duplicates(
-  const std::vector<std::string> &infiles,
-  const std::vector<file_info> &infos,
-  const std::uint32_t n_threads) -> std::vector<dups_init_t> {
+initialize_original_duplicates(const std::vector<std::string> &infiles,
+                               const std::vector<file_info> &infos,
+                               const std::uint32_t n_threads)
+  -> std::vector<dups_init_t> {
   assert(std::size(infiles) == std::size(infos));
   const std::uint64_t n_files = std::size(infiles);
   std::vector<dups_map_t> dups(n_files);
@@ -47,10 +150,10 @@ initialize_original_duplicates(
         }
         dups[file_id] =
           falco::is_mapped_reads(infos[file_id].format)
-            ? init_dups(infiles[file_id],
-                        duplication_results::max_reads_to_hash)
-            : init_dups_fq(infiles[file_id],
-                           duplication_results::max_reads_to_hash);
+            ? init_orig_dups_bam(infiles[file_id],
+                                 duplication_results::max_reads_to_hash)
+            : init_orig_dups_fq(infiles[file_id],
+                                duplication_results::max_reads_to_hash);
       }
     });
   std::ranges::for_each(workers, [](auto &w) { w.join(); });