Adding robustness to IUPAC degenerate nucleotides in reads

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)

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

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)

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

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

src/bamrec.cpp

@@ -37,23 +37,20 @@ bamrec::to_string() const -> std::string {
 template <class bidir_itr_t, class output_itr_t>
 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 <class bidir_itr_t, class output_itr_t>
 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;

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<FormatWithoutFlagDefaults>();
-    app.formatter(fmt);
+    const auto app_fmt = std::make_shared<FormatWithoutFlagDefaults>();
+    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()))

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<char> 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<fq_task_t>(task))
-              process_reads(res[file_id], std::get<fq_task_t>(task));
+              process_reads(res[file_id], std::get<fq_task_t>(task)
+#ifdef USE_MMAP_FOR_PLAIN_FASTQ
+                            , fq_buffer
+#endif  // USE_MMAP_FOR_PLAIN_FASTQ
+                            );
+            // clang-format on
+            else if (std::holds_alternative<fq_gz_task_t>(task))
+              process_reads(res[file_id], std::get<fq_gz_task_t>(task));
             else if (std::holds_alternative<bam_task_t>(task))
               process_reads(res[file_id], std::get<bam_task_t>(task));
             else if (std::holds_alternative<sam_task_t>(task))

src/falco_task.hpp

@@ -10,7 +10,11 @@
 
 #include <variant>
 
-using task_t =
-  std::variant<std::monostate, fq_task_t, bam_task_t, sam_task_t, bgzf_block_t>;
+using task_t = std::variant<std::monostate,
+                            fq_task_t,
+                            fq_gz_task_t,
+                            bam_task_t,
+                            sam_task_t,
+                            bgzf_block_t>;
 
 #endif  // SRC_FALCO_TASK_HPP_

src/falco_utils.hpp

@@ -20,6 +20,10 @@
 #include <type_traits>
 #include <vector>
 
+/// 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

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;