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() diff --git a/src/bam_file.cpp b/src/bam_file.cpp index 90533a8..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 @@ -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 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..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 @@ -129,10 +128,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/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 719e928..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 @@ -55,38 +53,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..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 @@ -128,10 +127,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, diff --git a/src/original_duplicates.cpp b/src/original_duplicates.cpp index d0b960b..44629d5 100644 --- a/src/original_duplicates.cpp +++ b/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 // IWYU pragma: keep +#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 { + // 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) + 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] = 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 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 &infiles, - const std::vector &infos, - const std::uint32_t n_threads) -> std::vector { +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); @@ -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(); });