diff --git a/CMakeLists.txt b/CMakeLists.txt index c3aca1c..69615ed 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -199,6 +199,9 @@ 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 @@ -206,6 +209,7 @@ target_link_libraries(falco PUBLIC adapter_matcher adapter_set duplication_results + original_duplicates kmer_counter fastq_file fastq_gz_file diff --git a/README.md b/README.md index aaf2b03..c8ef267 100644 --- a/README.md +++ b/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: diff --git a/cmake/tests.cmake b/cmake/tests.cmake index 71a087a..d5cdbfa 100644 --- a/cmake/tests.cmake +++ b/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) diff --git a/data/test_data/md5sum.txt b/data/test_data/md5sum.txt index 6e82c5f..f828598 100644 --- a/data/test_data/md5sum.txt +++ b/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 diff --git a/data/test_scripts/orig_dups.sh b/data/test_scripts/orig_dups.sh new file mode 100644 index 0000000..2d3958a --- /dev/null +++ b/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 diff --git a/src/bam_file.cpp b/src/bam_file.cpp index 8e24bbc..8c7a32c 100644 --- a/src/bam_file.cpp +++ b/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 @@ -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 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 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 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(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, // diff --git a/src/bam_file.hpp b/src/bam_file.hpp index 330c7d8..606df5d 100644 --- a/src/bam_file.hpp +++ b/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 #include @@ -119,6 +120,10 @@ 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 35a9b58..a3abf52 100644 --- a/src/duplication_results.cpp +++ b/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 #include #include @@ -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(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::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::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(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,7 +227,7 @@ 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( - 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), @@ -226,14 +235,12 @@ duplication_results::get_dups_summary() const -> dup_summary_t { std::ranges::to(); 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(); 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 diff --git a/src/duplication_results.hpp b/src/duplication_results.hpp index 6710f7a..bda01b2 100644 --- a/src/duplication_results.hpp +++ b/src/duplication_results.hpp @@ -8,17 +8,18 @@ #include "boost/boost_unordered.hpp" #include +#include +#include +#include #include #include -#ifdef ORIGINAL_DUPS -#include -#endif // ORIGINAL_DUPS - class run_mode; struct file_grades; struct file_info; +using dups_map_t = boost::unordered_flat_map; + struct overrep_t { falco_word w; std::uint64_t n_obs{}; // number of observations @@ -33,21 +34,35 @@ struct dup_summary_t { std::vector hist_dedup; }; +struct dups_init_t { + std::int64_t count_at_limit{}; + dups_map_t dups_zero; + dups_init_t() = default; + dups_init_t(const dups_map_t &dups) { + const auto vals = dups | std::views::values; + count_at_limit = static_cast( + std::reduce(std::cbegin(vals), std::cend(vals))); + for (const auto &fw : dups | std::views::keys) + dups_zero.emplace(fw, 0); + } + operator bool() const { return !dups_zero.empty(); } +}; + struct duplication_results { static constexpr auto max_n_reads_total{1'000'000}; -#ifdef ORIGINAL_DUPS static constexpr auto max_reads_to_hash{100'000}; -#endif // ORIGINAL_DUPS static constexpr auto default_read_skip{10}; static constexpr auto overrep_cutoff = 0.001; -#ifdef ORIGINAL_DUPS std::int64_t count_at_limit{}; -#endif // ORIGINAL_DUPS std::int64_t read_skip{default_read_skip}; std::int64_t read_idx{}; boost::unordered_flat_map dups; + auto + initialize(const run_mode &mode, const file_info &info, + const dups_init_t &dups_init) -> void; + auto initialize(const run_mode &mode, const file_info &info) -> void; @@ -60,6 +75,9 @@ struct duplication_results { boost::unordered_flat_map().swap(dups); } + [[nodiscard]] auto + get_dups_summary(const std::uint64_t n_reads) const -> dup_summary_t; + [[nodiscard]] auto get_dups_summary() const -> dup_summary_t; @@ -70,17 +88,13 @@ struct duplication_results { auto add_and_consume(duplication_results &rhs) -> void; -#ifdef ORIGINAL_DUPS auto - count_seqs(const auto seq_itr, const auto sz) { - ++read_idx; - const auto fw = falco_word(seq_itr, sz); - if (std::size(dups) < max_reads_to_hash || dups.contains(fw)) - ++dups[fw]; - else if (count_at_limit == 0) - count_at_limit = read_idx; + count_seqs_orig(const auto seq_itr, const auto sz) { + auto itr = dups.find(falco_word(seq_itr, sz)); + if (itr != std::cend(dups)) + ++itr->second; } -#else // NOT ORIGINAL_DUPS + auto count_seqs(const auto seq_itr, const auto sz) { if (read_idx-- == 0) [[unlikely]] { @@ -88,7 +102,6 @@ struct duplication_results { ++dups[falco_word(seq_itr, sz)]; } } -#endif // ORIGINAL_DUPS }; [[nodiscard]] auto diff --git a/src/falco.cpp b/src/falco.cpp index 8e22aac..042428b 100644 --- a/src/falco.cpp +++ b/src/falco.cpp @@ -41,6 +41,7 @@ Use these as templates. Copy and modify them to customize your analysis. #include "fastq_file.hpp" #include "fastq_gz_file.hpp" #include "get_binary_dir.hpp" +#include "original_duplicates.hpp" #include "quality_score.hpp" #include "reads_file.hpp" // IWYU pragma: keep #include "results_collector.hpp" @@ -83,16 +84,18 @@ write_file(const auto &filename, const auto &data) { } static auto -write_output(const run_mode &mode, std::vector &infos, - const std::vector &outdirs, - std::vector &results) { +write_output( + const run_mode &mode, std::vector &infos, + const std::vector &outdirs, + // NOLINTNEXTLINE(cppcoreguidelines-rvalue-reference-param-not-moved) + std::vector &&results) { static constexpr auto report_filename = "fastqc_data.txt"; static constexpr auto html_filename = "fastqc_report.html"; static constexpr auto summary_filename = "summary.txt"; for (const auto [result, info, outdir] : std::views::zip(results, infos, outdirs)) { const auto outdir_path = std::filesystem::path{outdir}; - const auto summary = results_summary(result, mode, info); + const auto summary = results_summary(std::move(result), mode, info); write_file(outdir_path / report_filename, summary.get_report()); write_file(outdir_path / html_filename, summary.get_html()); write_file(outdir_path / summary_filename, summary.get_summary()); @@ -216,8 +219,10 @@ main(int argc, char *argv[]) { int do_kmers{}; int do_dup_analysis{}; int do_adap{}; + int do_groups{}; int do_bisulfite{}; + int do_original_dups{}; std::uint32_t n_threads{1}; std::uint32_t max_read_length{}; @@ -299,6 +304,12 @@ main(int argc, char *argv[]) { app.add_flag("--bisulfite", do_bisulfite, "Assume bisulfite when grading sequence content") ->option_text(" "); + const auto orig_dups_opt = + app.add_flag("--orig-dups", [&](const auto x) { + do_original_dups = x; + do_dup_analysis = 1; + }, "Use original duplication mode (turns dups on)") + ->option_text(" "); app.add_flag("--groups", do_groups, "Group base positions in output") ->option_text(" "); app.add_flag("--tiles,!--no-tiles", do_tiles, @@ -306,6 +317,7 @@ main(int argc, char *argv[]) { ->option_text(" "); app.add_flag("--dups,!--no-dups", do_dup_analysis, "Toggle duplication/overrep analysis (default: on)") + ->excludes(orig_dups_opt) ->option_text(" "); app.add_flag("--adap,!--no-adap", do_adap, "Toggle adapter analysis (default: on)") @@ -340,6 +352,7 @@ main(int argc, char *argv[]) { mode.set_do_tiles(do_tiles); mode.set_do_groups(do_groups); mode.set_do_bisulfite(do_bisulfite); + mode.set_do_original_dups(do_original_dups); mode.set_unassigned(); const auto outdirs = make_outdirs(infiles, outdir); @@ -401,9 +414,14 @@ main(int argc, char *argv[]) { }); } - auto reads_files = make_reads_files(infos, infiles, buffer_size); - auto results = analyze(n_threads, mode, infos, std::move(reads_files)); - write_output(mode, infos, outdirs, results); + std::vector dups = + do_original_dups + ? initialize_original_duplicates(infiles, infos, n_threads) + : std::vector{}; + 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( diff --git a/src/falco_analyzer.cpp b/src/falco_analyzer.cpp index 1ed3132..549f9db 100644 --- a/src/falco_analyzer.cpp +++ b/src/falco_analyzer.cpp @@ -4,6 +4,7 @@ #include "bamrec.hpp" #include "bgzf_block.hpp" +#include "duplication_results.hpp" #include "file_info.hpp" #include "fqrec.hpp" #include "reads_file.hpp" @@ -26,21 +27,17 @@ class run_mode; [[nodiscard]] auto analyze(const std::uint32_t n_threads, const run_mode &mode, - std::vector &infos, std::vector reads_files) - -> std::vector { + std::vector &infos, std::vector reads_files, + std::vector dups_init) -> std::vector { assert(std::size(reads_files) == std::size(infos)); - const std::int32_t n_files = static_cast(std::size(infos)); + if (dups_init.empty()) + dups_init.resize(n_files); std::vector n_tasks(n_files); std::atomic_uint32_t n_active_files{static_cast(n_files)}; auto results = std::vector(n_threads, std::vector(n_files)); - // set per-file information used to do the analysis - for (auto &res : results) - for (const auto [file_id, info] : std::views::enumerate(infos)) - res[file_id].init(mode, info); - // add initial jobs so workers can work immediately task_queue tq; for (const auto file_id : std::views::iota(0, n_files)) @@ -51,6 +48,8 @@ analyze(const std::uint32_t n_threads, const run_mode &mode, for (const auto th_id : std::views::iota(0u, n_threads)) workers.emplace_back([&, n_threads, th_id] { auto &res = results[th_id]; + for (const auto [file_id, info] : std::views::enumerate(infos)) + res[file_id].init(mode, info, dups_init[file_id]); while (true) { auto tq_lock = tq.wait_and_acquire_lock(); if (tq.is_finished()) diff --git a/src/falco_analyzer.hpp b/src/falco_analyzer.hpp index 1ce1e6c..851d54c 100644 --- a/src/falco_analyzer.hpp +++ b/src/falco_analyzer.hpp @@ -3,18 +3,19 @@ #ifndef SRC_FALCO_ANALYZER_HPP_ #define SRC_FALCO_ANALYZER_HPP_ -#include "reads_file.hpp" // IWYU pragma: keep #include "results_collector.hpp" #include #include +class reads_file_t; class run_mode; +struct dups_init_t; struct file_info; [[nodiscard]] auto analyze(const std::uint32_t n_threads, const run_mode &mode, - std::vector &infos, std::vector reads_files) - -> std::vector; + std::vector &infos, std::vector reads_files, + std::vector dups_init) -> std::vector; #endif // SRC_FALCO_ANALYZER_HPP_ diff --git a/src/fastq_bgzf_file.cpp b/src/fastq_bgzf_file.cpp index 61222ef..cc8c9fb 100644 --- a/src/fastq_bgzf_file.cpp +++ b/src/fastq_bgzf_file.cpp @@ -1,12 +1,15 @@ // 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 #include @@ -48,6 +51,38 @@ 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; + std::unique_ptr f(bgzf_open(std::data(filename), "r"), + &bgzf_close); + if (!f) + 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(f.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, std::atomic_int32_t &n_tasks) -> void { diff --git a/src/fastq_bgzf_file.hpp b/src/fastq_bgzf_file.hpp index 38be33b..ad010c7 100644 --- a/src/fastq_bgzf_file.hpp +++ b/src/fastq_bgzf_file.hpp @@ -5,6 +5,7 @@ #include "bgzf_block.hpp" #include "bgzf_reader.hpp" +#include "duplication_results.hpp" #include #include @@ -121,6 +122,10 @@ 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 new file mode 100644 index 0000000..7aea3d3 --- /dev/null +++ b/src/original_duplicates.cpp @@ -0,0 +1,61 @@ +// SPDX-License-Identifier: MIT; Copyright 2026 Andrew D Smith + +#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" + +#include +#include +#include +#include +#include +#include +#include +#include + +[[nodiscard]] auto +initialize_original_duplicates( + const std::vector &infiles, + [[maybe_unused]] const std::vector &infos, + [[maybe_unused]] const std::uint32_t n_threads) -> std::vector { + assert(std::size(infiles) == std::size(infos)); + const auto n_files = std::size(infiles); + std::vector dups(n_files); + auto n_active_files = n_files; + std::vector workers; + const auto n_workers = + std::min(n_files, static_cast(n_threads)); + workers.reserve(n_workers); + std::mutex mtx; + // NOLINTNEXTLINE(clang-analyzer-deadcode.DeadStores) + for (const auto _ : std::views::iota(0u, n_workers)) + workers.emplace_back([&] { + std::uint64_t file_id{}; + while (true) { + { + std::scoped_lock l(mtx); + if (n_active_files == 0) + return; + file_id = --n_active_files; + } + 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); + } + }); + std::ranges::for_each(workers, [](auto &w) { w.join(); }); + std::vector ret; + ret.reserve(n_files); + for (const auto &d : dups) + ret.emplace_back(d); + return ret; +} diff --git a/src/original_duplicates.hpp b/src/original_duplicates.hpp new file mode 100644 index 0000000..1745a33 --- /dev/null +++ b/src/original_duplicates.hpp @@ -0,0 +1,19 @@ +// SPDX-License-Identifier: MIT; Copyright 2026 Andrew D Smith + +#ifndef SRC_ORIGINAL_DUPLICATES_HPP_ +#define SRC_ORIGINAL_DUPLICATES_HPP_ + +#include "duplication_results.hpp" + +#include +#include +#include + +struct file_info; + +[[nodiscard]] auto +initialize_original_duplicates( + const std::vector &infiles, const std::vector &infos, + const std::uint32_t n_threads) -> std::vector; + +#endif // SRC_ORIGINAL_DUPLICATES_HPP_ diff --git a/src/quality_score.cpp b/src/quality_score.cpp index a14e99f..3723228 100644 --- a/src/quality_score.cpp +++ b/src/quality_score.cpp @@ -7,7 +7,7 @@ #include #include -#include +#include // IWYU pragma: keep #include #include #include diff --git a/src/quality_score.hpp b/src/quality_score.hpp index 6ebe782..19db487 100644 --- a/src/quality_score.hpp +++ b/src/quality_score.hpp @@ -6,8 +6,11 @@ #include "nlohmann/json.hpp" #include +#include #include #include +#include // IWYU pragma: keep +#include // IWYU pragma: keep #include #include // IWYU pragma: keep diff --git a/src/results_collector.hpp b/src/results_collector.hpp index c2f0e4b..1f19f9a 100644 --- a/src/results_collector.hpp +++ b/src/results_collector.hpp @@ -48,6 +48,7 @@ struct alignas(assumed_page_size) results_collector { kmer_counter kc; bool do_tiles{}; bool do_kmers{}; + bool do_original_dups{}; results_collector() : lengths(1, 0) {} // in case all reads have length 0 @@ -63,8 +64,13 @@ struct alignas(assumed_page_size) results_collector { } auto - init(const run_mode &mode, const auto &info) { - dr.initialize(mode, info); + init(const run_mode &mode, const auto &info, const auto &dups_init = {}) { + if (mode.do_original_dups()) { + dr.initialize(mode, info, dups_init); + do_original_dups = true; + } + else + dr.initialize(mode, info); do_tiles = mode.do_tiles() && info.has_tiles; do_kmers = mode.do_kmers(); if (do_tiles) @@ -135,7 +141,10 @@ struct alignas(assumed_page_size) results_collector { count_ns(seq_itr, seq_end, n_counts); const auto tot = count_quals(get_qual(rec), get_qual_end(rec), qual_by_pos); ++qual_by_read[tot / read_len]; - dr.count_seqs(seq_itr, read_len); + if (do_original_dups) + dr.count_seqs_orig(seq_itr, read_len); + else + dr.count_seqs(seq_itr, read_len); am.match_adapters(seq_itr, read_len); if (do_tiles) tp(rec); diff --git a/src/results_summary.cpp b/src/results_summary.cpp index b11610a..6caeef4 100644 --- a/src/results_summary.cpp +++ b/src/results_summary.cpp @@ -56,15 +56,10 @@ results_summary::initialize() -> void { apply_groups(); // get summary structures -#ifdef ORIGINAL_DUPS - n_reads_for_dups = n_reads; - dup_summary = dr.get_dups_summary(); - dup_summary.n_reads = n_reads_for_dups; -#else // ORIGINAL_DUPS - n_reads_for_dups = dr.get_n_counted_reads(); - dup_summary = dr.get_dups_summary(); -#endif // ORIGINAL_DUPS - overrep = dr.get_overrepresented(n_reads_for_dups); + dup_summary = mode.do_original_dups() ? dr.get_dups_summary(n_reads) + : dr.get_dups_summary(); + + overrep = dr.get_overrepresented(dup_summary.n_reads); centered = tp.get_centered(); kmer_results = kc.get_kmer_results(); @@ -96,7 +91,7 @@ results_summary::assign_grades() -> void { if (mode.do_overrep()) grades.emplace("overrepresented", - get_grade_overrepresented(n_reads_for_dups, dr)); + get_grade_overrepresented(dup_summary.n_reads, dr)); if (mode.do_qual_base()) grades.emplace("quality_base", get_grade_quality_base(qual_by_pos)); diff --git a/src/results_summary.hpp b/src/results_summary.hpp index c300a0a..ee52bee 100644 --- a/src/results_summary.hpp +++ b/src/results_summary.hpp @@ -35,7 +35,6 @@ struct results_summary { std::vector qual_by_pos; falco::qual_array qual_by_read{}; duplication_results dr; - std::uint64_t n_reads_for_dups{}; std::vector overrep; dup_summary_t dup_summary; adapter_matcher am; @@ -49,7 +48,7 @@ struct results_summary { base_group_vec groups; // clang-format off - results_summary(results_collector &rc, const run_mode &mode, + results_summary(results_collector &&rc, const run_mode &mode, const file_info &info) : n_reads{rc.n_reads}, max_read_len{rc.max_read_len}, diff --git a/src/run_mode.cpp b/src/run_mode.cpp index 2337300..71664f3 100644 --- a/src/run_mode.cpp +++ b/src/run_mode.cpp @@ -83,7 +83,9 @@ 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; - // 'groups' not set in config file + // 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_original_dups_ == 0) do_original_dups_ = do_original_dups_default; // clang-format on } diff --git a/src/run_mode.hpp b/src/run_mode.hpp index b07d984..55d282b 100644 --- a/src/run_mode.hpp +++ b/src/run_mode.hpp @@ -19,9 +19,10 @@ class run_mode { } // clang-format off - // ADS: do_groups is not set in config file + // ADS: these first params are not set in config file [[nodiscard]] auto do_groups() const -> bool { return do_groups_ == 1; } [[nodiscard]] auto do_bisulfite() const -> bool { return do_bisulfite_ == 1; } + [[nodiscard]] auto do_original_dups() const -> bool { return do_original_dups_ == 1; } // [[nodiscard]] auto do_adap() const -> bool { return do_adap_ == 1; } [[nodiscard]] auto do_dups() const -> bool { return do_dups_ == 1; } @@ -43,6 +44,7 @@ class run_mode { // clang-format off auto set_do_groups(const int x) { if (x) do_groups_ = x; } auto set_do_bisulfite(const int x) { if (x) do_bisulfite_ = x; } + auto set_do_original_dups(const int x) { if (x) do_original_dups_ = x; } // auto set_do_adap(const int x) { if (x) do_adap_ = x; } auto set_do_dups(const int x) { if (x) do_dups_ = x; } @@ -65,8 +67,10 @@ class run_mode { private: // ADS: 1 is yes; -1 is no; 0 is not assigned - static constexpr auto do_groups_default = -1; // OFF - static constexpr auto do_bisulfite_default = -1; // OFF + // first settings are not in config file + static constexpr auto do_groups_default = -1; // OFF + static constexpr auto do_bisulfite_default = -1; // OFF + static constexpr auto do_original_dups_default = -1; // OFF static constexpr auto do_adap_default = 1; // affects processing static constexpr auto do_dups_default = 1; // affects processing @@ -80,8 +84,10 @@ class run_mode { static constexpr auto do_sequence_default = 1; static constexpr auto do_tiles_default = 1; // affects processing + // not in config file int do_groups_{}; int do_bisulfite_{}; + int do_original_dups_{}; int do_adap_{}; int do_dups_{};