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) 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 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) diff --git a/data/test_data/md5sum.txt b/data/test_data/md5sum.txt index ff3aa1f..d227693 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 +e71cef09fd4b6f4019f8fdb82ee2eaf3 fastq_iupac_out/fastq_1_iupac/fastqc_data.txt diff --git a/data/test_scripts/fastq_iupac.sh b/data/test_scripts/fastq_iupac.sh new file mode 100644 index 0000000..25e5f3f --- /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} -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 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; diff --git a/src/falco.cpp b/src/falco.cpp index 41627b5..a94d478 100644 --- a/src/falco.cpp +++ b/src/falco.cpp @@ -3,38 +3,39 @@ // 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. Get it with + these packages: tabix (apt), bgzip (pip and dnf) or htslib (brew and conda). +- Falco assumes nucleotides symbols ACGT and N (upper or lowercase). Any other + symbols are counted as N. 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 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 -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 +265,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 +338,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 +353,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( @@ -379,7 +382,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())) diff --git a/src/falco_analyzer.cpp b/src/falco_analyzer.cpp index 1b9bd88..bfda9e7 100644 --- a/src/falco_analyzer.cpp +++ b/src/falco_analyzer.cpp @@ -55,6 +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]); @@ -63,8 +69,16 @@ 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)); + 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)) process_reads(res[file_id], std::get(task)); else if (std::holds_alternative(task)) 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_ diff --git a/src/falco_utils.hpp b/src/falco_utils.hpp index aaf934d..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 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_file.cpp b/src/fastq_file.cpp index 2abd3e7..c634a1e 100644 --- a/src/fastq_file.cpp +++ b/src/fastq_file.cpp @@ -5,27 +5,33 @@ #include "fqrec.hpp" #include "task_queue.hpp" -#include -#include +#include // IWYU pragma: keep #include +#ifdef USE_MMAP_FOR_PLAIN_FASTQ +#include #include +#endif // USE_MMAP_FOR_PLAIN_FASTQ #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,15 +80,106 @@ 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"); +#ifndef __APPLE__ + posix_fadvise(fileno(in.get()), 0, 0, POSIX_FADV_SEQUENTIAL); +#endif + 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{ 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 @@ -166,10 +263,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 +276,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 +285,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..5aeb267 100644 --- a/src/fastq_file.hpp +++ b/src/fastq_file.hpp @@ -3,17 +3,91 @@ #ifndef SRC_FASTQ_FILE_HPP_ #define SRC_FASTQ_FILE_HPP_ +#ifdef USE_MMAP_FOR_PLAIN_FASTQ #include +#endif #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 +157,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; 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; diff --git a/src/fqrec.hpp b/src/fqrec.hpp index 6737415..7018539 100644 --- a/src/fqrec.hpp +++ b/src/fqrec.hpp @@ -57,10 +57,26 @@ 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 { + 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 +91,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); + while (itr != end_itr && *itr != '\n') + iupac(*itr++); if (itr++ == end_itr) return {}; const auto o = itr; 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 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); diff --git a/src/results_collector.hpp b/src/results_collector.hpp index 5e521fc..431bee8 100644 --- a/src/results_collector.hpp +++ b/src/results_collector.hpp @@ -196,8 +196,29 @@ accumulate_results(std::vector> &r, } } +// clang-format off inline auto -process_reads(results_collector &results, const fq_task_t &task) { +process_reads(results_collector &results, + 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) { results.process_reads_fq(task.beg, task.end); } 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; 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