IUPAC symbols for orig-dups mode: non-ACGT (case insensitive) are converted to N

CMakeLists.txt

@@ -130,6 +130,8 @@ target_link_libraries(falco_file_format PRIVATE HTSLIB::HTSLIB)
 target_link_libraries(tile_processor PRIVATE HTSLIB::HTSLIB)
 target_link_libraries(quality_score PRIVATE HTSLIB::HTSLIB)
 target_link_libraries(fastq_gz_file PRIVATE HTSLIB::HTSLIB)
+target_link_libraries(original_duplicates PRIVATE HTSLIB::HTSLIB)
+
 if(USE_ISAL)
   target_link_libraries(fastq_gz_file PUBLIC ISAL::ISAL)
 endif()

cmake/tests.cmake

@@ -38,3 +38,4 @@ 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)
+add_test(NAME "IUPAC and orig dups" COMMAND bash test_scripts/iupac_dups.sh)

data/test_data/md5sum.txt

@@ -1,6 +1,6 @@
 e4f4641559aadc9ec7b5c50cc4aa82e3  bam_kmers_out/bam_1/fastqc_data.txt
 090a74739499d287a7aeea689a78b4e9  bam_mult_out/bam_1/fastqc_data.txt
-7cc30a3a14796140f05940485ab1bed5  bam_mult_out/bam_2/fastqc_data.txt
+e3568c45f8cc358c263e51c4ece215cd  bam_mult_out/bam_2/fastqc_data.txt
 090a74739499d287a7aeea689a78b4e9  bam_orig_dups_out/bam_1/fastqc_data.txt
 090a74739499d287a7aeea689a78b4e9  bam_out/bam_1/fastqc_data.txt
 090a74739499d287a7aeea689a78b4e9  bam_tiles_out/bam_1/fastqc_data.txt
@@ -17,3 +17,5 @@ e786bd4d9e6065904ab7b45021e7190b  groups_out/bam_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
+c827b3a42f08404ac87d3cfb1dc51938  iupac_dups_out/iupac_dups/fastqc_data.txt
+7f06e30c51af849eb74c9066d087173b  iupac_dups_out/iupac_dups_bam/fastqc_data.txt

data/test_scripts/iupac_dups.sh

@@ -0,0 +1,30 @@
+#!/usr/bin/env bash
+# SPDX-License-Identifier: MIT
+
+# Test that IUPAC symbols are counted as N
+
+# Need samtools to check that this works for sam/bam formats
+if [[ ! -e $(type -P samtools) ]]; then
+    echo "samtools not found; skipping remaining tests";
+    exit 77;
+fi
+
+prog=./falco
+infile1=test_data/iupac_dups.fq.gz
+infile2=iupac_dups_bam.bam
+outdir=iupac_dups_out
+if [[ -e "${infile1}" ]]; then
+    ${prog} --orig-dups -o ${outdir} ${infile1}
+    samtools view -h -o ${outdir}/${infile2} ${infile1}
+    ${prog} --orig-dups -o ${outdir} ${outdir}/${infile2}
+    x=$(md5sum --ignore-missing -c test_data/md5sum.txt | \
+            grep "${outdir}" | \
+            grep -c "OK$")
+    if [[ "${x}" != "2" ]]; then
+        exit 1;
+    fi
+    rm -r ${outdir}
+else
+    echo "${infile1} not found; skipping remaining tests";
+    exit 77;
+fi

src/adapter_matcher.cpp

@@ -55,13 +55,17 @@ adapter_matcher::add_and_consume(adapter_matcher &&rhs) -> void {
 [[nodiscard]] auto
 adapter_matcher::get_grade(const std::uint64_t n_reads) const -> std::string {
   static constexpr auto label = "adapter";
-  const auto max_count =
-    adap_counts.empty()
-      ? 0
-      : std::ranges::max(adap_counts | std::views::transform(std::ranges::max));
+  // MRS (maximilian.stammnitz@crg.eu):
+  // The report shows, per adapter, the cumulative percentage of reads with the
+  // adapter starting at or before each position; as in FastQC, the grade is
+  // the maximum of that, i.e. the total count of each adapter over positions.
+  std::vector<std::uint64_t> totals(n_adapters, 0);
+  for (const auto &counts : adap_counts)
+    std::ranges::transform(totals, counts, std::begin(totals), std::plus{});
+  const auto max_count = totals.empty() ? 0 : std::ranges::max(totals);
   return grader_set::get_grade(
     label,
-    as_frac(max_count, std::max(static_cast<std::uint64_t>(1), n_reads)));
+    pct(as_frac(max_count, std::max(static_cast<std::uint64_t>(1), n_reads))));
 }
 
 [[nodiscard]] auto

