From 689769ed6a9abfe7b3ecfe612eef140b3b1a81cd Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Wed, 14 Oct 2026 15:42:38 -0700 Subject: [PATCH 01/40] src/get_binary_dir.cpp: added a warning message in case the config dir is not found (e.g., binary was moved relative to original install dir) --- src/get_binary_dir.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/get_binary_dir.cpp b/src/get_binary_dir.cpp index 3e47bb3..4b33761 100644 --- a/src/get_binary_dir.cpp +++ b/src/get_binary_dir.cpp @@ -61,7 +61,8 @@ get_share_dir() -> std::string { if (binary_dir.filename() == bin_name) binary_dir = binary_dir.parent_path(); const auto share_path = binary_dir / DATADIR / PROJECT_NAME; - return is_directory(share_path) ? share_path.string() : std::string{}; + return is_directory(share_path) ? share_path.string() : + std::string{"**config directory not found; reinstall if you need it**"}; } } // namespace falco From 83b7c7c0319e162777d17ccc2850a183d88a5271 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Wed, 14 Oct 2026 15:44:53 -0700 Subject: [PATCH 02/40] src/original_duplicates.cpp: removing unneeded maybe_unused --- src/original_duplicates.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/original_duplicates.cpp b/src/original_duplicates.cpp index 46585c6..d0b960b 100644 --- a/src/original_duplicates.cpp +++ b/src/original_duplicates.cpp @@ -23,8 +23,8 @@ [[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 { + 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 6e2478e13a791a7cf58e75a6737d5d485f9fc7c6 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Wed, 14 Oct 2026 16:11:28 -0700 Subject: [PATCH 03/40] src/falco.cpp: adding a note to the detailed help warning users about non ACGTN nucs in modified input files --- src/falco.cpp | 29 +++++++++++++++++------------ 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/src/falco.cpp b/src/falco.cpp index 41627b5..7219ca4 100644 --- a/src/falco.cpp +++ b/src/falco.cpp @@ -3,7 +3,11 @@ // clang-format off static constexpr auto about = R"(Falco v{})"; static constexpr auto description = - R"(Note: Always use bgzip when compressing to gz format. It comes with samtools. + R"(NOTES: + +- Prefer bgzip over gzip when compressing to gz format. It comes with samtools. +- Falco assumes nucleotides symbols ACGT and N. If you modify the input data + with other symbols, falco results might be inaccurate. EXAMPLES: @@ -29,12 +33,11 @@ Use configuration settings from 'my_limits.txt': Generate a file with duplication info for preseq analysis: $ falco --preseq -o results SRX081761_1.fq.gz -Output files include: -results/SRX081761_1/preseq_hist.txt +Output files will include: results/SRX081761_1/preseq_hist.txt Take input from a pipe (accepts FASTQ or SAM): $ samtools view SRX081761_1.bam | falco --stdin sam -o results SRX081761 -Output files will be created in "results/SRX081761". Reading from stdin +Output files will be created in 'results/SRX081761'. Reading from stdin disables tile analysis. To re-enable it, specify either 4 or 6 to indicate position of the tile id in read names. Example: "--stdin fq:4" @@ -264,12 +267,13 @@ make_outdirs(const auto &ins, const auto &outdir, const bool keep_extn = false) // ADS: in functions below, the input_record_multiplier is the relative space // for a record: ideally for FASTQ it would be 2x the (max) size of the read -// name, and 2x the (max) sequence length. For BAM it would be 0.5x the sequence -// length, the name length and the length of the quality score sequence, plus -// all other BAM fields. The buffer_multiplier is how much total buffer space is -// needed for a given record, which is 3 for BAM and BGZF because one for -// decompression, and then two inflated buffers that are swapped so once can be -// analyzed while the other is the target of other inflation work. +// name, and 2x the (max) sequence length. For BAM it would be 0.5x the +// sequence length, the name length and the length of the quality score +// sequence, plus all other BAM fields. The buffer_multiplier is how much +// total buffer space is needed for a given record, which is 3 for BAM and +// BGZF because one for decompression, and then two inflated buffers that are +// swapped so once can be analyzed while the other is the target of other +// inflation work. [[nodiscard]] static auto get_max_read_length(const auto buffer_size) { @@ -336,6 +340,7 @@ main(int argc, char *argv[]) { struct FormatWithoutFlagDefaults : public CLI::Formatter { FormatWithoutFlagDefaults() : Formatter() { CLI::FormatterBase::enable_default_flag_values_ = false; + CLI::FormatterBase::enable_footer_formatting_ = false; } }; @@ -350,8 +355,8 @@ main(int argc, char *argv[]) { }; CLI::App app{std::format(about, VERSION)}; - const auto fmt = std::make_shared(); - app.formatter(fmt); + const auto app_fmt = std::make_shared(); + app.formatter(app_fmt); argv = app.ensure_utf8(argv); app.usage( From 522c9d0a6dac16ad08705aacf3e9a5b29f410a78 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Wed, 21 Oct 2026 16:54:01 -0700 Subject: [PATCH 04/40] src/run_mode.hpp: adding option for robustness to iupac degen symbols --- src/run_mode.hpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/run_mode.hpp b/src/run_mode.hpp index 4f3de1f..fbbda2c 100644 --- a/src/run_mode.hpp +++ b/src/run_mode.hpp @@ -24,6 +24,7 @@ class run_mode { [[nodiscard]] auto do_original_dups() const -> bool { return do_original_dups_ == 1; } [[nodiscard]] auto do_original_gc() const -> bool { return do_original_gc_ == 1; } [[nodiscard]] auto do_quals() const -> bool { return do_quals_ == 1; } + [[nodiscard]] auto do_iupac() const -> bool { return do_iupac_ == 1; } // [[nodiscard]] auto do_adap() const -> bool { return do_adap_ == 1; } [[nodiscard]] auto do_dups() const -> bool { return do_dups_ == 1; } @@ -50,6 +51,7 @@ class run_mode { auto set_do_original_dups(const int x) { if (x) do_original_dups_ = x; } auto set_do_original_gc(const int x) { if (x) do_original_gc_ = x; } auto set_do_quals(const int x) { if (x) do_quals_ = x; } + auto set_do_iupac(const int x) { if (x) do_iupac_ = 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; } @@ -80,6 +82,7 @@ class run_mode { static constexpr auto do_original_dups_default = -1; // OFF static constexpr auto do_original_dups_gc = -1; // OFF static constexpr auto do_quals_default = 1; // ON + static constexpr auto do_iupac_default = -1; // OFF static constexpr auto do_adap_default = 1; // affects processing static constexpr auto do_dups_default = 1; // affects processing @@ -101,6 +104,7 @@ class run_mode { int do_original_dups_{}; int do_original_gc_{}; int do_quals_{}; + int do_iupac_{}; int do_adap_{}; int do_dups_{}; From 059b46b91ab727b595613acd43c2d6ffc7a057db Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Wed, 21 Oct 2026 16:55:44 -0700 Subject: [PATCH 05/40] src/results_summary.hpp: forgot to turn clang-format back on --- src/results_summary.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/results_summary.hpp b/src/results_summary.hpp index 2694f2b..98c1391 100644 --- a/src/results_summary.hpp +++ b/src/results_summary.hpp @@ -66,6 +66,7 @@ struct results_summary { info{info} { initialize(); } + // clang-format on auto apply_groups() -> void; From 9776897dea7b75a457a509a08e8b520ef42e4952 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Wed, 21 Oct 2026 16:56:12 -0700 Subject: [PATCH 06/40] src/falco_utils.hpp: adding basic functions for anlayzing sequences that might have iupac symbols --- src/falco_utils.hpp | 71 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 71 insertions(+) diff --git a/src/falco_utils.hpp b/src/falco_utils.hpp index aaf934d..90a94ed 100644 --- a/src/falco_utils.hpp +++ b/src/falco_utils.hpp @@ -148,6 +148,49 @@ is_gc(const char c) { return (c >> 1) & 1; } +/// For robustness to IUPAC degenerate nucleutide symbols in the input + +[[nodiscard]] inline constexpr auto +encode_iupac(const char c) { + static constexpr char mask = 31; // (32 - 1) + static constexpr char conversion[] = { + 3, 0, 3, 1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, // 16 + 3, 3, 3, 3, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, // 32 + }; + // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-constant-array-index) + return conversion[c & mask]; +} + +[[nodiscard]] inline constexpr auto +is_gc_iupac(const char c) { + static constexpr char mask = 31; + static constexpr char C_position = 3; + static constexpr char G_position = 7; + const auto low_bits = (c & mask); + return low_bits == C_position || low_bits == G_position; +} + +[[nodiscard]] inline constexpr auto +is_N_iupac(const char c) { + static constexpr auto mask = 31; + static constexpr auto N_encoding = 3; + static constexpr auto G_position = 7; + return encode_iupac(c) == N_encoding && (c & mask) != G_position; +} + +/// ADS: alternate version of is_N_iupac + +// [[nodiscard]] inline constexpr auto +// is_N_iupac(const char c) -> std::uint8_t { +// static constexpr char mask = 31; // (32 - 1) +// static constexpr std::uint8_t value_is_N[] = { +// 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, // 16 +// 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 32 +// }; +// // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-constant-array-index) +// return value_is_N[c & mask]; +// } + [[nodiscard]] inline constexpr auto encode_nibble(const char c) { static constexpr auto four_bit_mask = 15; // ADS: 15 is to keep 4 bits @@ -257,6 +300,15 @@ count_nucs(auto seq_itr, ++(*out_itr++)[encode(*seq_itr++)]; } +inline constexpr auto +count_nucs_iupac(auto seq_itr, + const auto seq_end, + auto &tab) { // cppcheck-suppress constParameterReference + auto out_itr = std::begin(tab); + while (seq_itr != seq_end) + ++(*out_itr++)[encode_iupac(*seq_itr++)]; +} + inline constexpr auto count_ns(auto seq_itr, const auto seq_end, @@ -266,6 +318,15 @@ count_ns(auto seq_itr, *out_itr++ += (*seq_itr++ == 'N'); } +inline constexpr auto +count_ns_iupac(auto seq_itr, + const auto seq_end, + auto &tab) { // cppcheck-suppress constParameterReference + auto out_itr = std::begin(tab); + while (seq_itr != seq_end) + *out_itr++ += is_N_iupac(*seq_itr++); +} + [[nodiscard]] inline auto count_gc(auto seq_itr, const auto seq_end) { auto gc = 0; @@ -276,6 +337,16 @@ count_gc(auto seq_itr, const auto seq_end) { return gc; } +[[nodiscard]] inline auto +count_gc_iupac(auto seq_itr, const auto seq_end) { + auto gc = 0; + while (seq_itr != seq_end) { + gc += is_gc_iupac(*seq_itr); + ++seq_itr; + } + return gc; +} + [[nodiscard]] inline auto tabular_dot(const auto &a) { typename std::remove_cvref_t::value_type total{}; From fe21db3e36fbee0d3c1b12ee4c8bd3798e29ab7a Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Wed, 21 Oct 2026 17:00:27 -0700 Subject: [PATCH 07/40] src/falco.cpp: adding args for robustness to iupac symbols --- src/falco.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/falco.cpp b/src/falco.cpp index 7219ca4..83e7e65 100644 --- a/src/falco.cpp +++ b/src/falco.cpp @@ -324,6 +324,7 @@ main(int argc, char *argv[]) { int do_original_dups{}; int do_original_gc{}; int do_quals{}; + int do_iupac{}; std::uint32_t n_threads{1}; std::uint32_t max_read_length{}; @@ -450,6 +451,9 @@ main(int argc, char *argv[]) { app.add_flag("--orig-gc", do_original_gc, "Use original GC smoothing mode (but no truncation)") ->option_text(" "); + app.add_flag("--iupac", do_iupac, + "For data modified to include non-canonical bases") + ->option_text(" "); app.add_flag("--groups", do_groups, "Group base positions in output (see help)") ->option_text(" "); @@ -508,6 +512,7 @@ main(int argc, char *argv[]) { mode.set_do_original_dups(do_original_dups); mode.set_do_quals(do_quals); mode.set_do_original_gc(do_original_gc); + mode.set_do_iupac(do_iupac); mode.set_unassigned(); const auto outdirs = make_outdirs(infiles, outdir, mode.do_stdin()); From 8c54686027e15ac94ffec39414de9f83b93e3537 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Wed, 21 Oct 2026 17:01:01 -0700 Subject: [PATCH 08/40] src/results_collector.hpp: adding iupac mode via template --- src/results_collector.hpp | 33 +++++++++++++++++++++------------ 1 file changed, 21 insertions(+), 12 deletions(-) diff --git a/src/results_collector.hpp b/src/results_collector.hpp index 5e521fc..366b65b 100644 --- a/src/results_collector.hpp +++ b/src/results_collector.hpp @@ -84,33 +84,37 @@ struct alignas(assumed_page_size) results_collector { base_counts[i][unknown_base_index] -= n_counts[i]; } + template auto process_reads_bam(bamrec::pos_t cursor, const bamrec::pos_t lim) -> void { bamrec rec{}; while (cursor < lim && get_next(cursor, lim, rec)) { - process_one_read(rec); + process_one_read(rec); ++n_reads; } } + template auto process_reads_sam(samrec::pos_t cursor, const samrec::pos_t lim) -> void { samrec rec{}; while (cursor < lim && get_next(cursor, lim, rec)) { - process_one_read(rec); + process_one_read(rec); ++n_reads; } } + template auto process_reads_fq(fqrec::pos_t cursor, const fqrec::pos_t lim) -> void { fqrec rec{}; while (cursor < lim && (rec = get_next(cursor, lim))) { - process_one_read(rec); + process_one_read(rec); ++n_reads; } } + template auto process_one_read(const auto &rec) -> void { // NOLINTBEGIN(cppcoreguidelines-pro-bounds-constant-array-index) @@ -137,7 +141,10 @@ struct alignas(assumed_page_size) results_collector { assert(static_cast(get_arr_idx(read_len)) < std::size(gc_content)); ++gc_content[get_arr_idx(read_len)][get_gc_idx(gc, read_len)]; - count_ns(seq_itr, seq_end, n_counts); + if constexpr (iupac) + count_ns_iupac(seq_itr, seq_end, n_counts); + else + 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]; if (do_original_dups) @@ -182,8 +189,7 @@ struct alignas(assumed_page_size) results_collector { }; inline auto -accumulate_results(std::vector> &r, - const std::int32_t file_id) { +accumulate_results(auto &r, const std::int32_t file_id) { auto j = 1u; while (j < std::size(r)) { std::vector workers; @@ -196,19 +202,22 @@ accumulate_results(std::vector> &r, } } +template inline auto -process_reads(results_collector &results, const fq_task_t &task) { - results.process_reads_fq(task.beg, task.end); +process_reads(auto &results, const fq_task_t &task) { + results.template process_reads_fq(task.beg, task.end); } +template inline auto -process_reads(results_collector &results, const bam_task_t &task) { - results.process_reads_bam(task.beg, task.end); +process_reads(auto &results, const bam_task_t &task) { + results.template process_reads_bam(task.beg, task.end); } +template inline auto -process_reads(results_collector &results, const sam_task_t &task) { - results.process_reads_sam(task.beg, task.end); +process_reads(auto &results, const sam_task_t &task) { + results.template process_reads_sam(task.beg, task.end); } #endif // SRC_RESULTS_COLLECTOR_HPP_ From 1548098c68d8f54799be05a70aaf1abe03f3b215 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Wed, 21 Oct 2026 17:02:20 -0700 Subject: [PATCH 09/40] src/falco_analyzer.cpp: wrapping the analyze implementation for switching to iupac robustness mode --- src/falco_analyzer.cpp | 31 +++++++++++++++++++++++-------- 1 file changed, 23 insertions(+), 8 deletions(-) diff --git a/src/falco_analyzer.cpp b/src/falco_analyzer.cpp index 1b9bd88..ff0d349 100644 --- a/src/falco_analyzer.cpp +++ b/src/falco_analyzer.cpp @@ -28,12 +28,14 @@ class run_mode; +template [[nodiscard]] auto -analyze(const std::uint32_t n_threads, - const run_mode &mode, - std::vector &infos, - std::vector reads_files, - std::vector dups_init) -> std::vector { +analyze_impl(const std::uint32_t n_threads, + const run_mode &mode, + std::vector &infos, + std::vector &&reads_files, + std::vector &&dups_init) + -> std::vector { assert(std::size(reads_files) == std::size(infos)); const auto n_files = static_cast(std::size(infos)); if (dups_init.empty()) @@ -64,11 +66,11 @@ analyze(const std::uint32_t n_threads, return; auto [file_id, task] = tq.pop_and_release_lock(tq_lock); if (std::holds_alternative(task)) - process_reads(res[file_id], std::get(task)); + process_reads(res[file_id], std::get(task)); else if (std::holds_alternative(task)) - process_reads(res[file_id], std::get(task)); + process_reads(res[file_id], std::get(task)); else if (std::holds_alternative(task)) - process_reads(res[file_id], std::get(task)); + process_reads(res[file_id], std::get(task)); else if (std::holds_alternative(task)) decompress(std::get(task)); else { // monostate means read more data @@ -122,3 +124,16 @@ analyze(const std::uint32_t n_threads, return finalized; } + +[[nodiscard]] auto +analyze(const std::uint32_t n_threads, + const run_mode &mode, + std::vector &infos, + std::vector reads_files, + std::vector dups_init) -> std::vector { + return mode.do_iupac() + ? analyze_impl(n_threads, mode, infos, std::move(reads_files), + std::move(dups_init)) + : analyze_impl(n_threads, mode, infos, std::move(reads_files), + std::move(dups_init)); +} From d0404f41ba6816c4f72b888f949e60ce6dbf2215 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Wed, 21 Oct 2026 18:37:46 -0700 Subject: [PATCH 10/40] cmake/tests.cmake: adding a test for iupac robustness --- cmake/tests.cmake | 1 + 1 file changed, 1 insertion(+) diff --git a/cmake/tests.cmake b/cmake/tests.cmake index 3937eeb..aac54a0 100644 --- a/cmake/tests.cmake +++ b/cmake/tests.cmake @@ -37,3 +37,4 @@ add_test(NAME "Preseq output" COMMAND bash test_scripts/preseq.sh) add_test(NAME "FASTQ from stdin" COMMAND bash test_scripts/fastq_stdin.sh) add_test(NAME "SAM from stdin" COMMAND bash test_scripts/sam_stdin.sh) add_test(NAME "Orig GC content" COMMAND bash test_scripts/fastq_orig_gc.sh) +add_test(NAME "IUPAC symbols" COMMAND bash test_scripts/fastq_iupac.sh) From 896db77824cb8bbbd3df5ce27a162d8be65ca7d6 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Wed, 21 Oct 2026 18:38:05 -0700 Subject: [PATCH 11/40] data/test_data/md5sum.txt: adding hash for iupac test --- data/test_data/md5sum.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/data/test_data/md5sum.txt b/data/test_data/md5sum.txt index ff3aa1f..0327217 100644 --- a/data/test_data/md5sum.txt +++ b/data/test_data/md5sum.txt @@ -16,3 +16,4 @@ e786bd4d9e6065904ab7b45021e7190b groups_out/bam_1/fastqc_data.txt 22573e9861f6c2fa3951712a73e4db4c sam_out/sam_1/fastqc_data.txt bcea1e9f86f0917496b14ae81b6ecbc5 sam_stdin_out/sam_stdin/fastqc_data.txt 088c04f706f340f058b663d6e0209dca fastq_orig_gc_out/fastq_bgzip_1/fastqc_data.txt +7da3b7b3cd98c53a1f7939151729b988 fastq_iupac_out/fastq_1_iupac/fastqc_data.txt From 5967c1dab5c2b44197340b28e8ee6200fd7167bd Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Wed, 21 Oct 2026 18:38:52 -0700 Subject: [PATCH 12/40] data/test_scripts/fastq_iupac.sh: adding test for robustness to iupac degen nucs --- data/test_scripts/fastq_iupac.sh | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) create mode 100644 data/test_scripts/fastq_iupac.sh diff --git a/data/test_scripts/fastq_iupac.sh b/data/test_scripts/fastq_iupac.sh new file mode 100644 index 0000000..0949362 --- /dev/null +++ b/data/test_scripts/fastq_iupac.sh @@ -0,0 +1,23 @@ +#!/usr/bin/env bash +# SPDX-License-Identifier: MIT + +prog=./falco +infile1=test_data/fastq_bgzip_1.fq.gz +infile2=fastq_1_iupac.fq +outdir=fastq_iupac_out +if [[ -e "${infile1}" ]]; then + mkdir -p ${outdir} + # ADS: on the first sequence line, modify all C to K + gunzip -c ${infile1}| sed -e "2s/C/K/g" > ${outdir}/${infile2} + ${prog} --iupac -o ${outdir} ${outdir}/${infile2} + 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 From cdfde1a1608a7671e2908dc4e670c1c7f96bde55 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Wed, 21 Oct 2026 19:25:51 -0700 Subject: [PATCH 13/40] src/run_mode.cpp: added more info related to the do_iupac run mode --- src/run_mode.cpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/run_mode.cpp b/src/run_mode.cpp index 499967d..c8babc9 100644 --- a/src/run_mode.cpp +++ b/src/run_mode.cpp @@ -45,6 +45,11 @@ // do_preseq (make another output file for input to preseq) // do_quals (controls any quality score analysis in case no qual scores present) // +// ADS (2026-10-21 edit): added do_iupac as a possibly temporary way for users +// to allow degenerate IUPAC symbols in their data. +// +// do_iupac (ensures robustness to IUPAC degenerate nucleotide symbols) +// // ADS (2026-09-08 edit): added do_stdin because otherwise it's a floating // almost global variable. // @@ -99,6 +104,7 @@ run_mode::set_unassigned() -> void { if (do_preseq_ == 0) do_preseq_ = do_preseq_default; if (do_original_dups_ == 0) do_original_dups_ = do_original_dups_default; if (do_quals_ == 0) do_quals_ = do_quals_default; + if (do_iupac_ == 0) do_iupac_ = do_iupac_default; // if (do_adap_ == 0) do_adap_ = do_adap_default; if (do_dups_ == 0) do_dups_ = do_dups_default; @@ -135,6 +141,7 @@ run_mode::string_verbose() const -> std::string { {"sequence", "sequence composition"}, {"sequence_length", "sequence length"}, {"tile", "per-tile quality"}, + // {"iupac", "ignore IUPAC degen nucs"}, }; static const auto fmt_mode = [&](const std::string &label, const auto the_mode) { @@ -155,6 +162,7 @@ run_mode::string_verbose() const -> std::string { lines.push_back(fmt_mode("quality_sequence", do_qual_seq())); lines.push_back(fmt_mode("sequence", do_sequence())); lines.push_back(fmt_mode("tile", do_tiles())); + // lines.push_back(fmt_mode("iupac", do_iupac())); const auto w = spacer + std::ranges::max(std::views::transform( std::views::values(descriptions), [](const auto &x) { return std::size(x); })); From 096303a5475f8872946c7f1bf8cac1c89c77fee4 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Wed, 21 Oct 2026 19:26:21 -0700 Subject: [PATCH 14/40] src/falco.cpp: some cleanup and better info --- src/falco.cpp | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/src/falco.cpp b/src/falco.cpp index 83e7e65..3e05ffb 100644 --- a/src/falco.cpp +++ b/src/falco.cpp @@ -5,35 +5,33 @@ static constexpr auto about = R"(Falco v{})"; static constexpr auto description = R"(NOTES: -- Prefer bgzip over gzip when compressing to gz format. It comes with samtools. +- Prefer bgzip over gzip when compressing to gz format. Get it with + these packages: tabix (apt), bgzip (pip and dnf) or htslib (brew and conda). - Falco assumes nucleotides symbols ACGT and N. If you modify the input data - with other symbols, falco results might be inaccurate. + with other symbols, use --iupac to ignore non-ACGTN symbols. EXAMPLES: -Use all defaults making output files in a dir named 'results/SRX081761_1': +Use all defaults settings: $ falco -o results SRX081761_1.fastq Output files will be: results/SRX081761_1/summary.txt results/SRX081761_1/fastqc_data.txt results/SRX081761_1/fastqc_results.html -Use 8 cores and analyze all files in 'project' with the fq suffix, and results -for each input file in their own subdirectory of 'results': +Multiple cores for multiple input files: $ falco -t 8 -o results project/*.fq - -Use 64 cores on multiple BAM files in the current directory: $ falco -t 64 -o results sample*.bam -Skip adapter content and tile-specific analysis: +No adapter or tile-specific analysis: $ falco --no-adap --no-tiles -o results SRX081761.bam -Use configuration settings from 'my_limits.txt': +Custom configuration file 'my_limits.txt': $ falco --config my_limits.txt -o results project_sample2.fq.gz -Generate a file with duplication info for preseq analysis: +Generate file for preseq analysis: $ falco --preseq -o results SRX081761_1.fq.gz -Output files will include: results/SRX081761_1/preseq_hist.txt +Output will include: results/SRX081761_1/preseq_hist.txt Take input from a pipe (accepts FASTQ or SAM): $ samtools view SRX081761_1.bam | falco --stdin sam -o results SRX081761 @@ -385,7 +383,7 @@ main(int argc, char *argv[]) { ->check(CLI::ExistingFile, "file_check"); app.add_option("-o,--output", outdir, "Output directory (required)") ->required() - ->option_text("DIR"); + ->option_text("OUTDIR"); app.add_option("-t,--threads", n_threads, std::format("Threads for analysis (this machine supports: {})", std::thread::hardware_concurrency())) From 6075942c0220afd16208366e9129c1efe9e1e2d3 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Wed, 21 Oct 2026 19:41:14 -0700 Subject: [PATCH 15/40] src/falco_analyzer.cpp: static analysis fixes --- src/falco_analyzer.cpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/falco_analyzer.cpp b/src/falco_analyzer.cpp index ff0d349..bb71972 100644 --- a/src/falco_analyzer.cpp +++ b/src/falco_analyzer.cpp @@ -10,6 +10,7 @@ #include "fqrec.hpp" #include "reads_file.hpp" #include "results_collector.hpp" +#include "run_mode.hpp" #include "samrec.hpp" #include "task_queue.hpp" @@ -26,15 +27,13 @@ #include #include -class run_mode; - template [[nodiscard]] auto analyze_impl(const std::uint32_t n_threads, const run_mode &mode, std::vector &infos, - std::vector &&reads_files, - std::vector &&dups_init) + std::vector reads_files, + std::vector dups_init) -> std::vector { assert(std::size(reads_files) == std::size(infos)); const auto n_files = static_cast(std::size(infos)); From 104a8d5fd662b69680fb1dd4c141d3651a2ea003 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 25 Oct 2026 13:59:36 -0700 Subject: [PATCH 16/40] src/bamrec.cpp: modifying how nucleotides are assigned from BAM codes for robustness to IUPAC degen nucs --- src/bamrec.cpp | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/src/bamrec.cpp b/src/bamrec.cpp index be29899..8a7f164 100644 --- a/src/bamrec.cpp +++ b/src/bamrec.cpp @@ -37,23 +37,20 @@ bamrec::to_string() const -> std::string { template static inline constexpr output_itr_t assign_sequence_revcomp(bidir_itr_t first, auto last, output_itr_t d_first) { - static constexpr std::span seq_nt16_str = "=ACMGRSVTWYHKDBN"; - constexpr auto complem = [](const auto x) { - return "TNGNNNCNNNNNNNNNNNNA"[x - 'A']; - }; + static constexpr std::span seq_nt16_str_complem = "NTGNCNNNANNNNNNN"; constexpr auto bam_seqi = [](const auto s, const auto i) -> int { constexpr auto low_nibble_on = 0xf; return s[i >> 1] >> ((~i & 1) << 2) & low_nibble_on; }; for (auto j = last; j != 0; ++d_first) - *d_first = complem(seq_nt16_str[bam_seqi(first, --j)]); + *d_first = seq_nt16_str_complem[bam_seqi(first, --j)]; return d_first; } template static inline constexpr output_itr_t assign_sequence(bidir_itr_t first, auto last, output_itr_t d_first) { - static constexpr std::span seq_nt16_str = "=ACMGRSVTWYHKDBN"; + static constexpr std::span seq_nt16_str = "NACNGNNNTNNNNNNN"; constexpr auto bam_seqi = [](const auto s, const auto i) -> int { constexpr auto low_nibble_on = 0xf; return s[i >> 1] >> ((~i & 1) << 2) & low_nibble_on; From d8d28ee2a44d94f606e7ee00e5b1523b8a060b2e Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 25 Oct 2026 13:59:51 -0700 Subject: [PATCH 17/40] src/samrec.cpp: modifying how nucleotides are assigned from BAM codes for robustness to IUPAC degen nucs --- src/samrec.cpp | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/samrec.cpp b/src/samrec.cpp index 98b4bb4..a79ac55 100644 --- a/src/samrec.cpp +++ b/src/samrec.cpp @@ -39,8 +39,11 @@ samrec::get_next(samrec::pos_t &cursor, const samrec::pos_t end_itr, samrec &rec) -> bool { static constexpr auto n_fields_to_skip = 7; - constexpr auto complem = [](const auto x) { - return "TNGNNNCNNNNNNNNNNNNA"[x - 'A']; + const auto iupac = [](const auto c) -> char { + return "NANCNNNGNNNNNNNNNNNNTNNNNNNNNNNN"[c & 31]; // NOLINT + }; + const auto iupac_cmpl = [](const auto c) -> char { + return "NTNGNNNCNNNNNNNNNNNNANNNNNNNNNNN"[c & 31]; // NOLINT }; static constexpr auto bam_is_rev = [](const auto flag) -> bool { static constexpr auto BAM_FREVERSE = 16; @@ -107,14 +110,14 @@ samrec::get_next(samrec::pos_t &cursor, if (bam_is_rev(flag)) { std::reverse_copy(seq_itr, seq_itr + seq_len, out_seq_itr); - std::transform(out_seq_itr, out_seq_itr + seq_len, out_seq_itr, complem); + std::transform(out_seq_itr, out_seq_itr + seq_len, out_seq_itr, iupac_cmpl); if (*qual_itr != qual_missing_code) std::reverse_copy(qual_itr, qual_itr + seq_len, out_qual_itr); else *out_qual_itr = qual_missing_code; } else { - std::copy(seq_itr, seq_itr + seq_len, out_seq_itr); + std::transform(seq_itr, seq_itr + seq_len, out_seq_itr, iupac); if (*qual_itr != qual_missing_code) std::copy(qual_itr, qual_itr + seq_len, out_qual_itr); else From 2de92b5e45cb334d7f43e84b387ae0c1cd43983f Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 25 Oct 2026 14:00:13 -0700 Subject: [PATCH 18/40] src/run_mode.hcpp: undoing the changes related to iupac mode --- src/run_mode.cpp | 8 -------- src/run_mode.hpp | 4 ---- 2 files changed, 12 deletions(-) diff --git a/src/run_mode.cpp b/src/run_mode.cpp index c8babc9..499967d 100644 --- a/src/run_mode.cpp +++ b/src/run_mode.cpp @@ -45,11 +45,6 @@ // do_preseq (make another output file for input to preseq) // do_quals (controls any quality score analysis in case no qual scores present) // -// ADS (2026-10-21 edit): added do_iupac as a possibly temporary way for users -// to allow degenerate IUPAC symbols in their data. -// -// do_iupac (ensures robustness to IUPAC degenerate nucleotide symbols) -// // ADS (2026-09-08 edit): added do_stdin because otherwise it's a floating // almost global variable. // @@ -104,7 +99,6 @@ run_mode::set_unassigned() -> void { if (do_preseq_ == 0) do_preseq_ = do_preseq_default; if (do_original_dups_ == 0) do_original_dups_ = do_original_dups_default; if (do_quals_ == 0) do_quals_ = do_quals_default; - if (do_iupac_ == 0) do_iupac_ = do_iupac_default; // if (do_adap_ == 0) do_adap_ = do_adap_default; if (do_dups_ == 0) do_dups_ = do_dups_default; @@ -141,7 +135,6 @@ run_mode::string_verbose() const -> std::string { {"sequence", "sequence composition"}, {"sequence_length", "sequence length"}, {"tile", "per-tile quality"}, - // {"iupac", "ignore IUPAC degen nucs"}, }; static const auto fmt_mode = [&](const std::string &label, const auto the_mode) { @@ -162,7 +155,6 @@ run_mode::string_verbose() const -> std::string { lines.push_back(fmt_mode("quality_sequence", do_qual_seq())); lines.push_back(fmt_mode("sequence", do_sequence())); lines.push_back(fmt_mode("tile", do_tiles())); - // lines.push_back(fmt_mode("iupac", do_iupac())); const auto w = spacer + std::ranges::max(std::views::transform( std::views::values(descriptions), [](const auto &x) { return std::size(x); })); diff --git a/src/run_mode.hpp b/src/run_mode.hpp index fbbda2c..4f3de1f 100644 --- a/src/run_mode.hpp +++ b/src/run_mode.hpp @@ -24,7 +24,6 @@ class run_mode { [[nodiscard]] auto do_original_dups() const -> bool { return do_original_dups_ == 1; } [[nodiscard]] auto do_original_gc() const -> bool { return do_original_gc_ == 1; } [[nodiscard]] auto do_quals() const -> bool { return do_quals_ == 1; } - [[nodiscard]] auto do_iupac() const -> bool { return do_iupac_ == 1; } // [[nodiscard]] auto do_adap() const -> bool { return do_adap_ == 1; } [[nodiscard]] auto do_dups() const -> bool { return do_dups_ == 1; } @@ -51,7 +50,6 @@ class run_mode { auto set_do_original_dups(const int x) { if (x) do_original_dups_ = x; } auto set_do_original_gc(const int x) { if (x) do_original_gc_ = x; } auto set_do_quals(const int x) { if (x) do_quals_ = x; } - auto set_do_iupac(const int x) { if (x) do_iupac_ = 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; } @@ -82,7 +80,6 @@ class run_mode { static constexpr auto do_original_dups_default = -1; // OFF static constexpr auto do_original_dups_gc = -1; // OFF static constexpr auto do_quals_default = 1; // ON - static constexpr auto do_iupac_default = -1; // OFF static constexpr auto do_adap_default = 1; // affects processing static constexpr auto do_dups_default = 1; // affects processing @@ -104,7 +101,6 @@ class run_mode { int do_original_dups_{}; int do_original_gc_{}; int do_quals_{}; - int do_iupac_{}; int do_adap_{}; int do_dups_{}; From b6a5daef74f80516fa198e5e34e3ed16d933c2c9 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 25 Oct 2026 14:01:25 -0700 Subject: [PATCH 19/40] src/fqrec.hpp: (1) assuming memory pointed to by 'cursor' is writable, and modifying the read sequence part to fix iupac degenerate nucleotides to N --- src/fqrec.hpp | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/src/fqrec.hpp b/src/fqrec.hpp index 6737415..f7c9bf7 100644 --- a/src/fqrec.hpp +++ b/src/fqrec.hpp @@ -59,8 +59,20 @@ struct fq_task_t { beg{std::to_address(beg_arg)}, end{std::to_address(end_arg)} {} }; +struct fq_gz_task_t { + fqrec::pos_t beg{}; + fqrec::pos_t end{}; + fq_gz_task_t(const fqrec::pos_t beg, const fqrec::pos_t end) : + beg{beg}, end{end} {} + fq_gz_task_t(const auto beg_arg, const auto end_arg) : + beg{std::to_address(beg_arg)}, end{std::to_address(end_arg)} {} +}; + [[nodiscard]] inline auto get_next(fqrec::pos_t &cursor, const fqrec::pos_t end_itr) -> fqrec { + const auto iupac = [](auto &c) { + c = "NANCNNNGNNNNNNNNNNNNTNNNNNNNNNNN"[c & 31]; // NOLINT + }; // ADS: need to make sure cursor < end_itr or we will move past const auto n = cursor; auto itr = n + 1; // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic) @@ -75,7 +87,8 @@ get_next(fqrec::pos_t &cursor, const fqrec::pos_t end_itr) -> fqrec { if (itr++ == end_itr) return {}; const auto r = itr; - next_newline(itr); + for (; itr != end_itr && *itr != '\n'; ++itr) + iupac(*itr); if (itr++ == end_itr) return {}; const auto o = itr; From 6218cad87ffd3f24a5e1c165f721bba53238f098 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 25 Oct 2026 14:04:45 -0700 Subject: [PATCH 20/40] src/falco_task.hpp: adding the fq_gz_task_t task type which allows for choice to copy the file's data buffer vs. modifying it directly. Direct modification can't work with memory mapped files using multiple cores because it invalidates all caches and just adding the write mode (so no segfault on mmap data modification) causes many-fold slowdown --- src/falco_task.hpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/falco_task.hpp b/src/falco_task.hpp index d2a5c35..78d01e6 100644 --- a/src/falco_task.hpp +++ b/src/falco_task.hpp @@ -10,7 +10,11 @@ #include -using task_t = - std::variant; +using task_t = std::variant; #endif // SRC_FALCO_TASK_HPP_ From 14f60f0948171526cdeb113ddc9281b3046c2ff6 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 25 Oct 2026 14:06:03 -0700 Subject: [PATCH 21/40] src/results_collector.hpp: mostly reverting templated results_collector but importantly having an overload to take fq_gz_task_t so that the fq_task_t can make a locally allocated buffer that can be modified --- src/results_collector.hpp | 40 +++++++++++++++++++-------------------- 1 file changed, 19 insertions(+), 21 deletions(-) diff --git a/src/results_collector.hpp b/src/results_collector.hpp index 366b65b..38c9133 100644 --- a/src/results_collector.hpp +++ b/src/results_collector.hpp @@ -84,37 +84,33 @@ struct alignas(assumed_page_size) results_collector { base_counts[i][unknown_base_index] -= n_counts[i]; } - template auto process_reads_bam(bamrec::pos_t cursor, const bamrec::pos_t lim) -> void { bamrec rec{}; while (cursor < lim && get_next(cursor, lim, rec)) { - process_one_read(rec); + process_one_read(rec); ++n_reads; } } - template auto process_reads_sam(samrec::pos_t cursor, const samrec::pos_t lim) -> void { samrec rec{}; while (cursor < lim && get_next(cursor, lim, rec)) { - process_one_read(rec); + process_one_read(rec); ++n_reads; } } - template auto process_reads_fq(fqrec::pos_t cursor, const fqrec::pos_t lim) -> void { fqrec rec{}; while (cursor < lim && (rec = get_next(cursor, lim))) { - process_one_read(rec); + process_one_read(rec); ++n_reads; } } - template auto process_one_read(const auto &rec) -> void { // NOLINTBEGIN(cppcoreguidelines-pro-bounds-constant-array-index) @@ -141,10 +137,7 @@ struct alignas(assumed_page_size) results_collector { assert(static_cast(get_arr_idx(read_len)) < std::size(gc_content)); ++gc_content[get_arr_idx(read_len)][get_gc_idx(gc, read_len)]; - if constexpr (iupac) - count_ns_iupac(seq_itr, seq_end, n_counts); - else - count_ns(seq_itr, seq_end, n_counts); + 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]; if (do_original_dups) @@ -189,7 +182,8 @@ struct alignas(assumed_page_size) results_collector { }; inline auto -accumulate_results(auto &r, const std::int32_t file_id) { +accumulate_results(std::vector> &r, + const std::int32_t file_id) { auto j = 1u; while (j < std::size(r)) { std::vector workers; @@ -202,22 +196,26 @@ accumulate_results(auto &r, const std::int32_t file_id) { } } -template inline auto -process_reads(auto &results, const fq_task_t &task) { - results.template process_reads_fq(task.beg, task.end); +process_reads(results_collector &results, const fq_task_t &task) { + std::vector local_data(task.beg, task.end); + results.process_reads_fq(std::to_address(std::begin(local_data)), + std::to_address(std::end(local_data))); +} + +inline auto +process_reads(results_collector &results, const fq_gz_task_t &task) { + results.process_reads_fq(task.beg, task.end); } -template inline auto -process_reads(auto &results, const bam_task_t &task) { - results.template process_reads_bam(task.beg, task.end); +process_reads(results_collector &results, const bam_task_t &task) { + results.process_reads_bam(task.beg, task.end); } -template inline auto -process_reads(auto &results, const sam_task_t &task) { - results.template process_reads_sam(task.beg, task.end); +process_reads(results_collector &results, const sam_task_t &task) { + results.process_reads_sam(task.beg, task.end); } #endif // SRC_RESULTS_COLLECTOR_HPP_ From 536e85f058dae095dbab058ce1c4f768df1b731a Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 25 Oct 2026 14:07:02 -0700 Subject: [PATCH 22/40] src/falco.cpp: reverting iupac options and modifying the comment in detailed help about iupac degen nucs --- src/falco.cpp | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/falco.cpp b/src/falco.cpp index 3e05ffb..198eb4a 100644 --- a/src/falco.cpp +++ b/src/falco.cpp @@ -8,7 +8,7 @@ static constexpr auto description = - Prefer bgzip over gzip when compressing to gz format. Get it with these packages: tabix (apt), bgzip (pip and dnf) or htslib (brew and conda). - Falco assumes nucleotides symbols ACGT and N. If you modify the input data - with other symbols, use --iupac to ignore non-ACGTN symbols. + with other symbols, they are treated as N, but ACGT are case-insensitive. EXAMPLES: @@ -322,7 +322,6 @@ main(int argc, char *argv[]) { int do_original_dups{}; int do_original_gc{}; int do_quals{}; - int do_iupac{}; std::uint32_t n_threads{1}; std::uint32_t max_read_length{}; @@ -449,9 +448,6 @@ main(int argc, char *argv[]) { app.add_flag("--orig-gc", do_original_gc, "Use original GC smoothing mode (but no truncation)") ->option_text(" "); - app.add_flag("--iupac", do_iupac, - "For data modified to include non-canonical bases") - ->option_text(" "); app.add_flag("--groups", do_groups, "Group base positions in output (see help)") ->option_text(" "); @@ -510,7 +506,6 @@ main(int argc, char *argv[]) { mode.set_do_original_dups(do_original_dups); mode.set_do_quals(do_quals); mode.set_do_original_gc(do_original_gc); - mode.set_do_iupac(do_iupac); mode.set_unassigned(); const auto outdirs = make_outdirs(infiles, outdir, mode.do_stdin()); From 630c930b7ca818212b87acdef54be478445c56cf Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 25 Oct 2026 14:07:38 -0700 Subject: [PATCH 23/40] src/falco_utils.hpp: reverting iupac code and adding a comment related to assumptions on functions that take the sequences --- src/falco_utils.hpp | 75 +++------------------------------------------ 1 file changed, 4 insertions(+), 71 deletions(-) diff --git a/src/falco_utils.hpp b/src/falco_utils.hpp index 90a94ed..b19b673 100644 --- a/src/falco_utils.hpp +++ b/src/falco_utils.hpp @@ -20,6 +20,10 @@ #include #include +/// ADS: (2026-10-24) intention currently is to be robust to IUPAC degenerate +/// nucleotide symbols, but current assumption is that this is taken care of +/// when constructing each record type (i.e., FASTQ, SAM, BAM, etc.) + // N (78)10 = (1001110)2 => 11 => 3 (same as G) // A (65)10 = (1000001)2 => 00 => 0 // C (67)10 = (1000011)2 => 01 => 1 @@ -148,49 +152,6 @@ is_gc(const char c) { return (c >> 1) & 1; } -/// For robustness to IUPAC degenerate nucleutide symbols in the input - -[[nodiscard]] inline constexpr auto -encode_iupac(const char c) { - static constexpr char mask = 31; // (32 - 1) - static constexpr char conversion[] = { - 3, 0, 3, 1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, // 16 - 3, 3, 3, 3, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, // 32 - }; - // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-constant-array-index) - return conversion[c & mask]; -} - -[[nodiscard]] inline constexpr auto -is_gc_iupac(const char c) { - static constexpr char mask = 31; - static constexpr char C_position = 3; - static constexpr char G_position = 7; - const auto low_bits = (c & mask); - return low_bits == C_position || low_bits == G_position; -} - -[[nodiscard]] inline constexpr auto -is_N_iupac(const char c) { - static constexpr auto mask = 31; - static constexpr auto N_encoding = 3; - static constexpr auto G_position = 7; - return encode_iupac(c) == N_encoding && (c & mask) != G_position; -} - -/// ADS: alternate version of is_N_iupac - -// [[nodiscard]] inline constexpr auto -// is_N_iupac(const char c) -> std::uint8_t { -// static constexpr char mask = 31; // (32 - 1) -// static constexpr std::uint8_t value_is_N[] = { -// 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, // 16 -// 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 32 -// }; -// // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-constant-array-index) -// return value_is_N[c & mask]; -// } - [[nodiscard]] inline constexpr auto encode_nibble(const char c) { static constexpr auto four_bit_mask = 15; // ADS: 15 is to keep 4 bits @@ -300,15 +261,6 @@ count_nucs(auto seq_itr, ++(*out_itr++)[encode(*seq_itr++)]; } -inline constexpr auto -count_nucs_iupac(auto seq_itr, - const auto seq_end, - auto &tab) { // cppcheck-suppress constParameterReference - auto out_itr = std::begin(tab); - while (seq_itr != seq_end) - ++(*out_itr++)[encode_iupac(*seq_itr++)]; -} - inline constexpr auto count_ns(auto seq_itr, const auto seq_end, @@ -318,15 +270,6 @@ count_ns(auto seq_itr, *out_itr++ += (*seq_itr++ == 'N'); } -inline constexpr auto -count_ns_iupac(auto seq_itr, - const auto seq_end, - auto &tab) { // cppcheck-suppress constParameterReference - auto out_itr = std::begin(tab); - while (seq_itr != seq_end) - *out_itr++ += is_N_iupac(*seq_itr++); -} - [[nodiscard]] inline auto count_gc(auto seq_itr, const auto seq_end) { auto gc = 0; @@ -337,16 +280,6 @@ count_gc(auto seq_itr, const auto seq_end) { return gc; } -[[nodiscard]] inline auto -count_gc_iupac(auto seq_itr, const auto seq_end) { - auto gc = 0; - while (seq_itr != seq_end) { - gc += is_gc_iupac(*seq_itr); - ++seq_itr; - } - return gc; -} - [[nodiscard]] inline auto tabular_dot(const auto &a) { typename std::remove_cvref_t::value_type total{}; From 6bcb1e37ec148aee6888c7176ac76303f6ee5128 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 25 Oct 2026 14:08:46 -0700 Subject: [PATCH 24/40] src/falco_analyzer.cpp: reverting iupac-specific (a template on analyze function implementation) and adding a condition for holding the new type fq_gz_task_t --- src/falco_analyzer.cpp | 36 ++++++++++++------------------------ 1 file changed, 12 insertions(+), 24 deletions(-) diff --git a/src/falco_analyzer.cpp b/src/falco_analyzer.cpp index bb71972..70105c4 100644 --- a/src/falco_analyzer.cpp +++ b/src/falco_analyzer.cpp @@ -10,7 +10,6 @@ #include "fqrec.hpp" #include "reads_file.hpp" #include "results_collector.hpp" -#include "run_mode.hpp" #include "samrec.hpp" #include "task_queue.hpp" @@ -27,14 +26,14 @@ #include #include -template +class run_mode; + [[nodiscard]] auto -analyze_impl(const std::uint32_t n_threads, - const run_mode &mode, - std::vector &infos, - std::vector reads_files, - std::vector dups_init) - -> std::vector { +analyze(const std::uint32_t n_threads, + const run_mode &mode, + std::vector &infos, + std::vector reads_files, + std::vector dups_init) -> std::vector { assert(std::size(reads_files) == std::size(infos)); const auto n_files = static_cast(std::size(infos)); if (dups_init.empty()) @@ -65,11 +64,13 @@ analyze_impl(const std::uint32_t n_threads, return; auto [file_id, task] = tq.pop_and_release_lock(tq_lock); if (std::holds_alternative(task)) - process_reads(res[file_id], std::get(task)); + process_reads(res[file_id], std::get(task)); + else if (std::holds_alternative(task)) + process_reads(res[file_id], std::get(task)); else if (std::holds_alternative(task)) - process_reads(res[file_id], std::get(task)); + process_reads(res[file_id], std::get(task)); else if (std::holds_alternative(task)) - process_reads(res[file_id], std::get(task)); + process_reads(res[file_id], std::get(task)); else if (std::holds_alternative(task)) decompress(std::get(task)); else { // monostate means read more data @@ -123,16 +124,3 @@ analyze_impl(const std::uint32_t n_threads, return finalized; } - -[[nodiscard]] auto -analyze(const std::uint32_t n_threads, - const run_mode &mode, - std::vector &infos, - std::vector reads_files, - std::vector dups_init) -> std::vector { - return mode.do_iupac() - ? analyze_impl(n_threads, mode, infos, std::move(reads_files), - std::move(dups_init)) - : analyze_impl(n_threads, mode, infos, std::move(reads_files), - std::move(dups_init)); -} From b8eef2731427d00793116e5019b4ef16b9bc0eb1 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 25 Oct 2026 14:10:37 -0700 Subject: [PATCH 25/40] src/fastq_gz_file.cpp and src/fastq_bgzf_file.cpp: changing the task type from fq_task_t to fq_gz_task_t which has the same behavior, as fq_task_t was just changed to make a copy of the buffer so it can be modified local to the thread's core --- src/fastq_bgzf_file.cpp | 2 +- src/fastq_gz_file.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/fastq_bgzf_file.cpp b/src/fastq_bgzf_file.cpp index edf6de2..91d30a7 100644 --- a/src/fastq_bgzf_file.cpp +++ b/src/fastq_bgzf_file.cpp @@ -163,7 +163,7 @@ fastq_bgzf_file::get_chunks(std::int64_t n_chunks, } ++n_tasks; // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic) - tq.push(file_id, fq_task_t(data + chunk_beg, data + chunk_end)); + tq.push(file_id, fq_gz_task_t(data + chunk_beg, data + chunk_end)); start_pos = stop_pos; } output_cursor = chunk_end; diff --git a/src/fastq_gz_file.cpp b/src/fastq_gz_file.cpp index e87e104..88bedb1 100644 --- a/src/fastq_gz_file.cpp +++ b/src/fastq_gz_file.cpp @@ -98,7 +98,7 @@ fastq_gz_file::get_chunks(std::int64_t n_chunks, } ++n_tasks; // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic) - tq.push(file_id, fq_task_t(data + chunk_beg, data + chunk_end)); + tq.push(file_id, fq_gz_task_t(data + chunk_beg, data + chunk_end)); start_pos = stop_pos; } cursor = chunk_end; From 3a444a7c1715462c92eb7c1e96dccc593c6af9ca Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 25 Oct 2026 14:13:40 -0700 Subject: [PATCH 26/40] data/test_data/md5sum.txt: updating test hash for iupac test falco no longer has --iupac arg --- data/test_data/md5sum.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/data/test_data/md5sum.txt b/data/test_data/md5sum.txt index 0327217..d227693 100644 --- a/data/test_data/md5sum.txt +++ b/data/test_data/md5sum.txt @@ -16,4 +16,4 @@ e786bd4d9e6065904ab7b45021e7190b groups_out/bam_1/fastqc_data.txt 22573e9861f6c2fa3951712a73e4db4c sam_out/sam_1/fastqc_data.txt bcea1e9f86f0917496b14ae81b6ecbc5 sam_stdin_out/sam_stdin/fastqc_data.txt 088c04f706f340f058b663d6e0209dca fastq_orig_gc_out/fastq_bgzip_1/fastqc_data.txt -7da3b7b3cd98c53a1f7939151729b988 fastq_iupac_out/fastq_1_iupac/fastqc_data.txt +e71cef09fd4b6f4019f8fdb82ee2eaf3 fastq_iupac_out/fastq_1_iupac/fastqc_data.txt From 826de0a8ec09af2049c14fb6ca89c5855f1764cf Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 25 Oct 2026 14:14:07 -0700 Subject: [PATCH 27/40] data/test_scripts/fastq_iupac.sh: updating test script iupac test as falco no longer has --iupac arg --- data/test_scripts/fastq_iupac.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/data/test_scripts/fastq_iupac.sh b/data/test_scripts/fastq_iupac.sh index 0949362..25e5f3f 100644 --- a/data/test_scripts/fastq_iupac.sh +++ b/data/test_scripts/fastq_iupac.sh @@ -9,7 +9,7 @@ if [[ -e "${infile1}" ]]; then mkdir -p ${outdir} # ADS: on the first sequence line, modify all C to K gunzip -c ${infile1}| sed -e "2s/C/K/g" > ${outdir}/${infile2} - ${prog} --iupac -o ${outdir} ${outdir}/${infile2} + ${prog} -o ${outdir} ${outdir}/${infile2} x=$(md5sum --ignore-missing -c test_data/md5sum.txt | \ grep "${outdir}" | \ grep -c "OK$") From d06cb363d422cc54b1df9ceb125ccca0601db982 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 25 Oct 2026 18:52:54 -0700 Subject: [PATCH 28/40] src/fqrec.hpp: adding a size function to fqrec --- src/fqrec.hpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/fqrec.hpp b/src/fqrec.hpp index f7c9bf7..ab91c92 100644 --- a/src/fqrec.hpp +++ b/src/fqrec.hpp @@ -57,6 +57,10 @@ struct fq_task_t { beg{beg}, end{end} {} fq_task_t(const auto beg_arg, const auto end_arg) : beg{std::to_address(beg_arg)}, end{std::to_address(end_arg)} {} + [[nodiscard]] auto + size() const -> std::size_t { + return static_cast(std::max(std::distance(beg, end), 0L)); + } }; struct fq_gz_task_t { From adabede04a7a8b32eceec161072a36eba96f05e2 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 25 Oct 2026 18:56:56 -0700 Subject: [PATCH 29/40] src/results_collector.hpp: modifying process_reads overload for fq_task_t so it accepts a buffer to copy mmap data into, so reallocates only as needed --- src/results_collector.hpp | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/src/results_collector.hpp b/src/results_collector.hpp index 38c9133..10a56f8 100644 --- a/src/results_collector.hpp +++ b/src/results_collector.hpp @@ -197,10 +197,15 @@ accumulate_results(std::vector> &r, } inline auto -process_reads(results_collector &results, const fq_task_t &task) { - std::vector local_data(task.beg, task.end); - results.process_reads_fq(std::to_address(std::begin(local_data)), - std::to_address(std::end(local_data))); +process_reads(results_collector &results, + const fq_task_t &task, + std::vector &fq_buffer) { + const auto task_size = std::size(task); + if (std::size(fq_buffer) < task_size) + fq_buffer.resize(task_size); + std::memcpy(std::data(fq_buffer), task.beg, task_size); + const auto beg = std::to_address(std::begin(fq_buffer)); + results.process_reads_fq(beg, beg + task_size); } inline auto From 4d0381f82c909340525e09536bbec804f4903c27 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 25 Oct 2026 19:01:54 -0700 Subject: [PATCH 30/40] src/falco_analyzer.cpp: adding fq_buffer to hold thread local data for mmap plain fastq file input and avoid reallocating; the process_reads for fq_task_t accepts and uses this buffer --- src/falco_analyzer.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/falco_analyzer.cpp b/src/falco_analyzer.cpp index 70105c4..b13f9b7 100644 --- a/src/falco_analyzer.cpp +++ b/src/falco_analyzer.cpp @@ -55,6 +55,10 @@ analyze(const std::uint32_t n_threads, for (const auto th_id : std::views::iota(0u, n_threads)) workers.emplace_back([&, n_threads, th_id] { try { + /// ADS: the 'fq_buffer' is used to make a thread-local copy of mmap + /// data so it can be modified without invalidating all + /// memories. Declare it here so we have one per thread. + std::vector fq_buffer; auto &res = results[th_id]; for (const auto [file_id, info] : falco::views::enumerate(infos)) res[file_id].init(mode, info, dups_init[file_id]); @@ -64,7 +68,7 @@ analyze(const std::uint32_t n_threads, return; auto [file_id, task] = tq.pop_and_release_lock(tq_lock); if (std::holds_alternative(task)) - process_reads(res[file_id], std::get(task)); + process_reads(res[file_id], std::get(task), fq_buffer); else if (std::holds_alternative(task)) process_reads(res[file_id], std::get(task)); else if (std::holds_alternative(task)) From 97f927e7beb07bfc62cffb9d81c5173972e96155 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Tue, 27 Oct 2026 16:04:56 -0700 Subject: [PATCH 31/40] src/fqrec.hpp: simplifying loop --- src/fqrec.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/fqrec.hpp b/src/fqrec.hpp index ab91c92..7018539 100644 --- a/src/fqrec.hpp +++ b/src/fqrec.hpp @@ -91,8 +91,8 @@ get_next(fqrec::pos_t &cursor, const fqrec::pos_t end_itr) -> fqrec { if (itr++ == end_itr) return {}; const auto r = itr; - for (; itr != end_itr && *itr != '\n'; ++itr) - iupac(*itr); + while (itr != end_itr && *itr != '\n') + iupac(*itr++); if (itr++ == end_itr) return {}; const auto o = itr; From 734a19758bf7622c9de36d5c9a7b42a814663488 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Tue, 27 Oct 2026 16:07:36 -0700 Subject: [PATCH 32/40] src/fastq_file.hcpp: merging the code for mmap fastq with the code for reading fastq using fread with a preprocessor guard to distinguish --- src/fastq_file.cpp | 107 ++++++++++++++++++++++++++++++++++++++++++--- src/fastq_file.hpp | 76 +++++++++++++++++++++++++++++++- 2 files changed, 176 insertions(+), 7 deletions(-) diff --git a/src/fastq_file.cpp b/src/fastq_file.cpp index 2abd3e7..08544c8 100644 --- a/src/fastq_file.cpp +++ b/src/fastq_file.cpp @@ -13,19 +13,24 @@ #include #include #include -#include #include #include #include #include #include -#include #include -#include #include #include +#ifndef USE_MMAP_FOR_PLAIN_FASTQ +#include +#include +#include #include #include +#else // USE_MMAP_FOR_PLAIN_FASTQ +#include +#include +#endif [[nodiscard]] auto estimate_n_reads_fastq(const std::string &filename) @@ -74,6 +79,89 @@ estimate_n_reads_fastq(const std::string &filename) return {static_cast(n_reads_est), readlen_est, filesize}; } +#ifndef USE_MMAP_FOR_PLAIN_FASTQ + +fastq_file::fastq_file(const std::string &filename, + const std::int64_t buf_size) : + buffer(buf_size + min_buf_size), + in(std::fopen(std::data(filename), "r"), &std::fclose) { + if (!in) + std::system_error(std::make_error_code(std::errc(errno)), + "failed to read file"); + cursor = std::begin(buffer); + last = std::begin(buffer); +} + +auto +fastq_file::shift_output_buffer() -> void { + if (cursor == std::cbegin(buffer)) // shifting here does nothing + return; + last = std::copy(cursor, last, std::begin(buffer)); + cursor = std::begin(buffer); +} + +auto +fastq_file::get_chunks(std::int64_t n_chunks, + const std::int32_t file_id, + task_queue &tq, + std::atomic_int32_t &n_tasks) -> void { + static constexpr std::int64_t min_chunk_size = 64 * 1024; + static constexpr auto rec_lines = 4; // FASTQ + assert(n_chunks > 0); + const auto beg_itr = std::begin(buffer); + const auto end_itr = last; + // clang-format off + const auto not_read_start = [&](const auto p) { + // ADS: could get confused if '+' lines have full name info + return *p != '@' || (p > beg_itr && *(p - 1) != '\n') || + (p > beg_itr + 2 && *(p - 2) == '+' && *(p - 3) == '\n'); + }; + const auto forward_to_start = [&](auto p) { + if (p == beg_itr) return p; + while (p < end_itr && not_read_start(p)) ++p; + return p; + }; + const auto rev_to_read_start = [&](auto p) { + while (p != beg_itr && (p == end_itr || not_read_start(p))) --p; + return p; + }; + // clang-format on + const std::int64_t n_bytes_available = std::distance(beg_itr, end_itr); + n_chunks = std::min((n_bytes_available + min_chunk_size - 1) / min_chunk_size, + n_chunks); + auto [chunk_size, remainder] = std::div(n_bytes_available, n_chunks); + assert(n_chunks > 0); + auto start_itr = beg_itr; + auto chunk_end = start_itr; + for (const auto chunk_idx : std::views::iota(0, n_chunks)) { + const auto chunk_beg = forward_to_start(start_itr); + const auto stop_itr = start_itr + chunk_size + (chunk_idx < remainder); + chunk_end = forward_to_start(stop_itr); + if (chunk_idx + 1 == n_chunks) { + const auto prev_start = rev_to_read_start(chunk_end); + const auto n_trailing = std::count(prev_start, end_itr, '\n'); + if (n_trailing < rec_lines) + chunk_end = prev_start; + } + ++n_tasks; + tq.push(file_id, fq_task_t(chunk_beg, chunk_end)); + start_itr = stop_itr; + } + cursor = chunk_end; +} + +auto +fastq_file::load_next() -> void { + const auto n_bytes = std::distance(last, std::end(buffer)); + last += static_cast( + std::fread(std::to_address(last), 1, n_bytes, in.get())); + if (std::ferror(in.get())) + std::system_error(std::make_error_code(std::errc(errno)), + "error reading FASTQ file"); +} + +#else // NOT USE_MMAP_FOR_PLAIN_FASTQ + fastq_file::fastq_file(const std::string &filename, const std::int64_t target_length) : target_length{ @@ -166,10 +254,12 @@ fastq_file::get_chunks(std::int64_t n_chunks, const auto chunk_beg = forward_to_start(start_pos); const auto stop_pos = start_pos + chunk_size + (chunk_idx < remainder); chunk_end = forward_to_start(stop_pos); - if (chunk_idx + 1 == n_chunks) - if (const auto prev = reverse_to_start(chunk_end); - std::count(prev, end_itr, '\n') < rec_lines) + if (chunk_idx + 1 == n_chunks) { + const auto prev = reverse_to_start(chunk_end); + const auto n_trailing = std::count(prev, end_itr, '\n'); + if (n_trailing < rec_lines) chunk_end = prev; + } ++n_tasks; tq.push(file_id, fq_task_t(chunk_beg, chunk_end)); start_pos = stop_pos; @@ -177,6 +267,8 @@ fastq_file::get_chunks(std::int64_t n_chunks, last = chunk_end; } +#endif // USE_MMAP_FOR_PLAIN_FASTQ + auto fastq_file::make_tasks(const std::int64_t n_threads, const std::int32_t file_id, @@ -184,6 +276,9 @@ fastq_file::make_tasks(const std::int64_t n_threads, std::atomic_int32_t &n_tasks) -> void { static constexpr auto n_chunks_per_thread = 8; const auto n_chunks = n_chunks_per_thread * n_threads; +#ifndef USE_MMAP_FOR_PLAIN_FASTQ + shift_output_buffer(); +#endif // USE_MMAP_FOR_PLAIN_FASTQ load_next(); get_chunks(n_chunks, file_id, tq, n_tasks); } diff --git a/src/fastq_file.hpp b/src/fastq_file.hpp index 19d5ec1..9d13d18 100644 --- a/src/fastq_file.hpp +++ b/src/fastq_file.hpp @@ -7,13 +7,85 @@ #include #include -#include #include #include #include +#ifndef USE_MMAP_FOR_PLAIN_FASTQ +#include +#include +#include +#include // IWYU pragma: keep +#include +#else // NOT USE_MMAP_FOR_PLAIN_FASTQ +#include +#endif struct task_queue; +#ifndef USE_MMAP_FOR_PLAIN_FASTQ + +class fastq_file { + static constexpr auto min_buf_size = 256 * 1024; + std::vector buffer; + std::vector::iterator cursor; + std::vector::iterator last; + std::unique_ptr in; + +public: + fastq_file(const std::string &filename, const std::int64_t buf_size); + operator bool() const { return cursor < last || !std::feof(in.get()); } + + // clang-format off + fastq_file(const fastq_file &) = delete; + auto operator=(const fastq_file &) -> fastq_file & = delete; + auto operator=(fastq_file &&) noexcept -> fastq_file & = delete; + fastq_file(fastq_file &&) noexcept = default; + ~fastq_file() = default; + // clang-format on + + friend auto + reset(fastq_file &reads_file) -> void; + + friend auto + make_tasks(fastq_file &reads_file, + const std::int64_t n_threads, + const std::int32_t file_id, + task_queue &tq, + std::atomic_int32_t &n_tasks) -> void; + +private: + auto + reset() -> void { + buffer.clear(); + buffer.shrink_to_fit(); + cursor = std::begin(buffer); + last = cursor; + } + + auto + get_chunks(const std::int64_t n_chunks, + const std::int32_t file_id, + task_queue &tq, + std::atomic_int32_t &n_tasks) -> void; + + auto + make_tasks(const std::int64_t n_chunks, + const std::int32_t file_id, + task_queue &tq, + std::atomic_int32_t &n_tasks) -> void; + + auto + shift_output_buffer() -> void; + + auto + load_next() -> void; + + [[nodiscard]] auto + skip_header() -> bool; +}; + +#else // USE_MMAP_FOR_PLAIN_FASTQ + struct fastq_file { static constexpr auto min_buf_size = 256 * 1024; std::int64_t target_length{}; @@ -83,6 +155,8 @@ struct fastq_file { std::atomic_int32_t &n_tasks) -> void; }; +#endif // USE_MMAP_FOR_PLAIN_FASTQ + [[nodiscard]] auto estimate_n_reads_fastq(const std::string &filename) -> std::tuple; From a73911db600e23209ae4095c6a78ed5bd3d120e8 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Tue, 27 Oct 2026 16:08:35 -0700 Subject: [PATCH 33/40] src/falco_analyzer.cpp: using a preprocessor guard to distinguish to distinguish between builds using mmapped files and not --- src/falco_analyzer.cpp | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/falco_analyzer.cpp b/src/falco_analyzer.cpp index b13f9b7..bfda9e7 100644 --- a/src/falco_analyzer.cpp +++ b/src/falco_analyzer.cpp @@ -55,10 +55,12 @@ analyze(const std::uint32_t n_threads, for (const auto th_id : std::views::iota(0u, n_threads)) workers.emplace_back([&, n_threads, th_id] { try { +#ifdef USE_MMAP_FOR_PLAIN_FASTQ /// ADS: the 'fq_buffer' is used to make a thread-local copy of mmap /// data so it can be modified without invalidating all /// memories. Declare it here so we have one per thread. std::vector fq_buffer; +#endif // USE_MMAP_FOR_PLAIN_FASTQ auto &res = results[th_id]; for (const auto [file_id, info] : falco::views::enumerate(infos)) res[file_id].init(mode, info, dups_init[file_id]); @@ -67,8 +69,14 @@ analyze(const std::uint32_t n_threads, if (tq.is_finished()) return; auto [file_id, task] = tq.pop_and_release_lock(tq_lock); + // clang-format off if (std::holds_alternative(task)) - process_reads(res[file_id], std::get(task), fq_buffer); + process_reads(res[file_id], std::get(task) +#ifdef USE_MMAP_FOR_PLAIN_FASTQ + , fq_buffer +#endif // USE_MMAP_FOR_PLAIN_FASTQ + ); + // clang-format on else if (std::holds_alternative(task)) process_reads(res[file_id], std::get(task)); else if (std::holds_alternative(task)) From dee904218f39884d86e64d4a852e67247436a767 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Tue, 27 Oct 2026 16:09:02 -0700 Subject: [PATCH 34/40] src/results_collector.hpp: using a preprocessor guard to distinguish to distinguish between builds using mmapped files and not --- src/results_collector.hpp | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/src/results_collector.hpp b/src/results_collector.hpp index 10a56f8..431bee8 100644 --- a/src/results_collector.hpp +++ b/src/results_collector.hpp @@ -196,17 +196,26 @@ accumulate_results(std::vector> &r, } } +// clang-format off inline auto process_reads(results_collector &results, - const fq_task_t &task, - std::vector &fq_buffer) { + const fq_task_t &task +#ifdef USE_MMAP_FOR_PLAIN_FASTQ + , std::vector &fq_buffer +#endif // USE_MMAP_FOR_PLAIN_FASTQ +) { +#ifdef USE_MMAP_FOR_PLAIN_FASTQ const auto task_size = std::size(task); if (std::size(fq_buffer) < task_size) fq_buffer.resize(task_size); std::memcpy(std::data(fq_buffer), task.beg, task_size); const auto beg = std::to_address(std::begin(fq_buffer)); results.process_reads_fq(beg, beg + task_size); +#else + results.process_reads_fq(task.beg, task.end); +#endif // USE_MMAP_FOR_PLAIN_FASTQ } +// clang-format on inline auto process_reads(results_collector &results, const fq_gz_task_t &task) { From 9f013e5a405e0615f0b4da6dd97a2ce5d86499f2 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Tue, 27 Oct 2026 16:09:31 -0700 Subject: [PATCH 35/40] README.md: adding explanation about building falco to use memory mapping for plain FASTQ files --- README.md | 34 ++++++++++++++++++++++++++++------ 1 file changed, 28 insertions(+), 6 deletions(-) diff --git a/README.md b/README.md index 4c789da..8b46b16 100644 --- a/README.md +++ b/README.md @@ -20,7 +20,7 @@ forward. You will be able to find binaries for Linux and macOS with the releases. Example: -``` +```console falco -o output input.fq ``` This generates 3 files in the directory named `output` (creating it if needed): @@ -35,7 +35,7 @@ This generates 3 files in the directory named `output` (creating it if needed): If you use [anaconda](https://anaconda.org) to manage your packages, and the `conda` binary is in your path, you can install the most recent release of `falco` by running -``` +```console conda install -c bioconda falco ``` @@ -56,7 +56,7 @@ conda install -c bioconda falco you need: `.github/workflows/linux-release.yml` and `.github/workflows/macos-release.yml` From the root of the repo if you have all dependencies: -``` +```console cmake -B build -DCMAKE_CXX_COMPILER=g++ -DCMAKE_BUILD_TYPE=Release cmake --build build -j8 # use -j for more cores when building ``` @@ -70,11 +70,11 @@ specified directly. ### Linux detailed instructions I'm explaining this via a clean Ubuntu instance in docker: -``` +```console docker run -it ubuntu:latest bash ``` Inside the docker: -``` +```console export DEBIAN_FRONTEND=noninteractive && apt-get update && apt-get install -y --no-install-recommends \ @@ -102,7 +102,7 @@ you can find them here: `.github/workflows/linux-release.yml` I don't have the same ability to test with clean OS images for macOS (suggestions welcome). The best I can do is use the GitHub macOS runners, which already have some of the dependencies installed. Here is what works: -``` +```console brew install libdeflate htslib samtools && # samtools for testing git clone https://github.com/smithlabcode/falco.git && cd falco && @@ -117,6 +117,28 @@ GitHub's macOS look Note: ISA-L was designed for Intel hardware. It works on Apple silicon, but only gives 1-2% speedup in my tests. +### Potential speedups for plain FASTQ file inputs + +For uncompressed FASTQ files, run times can be improved by building Falco to +memory map the input files. This only helps for *uncompressed FASTQ* files. The +speedup depends on whether your input files remain "hot" during your analysis. +This is determined by your pipeline, hardware, OS settings, etc. One reason a +file might stay "hot" is if it's stil in virtual memory from a previous step in +your pipeline. The speed difference can be almost 2x on some systems, but on +other systems you might see reduced performance. If you use Falco frequently on +plain FASTQ files, this might be worth trying. To enable memory mapped plain +FASTQ files, configure the build with `-DUSE_MMAP_FOR_PLAIN_FASTQ=on`: + +```console +cmake -B build -DCMAKE_CXX_COMPILER=g++ -DUSE_MMAP_FOR_PLAIN_FASTQ=on -DCMAKE_BUILD_TYPE=Release +``` + +If you build Falco from source, there's no harm in testing both ways and seeing +what works best for your pipeline. Just make sure to test using your entire +pipeline and with both orders of testing the two configurations. Otherwise you +risk seeing better performance for whichever is tested second. You can use +`vmtouch` (available through conda) to load or evict files while testing. + ## Changes in Falco v2.0 ### Tiles results From 6c7811659aaaf4a00d7dc2d11d144c7e17a7aebb Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Tue, 27 Oct 2026 16:22:45 -0700 Subject: [PATCH 36/40] CMakeLists.txt: adding configure argument that will turn on USE_MMAP_FOR_PLAIN_FASTQ --- CMakeLists.txt | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 827ae53..b4c0948 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -58,6 +58,14 @@ if(PACKAGE) include(cmake/package.cmake) endif() +option(USE_MMAP_FOR_PLAIN_FASTQ + "Enabling memory mapping for plain FASTQ files" + off +) +if(USE_MMAP_FOR_PLAIN_FASTQ) + add_compile_definitions(USE_MMAP_FOR_PLAIN_FASTQ) +endif() + find_package(LIBDEFLATE REQUIRED) find_package(HTSLIB 1.18 REQUIRED) find_package(ZLIB REQUIRED) From 9e1e8955f3b1250f3db574548684b5df0c0db532 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Tue, 27 Oct 2026 18:02:47 -0700 Subject: [PATCH 37/40] src/fastq_file.hcpp: adding posix_advise sequential for non-macos and swapping filesystem file size for st_size. And fixing headers for build conditions --- src/fastq_file.cpp | 20 ++++++++++++++++---- src/fastq_file.hpp | 2 ++ 2 files changed, 18 insertions(+), 4 deletions(-) diff --git a/src/fastq_file.cpp b/src/fastq_file.cpp index 08544c8..98694a2 100644 --- a/src/fastq_file.cpp +++ b/src/fastq_file.cpp @@ -5,10 +5,15 @@ #include "fqrec.hpp" #include "task_queue.hpp" -#include -#include +#ifndef __APPLE__ +// ADS: fcntl.h needed on macOS for O_RDONLY +#include // IWYU pragma: keep +#endif #include +#ifdef USE_MMAP_FOR_PLAIN_FASTQ +#include #include +#endif // USE_MMAP_FOR_PLAIN_FASTQ #include #include @@ -16,7 +21,6 @@ #include #include #include -#include #include #include #include @@ -88,6 +92,9 @@ fastq_file::fastq_file(const std::string &filename, if (!in) std::system_error(std::make_error_code(std::errc(errno)), "failed to read file"); +#ifndef __APPLE__ + posix_fadvise(fileno(in.get()), 0, 0, POSIX_FADV_SEQUENTIAL); +#endif cursor = std::begin(buffer); last = std::begin(buffer); } @@ -166,11 +173,16 @@ fastq_file::fastq_file(const std::string &filename, const std::int64_t target_length) : target_length{ std::max(target_length, static_cast(min_buf_size))}, - filesize{static_cast(std::filesystem::file_size(filename))}, fd{open(std::data(filename), O_RDONLY, 0)} { if (fd < 0) throw std::system_error(std::make_error_code(std::errc(errno)), "failed to open file: " + filename); +#ifndef __APPLE__ + posix_fadvise(fd, 0, 0, POSIX_FADV_SEQUENTIAL); +#endif + struct stat buf{}; + fstat(fd, &buf); + filesize = buf.st_size; } static inline auto diff --git a/src/fastq_file.hpp b/src/fastq_file.hpp index 9d13d18..5aeb267 100644 --- a/src/fastq_file.hpp +++ b/src/fastq_file.hpp @@ -3,7 +3,9 @@ #ifndef SRC_FASTQ_FILE_HPP_ #define SRC_FASTQ_FILE_HPP_ +#ifdef USE_MMAP_FOR_PLAIN_FASTQ #include +#endif #include #include From 3c11c97bd4dec868d787ecc6dc5d465d3732d841 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Tue, 27 Oct 2026 18:21:33 -0700 Subject: [PATCH 38/40] src/fastq_file.cpp: fixing prepro guard around fcntl.h on apple --- src/fastq_file.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fastq_file.cpp b/src/fastq_file.cpp index 98694a2..f374929 100644 --- a/src/fastq_file.cpp +++ b/src/fastq_file.cpp @@ -5,7 +5,7 @@ #include "fqrec.hpp" #include "task_queue.hpp" -#ifndef __APPLE__ +#ifdef __APPLE__ // ADS: fcntl.h needed on macOS for O_RDONLY #include // IWYU pragma: keep #endif From 6e6e398f177aeda21e61e68403487cdae7613a6f Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Tue, 27 Oct 2026 18:25:52 -0700 Subject: [PATCH 39/40] src/falco.cpp: attempting to make the message about ACGTN more clear --- src/falco.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/falco.cpp b/src/falco.cpp index 198eb4a..a94d478 100644 --- a/src/falco.cpp +++ b/src/falco.cpp @@ -7,8 +7,8 @@ static constexpr auto description = - Prefer bgzip over gzip when compressing to gz format. Get it with these packages: tabix (apt), bgzip (pip and dnf) or htslib (brew and conda). -- Falco assumes nucleotides symbols ACGT and N. If you modify the input data - with other symbols, they are treated as N, but ACGT are case-insensitive. +- Falco assumes nucleotides symbols ACGT and N (upper or lowercase). Any other + symbols are counted as N. EXAMPLES: From 0836bf1ef50559b3fc5a84f2d00c3f469ae74a90 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Tue, 27 Oct 2026 18:37:12 -0700 Subject: [PATCH 40/40] src/fastq_file.cpp: apparently fcntl.h is needed on systems besides mine, so including it generally --- src/fastq_file.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/fastq_file.cpp b/src/fastq_file.cpp index f374929..c634a1e 100644 --- a/src/fastq_file.cpp +++ b/src/fastq_file.cpp @@ -5,10 +5,7 @@ #include "fqrec.hpp" #include "task_queue.hpp" -#ifdef __APPLE__ -// ADS: fcntl.h needed on macOS for O_RDONLY #include // IWYU pragma: keep -#endif #include #ifdef USE_MMAP_FOR_PLAIN_FASTQ #include