From 6b70baf184d40d2db025807f1cb0f68644e7a41a Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Thu, 3 Sep 2026 07:29:32 -0700 Subject: [PATCH 01/25] src/reads_file.hpp: making a var name more sensible --- src/reads_file.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/reads_file.hpp b/src/reads_file.hpp index 641f5bd..8dff7aa 100644 --- a/src/reads_file.hpp +++ b/src/reads_file.hpp @@ -35,10 +35,10 @@ class reads_file_t { } friend auto - make_tasks(auto &reads_file, const std::int64_t n_chunks, + make_tasks(auto &reads_file, const std::int64_t n_threads, const std::int32_t file_id, task_queue &tq, std::atomic_int32_t &n_tasks) -> void { - reads_file.self_->make_tasks_(n_chunks, file_id, tq, n_tasks); + reads_file.self_->make_tasks_(n_threads, file_id, tq, n_tasks); } private: @@ -56,9 +56,9 @@ class reads_file_t { explicit model(T x) : data_(std::move(x)) {} auto - make_tasks_(const std::int64_t n_chunks, const std::int32_t file_id, + make_tasks_(const std::int64_t n_threads, const std::int32_t file_id, task_queue &tq, std::atomic_int32_t &n_tasks) -> void override { - make_tasks(data_, n_chunks, file_id, tq, n_tasks); + make_tasks(data_, n_threads, file_id, tq, n_tasks); } [[nodiscard]] auto From 31547166208ec7b733ed93bb181bf6597dd1e8ee Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Thu, 3 Sep 2026 07:30:20 -0700 Subject: [PATCH 02/25] src/sam_file.hcpp: adding load next and symmetry in code with fastq_file --- src/sam_file.cpp | 33 +++++++++++++++++++-------------- src/sam_file.hpp | 3 +++ 2 files changed, 22 insertions(+), 14 deletions(-) diff --git a/src/sam_file.cpp b/src/sam_file.cpp index 8fb4ced..6e91745 100644 --- a/src/sam_file.cpp +++ b/src/sam_file.cpp @@ -23,7 +23,7 @@ sam_file::sam_file(const std::string &filename, const std::int64_t buf_size) : std::system_error(std::make_error_code(std::errc(errno)), "failed to read file"); if (!skip_header()) - std::runtime_error("failed to validated SAM file header: " + filename); + std::runtime_error("failed to validate SAM header: " + filename); cursor = std::begin(buffer); last = std::begin(buffer); } @@ -77,14 +77,14 @@ sam_file::get_chunks(const std::int64_t n_chunks, // }; // clang-format on const auto n_bytes_available = std::distance(beg_itr, end_itr); - const auto chunk_size = (n_bytes_available + n_chunks - 1) / n_chunks; + const auto [chunk_size, remainder] = std::div(n_bytes_available, n_chunks); assert(n_chunks > 0); - auto start_pos = beg_itr; - auto chunk_end = beg_itr; + 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 = fwd_to_read_start(start_pos); - auto stop_pos = start_pos + chunk_size; - chunk_end = fwd_to_read_start(stop_pos); + const auto chunk_beg = fwd_to_read_start(start_itr); + const auto stop_itr = start_itr + chunk_size + (chunk_idx < remainder); + chunk_end = fwd_to_read_start(stop_itr); if (chunk_idx + 1 == n_chunks) { // final chunk has only full records const auto prev_start = rev_to_read_start(chunk_end); const auto n_trailing = std::count(prev_start, end_itr, '\n'); @@ -93,11 +93,21 @@ sam_file::get_chunks(const std::int64_t n_chunks, // } ++n_tasks; tq.push(file_id, sam_task_t(chunk_beg, chunk_end)); - start_pos = stop_pos; + start_itr = stop_itr; } cursor = chunk_end; } +auto +sam_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 SAM file"); +} + auto sam_file::make_tasks(const std::int64_t n_chunks, // const std::int32_t file_id, // @@ -105,11 +115,6 @@ sam_file::make_tasks(const std::int64_t n_chunks, // std::atomic_int32_t &n_tasks) -> void { n_tasks = 1; // for current task, which makes more tasks shift_output_buffer(); - const auto n_bytes = std::distance(last, std::end(buffer)); - const auto r = 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 SAM file"); - last += static_cast(r); + load_next(); get_chunks(n_chunks, file_id, tq, n_tasks); } diff --git a/src/sam_file.hpp b/src/sam_file.hpp index 03f547a..1cf562f 100644 --- a/src/sam_file.hpp +++ b/src/sam_file.hpp @@ -53,6 +53,9 @@ class sam_file { auto shift_output_buffer() -> void; + auto + load_next() -> void; + [[nodiscard]] auto skip_header() -> bool; }; From 0fc8941b7f834f56dc1dd08a5c83fe30a2114f7a Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Thu, 3 Sep 2026 07:35:26 -0700 Subject: [PATCH 03/25] CMakeLists.txt: adding fastq_stdin object --- CMakeLists.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 552d396..5c52a6a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -83,6 +83,7 @@ add_library(base_groups OBJECT src/base_groups.cpp) add_library(bgzf_block OBJECT src/bgzf_block.cpp) add_library(bgzf_reader OBJECT src/bgzf_reader.cpp) add_library(fastq_file OBJECT src/fastq_file.cpp) +add_library(fastq_stdin OBJECT src/fastq_stdin.cpp) add_library(fastq_gz_file OBJECT src/fastq_gz_file.cpp) add_library(fastq_bgzf_file OBJECT src/fastq_bgzf_file.cpp) add_library(bam_header OBJECT src/bam_header.cpp) @@ -131,6 +132,7 @@ target_link_libraries(falco PRIVATE original_duplicates kmer_counter fastq_file + fastq_stdin fastq_gz_file fastq_bgzf_file bam_file From 91d41494751cbab7adea9b159185a640dd9a8500 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Thu, 3 Sep 2026 07:54:46 -0700 Subject: [PATCH 04/25] src/fastq_stdin.hcpp: adding sources for reading fastq from stdin --- src/fastq_stdin.cpp | 111 ++++++++++++++++++++++++++++++++++++++++++++ src/fastq_stdin.hpp | 81 ++++++++++++++++++++++++++++++++ 2 files changed, 192 insertions(+) create mode 100644 src/fastq_stdin.cpp create mode 100644 src/fastq_stdin.hpp diff --git a/src/fastq_stdin.cpp b/src/fastq_stdin.cpp new file mode 100644 index 0000000..4056915 --- /dev/null +++ b/src/fastq_stdin.cpp @@ -0,0 +1,111 @@ +// SPDX-License-Identifier: MIT; Copyright 2026 Andrew D Smith + +#include "fastq_stdin.hpp" +#include "falco_utils.hpp" +#include "fqrec.hpp" +#include "task_queue.hpp" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +fastq_stdin::fastq_stdin(const std::int64_t buf_size) : + buffer(buf_size + min_buf_size), cursor{std::begin(buffer)}, + last{std::begin(buffer)} {} + +[[nodiscard]] auto +estimate_n_reads_fastq_stdin(const std::string &) + -> std::tuple { + static constexpr auto assumed_n_reads = 10'000'000; + static constexpr auto assumed_read_len = 150; + static constexpr auto assumed_filesize = 0; + return {assumed_n_reads, assumed_read_len, assumed_filesize}; +} + +auto +fastq_stdin::get_chunks(const std::int64_t n_chunks, const std::int32_t file_id, + task_queue &tq, std::atomic_int32_t &n_tasks) -> void { + static constexpr auto rec_lines = 4; // FASTQ + const auto beg_itr = std::begin(buffer); + const auto end_itr = last; + assert(n_chunks > 0); + // 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 fwd_to_read_start = [&](auto pos) { + if (pos == beg_itr) return pos; + while (pos < end_itr && not_read_start(pos)) ++pos; + return pos; + }; + const auto rev_to_read_start = [&](auto pos) { + while (pos > beg_itr && (pos == end_itr || not_read_start(pos))) --pos; + return pos; + }; + // clang-format on + const auto n_bytes_available = std::distance(beg_itr, last); + const auto [chunk_size, remainder] = std::div(n_bytes_available, n_chunks); + std::vector> chunks(n_chunks); + 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 = fwd_to_read_start(start_itr); + const auto stop_itr = start_itr + chunk_size + (chunk_idx < remainder); + chunk_end = fwd_to_read_start(stop_itr); + if (chunk_idx + 1 == n_chunks) + if (const auto prev = rev_to_read_start(chunk_end); + std::count(prev, std::end(buffer), '\n') < rec_lines) + chunk_end = prev; + ++n_tasks; + tq.push(file_id, + fq_task_t(std::to_address(chunk_beg), std::to_address(chunk_end))); + start_itr = stop_itr; + } + cursor = chunk_end; +} + +auto +fastq_stdin::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_stdin::load_next() -> void { + const auto n_bytes = std::distance(last, std::end(buffer)); + const auto r = read(0, std::to_address(last), n_bytes); + if (r == -1) + std::system_error(std::make_error_code(std::errc(errno)), + "error reading fastq from stdin"); + if (r == 0) + hit_eof = true; + last += static_cast(r); +} + +auto +fastq_stdin::make_tasks(const std::int64_t n_chunks, // + const std::int32_t file_id, // + task_queue &tq, // + std::atomic_int32_t &n_tasks) -> void { + n_tasks = 1; // for current task, which makes more tasks + shift_output_buffer(); + load_next(); + get_chunks(n_chunks, file_id, tq, n_tasks); +} diff --git a/src/fastq_stdin.hpp b/src/fastq_stdin.hpp new file mode 100644 index 0000000..2eac214 --- /dev/null +++ b/src/fastq_stdin.hpp @@ -0,0 +1,81 @@ +// SPDX-License-Identifier: MIT; Copyright 2026 Andrew D Smith + +#ifndef SRC_FASTQ_STDIN_HPP_ +#define SRC_FASTQ_STDIN_HPP_ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +struct task_queue; + +struct fastq_stdin { + static constexpr auto min_buf_size = 64 * 1024; + std::vector buffer; + std::vector::iterator cursor; + std::vector::iterator last; + bool hit_eof{}; + + fastq_stdin(const std::int64_t buf_size); + operator bool() const { return cursor < last || !hit_eof; } + + // clang-format off + fastq_stdin(const fastq_stdin &) = delete; + auto operator=(const fastq_stdin &) -> fastq_stdin & = delete; + auto operator=(fastq_stdin &&) noexcept -> fastq_stdin & = delete; + fastq_stdin(fastq_stdin &&) noexcept = default; + ~fastq_stdin() = default; + // clang-format on + + auto + reset() -> void { + buffer.clear(); + buffer.shrink_to_fit(); + cursor = std::begin(buffer); + last = std::begin(buffer); + } + + 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; + +private: + 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 + shift_output_buffer() -> void; + + auto + load_next() -> void; +}; + +[[nodiscard]] auto +estimate_n_reads_fastq(const std::string &filename) + -> std::tuple; + +inline auto +make_tasks(fastq_stdin &reads_file, // + const std::int64_t n_threads, // + const std::int32_t file_id, // + task_queue &tq, // + 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; + reads_file.make_tasks(n_chunks, file_id, tq, n_tasks); +} + +inline auto +reset(fastq_stdin &reads_file) -> void { + reads_file.reset(); +} + +#endif // SRC_FASTQ_STDIN_HPP_ From e6ae0ac247fdd3204834f2f1ff0a1113f3c84c92 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Thu, 3 Sep 2026 08:01:49 -0700 Subject: [PATCH 05/25] src/sam_file.cpp and src/fastq_stdin.cpp: preventing compiler from complaints about mismatching types in std::div --- src/fastq_stdin.cpp | 2 +- src/sam_file.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/fastq_stdin.cpp b/src/fastq_stdin.cpp index 4056915..2595e77 100644 --- a/src/fastq_stdin.cpp +++ b/src/fastq_stdin.cpp @@ -58,7 +58,7 @@ fastq_stdin::get_chunks(const std::int64_t n_chunks, const std::int32_t file_id, return pos; }; // clang-format on - const auto n_bytes_available = std::distance(beg_itr, last); + const std::int64_t n_bytes_available = std::distance(beg_itr, last); const auto [chunk_size, remainder] = std::div(n_bytes_available, n_chunks); std::vector> chunks(n_chunks); auto start_itr = beg_itr; diff --git a/src/sam_file.cpp b/src/sam_file.cpp index 6e91745..276d1ab 100644 --- a/src/sam_file.cpp +++ b/src/sam_file.cpp @@ -76,7 +76,7 @@ sam_file::get_chunks(const std::int64_t n_chunks, // return p; }; // clang-format on - const auto n_bytes_available = std::distance(beg_itr, end_itr); + const std::int64_t n_bytes_available = std::distance(beg_itr, end_itr); const auto [chunk_size, remainder] = std::div(n_bytes_available, n_chunks); assert(n_chunks > 0); auto start_itr = beg_itr; From 03661a20f5145a5c80881b7b087b87f31195e77a Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Thu, 3 Sep 2026 14:24:57 -0700 Subject: [PATCH 06/25] src/falco_file_format.hpp: format now uses names of file_format values --- src/falco_file_format.hpp | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/src/falco_file_format.hpp b/src/falco_file_format.hpp index 9762a42..8ff6118 100644 --- a/src/falco_file_format.hpp +++ b/src/falco_file_format.hpp @@ -23,6 +23,19 @@ enum class file_format : std::uint8_t { bam, }; +// clang-format off +static constexpr auto file_format_names_impl = std::array{ + "unknown", // + "fastq", // + "fastq_gz", // + "fastq_bgzf", // + "sam", // + "bam", // +}; +// clang-format on + +static constexpr std::span file_format_names = file_format_names_impl; + // NOLINTNEXTLINE NLOHMANN_JSON_SERIALIZE_ENUM( // file_format, // @@ -62,7 +75,7 @@ struct std::formatter : std::formatter { auto format(const falco::file_format &f, auto &ctx) const { return std::formatter::format( - std::to_string(std::to_underlying(f)), ctx); + falco::file_format_names[std::to_underlying(f)], ctx); } }; From 61c9e514572337df9a70620e20ffa02584ba3a9d Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Fri, 4 Sep 2026 14:45:53 -0700 Subject: [PATCH 07/25] src/fastq_stdin.cpp: using multiple reads when necessary to fill the buffer from stdin --- src/fastq_stdin.cpp | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/src/fastq_stdin.cpp b/src/fastq_stdin.cpp index 2595e77..1da7bfb 100644 --- a/src/fastq_stdin.cpp +++ b/src/fastq_stdin.cpp @@ -60,7 +60,6 @@ fastq_stdin::get_chunks(const std::int64_t n_chunks, const std::int32_t file_id, // clang-format on const std::int64_t n_bytes_available = std::distance(beg_itr, last); const auto [chunk_size, remainder] = std::div(n_bytes_available, n_chunks); - std::vector> chunks(n_chunks); auto start_itr = beg_itr; auto chunk_end = start_itr; for (const auto chunk_idx : std::views::iota(0, n_chunks)) { @@ -89,14 +88,16 @@ fastq_stdin::shift_output_buffer() -> void { auto fastq_stdin::load_next() -> void { - const auto n_bytes = std::distance(last, std::end(buffer)); - const auto r = read(0, std::to_address(last), n_bytes); - if (r == -1) - std::system_error(std::make_error_code(std::errc(errno)), - "error reading fastq from stdin"); - if (r == 0) - hit_eof = true; - last += static_cast(r); + auto space = std::distance(last, std::end(buffer)); + std::int64_t n{1}; + while (space > 0 && (n = read(0, std::to_address(last), space)) != 0) { + if (n == -1) + std::system_error(std::make_error_code(std::errc(errno)), + "error reading fastq from stdin"); + space -= n; + last += n; + } + hit_eof = (n == 0); } auto From c1fc0aa7e1b7f245b1359a4576b7c20e0483ba01 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Fri, 4 Sep 2026 14:46:32 -0700 Subject: [PATCH 08/25] src/tile_processor.cpp: when identifying tile position within read names, checking multiple reads, and ensuring they give integer values for tile ids --- src/tile_processor.cpp | 88 +++++++++++++++++++++++++++++++----------- 1 file changed, 65 insertions(+), 23 deletions(-) diff --git a/src/tile_processor.cpp b/src/tile_processor.cpp index fb3e9ed..a2cf8f8 100644 --- a/src/tile_processor.cpp +++ b/src/tile_processor.cpp @@ -35,33 +35,52 @@ tile_processor::init(const file_info &info) -> void { } [[nodiscard]] auto -get_name_fastq(const std::string &filename) -> std::string { +get_names_fastq(const std::string &filename) -> std::vector { + static constexpr auto n_records = 100; + static constexpr auto n_lines_per_record = 4; + static constexpr auto n_total_lines = n_records * n_lines_per_record; + static constexpr auto name_line = 0; std::unique_ptr in(bgzf_open(std::data(filename), "r"), &bgzf_close); if (!in) throw std::runtime_error("failed to open gz file: " + filename); - kstring_t str = KS_INITIALIZE; - const auto r = bgzf_getline(in.get(), '\n', &str); - if (r < 0) - throw std::runtime_error("failed to read line from: " + filename); - std::string line(str.s, str.l); - ks_free(&str); - return line; + std::vector names; + std::uint32_t line_count{}; + for (auto i = 0; i < n_total_lines; ++i) { + kstring_t str = KS_INITIALIZE; + const auto r = bgzf_getline(in.get(), '\n', &str); + if (r == -1) // EOF (htslib/bgzf.h) + break; + if (r < -1) // ERROR (htslib/bgzf.h) + throw std::runtime_error("failed to read line from: " + filename); + if (line_count % n_lines_per_record == name_line) + names.emplace_back(str.s, str.l); + ks_free(&str); + ++line_count; + } + return names; } [[nodiscard]] auto -get_name_bam(const std::string &filename) -> std::string { +get_names_bam(const std::string &filename) -> std::vector { + static constexpr auto n_records = 100; std::unique_ptr in( hts_open(std::data(filename), "r"), &hts_close); if (!in) throw std::runtime_error("failed to open BAM/SAM file: " + filename); std::unique_ptr h(sam_hdr_read(in.get()), &sam_hdr_destroy); - std::unique_ptr b(bam_init1(), &bam_destroy1); - const auto r = sam_read1(in.get(), h.get(), b.get()); // -1 on EOF - if (r < -1) - throw std::runtime_error("failed reading bam record"); - return bam_get_qname(b); + std::vector names; + for (auto i = 0; i < n_records; ++i) { + std::unique_ptr b(bam_init1(), &bam_destroy1); + const auto r = sam_read1(in.get(), h.get(), b.get()); // -1 on EOF + if (r == -1) // EOF (htslib/bgzf.h) + break; + if (r < -1) // ERROR (htslib/bgzf.h) + throw std::runtime_error("failed reading BAM/SAM record: " + filename); + names.emplace_back(bam_get_qname(b)); + } + return names; } [[nodiscard]] auto @@ -168,7 +187,7 @@ tile_processor::add_and_consume( const auto pair_plus = [](const auto &a, const auto &b) { return std::pair{a.first + b.first, a.second + b.second}; }; - for (auto &[rhs_tile_id, rhs_qual] : rhs.quals) { + for (auto &&[rhs_tile_id, rhs_qual] : rhs.quals) { const auto quals_itr = quals.find(rhs_tile_id); if (quals_itr != std::end(quals)) { auto &curr_qual = quals_itr->second; @@ -199,18 +218,41 @@ get_tile_info(const std::string &filename) -> std::uint32_t { throw std::runtime_error("failed to open file: " + filename); const auto hts_fmt = hts_get_format(fp.get()); - if (hts_fmt->format != fastq_format && hts_fmt->format != bam && + if (hts_fmt->format != fastq_format && // + hts_fmt->format != bam && // hts_fmt->format != sam) return 0; - const auto line = (hts_fmt->format == bam || hts_fmt->format == sam) - ? get_name_bam(filename) - : get_name_fastq(filename); + const auto names = (hts_fmt->format == bam || hts_fmt->format == sam) + ? get_names_bam(filename) + : get_names_fastq(filename); - const auto colons_found = std::ranges::count(line, ':'); - return colons_found >= colon_cutoff_1 - ? colon_cutoff_1_val - : (colons_found >= colon_cutoff_2 ? colon_cutoff_2_val : 0); + const auto n_colons = + std::ranges::min(names | std::views::transform([](const auto &x) { + return std::ranges::count(x, ':'); + })); + + const auto colon_cutoff_val = + (n_colons >= colon_cutoff_1) + ? colon_cutoff_1_val + : (n_colons >= colon_cutoff_2 ? colon_cutoff_2_val : 0); + // now verify that they are all valid + if (colon_cutoff_val > 0) + for (const auto &name : names) { + auto tile_itr = std::cbegin(name); + auto colon_count = 0; + while (colon_count < colon_cutoff_val && tile_itr != std::cend(name)) + colon_count += (*tile_itr++ == ':'); + std::uint32_t curr_tile_id{}; + const auto [_, ec] = + std::from_chars(std::to_address(tile_itr), + std::to_address(std::cend(name)), curr_tile_id); + if (ec != std::errc{}) + throw std::system_error( + std::make_error_code(ec), + "error identifying numerical tile id; rerun with --no-tiles"); + } + return colon_cutoff_val; } [[nodiscard]] auto From 7ece5ed9a8aef560b87ed301a9a3495501bffb46 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Fri, 4 Sep 2026 16:48:00 -0700 Subject: [PATCH 09/25] src/sam_stdin.hcpp: adding sources to read SAM from stdin --- src/sam_stdin.cpp | 130 ++++++++++++++++++++++++++++++++++++++++++++++ src/sam_stdin.hpp | 84 ++++++++++++++++++++++++++++++ 2 files changed, 214 insertions(+) create mode 100644 src/sam_stdin.cpp create mode 100644 src/sam_stdin.hpp diff --git a/src/sam_stdin.cpp b/src/sam_stdin.cpp new file mode 100644 index 0000000..56424cd --- /dev/null +++ b/src/sam_stdin.cpp @@ -0,0 +1,130 @@ +// SPDX-License-Identifier: MIT; Copyright 2026 Andrew D Smith + +#include "sam_stdin.hpp" +#include "samrec.hpp" +#include "task_queue.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +sam_stdin::sam_stdin(const std::int64_t buf_size) : + buffer(buf_size + min_buf_size), cursor{std::begin(buffer)}, + last{std::begin(buffer)} { + if (!skip_header()) + std::runtime_error("failed to validate SAM header: stdin"); +} + +[[nodiscard]] auto +estimate_n_reads_sam_stdin(const std::string &) + -> std::tuple { + static constexpr auto assumed_n_reads = 10'000'000; + static constexpr auto assumed_read_len = 150; + static constexpr auto assumed_filesize = 0; + return {assumed_n_reads, assumed_read_len, assumed_filesize}; +} + +[[nodiscard]] auto +sam_stdin::skip_header() -> bool { + bool pre_header = true; + while (cursor == std::cbegin(buffer)) { + load_next(); + if (pre_header && buffer[0] != '@') + break; + pre_header = false; + const auto lim = std::distance(std::begin(buffer), last); + for (auto i = 1; i < lim && cursor == std::begin(buffer); ++i) + if (buffer[i - 1] == '\n' && buffer[i] != '@') + cursor = std::begin(buffer) + i; + if (cursor == std::cbegin(buffer)) { + // corner case of '\n' and non-@ in different buffers + last = std::end(buffer); + cursor = std::prev(last); + shift_output_buffer(); // in prep for subsequent load-next + } + } + return true; +} + +auto +sam_stdin::get_chunks(const std::int64_t n_chunks, // + const std::int32_t file_id, // + task_queue &tq, // + std::atomic_int32_t &n_tasks) -> void { + assert(n_chunks > 0); + const auto beg_itr = std::begin(buffer); + const auto end_itr = last; + // clang-format off + const auto fwd_to_read_start = [&](auto p) { + if (p == beg_itr) return p; + while (p != end_itr && *p != '\n') ++p; + if (p != end_itr) ++p; + return p; + }; + const auto rev_to_read_start = [&](auto p) { + while (p != beg_itr && *(p - 1) != '\n') --p; + return p; + }; + // clang-format on + const std::int64_t n_bytes_available = std::distance(beg_itr, end_itr); + const 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 = fwd_to_read_start(start_itr); + const auto stop_itr = start_itr + chunk_size + (chunk_idx < remainder); + chunk_end = fwd_to_read_start(stop_itr); + if (chunk_idx + 1 == n_chunks) { // final chunk has only full records + const auto prev_start = rev_to_read_start(chunk_end); + const auto n_trailing = std::count(prev_start, end_itr, '\n'); + if (n_trailing == 0) + chunk_end = prev_start; + } + ++n_tasks; + tq.push(file_id, sam_task_t(chunk_beg, chunk_end)); + start_itr = stop_itr; + } + cursor = chunk_end; +} + +auto +sam_stdin::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 +sam_stdin::load_next() -> void { + auto space = std::distance(last, std::end(buffer)); + std::int64_t n{1}; + while (space > 0 && (n = read(0, std::to_address(last), space)) != 0) { + if (n == -1) + std::system_error(std::make_error_code(std::errc(errno)), + "error reading fastq from stdin"); + space -= n; + last += n; + } + hit_eof = (n == 0); +} + +auto +sam_stdin::make_tasks(const std::int64_t n_chunks, // + const std::int32_t file_id, // + task_queue &tq, // + std::atomic_int32_t &n_tasks) -> void { + n_tasks = 1; // for current task, which makes more tasks + shift_output_buffer(); + load_next(); + get_chunks(n_chunks, file_id, tq, n_tasks); +} diff --git a/src/sam_stdin.hpp b/src/sam_stdin.hpp new file mode 100644 index 0000000..9566f15 --- /dev/null +++ b/src/sam_stdin.hpp @@ -0,0 +1,84 @@ +// SPDX-License-Identifier: MIT; Copyright 2026 Andrew D Smith + +#ifndef SRC_SAM_STDIN_HPP_ +#define SRC_SAM_STDIN_HPP_ + +#include +#include +#include +#include +#include +#include +#include +#include + +struct task_queue; + +class sam_stdin { + static constexpr auto min_buf_size = 64 * 1024; + std::vector buffer; + std::vector::iterator cursor; + std::vector::iterator last; + bool hit_eof{}; + +public: + sam_stdin(const std::int64_t buf_size); + operator bool() const { return cursor < last || !hit_eof; } + + // clang-format off + sam_stdin(const sam_stdin &) = delete; + auto operator=(const sam_stdin &) -> sam_stdin & = delete; + auto operator=(sam_stdin &&) noexcept -> sam_stdin & = delete; + sam_stdin(sam_stdin &&) noexcept = default; + ~sam_stdin() = default; + // clang-format on + + 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 + reset() -> void { + buffer.clear(); + buffer.shrink_to_fit(); + cursor = std::begin(buffer); + last = cursor; + } + +private: + 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 + shift_output_buffer() -> void; + + auto + load_next() -> void; + + [[nodiscard]] auto + skip_header() -> bool; +}; + +[[nodiscard]] auto +estimate_n_reads_sam_stdin(const std::string &) + -> std::tuple; + +inline auto +make_tasks(sam_stdin &reads_file, // + const std::int64_t n_threads, // + const std::int32_t file_id, // + task_queue &tq, // + 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; + reads_file.make_tasks(n_chunks, file_id, tq, n_tasks); +} + +inline auto +reset(sam_stdin &reads_file) -> void { + reads_file.reset(); +} + +#endif // SRC_SAM_STDIN_HPP_ From d04ec7f9b20b72be922cd7f8926b9050249aa7f9 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Fri, 4 Sep 2026 16:48:17 -0700 Subject: [PATCH 10/25] CMakeLists.txt: adding sam_stdin target --- CMakeLists.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 5c52a6a..67251f3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -91,6 +91,7 @@ add_library(bamrec OBJECT src/bamrec.cpp) add_library(bam_file OBJECT src/bam_file.cpp) add_library(samrec OBJECT src/samrec.cpp) add_library(sam_file OBJECT src/sam_file.cpp) +add_library(sam_stdin OBJECT src/sam_stdin.cpp) add_library(falco_file_format OBJECT src/falco_file_format.cpp) add_library(quality_score OBJECT src/quality_score.cpp) add_library(tile_processor OBJECT src/tile_processor.cpp) @@ -138,6 +139,7 @@ target_link_libraries(falco PRIVATE bam_file bamrec sam_file + sam_stdin samrec bam_header falco_file_format From f21ae23f80a0f885b07e3677e72042b1408ce48c Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Fri, 4 Sep 2026 19:35:11 -0700 Subject: [PATCH 11/25] src/fastq_stdin.cpp: added validation for the file format that looks at the first 16k bytes of each buffer load; this probably should be done differently --- src/fastq_stdin.cpp | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/src/fastq_stdin.cpp b/src/fastq_stdin.cpp index 1da7bfb..7e92227 100644 --- a/src/fastq_stdin.cpp +++ b/src/fastq_stdin.cpp @@ -86,6 +86,28 @@ fastq_stdin::shift_output_buffer() -> void { cursor = std::begin(buffer); } +[[nodiscard]] static auto +validate_fastq(const auto &buffer) { + static constexpr auto n_bytes_to_validate = 16 * 1024; + // verify that no two consecutive newlines are followed by a '@' + bool prev_was_ampersand{false}; + auto n_bytes = 0; + for (auto itr = std::cbegin(buffer); itr + 1 != std::cend(buffer); ++itr) { + if (*itr == '\n') { + if (*(itr + 1) == '@') { + if (prev_was_ampersand) + return false; + prev_was_ampersand = true; + } + else + prev_was_ampersand = false; + } + if (++n_bytes == n_bytes_to_validate) + return true; + } + return true; +} + auto fastq_stdin::load_next() -> void { auto space = std::distance(last, std::end(buffer)); @@ -98,6 +120,8 @@ fastq_stdin::load_next() -> void { last += n; } hit_eof = (n == 0); + if (!validate_fastq(buffer)) + throw std::runtime_error("input appears not to be FASTQ"); } auto From 7c670f38a4c4cc298dc6e4c14d60ecf43899b4de Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Fri, 4 Sep 2026 19:39:28 -0700 Subject: [PATCH 12/25] src/falco.cpp: added arguments and logic to enable reading data from stdin, with user-specified format as either fastq or sam, and with user-specified colon position for tile id within read names --- src/falco.cpp | 120 +++++++++++++++++++++++++++++++++++++++----------- 1 file changed, 94 insertions(+), 26 deletions(-) diff --git a/src/falco.cpp b/src/falco.cpp index 6a26b3f..6bbf75d 100644 --- a/src/falco.cpp +++ b/src/falco.cpp @@ -3,7 +3,7 @@ // clang-format off static constexpr auto about = R"(Falco v{})"; static constexpr auto description = - R"(Note: Always use bgzip when compressing files to gz format. It comes with samtools. + R"(Note: Always use bgzip when compressing to gz format. It comes with samtools. EXAMPLES: @@ -32,6 +32,12 @@ Generate a file with duplication info for preseq analysis: Output files 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 +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" + Default configuration files can be found here: {} Use these as templates. Copy and modify them to customize your analysis. @@ -50,6 +56,7 @@ Use these as templates. Copy and modify them to customize your analysis. #include "fastq_bgzf_file.hpp" #include "fastq_file.hpp" #include "fastq_gz_file.hpp" +#include "fastq_stdin.hpp" #include "file_info.hpp" #include "get_binary_dir.hpp" #include "original_duplicates.hpp" @@ -57,6 +64,7 @@ Use these as templates. Copy and modify them to customize your analysis. #include "results_summary.hpp" #include "run_mode.hpp" #include "sam_file.hpp" +#include "sam_stdin.hpp" #include "tile_processor.hpp" #include "CLI11/CLI11.hpp" @@ -142,6 +150,38 @@ make_reads_files(const std::vector &infos, return reads_files; } +[[nodiscard]] static auto +make_reads_file_stdin(const std::vector &infos, + const std::int64_t buf_size) + -> std::vector { + std::vector reads_files; + switch (infos.front().format) { + case falco::file_format::fastq: + reads_files.emplace_back(fastq_stdin(buf_size)); + break; + case falco::file_format::sam: + reads_files.emplace_back(sam_stdin(buf_size)); + break; + default: + throw std::runtime_error("unsupported stdin file format"); + } + return reads_files; +} + +[[nodiscard]] static auto +get_file_info_stdin(const std::vector &names, + const std::pair &ft_tile) + -> std::vector { + file_info info; + info.name = names.front(); + info.format = ft_tile.first; + info.description = std::format("{} from standard input", info.format); + info.size = 0; + info.has_tiles = (ft_tile.second != 0); + info.tile_id_position = ft_tile.second; + return std::vector(1, info); +} + [[nodiscard]] static auto get_file_info(const auto &infiles) { std::vector infos; @@ -150,17 +190,17 @@ get_file_info(const auto &infiles) { const auto tile_id_position = get_tile_info(infile); const bool has_tiles = (tile_id_position != 0); const auto [n_reads_est, read_len_est, filesize] = [&] { - if (input_format == falco::file_format::bam) - return estimate_n_reads_bam(infile); - if (input_format == falco::file_format::sam) - return estimate_n_reads_bam(infile); - if (input_format == falco::file_format::fastq_bgzf) - return estimate_n_reads_fastq_bgzf(infile); - if (input_format == falco::file_format::fastq_gz) - return estimate_n_reads_fastq_gz(infile); - if (input_format == falco::file_format::fastq) - return estimate_n_reads_fastq(infile); - throw std::runtime_error("invalid reads file format"); + // clang-format off + using falco::file_format; + switch (input_format) { + case file_format::bam: return estimate_n_reads_bam(infile); + case file_format::sam: return estimate_n_reads_bam(infile); + case file_format::fastq_bgzf: return estimate_n_reads_fastq_bgzf(infile); + case file_format::fastq_gz: return estimate_n_reads_fastq_gz(infile); + case file_format::fastq: return estimate_n_reads_fastq(infile); + default: throw std::runtime_error("invalid reads file format"); + } + // clang-format on }(); infos.push_back({ .name = std::filesystem::path{infile}.filename().string(), @@ -256,17 +296,22 @@ main(int argc, char *argv[]) { {"k"s, kilobytes}, }); - const auto license_callback = [&](auto) { - std::print("{}", license_text); - throw CLI::Success(); - }; - struct FormatWithoutFlagDefaults : public CLI::Formatter { FormatWithoutFlagDefaults() : Formatter() { CLI::FormatterBase::enable_default_flag_values_ = false; } }; + // Related to reading data from stdin + const auto format_name_map = std::map{ + std::pair{"fq"s, falco::file_format::fastq}, + {"sam"s, falco::file_format::sam}, + }; + auto stdin_info = std::pair{ + falco::file_format::unknown, + 0U, + }; + CLI::App app{std::format(about, VERSION)}; const auto fmt = std::make_shared(); app.formatter(fmt); @@ -277,18 +322,20 @@ main(int argc, char *argv[]) { if (argc >= 2) app.footer(std::format(description, falco::get_share_dir())); + // clang-format off // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers) app.get_formatter()->long_option_alignment_ratio(0.2); app.set_help_flag("-h,--help", "Print more detailed help"); app.set_version_flag("--version", VERSION, "Print program version"); - // clang-format off - app.add_flag("--license", license_callback, "Print full license") + app.add_flag("--license", [&](auto) { std::print("{}", license_text); throw CLI::Success(); }, + "Print full license") ->callback_priority(CLI::CallbackPriority::PreRequirementsCheck); - app.add_option("INFILES", infiles, - "FASTQ (plain, GZIP or BGZF) or BAM/SAM") - ->required() + auto infiles_opt = + app.add_option("INFILES", infiles, + "FASTQ (plain, GZIP or BGZF) or BAM/SAM") ->option_text(" ") - ->check(CLI::ExistingFile); + ->required() + ->check(CLI::ExistingFile, "file_check"); app.add_option("-o,--output", outdir, "Output directory (required)") ->required() ->option_text("DIR"); @@ -321,6 +368,22 @@ main(int argc, char *argv[]) { ->option_text(" ") ->capture_default_str() ->transform(size_from_units); + app.add_option_function>( + "--stdin", + [&](const auto &arg) { // callback is to allow trailing arg to be name + stdin_info = arg; + infiles_opt->get_validator("file_check")->active(false); + infiles_opt->expected(1); + }, + "Read from stdin assuming given format (see help)") + ->option_text("fq|sam[:{4,6}]") + ->delimiter(':') + ->allow_extra_args(false) + ->type_size(1, 2) + ->transform(CLI::CheckedTransformer(format_name_map, CLI::ignore_case) + .application_index(0)) + ->check(CLI::IsMember({4, 6}).application_index(1)) + ->callback_priority(CLI::CallbackPriority::PreRequirementsCheck); app.add_flag("--bisulfite", do_bisulfite, "Assume bisulfite when grading sequence content") ->option_text(" "); @@ -358,6 +421,8 @@ main(int argc, char *argv[]) { } CLI11_PARSE(app, argc, argv); + const bool do_stdin = stdin_info.first != falco::file_format::unknown; + run_mode mode; // declare mode here so we can assign from config file if (!config_file.empty()) load_config_and_set_graders(config_file, mode); @@ -401,7 +466,8 @@ main(int argc, char *argv[]) { adapters_file, adapter_set::n_adapters()); // not const because infos will change later when we can deduce the encoding - auto infos = get_file_info(infiles); + auto infos = do_stdin ? get_file_info_stdin(infiles, stdin_info) + : get_file_info(infiles); // restrict buffer size to avoid using a possibly harmful amount of memory const auto get_sz = [](const auto &i) { return i.size; }; @@ -440,9 +506,11 @@ main(int argc, char *argv[]) { auto dups = do_original_dups ? initialize_original_duplicates(infiles, infos, n_threads) : std::vector{}; + auto reads_files = do_stdin ? make_reads_file_stdin(infos, buffer_size) + : make_reads_files(infos, infiles, buffer_size); auto results = - analyze(n_threads, mode, infos, - make_reads_files(infos, infiles, buffer_size), std::move(dups)); + analyze(n_threads, mode, infos, std::move(reads_files), std::move(dups)); + write_output(mode, infos, outdirs, std::move(results)); if (verbose) From c0d19a61bec0b2c9b3c01b38a73912104406666a Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Fri, 4 Sep 2026 19:41:34 -0700 Subject: [PATCH 13/25] src/fastq_file.hcpp: wrapping the memory mapped buffer in a span so the code to work with the fastq records can be more similar to other file formats. Tests show memory mapping does help with speed, but with enough files analyzed at the same time, memory mapping is slower. These changes will make it easier to remove the memory mapping later. --- src/fastq_file.cpp | 124 +++++++++++++++++++++++++-------------------- src/fastq_file.hpp | 38 ++++++-------- 2 files changed, 84 insertions(+), 78 deletions(-) diff --git a/src/fastq_file.cpp b/src/fastq_file.cpp index cd99daf..be447ea 100644 --- a/src/fastq_file.cpp +++ b/src/fastq_file.cpp @@ -27,47 +27,59 @@ estimate_n_reads_fastq(const std::string &filename) -> std::tuple { static constexpr auto fastq_lines_per_read = 4; - static constexpr auto n_parts = 10; + static constexpr auto n_parts = 10L; static constexpr auto max_part_size = 1024 * 1024; - static const auto page_mask = ~(sysconf(_SC_PAGESIZE) - 1); - const int fd = open(std::data(filename), O_RDONLY, 0); - if (fd < 0) + std::vector buffer(max_part_size); + + std::unique_ptr in( + std::fopen(std::data(filename), "r"), &std::fclose); + if (!in) throw std::system_error(std::make_error_code(std::errc(errno)), "failed to open file: " + filename); struct stat buf{}; - fstat(fd, &buf); + fstat(fileno(in.get()), &buf); const auto filesize = buf.st_size; if (filesize < n_parts) return {{}, {}, filesize}; - const auto part_size = - filesize < n_parts * max_part_size ? filesize / n_parts : max_part_size; - - auto total_newlines = 0ul; - auto readlen_est = 0ul; + const auto [part_size, remainder] = (filesize < n_parts * max_part_size) + ? std::div(filesize, n_parts) + : std::ldiv_t{max_part_size, 0}; + auto n_lines = 0LU; + auto readlen_est = 0LU; + auto offset = 0L; for (auto i = 0; i < n_parts; ++i) { - const auto offset = (i * part_size) & page_mask; - auto raw = mmap(nullptr, part_size, PROT_READ, MAP_PRIVATE, fd, offset); - if (raw == MAP_FAILED) - throw std::system_error(std::make_error_code(std::errc(errno)), - "failed to mmap file"); - const auto data = std::span(static_cast(raw), part_size); - total_newlines += std::ranges::count(data, '\n'); - readlen_est += estimate_read_length_fastq_chunk(std::data(data), part_size); - if (munmap(raw, part_size)) - throw std::system_error(std::make_error_code(std::errc(errno)), - "failed to unmap memory"); + if (std::fseek(in.get(), offset, SEEK_SET)) + std::system_error(std::make_error_code(std::errc(errno)), + "error reading fastq file: " + filename); + const auto r = std::fread(std::data(buffer), 1, part_size, in.get()); + if (std::ferror(in.get())) + std::system_error(std::make_error_code(std::errc(errno)), + "error reading fastq file: " + filename); + n_lines += std::ranges::count(std::span{std::cbegin(buffer), r}, '\n'); + readlen_est += estimate_read_length_fastq_chunk(buffer, r); + offset += part_size + (i < remainder); } - close(fd); readlen_est /= n_parts; const auto n_reads_est = - as_frac(total_newlines, fastq_lines_per_read) * + as_frac(n_lines, fastq_lines_per_read) * as_frac(static_cast(filesize), (part_size * n_parts)); return {static_cast(n_reads_est), readlen_est, filesize}; } +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); +} + static inline auto mmap_fastq(const int fd, const std::int64_t offset, const std::int64_t length, auto &data) { @@ -80,33 +92,34 @@ mmap_fastq(const int fd, const std::int64_t offset, const std::int64_t length, } static inline auto -cleanup_mmap_fastq(auto &buffer, std::int64_t &buffer_size) { - if (buffer == nullptr) +cleanup_mmap_fastq(auto &data, std::int64_t &length) { + if (data == nullptr) return; - munmap(static_cast(buffer), buffer_size); - buffer = nullptr; - buffer_size = 0; + munmap(static_cast(data), length); + data = nullptr; + length = 0; } auto fastq_file::reset() -> void { - cleanup_mmap_fastq(buffer, buffer_size); + cleanup_mmap_fastq(mmap_data, length); } auto fastq_file::load_next() -> void { // memory mapped data is page aligned but the data we need is not + std::int64_t offset_in_buf = std::distance(mmap_data, std::to_address(last)); static const auto page_mask = sysconf(_SC_PAGESIZE) - 1; - std::tie(start_in_file, cursor) = [&] { - const auto pos_in_file = start_in_file + cursor; + std::tie(start_in_file, offset_in_buf) = [&] { + const auto pos_in_file = start_in_file + offset_in_buf; return std::tuple(pos_in_file & (~page_mask), pos_in_file & page_mask); }(); - stop_in_file = std::min(filesize, start_in_file + target_buffer_size); - if (buffer_size > 0) - cleanup_mmap_fastq(buffer, buffer_size); + cleanup_mmap_fastq(mmap_data, length); + stop_in_file = std::min(filesize, start_in_file + target_length); if (start_in_file < stop_in_file) { // this exist for empty files - buffer_size = stop_in_file - start_in_file; - mmap_fastq(fd, start_in_file, buffer_size, buffer); + length = stop_in_file - start_in_file; + mmap_fastq(fd, start_in_file, length, mmap_data); + buffer = std::span(mmap_data + offset_in_buf, length - offset_in_buf); } } @@ -115,40 +128,39 @@ fastq_file::get_chunks(const std::int64_t n_chunks, const std::int32_t file_id, task_queue &tq, std::atomic_int32_t &n_tasks) -> void { static constexpr auto rec_lines = 4; // FASTQ assert(n_chunks > 0); - const auto buf = std::span(buffer, buffer_size); + const auto beg_itr = std::begin(buffer); + const auto end_itr = std::cend(buffer); // clang-format off const auto not_read_start = [&](const auto p) { // ADS: could get confused if '+' lines have full name info - return buf[p] != '@' || (p > 0 && buf[p - 1] != '\n') || - (p > 2 && buf[p - 2] == '+' && buf[p - 3] == '\n'); + return *p != '@' || (p > beg_itr && *(p - 1) != '\n') || + (p > beg_itr + 2 && *(p - 2) == '+' && *(p - 3) == '\n'); }; - const auto fwd_to_read_start = [&](auto pos) { - if (pos == 0) return pos; - while (pos < buffer_size && not_read_start(pos)) ++pos; + const auto forward_to_start = [&](auto pos) { + if (pos == beg_itr) return pos; + while (pos < end_itr && not_read_start(pos)) ++pos; return pos; }; - const auto rev_to_read_start = [&](auto pos) { - while (pos > 0 && (pos == buffer_size || not_read_start(pos))) --pos; + const auto reverse_to_start = [&](auto pos) { + while (pos > beg_itr && (pos == end_itr || not_read_start(pos))) --pos; return pos; }; // clang-format on - const auto n_bytes_available = buffer_size - cursor; - std::int64_t start_pos = cursor; - const auto [chunk_size, remainder] = std::div(n_bytes_available, n_chunks); - std::vector> chunks(n_chunks); - std::int64_t chunk_end{start_pos}; + const auto [chunk_size, remainder] = std::div(std::ssize(buffer), n_chunks); + auto start_pos = beg_itr; + auto chunk_end = start_pos; for (const auto chunk_idx : std::views::iota(0, n_chunks)) { - const auto chunk_beg = fwd_to_read_start(start_pos); + const auto chunk_beg = forward_to_start(start_pos); const auto stop_pos = start_pos + chunk_size + (chunk_idx < remainder); - chunk_end = fwd_to_read_start(stop_pos); + chunk_end = forward_to_start(stop_pos); if (chunk_idx + 1 == n_chunks) - if (const auto prev = rev_to_read_start(chunk_end); - std::count(std::cbegin(buf) + prev, std::cend(buf), '\n') < rec_lines) + if (const auto prev = reverse_to_start(chunk_end); + std::count(prev, end_itr, '\n') < rec_lines) chunk_end = prev; ++n_tasks; - // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic) - tq.push(file_id, fq_task_t(buffer + chunk_beg, buffer + chunk_end)); + tq.push(file_id, + fq_task_t(std::to_address(chunk_beg), std::to_address(chunk_end))); start_pos = stop_pos; } - cursor = chunk_end; + last = chunk_end; } diff --git a/src/fastq_file.hpp b/src/fastq_file.hpp index 624bd06..678b82b 100644 --- a/src/fastq_file.hpp +++ b/src/fastq_file.hpp @@ -12,6 +12,7 @@ #include #include #include +#include #include #include #include @@ -21,25 +22,17 @@ struct task_queue; struct fastq_file { static constexpr auto min_buf_size = 65536; - std::int64_t target_buffer_size{}; + std::int64_t target_length{}; std::int64_t filesize{}; - char *buffer{}; - std::int64_t buffer_size{}; + char *mmap_data{}; + std::int64_t length{}; std::int64_t start_in_file{}; std::int64_t stop_in_file{}; - std::int64_t cursor{}; int fd{}; + std::span buffer{}; + std::span::iterator last{}; - fastq_file(const std::string &filename, - const std::int64_t target_buffer_size) : - target_buffer_size{ - std::max(target_buffer_size, 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); - } + fastq_file(const std::string &filename, const std::int64_t target_length); // clang-format off fastq_file(const fastq_file &) = delete; @@ -48,14 +41,15 @@ struct fastq_file { // clang-format on fastq_file(fastq_file &&src) noexcept : - target_buffer_size{src.target_buffer_size}, // - filesize{src.filesize}, // - buffer{src.buffer}, // - buffer_size{src.buffer_size}, // - start_in_file{src.start_in_file}, // - stop_in_file{src.stop_in_file}, // - cursor{src.cursor}, // - fd{dup(src.fd)} // <- LOOK + target_length{src.target_length}, // + filesize{src.filesize}, // + mmap_data{src.mmap_data}, // + length{src.length}, // + start_in_file{src.start_in_file}, // + stop_in_file{src.stop_in_file}, // + fd{dup(src.fd)}, // <- LOOK + buffer{src.buffer}, // + last{src.last} // {} auto From 46bbaa43a0459622852a586b378466550d341105 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Fri, 4 Sep 2026 19:43:31 -0700 Subject: [PATCH 14/25] src/tile_processor.hpp: added upper limit on the number of tiles that can be identified, following FastQC, mostly as a precaution in case the user mis-specifies the tile id position within reads when data is read from stdin --- src/tile_processor.hpp | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/src/tile_processor.hpp b/src/tile_processor.hpp index 839144d..1369b0e 100644 --- a/src/tile_processor.hpp +++ b/src/tile_processor.hpp @@ -23,7 +23,7 @@ struct file_info; // Notes // -// - The number of tiles should never exceed 2500? +// - The number of tiles should not exceed 2500 // - Among the first 10k reads, all should contribute to tiles? class tile_processor { @@ -33,6 +33,7 @@ class tile_processor { using qual_vec = std::vector>; using tile_qual_map_t = boost::unordered_flat_map; static constexpr auto read_skip = 10 - 1; + static constexpr auto max_allowed_tiles = 2500; // from FastQC std::uint32_t tile_id_position{}; std::int32_t read_idx{}; @@ -61,8 +62,8 @@ class tile_processor { if (read_idx-- == 0) [[unlikely]] { read_idx = read_skip; update_tile_id(get_name(rec), get_name_end(rec)); - const auto curr_len = static_cast(get_seq_size(rec)); - if (curr_len > max_read_len) + if (const auto curr_len = static_cast(get_seq_size(rec)); + curr_len > max_read_len) resize(curr_len); count_quals_itr(get_qual(rec), get_qual_end(rec), qual); } @@ -99,6 +100,8 @@ class tile_processor { auto update_tile_id(const auto name_beg, const auto name_end) { + static constexpr auto too_many_tiles_msg = + R"(too many tiles observed; likely misidentified tile id position)"; auto tile_itr = name_beg; auto colon_count = 0u; while (colon_count < tile_id_position && tile_itr != name_end) @@ -107,15 +110,16 @@ class tile_processor { const auto [_, ec] = std::from_chars( std::to_address(tile_itr), std::to_address(name_end), curr_tile_id); if (ec != std::errc{}) - throw std::system_error(std::make_error_code(ec), - "failed to parse tile id"); + throw std::system_error(std::make_error_code(ec), "parsing tile id"); if (curr_tile_id != tile_id) { tile_id = curr_tile_id; - auto tile_id_itr = quals.find(tile_id); - if (tile_id_itr == std::cend(quals)) - tile_id_itr = - quals.emplace(tile_id, qual_vec(max_read_len, {0, 0})).first; - qual = std::begin(tile_id_itr->second); + auto id_itr = quals.find(tile_id); + if (id_itr == std::cend(quals)) { + id_itr = quals.emplace(tile_id, qual_vec(max_read_len, {0, 0})).first; + if (std::size(quals) > max_allowed_tiles) + throw std::runtime_error(too_many_tiles_msg); + } + qual = std::begin(id_itr->second); } } }; From 85c78429527a91b5d4d627d3437310ae794e3f0e Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Fri, 4 Sep 2026 20:27:30 -0700 Subject: [PATCH 15/25] static analysis --- src/falco.cpp | 3 ++- src/falco_file_format.hpp | 2 ++ src/fastq_file.cpp | 13 +++++++++---- src/fastq_file.hpp | 6 ------ src/fastq_stdin.cpp | 5 +---- src/fastq_stdin.hpp | 5 +---- src/sam_file.cpp | 1 + src/sam_file.hpp | 2 +- src/sam_stdin.cpp | 19 ++++++++++--------- src/sam_stdin.hpp | 9 ++++----- src/tile_processor.hpp | 1 + 11 files changed, 32 insertions(+), 34 deletions(-) diff --git a/src/falco.cpp b/src/falco.cpp index 6bbf75d..9389bf6 100644 --- a/src/falco.cpp +++ b/src/falco.cpp @@ -370,7 +370,7 @@ main(int argc, char *argv[]) { ->transform(size_from_units); app.add_option_function>( "--stdin", - [&](const auto &arg) { // callback is to allow trailing arg to be name + [&](const auto &arg) { // callback is to allow trailing arg to be name stdin_info = arg; infiles_opt->get_validator("file_check")->active(false); infiles_opt->expected(1); @@ -382,6 +382,7 @@ main(int argc, char *argv[]) { ->type_size(1, 2) ->transform(CLI::CheckedTransformer(format_name_map, CLI::ignore_case) .application_index(0)) + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers) ->check(CLI::IsMember({4, 6}).application_index(1)) ->callback_priority(CLI::CallbackPriority::PreRequirementsCheck); app.add_flag("--bisulfite", do_bisulfite, diff --git a/src/falco_file_format.hpp b/src/falco_file_format.hpp index 8ff6118..b8c0c64 100644 --- a/src/falco_file_format.hpp +++ b/src/falco_file_format.hpp @@ -5,8 +5,10 @@ #include "nlohmann/json.hpp" +#include #include #include +#include #include #include #include diff --git a/src/fastq_file.cpp b/src/fastq_file.cpp index be447ea..add9e7f 100644 --- a/src/fastq_file.cpp +++ b/src/fastq_file.cpp @@ -10,17 +10,21 @@ #include #include -#include // IWYU pragma: keep +#include #include +#include +#include #include +#include #include +#include +#include +#include #include #include #include #include -#include // IWYU pragma: keep -#include -#include +#include #include [[nodiscard]] auto @@ -119,6 +123,7 @@ fastq_file::load_next() -> void { if (start_in_file < stop_in_file) { // this exist for empty files length = stop_in_file - start_in_file; mmap_fastq(fd, start_in_file, length, mmap_data); + // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic) buffer = std::span(mmap_data + offset_in_buf, length - offset_in_buf); } } diff --git a/src/fastq_file.hpp b/src/fastq_file.hpp index 678b82b..4a8a126 100644 --- a/src/fastq_file.hpp +++ b/src/fastq_file.hpp @@ -3,18 +3,12 @@ #ifndef SRC_FASTQ_FILE_HPP_ #define SRC_FASTQ_FILE_HPP_ -#include #include -#include #include -#include #include -#include -#include #include #include -#include #include #include diff --git a/src/fastq_stdin.cpp b/src/fastq_stdin.cpp index 7e92227..946d9aa 100644 --- a/src/fastq_stdin.cpp +++ b/src/fastq_stdin.cpp @@ -1,7 +1,6 @@ // SPDX-License-Identifier: MIT; Copyright 2026 Andrew D Smith #include "fastq_stdin.hpp" -#include "falco_utils.hpp" #include "fqrec.hpp" #include "task_queue.hpp" @@ -11,15 +10,13 @@ #include #include #include -#include #include #include #include -#include +#include #include #include #include -#include #include fastq_stdin::fastq_stdin(const std::int64_t buf_size) : diff --git a/src/fastq_stdin.hpp b/src/fastq_stdin.hpp index 2eac214..cdf5a85 100644 --- a/src/fastq_stdin.hpp +++ b/src/fastq_stdin.hpp @@ -6,12 +6,9 @@ #include #include #include -#include #include -#include #include #include -#include #include struct task_queue; @@ -23,7 +20,7 @@ struct fastq_stdin { std::vector::iterator last; bool hit_eof{}; - fastq_stdin(const std::int64_t buf_size); + explicit fastq_stdin(const std::int64_t buf_size); operator bool() const { return cursor < last || !hit_eof; } // clang-format off diff --git a/src/sam_file.cpp b/src/sam_file.cpp index 276d1ab..4624c5c 100644 --- a/src/sam_file.cpp +++ b/src/sam_file.cpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include diff --git a/src/sam_file.hpp b/src/sam_file.hpp index 1cf562f..7832598 100644 --- a/src/sam_file.hpp +++ b/src/sam_file.hpp @@ -8,7 +8,7 @@ #include #include #include -#include +#include // IWYU pragma: keep #include #include diff --git a/src/sam_stdin.cpp b/src/sam_stdin.cpp index 56424cd..f9465c5 100644 --- a/src/sam_stdin.cpp +++ b/src/sam_stdin.cpp @@ -4,23 +4,25 @@ #include "samrec.hpp" #include "task_queue.hpp" +#include + #include #include #include #include #include -#include +#include #include +#include #include -#include #include #include +#include // IWYU pragma: keep sam_stdin::sam_stdin(const std::int64_t buf_size) : buffer(buf_size + min_buf_size), cursor{std::begin(buffer)}, last{std::begin(buffer)} { - if (!skip_header()) - std::runtime_error("failed to validate SAM header: stdin"); + skip_header(); } [[nodiscard]] auto @@ -32,8 +34,8 @@ estimate_n_reads_sam_stdin(const std::string &) return {assumed_n_reads, assumed_read_len, assumed_filesize}; } -[[nodiscard]] auto -sam_stdin::skip_header() -> bool { +auto +sam_stdin::skip_header() -> void { bool pre_header = true; while (cursor == std::cbegin(buffer)) { load_next(); @@ -51,7 +53,6 @@ sam_stdin::skip_header() -> bool { shift_output_buffer(); // in prep for subsequent load-next } } - return true; } auto @@ -110,8 +111,8 @@ sam_stdin::load_next() -> void { std::int64_t n{1}; while (space > 0 && (n = read(0, std::to_address(last), space)) != 0) { if (n == -1) - std::system_error(std::make_error_code(std::errc(errno)), - "error reading fastq from stdin"); + throw std::system_error(std::make_error_code(std::errc(errno)), + "error reading fastq from stdin"); space -= n; last += n; } diff --git a/src/sam_stdin.hpp b/src/sam_stdin.hpp index 9566f15..1c84316 100644 --- a/src/sam_stdin.hpp +++ b/src/sam_stdin.hpp @@ -6,10 +6,9 @@ #include #include #include -#include #include -#include #include +#include #include struct task_queue; @@ -22,7 +21,7 @@ class sam_stdin { bool hit_eof{}; public: - sam_stdin(const std::int64_t buf_size); + explicit sam_stdin(const std::int64_t buf_size); operator bool() const { return cursor < last || !hit_eof; } // clang-format off @@ -56,8 +55,8 @@ class sam_stdin { auto load_next() -> void; - [[nodiscard]] auto - skip_header() -> bool; + auto + skip_header() -> void; }; [[nodiscard]] auto diff --git a/src/tile_processor.hpp b/src/tile_processor.hpp index 1369b0e..ac768a6 100644 --- a/src/tile_processor.hpp +++ b/src/tile_processor.hpp @@ -14,6 +14,7 @@ #include #include #include +#include #include #include #include From 5dbf574d05cd1eef4ae24ab08d175d7a2615100b Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sat, 5 Sep 2026 13:16:02 -0700 Subject: [PATCH 16/25] src/fastq_stdin.cpp: fixing wrong algorithm for validating fastq format from stdin --- src/fastq_stdin.cpp | 28 ++++++++++++++++------------ 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/src/fastq_stdin.cpp b/src/fastq_stdin.cpp index 946d9aa..5c35d71 100644 --- a/src/fastq_stdin.cpp +++ b/src/fastq_stdin.cpp @@ -86,19 +86,23 @@ fastq_stdin::shift_output_buffer() -> void { [[nodiscard]] static auto validate_fastq(const auto &buffer) { static constexpr auto n_bytes_to_validate = 16 * 1024; - // verify that no two consecutive newlines are followed by a '@' - bool prev_was_ampersand{false}; + static constexpr auto name_line_symbol = '@'; + static constexpr auto plus_line_symbol = '+'; + static constexpr auto lines_per_rec = 4; + static constexpr auto name_line = 0; + static constexpr auto plus_line = 2; + bool prev_nl{true}; auto n_bytes = 0; - for (auto itr = std::cbegin(buffer); itr + 1 != std::cend(buffer); ++itr) { - if (*itr == '\n') { - if (*(itr + 1) == '@') { - if (prev_was_ampersand) - return false; - prev_was_ampersand = true; - } - else - prev_was_ampersand = false; - } + auto n_lines = 0; + for (const auto c : buffer) { + if (prev_nl && (n_lines % lines_per_rec == name_line) && + c != name_line_symbol) + return false; + if (prev_nl && (n_lines % lines_per_rec == plus_line) && + c != plus_line_symbol) + return false; + prev_nl = (c == '\n'); + n_lines += prev_nl; if (++n_bytes == n_bytes_to_validate) return true; } From b0921a740bf803982c221f32f8fa0828d2ea0d11 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sat, 5 Sep 2026 13:21:35 -0700 Subject: [PATCH 17/25] cmake/tests.cmake: adding tests for data coming from stdin --- cmake/tests.cmake | 2 ++ 1 file changed, 2 insertions(+) diff --git a/cmake/tests.cmake b/cmake/tests.cmake index 59d3076..b2ea0f9 100644 --- a/cmake/tests.cmake +++ b/cmake/tests.cmake @@ -34,3 +34,5 @@ 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) 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) From ae998e251162d2f9d09ed6a444b6115eb57a443c Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sat, 5 Sep 2026 13:22:34 -0700 Subject: [PATCH 18/25] data/test_data/md5sum.txt: adding hashes for tests with data from stdin --- data/test_data/md5sum.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/data/test_data/md5sum.txt b/data/test_data/md5sum.txt index 22053d6..d90b2a1 100644 --- a/data/test_data/md5sum.txt +++ b/data/test_data/md5sum.txt @@ -13,3 +13,5 @@ bc8fd6e40a0ca55cb00634bc306a896c groups_out/bam_1/fastqc_data.txt a1952ae366bd7c7207f40db833b1d16b sam_out/sam_1/fastqc_data.txt 3fe18e2ee85e3912ede8c57559bf9f88 preseq_out/bam_1/preseq_hist.txt 045639456fbe4a81b5de2789ceb95ddd bam_kmers_out/bam_1/fastqc_data.txt +a23fa58eafe340d2a9949374187da4bc fastq_stdin_out/fastq_stdin/fastqc_data.txt +d69f1b32d3e09772667c48031ee1a915 sam_stdin_out/sam_stdin/fastqc_data.txt From 29bb45be4fdbfe2c60fa88c5542100624445a977 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sat, 5 Sep 2026 13:24:19 -0700 Subject: [PATCH 19/25] data/test_scripts/fastq_stdin.sh data/test_scripts/sam_stdin.sh: adding tests for stdin --- data/test_scripts/fastq_stdin.sh | 24 ++++++++++++++++++++++++ data/test_scripts/sam_stdin.sh | 24 ++++++++++++++++++++++++ 2 files changed, 48 insertions(+) create mode 100644 data/test_scripts/fastq_stdin.sh create mode 100644 data/test_scripts/sam_stdin.sh diff --git a/data/test_scripts/fastq_stdin.sh b/data/test_scripts/fastq_stdin.sh new file mode 100644 index 0000000..6ddda81 --- /dev/null +++ b/data/test_scripts/fastq_stdin.sh @@ -0,0 +1,24 @@ +#!/usr/bin/env bash +# SPDX-License-Identifier: MIT + +# Test for fastq input from stdin + +prog=./falco +infile1=test_data/fastq_bgzip_1.fq.gz +name=fastq_stdin +outdir=fastq_stdin_out +if [[ -e "${infile1}" ]]; then + mkdir -p ${outdir} + gunzip -c ${infile1} | \ + ${prog} -o ${outdir} --stdin fq ${name} + 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/data/test_scripts/sam_stdin.sh b/data/test_scripts/sam_stdin.sh new file mode 100644 index 0000000..30f1c78 --- /dev/null +++ b/data/test_scripts/sam_stdin.sh @@ -0,0 +1,24 @@ +#!/usr/bin/env bash +# SPDX-License-Identifier: MIT + +# Test for SAM input from stdin + +prog=./falco +infile1=test_data/bam_1.bam +name=sam_stdin +outdir=sam_stdin_out +if [[ -e "${infile1}" ]]; then + mkdir -p ${outdir} + samtools view -h ${infile1} | \ + ${prog} -o ${outdir} --stdin sam ${name} + 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 95dac0c193c0cd8e339b5a31773d6622ede60275 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sat, 5 Sep 2026 13:44:55 -0700 Subject: [PATCH 20/25] src/fastq_file.cpp: fixing std:div to std::ldiv because macos complains --- src/fastq_file.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/fastq_file.cpp b/src/fastq_file.cpp index add9e7f..ba9b1e2 100644 --- a/src/fastq_file.cpp +++ b/src/fastq_file.cpp @@ -47,7 +47,7 @@ estimate_n_reads_fastq(const std::string &filename) return {{}, {}, filesize}; const auto [part_size, remainder] = (filesize < n_parts * max_part_size) - ? std::div(filesize, n_parts) + ? std::ldiv(filesize, n_parts) : std::ldiv_t{max_part_size, 0}; auto n_lines = 0LU; auto readlen_est = 0LU; @@ -151,7 +151,7 @@ fastq_file::get_chunks(const std::int64_t n_chunks, const std::int32_t file_id, return pos; }; // clang-format on - const auto [chunk_size, remainder] = std::div(std::ssize(buffer), n_chunks); + const auto [chunk_size, remainder] = std::ldiv(std::ssize(buffer), n_chunks); auto start_pos = beg_itr; auto chunk_end = start_pos; for (const auto chunk_idx : std::views::iota(0, n_chunks)) { From 6ba878e59afc51bb6f8938844ace4f5c43e9a31d Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sat, 5 Sep 2026 17:01:16 -0700 Subject: [PATCH 21/25] src/report.cpp: fixing constant array index --- src/report.cpp | 30 ++++++++++++++---------------- 1 file changed, 14 insertions(+), 16 deletions(-) diff --git a/src/report.cpp b/src/report.cpp index 0dc8549..9f576a5 100644 --- a/src/report.cpp +++ b/src/report.cpp @@ -118,19 +118,17 @@ quality_sequence_report(const falco::qual_array &qual_by_read, r += header; // output quality values between first and last non-zero const auto gt0 = [&](const auto x) { return x > 0; }; - const auto begin_obs_itr = std::ranges::find_if(qual_by_read, gt0); - if (begin_obs_itr == std::cend(qual_by_read)) + const auto qbr = std::span(qual_by_read); + const auto beg_obs_itr = std::ranges::find_if(qbr, gt0); + const auto end_obs_itr = std::cbegin(std::ranges::find_last_if(qbr, gt0)); + if (beg_obs_itr == std::cend(qbr) && end_obs_itr == std::cend(qbr)) throw std::runtime_error("error finding quality scores generating report"); - const std::int64_t begin_obs = - std::distance(std::cbegin(qual_by_read), begin_obs_itr); - const auto end_obs_subrange = std::ranges::find_last_if(qual_by_read, gt0); - const std::int64_t end_obs = - std::ssize(qual_by_read) - std::ssize(end_obs_subrange) + 1; - assert(begin_obs >= 0 && end_obs <= falco::max_qual_val); - for (const auto q : std::views::iota(begin_obs, end_obs)) - // cppcheck-suppress useStlAlgorithm - // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-constant-array-index) - r += std::format("{}\t{}\n", q, qual_by_read[q]); + const auto beg_idx = std::distance(std::cbegin(qbr), beg_obs_itr); + const auto end_idx = std::distance(std::cbegin(qbr), end_obs_itr); + assert(end_idx <= falco::max_qual_val); + std::ranges::for_each( + std::views::iota(beg_idx, end_idx + 1), + [&](const auto q) { r += std::format("{}\t{}\n", q, qbr[q]); }); r += end_module_tag; return r; } @@ -208,8 +206,8 @@ basic_stats_report(const file_info &info, const std::uint64_t n_reads, [[nodiscard]] auto tile_report(const tile_processor::tiles_centered_t ¢ered, - const std::vector &groups, const file_grades &grades) - -> std::string { + const std::vector &groups, + const file_grades &grades) -> std::string { static constexpr auto label = "tile"; static constexpr auto max_precision{std::numeric_limits::digits10}; static constexpr auto start_tag = ">>Per tile sequence quality\t{}\n"; @@ -228,8 +226,8 @@ tile_report(const tile_processor::tiles_centered_t ¢ered, } [[nodiscard]] auto -kmer_report(const std::vector &results, const file_grades &grades) - -> std::string { +kmer_report(const std::vector &results, + const file_grades &grades) -> std::string { static constexpr auto label = "kmer"; static constexpr auto start_tag = ">>Kmer Content\t{}\n"; static constexpr auto header = "#Sequence\t" From ef471bd39e17ee9fcbd9c1702c24da0b026cc00a Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sat, 5 Sep 2026 17:21:41 -0700 Subject: [PATCH 22/25] src/fastq_stdin.hpp: fixing wrong function name in declaration --- src/fastq_stdin.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fastq_stdin.hpp b/src/fastq_stdin.hpp index cdf5a85..6cf2fa8 100644 --- a/src/fastq_stdin.hpp +++ b/src/fastq_stdin.hpp @@ -56,7 +56,7 @@ struct fastq_stdin { }; [[nodiscard]] auto -estimate_n_reads_fastq(const std::string &filename) +estimate_n_reads_fastq_stdin(const std::string &filename) -> std::tuple; inline auto From 6cdd2b22824ae7fa191c118918adec91993c7b9d Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sat, 5 Sep 2026 17:23:32 -0700 Subject: [PATCH 23/25] src/duplication_results.hcpp: added a minimum absolute count for the overrep cutoff so that on very small data and stdin the overrep wont trigger at 1 observation --- src/duplication_results.cpp | 4 +++- src/duplication_results.hpp | 3 ++- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/src/duplication_results.cpp b/src/duplication_results.cpp index 2e3dd17..05e79d8 100644 --- a/src/duplication_results.cpp +++ b/src/duplication_results.cpp @@ -109,7 +109,9 @@ duplication_results::get_preseq_hist() const [[nodiscard]] auto duplication_results::get_overrepresented(const std::uint64_t n_reads) const -> std::vector { - const auto cutoff = static_cast(n_reads) * overrep_cutoff; + const auto cutoff = + std::max(static_cast(n_reads) * overrep_frac_cutoff, + static_cast(overrep_count_cutoff)); const auto gte_cutoff = [&](const auto p) { return p.second >= cutoff; }; const auto rev_p = [&](const auto p) { return std::pair{p.second, p.first}; }; auto overrep = dups | std::views::filter(gte_cutoff) | diff --git a/src/duplication_results.hpp b/src/duplication_results.hpp index 75e3f55..322b583 100644 --- a/src/duplication_results.hpp +++ b/src/duplication_results.hpp @@ -53,7 +53,8 @@ struct duplication_results { static constexpr auto max_n_reads_total{1'000'000}; static constexpr auto max_reads_to_hash{100'000}; static constexpr auto default_read_skip{10}; - static constexpr auto overrep_cutoff = 0.001; + static constexpr auto overrep_frac_cutoff = 0.001; + static constexpr auto overrep_count_cutoff = 5; std::int64_t count_at_limit{}; std::int64_t read_skip{default_read_skip}; From 29ed2e7f7b6f491dbac34cc77ad6f7e6d7e2cc77 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sat, 5 Sep 2026 17:28:41 -0700 Subject: [PATCH 24/25] data/test_data/md5sum.txt: updating test hashes for sam and fastq from stdin because of changes to output for dups and overrep --- data/test_data/md5sum.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/data/test_data/md5sum.txt b/data/test_data/md5sum.txt index d90b2a1..dd86140 100644 --- a/data/test_data/md5sum.txt +++ b/data/test_data/md5sum.txt @@ -13,5 +13,5 @@ bc8fd6e40a0ca55cb00634bc306a896c groups_out/bam_1/fastqc_data.txt a1952ae366bd7c7207f40db833b1d16b sam_out/sam_1/fastqc_data.txt 3fe18e2ee85e3912ede8c57559bf9f88 preseq_out/bam_1/preseq_hist.txt 045639456fbe4a81b5de2789ceb95ddd bam_kmers_out/bam_1/fastqc_data.txt -a23fa58eafe340d2a9949374187da4bc fastq_stdin_out/fastq_stdin/fastqc_data.txt -d69f1b32d3e09772667c48031ee1a915 sam_stdin_out/sam_stdin/fastqc_data.txt +a1aa02211d63f507646a9e12122f0226 sam_stdin_out/sam_stdin/fastqc_data.txt +c61e37d2c8d64b6950a5dad300e5fd60 fastq_stdin_out/fastq_stdin/fastqc_data.txt From b3dd4b54ea1eb7cf3e82d018fadff04022e1cab7 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sat, 5 Sep 2026 17:50:33 -0700 Subject: [PATCH 25/25] src/falco.cpp: adding cli logic for reading sam or fastq from stdin, including handling tiles --- src/falco.cpp | 25 ++++++++++++++++++++----- 1 file changed, 20 insertions(+), 5 deletions(-) diff --git a/src/falco.cpp b/src/falco.cpp index 9389bf6..420d3a5 100644 --- a/src/falco.cpp +++ b/src/falco.cpp @@ -82,7 +82,6 @@ Use these as templates. Copy and modify them to customize your analysis. #include #include #include -#include #include #include #include @@ -172,6 +171,10 @@ make_reads_file_stdin(const std::vector &infos, get_file_info_stdin(const std::vector &names, const std::pair &ft_tile) -> std::vector { + const auto [n_reads_est, read_len_est, _] = + ft_tile.first == falco::file_format::fastq + ? estimate_n_reads_fastq_stdin(names.front()) + : estimate_n_reads_sam_stdin(names.front()); file_info info; info.name = names.front(); info.format = ft_tile.first; @@ -179,6 +182,8 @@ get_file_info_stdin(const std::vector &names, info.size = 0; info.has_tiles = (ft_tile.second != 0); info.tile_id_position = ft_tile.second; + info.n_reads_est = n_reads_est; + info.read_len_est = read_len_est; return std::vector(1, info); } @@ -267,6 +272,7 @@ main(int argc, char *argv[]) { static constexpr auto buffer_size_default = 256 * 1024 * 1024; static constexpr std::int64_t min_buf_size = 1024 * 1024; + static constexpr std::int64_t max_buf_size = 1024L * 1024L * 1024L * 1024L; std::vector infiles; std::string contam_file; std::string config_file; @@ -327,8 +333,9 @@ main(int argc, char *argv[]) { app.get_formatter()->long_option_alignment_ratio(0.2); app.set_help_flag("-h,--help", "Print more detailed help"); app.set_version_flag("--version", VERSION, "Print program version"); - app.add_flag("--license", [&](auto) { std::print("{}", license_text); throw CLI::Success(); }, - "Print full license") + app.add_flag("--license", [&](auto) { + std::print("{}", license_text); throw CLI::Success(); }, + "Print full license") ->callback_priority(CLI::CallbackPriority::PreRequirementsCheck); auto infiles_opt = app.add_option("INFILES", infiles, @@ -345,7 +352,7 @@ main(int argc, char *argv[]) { ->option_text(std::format("[{}]", n_threads)); app.add_option("-m,--mem", buffer_size, "Input memory buffer size (G/M/K units ok)") - ->check(CLI::Range(min_buf_size, std::numeric_limits::max())) + ->check(CLI::Range(min_buf_size, max_buf_size)) ->option_text(std::format("[{}]", size_to_units(buffer_size_default))) ->capture_default_str() ->transform(size_from_units); @@ -423,6 +430,14 @@ main(int argc, char *argv[]) { CLI11_PARSE(app, argc, argv); const bool do_stdin = stdin_info.first != falco::file_format::unknown; + if (do_stdin) { + const auto deduced_do_tiles = stdin_info.second ? 1 : -1; + if (do_tiles && do_tiles != deduced_do_tiles) { + std::println("inconsistent tile analysis args for data from stdin"); + return EXIT_FAILURE; + } + do_tiles = deduced_do_tiles; + } run_mode mode; // declare mode here so we can assign from config file if (!config_file.empty()) @@ -473,7 +488,7 @@ main(int argc, char *argv[]) { // restrict buffer size to avoid using a possibly harmful amount of memory const auto get_sz = [](const auto &i) { return i.size; }; const auto max_sz = std::ranges::max(std::views::transform(infos, get_sz)); - buffer_size = buffer_size < max_sz ? buffer_size : max_sz; + buffer_size = buffer_size < max_sz ? buffer_size : min_buf_size; const auto min_buffer_size = get_min_buffer_size(max_read_length); if (min_buffer_size > buffer_size) {