src/bam_file.cpp

@@ -2,9 +2,7 @@
 
 #include "bam_file.hpp"
 #include "bamrec.hpp"
-#include "duplication_results.hpp"
 #include "falco_utils.hpp"
-#include "falco_word.hpp"
 #include "task_queue.hpp"
 
 #include <htslib/bgzf.h>  // IWYU pragma: keep
@@ -65,52 +63,6 @@ estimate_n_reads_bam(const std::string &filename)
   return {n_reads_estimate, read_len_estimate, filesize};
 }
 
-[[nodiscard]] auto
-init_dups(const std::string &filename, const std::uint64_t n_unique)
-  -> dups_map_t {
-  static constexpr auto complem = [](const auto x) {
-    return "TNGNNNCNNNNNNNNNNNNA"[x - 'A'];
-  };
-  static constexpr auto n_unique_multiplier = 10;
-  const std::unique_ptr<htsFile, int (*)(htsFile *)> in(
-    hts_open(std::data(filename), "r"), &hts_close);
-  if (!in)
-    throw std::system_error(std::make_error_code(std::errc(errno)),
-                            "failed to open file: " + filename);
-  const std::unique_ptr<sam_hdr_t, void (*)(sam_hdr_t *)> h(
-    sam_hdr_read(in.get()), &sam_hdr_destroy);
-  if (!h)
-    throw std::system_error(std::make_error_code(std::errc(errno)),
-                            "failed to read header: " + filename);
-  // NOLINTNEXTLINE(cppcoreguidelines-pro-type-union-access)
-  std::unique_ptr<bam1_t, void (*)(bam1_t *)> rec(bam_init1(), &bam_destroy1);
-
-  dups_map_t dups;
-  std::vector<char> buffer;
-  int r{};
-  std::uint64_t n_reads{};
-  const auto max_n_reads = n_unique_multiplier * n_unique;
-  while (n_reads++ < max_n_reads && std::size(dups) < n_unique &&
-         (r = sam_read1(in.get(), h.get(), rec.get())) >= 0) {
-    const auto l_qseq = rec->core.l_qseq;
-    const auto seq = bam_get_seq(rec.get());
-    if (std::ssize(buffer) < l_qseq)
-      buffer.resize(l_qseq);
-    for (auto i = 0; i < l_qseq; ++i)
-      buffer[i] = seq_nt16_str[bam_seqi(seq, i)];
-    if (bam_is_rev(rec.get())) {
-      std::reverse(std::begin(buffer), std::begin(buffer) + l_qseq);
-      std::transform(std::cbegin(buffer), std::cbegin(buffer) + l_qseq,
-                     std::begin(buffer), complem);
-    }
-    ++dups[falco_word(std::data(buffer), l_qseq)];
-  }
-  if (r < -1)  // error
-    throw std::system_error(std::make_error_code(std::errc(errno)),
-                            "error reading bam record from: " + filename);
-  return dups;
-}
-
 auto
 bam_file::load_next(const std::int32_t file_id,
                     task_queue &tq,

src/bam_file.hpp

@@ -6,7 +6,6 @@
 #include "bam_header.hpp"
 #include "bgzf_block.hpp"
 #include "bgzf_reader.hpp"
-#include "duplication_results.hpp"
 
 #include <atomic>
 #include <cstdint>
@@ -129,10 +128,6 @@ class bam_file {
 estimate_n_reads_bam(const std::string &filename)
   -> std::tuple<std::uint64_t, std::uint64_t, std::int64_t>;
 
-[[nodiscard]] auto
-init_dups(const std::string &filename, const std::uint64_t n_unique)
-  -> dups_map_t;
-
 inline auto
 make_tasks(bam_file &reads_file,
            const std::int64_t n_threads,

src/duplication_results.cpp

@@ -27,6 +27,7 @@
 #include <numeric>
 #include <ranges>
 #include <span>
+#include <stdexcept>
 #include <string>
 #include <tuple>
 #include <unordered_map>
@@ -157,11 +158,13 @@ duplication_results::add_and_consume(duplication_results &rhs) -> void {
 get_grade_overrepresented(const std::uint64_t n_reads,
                           const duplication_results &dr) -> std::string {
   static constexpr auto label = "overrepresented";
-  const auto max_n_obs =
-    dr.dups.empty() ? 0LU : std::ranges::max(std::views::values(dr.dups));
-  return grader_set::get_grade(
-    label,
-    as_frac(max_n_obs, std::max(static_cast<std::uint64_t>(1), n_reads)));
+  // MRS (maximilian.stammnitz@crg.eu):
+  // As in FastQC, grade on the list the report shows: the percentage of the
+  // most frequent listed sequence (in percent, like the cutoffs), or zero if
+  // no sequence passes the list's cutoffs.
+  const auto overrep = dr.get_overrepresented(n_reads);
+  return grader_set::get_grade(label,
+                               overrep.empty() ? 0.0 : overrep.front().pct_val);
 }
 
 [[nodiscard]] auto
@@ -236,6 +239,8 @@ duplication_results::get_dups_summary(const std::uint64_t n_reads) const
   std::vector<std::uint64_t> hist_dedup(max_dup + 1);
   for (const auto n_copies : std::views::values(dups))
     ++hist_dedup[n_copies];
+  if (hist_dedup[0] > 0)
+    throw std::runtime_error("Non-zero count of indexed reads seen 0 times");
   for (auto [idx, val] : falco::views::enumerate(hist_dedup))
     val = static_cast<std::uint64_t>(
       get_corrected_count(count_at_limit, n_reads, idx, val));
@@ -283,7 +288,7 @@ get_grade_duplication(const dup_summary_t &summary) -> std::string {
   };
   const auto frac_dedup =
     as_frac(reduce(summary.hist_dedup), reduce(summary.hist_mass));
-  return grader_set::get_grade(label, frac_dedup);
+  return grader_set::get_grade(label, pct(frac_dedup));
 }
 
 [[nodiscard]] auto

src/falco_grade.cpp

@@ -24,15 +24,15 @@ static const auto default_graders = std::array{  // NOLINT(cert-err58-cpp)
   grader{.name = "quality_base_median", .warn = 25.0, .fail = 20.0},
   grader{.name = "quality_base_lower",  .warn = 10.0, .fail = 5.00},
   grader{.name = "quality_sequence",    .warn = 27.0, .fail = 20.0},
-  grader{.name = "n_content",           .warn = 0.05, .fail = 0.20},
+  grader{.name = "n_content",           .warn = 5.00, .fail = 20.0},
   grader{.name = "sequence",            .warn = 10.0, .fail = 20.0},
   grader{.name = "sequence_length",     .warn = 1.00, .fail = 1.00},
   grader{.name = "gc_sequence",         .warn = 15.0, .fail = 30.0},
-  grader{.name = "duplication",         .warn = 0.70, .fail = 0.50},
+  grader{.name = "duplication",         .warn = 70.0, .fail = 50.0},
   grader{.name = "overrepresented",     .warn = 0.10, .fail = 1.00},
   grader{.name = "kmer",                .warn = 2.00, .fail = 5.00},
   grader{.name = "tile",                .warn = 5.00, .fail = 10.0},
-  grader{.name = "adapter",             .warn = 0.05, .fail = 0.10},
+  grader{.name = "adapter",             .warn = 5.00, .fail = 10.0},
 };
 // clang-format on
 
@@ -74,9 +74,12 @@ file_grades::get_title(const std::string &name) -> std::string {
 grader::identify_grade(const double value) const -> std::string {
   using std::string_literals::operator""s;
   if (warn < fail) {
-    if (value < warn)
+    // MRS (maximilian.stammnitz@crg.eu):
+    // higher is worse: as in FastQC, a value equal to a cutoff does not
+    // reach it (e.g., 0.1% overrepresented passes, more than 0.1% warns)
+    if (value <= warn)
       return "pass"s;
-    if (value < fail)
+    if (value <= fail)
       return "warn"s;
     return "fail"s;
   }
@@ -198,7 +201,7 @@ get_grade_n_content(const std::vector<std::uint64_t> &n_counts,
     std::reduce(std::cbegin(nucs[max_idx]), std::cend(nucs[max_idx]));
   const auto max_n_frac =
     as_frac(n_counts[max_idx], n_counts[max_idx] + total_non_n);
-  return grader_set::get_grade(label, max_n_frac);
+  return grader_set::get_grade(label, pct(max_n_frac));
 }
 
 [[nodiscard]] auto

src/fastq_bgzf_file.cpp

@@ -1,15 +1,13 @@
 // SPDX-License-Identifier: MIT; Copyright 2026 Andrew D Smith
 
 #include "fastq_bgzf_file.hpp"
-#include "duplication_results.hpp"
 #include "falco_utils.hpp"
-#include "falco_word.hpp"
 #include "fqrec.hpp"
 #include "task_queue.hpp"
 
 #include <htslib/bgzf.h>
 #include <htslib/hfile.h>
-#include <htslib/sam.h>
+#include <htslib/sam.h>  // IWYU pragma: keep
 
 #include <algorithm>
 #include <cassert>
@@ -55,38 +53,6 @@ estimate_n_reads_fastq_bgzf(const std::string &filename)
   return {static_cast<std::uint64_t>(n_reads_est), read_len_est, filesize};
 }
 
-[[nodiscard]] auto
-init_dups_fq(const std::string &filename, const std::uint64_t n_unique)
-  -> dups_map_t {
-  static constexpr auto n_unique_multiplier = 10;
-  const std::unique_ptr<BGZF, int (*)(BGZF *)> in(
-    bgzf_open(std::data(filename), "r"), &bgzf_close);
-  if (!in)
-    throw std::system_error(std::make_error_code(std::errc(errno)),
-                            "failed to open file: " + filename);
-  dups_map_t dups;
-  kstring_t s = KS_INITIALIZE;
-  int r{};
-  std::uint64_t n_lines{};
-  std::uint64_t n_reads{};
-  const auto max_n_reads = n_unique_multiplier * n_unique;
-  while (n_reads < max_n_reads && std::size(dups) < n_unique &&
-         (r = bgzf_getline(in.get(), '\n', &s)) >= 0) {
-    if (n_lines % 4 == 1) {
-      const auto l_qseq = ks_len(&s);
-      const auto seq = ks_c_str(&s);
-      ++dups[falco_word(seq, l_qseq)];
-      ++n_reads;
-    }
-    ++n_lines;
-  }
-  if (r < -1)  // error
-    throw std::system_error(std::make_error_code(std::errc(errno)),
-                            "error reading bam record from: " + filename);
-  ks_free(&s);
-  return dups;
-}
-
 auto
 fastq_bgzf_file::load_next(const std::int32_t file_id,
                            task_queue &tq,

src/fastq_bgzf_file.hpp

@@ -5,7 +5,6 @@
 
 #include "bgzf_block.hpp"
 #include "bgzf_reader.hpp"
-#include "duplication_results.hpp"
 
 #include <atomic>
 #include <cstdint>
@@ -128,10 +127,6 @@ class fastq_bgzf_file {
 estimate_n_reads_fastq_bgzf(const std::string &filename)
   -> std::tuple<std::uint64_t, std::uint64_t, std::int64_t>;
 
-[[nodiscard]] auto
-init_dups_fq(const std::string &filename, const std::uint64_t n_unique)
-  -> dups_map_t;
-
 inline auto
 make_tasks(fastq_bgzf_file &reads_file,
            const std::int64_t n_threads,