From 09ef52395f683eff003386a34026b9c53314fb9e Mon Sep 17 00:00:00 2001 From: Maximilian Stammnitz Date: Mon, 2 Nov 2026 14:49:55 +0100 Subject: [PATCH] src: grade adapter, overrepresented, N content and duplication as FastQC Four modules were graded against cutoffs in different units from the values compared with them, or on a different quantity from the one the report shows: - overrepresented: the value was a fraction of reads, the default cutoffs (0.1, 1) are percent, so the effective cutoffs were 10% and 100%. The grade is now the percentage of the top sequence in the reported list (zero if the list is empty), as in FastQC, so a sequence below the list's own cutoffs can no longer raise it. - adapter: the value was the largest per-position count of an adapter start, while the report (and FastQC) shows the cumulative percentage of reads with the adapter at or before each position. The grade is now the maximum of that cumulative percentage. - n_content, duplication: fractions compared with fractional default cutoffs, which worked, but a config file in the FastQC limits.txt format (percent: n_content 5/20, duplication 70/50) made N content unable to warn and duplication always fail. Both are now graded in percent and the defaults are given in percent, like every other module and like Configuration/limits.txt. Also, for modules where higher is worse, a value equal to a cutoff no longer reaches it (FastQC tests value > cutoff), e.g. a sequence at exactly 0.1% passes. data/test_data/md5sum.txt: bam_2 now warns for overrepresented sequences (top sequence 0.12% of reads), as FastQC 0.13.0 does. --- data/test_data/md5sum.txt | 2 +- src/adapter_matcher.cpp | 13 ++++++++----- src/duplication_results.cpp | 13 +++++++------ src/falco_grade.cpp | 14 ++++++++------ 4 files changed, 24 insertions(+), 18 deletions(-) diff --git a/data/test_data/md5sum.txt b/data/test_data/md5sum.txt index d227693..0987ae9 100644 --- a/data/test_data/md5sum.txt +++ b/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 diff --git a/src/adapter_matcher.cpp b/src/adapter_matcher.cpp index bc9552c..e088e15 100644 --- a/src/adapter_matcher.cpp +++ b/src/adapter_matcher.cpp @@ -55,13 +55,16 @@ 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)); + // 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 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(1), n_reads))); + pct(as_frac(max_count, std::max(static_cast(1), n_reads)))); } [[nodiscard]] auto diff --git a/src/duplication_results.cpp b/src/duplication_results.cpp index 0a33242..bfdb405 100644 --- a/src/duplication_results.cpp +++ b/src/duplication_results.cpp @@ -157,11 +157,12 @@ 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(1), n_reads))); + // 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 @@ -283,7 +284,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 diff --git a/src/falco_grade.cpp b/src/falco_grade.cpp index 28896e0..7c558e7 100644 --- a/src/falco_grade.cpp +++ b/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,11 @@ 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) + // 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 +200,7 @@ get_grade_n_content(const std::vector &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