diff --git a/CMakeLists.txt b/CMakeLists.txt index b4c0948..c7a7d41 100644 --- a/CMakeLists.txt +++ b/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() diff --git a/cmake/tests.cmake b/cmake/tests.cmake index aac54a0..3097ff0 100644 --- a/cmake/tests.cmake +++ b/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) diff --git a/data/test_data/iupac_dups.fq.gz b/data/test_data/iupac_dups.fq.gz new file mode 100644 index 0000000..4c04c35 Binary files /dev/null and b/data/test_data/iupac_dups.fq.gz differ diff --git a/data/test_data/md5sum.txt b/data/test_data/md5sum.txt index d227693..7fcce8d 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 @@ -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 diff --git a/data/test_scripts/iupac_dups.sh b/data/test_scripts/iupac_dups.sh new file mode 100644 index 0000000..1b34708 --- /dev/null +++ b/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 diff --git a/src/adapter_matcher.cpp b/src/adapter_matcher.cpp index bc9552c..bc0ea1e 100644 --- a/src/adapter_matcher.cpp +++ b/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 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/bam_file.cpp b/src/bam_file.cpp index 90533a8..b8a28f7 100644 --- a/src/bam_file.cpp +++ b/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 // 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 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 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 rec(bam_init1(), &bam_destroy1); - - dups_map_t dups; - std::vector 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, diff --git a/src/bam_file.hpp b/src/bam_file.hpp index 97a2473..f781c89 100644 --- a/src/bam_file.hpp +++ b/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 #include @@ -129,10 +128,6 @@ class bam_file { estimate_n_reads_bam(const std::string &filename) -> std::tuple; -[[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, diff --git a/src/duplication_results.cpp b/src/duplication_results.cpp index 0a33242..ae1f616 100644 --- a/src/duplication_results.cpp +++ b/src/duplication_results.cpp @@ -27,6 +27,7 @@ #include #include #include +#include #include #include #include @@ -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(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 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( 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 diff --git a/src/falco_grade.cpp b/src/falco_grade.cpp index 28896e0..d2133a9 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,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 &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 diff --git a/src/fastq_bgzf_file.cpp b/src/fastq_bgzf_file.cpp index 719e928..68536e5 100644 --- a/src/fastq_bgzf_file.cpp +++ b/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 #include -#include +#include // IWYU pragma: keep #include #include @@ -55,38 +53,6 @@ estimate_n_reads_fastq_bgzf(const std::string &filename) return {static_cast(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 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, diff --git a/src/fastq_bgzf_file.hpp b/src/fastq_bgzf_file.hpp index 0023291..8656e2d 100644 --- a/src/fastq_bgzf_file.hpp +++ b/src/fastq_bgzf_file.hpp @@ -5,7 +5,6 @@ #include "bgzf_block.hpp" #include "bgzf_reader.hpp" -#include "duplication_results.hpp" #include #include @@ -128,10 +127,6 @@ class fastq_bgzf_file { estimate_n_reads_fastq_bgzf(const std::string &filename) -> std::tuple; -[[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, diff --git a/src/original_duplicates.cpp b/src/original_duplicates.cpp index d0b960b..44629d5 100644 --- a/src/original_duplicates.cpp +++ b/src/original_duplicates.cpp @@ -2,29 +2,132 @@ #include "original_duplicates.hpp" -#include "bam_file.hpp" #include "duplication_results.hpp" #include "falco_file_format.hpp" -#include "fastq_bgzf_file.hpp" +#include "falco_word.hpp" #include "file_info.hpp" #include "boost/boost_unordered.hpp" +#include // IWYU pragma: keep +#include + #include #include +#include #include +#include #include #include +#include #include +#include #include #include // for std::move #include +namespace { +[[nodiscard]] auto +init_orig_dups_bam(const std::string &filename, const std::uint64_t n_unique) + -> dups_map_t { + // ADS: n_unique_multiplier: contributes to max reads to examine, multiplying + // the n_unique that is the target; impossible to know how many needed to hit + // the target, so try to overshoot. + static constexpr auto n_unique_multiplier = 10; + static constexpr auto complem = [](char &x) { + assert(x >= 'A' && x <= 'T'); + x = "TNGNNNCNNNNNNNNNNNNA"[x - 'A']; + }; + static constexpr auto to_acgtn = [](const char x) { + static constexpr char mask = 15; + static constexpr std::span seq_nt16_str_acgtn = "NACNGNNNTNNNNNNN"; + return seq_nt16_str_acgtn[x & mask]; + }; + const std::unique_ptr 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 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 rec(bam_init1(), &bam_destroy1); + + dups_map_t dups; + std::vector 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] = to_acgtn(bam_seqi(seq, i)); + if (bam_is_rev(rec.get())) { + std::reverse(std::begin(buffer), std::begin(buffer) + l_qseq); + std::for_each_n(std::begin(buffer), l_qseq, 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; +} + +[[nodiscard]] auto +init_orig_dups_fq(const std::string &filename, const std::uint64_t n_unique) + -> dups_map_t { + // ADS: n_unique_multiplier: contributes to max reads to examine, multiplying + // the n_unique that is the target; impossible to know how many needed to hit + // the target, so try to overshoot. + static constexpr auto n_unique_multiplier = 10; + static constexpr auto to_acgtn = [](const auto c) { + static constexpr auto mask = 31; + static constexpr std::span to_acgtn_sp = "NANCNNNGNNNNNNNNNNNNTNNNNNNNNNNN"; + return to_acgtn_sp[c & mask]; + }; + const std::unique_ptr 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); + std::ranges::transform(std::span(s.s, l_qseq), s.s, to_acgtn); + 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; +} +} // namespace + [[nodiscard]] auto -initialize_original_duplicates( - const std::vector &infiles, - const std::vector &infos, - const std::uint32_t n_threads) -> std::vector { +initialize_original_duplicates(const std::vector &infiles, + 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); @@ -47,10 +150,10 @@ initialize_original_duplicates( } dups[file_id] = falco::is_mapped_reads(infos[file_id].format) - ? init_dups(infiles[file_id], - duplication_results::max_reads_to_hash) - : init_dups_fq(infiles[file_id], - duplication_results::max_reads_to_hash); + ? init_orig_dups_bam(infiles[file_id], + duplication_results::max_reads_to_hash) + : init_orig_dups_fq(infiles[file_id], + duplication_results::max_reads_to_hash); } }); std::ranges::for_each(workers, [](auto &w) { w.join(); });