diff --git a/src/analysis/cpgbins.cpp b/src/analysis/cpgbins.cpp index 7c2ecdb4..949204d0 100644 --- a/src/analysis/cpgbins.cpp +++ b/src/analysis/cpgbins.cpp @@ -37,25 +37,11 @@ #include #include -using std::cerr; -using std::cout; -using std::distance; -using std::endl; -using std::ifstream; -using std::ostream; -using std::runtime_error; -using std::size; -using std::string; -using std::uint32_t; -using std::uint64_t; -using std::unordered_map; -using std::vector; - using bamxx::bgzf_file; namespace fs = std::filesystem; -static string +static std::string format_levels_counter(const LevelsCounter &lc) { // ... // (7) weighted mean methylation @@ -73,53 +59,57 @@ 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.total_c + lc.total_t) << '\t' + << lc.total_meth << '\t' + << lc.sites_covered; // clang-format on return oss.str(); } -static unordered_map -get_chrom_sizes(const string &chrom_sizes_file) { - unordered_map chrom_sizes; +static std::unordered_map +get_chrom_sizes(const std::string &chrom_sizes_file) { + std::unordered_map chrom_sizes; - ifstream in(chrom_sizes_file); + std::ifstream in(chrom_sizes_file); if (!in) - throw runtime_error("failed to open file: " + chrom_sizes_file); + throw std::runtime_error("failed to open file: " + chrom_sizes_file); - string line; + std::string line; while (getline(in, line)) { std::istringstream iss(line); - string chrom_name{}; - uint64_t chrom_size{}; + std::string chrom_name{}; + std::uint64_t chrom_size{}; if (!(iss >> chrom_name >> chrom_size)) - throw runtime_error("bad line in " + chrom_sizes_file + ":\n" + line); - string dummy; + throw std::runtime_error("bad line in " + chrom_sizes_file + ":\n" + + line); + std::string dummy; if (iss >> dummy) - throw runtime_error("too many columns: " + line); + throw std::runtime_error("too many columns: " + line); - if (chrom_sizes.find(chrom_name) != cend(chrom_sizes)) - throw runtime_error("repeated entry " + chrom_name + " in " + - chrom_sizes_file); + if (chrom_sizes.find(chrom_name) != std::cend(chrom_sizes)) + throw std::runtime_error("repeated entry " + chrom_name + " in " + + chrom_sizes_file); chrom_sizes[chrom_name] = chrom_size; } return chrom_sizes; } -static vector -get_chrom_names(const string &chrom_sizes_file) { - ifstream in(chrom_sizes_file); +static std::vector +get_chrom_names(const std::string &chrom_sizes_file) { + std::ifstream in(chrom_sizes_file); if (!in) - throw runtime_error("failed to open file: " + chrom_sizes_file); + throw std::runtime_error("failed to open file: " + chrom_sizes_file); - vector chrom_names; + std::vector chrom_names; - string line; + std::string line; while (getline(in, line)) { std::istringstream iss(line); - string chrom_name{}; + std::string chrom_name{}; if (!(iss >> chrom_name)) - throw runtime_error("bad line in " + chrom_sizes_file + ":\n" + line); + throw std::runtime_error("bad line in " + chrom_sizes_file + ":\n" + + line); chrom_names.push_back(chrom_name); } @@ -128,13 +118,13 @@ get_chrom_names(const string &chrom_sizes_file) { static void update(LevelsCounter &lc, const xcounts_entry &xce) { - const uint64_t n_reads = xce.n_meth + xce.n_unmeth; + const std::uint64_t n_reads = xce.n_reads(); if (n_reads > 0) { ++lc.sites_covered; lc.max_depth = std::max(lc.max_depth, n_reads); lc.total_c += xce.n_meth; lc.total_t += xce.n_unmeth; - const auto meth = static_cast(xce.n_unmeth) / n_reads; + const auto meth = xce.frac(); lc.total_meth += meth; double lower = 0.0, upper = 0.0; wilson_ci_for_binomial(lc.alpha, n_reads, meth, lower, upper); @@ -144,20 +134,33 @@ update(LevelsCounter &lc, const xcounts_entry &xce) { ++lc.total_sites; } +[[nodiscard]] static std::string +get_name(const bool report_more_info, const char level_code, + const LevelsCounter &lc) { + static const std::string tag = "CpG"; + return tag + (report_more_info + ? "" + : "_" + std::to_string( + (level_code == 'w' + ? lc.coverage() + : (level_code == 'u' ? lc.sites_covered + : lc.total_called())))); +} + static void process_chrom(const bool report_more_info, const char level_code, - const string &chrom_name, const uint64_t chrom_size, - const uint64_t bin_size, const vector &sites, - ostream &out) { + const std::string &chrom_name, const std::uint64_t chrom_size, + const std::uint64_t bin_size, + const std::vector &sites, std::ostream &out) { GenomicRegion r(chrom_name, 0, 0, "CpG", 0.0, '+'); - uint64_t j = 0; + std::uint64_t j = 0; for (auto i = 0ul; i < chrom_size; i += bin_size) { - while (j < size(sites) && sites[j].pos < i) + while (j < std::size(sites) && sites[j].pos < i) ++j; LevelsCounter lc; - while (j < size(sites) && sites[j].pos < i + bin_size) + while (j < std::size(sites) && sites[j].pos < i + bin_size) update(lc, sites[j++]); r.set_start(i); @@ -165,11 +168,9 @@ process_chrom(const bool report_more_info, const char level_code, r.set_score(level_code == 'w' ? lc.mean_meth_weighted() : (level_code == 'u' ? lc.mean_meth() : lc.fractional_meth())); - r.set_name("CpG_" + - std::to_string((level_code == 'w' - ? lc.coverage() - : (level_code == 'u' ? lc.sites_covered - : lc.total_called())))); + + r.set_name(get_name(report_more_info, level_code, lc)); + out << r; if (report_more_info) out << '\t' << format_levels_counter(lc); @@ -178,12 +179,12 @@ process_chrom(const bool report_more_info, const char level_code, } static void -process_chrom(const bool report_more_info, const string &chrom_name, - const uint64_t chrom_size, const uint64_t bin_size, - ostream &out) { +process_chrom(const bool report_more_info, const std::string &chrom_name, + const std::uint64_t chrom_size, const std::uint64_t bin_size, + std::ostream &out) { GenomicRegion r(chrom_name, 0, 0, "CpG_0", 0.0, '+'); LevelsCounter lc; - const string lc_formatted = format_levels_counter(lc); + const std::string lc_formatted = format_levels_counter(lc); for (auto i = 0ul; i < chrom_size; i += bin_size) { r.set_start(i); r.set_end(std::min(i + bin_size, chrom_size)); @@ -197,7 +198,7 @@ process_chrom(const bool report_more_info, const string &chrom_name, int main_cpgbins(int argc, char *argv[]) { try { - static const string description = R"""( + static const std::string description = R"""( Compute average site methylation levels in each non-overlapping genomic bin of the specified size. The 5th column (the "score" column in BED format) is determined by the '-l' or '-level' argument. @@ -212,22 +213,22 @@ 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 report_more_info = false; - uint32_t n_threads = 1; - // uint64_t bin_size = 1000; + std::uint32_t n_threads = 1; + // std::uint64_t bin_size = 1000; size_t bin_size = 1000; // ADS: for macOS gcc-14.2.0 - string level_code = "w"; - string outfile; + std::string level_code = "w"; + std::string outfile; /****************** COMMAND LINE OPTIONS ********************/ OptionParser opt_parse(strip_path(argv[0]), description, " "); opt_parse.set_show_defaults(); opt_parse.add_opt("output", 'o', "name of output file (default: stdout)", - false, outfile); + true, outfile); opt_parse.add_opt("bin", 'b', "bin size in base pairs", false, bin_size); opt_parse.add_opt("threads", 't', "number of threads", false, n_threads); opt_parse.add_opt("level", 'l', @@ -237,54 +238,51 @@ Columns (beyond the first 6) in the BED format output: opt_parse.add_opt("more-levels", 'M', "report more methylation information", false, report_more_info); 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() || opt_parse.about_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.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 chrom_sizes_file = leftover_args.front(); - const string xcounts_file = leftover_args.back(); + const std::string chrom_sizes_file = leftover_args.front(); + const std::string xcounts_file = leftover_args.back(); /****************** END COMMAND LINE OPTIONS *****************/ if (!fs::is_regular_file(chrom_sizes_file)) - throw runtime_error("chromosome sizes file not a regular file: " + - chrom_sizes_file); + throw std::runtime_error("chromosome sizes file not a regular file: " + + chrom_sizes_file); if (!fs::is_regular_file(xcounts_file)) - throw runtime_error("xsym file not a regular file: " + xcounts_file); + throw std::runtime_error("xsym file not a regular file: " + xcounts_file); const auto sites_by_chrom = read_xcounts_by_chrom(n_threads, xcounts_file); const auto chrom_names = get_chrom_names(chrom_sizes_file); const auto chrom_sizes = get_chrom_sizes(chrom_sizes_file); - std::ofstream of; - if (!outfile.empty()) { - of.open(outfile); - if (!of) - throw runtime_error("failed to open outfile: " + outfile); - } - std::ostream out(outfile.empty() ? cout.rdbuf() : of.rdbuf()); + std::ofstream out(outfile); + if (!out) + throw std::runtime_error("failed to open outfile: " + outfile); for (const auto &chrom_name : chrom_names) { const auto sites = sites_by_chrom.find(chrom_name); const auto chrom_size = chrom_sizes.find(chrom_name); - if (chrom_size == cend(chrom_sizes)) - throw runtime_error("failed chrom size lookup"); - if (sites != cend(sites_by_chrom)) + if (chrom_size == std::cend(chrom_sizes)) + throw std::runtime_error("failed chrom size lookup"); + if (sites != std::cend(sites_by_chrom)) process_chrom(report_more_info, level_code[0], chrom_name, chrom_size->second, bin_size, sites->second, out); else @@ -293,7 +291,7 @@ Columns (beyond the first 6) in the BED format output: } } catch (const std::exception &e) { - cerr << e.what() << endl; + std::cerr << e.what() << '\n'; return EXIT_FAILURE; } return EXIT_SUCCESS; diff --git a/src/common/LevelsCounter.hpp b/src/common/LevelsCounter.hpp index fe44fbb9..b264a66c 100644 --- a/src/common/LevelsCounter.hpp +++ b/src/common/LevelsCounter.hpp @@ -2,15 +2,15 @@ * * 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. */ #ifndef LEVELS_COUNTER_HPP @@ -28,50 +28,49 @@ struct LevelsCounter { // tools provide more information about this context string. std::string context; - // total_sites is the total number of sites summarized, all of which - // are consistent with the context string. If the context string is - // empty or not used, then any constraints on which sites are - // counted among this total is defined separately. + // total_sites is the total number of sites summarized, all of which are + // consistent with the context string. If the context string is empty or not + // used, then any constraints on which sites are counted among this total is + // defined separately. std::uint64_t total_sites{}; // sites_covered is the number of reference genome sites analyzed // that have at least one mapped read covering them. std::uint64_t sites_covered{}; - // max_depth is the number of reads covering the reference genome - // site having the most mapped reads covering it. + // max_depth is the number of reads covering the reference genome site + // having the most mapped reads covering it. std::uint64_t max_depth{}; - // mutations is the number of sites tagged as having been mutated in - // the sample based on the mapped reads in comparisin with the - // underlying reference genome site. This may be indicated in - // various ways and defined in various ways. + // mutations is the number of sites tagged as having been mutated in the + // sample based on the mapped reads in comparisin with the underlying + // reference genome site. This may be indicated in various ways and defined + // in various ways. std::uint64_t mutations{}; - // total_c is the number of cytosines in mapped reads covering - // relevant cytosines in the reference genome. + // total_c is the number of cytosines in mapped reads covering relevant + // cytosines in the reference genome. std::uint64_t total_c{}; - // total_t is the number of thymines in mapped reads covering - // relevant cytosines in the reference genome. + // total_t is the number of thymines in mapped reads covering relevant + // cytosines in the reference genome. std::uint64_t total_t{}; - // called_meth is the number of cytosines in the reference genome - // "called" as being methylated using some statistical criteria. + // called_meth is the number of cytosines in the reference genome "called" + // as being methylated using some statistical criteria. std::uint64_t called_meth{}; - // called_unmeth is the number of cytosines in the reference genome - // "called" as being unmethylated using some statistical criteria. + // called_unmeth is the number of cytosines in the reference genome "called" + // as being unmethylated using some statistical criteria. std::uint64_t called_unmeth{}; - // total_meth is the sum over all relevant cytosines of the - // methylation level estimated at that cytosine using the mapped - // reads. This value contributes to the unweighted mean methylation. + // total_meth is the sum over all relevant cytosines of the methylation + // level estimated at that cytosine using the mapped reads. This value + // contributes to the unweighted mean methylation. double total_meth{}; - // coverage is equal to total_c plus total_t. These are the counts - // that contribute towards estimates of the various methylation - // levels. + // coverage is equal to total_c plus total_t. These are the counts that + // contribute towards estimates of the various methylation levels. [[nodiscard]] std::uint64_t coverage() const { return total_c + total_t; @@ -92,24 +91,23 @@ struct LevelsCounter { std::max(coverage(), static_cast(1)); } + // mean_meth is the unweighted mean methylation level. This is the ratio of + // total_meth divided by sites_covered. Value is always between 0 and 1. + [[nodiscard]] double + mean_meth() const { + return static_cast(total_meth) / + std::max(sites_covered, static_cast(1)); + } + // fractional_meth is the fraction of "called" sites that are called - // methylated. It is the ratio of called_meth divided by - // total_called. This value is always between 0 and 1. + // methylated. It is the ratio of called_meth divided by total_called. This + // value is always between 0 and 1. [[nodiscard]] double fractional_meth() const { return static_cast(called_meth) / std::max(total_called(), static_cast(1)); } - // mean_meth is the unweighted mean methylation level. This is the - // ratio of total_meth divided by sites_covered. This value is - // always between 0 and 1. - [[nodiscard]] double - mean_meth() const { - return static_cast(total_meth) / - std::max(sites_covered, static_cast(1)); - } - explicit LevelsCounter(const std::string &c) : context{c} {} LevelsCounter() = default; diff --git a/src/common/bsutils.cpp b/src/common/bsutils.cpp index 867be23c..7d0a661e 100644 --- a/src/common/bsutils.cpp +++ b/src/common/bsutils.cpp @@ -1,52 +1,54 @@ -/* - * Copyright (C) 2018-2022 Andrew D. Smith +/* Copyright (C) 2018-2025 Andrew D. Smith * - * Authors: Andrew D. Smith + * Author: Andrew D. Smith * - * 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 "bsutils.hpp" - -#include -#include -#include -#include +#include "dnmtools_gaussinv.hpp" #include #include -#include "dnmtools_gaussinv.hpp" -//// CONFIDENCE INTERVALS //**************//////////////////////// +#include +#include +#include +#include + void -wilson_ci_for_binomial(const double alpha, const double n, - const double p_hat, double &lower, double &upper) { - const double z = dnmt_gsl_cdf_ugaussian_Pinv(1 - alpha/2); - const double denom = 1 + z*z/n; - const double first_term = p_hat + z*z/(2*n); - const double discriminant = p_hat*(1 - p_hat)/n + z*z/(4*n*n); - lower = std::max(0.0, (first_term - z*std::sqrt(discriminant))/denom); - upper = std::min(1.0, (first_term + z*std::sqrt(discriminant))/denom); +wilson_ci_for_binomial(const double alpha, const double n, const double p_hat, + double &lower, double &upper) { + if (n <= 0.0) { // protection + lower = 0.0; + upper = 1.0; + return; + } + const double z = dnmt_gsl_cdf_ugaussian_Pinv(1 - alpha / 2); + const double denom = 1 + z * z / n; + const double first_term = p_hat + z * z / (2 * n); + const double discriminant = + std::max(0.0, p_hat * (1 - p_hat) / n + z * z / (4 * n * n)); + lower = std::max(0.0, (first_term - z * std::sqrt(discriminant)) / denom); + upper = std::min(1.0, (first_term + z * std::sqrt(discriminant)) / denom); } -//////////////////////////**************//////////////////////// - void -adjust_region_ends(const std::vector > &clusters, +adjust_region_ends(const std::vector> &clusters, std::vector ®ions) { - assert(clusters.size() == regions.size()); - for (size_t i = 0; i < regions.size(); ++i) { - size_t max_pos = regions[i].get_end(); - size_t min_pos = regions[i].get_start(); - for (size_t j = 0; j < clusters[i].size(); ++j) { + assert(std::size(clusters) == std::size(regions)); + for (std::size_t i = 0; i < std::size(regions); ++i) { + std::size_t max_pos = regions[i].get_end(); + std::size_t min_pos = regions[i].get_start(); + for (std::size_t j = 0; j < std::size(clusters[i]); ++j) { max_pos = std::max(clusters[i][j].get_end(), max_pos); min_pos = std::min(clusters[i][j].get_start(), min_pos); } @@ -55,20 +57,18 @@ adjust_region_ends(const std::vector > &clusters, } } - void relative_sort(const std::vector &mapped_locations, const std::vector &names, - std::vector &lookup) { + std::vector &lookup) { - std::unordered_map names_map; - for (size_t i = 0; i < names.size(); ++i) + std::unordered_map names_map; + for (std::size_t i = 0; i < std::size(names); ++i) names_map[names[i]] = i; - for (size_t i = 0; i < mapped_locations.size(); ++i) { - const std::unordered_map::const_iterator - j(names_map.find(mapped_locations[i].get_name())); - if (j == names_map.end()) + for (std::size_t i = 0; i < std::size(mapped_locations); ++i) { + const auto j = names_map.find(mapped_locations[i].get_name()); + if (j == std::cend(names_map)) throw std::runtime_error("read sequence not found for: " + names[i]); lookup.push_back(j->second); } diff --git a/src/common/xcounts_utils.hpp b/src/common/xcounts_utils.hpp index b4a14a49..91f0ae61 100644 --- a/src/common/xcounts_utils.hpp +++ b/src/common/xcounts_utils.hpp @@ -19,19 +19,36 @@ #ifndef XCOUNTS_UTILS_HPP #define XCOUNTS_UTILS_HPP -#include -#include #include +#include +#include #include +#include struct xcounts_entry { - uint64_t pos{}; // absolute position - uint32_t n_meth{}; - uint32_t n_unmeth{}; + std::uint64_t pos{}; // absolute position + std::uint32_t n_meth{}; + std::uint32_t n_unmeth{}; + + [[nodiscard]] std::uint32_t + n_reads() const { + return n_meth + n_unmeth; + } + + [[nodiscard]] double + frac() const { + return static_cast(n_meth) / n_reads(); + } }; +inline std::ostream & +operator<<(std::ostream &o, const xcounts_entry &e) { + return o << e.pos << '\t' << e.n_meth << '\t' << e.n_unmeth; +} + std::unordered_map> -read_xcounts_by_chrom(const uint32_t n_threads, const std::string &xcounts_file); +read_xcounts_by_chrom(const std::uint32_t n_threads, + const std::string &xcounts_file); bool get_is_xcounts_file(const std::string &filename);