diff --git a/README.md b/README.md index b3d3747..24991b9 100644 --- a/README.md +++ b/README.md @@ -41,10 +41,9 @@ conda install -c bioconda falco ## Building -- Falco uses features of C++23. To compile on Linux: GCC >= 14.2.0 or LLVM-Clang >= 20.0.0. - On macOS, GCC >= 15. Unfortunately GCC 16.1 has bugs that impact - falco, and 16.1 is the most recent on macOS through Homebrew. These seem to be - fixed in GCC 16.2, and GCC 15 can also be obtained through Homebrew. +- Falco uses features of C++23. To compile on Linux: + GCC >= 14.2.0 or LLVM-Clang >= 20.0.0. On macOS, GCC >= 15 or GCC 16.2, which + is the most recent (not GCC 16.1, which has a bug). - Falco uses the cmake build system. - Dependencies: * [HTSLib](https://github.com/samtools/htslib): used for identifying file formats. @@ -72,7 +71,6 @@ specified directly. I'm explaining this via a clean Ubuntu instance in docker: ``` -docker pull ubuntu:latest docker run -it ubuntu:latest bash ``` Inside the docker: @@ -80,7 +78,6 @@ Inside the docker: export DEBIAN_FRONTEND=noninteractive && apt-get update && apt-get install -y --no-install-recommends \ - libssl-dev \ zlib1g-dev \ libdeflate-dev \ libisal-dev \ @@ -106,10 +103,10 @@ 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: ``` -brew install libdeflate isa-l htslib samtools && # samtools for testing +brew install libdeflate htslib samtools && # samtools for testing git clone https://github.com/smithlabcode/falco.git && cd falco && -cmake -B build -DUSE_ISAL=on -DCMAKE_CXX_COMPILER=g++-15 -DCMAKE_BUILD_TYPE=Release && +cmake -B build -DCMAKE_CXX_COMPILER=g++-15 -DCMAKE_BUILD_TYPE=Release && cmake --build build -j8 && ctest --test-dir build ``` @@ -117,6 +114,8 @@ ZLib is already installed on macOS, HTSLib installs libdeflate as a dependency and samtools installs both as dependency. To see what's already installed on GitHub's macOS look [here](https://github.com/actions/runner-images/blob/main/images/macos/macos-26-Readme.md). +Note: ISA-L was designed for Intel hardware. It works on Apple silicon, but only +gives 1-2% speedup in my tests. ## Changes in Falco v2.0 @@ -124,13 +123,14 @@ GitHub's macOS look I found that the method for tile analysis is a bit unstable, and the tile grade can be slightly unstable. The only way to notice this is to process reads from -the same input file in different orders. This happens as a side effect of +the same input file in different orders, which happens as a side effect of analyzing reads concurrently with threads. Here is my understanding of how -FastQC works, and how I implemented falco v2.0. Please comment if you see -anything incorrect. +FastQC works, and how I implemented tile analysis in falco v2. Please comment if +you see anything incorrect. - Tile analysis is done for 1/10 of the reads (though FastQC includes all among the first 10k reads). +- If more than 2500 tiles are identified, an error is assumed. - Accumulating results: For each counted read, for each position in the read, the quality score contributes to that tile's mean for the given position. - Summarizing tile results: For each read position, the mean over tiles' quality @@ -146,6 +146,11 @@ runs. I've noticed that this can lead to differences between runs, and in some cases this has changed the grade between pass/warn and warn/fail. So it is possible the grade can differ between runs for the same data. +About input from stdin: If you are using standard input tile analysis is diabled +unless you specify where to find the tile info in the read name, e.g., +`cat file.fq | falco --stdin fq:4 outdir`, where the 4 indicates that the tile is +after the 4th colon in the read name. The only supported positions are 4 and 6. + ### Duplication results I changed how falco evaluates "duplcation". Although the format of the output is diff --git a/src/falco.cpp b/src/falco.cpp index 420d3a5..a9ddea3 100644 --- a/src/falco.cpp +++ b/src/falco.cpp @@ -167,6 +167,15 @@ make_reads_file_stdin(const std::vector &infos, return reads_files; } +[[nodiscard]] static auto +make_reads_files(const run_mode &mode, const std::vector &infos, + const std::vector &infiles, + const std::int64_t buffer_size) -> std::vector { + // ADS: need to do this differently for stdin + return mode.do_stdin() ? make_reads_file_stdin(infos, buffer_size) + : make_reads_files(infos, infiles, buffer_size); +} + [[nodiscard]] static auto get_file_info_stdin(const std::vector &names, const std::pair &ft_tile) @@ -222,12 +231,14 @@ get_file_info(const auto &infiles) { } [[nodiscard]] static auto -make_outdirs(const auto &ins, const auto &outdir) -> std::vector { +make_outdirs(const auto &ins, const auto &outdir, + const bool keep_extn = false) -> std::vector { namespace fs = std::filesystem; fs::create_directory(outdir); const auto compose_dirname = [&](const auto &fname) { - const auto without_path = fs::path{fname}.filename(); - return (fs::path{outdir} / remove_extension(without_path)).string(); + const auto without_path = fs::path{fname}.filename().string(); + const auto name = keep_extn ? without_path : remove_extension(without_path); + return (fs::path{outdir} / name).string(); }; const auto dnames = ins | std::views::transform(compose_dirname); std::ranges::for_each(dnames, [](const auto &d) { fs::create_directory(d); }); @@ -285,6 +296,7 @@ main(int argc, char *argv[]) { int do_dup_analysis{}; int do_adap{}; + int do_stdin{}; int do_groups{}; int do_bisulfite{}; int do_preseq{}; @@ -309,8 +321,8 @@ main(int argc, char *argv[]) { }; // Related to reading data from stdin - const auto format_name_map = std::map{ - std::pair{"fq"s, falco::file_format::fastq}, + const auto format_name_map = std::map{ + {"fq"s, falco::file_format::fastq}, {"sam"s, falco::file_format::sam}, }; auto stdin_info = std::pair{ @@ -378,6 +390,7 @@ main(int argc, char *argv[]) { app.add_option_function>( "--stdin", [&](const auto &arg) { // callback is to allow trailing arg to be name + do_stdin = 1; stdin_info = arg; infiles_opt->get_validator("file_check")->active(false); infiles_opt->expected(1); @@ -429,8 +442,15 @@ main(int argc, char *argv[]) { } CLI11_PARSE(app, argc, argv); - const bool do_stdin = stdin_info.first != falco::file_format::unknown; - if (do_stdin) { + 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); + else + // if no config file, use default graders + grader_set::instance(); + + mode.set_do_stdin(do_stdin); + if (mode.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"); @@ -439,13 +459,6 @@ main(int argc, char *argv[]) { do_tiles = deduced_do_tiles; } - 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); - else - // if no config file, use default graders - grader_set::instance(); - // now set run mode values to take priority over config file if (!adapters_file.empty()) do_adap = 1; @@ -460,7 +473,7 @@ main(int argc, char *argv[]) { mode.set_do_original_dups(do_original_dups); mode.set_unassigned(); - const auto outdirs = make_outdirs(infiles, outdir); + const auto outdirs = make_outdirs(infiles, outdir, mode.do_stdin()); if (!contam_file.empty()) { load_contaminants(contam_file); @@ -482,8 +495,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 = do_stdin ? get_file_info_stdin(infiles, stdin_info) - : get_file_info(infiles); + auto infos = mode.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; }; @@ -522,8 +535,9 @@ 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 reads_files = make_reads_files(mode, infos, infiles, buffer_size); + auto results = analyze(n_threads, mode, infos, std::move(reads_files), std::move(dups)); diff --git a/src/fastq_file.cpp b/src/fastq_file.cpp index ba9b1e2..c1fa3a7 100644 --- a/src/fastq_file.cpp +++ b/src/fastq_file.cpp @@ -163,8 +163,7 @@ fastq_file::get_chunks(const std::int64_t n_chunks, const std::int32_t file_id, std::count(prev, end_itr, '\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))); + tq.push(file_id, fq_task_t(chunk_beg, chunk_end)); start_pos = stop_pos; } last = chunk_end; diff --git a/src/fastq_stdin.cpp b/src/fastq_stdin.cpp index 5c35d71..6c3c8fb 100644 --- a/src/fastq_stdin.cpp +++ b/src/fastq_stdin.cpp @@ -68,8 +68,7 @@ fastq_stdin::get_chunks(const std::int64_t n_chunks, const std::int32_t file_id, 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))); + tq.push(file_id, fq_task_t(chunk_beg, chunk_end)); start_itr = stop_itr; } cursor = chunk_end; @@ -84,7 +83,7 @@ fastq_stdin::shift_output_buffer() -> void { } [[nodiscard]] static auto -validate_fastq(const auto &buffer) { +validate_fastq(const auto &buffer) -> bool { static constexpr auto n_bytes_to_validate = 16 * 1024; static constexpr auto name_line_symbol = '@'; static constexpr auto plus_line_symbol = '+'; @@ -121,8 +120,6 @@ 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 @@ -133,5 +130,7 @@ fastq_stdin::make_tasks(const std::int64_t n_chunks, // n_tasks = 1; // for current task, which makes more tasks shift_output_buffer(); load_next(); + if (!validate_fastq(buffer)) + throw std::runtime_error("input appears not to be FASTQ"); get_chunks(n_chunks, file_id, tq, n_tasks); } diff --git a/src/fqrec.hpp b/src/fqrec.hpp index 8bf0357..a648a28 100644 --- a/src/fqrec.hpp +++ b/src/fqrec.hpp @@ -51,6 +51,10 @@ struct fqrec { struct fq_task_t { fqrec::pos_t beg{}; fqrec::pos_t end{}; + fq_task_t(const fqrec::pos_t beg, const fqrec::pos_t end) : + 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]] inline auto diff --git a/src/run_mode.cpp b/src/run_mode.cpp index 5993cba..1713ad8 100644 --- a/src/run_mode.cpp +++ b/src/run_mode.cpp @@ -44,6 +44,11 @@ // do_original_dups (use the original duplication mode from FastQC and Falco v1) // do_preseq (make another output file for input to preseq) // +// ADS (2026-09-08 edit): added do_stdin because otherwise it's a floating +// almost global variable. +// +// do_stdin (multiple effects; allows input of FASTQ or SAM from stdin) +// // clang-format on // clang-format off @@ -64,6 +69,7 @@ std::vector run_mode::labels{ // NOLINT(cert-err58-cpp) auto run_mode::assign(const std::unordered_map &modes) -> void { + // ADS: only assigns settings that could be from the config file static const auto set_mode = [&](const std::string &label, auto &the_mode) { const auto itr = modes.find(label); if (itr != std::cend(modes)) @@ -86,6 +92,7 @@ auto run_mode::set_unassigned() -> void { // clang-format off // settings below are not in config file + if (do_stdin_ == 0) do_stdin_ = do_stdin_default; if (do_groups_ == 0) do_groups_ = do_groups_default; if (do_bisulfite_ == 0) do_bisulfite_ = do_bisulfite_default; if (do_preseq_ == 0) do_preseq_ = do_preseq_default; diff --git a/src/run_mode.hpp b/src/run_mode.hpp index ec63883..4e838a1 100644 --- a/src/run_mode.hpp +++ b/src/run_mode.hpp @@ -17,6 +17,7 @@ class run_mode { // clang-format off // ADS: these first params are not set in config file + [[nodiscard]] auto do_stdin() const -> bool { return do_stdin_ == 1; } [[nodiscard]] auto do_groups() const -> bool { return do_groups_ == 1; } [[nodiscard]] auto do_bisulfite() const -> bool { return do_bisulfite_ == 1; } [[nodiscard]] auto do_preseq() const -> bool { return do_preseq_ == 1; } @@ -40,6 +41,7 @@ class run_mode { // clang-format on // clang-format off + auto set_do_stdin(const int x) { if (x) do_stdin_ = x; } auto set_do_groups(const int x) { if (x) do_groups_ = x; } auto set_do_bisulfite(const int x) { if (x) do_bisulfite_ = x; } auto set_do_preseq(const int x) { if (x) do_preseq_ = x; } @@ -67,6 +69,7 @@ class run_mode { private: // ADS: 1 is yes; -1 is no; 0 is not assigned // first settings are not in config file + static constexpr auto do_stdin_default = -1; // OFF static constexpr auto do_groups_default = -1; // OFF static constexpr auto do_bisulfite_default = -1; // OFF static constexpr auto do_preseq_default = -1; // OFF @@ -85,6 +88,7 @@ class run_mode { static constexpr auto do_tiles_default = 1; // affects processing // not in config file + int do_stdin_{}; int do_groups_{}; int do_bisulfite_{}; int do_preseq_{}; diff --git a/src/sam_stdin.cpp b/src/sam_stdin.cpp index f9465c5..aa8c8e1 100644 --- a/src/sam_stdin.cpp +++ b/src/sam_stdin.cpp @@ -15,6 +15,7 @@ #include #include #include +#include #include #include #include // IWYU pragma: keep @@ -119,6 +120,20 @@ sam_stdin::load_next() -> void { hit_eof = (n == 0); } +[[nodiscard]] static auto +validate_sam(const auto &buffer) -> bool { + // SAM format has 11+ fields, tab separated and the docs give a regex for each + // field. We are only checking the first char of the first field. + static constexpr auto n_bytes_to_validate = 16L * 1024; + assert(std::size(buffer)); + const auto n_bytes = std::min(n_bytes_to_validate, std::ssize(buffer)); + const auto buf_end = std::cbegin(buffer) + n_bytes - (n_bytes >= 1L); + for (auto itr = std::cbegin(buffer); itr != buf_end; ++itr) + if (*itr == '\n' && *(itr + 1) == '@') + return false; + return true; +} + auto sam_stdin::make_tasks(const std::int64_t n_chunks, // const std::int32_t file_id, // @@ -127,5 +142,7 @@ sam_stdin::make_tasks(const std::int64_t n_chunks, // n_tasks = 1; // for current task, which makes more tasks shift_output_buffer(); load_next(); + if (!validate_sam(buffer)) + throw std::runtime_error("input appears not to be SAM"); get_chunks(n_chunks, file_id, tq, n_tasks); }