From 2708965917d3f359e67fad3e5990707669b2c1ed Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Fri, 26 Sep 2025 17:20:55 -0700 Subject: [PATCH 01/11] src/utils/uniq.cpp: computing the derived values when they are printed and only keeping track of the fundamental counts; this is so the db will be in normal form --- src/utils/uniq.cpp | 46 ++++++++++++++++++++++++++-------------------- 1 file changed, 26 insertions(+), 20 deletions(-) diff --git a/src/utils/uniq.cpp b/src/utils/uniq.cpp index 44f11ffb..9339ba75 100644 --- a/src/utils/uniq.cpp +++ b/src/utils/uniq.cpp @@ -73,18 +73,10 @@ struct rd_stats { // keep track of good bases/reads in and out struct uniq_summary { uniq_summary(const rd_stats &rs_in, const rd_stats &rs_out, - const std::size_t reads_duped) { - total_reads = rs_in.reads; - total_bases = rs_in.bases; - unique_reads = rs_out.reads; - unique_read_bases = rs_out.bases; - reads_removed = rs_in.reads - rs_out.reads; - non_duplicate_fraction = static_cast(rs_out.reads - reads_duped) / - std::max(1ul, rs_in.reads); - duplication_rate = static_cast(reads_removed + reads_duped) / - std::max(1ul, reads_duped); - duplicate_reads = reads_duped; - } + const std::size_t reads_duped) : + total_reads(rs_in.reads), total_bases(rs_in.bases), + unique_reads(rs_out.reads), unique_read_bases(rs_out.bases), + duplicate_reads(reads_duped) {} // total_reads is the number of input reads std::size_t total_reads{}; @@ -94,16 +86,30 @@ struct uniq_summary { std::size_t unique_reads{}; // unique_read_bases is the total number of bases for the unique reads std::size_t unique_read_bases{}; - // non_duplicate_fraction is the ratio of the number of unique reads with - // no duplicates to that of the input reads - double non_duplicate_fraction{}; // duplicate_reads is the number of unique reads with at least one duplicate std::size_t duplicate_reads{}; + // reads_removed is the number of duplicate reads that have been removed - std::size_t reads_removed{}; + [[nodiscard]] std::size_t + reads_removed() const { + return total_reads - unique_reads; + } + + // non_duplicate_fraction is the ratio of the number of unique reads with + // no duplicates to that of the input reads + [[nodiscard]] double + non_duplicate_fraction() const { + return static_cast(unique_reads - duplicate_reads) / + std::max(1ul, total_reads); + } + // duplication_rate is the average number of duplicates for the reads with // at least one duplicate (>1 by definition) - double duplication_rate{}; + [[nodiscard]] double + duplication_rate() const { + return static_cast(reads_removed() + duplicate_reads) / + std::max(1ul, duplicate_reads); + } std::string to_string() { @@ -112,10 +118,10 @@ struct uniq_summary { << "total_bases: " << total_bases << "\n" << "unique_reads: " << unique_reads << "\n" << "unique_read_bases: " << unique_read_bases << "\n" - << "non_duplicate_fraction: " << non_duplicate_fraction << "\n" + << "non_duplicate_fraction: " << non_duplicate_fraction() << "\n" << "duplicate_reads: " << duplicate_reads << "\n" - << "reads_removed: " << reads_removed << "\n" - << "duplication_rate: " << duplication_rate; + << "reads_removed: " << reads_removed() << "\n" + << "duplication_rate: " << duplication_rate(); return oss.str(); } }; From 7f30aad43314006f37313188f6a7c5fd96749c12 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Fri, 26 Sep 2025 17:21:17 -0700 Subject: [PATCH 02/11] src/analysis/bsrate.cpp: computing the derived values when they are printed and only keeping track of the fundamental counts; this is so the db will be in normal form --- src/analysis/bsrate.cpp | 451 +++++++++++++++++++++------------------- 1 file changed, 235 insertions(+), 216 deletions(-) diff --git a/src/analysis/bsrate.cpp b/src/analysis/bsrate.cpp index e0324bff..e5bbc6e6 100644 --- a/src/analysis/bsrate.cpp +++ b/src/analysis/bsrate.cpp @@ -31,78 +31,109 @@ #include #include #include +#include #include #include -using std::accumulate; -using std::cerr; -using std::cout; -using std::endl; -using std::for_each; -using std::max; -using std::numeric_limits; -using std::pair; -using std::runtime_error; -using std::string; -using std::unordered_map; -using std::unordered_set; -using std::vector; - -using bamxx::bam_rec; - struct bsrate_summary { // converted_count_positive is the number of nucleotides covering a // cytosine in the reference that show a thymine in the read, and // for reads mapping to the positive strand. - uint64_t converted_count_positive{}; + std::uint64_t converted_count_positive{}; + // total_count_positive is the number of nucleotides covering a // cytosine in the reference that show either a cytosine or a // thymine in the read, and for reads mapping to the positive // strand. - uint64_t total_count_positive{}; - // bisulfite_conversion_rate_positive is equal to - // converted_count_positive divided by total_count_positive, a value - // that is always between 0 and 1. When total_count_positive is 0, - // then bisulfite_conversion_rate_positive is given a value of 0. - double bisulfite_conversion_rate_positive{}; + std::uint64_t total_count_positive{}; // converted_count_negative is the number of nucleotides covering a // cytosine in the reference that show a thymine in the read, and // for reads mapping to the negative strand. - uint64_t converted_count_negative{}; + std::uint64_t converted_count_negative{}; + // total_count_negative is the number of nucleotides covering a // cytosine in the reference that show either a cytosine or a // thymine in the read, and for reads mapping to the negative // strand. - uint64_t total_count_negative{}; - // bisulfite_conversion_rate_negative is equal to - // converted_count_negative divided by total_count_negative, a value - // that is always between 0 and 1. When total_count_negative is 0, - // then bisulfite_conversion_rate_negative is given a value of 0. - double bisulfite_conversion_rate_negative{}; + std::uint64_t total_count_negative{}; + + // error_count_positive is the number of nucleotides covering a cytosine in + // the reference shows either an A or a G in a read mapping on the positive + // strand. + std::uint64_t error_count_positive{}; + + // error_count_negative is the number of nucleotides covering a cytosine in + // the reference shows either an A or a G in a read mapping on the negative + // strand. + std::uint64_t error_count_negative{}; - // converted_count is equal to the sum of converted_count_positive - // and converted_count_negative. - uint64_t converted_count{}; // total_count is equal to the sum of total_count_positive and // total_count_negative - uint64_t total_count{}; + [[nodiscard]] std::uint64_t + total_count() const { + return total_count_positive + total_count_negative; + } + + // converted_count is equal to the sum of converted_count_positive + // and converted_count_negative. + [[nodiscard]] std::uint64_t + converted_count() const { + return converted_count_positive + converted_count_negative; + } + + // bisulfite_conversion_rate_positive is equal to converted_count_positive + // divided by total_count_positive, a value that is always between 0 and + // 1. When total_count_positive is 0, then + // bisulfite_conversion_rate_positive is given a value of 0. + [[nodiscard]] double + bisulfite_conversion_rate_positive() const { + return static_cast(converted_count_positive) / + std::max(total_count_positive, static_cast(1)); + } + + // bisulfite_conversion_rate_negative is equal to converted_count_negative + // divided by total_count_negative, a value that is always between 0 and + // 1. When total_count_negative is 0, then + // bisulfite_conversion_rate_negative is given a value of 0. + [[nodiscard]] double + bisulfite_conversion_rate_negative() const { + return static_cast(converted_count_negative) / + std::max(total_count_negative, static_cast(1)); + } + // bisulfite_conversion_rate is equal to converted_count divided by // total_count, a value that is always between 0 and 1. When // total_count is 0, then bisulfite_conversion_rate is given a value // of 0. - double bisulfite_conversion_rate{}; + [[nodiscard]] double + bisulfite_conversion_rate() const { + return static_cast(converted_count()) / + std::max(total_count(), static_cast(1)); + } + + // error_count is equal to the sum of error_count_positive and + // error_count_negative + [[nodiscard]] std::uint64_t + error_count() const { + return error_count_positive + error_count_negative; + } - // error_count is the number of nucleotides covering a cytosine in - // the reference shows either an A or a G in the read. - uint64_t error_count{}; // valid_count is the number of nucleotides covering a cytosine in // the reference shows any nucleotide that is not an N in the read. - uint64_t valid_count{}; + [[nodiscard]] std::uint64_t + valid_count() const { + return total_count() + error_count(); + } + // error_rate is equal to error_count divided by valid_count, and is // a value that is always between 0 and 1. When valid_count is 0, // then error_rate is given a value of 0. - double error_rate{}; + [[nodiscard]] double + error_rate() const { + return static_cast(error_count()) / + std::max(valid_count(), static_cast(1)); + } void update_pos(const char nt) { @@ -111,7 +142,7 @@ struct bsrate_summary { converted_count_positive += (nt == 'T'); } else if (nt != 'N') - ++error_count; + ++error_count_positive; } void @@ -121,55 +152,21 @@ struct bsrate_summary { converted_count_negative += (nt == 'T'); } else if (nt != 'N') - ++error_count; + ++error_count_negative; } bsrate_summary & operator+=(const bsrate_summary &rhs) { - // ADS: the "rates" are set to 0.0 here to ensure that they are - // computed properly using the full data after any accumulation of - // integer values has completed. converted_count_positive += rhs.converted_count_positive; total_count_positive += rhs.total_count_positive; - bisulfite_conversion_rate_positive = 0.0; - converted_count_negative += rhs.converted_count_negative; total_count_negative += rhs.total_count_negative; - bisulfite_conversion_rate_negative = 0.0; - - converted_count += rhs.converted_count; - total_count += rhs.total_count; - bisulfite_conversion_rate = 0.0; - - error_count += rhs.error_count; - valid_count += rhs.valid_count; - error_rate = 0.0; + error_count_positive += rhs.error_count_positive; + error_count_negative += rhs.error_count_negative; return *this; } - void - set_values() { - bisulfite_conversion_rate_positive = - static_cast(converted_count_positive) / - max(total_count_positive, static_cast(1)); - - bisulfite_conversion_rate_negative = - static_cast(converted_count_negative) / - max(total_count_negative, static_cast(1)); - - converted_count = converted_count_positive + converted_count_negative; - total_count = total_count_positive + total_count_negative; - - bisulfite_conversion_rate = static_cast(converted_count) / - max(total_count, static_cast(1)); - - valid_count = total_count + error_count; - - error_rate = static_cast(error_count) / - max(valid_count, static_cast(1)); - } - - string + std::string tostring_as_row() const { static constexpr auto precision_val = 5u; std::ostringstream oss; @@ -178,22 +175,22 @@ struct bsrate_summary { // clang-format off oss << total_count_positive << '\t' << converted_count_positive << '\t' - << bisulfite_conversion_rate_positive << '\t'; + << bisulfite_conversion_rate_positive() << '\t'; oss << total_count_negative << '\t' << converted_count_negative << '\t' - << bisulfite_conversion_rate_negative << '\t'; - oss << total_count << '\t' - << converted_count << '\t' - << bisulfite_conversion_rate << '\t'; - oss << error_count << '\t' - << valid_count << '\t' - << error_rate; + << bisulfite_conversion_rate_negative() << '\t'; + oss << total_count() << '\t' + << converted_count() << '\t' + << bisulfite_conversion_rate() << '\t'; + oss << error_count() << '\t' + << valid_count() << '\t' + << error_rate(); // clang-format on return oss.str(); } - string - tostring_as_yaml_list(const uint32_t position) const { + std::string + tostring_as_yaml_list(const std::uint32_t position) const { static constexpr auto precision_val = 5u; std::ostringstream oss; oss.precision(precision_val); @@ -202,21 +199,21 @@ struct bsrate_summary { oss << " - base: " << position << '\n' << " ptot: " << total_count_positive << '\n' << " pconv: " << converted_count_positive << '\n' - << " prate: " << bisulfite_conversion_rate_positive << '\n' + << " prate: " << bisulfite_conversion_rate_positive() << '\n' << " ntot: " << total_count_negative << '\n' << " nconv: " << converted_count_negative << '\n' - << " nrate: " << bisulfite_conversion_rate_negative << '\n' - << " bthtot: " << total_count << '\n' - << " bthconv: " << converted_count << '\n' - << " bthrate: " << bisulfite_conversion_rate << '\n' - << " err: " << error_count << '\n' - << " all: " << valid_count << '\n' - << " errrate: " << error_rate << '\n'; + << " nrate: " << bisulfite_conversion_rate_negative() << '\n' + << " bthtot: " << total_count() << '\n' + << " bthconv: " << converted_count() << '\n' + << " bthrate: " << bisulfite_conversion_rate() << '\n' + << " err: " << error_count() << '\n' + << " all: " << valid_count() << '\n' + << " errrate: " << error_rate() << '\n'; // clang-format on return oss.str(); } - string + std::string tostring() const { static constexpr auto sep = ": "; std::ostringstream oss; @@ -226,7 +223,7 @@ struct bsrate_summary { oss << "total_count_positive" << sep << total_count_positive << '\n'; oss << "bisulfite_conversion_rate_positive" << sep - << bisulfite_conversion_rate_positive << '\n'; + << bisulfite_conversion_rate_positive() << '\n'; oss << "converted_count_negative" << sep << converted_count_negative << '\n'; @@ -234,17 +231,17 @@ struct bsrate_summary { oss << "total_count_negative" << sep << total_count_negative << '\n'; oss << "bisulfite_conversion_rate_negative" << sep - << bisulfite_conversion_rate_negative << '\n'; + << bisulfite_conversion_rate_negative() << '\n'; - oss << "converted_count" << sep << converted_count << '\n'; - oss << "total_count" << sep << total_count << '\n'; + oss << "converted_count" << sep << converted_count() << '\n'; + oss << "total_count" << sep << total_count() << '\n'; - oss << "bisulfite_conversion_rate" << sep << bisulfite_conversion_rate + oss << "bisulfite_conversion_rate" << sep << bisulfite_conversion_rate() << '\n'; - oss << "error_count" << sep << error_count << '\n'; - oss << "valid_count" << sep << valid_count << '\n'; - oss << "error_rate" << sep << error_rate; + oss << "error_count" << sep << error_count() << '\n'; + oss << "valid_count" << sep << valid_count() << '\n'; + oss << "error_rate" << sep << error_rate(); return oss.str(); } @@ -256,11 +253,11 @@ operator+(bsrate_summary lhs, const bsrate_summary &rhs) { return lhs; } -static pair -count_states_pos(const bool INCLUDE_CPGS, const string &chrom, - const bam_rec &aln, vector &summaries, - size_t &hanging) { - uint32_t n_conv = 0, n_uconv = 0; +static std::pair +count_states_pos(const bool INCLUDE_CPGS, const std::string &chrom, + const bamxx::bam_rec &aln, + std::vector &summaries, std::size_t &hanging) { + std::uint32_t n_conv = 0, n_uconv = 0; /* iterate through reference, query/read and fragment */ const auto seq = bam_get_seq(aln); @@ -278,7 +275,8 @@ count_states_pos(const bool INCLUDE_CPGS, const string &chrom, if (cigar_eats_ref(op) && cigar_eats_query(op)) { const decltype(qpos) end_qpos = qpos + n; for (; qpos < end_qpos; ++qpos, ++rpos, ++fpos) { - if (rpos > chrom_lim) ++hanging; + if (rpos > chrom_lim) + ++hanging; if (is_cytosine(chrom[rpos]) && (rpos >= chrom_lim || !is_guanine(chrom[rpos + 1]) || INCLUDE_CPGS)) { @@ -290,20 +288,23 @@ count_states_pos(const bool INCLUDE_CPGS, const string &chrom, } } else { - if (cigar_eats_query(op)) qpos += n; - if (cigar_eats_ref(op)) rpos += n; - if (cigar_eats_frag(op)) fpos += n; + if (cigar_eats_query(op)) + qpos += n; + if (cigar_eats_ref(op)) + rpos += n; + if (cigar_eats_frag(op)) + fpos += n; } } assert(qpos == get_l_qseq(aln)); return {n_conv, n_conv + n_uconv}; } -static pair -count_states_neg(const bool INCLUDE_CPGS, const string &chrom, - const bam_rec &aln, vector &summaries, - size_t &hanging) { - uint32_t n_conv = 0, n_uconv = 0; +static std::pair +count_states_neg(const bool INCLUDE_CPGS, const std::string &chrom, + const bamxx::bam_rec &aln, + std::vector &summaries, std::size_t &hanging) { + std::uint32_t n_conv = 0, n_uconv = 0; /* iterate backward over query/read positions but forward over reference and fragment positions */ @@ -322,7 +323,8 @@ count_states_neg(const bool INCLUDE_CPGS, const string &chrom, if (cigar_eats_ref(op) && cigar_eats_query(op)) { const decltype(qpos) end_qpos = qpos - n; for (; qpos > end_qpos; --qpos, ++rpos, ++fpos) { - if (rpos > chrom_lim) ++hanging; + if (rpos > chrom_lim) + ++hanging; if (is_guanine(chrom[rpos]) && (rpos == 0 || !is_cytosine(chrom[rpos - 1]) || INCLUDE_CPGS)) { const auto nt = seq_nt16_str[bam_seqi(seq, qpos - 1)]; @@ -333,9 +335,12 @@ count_states_neg(const bool INCLUDE_CPGS, const string &chrom, } } else { - if (cigar_eats_query(op)) qpos -= n; - if (cigar_eats_ref(op)) rpos += n; - if (cigar_eats_frag(op)) fpos += n; + if (cigar_eats_query(op)) + qpos -= n; + if (cigar_eats_ref(op)) + rpos += n; + if (cigar_eats_frag(op)) + fpos += n; } } assert(qpos == 0); @@ -343,22 +348,25 @@ count_states_neg(const bool INCLUDE_CPGS, const string &chrom, } static void -write_output(const string &outfile, const vector &summaries) { +write_output(const std::string &outfile, + const std::vector &summaries) { std::ofstream of; - if (!outfile.empty()) of.open(outfile.c_str()); + if (!outfile.empty()) + of.open(outfile.c_str()); std::ostream out(outfile.empty() ? std::cout.rdbuf() : of.rdbuf()); - if (!out) throw dnmt_error("failed to open output file"); + if (!out) + throw dnmt_error("failed to open output file"); - bsrate_summary overall_summary = reduce(cbegin(summaries), cend(summaries)); - overall_summary.set_values(); + const bsrate_summary overall_summary = + std::reduce(std::cbegin(summaries), std::cend(summaries)); out << "OVERALL CONVERSION RATE = " - << overall_summary.bisulfite_conversion_rate << '\n' + << overall_summary.bisulfite_conversion_rate() << '\n' << "POS CONVERSION RATE = " - << overall_summary.bisulfite_conversion_rate_positive << '\t' + << overall_summary.bisulfite_conversion_rate_positive() << '\t' << overall_summary.total_count_positive << '\n' << "NEG CONVERSION RATE = " - << overall_summary.bisulfite_conversion_rate_negative << '\t' + << overall_summary.bisulfite_conversion_rate_negative() << '\t' << overall_summary.total_count_negative << '\n'; // clang-format off @@ -374,12 +382,12 @@ write_output(const string &outfile, const vector &summaries) { << "BTHRATE" << '\t' << "ERR" << '\t' << "ALL" << '\t' - << "ERRRATE" << endl; + << "ERRRATE" << '\n'; // clang-format on // figure out how many positions to print in the output, capped at 1000 auto output_len = std::min(std::size(summaries), 1000ul); - while (output_len > 0 && summaries[output_len - 1].total_count == 0) + while (output_len > 0 && summaries[output_len - 1].total_count() == 0) --output_len; for (auto i = 0u; i < output_len; ++i) @@ -387,27 +395,30 @@ write_output(const string &outfile, const vector &summaries) { } static void -write_output_yaml(const string &outfile, - const vector &summaries) { +write_output_yaml(const std::string &outfile, + const std::vector &summaries) { std::ofstream of; - if (!outfile.empty()) of.open(outfile.c_str()); + if (!outfile.empty()) + of.open(outfile.c_str()); std::ostream out(outfile.empty() ? std::cout.rdbuf() : of.rdbuf()); - if (!out) throw dnmt_error("failed to open output file"); + if (!out) + throw dnmt_error("failed to open output file"); - bsrate_summary overall_summary = reduce(cbegin(summaries), cend(summaries)); - overall_summary.set_values(); + const bsrate_summary overall_summary = + std::reduce(std::cbegin(summaries), std::cend(summaries)); + // overall_summary.set_values(); out << "overall_conversion_rate: " - << overall_summary.bisulfite_conversion_rate << '\n' + << overall_summary.bisulfite_conversion_rate() << '\n' << "pos_conversion_rate: " - << overall_summary.bisulfite_conversion_rate_positive << '\n' + << overall_summary.bisulfite_conversion_rate_positive() << '\n' << "pos_count: " << overall_summary.total_count_positive << '\n' << "neg_conversion_rate: " - << overall_summary.bisulfite_conversion_rate_negative << '\n' + << overall_summary.bisulfite_conversion_rate_negative() << '\n' << "neg_count: " << overall_summary.total_count_negative << '\n'; // figure out how many positions to print in the output, capped at 1000 auto output_len = std::min(std::size(summaries), 1000ul); - while (output_len > 0 && summaries[output_len - 1].total_count == 0) + while (output_len > 0 && summaries[output_len - 1].total_count() == 0) --output_len; out << "rows:\n"; @@ -416,28 +427,31 @@ write_output_yaml(const string &outfile, } static void -write_summary(const string &summary_file, vector &summaries) { - auto s = reduce(cbegin(summaries), cend(summaries)); - s.set_values(); - +write_summary(const std::string &summary_file, + std::vector &summaries) { + const auto s = std::reduce(std::cbegin(summaries), std::cend(summaries)); std::ofstream out(summary_file); - if (!out) throw dnmt_error("failed to open file: " + summary_file); + if (!out) + throw dnmt_error("failed to open file: " + summary_file); out << s.tostring() << '\n'; } -template static inline void -update_per_read_stats(const pair &x, vector> &tab) { - if (x.second < std::size(tab)) ++tab[x.second][x.first]; +template +static inline void +update_per_read_stats(const std::pair &x, + std::vector> &tab) { + if (x.second < std::size(tab)) + ++tab[x.second][x.first]; } -static inline vector -format_histogram(const vector> &tab, - const size_t n_hist_bins) { +static inline std::vector +format_histogram(const std::vector> &tab, + const std::size_t n_hist_bins) { static constexpr auto epsilon = 1e-6; - vector hist(n_hist_bins, 0.0); - for (size_t i = 1; i < tab.size(); ++i) { + std::vector hist(n_hist_bins, 0.0); + for (std::size_t i = 1; i < tab.size(); ++i) { const double denom = i + epsilon; - for (size_t j = 0; j <= i; ++j) { + for (std::size_t j = 0; j <= i; ++j) { const double frac = j / denom; const auto bin_id = std::floor(frac * n_hist_bins); assert(bin_id < hist.size()); @@ -448,48 +462,49 @@ format_histogram(const vector> &tab, } static void -write_hanging_read_message(std::ostream &out, const size_t n_hanging) { +write_hanging_read_message(std::ostream &out, const std::size_t n_hanging) { out << "Warning: hanging reads detected at chromosome ends " - << "(N=" << n_hanging << ")." << endl - << "High numbers of hanging reads suggest inconsistent " << endl - << "reference genomes between stages of analysis." << endl - << "This is likely to result in analysis errors." << endl; + << "(N=" << n_hanging << ")." << '\n' + << "High numbers of hanging reads suggest inconsistent " << '\n' + << "reference genomes between stages of analysis." << '\n' + << "This is likely to result in analysis errors." << '\n'; } -template static void -write_per_read_histogram(const vector> &tab, const size_t n_hist_bins, - std::ostream &out) { +template +static void +write_per_read_histogram(const std::vector> &tab, + const std::size_t n_hist_bins, std::ostream &out) { const auto hist = format_histogram(tab, n_hist_bins); out << std::fixed; for (auto i = 0.0; i < hist.size(); ++i) out << std::setprecision(3) << i / hist.size() << '\t' << std::setprecision(3) << (i + 1) / hist.size() << '\t' - << std::setprecision(0) << hist[i] << endl; + << std::setprecision(0) << hist[i] << '\n'; } int main_bsrate(int argc, char *argv[]) { try { // assumed maximum length of a fragment - static constexpr const size_t output_size = 10000; + static constexpr const std::size_t output_size = 10000; // Assumed maximum cytosines per fragment. Currently the per-read // information is collected as counts in a 2D array, so that later // features may be able to fit distributions based these counts. - static constexpr const size_t max_cytosine_per_frag = 1000; + static constexpr const std::size_t max_cytosine_per_frag = 1000; bool VERBOSE = false; bool INCLUDE_CPGS = false; bool output_yaml = false; bool reads_are_a_rich = false; bool report_per_read = false; - size_t n_threads = 1; - size_t n_hist_bins = 20; + std::size_t n_threads = 1; + std::size_t n_hist_bins = 20; - string chroms_file; - string summary_file; - string outfile; - string seq_to_use; // use only this chrom/sequence in the analysis + std::string chroms_file; + std::string summary_file; + std::string outfile; + std::string seq_to_use; // use only this chrom/sequence in the analysis /****************** COMMAND LINE OPTIONS ********************/ OptionParser opt_parse(strip_path(argv[0]), @@ -515,82 +530,88 @@ main_bsrate(int argc, char *argv[]) { opt_parse.add_opt("yaml", 'y', "output in yaml format", false, output_yaml); opt_parse.add_opt("threads", 't', "number of threads", false, n_threads); opt_parse.add_opt("verbose", 'v', "print more run info", false, VERBOSE); - vector leftover_args; + std::vector leftover_args; opt_parse.parse(argc, argv, leftover_args); if (argc == 1 || opt_parse.help_requested()) { - cerr << opt_parse.help_message() << endl - << opt_parse.about_message() << endl; + std::cerr << opt_parse.help_message() << '\n' + << opt_parse.about_message() << '\n'; return EXIT_SUCCESS; } if (opt_parse.about_requested()) { - cerr << opt_parse.about_message() << endl; + std::cerr << opt_parse.about_message() << '\n'; return EXIT_SUCCESS; } if (opt_parse.option_missing()) { - cerr << opt_parse.option_missing_message() << endl; + std::cerr << opt_parse.option_missing_message() << '\n'; return EXIT_SUCCESS; } if (leftover_args.size() != 1) { - cerr << opt_parse.help_message() << endl; + std::cerr << opt_parse.help_message() << '\n'; return EXIT_SUCCESS; } - const string bam_file = leftover_args.front(); + const std::string bam_file = leftover_args.front(); /****************** END COMMAND LINE OPTIONS *****************/ - vector chroms; - vector names; + std::vector chroms; + std::vector names; read_fasta_file_short_names(chroms_file, names, chroms); for (auto &&chrom : chroms) - std::for_each(begin(chrom), end(chrom), + std::for_each(std::begin(chrom), std::end(chrom), [](const char c) { return std::toupper(c); }); if (VERBOSE) - cerr << "[n chroms in reference: " << chroms.size() << "]" << endl; + std::cerr << "[n chroms in reference: " << chroms.size() << "]" << '\n'; bamxx::bam_tpool tp(n_threads); bamxx::bam_in hts(bam_file); - if (!hts) throw dnmt_error("failed to open input file: " + bam_file); + if (!hts) + throw dnmt_error("failed to open input file: " + bam_file); bamxx::bam_header hdr(hts); - if (!hdr) throw dnmt_error("failed to read header"); + if (!hdr) + throw dnmt_error("failed to read header"); - if (n_threads > 1) tp.set_io(hts); + if (n_threads > 1) + tp.set_io(hts); // map the bam header index for each "target" to a sequence in the // reference genome - unordered_map chrom_lookup; - size_t chrom_idx_to_use = numeric_limits::max(); + std::unordered_map chrom_lookup; + std::size_t chrom_idx_to_use = std::numeric_limits::max(); for (auto i = 0u; i < std::size(chroms); ++i) { - if (names[i] == seq_to_use) chrom_idx_to_use = i; + if (names[i] == seq_to_use) + chrom_idx_to_use = i; chrom_lookup.emplace(sam_hdr_name2tid(hdr.h, names[i].data()), i); } - vector> per_read_counts( - max_cytosine_per_frag, vector(max_cytosine_per_frag, 0)); + std::vector> per_read_counts( + max_cytosine_per_frag, + std::vector(max_cytosine_per_frag, 0)); - vector summaries(output_size); + std::vector summaries(output_size); - int32_t current_tid = -1; - size_t chrom_idx = numeric_limits::max(); - size_t hanging = 0; + std::int32_t current_tid = -1; + std::size_t chrom_idx = std::numeric_limits::max(); + std::size_t hanging = 0; bool use_this_chrom = seq_to_use.empty(); - bam_rec aln; - unordered_set chroms_seen; + bamxx::bam_rec aln; + std::unordered_set chroms_seen; while (hts.read(hdr, aln)) { - const int32_t the_tid = get_tid(aln); + const std::int32_t the_tid = get_tid(aln); if (the_tid == -1) continue; - if (reads_are_a_rich) flip_conversion(aln); + if (reads_are_a_rich) + flip_conversion(aln); // get the correct chrom if it has changed if (the_tid != current_tid) { // make sure all reads from same chrom are contiguous in the file if (chroms_seen.find(the_tid) != end(chroms_seen)) - throw runtime_error("chroms out of order in mapped reads file"); + throw std::runtime_error("chroms out of order in mapped reads file"); current_tid = the_tid; @@ -598,12 +619,13 @@ main_bsrate(int argc, char *argv[]) { auto chrom_itr = chrom_lookup.find(the_tid); if (chrom_itr == end(chrom_lookup)) - throw runtime_error("could not find chrom: " + - std::to_string(the_tid)); + throw std::runtime_error("could not find chrom: " + + std::to_string(the_tid)); chrom_idx = chrom_itr->second; - if (VERBOSE) cerr << "processing " << names[chrom_idx] << endl; + if (VERBOSE) + std::cerr << "processing " << names[chrom_idx] << '\n'; use_this_chrom = seq_to_use.empty() || chrom_idx == chrom_idx_to_use; } @@ -620,25 +642,22 @@ main_bsrate(int argc, char *argv[]) { } } - // before writing the output we need to ensure the derived - // quantities in bsrate_summary have been calculated - for_each(begin(summaries), end(summaries), - [](bsrate_summary &s) { s.set_values(); }); - if (output_yaml) write_output_yaml(outfile, summaries); else write_output(outfile, summaries); - if (hanging > 0) write_hanging_read_message(cerr, hanging); + if (hanging > 0) + write_hanging_read_message(std::cerr, hanging); - if (!summary_file.empty()) write_summary(summary_file, summaries); + if (!summary_file.empty()) + write_summary(summary_file, summaries); if (report_per_read) - write_per_read_histogram(per_read_counts, n_hist_bins, cout); + write_per_read_histogram(per_read_counts, n_hist_bins, std::cout); } - catch (const runtime_error &e) { - cerr << e.what() << endl; + catch (const std::runtime_error &e) { + std::cerr << e.what() << '\n'; return EXIT_FAILURE; } return EXIT_SUCCESS; From 90d19ce234dd97d83594db3468c030ced090e378 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Tue, 7 Oct 2025 13:05:16 -0700 Subject: [PATCH 03/11] src/common/LevelsCounter.hcpp: made functions for mean_depth, mean_depth_covered and sites_covered_fraction --- src/common/LevelsCounter.cpp | 57 ++++++++++++++++++------------------ src/common/LevelsCounter.hpp | 21 +++++++++++++ 2 files changed, 49 insertions(+), 29 deletions(-) diff --git a/src/common/LevelsCounter.cpp b/src/common/LevelsCounter.cpp index 077051de..373d5403 100644 --- a/src/common/LevelsCounter.cpp +++ b/src/common/LevelsCounter.cpp @@ -16,13 +16,13 @@ #include "LevelsCounter.hpp" #include "bsutils.hpp" -#include #include #include +#include +using std::runtime_error; using std::string; using std::to_string; -using std::runtime_error; void LevelsCounter::update(const MSite &s) { @@ -43,9 +43,9 @@ LevelsCounter::update(const MSite &s) { ++total_sites; } -string +std::string LevelsCounter::tostring() const { - static const string indent = string(2, ' '); + static const std::string indent = std::string(2, ' '); const bool good = (sites_covered != 0); std::ostringstream oss; // directly counted values @@ -62,22 +62,20 @@ LevelsCounter::tostring() const { // derived values oss << indent << "coverage: " << coverage() << '\n' - << indent << "sites_covered_fraction: " - << static_cast(sites_covered)/total_sites << '\n' - << indent << "mean_depth: " - << static_cast(coverage())/total_sites << '\n' - << indent << "mean_depth_covered: " - << static_cast(coverage())/sites_covered << '\n' - << indent << "mean_meth: " - << (good ? to_string(mean_meth()) : "NA") << '\n' + << indent << "sites_covered_fraction: " << sites_covered_fraction() + << '\n' + << indent << "mean_depth: " << mean_depth() << '\n' + << indent << "mean_depth_covered: " << mean_depth_covered() << '\n' + << indent << "mean_meth: " << (good ? std::to_string(mean_meth()) : "NA") + << '\n' << indent << "mean_meth_weighted: " - << (good ? to_string(mean_meth_weighted()) : "NA") << '\n' + << (good ? std::to_string(mean_meth_weighted()) : "NA") << '\n' << indent << "fractional_meth: " - << (good ? to_string(fractional_meth()) : "NA"); + << (good ? std::to_string(fractional_meth()) : "NA"); return oss.str(); } -string +std::string LevelsCounter::format_row() const { const bool good = (sites_covered > 0); std::ostringstream oss; @@ -104,7 +102,7 @@ LevelsCounter::format_row() const { return oss.str(); } -string +std::string LevelsCounter::format_header() { std::ostringstream oss; // directly counted values @@ -136,42 +134,43 @@ operator<<(std::ostream &out, const LevelsCounter &cs) { } static void -check_label(const string &observed, const string expected) { +check_label(const std::string &observed, const std::string expected) { if (observed != expected) - throw runtime_error("bad levels format [" + observed + "," + expected + "]"); + throw runtime_error("bad levels format [" + observed + "," + expected + + "]"); } std::istream & operator>>(std::istream &in, LevelsCounter &cs) { - in >> cs.context; // get the context + in >> cs.context; // get the context cs.context = cs.context.substr(0, cs.context.find_first_of(":")); - string label; - in >> label >> cs.total_sites; // the total sites + std::string label; + in >> label >> cs.total_sites; // the total sites check_label(label, "total_sites:"); - in >> label >> cs.sites_covered; // the sites covered + in >> label >> cs.sites_covered; // the sites covered check_label(label, "sites_covered:"); - in >> label >> cs.total_c; // the total c + in >> label >> cs.total_c; // the total c check_label(label, "total_c:"); - in >> label >> cs.total_t; // the total t + in >> label >> cs.total_t; // the total t check_label(label, "total_t:"); - in >> label >> cs.max_depth; // the max depth + in >> label >> cs.max_depth; // the max depth check_label(label, "max_depth:"); - in >> label >> cs.mutations; // the number of mutations + in >> label >> cs.mutations; // the number of mutations check_label(label, "mutations:"); - in >> label >> cs.called_meth; // the number of sites called methylated + in >> label >> cs.called_meth; // the number of sites called methylated check_label(label, "called_meth:"); - in >> label >> cs.called_unmeth; // the number of sites called unmethylated + in >> label >> cs.called_unmeth; // the number of sites called unmethylated check_label(label, "called_unmeth:"); - in >> label >> cs.total_meth; // the mean aggregate + in >> label >> cs.total_meth; // the mean aggregate check_label(label, "total_meth:"); return in; diff --git a/src/common/LevelsCounter.hpp b/src/common/LevelsCounter.hpp index b264a66c..512c5794 100644 --- a/src/common/LevelsCounter.hpp +++ b/src/common/LevelsCounter.hpp @@ -76,6 +76,27 @@ struct LevelsCounter { return total_c + total_t; } + // mean_depth is the average number of reads contributing to methylation + // calls over all sites. + [[nodiscard]] double + mean_depth() const { + return static_cast(coverage()) / total_sites; + } + + // mean_depth_covered is the average number of reads contributing to + // methylation calls over all sites that are covered at least once. + [[nodiscard]] double + mean_depth_covered() const { + return static_cast(coverage()) / sites_covered; + } + + // sites_covered_fraction is the fraction of all sites that are covered at + // least once. + [[nodiscard]] double + sites_covered_fraction() const { + return static_cast(sites_covered) / total_sites; + } + // total_called is equal to called_meth plus called_unmeth [[nodiscard]] std::uint64_t total_called() const { From cfddf690112171ae3dc3cb307d26d5ab3baf7b54 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Tue, 7 Oct 2025 13:06:57 -0700 Subject: [PATCH 04/11] src/abismal: updating submodule --- src/abismal | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/abismal b/src/abismal index ad4832c4..40548a0a 160000 --- a/src/abismal +++ b/src/abismal @@ -1 +1 @@ -Subproject commit ad4832c43ed8c994e1de3b42c702ecffe7515a7f +Subproject commit 40548a0a6a1a38fde9e66d67a6f84225c6f1cd30 From bf18742d41d5edee416990425bdab2cb22d12a78 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Tue, 7 Oct 2025 13:09:45 -0700 Subject: [PATCH 05/11] src/common/LevelsCounter.hcpp: renamed sites_covered_fraction to sites_covered_frac to match the colname in the view in methbase db --- src/common/LevelsCounter.cpp | 3 +-- src/common/LevelsCounter.hpp | 6 +++--- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/src/common/LevelsCounter.cpp b/src/common/LevelsCounter.cpp index 373d5403..de9d7109 100644 --- a/src/common/LevelsCounter.cpp +++ b/src/common/LevelsCounter.cpp @@ -62,8 +62,7 @@ LevelsCounter::tostring() const { // derived values oss << indent << "coverage: " << coverage() << '\n' - << indent << "sites_covered_fraction: " << sites_covered_fraction() - << '\n' + << indent << "sites_covered_frac: " << sites_covered_frac() << '\n' << indent << "mean_depth: " << mean_depth() << '\n' << indent << "mean_depth_covered: " << mean_depth_covered() << '\n' << indent << "mean_meth: " << (good ? std::to_string(mean_meth()) : "NA") diff --git a/src/common/LevelsCounter.hpp b/src/common/LevelsCounter.hpp index 512c5794..2ce520fc 100644 --- a/src/common/LevelsCounter.hpp +++ b/src/common/LevelsCounter.hpp @@ -90,10 +90,10 @@ struct LevelsCounter { return static_cast(coverage()) / sites_covered; } - // sites_covered_fraction is the fraction of all sites that are covered at - // least once. + // sites_covered_frac is the fraction of all sites that are covered at least + // once. [[nodiscard]] double - sites_covered_fraction() const { + sites_covered_frac() const { return static_cast(sites_covered) / total_sites; } From c0f37284c22a98e47b73918bc414036530171d9e Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Tue, 7 Oct 2025 14:12:51 -0700 Subject: [PATCH 06/11] src/analysis/hmr.cpp: cleanup --- src/analysis/hmr.cpp | 363 +++++++++++++++++++++---------------------- 1 file changed, 178 insertions(+), 185 deletions(-) diff --git a/src/analysis/hmr.cpp b/src/analysis/hmr.cpp index dda0c230..9fac1574 100644 --- a/src/analysis/hmr.cpp +++ b/src/analysis/hmr.cpp @@ -14,6 +14,17 @@ * GNU General Public License for more details. */ +#include "MSite.hpp" +#include "TwoStateHMM.hpp" +#include "counts_header.hpp" + +#include "GenomicRegion.hpp" +#include "OptionParser.hpp" +#include "smithlab_os.hpp" +#include "smithlab_utils.hpp" + +#include + #include #include // for [u]int[0-9]+_t #include @@ -24,57 +35,29 @@ #include #include -#include - -#include "GenomicRegion.hpp" -#include "OptionParser.hpp" -#include "smithlab_os.hpp" -#include "smithlab_utils.hpp" - -#include "MSite.hpp" -#include "TwoStateHMM.hpp" -#include "counts_header.hpp" - -using std::accumulate; -using std::cerr; -using std::cout; -using std::endl; -using std::make_pair; -using std::max; -using std::min; -using std::numeric_limits; -using std::pair; -using std::reduce; -using std::runtime_error; -using std::string; -using std::to_string; -using std::unordered_set; -using std::vector; - -using bamxx::bgzf_file; - struct hmr_summary { - hmr_summary(const vector &hmrs) { + hmr_summary(const std::vector &hmrs) { hmr_count = hmrs.size(); - hmr_total_size = accumulate(cbegin(hmrs), cend(hmrs), 0ul, - [](const uint64_t t, const GenomicRegion &p) { - return t + p.get_width(); - }); + hmr_total_size = + std::accumulate(std::cbegin(hmrs), std::cend(hmrs), 0ul, + [](const std::uint64_t t, const GenomicRegion &p) { + return t + p.get_width(); + }); hmr_mean_size = static_cast(hmr_total_size) / - std::max(hmr_count, static_cast(1)); + std::max(hmr_count, static_cast(1)); } // hmr_count is the number of identified HMRs. - uint64_t hmr_count{}; + std::uint64_t hmr_count{}; // total_hmr_size is the sum of the sizes of the identified HMRs - uint64_t hmr_total_size{}; + std::uint64_t hmr_total_size{}; // mean_hmr_size is the mean size of the identified HMRs double hmr_mean_size{}; - string + std::string tostring() { std::ostringstream oss; - oss << "hmr_count: " << hmr_count << endl - << "hmr_total_size: " << hmr_total_size << endl + oss << "hmr_count: " << hmr_count << '\n' + << "hmr_total_size: " << hmr_total_size << '\n' << "hmr_mean_size: " << std::fixed << std::setprecision(2) << hmr_mean_size; @@ -87,36 +70,37 @@ as_gen_rgn(const MSite &s) { return GenomicRegion(s.chrom, s.pos, s.pos + 1); } -static string -format_cpg_meth_tag(const pair &m) { - return "CpG:" + to_string(static_cast(m.first)) + ":" + - to_string(static_cast(m.second)); +static std::string +format_cpg_meth_tag(const std::pair &m) { + return "CpG:" + std::to_string(static_cast(m.first)) + ":" + + std::to_string(static_cast(m.second)); } static double -get_stepup_cutoff(vector scores, const double cutoff) { +get_stepup_cutoff(std::vector scores, const double cutoff) { if (cutoff <= 0) - return numeric_limits::max(); + return std::numeric_limits::max(); else if (cutoff > 1) - return numeric_limits::min(); + return std::numeric_limits::min(); - const size_t n = scores.size(); - std::sort(begin(scores), end(scores)); - size_t i = 1; + const std::size_t n = scores.size(); + std::sort(begin(scores), std::end(scores)); + std::size_t i = 1; while (i < n && scores[i - 1] < (cutoff * i) / n) ++i; return scores[i - 1]; } static void -get_domain_scores(const vector &state_ids, - const vector> &meth, - const vector &reset_points, vector &scores) { +get_domain_scores(const std::vector &state_ids, + const std::vector> &meth, + const std::vector &reset_points, + std::vector &scores) { - size_t reset_idx = 1; + std::size_t reset_idx = 1; bool in_domain = false; double score = 0; - for (size_t i = 0; i < state_ids.size(); ++i) { + for (std::size_t i = 0; i < state_ids.size(); ++i) { if (reset_points[reset_idx] == i) { if (in_domain) { in_domain = false; @@ -141,12 +125,14 @@ get_domain_scores(const vector &state_ids, } static void -build_domains(const vector &cpgs, const vector &reset_points, - const vector &state_ids, vector &domains) { +build_domains(const std::vector &cpgs, + const std::vector &reset_points, + const std::vector &state_ids, + std::vector &domains) { - size_t n_cpgs = 0, reset_idx = 1, prev_end = 0; + std::size_t n_cpgs = 0, reset_idx = 1, prev_end = 0; bool in_domain = false; - for (size_t i = 0; i < state_ids.size(); ++i) { + for (std::size_t i = 0; i < state_ids.size(); ++i) { if (reset_points[reset_idx] == i) { if (in_domain) { in_domain = false; @@ -160,7 +146,7 @@ build_domains(const vector &cpgs, const vector &reset_points, if (!in_domain) { in_domain = true; domains.push_back(as_gen_rgn(cpgs[i])); - domains.back().set_name("HYPO" + to_string(domains.size())); + domains.back().set_name("HYPO" + std::to_string(domains.size())); } ++n_cpgs; } @@ -180,38 +166,39 @@ build_domains(const vector &cpgs, const vector &reset_points, template static void -separate_regions(const bool verbose, const size_t desert_size, - vector &cpgs, vector &meth, vector &reads, - vector &reset_points) { +separate_regions(const bool verbose, const std::size_t desert_size, + std::vector &cpgs, std::vector &meth, + std::vector &reads, + std::vector &reset_points) { if (verbose) - cerr << "[separating by cpg desert]" << endl; + std::cerr << "[separating by cpg desert]" << '\n'; // eliminate the zero-read cpgs - size_t j = 0; - for (size_t i = 0; i < cpgs.size(); ++i) + std::size_t j = 0; + for (std::size_t i = 0; i < cpgs.size(); ++i) if (reads[i] > 0) { cpgs[j] = cpgs[i]; meth[j] = meth[i]; reads[j] = reads[i]; ++j; } - cpgs.erase(begin(cpgs) + j, end(cpgs)); - meth.erase(begin(meth) + j, end(meth)); - reads.erase(begin(reads) + j, end(reads)); + cpgs.erase(std::begin(cpgs) + j, std::end(cpgs)); + meth.erase(std::begin(meth) + j, std::end(meth)); + reads.erase(std::begin(reads) + j, std::end(reads)); double total_bases = 0; double bases_in_deserts = 0; // segregate cpgs - size_t prev_pos = 0; - for (size_t i = 0; i < cpgs.size(); ++i) { - const size_t dist = (i > 0 && cpgs[i].chrom == cpgs[i - 1].chrom) - ? cpgs[i].pos - prev_pos - : numeric_limits::max(); + std::size_t prev_pos = 0; + for (std::size_t i = 0; i < cpgs.size(); ++i) { + const std::size_t dist = (i > 0 && cpgs[i].chrom == cpgs[i - 1].chrom) + ? cpgs[i].pos - prev_pos + : std::numeric_limits::max(); if (dist > desert_size) { reset_points.push_back(i); - if (dist < numeric_limits::max()) + if (dist < std::numeric_limits::max()) bases_in_deserts += dist; } - if (dist < numeric_limits::max()) + if (dist < std::numeric_limits::max()) total_bases += dist; prev_pos = cpgs[i].pos; @@ -219,10 +206,10 @@ separate_regions(const bool verbose, const size_t desert_size, reset_points.push_back(cpgs.size()); if (verbose) - cerr << "[cpgs retained: " << cpgs.size() << "]" << endl - << "[deserts removed: " << reset_points.size() - 2 << "]" << endl - << "[genome fraction covered: " - << 1.0 - (bases_in_deserts / total_bases) << "]" << endl; + std::cerr << "[cpgs retained: " << cpgs.size() << "]" << '\n' + << "[deserts removed: " << reset_points.size() - 2 << "]" << '\n' + << "[genome fraction covered: " + << 1.0 - (bases_in_deserts / total_bases) << "]" << '\n'; } /* function to "fold" the methylation profile so that the middle @@ -230,9 +217,9 @@ separate_regions(const bool verbose, const size_t desert_size, * methylation become high. this method actually seems to work. */ static void -make_partial_meth(const vector &reads, - vector> &meth) { - for (size_t i = 0; i < reads.size(); ++i) { +make_partial_meth(const std::vector &reads, + std::vector> &meth) { + for (std::size_t i = 0; i < reads.size(); ++i) { double m = meth[i].first / reads[i]; m = (m <= 0.5) ? (1.0 - 2 * m) : (1.0 - 2 * (1.0 - m)); meth[i].first = reads[i] * m; @@ -241,58 +228,60 @@ make_partial_meth(const vector &reads, } static void -shuffle_cpgs(const size_t rng_seed, const TwoStateHMM &hmm, - vector> meth, vector reset_points, - const double p_fb, const double p_bf, const double fg_alpha, - const double fg_beta, const double bg_alpha, const double bg_beta, - vector &domain_scores) { +shuffle_cpgs(const std::size_t rng_seed, const TwoStateHMM &hmm, + std::vector> meth, + std::vector reset_points, const double p_fb, + const double p_bf, const double fg_alpha, const double fg_beta, + const double bg_alpha, const double bg_beta, + std::vector &domain_scores) { auto eng = std::default_random_engine(rng_seed); - std::shuffle(begin(meth), end(meth), eng); + std::shuffle(std::begin(meth), std::end(meth), eng); - vector state_ids; - vector scores; + std::vector state_ids; + std::vector scores; hmm.PosteriorDecoding(meth, reset_points, p_fb, p_bf, fg_alpha, fg_beta, bg_alpha, bg_beta, state_ids, scores); get_domain_scores(state_ids, meth, reset_points, domain_scores); - sort(begin(domain_scores), end(domain_scores)); + sort(std::begin(domain_scores), std::end(domain_scores)); } static void -assign_p_values(const vector &random_scores, - const vector &observed_scores, - vector &p_values) { +assign_p_values(const std::vector &random_scores, + const std::vector &observed_scores, + std::vector &p_values) { const double n_randoms = random_scores.empty() ? 1 : random_scores.size(); - for (size_t i = 0; i < observed_scores.size(); ++i) - p_values.push_back((end(random_scores) - upper_bound(begin(random_scores), - end(random_scores), - observed_scores[i])) / - n_randoms); + for (std::size_t i = 0; i < observed_scores.size(); ++i) + p_values.push_back( + (std::end(random_scores) - upper_bound(std::begin(random_scores), + std::end(random_scores), + observed_scores[i])) / + n_randoms); } static void -read_params_file(const bool verbose, const string ¶ms_file, +read_params_file(const bool verbose, const std::string ¶ms_file, double &fg_alpha, double &fg_beta, double &bg_alpha, double &bg_beta, double &p_fb, double &p_bf, double &domain_score_cutoff) { - string jnk; + std::string jnk; std::ifstream in(params_file); if (!in) - throw runtime_error("failed to parse params file: " + params_file); + throw std::runtime_error("failed to parse params file: " + params_file); in >> jnk >> fg_alpha >> jnk >> fg_beta >> jnk >> bg_alpha >> jnk >> bg_beta >> jnk >> p_fb >> jnk >> p_bf >> jnk >> domain_score_cutoff; if (verbose) - cerr << "FG_ALPHA\t" << fg_alpha << endl - << "FG_BETA\t" << fg_beta << endl - << "BG_ALPHA\t" << bg_alpha << endl - << "BG_BETA\t" << bg_beta << endl - << "F_B\t" << p_fb << endl - << "B_F\t" << p_bf << endl - << "DOMAIN_SCORE_CUTOFF\t" << domain_score_cutoff << endl; + std::cerr << "FG_ALPHA\t" << fg_alpha << '\n' + << "FG_BETA\t" << fg_beta << '\n' + << "BG_ALPHA\t" << bg_alpha << '\n' + << "BG_BETA\t" << bg_beta << '\n' + << "F_B\t" << p_fb << '\n' + << "B_F\t" << p_bf << '\n' + << "DOMAIN_SCORE_CUTOFF\t" << domain_score_cutoff << '\n'; } static void -write_params_file(const string &outfile, const double fg_alpha, +write_params_file(const std::string &outfile, const double fg_alpha, const double fg_beta, const double bg_alpha, const double bg_beta, const double p_fb, const double p_bf, const double domain_score_cutoff) { @@ -303,13 +292,13 @@ write_params_file(const string &outfile, const double fg_alpha, std::ostream out(outfile.empty() ? std::cout.rdbuf() : of.rdbuf()); out.precision(30); - out << "FG_ALPHA\t" << fg_alpha << endl - << "FG_BETA\t" << fg_beta << endl - << "BG_ALPHA\t" << bg_alpha << endl - << "BG_BETA\t" << bg_beta << endl - << "F_B\t" << p_fb << endl - << "B_F\t" << p_bf << endl - << "DOMAIN_SCORE_CUTOFF\t" << domain_score_cutoff << endl; + out << "FG_ALPHA\t" << fg_alpha << '\n' + << "FG_BETA\t" << fg_beta << '\n' + << "BG_ALPHA\t" << bg_alpha << '\n' + << "BG_BETA\t" << bg_beta << '\n' + << "F_B\t" << p_fb << '\n' + << "B_F\t" << p_bf << '\n' + << "DOMAIN_SCORE_CUTOFF\t" << domain_score_cutoff << '\n'; } [[nodiscard]] @@ -325,12 +314,13 @@ get_n_unmeth(const MSite &s) { } static void -load_cpgs(const string &cpgs_file, vector &cpgs, - vector> &meth, vector &reads) { +load_cpgs(const std::string &cpgs_file, std::vector &cpgs, + std::vector> &meth, + std::vector &reads) { - bgzf_file in(cpgs_file, "r"); + bamxx::bgzf_file in(cpgs_file, "r"); if (!in) - throw runtime_error("failed opening file: " + cpgs_file); + throw std::runtime_error("failed opening file: " + cpgs_file); if (get_has_counts_header(cpgs_file)) skip_counts_header(in); @@ -338,10 +328,12 @@ load_cpgs(const string &cpgs_file, vector &cpgs, MSite prev_site, the_site; while (read_site(in, the_site)) { if (!the_site.is_cpg() || distance(prev_site, the_site) < 2) - throw runtime_error("error: input is not symmetric-CpGs: " + cpgs_file); + throw std::runtime_error("error: input is not symmetric-CpGs: " + + cpgs_file); cpgs.push_back(the_site); reads.push_back(the_site.n_reads); - meth.push_back(make_pair(get_n_meth(the_site), get_n_unmeth(the_site))); + meth.push_back( + std::make_pair(get_n_meth(the_site), get_n_unmeth(the_site))); prev_site = the_site; } } @@ -349,7 +341,7 @@ load_cpgs(const string &cpgs_file, vector &cpgs, template static auto get_mean(InputIterator first, InputIterator last) -> T { - return accumulate(first, last, static_cast(0)) / + return std::accumulate(first, last, static_cast(0)) / std::distance(first, last); } @@ -361,15 +353,15 @@ check_sorted_within_chroms(T first, const T last) { if (first == last || first + 1 == last) return; - unordered_set chroms_seen; - string cur_chrom = ""; + std::unordered_set chroms_seen; + std::string cur_chrom = ""; for (auto fast = first + 1; fast != last; ++fast, ++first) { if (fast->chrom != cur_chrom) { - if (chroms_seen.find(fast->chrom) != end(chroms_seen)) { - throw runtime_error("input not grouped by chromosomes. " - "Error in the following line:\n" + - fast->tostring()); + if (chroms_seen.find(fast->chrom) != std::end(chroms_seen)) { + throw std::runtime_error("input not grouped by chromosomes. " + "Error in the following line:\n" + + fast->tostring()); } cur_chrom = fast->chrom; chroms_seen.insert(cur_chrom); @@ -377,9 +369,9 @@ check_sorted_within_chroms(T first, const T last) { } else { // first and fast are in the same chrom if (first->pos >= fast->pos) { - throw runtime_error("input file not sorted properly. " - "Error in the following lines:\n" + - first->tostring() + "\n" + fast->tostring()); + throw std::runtime_error("input file not sorted properly. " + "Error in the following lines:\n" + + first->tostring() + "\n" + fast->tostring()); } } } @@ -392,13 +384,13 @@ main_hmr(int argc, char *argv[]) { constexpr double min_coverage = 1.0; - string outfile; - string hypo_post_outfile; - string meth_post_outfile; + std::string outfile; + std::string hypo_post_outfile; + std::string meth_post_outfile; - size_t desert_size = 1000; - size_t max_iterations = 10; - size_t rng_seed = 408; + std::size_t desert_size = 1000; + std::size_t max_iterations = 10; + std::size_t rng_seed = 408; // run mode flags bool verbose = false; @@ -408,12 +400,12 @@ main_hmr(int argc, char *argv[]) { // corrections for small values const double tolerance = 1e-10; - string summary_file; + std::string summary_file; - string params_in_file; - string params_out_file; + std::string params_in_file; + std::string params_out_file; - const string description = + const std::string description = "Identify HMRs in methylomes. Methylation must be provided in the \ methcounts format (chrom, position, strand, context, \ methylation, reads). See the methcounts documentation for \ @@ -455,71 +447,71 @@ main_hmr(int argc, char *argv[]) { "input has extra fields (used for nanopore)", false, allow_extra_fields); opt_parse.set_show_defaults(); - vector leftover_args; + std::vector leftover_args; opt_parse.parse(argc, argv, leftover_args); if (argc == 1 || opt_parse.help_requested()) { - cerr << opt_parse.help_message() << endl - << opt_parse.about_message() << endl; + std::cerr << opt_parse.help_message() << '\n' + << opt_parse.about_message() << '\n'; return EXIT_SUCCESS; } if (opt_parse.about_requested()) { - cerr << opt_parse.about_message() << endl; + std::cerr << opt_parse.about_message() << '\n'; return EXIT_SUCCESS; } if (opt_parse.option_missing()) { - cerr << opt_parse.option_missing_message() << endl; + std::cerr << opt_parse.option_missing_message() << '\n'; return EXIT_SUCCESS; } if (leftover_args.empty()) { - cerr << opt_parse.help_message() << endl; + std::cerr << opt_parse.help_message() << '\n'; return EXIT_SUCCESS; } - const string cpgs_file = leftover_args.front(); + const std::string cpgs_file = leftover_args.front(); /****************** END COMMAND LINE OPTIONS *****************/ MSite::no_extra_fields = (allow_extra_fields == false); if (!is_msite_file(cpgs_file)) - throw runtime_error("malformed counts file: " + cpgs_file); + throw std::runtime_error("malformed counts file: " + cpgs_file); // separate the regions by chrom and by desert - vector cpgs; - vector> meth; - vector reads; + std::vector cpgs; + std::vector> meth; + std::vector reads; if (verbose) - cerr << "[reading methylation levels]" << endl; + std::cerr << "[reading methylation levels]\n"; load_cpgs(cpgs_file, cpgs, meth, reads); if (verbose) - cerr << "[checking if input is properly formatted]" << endl; - check_sorted_within_chroms(begin(cpgs), end(cpgs)); + std::cerr << "[checking if input is properly formatted]\n"; + check_sorted_within_chroms(std::begin(cpgs), std::end(cpgs)); if (PARTIAL_METH) make_partial_meth(reads, meth); - const auto mean_coverage = get_mean(cbegin(reads), cend(reads)); + const auto mean_coverage = get_mean(std::cbegin(reads), std::cend(reads)); if (verbose) - cerr << "[total_cpgs=" << size(cpgs) << "]" << '\n' - << "[mean_coverage=" << mean_coverage << "]" << endl; + std::cerr << "[total_cpgs=" << size(cpgs) << "]" << '\n' + << "[mean_coverage=" << mean_coverage << "]" << '\n'; // check for sufficient data if (mean_coverage < static_cast(min_coverage)) { if (verbose) - cerr << "error: insufficient data" - << " mean_coverage=" << mean_coverage - << " min_coverage=" << min_coverage << endl; + std::cerr << "error: insufficient data" + << " mean_coverage=" << mean_coverage + << " min_coverage=" << min_coverage << '\n'; if (!summary_file.empty()) { std::ofstream summary_out(summary_file); if (!summary_out) - throw runtime_error("failed to open: " + summary_file); - summary_out << hmr_summary({}).tostring() << endl; + throw std::runtime_error("failed to open: " + summary_file); + summary_out << hmr_summary({}).tostring() << '\n'; } return EXIT_SUCCESS; } // separate the regions by chrom and by desert, and eliminate // those isolated CpGs - vector reset_points; + std::vector reset_points; separate_regions(verbose, desert_size, cpgs, meth, reads, reset_points); const TwoStateHMM hmm(tolerance, max_iterations, verbose); @@ -537,7 +529,7 @@ main_hmr(int argc, char *argv[]) { max_iterations = 0; } else { - const double n_reads = get_mean(begin(reads), end(reads)); + const double n_reads = get_mean(std::cbegin(reads), std::cend(reads)); fg_alpha = 0.33 * n_reads; fg_beta = 0.67 * n_reads; bg_alpha = 0.67 * n_reads; @@ -548,22 +540,22 @@ main_hmr(int argc, char *argv[]) { hmm.BaumWelchTraining(meth, reset_points, p_fb, p_bf, fg_alpha, fg_beta, bg_alpha, bg_beta); // DECODE THE DOMAINS - vector state_ids; - vector posteriors; + std::vector state_ids; + std::vector posteriors; hmm.PosteriorDecoding(meth, reset_points, p_fb, p_bf, fg_alpha, fg_beta, bg_alpha, bg_beta, state_ids, posteriors); - vector domain_scores; + std::vector domain_scores; get_domain_scores(state_ids, meth, reset_points, domain_scores); - vector random_scores; + std::vector random_scores; shuffle_cpgs(rng_seed, hmm, meth, reset_points, p_fb, p_bf, fg_alpha, fg_beta, bg_alpha, bg_beta, random_scores); - vector p_values; + std::vector p_values; assign_p_values(random_scores, domain_scores, p_values); - if (domain_score_cutoff == numeric_limits::max() && + if (domain_score_cutoff == std::numeric_limits::max() && !domain_scores.empty()) domain_score_cutoff = get_stepup_cutoff(p_values, 0.01); @@ -572,7 +564,7 @@ main_hmr(int argc, char *argv[]) { write_params_file(params_out_file, fg_alpha, fg_beta, bg_alpha, bg_beta, p_fb, p_bf, domain_score_cutoff); - vector domains; + std::vector domains; if (!domain_scores.empty()) build_domains(cpgs, reset_points, state_ids, domains); @@ -581,11 +573,12 @@ main_hmr(int argc, char *argv[]) { of.open(outfile); std::ostream out(outfile.empty() ? std::cout.rdbuf() : of.rdbuf()); - size_t good_hmr_count = 0; + std::size_t good_hmr_count = 0; for (auto i = 0u; i < size(domains); ++i) if (p_values[i] < domain_score_cutoff) { domains[good_hmr_count] = domains[i]; - domains[good_hmr_count].set_name("HYPO" + to_string(good_hmr_count)); + domains[good_hmr_count].set_name("HYPO" + + std::to_string(good_hmr_count)); ++good_hmr_count; } domains.resize(good_hmr_count); @@ -595,9 +588,9 @@ main_hmr(int argc, char *argv[]) { if (!hypo_post_outfile.empty()) { if (verbose) - cerr << "[writing=" << hypo_post_outfile << "]" << endl; + std::cerr << "[writing=" << hypo_post_outfile << "]\n"; std::ofstream out(hypo_post_outfile); - for (size_t i = 0; i < cpgs.size(); ++i) { + for (std::size_t i = 0; i < cpgs.size(); ++i) { GenomicRegion cpg(as_gen_rgn(cpgs[i])); cpg.set_name(format_cpg_meth_tag(meth[i])); cpg.set_score(posteriors[i]); @@ -608,8 +601,8 @@ main_hmr(int argc, char *argv[]) { if (!meth_post_outfile.empty()) { std::ofstream out(meth_post_outfile); if (verbose) - cerr << "[writing=" << meth_post_outfile << "]" << endl; - for (size_t i = 0; i < cpgs.size(); ++i) { + std::cerr << "[writing=" << meth_post_outfile << "]\n"; + for (std::size_t i = 0; i < cpgs.size(); ++i) { GenomicRegion cpg(as_gen_rgn(cpgs[i])); cpg.set_name(format_cpg_meth_tag(meth[i])); cpg.set_score(1.0 - posteriors[i]); @@ -619,12 +612,12 @@ main_hmr(int argc, char *argv[]) { if (!summary_file.empty()) { std::ofstream summary_out(summary_file); if (!summary_out) - throw runtime_error("failed to open: " + summary_file); - summary_out << hmr_summary(domains).tostring() << endl; + throw std::runtime_error("failed to open: " + summary_file); + summary_out << hmr_summary(domains).tostring() << '\n'; } } - catch (runtime_error &e) { - cerr << e.what() << endl; + catch (std::exception &e) { + std::cerr << e.what() << '\n'; return EXIT_FAILURE; } return EXIT_SUCCESS; From b9e626249e9a4755be81203249b448f9486aa48e Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Tue, 7 Oct 2025 14:13:21 -0700 Subject: [PATCH 07/11] src/utils/format-reads.cpp: now requiring that the input be given as strict SAM format --- src/utils/format-reads.cpp | 57 ++++++++++++++++---------------------- 1 file changed, 24 insertions(+), 33 deletions(-) diff --git a/src/utils/format-reads.cpp b/src/utils/format-reads.cpp index 31b9c612..2d7a75a3 100644 --- a/src/utils/format-reads.cpp +++ b/src/utils/format-reads.cpp @@ -27,6 +27,15 @@ * A-rich and then changing that format to indicate T-rich. */ +// from dnmtools +#include "bam_record_utils.hpp" +#include "dnmt_error.hpp" + +// from smithlab +#include "OptionParser.hpp" +#include "smithlab_os.hpp" +#include "smithlab_utils.hpp" + #include #include @@ -40,21 +49,10 @@ #include #include -// from smithlab -#include "OptionParser.hpp" -#include "smithlab_os.hpp" -#include "smithlab_utils.hpp" - -// from dnmtools -#include "bam_record_utils.hpp" -#include "dnmt_error.hpp" - -using bamxx::bam_rec; - static std::int32_t -merge_mates(bam_rec &one, bam_rec &two, bam_rec &merged) { +merge_ends(bamxx::bam_rec &one, bamxx::bam_rec &two, bamxx::bam_rec &merged) { if (!are_mates(one, two)) - return -std::numeric_limits::max(); + return -std::numeric_limits::min(); // arithmetic easier using base 0 so subtracting 1 from pos const int one_s = get_pos(one); @@ -159,7 +157,7 @@ load_read_names(const std::string &inputfile, const std::size_t n_reads) { if (!hdr) throw dnmt_error("failed to get header: " + inputfile); - bam_rec aln; + bamxx::bam_rec aln; std::vector names; std::size_t count = 0; while (hts.read(hdr, aln) && count++ < n_reads) @@ -303,7 +301,7 @@ check_format_in_header(const std::string &input_format, } static inline bool -same_name(const bamxx::bam_rec &a, const bam_rec &b, +same_name(const bamxx::bam_rec &a, const bamxx::bam_rec &b, const std::size_t suff_len) { // "+ 1" below: extranul counts *extras*; we don't want *any* nulls const std::uint16_t a_l = a.b->core.l_qname - (a.b->core.l_extranul + 1); @@ -315,16 +313,15 @@ same_name(const bamxx::bam_rec &a, const bam_rec &b, } static inline void -swap(bam_rec &a, bam_rec &b) { +swap(bamxx::bam_rec &a, bamxx::bam_rec &b) { std::swap(a.b, b.b); } static void -format(const bool sam_orientation, const std::string &cmd, - const std::size_t n_threads, const std::string &inputfile, - const std::string &outfile, const bool bam_format, - const std::string &input_format, const std::size_t suff_len, - const std::int32_t max_frag_len) { +format(const std::string &cmd, const std::size_t n_threads, + const std::string &inputfile, const std::string &outfile, + const bool bam_format, const std::string &input_format, + const std::size_t suff_len, const std::int32_t max_frag_len) { static const dnmt_error bam_write_err{"error writing bam"}; bamxx::bam_tpool tp(n_threads); @@ -348,7 +345,7 @@ format(const bool sam_orientation, const std::string &cmd, tp.set_io(out); } - bam_rec aln, prev_aln, merged; + bamxx::bam_rec aln, prev_aln, merged; bool previous_was_merged = false; const bool empty_reads_file = !hts.read(hdr, aln); @@ -359,7 +356,7 @@ format(const bool sam_orientation, const std::string &cmd, // in the flags has been reverse complemented and is not the original // sequence from the FASTQ file. So we need to undo the reverse // complementing. - if (sam_orientation && bam_is_rev(aln)) + if (bam_is_rev(aln)) revcomp_qseq(aln); standardize_format(input_format, aln); @@ -371,7 +368,7 @@ format(const bool sam_orientation, const std::string &cmd, if (get_tid(aln) == -1) continue; - if (sam_orientation && bam_is_rev(aln)) + if (bam_is_rev(aln)) revcomp_qseq(aln); standardize_format(input_format, aln); @@ -379,7 +376,7 @@ format(const bool sam_orientation, const std::string &cmd, // below: essentially check for dovetail if (!bam_is_rev(aln)) swap(prev_aln, aln); - const auto frag_len = merge_mates(prev_aln, aln, merged); + const auto frag_len = merge_ends(prev_aln, aln, merged); if (frag_len > 0 && frag_len < max_frag_len) { if (is_a_rich(merged)) flip_conversion(merged); @@ -446,7 +443,6 @@ main_format(int argc, char *argv[]) { std::size_t suff_len = 0; bool single_end{false}; bool verbose{false}; - bool sam_orientation{false}; bool force{false}; std::size_t n_threads{1}; @@ -471,9 +467,6 @@ main_format(int argc, char *argv[]) { opt_parse.add_opt("check", 'c', "number of reads to confirm read name suffix", false, n_reads_to_check); - opt_parse.add_opt("sam", 'S', - "negative strand mappers are reverse complemented", false, - sam_orientation); opt_parse.add_opt("force", 'F', "force formatting for " "mixed single and paired reads", @@ -526,8 +519,6 @@ main_format(int argc, char *argv[]) { << "[output type: " << (bam_format ? "B" : "S") << "AM]\n" << "[force formatting: " << (force ? "yes" : "no") << "]\n" << "[threads requested: " << n_threads << "]\n" - << "[SAM orientation: " << std::boolalpha << sam_orientation - << "]\n" << "[command line: \"" << command << "\"]\n"; check_input_file(infile); @@ -559,8 +550,8 @@ main_format(int argc, char *argv[]) { if (verbose && !single_end) std::cerr << "[readname suffix length: " << suff_len << "]" << '\n'; - format(sam_orientation, command, n_threads, infile, outfile, bam_format, - input_format, suff_len, max_frag_len); + format(command, n_threads, infile, outfile, bam_format, input_format, + suff_len, max_frag_len); } catch (const std::exception &e) { std::cerr << e.what() << '\n'; From 5fba078f6a9c8b6962aa9904f22c46ceb6c2ea20 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Tue, 7 Oct 2025 14:18:31 -0700 Subject: [PATCH 08/11] src/utils/format-reads.cpp: fixing a bug of negative of negative --- src/utils/format-reads.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils/format-reads.cpp b/src/utils/format-reads.cpp index 2d7a75a3..f2b674d9 100644 --- a/src/utils/format-reads.cpp +++ b/src/utils/format-reads.cpp @@ -52,7 +52,7 @@ static std::int32_t merge_ends(bamxx::bam_rec &one, bamxx::bam_rec &two, bamxx::bam_rec &merged) { if (!are_mates(one, two)) - return -std::numeric_limits::min(); + return std::numeric_limits::min(); // arithmetic easier using base 0 so subtracting 1 from pos const int one_s = get_pos(one); From c2ae4878aed042b555388ca7affa9ada0cd88ca4 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Tue, 7 Oct 2025 14:18:54 -0700 Subject: [PATCH 09/11] src/common/LevelsCounter.cpp: cleanup --- src/common/LevelsCounter.cpp | 24 ++++++++++-------------- 1 file changed, 10 insertions(+), 14 deletions(-) diff --git a/src/common/LevelsCounter.cpp b/src/common/LevelsCounter.cpp index de9d7109..7931a613 100644 --- a/src/common/LevelsCounter.cpp +++ b/src/common/LevelsCounter.cpp @@ -1,16 +1,16 @@ -/* Copyright (C) 2018-2022 Andrew D. Smith +/* Copyright (C) 2018-2025 Andrew D. Smith * * Authors: Andrew D. Smith * - * This is free software: you can redistribute it and/or modify it - * under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. + * This is free software: you can redistribute it and/or modify it under the + * terms of the GNU General Public License as published by the Free Software + * Foundation, either version 3 of the License, or (at your option) any later + * version. * * This software is distributed in the hope that it will be useful, but - * WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * General Public License for more details. + * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY + * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. */ #include "LevelsCounter.hpp" @@ -20,10 +20,6 @@ #include #include -using std::runtime_error; -using std::string; -using std::to_string; - void LevelsCounter::update(const MSite &s) { if (s.is_mutated()) { @@ -135,8 +131,8 @@ operator<<(std::ostream &out, const LevelsCounter &cs) { static void check_label(const std::string &observed, const std::string expected) { if (observed != expected) - throw runtime_error("bad levels format [" + observed + "," + expected + - "]"); + throw std::runtime_error("bad levels format [" + observed + "," + expected + + "]"); } std::istream & From 028716d95fbb869cf23546a335fa01db762ae19f Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Tue, 7 Oct 2025 14:19:35 -0700 Subject: [PATCH 10/11] src/analysis/roimethstat.cpp: allowing an option for extra cols in output, which is useful for current counts format output by counts-nano --- src/analysis/roimethstat.cpp | 407 ++++++++++++++++++----------------- 1 file changed, 207 insertions(+), 200 deletions(-) diff --git a/src/analysis/roimethstat.cpp b/src/analysis/roimethstat.cpp index 1f7fdb6e..4555e816 100644 --- a/src/analysis/roimethstat.cpp +++ b/src/analysis/roimethstat.cpp @@ -1,34 +1,33 @@ /* roimethstat: average methylation in each of a set of regions * - * Copyright (C) 2014-2024 Andrew D. Smith + * Copyright (C) 2014-2025 Andrew D. Smith * * Authors: Andrew D. Smith and Masaru Nakajima * - * This program is free software: you can redistribute it and/or - * modify it under the terms of the GNU General Public License as - * published by the Free Software Foundation, either version 3 of the - * License, or (at your option) any later version. + * This program is free software: you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) + * any later version. * - * This program is distributed in the hope that it will be useful, but - * WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * General Public License for more details. + * This program is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for + * more details. */ #include #include +#include +#include #include #include #include #include #include #include -#include #include #include -#include -#include #include "GenomicRegion.hpp" #include "LevelsCounter.hpp" @@ -38,31 +37,7 @@ #include "smithlab_utils.hpp" #include "xcounts_utils.hpp" -using std::cerr; -using std::cout; -using std::distance; -using std::endl; -using std::ifstream; -using std::is_sorted; -using std::pair; -using std::runtime_error; -using std::string; -using std::to_string; -using std::unordered_map; -using std::unordered_set; -using std::vector; -using std::from_chars; -using std::size; -using std::cend; -using std::ostream; -using std::size; - -using bamxx::bgzf_file; - -namespace fs = std::filesystem; - - -static string +[[nodiscard]] static std::string format_levels_counter(const LevelsCounter &lc) { // ... // (7) weighted mean methylation @@ -80,22 +55,20 @@ format_levels_counter(const LevelsCounter &lc) { << lc.total_sites << '\t' << lc.sites_covered << '\t' << lc.total_c << '\t' - << (lc.total_c + lc.total_t); + << lc.coverage(); // clang-format on return oss.str(); } - struct genomic_interval { - string chrom{}; - uint64_t start_pos{}; - uint64_t end_pos{}; + std::string chrom{}; + std::uint64_t start_pos{}; + std::uint64_t end_pos{}; }; - static void update(LevelsCounter &lc, const xcounts_entry &xse) { - const uint64_t n_reads = xse.n_meth + xse.n_unmeth; + const std::uint64_t n_reads = xse.n_meth + xse.n_unmeth; if (n_reads > 0) { ++lc.sites_covered; lc.max_depth = std::max(lc.max_depth, n_reads); @@ -113,62 +86,61 @@ update(LevelsCounter &lc, const xcounts_entry &xse) { static void process_chrom(const bool report_more_info, const char level_code, - const vector &intervals, - const vector &sites, - ostream &out) { + const std::vector &intervals, + const std::vector &sites, std::ostream &out) { - uint64_t j = 0; - for (auto i = 0ul; i < intervals.size(); ++i) { - while (j < size(sites) && sites[j].pos < intervals[i].get_start()) ++j; + std::uint64_t j = 0; + for (auto i = 0ul; i < std::size(intervals); ++i) { + while (j < std::size(sites) && sites[j].pos < intervals[i].get_start()) + ++j; LevelsCounter lc; - while (j < size(sites) && sites[j].pos < intervals[i].get_end()) + while (j < std::size(sites) && sites[j].pos < intervals[i].get_end()) update(lc, sites[j++]); GenomicRegion r(intervals[i]); r.set_score(level_code == 'w' ? lc.mean_meth_weighted() - : (level_code == 'u' ? lc.mean_meth() - : lc.fractional_meth())); + : (level_code == 'u' ? lc.mean_meth() + : lc.fractional_meth())); r.set_name(r.get_name() + "_" + std::to_string((level_code == 'w' - ? lc.coverage() - : (level_code == 'u' ? lc.sites_covered - : lc.total_called())))); + ? lc.coverage() + : (level_code == 'u' ? lc.sites_covered + : lc.total_called())))); out << r; - if (report_more_info) out << '\t' << format_levels_counter(lc); + if (report_more_info) + out << '\t' << format_levels_counter(lc); out << '\n'; } } static void process_chrom(const bool report_more_info, - const vector &intervals, - ostream &out) { + const std::vector &intervals, std::ostream &out) { LevelsCounter lc; - const string lc_formatted = format_levels_counter(lc); - for (const auto &r: intervals) { + const std::string lc_formatted = format_levels_counter(lc); + for (const auto &r : intervals) { out << r; - if (report_more_info) out << '\t' << lc_formatted; + if (report_more_info) + out << '\t' << lc_formatted; out << '\n'; } } - static void -process_from_xcounts(const uint32_t n_threads, - const bool report_more_info, - const char level_code, const string &xsym_file, - const vector &intervals, - ostream &out) { +process_from_xcounts(const std::uint32_t n_threads, const bool report_more_info, + const char level_code, const std::string &xsym_file, + const std::vector &intervals, + std::ostream &out) { const auto sites_by_chrom = read_xcounts_by_chrom(n_threads, xsym_file); // const auto intervals = get_GenomicRegions(intervals_file); - vector> intervals_by_chrom; - string prev_chrom; - for (auto i = 0u; i < size(intervals); ++i) { + std::vector> intervals_by_chrom; + std::string prev_chrom; + for (auto i = 0u; i < std::size(intervals); ++i) { if (intervals[i].get_chrom() != prev_chrom) { - intervals_by_chrom.push_back(vector()); + intervals_by_chrom.push_back(std::vector()); prev_chrom = intervals[i].get_chrom(); } intervals_by_chrom.back().push_back(intervals[i]); @@ -177,76 +149,85 @@ process_from_xcounts(const uint32_t n_threads, for (const auto &intervals : intervals_by_chrom) { const auto chrom_name = intervals.front().get_chrom(); const auto sites = sites_by_chrom.find(chrom_name); - if (sites != cend(sites_by_chrom)) - process_chrom(report_more_info, level_code, - intervals, sites->second, out); + if (sites != std::cend(sites_by_chrom)) + process_chrom(report_more_info, level_code, intervals, sites->second, + out); else process_chrom(report_more_info, intervals, out); } } - - -bool +[[nodiscard]] bool cmp_within_chrom(const GenomicRegion &r1, const GenomicRegion &r2) { - return (r1.get_start() < r2.get_start() || - (r1.get_start() == r2.get_start() && - (r1.get_end() < r2.get_end() || - (r1.get_end() == r2.get_end() && - (r1.get_strand() < r2.get_strand()))))); + return r1.get_start() < r2.get_start() || + (r1.get_start() == r2.get_start() && + (r1.get_end() < r2.get_end() || + (r1.get_end() == r2.get_end() && + (r1.get_strand() < r2.get_strand())))); } -bool +[[nodiscard]] bool cmp_within_chrom(const MSite &s1, const MSite &s2) { return s1.pos < s2.pos; } -template typename T::const_iterator +template +[[nodiscard]] typename T::const_iterator get_chrom_order(const T &order, const MSite &s) { return order.find(s.chrom); } -template typename T::const_iterator +template +[[nodiscard]] typename T::const_iterator get_chrom_order(const T &order, const GenomicRegion &r) { return order.find(r.get_chrom()); } struct cmp_chrom_order { - cmp_chrom_order(const unordered_map &m): order{m} {} + cmp_chrom_order(const std::unordered_map &m) : + order{m} {} - template bool operator()(const T &a, const T &b) const { + template + [[nodiscard]] bool + operator()(const T &a, const T &b) const { const auto ac = get_chrom_order(order, a); - if (ac == cend(order)) return false; + if (ac == std::cend(order)) + return false; const auto bc = get_chrom_order(order, b); - if (bc == cend(order)) return true; + if (bc == std::cend(order)) + return true; const auto c = static_cast(ac->second) - static_cast(bc->second); return c < 0 || (c == 0 && cmp_within_chrom(a, b)); } - const unordered_map ℴ + const std::unordered_map ℴ }; -static vector -get_chroms(const string &filename) { +[[nodiscard]] static std::vector +get_chroms(const std::string &filename) { // this function might look slow but it seems not to bottleneck - constexpr auto not_in_chrom = [](const uint8_t c) { + constexpr auto not_in_chrom = [](const std::uint8_t c) { return c == ' ' || c == '\t' || c == '\n'; }; - bgzf_file in(filename, "r"); - if (!in) throw runtime_error("cannot open file: " + filename); + bamxx::bgzf_file in(filename, "r"); + if (!in) + throw std::runtime_error("cannot open file: " + filename); - string prev_chrom; - string line; - vector chroms_seen; + std::string prev_chrom; + std::string line; + std::vector chroms_seen; while (getline(in, line)) { - const auto chrom_end = find_if(cbegin(line), cend(line), not_in_chrom); - if (chrom_end == cend(line)) - throw runtime_error("bad counts line: " + line); - if (!equal(cbegin(prev_chrom), cend(prev_chrom), cbegin(line), chrom_end)) { - const string chrom{cbegin(line), chrom_end}; - auto x = find(cbegin(chroms_seen), cend(chroms_seen), chrom); - if (x != cend(chroms_seen)) throw runtime_error("chroms not sorted"); + const auto chrom_end = + std::find_if(std::cbegin(line), std::cend(line), not_in_chrom); + if (chrom_end == std::cend(line)) + throw std::runtime_error("bad counts line: " + line); + if (!std::equal(std::cbegin(prev_chrom), std::cend(prev_chrom), + std::cbegin(line), chrom_end)) { + const std::string chrom{std::cbegin(line), chrom_end}; + auto x = find(std::cbegin(chroms_seen), std::cend(chroms_seen), chrom); + if (x != std::cend(chroms_seen)) + throw std::runtime_error("chroms not sorted"); chroms_seen.push_back(chrom); prev_chrom = std::move(chrom); } @@ -255,17 +236,17 @@ get_chroms(const string &filename) { } template -static inline pair -region_bounds(const unordered_map &chrom_order, +[[nodiscard]] static inline std::pair +region_bounds(const std::unordered_map &chrom_order, ForwardIt first, ForwardIt last, const T ®ion) { // ADS: in the sites below, a and b, the strand is set to '+' // always. Otherwise the ordering on MSite would force strand to be - // used. The consequence is that sites on the positive strand would - // not be included if they begin at the same nucleotide where a - // region starts if that region is on the negative strand - // const char strand(region.get_strand()); + // used. The consequence is that sites on the positive strand would not be + // included if they begin at the same nucleotide where a region starts if + // that region is on the negative strand const char + // strand(region.get_strand()); - const string chrom(region.get_chrom()); + const std::string chrom(region.get_chrom()); const MSite a(chrom, region.get_start(), '+', "", 0, 0); const MSite b(chrom, region.get_end(), '+', "", 0, 0); @@ -277,53 +258,59 @@ region_bounds(const unordered_map &chrom_order, return {lower_bound(first, last, a, cmp), lower_bound(first, last, b, cmp)}; } - -static bool -is_sorted_within_chrom(const vector &sites) { +[[nodiscard]] static bool +is_sorted_within_chrom(const std::vector &sites) { constexpr auto pos_cmp = [](const MSite &s1, const MSite &s2) { return s1.pos < s2.pos; }; - auto a = cbegin(sites); - while (a != cend(sites)) { + auto a = std::cbegin(sites); + while (a != std::cend(sites)) { const auto a_chrom = a->chrom; - const auto b = find_if(a, cend(sites), [&a_chrom](const MSite &s) { + const auto b = find_if(a, std::cend(sites), [&a_chrom](const MSite &s) { return s.chrom != a_chrom; }); - if (!is_sorted(a, b, pos_cmp)) return false; + if (!std::is_sorted(a, b, pos_cmp)) + return false; a = b; } return true; } -static vector -read_sites(const string &filename) { - vector sites; - bgzf_file in(filename, "r"); +[[nodiscard]] static std::vector +read_sites(const std::string &filename) { + std::vector sites; + bamxx::bgzf_file in(filename, "r"); if (in) { MSite s; - while (read_site(in, s)) sites.push_back(s); + while (read_site(in, s)) + sites.push_back(s); } return sites; } static void -process_preloaded(const bool VERBOSE, const bool report_more_info, - const char level_code, const string &sites_file, - const unordered_map &chrom_order, - const vector ®ions, ostream &out) { +process_preloaded( + const bool VERBOSE, const bool report_more_info, const char level_code, + const std::string &sites_file, + const std::unordered_map &chrom_order, + const std::vector ®ions, std::ostream &out) { const auto sites = read_sites(sites_file); - if (sites.empty()) throw runtime_error("failed to read sites: " + sites_file); + if (sites.empty()) + throw std::runtime_error("failed to read sites: " + sites_file); if (!is_sorted_within_chrom(sites)) - throw runtime_error("sites not sorted: " + sites_file); + throw std::runtime_error("sites not sorted: " + sites_file); - if (VERBOSE) cerr << "[n_sites=" << sites.size() << "]" << endl; + if (VERBOSE) + std::cerr << "[n_sites=" << std::size(sites) << "]\n"; for (auto &r : regions) { LevelsCounter lc; - const auto b = region_bounds(chrom_order, cbegin(sites), cend(sites), r); - for (auto c = b.first; c != b.second; ++c) lc.update(*c); + const auto b = + region_bounds(chrom_order, std::cbegin(sites), std::cend(sites), r); + for (auto c = b.first; c != b.second; ++c) + lc.update(*c); const double score = level_code == 'w' @@ -338,30 +325,33 @@ process_preloaded(const bool VERBOSE, const bool report_more_info, } } -template static inline bool -not_after(const T &chrom_order, const MSite &site, const uint32_t &chrom_idx, - const size_t end_pos) { +template +[[nodiscard]] static inline bool +not_after(const T &chrom_order, const MSite &site, + const std::uint32_t &chrom_idx, const std::size_t end_pos) { const auto sc_itr = chrom_order.find(site.chrom); - if (sc_itr == cend(chrom_order)) - throw runtime_error("lookup failure: " + site.chrom); + if (sc_itr == std::cend(chrom_order)) + throw std::runtime_error("lookup failure: " + site.chrom); const auto sc = sc_itr->second; return (sc == chrom_idx && site.pos < end_pos) || sc < chrom_idx; } -static inline bool -at_or_after(const MSite &site, const string &chrom, const size_t start_pos) { +[[nodiscard]] static inline bool +at_or_after(const MSite &site, const std::string &chrom, + const std::size_t start_pos) { // not using dictionary in this function because equality of chrom // can be tested using the strings regardless of index return (start_pos <= site.pos && site.chrom == chrom); } -static LevelsCounter -calc_site_stats(ifstream &sites_in, const GenomicRegion ®ion, - const unordered_map &chrom_order) { - const string chrom{region.get_chrom()}; +[[nodiscard]] static LevelsCounter +calc_site_stats( + std::ifstream &sites_in, const GenomicRegion ®ion, + const std::unordered_map &chrom_order) { + const std::string chrom{region.get_chrom()}; const auto chrom_itr = chrom_order.find(chrom); - if (chrom_itr == cend(chrom_order)) - throw runtime_error("lookup failure: " + chrom); + if (chrom_itr == std::cend(chrom_order)) + throw std::runtime_error("lookup failure: " + chrom); const auto chrom_idx = chrom_itr->second; const auto start_pos = region.get_start(); const auto end_pos = region.get_end(); @@ -371,7 +361,8 @@ calc_site_stats(ifstream &sites_in, const GenomicRegion ®ion, MSite site; while (sites_in >> site && not_after(chrom_order, site, chrom_idx, end_pos)) - if (at_or_after(site, chrom, start_pos)) lc.update(site); + if (at_or_after(site, chrom, start_pos)) + lc.update(site); sites_in.clear(); // clear here otherwise an eof means we skip later // intervals, which could be valid @@ -379,12 +370,14 @@ calc_site_stats(ifstream &sites_in, const GenomicRegion ®ion, } static void -process_on_disk(const bool report_more_info, const char level_code, - const string &sites_file, - const unordered_map &chrom_order, - const vector ®ions, ostream &out) { - ifstream in(sites_file); - if (!in) throw runtime_error("failed to open file: " + sites_file); +process_on_disk( + const bool report_more_info, const char level_code, + const std::string &sites_file, + const std::unordered_map &chrom_order, + const std::vector ®ions, std::ostream &out) { + std::ifstream in(sites_file); + if (!in) + throw std::runtime_error("failed to open file: " + sites_file); for (auto ®ion : regions) { const auto lc = calc_site_stats(in, region, chrom_order); @@ -401,27 +394,28 @@ process_on_disk(const bool report_more_info, const char level_code, } } -static size_t -get_bed_columns(const string ®ions_file) { - ifstream in(regions_file); - if (!in) throw runtime_error("failed to open file: " + regions_file); +[[nodiscard]] static std::size_t +get_bed_columns(const std::string ®ions_file) { + std::ifstream in(regions_file); + if (!in) + throw std::runtime_error("failed to open file: " + regions_file); - string line; - getline(in, line); + std::string line; + std::getline(in, line); std::istringstream iss(line); - string token; - size_t n_columns = 0; - while (iss >> token) ++n_columns; + std::string token; + std::size_t n_columns = 0; + while (iss >> token) + ++n_columns; return n_columns; } - int main_roimethstat(int argc, char *argv[]) { try { - static const string description = R"""( + static const std::string description = R"""( Compute average site methylation levels in each interval from a given set of genomic intervals. The 5th column (the "score" column in BED format) is determined by the '-l' or '-level' argument. @@ -436,18 +430,19 @@ Columns (beyond the first 6) in the BED format output: (13) total number of observations from reads in the region )"""; - static const string default_name_prefix = "X"; + static const std::string default_name_prefix = "X"; bool VERBOSE = false; bool print_numeric_only = false; bool preload = false; bool report_more_info = false; bool sort_data_if_needed = false; - uint32_t n_threads = 1; + bool allow_extra_fields = false; + std::uint32_t n_threads = 1; - string level_code = "w"; + std::string level_code = "w"; - string outfile; + std::string outfile; /****************** COMMAND LINE OPTIONS ********************/ OptionParser opt_parse(strip_path(argv[0]), description, @@ -468,85 +463,97 @@ Columns (beyond the first 6) in the BED format output: false, report_more_info); opt_parse.add_opt("threads", 't', "threads to use (if input compressed)", false, n_threads); + opt_parse.add_opt("relaxed", '\0', + "input has extra fields (used for nanopore)", false, + allow_extra_fields); opt_parse.add_opt("verbose", 'v', "print more run info", false, VERBOSE); - vector leftover_args; + std::vector leftover_args; opt_parse.parse(argc, argv, leftover_args); if (argc == 1 || opt_parse.help_requested()) { - cerr << opt_parse.help_message() - << opt_parse.about_message_raw() << endl; + std::cerr << opt_parse.help_message() << opt_parse.about_message_raw() + << '\n'; return EXIT_SUCCESS; } if (opt_parse.about_requested()) { - cerr << opt_parse.about_message_raw() << endl; + std::cerr << opt_parse.about_message_raw() << '\n'; return EXIT_SUCCESS; } if (opt_parse.option_missing()) { - cerr << opt_parse.option_missing_message() << endl; + std::cerr << opt_parse.option_missing_message() << '\n'; return EXIT_FAILURE; } - if (leftover_args.size() != 2) { - cerr << opt_parse.help_message() << endl; + if (std::size(leftover_args) != 2) { + std::cerr << opt_parse.help_message() << '\n'; return EXIT_FAILURE; } if (level_code != "w" && level_code != "u" && level_code != "f") { - cerr << "selected level must be in {w, u, f}" << endl; + std::cerr << "selected level must be in {w, u, f}\n"; return EXIT_FAILURE; } - const string regions_file = leftover_args.front(); - const string sites_file = leftover_args.back(); + const std::string regions_file = leftover_args.front(); + const std::string sites_file = leftover_args.back(); /****************** END COMMAND LINE OPTIONS *****************/ + MSite::no_extra_fields = (allow_extra_fields == false); + const bool is_xcounts = get_is_xcounts_file(sites_file); if (!is_msite_file(sites_file) && !is_xcounts) - throw runtime_error("dnmtools counts or xcounts format required: " + - sites_file); + throw std::runtime_error("dnmtools counts or xcounts format required: " + + sites_file); // make a map that specifies their order; otherwise we can't // ensure regions are sorted in the same way - unordered_map chrom_order; + std::unordered_map chrom_order; if (!is_xcounts) for (auto &i : get_chroms(sites_file)) - chrom_order.emplace(i, chrom_order.size()); + chrom_order.emplace(i, std::size(chrom_order)); - if (VERBOSE) cerr << "loading regions" << endl; + if (VERBOSE) + std::cerr << "loading regions" << '\n'; - if (!fs::is_regular_file(regions_file)) + if (!std::filesystem::is_regular_file(regions_file)) // otherwise we could not read the file twice - throw runtime_error("regions file must be regular file"); + throw std::runtime_error("regions file must be regular file"); // MAGIC: below allow for ==3 or >=6 columns in the bed format const auto n_columns = get_bed_columns(regions_file); if (n_columns != 3 && n_columns < 6) - throw runtime_error("format must be 3 or 6+ column bed: " + regions_file); + throw std::runtime_error("format must be 3 or 6+ column bed: " + + regions_file); if (is_msite_file(regions_file)) - throw runtime_error("seems to be a counts file: " + regions_file + "\n" + - "check order of input arguments"); + throw std::runtime_error("seems to be a counts file: " + regions_file + + "\ncheck order of input arguments"); - vector regions; + std::vector regions; ReadBEDFile(regions_file, regions); - if (!is_sorted(begin(regions), end(regions))) { + if (!std::is_sorted(std::begin(regions), std::end(regions))) { if (sort_data_if_needed) { - if (VERBOSE) cerr << "sorting regions" << endl; - sort(begin(regions), end(regions), cmp_chrom_order(chrom_order)); + if (VERBOSE) + std::cerr << "sorting regions\n"; + std::sort(std::begin(regions), std::end(regions), + cmp_chrom_order(chrom_order)); } else - throw runtime_error("not sorted: " + regions_file + " (consider -s)"); + throw std::runtime_error("not sorted: " + regions_file + + " (consider -s)"); } // if columns are 3 we name the regions otherwise that column will // be empty and the output format will be wrong if (n_columns == 3) - for (auto i = 0u; i < regions.size(); ++i) - regions[i].set_name(default_name_prefix + to_string(i)); + for (auto i = 0u; i < std::size(regions); ++i) + regions[i].set_name(default_name_prefix + std::to_string(i)); - if (VERBOSE) cerr << "[n_regions=" << regions.size() << "]" << endl; + if (VERBOSE) + std::cerr << "[n_regions=" << std::size(regions) << "]\n"; std::ofstream of; if (!outfile.empty()) { of.open(outfile); - if (!of) throw runtime_error("failed to open outfile: " + outfile); + if (!of) + throw std::runtime_error("failed to open outfile: " + outfile); } - ostream out(outfile.empty() ? cout.rdbuf() : of.rdbuf()); + std::ostream out(outfile.empty() ? std::cout.rdbuf() : of.rdbuf()); if (is_xcounts) process_from_xcounts(n_threads, report_more_info, level_code[0], @@ -559,7 +566,7 @@ Columns (beyond the first 6) in the BED format output: regions, out); } catch (const std::exception &e) { - cerr << e.what() << endl; + std::cerr << e.what() << '\n'; return EXIT_FAILURE; } return EXIT_SUCCESS; From 8dff2765a563ad99b62ac409fd930964c97ca689 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Tue, 7 Oct 2025 14:21:09 -0700 Subject: [PATCH 11/11] data/md5sum.txt: updates for changes to abismal submodule and very minor change to a tag in the output of levels. Might break. --- data/md5sum.txt | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/data/md5sum.txt b/data/md5sum.txt index 738ab53d..a66a8504 100644 --- a/data/md5sum.txt +++ b/data/md5sum.txt @@ -7,8 +7,6 @@ e73facd597c3b903cbfe29afa9f58371 tests/reads.counts 0f72560aa101e85679783a1ecaf80615 tests/reads.epiread 9dbd476424d48a8d0f043dfc00af0d23 tests/reads.fmt.srt.sam 4085cc74b003a918b4a4743fca7922a4 tests/reads.hmr -5c6fbbc6e565f23afa4f597f2483acb0 tests/reads.levels -a389e4ee1687895517b0b0d26e9338ee tests/reads.mstats d8856f9731af76b8a4ab3cc7d667cdb2 tests/reads.ustats bcbf01be810cbf4051292813eb6b9225 tests/tRex1.idx ec6a686617cad31e9f7a37a3d378e6ed tests/two_epialleles.states @@ -17,8 +15,10 @@ d947fe3d61ef7b1564558a69608f0e64 tests/methylome.pmd d41d8cd98f00b204e9800998ecf8427e tests/two_epialleles.amr 001b9d966f62fa439b24cf2198cc3de5 tests/reads.counts.sym 2b8a0406015458be51b8b1c9e58b3602 tests/tRex1_promoters.roi.bed -9bdb361091d1c0626df30be8ba2c408c tests/reads.sam 33640b24cb64ad3179f364af5a887f95 tests/reads.fmt.sam b5a63997c57dcde5c3a6635f7beb2cce tests/reads.fmt.srt.uniq.sam 3ac2b51545740bafd1a548ba7f73e739 tests/reads.xcounts 054fe804a32063c80862fbee30f74579 tests/reads.unxcounts +830157684f1ddbf1f1c37d354188dc2b tests/reads.sam +8dbcdabecb6cfe6aebb73c94c605f696 tests/reads.mstats +490723e9af084c8f957f5a265cf02994 tests/reads.levels