From b97a8762fa9ccba1f046e0c68128a8166f8d4f85 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 1 Nov 2026 16:45:28 -0800 Subject: [PATCH 1/7] src/original_duplicates.cpp: moving code to initialize original duplicates analysis here --- src/original_duplicates.cpp | 88 ++++++++++++++++++++++++++++++++++--- 1 file changed, 82 insertions(+), 6 deletions(-) diff --git a/src/original_duplicates.cpp b/src/original_duplicates.cpp index d0b960b..6fd00e8 100644 --- a/src/original_duplicates.cpp +++ b/src/original_duplicates.cpp @@ -2,10 +2,8 @@ #include "original_duplicates.hpp" -#include "bam_file.hpp" #include "duplication_results.hpp" #include "falco_file_format.hpp" -#include "fastq_bgzf_file.hpp" #include "file_info.hpp" #include "boost/boost_unordered.hpp" @@ -21,10 +19,88 @@ #include [[nodiscard]] auto -initialize_original_duplicates( - const std::vector &infiles, - const std::vector &infos, - const std::uint32_t n_threads) -> std::vector { +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 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 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 rec(bam_init1(), &bam_destroy1); + + dups_map_t dups; + std::vector 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; +} + +[[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 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; +} + +[[nodiscard]] auto +initialize_original_duplicates(const std::vector &infiles, + const std::vector &infos, + const std::uint32_t n_threads) + -> std::vector { assert(std::size(infiles) == std::size(infos)); const std::uint64_t n_files = std::size(infiles); std::vector dups(n_files); From 0e0bfc3c5c9afa1147b0e5d4cef298d55df0f38f Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 1 Nov 2026 16:45:49 -0800 Subject: [PATCH 2/7] src/fastq_bgzf_file.hcpp and src/bam_file.hcpp: moving code for initializing original duplication data structures to original_duplicates.cpp --- src/bam_file.cpp | 46 ----------------------------------------- src/bam_file.hpp | 4 ---- src/fastq_bgzf_file.cpp | 32 ---------------------------- src/fastq_bgzf_file.hpp | 4 ---- 4 files changed, 86 deletions(-) diff --git a/src/bam_file.cpp b/src/bam_file.cpp index 90533a8..d160a56 100644 --- a/src/bam_file.cpp +++ b/src/bam_file.cpp @@ -65,52 +65,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 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 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 rec(bam_init1(), &bam_destroy1); - - dups_map_t dups; - std::vector 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, diff --git a/src/bam_file.hpp b/src/bam_file.hpp index 97a2473..10c5458 100644 --- a/src/bam_file.hpp +++ b/src/bam_file.hpp @@ -129,10 +129,6 @@ class bam_file { estimate_n_reads_bam(const std::string &filename) -> std::tuple; -[[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, diff --git a/src/fastq_bgzf_file.cpp b/src/fastq_bgzf_file.cpp index 719e928..0c3605c 100644 --- a/src/fastq_bgzf_file.cpp +++ b/src/fastq_bgzf_file.cpp @@ -55,38 +55,6 @@ estimate_n_reads_fastq_bgzf(const std::string &filename) return {static_cast(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 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, diff --git a/src/fastq_bgzf_file.hpp b/src/fastq_bgzf_file.hpp index 0023291..26a2169 100644 --- a/src/fastq_bgzf_file.hpp +++ b/src/fastq_bgzf_file.hpp @@ -128,10 +128,6 @@ class fastq_bgzf_file { estimate_n_reads_fastq_bgzf(const std::string &filename) -> std::tuple; -[[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, From 754dfd79119f7d34ca14df51806001b81a0953b6 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 1 Nov 2026 16:47:34 -0800 Subject: [PATCH 3/7] src/original_duplicates.cpp: adding forgotten headers --- src/original_duplicates.cpp | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/src/original_duplicates.cpp b/src/original_duplicates.cpp index 6fd00e8..7f369a3 100644 --- a/src/original_duplicates.cpp +++ b/src/original_duplicates.cpp @@ -8,6 +8,10 @@ #include "boost/boost_unordered.hpp" +#include // IWYU pragma: keep +#include +#include + #include #include #include @@ -19,7 +23,7 @@ #include [[nodiscard]] auto -init_dups(const std::string &filename, const std::uint64_t n_unique) +init_orig_dups_bam(const std::string &filename, const std::uint64_t n_unique) -> dups_map_t { static constexpr auto complem = [](const auto x) { return "TNGNNNCNNNNNNNNNNNNA"[x - 'A']; @@ -65,7 +69,7 @@ init_dups(const std::string &filename, const std::uint64_t n_unique) } [[nodiscard]] auto -init_dups_fq(const std::string &filename, const std::uint64_t n_unique) +init_orig_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 in( @@ -123,10 +127,10 @@ initialize_original_duplicates(const std::vector &infiles, } 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(); }); From feccee05cb99a4317c1a8c4d62cd5466fcbf9bc4 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 1 Nov 2026 16:55:14 -0800 Subject: [PATCH 4/7] src/original_duplicates.cpp: adding to_acgtn in both the init_orig_dups_bam and init_orig_dups_fq to convert all non-ACGT letters to N prior to forming a falco_word when initializing the unique 'reads' for original duplicates analysis --- src/original_duplicates.cpp | 26 ++++++++++++++++++++++---- 1 file changed, 22 insertions(+), 4 deletions(-) diff --git a/src/original_duplicates.cpp b/src/original_duplicates.cpp index 7f369a3..f25d54a 100644 --- a/src/original_duplicates.cpp +++ b/src/original_duplicates.cpp @@ -25,10 +25,19 @@ [[nodiscard]] auto init_orig_dups_bam(const std::string &filename, const std::uint64_t n_unique) -> dups_map_t { - static constexpr auto complem = [](const auto x) { - return "TNGNNNCNNNNNNNNNNNNA"[x - 'A']; - }; + // 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 in( hts_open(std::data(filename), "r"), &hts_close); if (!in) @@ -54,7 +63,7 @@ init_orig_dups_bam(const std::string &filename, const std::uint64_t n_unique) 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)]; + 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::transform(std::cbegin(buffer), std::cbegin(buffer) + l_qseq, @@ -71,7 +80,15 @@ init_orig_dups_bam(const std::string &filename, const std::uint64_t n_unique) [[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 in( bgzf_open(std::data(filename), "r"), &bgzf_close); if (!in) @@ -87,6 +104,7 @@ init_orig_dups_fq(const std::string &filename, const std::uint64_t 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; From 03bcc802d59099e7cf25ace004f1258dc25cea58 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 1 Nov 2026 16:56:38 -0800 Subject: [PATCH 5/7] CMakeLists.txt: original_duplicates object links to HTSLIB so that the htslib headers will be visibile when making original_duplicates.o --- CMakeLists.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index b4c0948..c7a7d41 100644 --- a/CMakeLists.txt +++ b/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() From 91e83cbf26282c99c4cf41d25b15f62fd83e6b4d Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 1 Nov 2026 17:00:44 -0800 Subject: [PATCH 6/7] src/original_duplicates.cpp: fix changing transform to for_each_n when complementing the already reversed read sequence for revcomp mappers --- src/original_duplicates.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/original_duplicates.cpp b/src/original_duplicates.cpp index f25d54a..66c71da 100644 --- a/src/original_duplicates.cpp +++ b/src/original_duplicates.cpp @@ -66,8 +66,7 @@ init_orig_dups_bam(const std::string &filename, const std::uint64_t n_unique) 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::transform(std::cbegin(buffer), std::cbegin(buffer) + l_qseq, - std::begin(buffer), complem); + std::for_each_n(std::begin(buffer), l_qseq, complem); } ++dups[falco_word(std::data(buffer), l_qseq)]; } From 8a2a40e2a241b70175d7e5feb92547a1baabc240 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 1 Nov 2026 17:43:16 -0800 Subject: [PATCH 7/7] static analysis --- src/bam_file.cpp | 2 -- src/bam_file.hpp | 1 - src/duplication_results.cpp | 3 +++ src/fastq_bgzf_file.cpp | 4 +--- src/fastq_bgzf_file.hpp | 1 - src/original_duplicates.cpp | 8 +++++++- 6 files changed, 11 insertions(+), 8 deletions(-) diff --git a/src/bam_file.cpp b/src/bam_file.cpp index d160a56..b8a28f7 100644 --- a/src/bam_file.cpp +++ b/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 // IWYU pragma: keep diff --git a/src/bam_file.hpp b/src/bam_file.hpp index 10c5458..f781c89 100644 --- a/src/bam_file.hpp +++ b/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 #include diff --git a/src/duplication_results.cpp b/src/duplication_results.cpp index 0a33242..ff269dc 100644 --- a/src/duplication_results.cpp +++ b/src/duplication_results.cpp @@ -27,6 +27,7 @@ #include #include #include +#include #include #include #include @@ -236,6 +237,8 @@ duplication_results::get_dups_summary(const std::uint64_t n_reads) const std::vector 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( get_corrected_count(count_at_limit, n_reads, idx, val)); diff --git a/src/fastq_bgzf_file.cpp b/src/fastq_bgzf_file.cpp index 0c3605c..68536e5 100644 --- a/src/fastq_bgzf_file.cpp +++ b/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 #include -#include +#include // IWYU pragma: keep #include #include diff --git a/src/fastq_bgzf_file.hpp b/src/fastq_bgzf_file.hpp index 26a2169..8656e2d 100644 --- a/src/fastq_bgzf_file.hpp +++ b/src/fastq_bgzf_file.hpp @@ -5,7 +5,6 @@ #include "bgzf_block.hpp" #include "bgzf_reader.hpp" -#include "duplication_results.hpp" #include #include diff --git a/src/original_duplicates.cpp b/src/original_duplicates.cpp index 66c71da..44629d5 100644 --- a/src/original_duplicates.cpp +++ b/src/original_duplicates.cpp @@ -4,24 +4,29 @@ #include "duplication_results.hpp" #include "falco_file_format.hpp" +#include "falco_word.hpp" #include "file_info.hpp" #include "boost/boost_unordered.hpp" #include // IWYU pragma: keep -#include #include #include #include +#include #include +#include #include #include +#include #include +#include #include #include // for std::move #include +namespace { [[nodiscard]] auto init_orig_dups_bam(const std::string &filename, const std::uint64_t n_unique) -> dups_map_t { @@ -116,6 +121,7 @@ init_orig_dups_fq(const std::string &filename, const std::uint64_t n_unique) ks_free(&s); return dups; } +} // namespace [[nodiscard]] auto initialize_original_duplicates(const std::vector &infiles,