Adding option to take input from stdin

CMakeLists.txt

@@ -83,13 +83,15 @@ 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)
 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)
@@ -131,11 +133,13 @@ target_link_libraries(falco PRIVATE
   original_duplicates
   kmer_counter
   fastq_file
+  fastq_stdin
   fastq_gz_file
   fastq_bgzf_file
   bam_file
   bamrec
   sam_file
+  sam_stdin
   samrec
   bam_header
   falco_file_format

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)

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
+a1aa02211d63f507646a9e12122f0226  sam_stdin_out/sam_stdin/fastqc_data.txt
+c61e37d2c8d64b6950a5dad300e5fd60  fastq_stdin_out/fastq_stdin/fastqc_data.txt

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

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

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<overrep_t> {
-  const auto cutoff = static_cast<double>(n_reads) * overrep_cutoff;
+  const auto cutoff =
+    std::max(static_cast<double>(n_reads) * overrep_frac_cutoff,
+             static_cast<double>(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) |

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};

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,13 +56,15 @@ 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"
 #include "reads_file.hpp"  // IWYU pragma: keep
 #include "results_summary.hpp"
 #include "run_mode.hpp"
 #include "sam_file.hpp"
+#include "sam_stdin.hpp"
 #include "tile_processor.hpp"
 
 #include "CLI11/CLI11.hpp"
@@ -74,7 +82,6 @@ Use these as templates. Copy and modify them to customize your analysis.
 #include <fstream>
 #include <functional>
 #include <iterator>
-#include <limits>
 #include <map>
 #include <memory>
 #include <print>
@@ -142,6 +149,44 @@ make_reads_files(const std::vector<file_info> &infos,
   return reads_files;
 }
 
+[[nodiscard]] static auto
+make_reads_file_stdin(const std::vector<file_info> &infos,
+                      const std::int64_t buf_size)
+  -> std::vector<reads_file_t> {
+  std::vector<reads_file_t> 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<std::string> &names,
+                    const std::pair<falco::file_format, std::uint32_t> &ft_tile)
+  -> std::vector<file_info> {
+  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;
+  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;
+  info.n_reads_est = n_reads_est;
+  info.read_len_est = read_len_est;
+  return std::vector<file_info>(1, info);
+}
+
 [[nodiscard]] static auto
 get_file_info(const auto &infiles) {
   std::vector<file_info> infos;
@@ -150,17 +195,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(),
@@ -227,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<std::string> infiles;
     std::string contam_file;
     std::string config_file;
@@ -256,17 +302,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<FormatWithoutFlagDefaults>();
     app.formatter(fmt);
@@ -277,18 +328,21 @@ 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");
@@ -298,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<std::int64_t>::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);
@@ -321,6 +375,23 @@ main(int argc, char *argv[]) {
       ->option_text(" ")
       ->capture_default_str()
       ->transform(size_from_units);
+    app.add_option_function<std::pair<falco::file_format, std::uint32_t>>(
+      "--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))
+      // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers)
+      ->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 +429,16 @@ 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())
       load_config_and_set_graders(config_file, mode);
@@ -401,12 +482,13 @@ 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; };
     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) {
@@ -440,9 +522,11 @@ main(int argc, char *argv[]) {
     auto dups = do_original_dups
                   ? initialize_original_duplicates(infiles, infos, n_threads)
                   : std::vector<dups_init_t>{};
+    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)