From 7d1738bf91389f08a1c64b6a1b72e43caf69c975 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 15 Jun 2025 17:12:35 -0700 Subject: [PATCH 1/2] docs/content/counts-nano.md: adding first docs for counts-nano --- docs/content/counts-nano.md | 104 ++++++++++++++++++++++++++++++++++++ 1 file changed, 104 insertions(+) create mode 100644 docs/content/counts-nano.md diff --git a/docs/content/counts-nano.md b/docs/content/counts-nano.md new file mode 100644 index 00000000..b071ba77 --- /dev/null +++ b/docs/content/counts-nano.md @@ -0,0 +1,104 @@ +# counts-nano - compute single-site methylation from nanopore data + +## Synopsis +```console +$ dnmtools counts-nano [OPTIONS] -c +``` + +## Description + +The `counts-nano` command introduced in v1.5.0 is designed specifically to +generate DNMTools [counts](../counts) format files from nanopore data called +for the `5mCG_5hmCG` modification. Currently this is only supported for +methylation and hydroxymethylation called at CpG sites. + +More documentation will come as this tool evolves, but for now: + +- Most behavior is very similar to what you will find from [counts](../counts). +- Mutation information is not estimated by `nano-counts`. +- Currently this only works for CpG sites and when the only modified sites are + marked as `C+m?` or `C+h?` in the `MM` field of each BAM/SAM read record. +- The first 6 columns of the output are the same as explained in the + [counts](../counts) format, except the fraction for the 5th column is both + 5mC and 5hmC. The 7th column is for 5hmC alone and the 8th is for 5mC alone. +- The methylation levels will not result in integer values when multiplied by + the number of reads because probabilities on modifications are used, so + methylation levels for each site are expected values (the best estimates we + can make), and do not use arbitrary cutoffs. +- Other commands in DNMTools have been modified to use this form of expected + methylation level, and behave as previously for bisulfite sequencing data, + but have updated behavior when the data is from nanopore. The user does not + need to specify the technology used. +- Some commands need to use a `-relaxed` flag to work with the additional + columns in the output from `counts-nano` compared with `counts`. For + commands without this option, simply do `cut -f1-6` on the output of + `counts-nano` to remove those. + +## Options + +```txt +-o, -output +``` +Output file name. The default is to write output to the terminal, +which is not useful unless commands are piped. + +```txt +-c, -chrom +``` +Reference genome file, which must be in FASTA format. This is +required. + +```txt +-t, -threads +``` + +The number of threads to use. This is only really helpful if the input is BAM +(not helpful for SAM), and the output is to be zipped (see `-z` below). These +threads are used to decompress BAM input and compress gzip output. If only one +of these conditions holds, using more threads can still help. Because most +computation in `counts-nano` is processing reads sequentially, using too many +threads will have decreasing returns. + +```txt +-z, -zip +``` + +The output should be zipped (in gzip format). This is not deduced by the +filename, but specifying this argument should be accompanied by using a `.gz` +filename suffix for the output. + +```txt +-n, -cpg-only +``` + +Print only CpG context cytosines. This significantly reduces the output size +in most genomes. Note that using this option does not merge data as symmetric +CpGs. + +```txt +-sym +``` + +This will turn on `-n, -cpg-only` automatically and will output symmetric CpG +sites, with each level including all counts and methylation levels as a +(weighted) average of both strands. + +```txt +-H, -header +``` + +Add a header to the output file to identify the reference genome. This will be +in the form of "comment" lines beginning with `#`. This is not required for most +downstream processing, but is used by commands that check for consistency with +a reference genome. + +```txt +-v, -verbose +``` + +Report more information while the program is running. + +```txt +-progress +``` +Show progress while the program is running. From 5d4554598384bdc5a525aad187da13cba3ec43db Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Mon, 16 Jun 2025 20:43:52 -0700 Subject: [PATCH 2/2] src/radmeth/methdiff.cpp: updates to how output is formatted for writing, which should gain some speed, but that was likely not a bottleneck. Also cleaned up code stype --- src/radmeth/methdiff.cpp | 251 +++++++++++++++++++-------------------- 1 file changed, 125 insertions(+), 126 deletions(-) diff --git a/src/radmeth/methdiff.cpp b/src/radmeth/methdiff.cpp index d7421637..9b2a5cc9 100644 --- a/src/radmeth/methdiff.cpp +++ b/src/radmeth/methdiff.cpp @@ -1,8 +1,8 @@ -/* methdiff: Computes probability that individual CpGs have higher - * methylation in file A than in file B, where files A and B - * are specified on the command line. +/* methdiff: Computes probability that individual CpGs have higher methylation + * in file A than in file B, where files A and B are specified on + * the command line. * - * Copyright (C) 2011-2023 Andrew D Smith + * Copyright (C) 2011-2025 Andrew D Smith * * Author: Andrew D. Smith * @@ -17,13 +17,13 @@ * General Public License for more details. */ +#include #include #include #include #include #include #include -#include #include "OptionParser.hpp" #include "smithlab_os.hpp" @@ -31,129 +31,126 @@ #include "MSite.hpp" -using std::cerr; -using std::cout; -using std::endl; -using std::min; -using std::runtime_error; -using std::string; -using std::vector; - -using std::ofstream; -using std::ostream_iterator; - using bamxx::bgzf_file; -struct quick_buf : public std::ostringstream, - public std::basic_stringbuf { - quick_buf() { static_cast &>(*this).rdbuf(this); } - void clear() { setp(pbase(), pbase()); } - char const *c_str() { - *pptr() = '\0'; - return pbase(); - } -}; - static inline double log_sum_log(const double p, const double q) { - if (p == 0) { return q; } - else if (q == 0) { return p; } + if (p == 0) + return q; + if (q == 0) + return p; const double larger = (p > q) ? p : q; const double smaller = (p > q) ? q : p; - return larger + log(1.0 + exp(smaller - larger)); + return larger + std::log1p(std::exp(smaller - larger)); } static inline double lnchoose(const unsigned int n, unsigned int m) { - if (m == n || m == 0) return 0; - if (m * 2 > n) m = n - m; - using std::lgamma; - return lgamma(n + 1.0) - lgamma(m + 1.0) - lgamma((n - m) + 1.0); + if (m == n || m == 0) + return 0; + if (m * 2 > n) + m = n - m; + return std::lgamma(n + 1.0) - std::lgamma(m + 1.0) - std::lgamma(n - m + 1.0); } static inline double -log_hyper_g_greater(const size_t meth_a, const size_t unmeth_a, - const size_t meth_b, const size_t unmeth_b, size_t k) { - return ( - lnchoose(meth_b + unmeth_b - 1, k) + - lnchoose(meth_a + unmeth_a - 1, meth_a + meth_b - 1 - k) - - lnchoose(meth_a + unmeth_a + meth_b + unmeth_b - 2, meth_a + meth_b - 1)); +log_hyper_g_greater(const std::size_t meth_a, const std::size_t unmeth_a, + const std::size_t meth_b, const std::size_t unmeth_b, + std::size_t k) { + return lnchoose(meth_b + unmeth_b - 1, k) + + lnchoose(meth_a + unmeth_a - 1, meth_a + meth_b - 1 - k) - + lnchoose(meth_a + unmeth_a + meth_b + unmeth_b - 2, + meth_a + meth_b - 1); } static inline double -test_greater_population(const size_t meth_a, const size_t unmeth_a, - const size_t meth_b, const size_t unmeth_b) { +test_greater_population(const std::size_t meth_a, const std::size_t unmeth_a, + const std::size_t meth_b, const std::size_t unmeth_b) { double p = 0; - for (size_t k = (meth_b > unmeth_a) ? meth_b - unmeth_a : 0; k < meth_b; ++k) + for (std::size_t k = (meth_b > unmeth_a) ? meth_b - unmeth_a : 0; k < meth_b; + ++k) p = log_sum_log(p, log_hyper_g_greater(meth_a, unmeth_a, meth_b, unmeth_b, k)); - return exp(p); + return std::exp(p); } -template T & +template +T & write_methdiff_site(T &out, const MSite &a, const MSite &b, const double diffscore) { + // static constexpr auto out_fmt = + // "%s\t%ld\t%c\t%s\t%.6g\t%ld\t%ld\t%ld\t%ld\n"; // clang-format off - return out << a.chrom << '\t' - << a.pos << '\t' - << a.strand << '\t' - << a.context << '\t' - << diffscore << '\t' - << a.n_meth() << '\t' - << a.n_unmeth() << '\t' - << b.n_meth() << '\t' - << b.n_unmeth() << '\n'; + static constexpr auto out_fmt = + "%s" // chrom + "\t" + "%ld" // pos + "\t" + "%c" // strand + "\t" + "%s" // context + "\t" + "%.6g" // diffscore + "\t" + "%ld" // a.n_meth() + "\t" + "%ld" // a.n_unmeth() + "\t" + "%ld" // b.n_meth() + "\t" + "%ld" // b.n_unmeth() + "\n"; // clang-format on -} + static constexpr auto buf_size = 1024; + static char buffer[buf_size]; -bgzf_file & -write_methdiff_site(bgzf_file &out, const MSite &a, const MSite &b, - const double diffscore) { - quick_buf buf; // clang-format off - buf << a.chrom << '\t' - << a.pos << '\t' - << a.strand << '\t' - << a.context << '\t' - << diffscore << '\t' - << a.n_meth() << '\t' - << a.n_unmeth() << '\t' - << b.n_meth() << '\t' - << b.n_unmeth() << '\n'; + const int r = std::sprintf(buffer, out_fmt, + a.chrom.data(), + a.pos, + a.strand, + a.context.data(), + diffscore, + a.n_meth(), + a.n_unmeth(), + b.n_meth(), + b.n_unmeth()); // clang-format on - out.write(buf.c_str(), buf.tellp()); + if (r < 0) + throw std::runtime_error("failed to write to output buffer"); + out.write(buffer, r); return out; } static inline bool -site_precedes(const size_t lhs_chrom_idx, const size_t lhs_pos, - const size_t rhs_chrom_idx, const size_t rhs_pos) { +site_precedes(const std::size_t lhs_chrom_idx, const std::size_t lhs_pos, + const std::size_t rhs_chrom_idx, const std::size_t rhs_pos) { return (lhs_chrom_idx < rhs_chrom_idx || (lhs_chrom_idx == rhs_chrom_idx && (lhs_pos < rhs_pos))); } -static size_t -get_chrom_id(std::unordered_map &chrom_order, - std::unordered_set &chroms_seen, const MSite &s) { - if (chroms_seen.find(s.chrom) != end(chroms_seen)) - throw runtime_error("unsorted chromosomes found"); +static std::size_t +get_chrom_id(std::unordered_map &chrom_order, + std::unordered_set &chroms_seen, const MSite &s) { + if (chroms_seen.find(s.chrom) != std::end(chroms_seen)) + throw std::runtime_error("unsorted chromosomes found"); chroms_seen.emplace(s.chrom); - auto x = chrom_order.find(s.chrom); - if (x == end(chrom_order)) { - const auto idx = chrom_order.size(); + const auto itr = chrom_order.find(s.chrom); + if (itr == std::cend(chrom_order)) { + const auto idx = std::size(chrom_order); chrom_order.emplace(s.chrom, idx); return idx; } - else return x->second; + return itr->second; } -static string -bad_order(const std::unordered_map &chrom_order, - const string prev_chrom, const size_t prev_pos, const string chrom, - const size_t pos) { +static std::string +bad_order(const std::unordered_map &chrom_order, + const std::string prev_chrom, const std::size_t prev_pos, + const std::string chrom, const std::size_t pos) { std::ostringstream oss; - const size_t chrom_id = chrom_order.find(chrom)->second; - const size_t prev_chrom_id = chrom_order.find(prev_chrom)->second; + const std::size_t chrom_id = chrom_order.find(chrom)->second; + const std::size_t prev_chrom_id = chrom_order.find(prev_chrom)->second; // clang-format off oss << "bad order:\n" << "chrom=" << chrom << " [id=" << chrom_id << "] " @@ -165,20 +162,21 @@ bad_order(const std::unordered_map &chrom_order, return oss.str(); } -template static void +template +static void process_sites(const bool show_progress, bgzf_file &in_a, bgzf_file &in_b, const bool allow_uncovered, const double pseudocount, T &out) { // chromosome order in the files - std::unordered_map chrom_order; - std::unordered_set chroms_seen_a, chroms_seen_b; + std::unordered_map chrom_order; + std::unordered_set chroms_seen_a, chroms_seen_b; MSite a, b; // ADS: default "pos" might be num lim max a.pos = 0; b.pos = 0; - string prev_chrom_a, prev_chrom_b; - size_t chrom_id_a = 0, chrom_id_b = 0; - size_t prev_chrom_id_a = 0, prev_chrom_id_b = 0; - size_t prev_pos_a = 0, prev_pos_b = 0; + std::string prev_chrom_a, prev_chrom_b; + std::size_t chrom_id_a = 0, chrom_id_b = 0; + std::size_t prev_chrom_id_a = 0, prev_chrom_id_b = 0; + std::size_t prev_pos_a = 0, prev_pos_b = 0; bool advance_a = true; bool advance_b = true; @@ -189,11 +187,12 @@ process_sites(const bool show_progress, bgzf_file &in_a, bgzf_file &in_b, if (prev_chrom_a.compare(a.chrom) != 0) { prev_chrom_id_a = chrom_id_a; chrom_id_a = get_chrom_id(chrom_order, chroms_seen_a, a); - if (show_progress) cerr << "processing " << a.chrom << endl; + if (show_progress) + std::cerr << "processing " << a.chrom << std::endl; prev_chrom_a = a.chrom; } if (site_precedes(chrom_id_a, a.pos, prev_chrom_id_a, prev_pos_a)) - throw runtime_error( + throw std::runtime_error( bad_order(chrom_order, prev_chrom_a, prev_pos_a, a.chrom, a.pos)); advance_a = site_precedes(chrom_id_a, a.pos, chrom_id_b, b.pos); prev_pos_a = a.pos; @@ -206,20 +205,21 @@ process_sites(const bool show_progress, bgzf_file &in_a, bgzf_file &in_b, prev_chrom_b = b.chrom; } if (site_precedes(chrom_id_b, b.pos, prev_chrom_id_b, prev_pos_b)) - throw runtime_error( + throw std::runtime_error( bad_order(chrom_order, prev_chrom_b, prev_pos_b, b.chrom, b.pos)); advance_b = site_precedes(chrom_id_b, b.pos, chrom_id_a, a.pos); prev_pos_b = b.pos; } - if (!in_a || !in_b) break; + if (!in_a || !in_b) + break; if (chrom_id_a == chrom_id_b && a.pos == b.pos) { - if (allow_uncovered || min(a.n_reads, b.n_reads) > 0) { - const size_t meth_a = a.n_meth() + pseudocount; - const size_t unmeth_a = a.n_unmeth() + pseudocount; - const size_t meth_b = b.n_meth() + pseudocount; - const size_t unmeth_b = b.n_unmeth() + pseudocount; + if (allow_uncovered || std::min(a.n_reads, b.n_reads) > 0) { + const std::size_t meth_a = a.n_meth() + pseudocount; + const std::size_t unmeth_a = a.n_unmeth() + pseudocount; + const std::size_t meth_b = b.n_meth() + pseudocount; + const std::size_t unmeth_b = b.n_unmeth() + pseudocount; const double diffscore = test_greater_population(meth_b, unmeth_b, meth_a, unmeth_a); @@ -239,12 +239,12 @@ process_sites(const bool show_progress, bgzf_file &in_a, bgzf_file &in_b, int main_methdiff(int argc, const char **argv) { try { - string outfile; + std::string outfile; double pseudocount = 1.0; // run mode flags bool allow_uncovered = true; - bool VERBOSE = false; + bool verbose = false; /****************** COMMAND LINE OPTIONS ********************/ OptionParser opt_parse(strip_path(argv[0]), @@ -257,54 +257,53 @@ main_methdiff(int argc, const char **argv) { "process only sites with coveage in both samples", false, allow_uncovered); opt_parse.add_opt("out", 'o', "output file", true, outfile); - opt_parse.add_opt("verbose", 'v', "print more run info", false, VERBOSE); - vector leftover_args; + opt_parse.add_opt("verbose", 'v', "print more run info", false, verbose); + 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() << std::endl + << opt_parse.about_message() << std::endl; return EXIT_SUCCESS; } if (opt_parse.about_requested()) { - cerr << opt_parse.about_message() << endl; + std::cerr << opt_parse.about_message() << std::endl; return EXIT_SUCCESS; } if (opt_parse.option_missing()) { - cerr << opt_parse.option_missing_message() << endl; + std::cerr << opt_parse.option_missing_message() << std::endl; return EXIT_SUCCESS; } if (leftover_args.size() != 2) { - cerr << opt_parse.help_message() << endl; + std::cerr << opt_parse.help_message() << std::endl; return EXIT_SUCCESS; } - const string cpgs_file_a = leftover_args[0]; - const string cpgs_file_b = leftover_args[1]; + const std::string cpgs_file_a = leftover_args[0]; + const std::string cpgs_file_b = leftover_args[1]; /****************** END COMMAND LINE OPTIONS *****************/ - if (VERBOSE) - cerr << "[opening counts file: " << cpgs_file_a << "]" << endl; + if (verbose) + std::cerr << "[opening counts file: " << cpgs_file_a << "]" << std::endl; bgzf_file in_a(cpgs_file_a, "r"); - if (!in_a) throw runtime_error("cannot open file: " + cpgs_file_a); + if (!in_a) + throw std::runtime_error("cannot open file: " + cpgs_file_a); - if (VERBOSE) - cerr << "[opening counts file: " << cpgs_file_b << "]" << endl; + if (verbose) + std::cerr << "[opening counts file: " << cpgs_file_b << "]" << std::endl; bgzf_file in_b(cpgs_file_b, "r"); - if (!in_b) throw runtime_error("cannot open file: " + cpgs_file_b); + if (!in_b) + throw std::runtime_error("cannot open file: " + cpgs_file_b); if (outfile.empty() || !has_gz_ext(outfile)) { - std::ofstream of; - if (!outfile.empty()) of.open(outfile); - std::ostream out(outfile.empty() ? cout.rdbuf() : of.rdbuf()); - - process_sites(VERBOSE, in_a, in_b, allow_uncovered, pseudocount, out); + std::ofstream out(outfile); + process_sites(verbose, in_a, in_b, allow_uncovered, pseudocount, out); } else { bgzf_file out(outfile, "w"); - process_sites(VERBOSE, in_a, in_b, allow_uncovered, pseudocount, out); + process_sites(verbose, in_a, in_b, allow_uncovered, pseudocount, out); } } - catch (const runtime_error &e) { - cerr << e.what() << endl; + catch (const std::exception &e) { + std::cerr << e.what() << std::endl; return EXIT_FAILURE; } return EXIT_SUCCESS;