diff --git a/Makefile.am b/Makefile.am index af091f3c..58a1c915 100644 --- a/Makefile.am +++ b/Makefile.am @@ -189,8 +189,12 @@ libdnmtools_a_SOURCES += \ # ADS: additional radmeth sources to help isolate the parts using GSL for # optimization libdnmtools_a_SOURCES += \ - src/radmeth/radmeth_optimize.hpp \ - src/radmeth/radmeth_model.hpp + src/radmeth/radmeth_optimize_series.hpp \ + src/radmeth/radmeth_optimize_gamma.hpp \ + src/radmeth/radmeth_optimize_params.hpp \ + src/radmeth/radmeth_model.hpp \ + src/radmeth/radmeth_utils.hpp \ + src/radmeth/radmeth_design.hpp LDADD = libdnmtools.a src/abismal/libabismal.a src/smithlab_cpp/libsmithlab_cpp.a @@ -236,9 +240,12 @@ dnmtools_SOURCES += src/amrfinder/amrtester.cpp dnmtools_SOURCES += src/radmeth/dmr.cpp dnmtools_SOURCES += src/radmeth/methdiff.cpp -dnmtools_SOURCES += src/radmeth/radmeth_optimize.cpp -dnmtools_SOURCES += src/radmeth/radmeth_model.cpp +dnmtools_SOURCES += src/radmeth/radmeth_utils.cpp +dnmtools_SOURCES += src/radmeth/radmeth_optimize_series.cpp +dnmtools_SOURCES += src/radmeth/radmeth_optimize_gamma.cpp +dnmtools_SOURCES += src/radmeth/radmeth_design.cpp dnmtools_SOURCES += src/radmeth/radmeth.cpp +dnmtools_SOURCES += src/radmeth/radmeth_nano.cpp dnmtools_SOURCES += src/radmeth/radmeth-adjust.cpp dnmtools_SOURCES += src/radmeth/radmeth-merge.cpp diff --git a/src/analysis/nanopore.cpp b/src/analysis/nanopore.cpp index 6eb3ba53..cdf29b4c 100644 --- a/src/analysis/nanopore.cpp +++ b/src/analysis/nanopore.cpp @@ -845,7 +845,8 @@ struct read_processor { if (chrom_name == nullptr) throw std::runtime_error("failed to identify chrom for tid: " + std::to_string(tid)); - const auto chrom_name_offset = std::sprintf(buffer, "%s\t", chrom_name); + const auto chrom_name_offset = + std::snprintf(buffer, buf_size, "%s\t", chrom_name); if (chrom_name_offset < 0) throw std::runtime_error("failed to write to output buffer"); auto buffer_after_chrom = buffer + chrom_name_offset; @@ -869,14 +870,16 @@ struct read_processor { const double mods = counts[chrom_posn].get_mods(is_c) / denom; // clang-format off - const int r = std::sprintf(buffer_after_chrom, out_fmt, - chrom_posn, - is_c ? '+' : '-', - tag_values[the_tag], - mods, - n_reads, - hydroxy, - methyl); + const int r = std::snprintf(buffer_after_chrom, + buf_size - chrom_name_offset, + out_fmt, + chrom_posn, + is_c ? '+' : '-', + tag_values[the_tag], + mods, + n_reads, + hydroxy, + methyl); // clang-format on if (r < 0) @@ -900,7 +903,8 @@ struct read_processor { if (chrom_name == nullptr) throw std::runtime_error("failed to identify chrom for tid: " + std::to_string(tid)); - const auto chrom_name_offset = std::sprintf(buffer, "%s\t", chrom_name); + const auto chrom_name_offset = + std::snprintf(buffer, buf_size, "%s\t", chrom_name); if (chrom_name_offset < 0) throw std::runtime_error("failed to write to output buffer"); auto buffer_after_chrom = buffer + chrom_name_offset; @@ -935,12 +939,14 @@ struct read_processor { const auto mods = hydroxy + methyl; // clang-format off - const int r = std::sprintf(buffer_after_chrom, out_fmt, - chrom_posn - 1, // for previous position - mods, - n_reads, - hydroxy, - methyl); + const int r = std::snprintf(buffer_after_chrom, + buf_size - chrom_name_offset, + out_fmt, + chrom_posn - 1, // for previous position + mods, + n_reads, + hydroxy, + methyl); // clang-format on if (r < 0) diff --git a/src/dnmtools.cpp b/src/dnmtools.cpp index b4bb81fd..4c56673c 100644 --- a/src/dnmtools.cpp +++ b/src/dnmtools.cpp @@ -111,6 +111,8 @@ main_radmeth_adjust(int argc, char *argv[]); int main_radmeth(int argc, char *argv[]); int +main_radmeth_nano(int argc, char *argv[]); +int main_radmeth_merge(int argc, char *argv[]); int main_clean_hairpins(int argc, char *argv[]); @@ -204,6 +206,7 @@ main(int argc, char *argv[]) { {{{"dmr", "identify DMRs from genomic intervals and single-CpG DM probabilities", main_dmr}, {"diff", "compute single-CpG DM probability between two methylomes", main_methdiff}, {"radmeth", "compute DM probabilities for each CpG using multiple methylomes", main_radmeth}, + {"radmeth-nano", "radmeth designed for nanopore data", main_radmeth_nano}, {"radadjust", "adjust p-values from radmeth output", main_radmeth_adjust}, {"radmerge", "merge significant CpGs in radmeth output", main_radmeth_merge}}}}, diff --git a/src/radmeth/methdiff.cpp b/src/radmeth/methdiff.cpp index 2fe12816..33d4d474 100644 --- a/src/radmeth/methdiff.cpp +++ b/src/radmeth/methdiff.cpp @@ -105,16 +105,16 @@ write_methdiff_site(T &out, const MSite &a, const MSite &b, static char buffer[buf_size]; // clang-format off - 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()); + const int r = std::snprintf(buffer, buf_size, 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 if (r < 0) throw std::runtime_error("failed to write to output buffer"); diff --git a/src/radmeth/radmeth.cpp b/src/radmeth/radmeth.cpp index fcf9bae1..68eab100 100644 --- a/src/radmeth/radmeth.cpp +++ b/src/radmeth/radmeth.cpp @@ -14,8 +14,11 @@ * more details. */ +#include "radmeth_design.hpp" #include "radmeth_model.hpp" -#include "radmeth_optimize.hpp" +#include "radmeth_optimize_params.hpp" +#include "radmeth_optimize_series.hpp" +#include "radmeth_utils.hpp" // smithlab headers #include "GenomicRegion.hpp" @@ -26,7 +29,6 @@ #include #include #include -#include #include #include #include @@ -36,169 +38,25 @@ #include #include -[[nodiscard]] static std::string -format_duration(const std::chrono::duration elapsed) { - static constexpr auto s_per_h = 3600; - static constexpr auto s_per_m = 60; - const double tot_s = elapsed.count(); - - // break down into hours, minutes, seconds - const std::uint32_t hours = tot_s / 3600; - const std::uint32_t minutes = (static_cast(tot_s) % s_per_h) / s_per_m; - const double seconds = tot_s - (hours * s_per_h) - (minutes * s_per_m); - - std::ostringstream oss; - oss << std::setfill('0') << std::setw(2) << hours << ":" << std::setfill('0') - << std::setw(2) << minutes << ":" << std::fixed << std::setprecision(2) - << std::setw(5) << seconds; - return oss.str(); -} - -struct file_progress { - double one_thousand_over_filesize{}; - std::size_t prev_offset{}; - explicit file_progress(const std::string &filename) : - one_thousand_over_filesize{1000.0 / std::filesystem::file_size(filename)} {} - void - operator()(std::ifstream &in) { - const std::size_t curr_offset = - in.eof() ? 1000 : in.tellg() * one_thousand_over_filesize; - if (curr_offset <= prev_offset) - return; - std::ios old_state(nullptr); - old_state.copyfmt(std::cerr); - std::cerr << "\r[progress: " << std::setw(5) << std::fixed - << std::setprecision(1) << (curr_offset / 10.0) - << (curr_offset == 1000 ? "%]\n" : "%]"); - std::cerr.copyfmt(old_state); - prev_offset = (curr_offset == 1000) - ? std::numeric_limits::max() - : curr_offset; - } -}; - +template [[nodiscard]] static bool -consistent_sample_names(const Regression ®, const std::string &header) { - std::istringstream iss(header); - auto nm_itr(std::begin(reg.design.sample_names)); - const auto nm_end(std::end(reg.design.sample_names)); - std::string token; - while (iss >> token && nm_itr != nm_end) - if (token != *nm_itr++) - return false; - return true; -} - -[[nodiscard]] static std::vector -get_sample_names_from_header(const std::string &header) { - std::istringstream iss(header); - std::string token; - std::vector sample_names; - while (iss >> token) - sample_names.push_back(token); - return sample_names; -} - -// Series representation for the lower incomplete gamma P(a,x) -[[nodiscard]] static double -gamma_p_series(const double a, const double x) { - constexpr auto eps = std::numeric_limits::epsilon(); - constexpr auto max_iter = 100; - double sum = 1.0 / a; - double term = sum; - for (auto n = 1; n < max_iter; ++n) { - term *= x / (a + n); - sum += term; - if (term < eps * sum) - break; - } - return sum * std::exp(-x + a * std::log(x) - std::lgamma(a)); -} - -[[nodiscard]] static inline double -safe_floor(const double x, const double floor_val) { - return std::abs(x) < floor_val ? floor_val : x; -} - -// Continued fraction representation for the upper incomplete gamma Q(a,x) -[[nodiscard]] static double -gamma_q_contfrac(const double a, const double x) { - constexpr auto epsilon = std::numeric_limits::epsilon(); - constexpr auto fpmin = std::numeric_limits::min() / epsilon; - constexpr auto max_iter = 100; - double b = x + 1 - a; - double c = 1.0 / fpmin; - double d = 1.0 / b; - double h = d; - for (auto i = 1; i < max_iter; ++i) { - const double an = -i * (i - a); - b += 2; - d = safe_floor(an * d + b, fpmin); - c = safe_floor(b + an / c, fpmin); - d = 1.0 / d; - const double delta = d * c; - h *= delta; - if (std::abs(delta - 1.0) < epsilon) - break; - } - return std::exp(-x + a * std::log(x) - std::lgamma(a)) * h; -} - -// Regularized lower incomplete gamma P(a,x) -[[nodiscard]] static double -gamma_p(const double a, const double x) { - if (x < 0 || a <= 0) - return 0.0; - if (x == 0) - return 0.0; - if (x < a + 1.0) - return gamma_p_series(a, x); - return 1.0 - gamma_q_contfrac(a, x); -} - -// chi-square CDF: P(k/2, x/2) -[[nodiscard]] static double -chi_square_cdf(const double x, const double k) { - return gamma_p(k * 0.5, x * 0.5); -} - -// Given the maximum likelihood estimates of the full and reduced models, the -// function outputs the p-value of the log-likelihood ratio. *Note* that it is -// assumed that the reduced model has one fewer factor than the reduced model. -[[nodiscard]] static double -llr_test(const double null_loglik, const double full_loglik) { - // The log-likelihood ratio statistic. - const double log_lik_stat = -2 * (null_loglik - full_loglik); - - // It is assumed that null model has one fewer factor than the full - // model. Hence the number of degrees of freedom is 1. - const std::size_t degrees_of_freedom = 1; - - // Log-likelihood ratio statistic has a chi-sqare distribution. - const double chisq_p = chi_square_cdf(log_lik_stat, degrees_of_freedom); - - const double p_value = 1.0 - chisq_p; - - return p_value; -} - -static bool -has_low_coverage(const Regression ®, const std::size_t test_factor) { +has_low_coverage(const RegressionType ®, const std::size_t test_factor) { bool cvrd_in_test_fact_smpls = false; const auto &tcol = reg.design.tmatrix[test_factor]; for (std::size_t i = 0; i < reg.n_samples() && !cvrd_in_test_fact_smpls; ++i) - cvrd_in_test_fact_smpls = (tcol[i] == 1 && reg.props.mc[i].n_reads != 0); + cvrd_in_test_fact_smpls = (tcol[i] == 1 && reg.mc[i].n_reads != 0); bool cvrd_in_other_smpls = false; for (std::size_t i = 0; i < reg.n_samples() && !cvrd_in_other_smpls; ++i) - cvrd_in_other_smpls = (tcol[i] != 1 && reg.props.mc[i].n_reads != 0); + cvrd_in_other_smpls = (tcol[i] != 1 && reg.mc[i].n_reads != 0); return !cvrd_in_test_fact_smpls || !cvrd_in_other_smpls; } +template [[nodiscard]] static bool -has_extreme_counts(const Regression ®) { - const auto &mc = reg.props.mc; +has_extreme_counts(const RegressionType ®) { + const auto &mc = reg.mc; bool full_meth = true; for (auto i = std::cbegin(mc); i != std::cend(mc) && full_meth; ++i) @@ -211,42 +69,6 @@ has_extreme_counts(const Regression ®) { return full_meth || no_meth; } -[[nodiscard]] static bool -has_two_values(const Regression ®, const std::size_t test_factor) { - const auto &tcol = reg.design.tmatrix[test_factor]; - for (const auto x : tcol) - if (x != tcol[0]) - return true; - return false; -} - -[[nodiscard]] static std::uint32_t -get_test_factor_idx(const Regression &model, const std::string &test_factor) { - const auto &factors = model.design.factor_names; - const auto itr = - std::find(std::cbegin(factors), std::cend(factors), test_factor); - - if (itr == std::cend(factors)) - throw std::runtime_error("factor not part of design: " + test_factor); - - return std::distance(std::cbegin(factors), itr); -} - -[[nodiscard]] static Design -read_design(const std::string &design_filename) { - std::ifstream design_file(design_filename); - if (!design_file) - throw std::runtime_error("could not open file: " + design_filename); - Design design; - design_file >> design; - return design; -} - -[[nodiscard]] static inline double -overdispersion_factor(const std::uint32_t n_samples, const double dispersion) { - return (n_samples - 1) / (dispersion + 1); -} - enum class row_status : std::uint8_t { ok, na, @@ -254,11 +76,12 @@ enum class row_status : std::uint8_t { na_extreme_cnt, }; +template static void radmeth(const bool show_progress, const bool more_na_info, const std::uint32_t n_threads, const std::string &table_filename, - const std::string &outfile, const Regression &alt_model, - const Regression &null_model, const std::uint32_t test_factor_idx) { + const std::string &outfile, const Regression &alt_model, + const Regression &null_model, const std::uint32_t test_factor_idx) { static constexpr auto buf_size = 1024; static constexpr auto max_lines = 16384; @@ -274,7 +97,7 @@ radmeth(const bool show_progress, const bool more_na_info, std::getline(table_file, sample_names_header); const auto sample_names = get_sample_names_from_header(sample_names_header); - if (!consistent_sample_names(alt_model, sample_names_header)) { + if (!consistent_sample_names(alt_model.design, sample_names_header)) { // clang-format off const auto message = R"( @@ -298,8 +121,8 @@ that the design matrix and the proportion table are correctly formatted. std::vector n_bytes(max_lines, 0); std::vector lines(max_lines); - std::vector alt_models(n_threads, alt_model); - std::vector null_models(n_threads, null_model); + std::vector> alt_models(n_threads, alt_model); + std::vector> null_models(n_threads, null_model); const auto n_groups = alt_model.n_groups(); @@ -328,7 +151,7 @@ that the design matrix and the proportion table are correctly formatted. // balances work better. if (b % n_threads != thread_id) continue; - t_alt_model.props.parse(lines[b]); + t_alt_model.parse(lines[b]); if (t_alt_model.props_size() != n_samples) throw std::runtime_error("found row with wrong number of columns"); @@ -344,7 +167,8 @@ that the design matrix and the proportion table are correctly formatted. fit_regression_model(t_alt_model, p_estim_alt, phi_estim_alt); - t_null_model.props = t_alt_model.props; + t_null_model.mc = t_alt_model.mc; + t_null_model.rowname = t_alt_model.rowname; fit_regression_model(t_null_model, p_estim_null, phi_estim_null); @@ -359,48 +183,54 @@ that the design matrix and the proportion table are correctly formatted. const auto status = std::get<1>(p_val_status); n_bytes[b] = [&] { + auto bufsize = std::size(bufs[b]); // clang-format off const int n_prefix_bytes = - std::sprintf(bufs[b].data(), "%s\t", t_alt_model.rowname().data()); + std::snprintf(bufs[b].data(), bufsize, + "%s\t", t_alt_model.rowname.data()); // clang-format on if (n_prefix_bytes < 0) return n_prefix_bytes; + bufsize -= n_prefix_bytes; auto cursor = bufs[b].data() + n_prefix_bytes; const int n_pval_bytes = [&] { if (status == row_status::ok) - return std::sprintf(cursor, "%.6g", p_val); + return std::snprintf(cursor, bufsize, "%.6g", p_val); if (!more_na_info || status == row_status::na) - return std::sprintf(cursor, "NA"); + return std::snprintf(cursor, bufsize, "NA"); if (status == row_status::na_extreme_cnt) - return std::sprintf(cursor, "NA_EXTREME_CNT"); + return std::snprintf(cursor, bufsize, "NA_EXTREME_CNT"); // if (status == row_status::na_low_cov) - return std::sprintf(cursor, "NA_LOW_COV"); + return std::snprintf(cursor, bufsize, "NA_LOW_COV"); }(); if (n_pval_bytes < 0) return n_pval_bytes; + bufsize -= n_pval_bytes; cursor += n_pval_bytes; const int n_param_bytes = [&] { - std::int32_t n_bytes = 0; + std::int32_t n_param_bytes = 0; if (p_estim_alt.empty()) p_estim_alt.resize(n_groups, 0.0); for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) { - const int n = std::sprintf(cursor, "\t%f", p_estim_alt[g_idx]); + const int n = + std::snprintf(cursor, bufsize, "\t%f", p_estim_alt[g_idx]); + bufsize -= n; cursor += n; if (n < 0) return n; - n_bytes += n; + n_param_bytes += n; } const auto od = overdispersion_factor(n_samples, phi_estim_alt); - const int n = std::sprintf(cursor, "\t%f\n", od); + const int n = std::snprintf(cursor, bufsize, "\t%f\n", od); if (n < 0) return n; - n_bytes += n; - return n_bytes; + n_param_bytes += n; + return n_param_bytes; }(); if (n_param_bytes < 0) @@ -432,19 +262,6 @@ that the design matrix and the proportion table are correctly formatted. progress(table_file); } -static void -ensure_sample_order(const std::string &table_filename, Regression &alt_model, - Regression &null_model) { - std::ifstream table_file(table_filename); - if (!table_file) - throw std::runtime_error("could not open file: " + table_filename); - std::string header; - std::getline(table_file, header); - const auto sample_names = get_sample_names_from_header(header); - alt_model.order_samples(sample_names); - null_model.order_samples(sample_names); -} - int main_radmeth(int argc, char *argv[]) { try { @@ -476,76 +293,71 @@ main_radmeth(int argc, char *argv[]) { opt_parse.add_opt("factor", 'f', "a factor to test", true, test_factor); opt_parse.add_opt("tolerance", '\0', "numerical tolerance to test convergence", false, - Regression::tolerance); + radmeth_optimize_params::tolerance); opt_parse.add_opt("max-iter", '\0', "max iterations when estimating parameters", false, - Regression::max_iter); - + radmeth_optimize_params::max_iter); std::vector leftover_args; opt_parse.parse(argc, argv, leftover_args); if (argc == 1 || opt_parse.help_requested()) { - std::cerr << opt_parse.help_message() << std::endl - << opt_parse.about_message() << std::endl; + std::cerr << opt_parse.help_message() << '\n' + << opt_parse.about_message() << '\n'; return EXIT_SUCCESS; } if (opt_parse.about_requested()) { - std::cerr << opt_parse.about_message() << std::endl; + std::cerr << opt_parse.about_message() << '\n'; return EXIT_SUCCESS; } if (opt_parse.option_missing()) { - std::cerr << opt_parse.option_missing_message() << std::endl; + std::cerr << opt_parse.option_missing_message() << '\n'; return EXIT_SUCCESS; } if (leftover_args.size() != 2) { - std::cerr << opt_parse.help_message() << std::endl; + std::cerr << opt_parse.help_message() << '\n'; return EXIT_SUCCESS; } const std::string design_filename(leftover_args.front()); const std::string table_filename(leftover_args.back()); /****************** END COMMAND LINE OPTIONS *****************/ - if (verbose) - std::cerr << "design table filename: " << design_filename << "\n\n"; - - // initialize full design matrix from file - Regression alt_model; - alt_model.design = read_design(design_filename); - - if (verbose) - std::cerr << "Alternate model:\n" << alt_model.design << '\n'; + Design design = Design::read_design(design_filename); // Check that provided test factor name exists and find its index. - const auto test_factor_idx = get_test_factor_idx(alt_model, test_factor); + const auto test_factor_idx = design.get_test_factor_idx(test_factor); + + ensure_sample_order(table_filename, design); // verify that the design includes more than one level for the // test factor - if (!has_two_values(alt_model, test_factor_idx)) { - const auto first_level = alt_model.design.matrix[0][test_factor_idx]; + if (!design.has_two_values(test_factor_idx)) { + const auto first_level = design.matrix[0][test_factor_idx]; throw std::runtime_error("test factor only one level: " + test_factor + ", " + std::to_string(first_level)); } - Regression null_model = alt_model; - null_model.design = alt_model.design.drop_factor(test_factor_idx); - if (verbose) - std::cerr << "Null model:\n" << null_model.design << '\n'; - - ensure_sample_order(table_filename, alt_model, null_model); + const Design null_design = design.drop_factor(test_factor_idx); // clang-format off if (verbose) - std::cerr << "Output columns:\n" + std::cerr << "design table filename: " << design_filename << "\n\n" + << "Alternate model:\n" << design << '\n' + << "Null model:\n" << null_design << '\n' + << "Output columns:\n" << "(1) chrom\n" << "(2) position\n" << "(3) strand\n" << "(4) cytosine context\n" << "(5) p-value\n" - << "(6) estimated methylation 0\n" - << "(7) estimated methylation 1\n" - << "(8) overdispersion factor (variance above binomial)\n" - << "estimated methylation is for test factor value (0 or 1)\n" - << '\n'; + << "(6) overdispersion factor (variance above binomial)\n" + << "(7-) estimated methylation for each unique group\n" + << "groups are defined by combinations of factor levels\n\n"; // clang-format on + + Regression alt_model; + alt_model.design = design; + Regression null_model; + null_model.design = null_design; + const auto start_time = std::chrono::steady_clock::now(); radmeth(show_progress, more_na_info, n_threads, table_filename, outfile, alt_model, null_model, test_factor_idx); diff --git a/src/radmeth/radmeth_model.cpp b/src/radmeth/radmeth_design.cpp similarity index 71% rename from src/radmeth/radmeth_model.cpp rename to src/radmeth/radmeth_design.cpp index 1aec01fe..f08dd9ea 100644 --- a/src/radmeth/radmeth_model.cpp +++ b/src/radmeth/radmeth_design.cpp @@ -2,22 +2,23 @@ * * Authors: 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 "radmeth_model.hpp" +#include "radmeth_design.hpp" #include #include #include +#include #include #include #include @@ -25,60 +26,6 @@ #include #include -double Regression::tolerance = 1e-4; -double Regression::stepsize = 0.01; -std::uint32_t Regression::max_iter = 250; - -void -SiteProp::parse(const std::string &line) { - const auto first_ws = line.find_first_of(" \t"); - - // Parse the row name (must be like: "chr:position:strand:context") - bool failed = (first_ws == std::string::npos); - - auto field_s = line.data(); - auto field_e = line.data() + first_ws; - rowname = std::string{field_s, field_e}; - std::replace(std::begin(rowname), std::end(rowname), ':', '\t'); - if (failed) - throw std::runtime_error("failed to parse label from:\n" + line); - - // Parse the counts of total reads and methylated reads - const auto is_sep = [](const auto x) { return std::isspace(x); }; - const auto not_sep = [](const auto x) { return std::isdigit(x); }; - - const auto line_end = line.data() + std::size(line); - mcounts mc1{}; - mc.clear(); - while (field_e != line_end) { - // get the total count - field_s = std::find_if(field_e + 1, line_end, not_sep); - field_e = std::find_if(field_s + 1, line_end, is_sep); - { - const auto [ptr, ec] = std::from_chars(field_s, field_e, mc1.n_reads); - failed = failed || (ec != std::errc{}); - } - - // get the methylated count - field_s = std::find_if(field_e + 1, line_end, not_sep); - field_e = std::find_if(field_s + 1, line_end, is_sep); - { -#ifdef __APPLE__ - const int ret = std::sscanf(field_s, "%lf", &mc1.n_meth); - failed = failed || (ret < 1); -#else - const auto [ptr, ec] = std::from_chars(field_s, field_e, mc1.n_meth); - failed = failed || (ec != std::errc{}); -#endif - } - - mc.push_back(mc1); - } - - if (failed) - throw std::runtime_error("failed to parse counts from:\n" + line); -} - static void make_groups(Design &design) { auto s = design.matrix; @@ -252,3 +199,64 @@ operator<<(std::ostream &out, const Design &design) { return out; } + +void +ensure_sample_order(const std::string &table_filename, Design &design) { + std::ifstream table_file(table_filename); + if (!table_file) + throw std::runtime_error("could not open file: " + table_filename); + std::string header; + std::getline(table_file, header); + const auto sample_names = get_sample_names_from_header(header); + design.order_samples(sample_names); +} + +[[nodiscard]] std::vector +get_sample_names_from_header(const std::string &header) { + std::istringstream iss(header); + std::string token; + std::vector sample_names; + while (iss >> token) + sample_names.push_back(token); + return sample_names; +} + +[[nodiscard]] bool +consistent_sample_names(const Design &design, const std::string &header) { + std::istringstream iss(header); + auto nm_itr(std::begin(design.sample_names)); + const auto nm_end(std::end(design.sample_names)); + std::string token; + while (iss >> token && nm_itr != nm_end) + if (token != *nm_itr++) + return false; + return true; +} + +[[nodiscard]] Design +Design::read_design(const std::string &design_filename) { + std::ifstream design_file(design_filename); + if (!design_file) + throw std::runtime_error("could not open file: " + design_filename); + Design design; + design_file >> design; + return design; +} + +[[nodiscard]] std::uint32_t +Design::get_test_factor_idx(const std::string &test_factor) const { + const auto itr = + std::find(std::cbegin(factor_names), std::cend(factor_names), test_factor); + if (itr == std::cend(factor_names)) + throw std::runtime_error("factor not part of design: " + test_factor); + return std::distance(std::cbegin(factor_names), itr); +} + +[[nodiscard]] bool +Design::has_two_values(const std::size_t test_factor) const { + const auto &tcol = tmatrix[test_factor]; + for (const auto x : tcol) + if (x != tcol[0]) + return true; + return false; +} diff --git a/src/radmeth/radmeth_design.hpp b/src/radmeth/radmeth_design.hpp new file mode 100644 index 00000000..38d37cc6 --- /dev/null +++ b/src/radmeth/radmeth_design.hpp @@ -0,0 +1,79 @@ +/* Copyright (C) 2025 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 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. + */ + +#ifndef RADMETH_DESIGN_HPP +#define RADMETH_DESIGN_HPP + +#include +#include +#include +#include +#include + +struct Design { + std::vector factor_names; + std::vector sample_names; + std::vector> matrix; // samples=rows, factors=cols + std::vector> tmatrix; // factors=rows, samples=cols + std::vector> groups; // combs of fact levels + std::vector group_id; // assign group to sample + + [[nodiscard]] static Design + read_design(const std::string &design_filename); + + [[nodiscard]] std::size_t + n_factors() const { + return std::size(factor_names); + } + + [[nodiscard]] std::size_t + n_groups() const { + return std::size(groups); + } + + [[nodiscard]] std::size_t + n_samples() const { + return std::size(sample_names); + } + + [[nodiscard]] Design + drop_factor(const std::size_t factor_idx); + + void + order_samples(const std::vector &ordered_names); + + [[nodiscard]] std::uint32_t + get_test_factor_idx(const std::string &test_factor) const; + + [[nodiscard]] bool + has_two_values(const std::size_t test_factor) const; +}; + +std::istream & +operator>>(std::istream &is, Design &design); + +std::ostream & +operator<<(std::ostream &os, const Design &design); + +void +ensure_sample_order(const std::string &table_filename, Design &design); + +[[nodiscard]] std::vector +get_sample_names_from_header(const std::string &header); + +[[nodiscard]] bool +consistent_sample_names(const Design &design, const std::string &header); + +#endif // RADMETH_DESIGN_HPP diff --git a/src/radmeth/radmeth_model.hpp b/src/radmeth/radmeth_model.hpp index 013e51f4..9be0271a 100644 --- a/src/radmeth/radmeth_model.hpp +++ b/src/radmeth/radmeth_model.hpp @@ -1,82 +1,46 @@ -/* Copyright (C) 2013-2023 University of Southern California and - * Egor Dolzhenko - * Andrew D Smith - * Guilherme Sena +/* Copyright (C) 2025 Andrew D Smith * - * Authors: Andrew D. Smith and Egor Dolzhenko and Guilherme Sena + * 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. */ #ifndef RADMETH_MODEL_HPP #define RADMETH_MODEL_HPP +#include "radmeth_design.hpp" + +#include +#include #include #include #include #include #include -struct Design { - std::vector factor_names; - std::vector sample_names; - std::vector> matrix; // samples=rows, factors=cols - std::vector> tmatrix; // factors=rows, samples=cols - std::vector> groups; // combs of fact levels - std::vector group_id; // assign group to sample - - [[nodiscard]] std::size_t - n_factors() const { - return std::size(factor_names); - } - - [[nodiscard]] std::size_t - n_groups() const { - return std::size(groups); - } - - [[nodiscard]] std::size_t - n_samples() const { - return std::size(sample_names); - } - - [[nodiscard]] Design - drop_factor(const std::size_t factor_idx); - - void - order_samples(const std::vector &ordered_names); -}; - -std::istream & -operator>>(std::istream &is, Design &design); - -std::ostream & -operator<<(std::ostream &os, const Design &design); - -struct mcounts { +template struct mcounts { std::uint32_t n_reads{}; - double n_meth{}; + T n_meth{}; }; +template [[nodiscard]] inline std::istream & -operator>>(std::istream &in, mcounts &rm) { +operator>>(std::istream &in, mcounts &rm) { return in >> rm.n_reads >> rm.n_meth; } -struct SiteProp { - std::string rowname; - std::vector mc; - - void - parse(const std::string &line); +struct cumul_counts { + std::vector m_counts; + std::vector r_counts; + std::vector d_counts; }; struct vars_cache { @@ -91,16 +55,26 @@ struct vars_cache { double digamma_a_b{}; }; -struct Regression { +template struct Regression { static double tolerance; // 1e-3; static double stepsize; // 0.001; static std::uint32_t max_iter; // 250; Design design; - SiteProp props; + std::string rowname; + std::vector> mc; double max_loglik{}; + // scratch space + std::vector cumul; + std::vector p_v; + std::vector cache_log1p_factors; + std::vector cache_dispersion_effect; std::vector cache; // scratch space + std::uint32_t max_r_count{}; + + void + parse(const std::string &line); [[nodiscard]] std::size_t n_factors() const { @@ -119,23 +93,64 @@ struct Regression { [[nodiscard]] std::size_t props_size() const { - return std::size(props.mc); + return std::size(mc); } [[nodiscard]] std::size_t n_samples() const { return design.n_samples(); } +}; - [[nodiscard]] std::string - rowname() const { - return props.rowname; - } +template +void +Regression::parse(const std::string &line) { + const auto first_ws = line.find_first_of(" \t"); + + // Parse the row name (must be like: "chr:position:strand:context") + bool failed = (first_ws == std::string::npos); + + auto field_s = line.data(); + auto field_e = line.data() + first_ws; + rowname = std::string{field_s, field_e}; + std::replace(std::begin(rowname), std::end(rowname), ':', '\t'); + if (failed) + throw std::runtime_error("failed to parse label from:\n" + line); + + // Parse the counts of total reads and methylated reads + const auto is_sep = [](const auto x) { return std::isspace(x); }; + const auto not_sep = [](const auto x) { return std::isdigit(x); }; + + const auto line_end = line.data() + std::size(line); + mcounts mc1{}; + mc.clear(); + while (field_e != line_end) { + // get the total count + field_s = std::find_if(field_e + 1, line_end, not_sep); + field_e = std::find_if(field_s + 1, line_end, is_sep); + { + const auto [ptr, ec] = std::from_chars(field_s, field_e, mc1.n_reads); + failed = failed || (ptr != field_e); + } + + field_s = std::find_if(field_e + 1, line_end, not_sep); + field_e = std::find_if(field_s + 1, line_end, is_sep); + + { +#ifdef __APPLE__ + const int ret = std::sscanf(field_s, "%lf", &mc1.n_meth); + failed = failed || (ret < 1); +#else + const auto [ptr, ec] = std::from_chars(field_s, field_e, mc1.n_meth); + failed = failed || (ec != std::errc{}); +#endif + } - void - order_samples(const std::vector &ordered_names) { - design.order_samples(ordered_names); + mc.push_back(mc1); } -}; -#endif + if (failed) + throw std::runtime_error("failed to parse counts from:\n" + line); +} + +#endif // RADMETH_MODEL_HPP diff --git a/src/radmeth/radmeth_nano.cpp b/src/radmeth/radmeth_nano.cpp new file mode 100644 index 00000000..0932cb84 --- /dev/null +++ b/src/radmeth/radmeth_nano.cpp @@ -0,0 +1,374 @@ +/* Copyright (C) 2025 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 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 "radmeth_design.hpp" +#include "radmeth_model.hpp" +#include "radmeth_optimize_gamma.hpp" +#include "radmeth_optimize_params.hpp" +#include "radmeth_utils.hpp" + +// smithlab headers +#include "GenomicRegion.hpp" +#include "OptionParser.hpp" +#include "smithlab_os.hpp" +#include "smithlab_utils.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +[[nodiscard]] static bool +has_low_coverage(const Regression ®, const std::size_t test_factor) { + bool cvrd_in_test_fact_smpls = false; + const auto &tcol = reg.design.tmatrix[test_factor]; + for (std::size_t i = 0; i < reg.n_samples() && !cvrd_in_test_fact_smpls; ++i) + cvrd_in_test_fact_smpls = (tcol[i] == 1 && reg.mc[i].n_reads != 0); + + bool cvrd_in_other_smpls = false; + for (std::size_t i = 0; i < reg.n_samples() && !cvrd_in_other_smpls; ++i) + cvrd_in_other_smpls = (tcol[i] != 1 && reg.mc[i].n_reads != 0); + + return !cvrd_in_test_fact_smpls || !cvrd_in_other_smpls; +} + +[[nodiscard]] static bool +has_extreme_counts(const Regression ®) { + const auto &mc = reg.mc; + + bool full_meth = true; + for (auto i = std::cbegin(mc); i != std::cend(mc) && full_meth; ++i) + full_meth = (i->n_reads == i->n_meth); + + bool no_meth = true; + for (auto i = std::cbegin(mc); i != std::cend(mc) && no_meth; ++i) + no_meth = (i->n_meth == 0.0); + + return full_meth || no_meth; +} + +enum class row_status : std::uint8_t { + ok, + na, + na_low_cov, + na_extreme_cnt, +}; + +template +static void +radmeth(const bool show_progress, const bool more_na_info, + const std::uint32_t n_threads, const std::string &table_filename, + const std::string &outfile, const Regression &alt_model, + const Regression &null_model, const std::uint32_t test_factor_idx) { + static constexpr auto buf_size = 1024; + static constexpr auto max_lines = 16384; + + // ADS: open the data table file + std::ifstream table_file(table_filename); + if (!table_file) + throw std::runtime_error("could not open file: " + table_filename); + + // Make sure that the first line of the proportion table file contains + // names of the samples. Throw an exception if the names or their order + // in the proportion table does not match those in the full design matrix. + std::string sample_names_header; + std::getline(table_file, sample_names_header); + + const auto sample_names = get_sample_names_from_header(sample_names_header); + if (!consistent_sample_names(alt_model.design, sample_names_header)) { + // clang-format off + const auto message = + R"( +does not match factor names or their order in the design matrix. Check +that the design matrix and the proportion table are correctly formatted. +)"; + // clang-format on + throw std::runtime_error("header:\n" + sample_names_header + message); + } + + const std::size_t n_samples = alt_model.n_samples(); + + std::ofstream out(outfile); + if (!out) + throw std::runtime_error("failed to open output file: " + outfile); + + file_progress progress{table_filename}; + + std::vector> bufs(max_lines, + std::vector(buf_size, 0)); + std::vector n_bytes(max_lines, 0); + std::vector lines(max_lines); + + std::vector> alt_models(n_threads, alt_model); + std::vector> null_models(n_threads, null_model); + + const auto n_groups = alt_model.n_groups(); + + // Iterate over rows in the file and do llr test on proportions from + // each. Do this in sets of rows to avoid having to spawn too many threads. + while (true) { + std::uint32_t n_lines = 0; + while (n_lines < max_lines && std::getline(table_file, lines[n_lines])) + ++n_lines; + if (n_lines == 0) + break; + + std::vector threads; + for (auto thread_id = 0u; thread_id < n_threads; ++thread_id) { + threads.emplace_back([&, thread_id] { + std::vector p_estim_alt; + std::vector p_estim_null; + double phi_estim_alt{}; + double phi_estim_null{}; + + auto &t_alt_model = alt_models[thread_id]; + auto &t_null_model = null_models[thread_id]; + for (auto b = 0u; b < n_lines; ++b) { + // ADS: rows done by different threads are interleaved because the + // difficult (e.g., high-coverage) rows can be consecutive and this + // balances work better. + if (b % n_threads != thread_id) + continue; + t_alt_model.parse(lines[b]); + if (t_alt_model.props_size() != n_samples) + throw std::runtime_error("found row with wrong number of columns"); + + const auto p_val_status = [&]() -> std::tuple { + // Skip the test if (1) no coverage in all cases or in all + // controls, or (2) the site is completely methylated or + // completely unmethylated across all samples. + if (has_low_coverage(t_alt_model, test_factor_idx)) + return std::tuple{1.0, row_status::na_low_cov}; + + if (has_extreme_counts(t_alt_model)) + return std::tuple{1.0, row_status::na_extreme_cnt}; + + fit_regression_model_gamma_fdf(t_alt_model, p_estim_alt, + phi_estim_alt); + + t_null_model.mc = t_alt_model.mc; + t_null_model.rowname = t_alt_model.rowname; + + fit_regression_model_gamma_fdf(t_null_model, p_estim_null, + phi_estim_null); + + const double p_value = + llr_test(t_null_model.max_loglik, t_alt_model.max_loglik); + + return (p_value != p_value) ? std::tuple{1.0, row_status::na} + : std::tuple{p_value, row_status::ok}; + }(); + // ADS: avoid capture structured binding in C++17 + const auto p_val = std::get<0>(p_val_status); + const auto status = std::get<1>(p_val_status); + + n_bytes[b] = [&] { + auto bufsize = std::size(bufs[b]); + // clang-format off + const int n_prefix_bytes = + std::snprintf(bufs[b].data(), bufsize, + "%s\t", t_alt_model.rowname.data()); + // clang-format on + if (n_prefix_bytes < 0) + return n_prefix_bytes; + + bufsize -= n_prefix_bytes; + auto cursor = bufs[b].data() + n_prefix_bytes; + + const int n_pval_bytes = [&] { + if (status == row_status::ok) + return std::snprintf(cursor, bufsize, "%.6g", p_val); + if (!more_na_info || status == row_status::na) + return std::snprintf(cursor, bufsize, "NA"); + if (status == row_status::na_extreme_cnt) + return std::snprintf(cursor, bufsize, "NA_EXTREME_CNT"); + // if (status == row_status::na_low_cov) + return std::snprintf(cursor, bufsize, "NA_LOW_COV"); + }(); + + if (n_pval_bytes < 0) + return n_pval_bytes; + + bufsize -= n_pval_bytes; + cursor += n_pval_bytes; + + const int n_param_bytes = [&] { + std::int32_t n_param_bytes = 0; + if (p_estim_alt.empty()) + p_estim_alt.resize(n_groups, 0.0); + for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) { + const int n = + std::snprintf(cursor, bufsize, "\t%f", p_estim_alt[g_idx]); + bufsize -= n; + cursor += n; + if (n < 0) + return n; + n_param_bytes += n; + } + const auto od = overdispersion_factor(n_samples, phi_estim_alt); + const int n = std::snprintf(cursor, bufsize, "\t%f\n", od); + if (n < 0) + return n; + n_param_bytes += n; + return n_param_bytes; + }(); + + if (n_param_bytes < 0) + return n_param_bytes; + + return n_prefix_bytes + n_pval_bytes + n_param_bytes; + }(); + } + }); + } + + for (auto &thread : threads) + thread.join(); + + for (auto i = 0u; i < n_lines; ++i) { + if (n_bytes[i] < 0) + throw std::runtime_error("failed to write to output buffer"); + out.write(bufs[i].data(), n_bytes[i]); + } + + if (show_progress) + progress(table_file); + + if (n_lines < n_threads) + break; + } + + if (show_progress) + progress(table_file); +} + +int +main_radmeth_nano(int argc, char *argv[]) { + try { + static const std::string description = + "calculate differential methylation scores for nanopore data"; + + std::string outfile; + std::string test_factor; + bool verbose{false}; + bool show_progress{false}; + bool more_na_info{false}; + std::uint32_t n_threads{1}; + + /****************** COMMAND LINE OPTIONS ********************/ + OptionParser opt_parse(strip_path(argv[0]), description, + " "); + opt_parse.set_show_defaults(); + opt_parse.add_opt("out", 'o', "output file", true, outfile); + opt_parse.add_opt("threads", 't', "number of threads to use", false, + n_threads); + opt_parse.add_opt("verbose", 'v', "print more run info", false, verbose); + opt_parse.add_opt("progress", '\0', "show progress", false, show_progress); + opt_parse.add_opt( + "na-info", 'n', + "if a p-value is not calculated, print NAs in more detail: " + "low count (NA_LOW_COV) extreme values (NA_EXTREME_CNT) or " + "numerical errors in likelihood ratios (NA)", + false, more_na_info); + opt_parse.add_opt("factor", 'f', "a factor to test", true, test_factor); + opt_parse.add_opt("tolerance", '\0', + "numerical tolerance to test convergence", false, + radmeth_optimize_params::tolerance); + opt_parse.add_opt("max-iter", '\0', + "max iterations when estimating parameters", false, + radmeth_optimize_params::max_iter); + std::vector leftover_args; + opt_parse.parse(argc, argv, leftover_args); + if (argc == 1 || opt_parse.help_requested()) { + std::cerr << opt_parse.help_message() << '\n' + << opt_parse.about_message() << '\n'; + return EXIT_SUCCESS; + } + if (opt_parse.about_requested()) { + std::cerr << opt_parse.about_message() << '\n'; + return EXIT_SUCCESS; + } + if (opt_parse.option_missing()) { + std::cerr << opt_parse.option_missing_message() << '\n'; + return EXIT_SUCCESS; + } + if (leftover_args.size() != 2) { + std::cerr << opt_parse.help_message() << '\n'; + return EXIT_SUCCESS; + } + const std::string design_filename(leftover_args.front()); + const std::string table_filename(leftover_args.back()); + /****************** END COMMAND LINE OPTIONS *****************/ + + Design design = Design::read_design(design_filename); + + // Check that provided test factor name exists and find its index. + const auto test_factor_idx = design.get_test_factor_idx(test_factor); + + ensure_sample_order(table_filename, design); + + // verify that the design includes more than one level for the + // test factor + if (!design.has_two_values(test_factor_idx)) { + const auto first_level = design.matrix[0][test_factor_idx]; + throw std::runtime_error("test factor only one level: " + test_factor + + ", " + std::to_string(first_level)); + } + + const Design null_design = design.drop_factor(test_factor_idx); + + // clang-format off + if (verbose) + std::cerr << "design table filename: " << design_filename << "\n\n" + << "Alternate model:\n" << design << '\n' + << "Null model:\n" << null_design << '\n' + << "Output columns:\n" + << "(1) chrom\n" + << "(2) position\n" + << "(3) strand\n" + << "(4) cytosine context\n" + << "(5) p-value\n" + << "(6) overdispersion factor (variance above binomial)\n" + << "(7-) estimated methylation for each unique group\n" + << "groups are defined by combinations of factor levels\n\n"; + // clang-format on + + Regression alt_model; + alt_model.design = design; + Regression null_model; + null_model.design = null_design; + + const auto start_time = std::chrono::steady_clock::now(); + radmeth(show_progress, more_na_info, n_threads, table_filename, outfile, + alt_model, null_model, test_factor_idx); + const auto stop_time = std::chrono::steady_clock::now(); + + if (verbose) + std::cerr << "[total time: " << format_duration(stop_time - start_time) + << "]\n"; + } + catch (const std::exception &e) { + std::cerr << e.what() << '\n'; + return EXIT_FAILURE; + } + return EXIT_SUCCESS; +} diff --git a/src/radmeth/radmeth_optimize.cpp b/src/radmeth/radmeth_optimize.cpp deleted file mode 100644 index 6c6e2621..00000000 --- a/src/radmeth/radmeth_optimize.cpp +++ /dev/null @@ -1,250 +0,0 @@ -/* Copyright (C) 2013-2025 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 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 "radmeth_optimize.hpp" -#include "radmeth_model.hpp" - -#include -#include // for gsl_sf_psi (digamma) -#include - -#include -#include -#include -#include - -[[nodiscard]] static inline double -logistic(const double x) { - return 1.0 / (1.0 / std::exp(x) + 1.0); -} - -template -[[nodiscard]] static double -get_p(const std::vector &v, const gsl_vector *params) { - const auto a = v.data(); - return logistic(std::inner_product(a, a + std::size(v), params->data, 0.0)); -} - -static void -get_cache_lgamma(const std::vector &group, - const gsl_vector *params, const double phi, - vars_cache &cache) { - const auto p = get_p(group, params); - const auto a = p * phi; - const auto b = (1.0 - p) * phi; - cache.p = p; - cache.a = a; - cache.b = b; - cache.lgamma_a = std::lgamma(a); - cache.lgamma_b = std::lgamma(b); - cache.lgamma_a_b = std::lgamma(a + b); -} - -[[nodiscard]] static double -log_likelihood(const gsl_vector *params, Regression ®) { - const auto phi = 1.0 / std::exp(gsl_vector_get(params, reg.n_factors())); - - const auto n_groups = reg.n_groups(); - const auto &groups = reg.design.groups; - auto &cache = reg.cache; - for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) - get_cache_lgamma(groups[g_idx], params, phi, cache[g_idx]); - - const auto &group_id = reg.design.group_id; - const auto &mc = reg.props.mc; - double ll = 0.0; - const auto n_samples = reg.n_samples(); - for (auto i = 0u; i < n_samples; ++i) { - const auto y = mc[i].n_meth; - const auto n = mc[i].n_reads; - - const auto &c = cache[group_id[i]]; - const auto a = c.a; - const auto b = c.b; - - // clang-format off - ll += ((std::lgamma(y + a) - c.lgamma_a) + - (std::lgamma(n - y + b) - c.lgamma_b) + - (c.lgamma_a_b - std::lgamma(n + a + b))); - // clang-format on - } - - return ll; -} - -[[nodiscard]] static inline double -digamma(const double x) { - return gsl_sf_psi(x); -} - -static void -get_cache_digamma(const std::vector &group, - const gsl_vector *params, const double phi, - vars_cache &cache) { - const auto p = get_p(group, params); - const auto a = p * phi; - const auto b = (1.0 - p) * phi; - cache.p = p; - cache.a = a; - cache.b = b; - cache.digamma_a = digamma(a); - cache.digamma_b = digamma(b); - cache.digamma_a_b = digamma(a + b); -} - -static void -gradient(const gsl_vector *params, Regression ®, gsl_vector *output) { - const auto n_factors = reg.n_factors(); - const auto phi = 1.0 / std::exp(gsl_vector_get(params, n_factors)); - - const auto &mc = reg.props.mc; - const auto &matrix = reg.design.matrix; - - const auto n_groups = reg.n_groups(); - const auto &groups = reg.design.groups; - auto &cache = reg.cache; // ADS: reusing scratch space - for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) - get_cache_digamma(groups[g_idx], params, phi, cache[g_idx]); - - gsl_vector_set_all(output, 0.0); // init output to zero for all factors - auto &data = output->data; - - const auto &group_id = reg.design.group_id; - double grad_phi = 0.0; - const auto n_samples = reg.n_samples(); - for (auto i = 0u; i < n_samples; ++i) { - const auto y = mc[i].n_meth; - const auto n = mc[i].n_reads; - - const auto &c = cache[group_id[i]]; - const auto p = c.p; - const auto a = c.a; - const auto b = c.b; - - const auto digamma_a = c.digamma_a; - const auto digamma_b = c.digamma_b; - const auto digamma_y_a = digamma(y + a); - const auto digamma_n_y_b = digamma(n - y + b); - - // grad wrt p - const auto dlogl_dp = - phi * (digamma_y_a - digamma_a - digamma_n_y_b + digamma_b); - const auto dp_delta = p * (1.0 - p); // chain rule: d p / d param - const auto dlogl_delta = dlogl_dp * dp_delta; // grad wrt params - - auto matrix_itr = std::cbegin(matrix[i]); - const auto data_end = data + n_factors; - for (auto data_itr = data; data_itr != data_end; ++data_itr) - *data_itr += (*matrix_itr++) * dlogl_delta; - - const auto digamma_delta = c.digamma_a_b - digamma(n + a + b); - const auto dphi_term1 = p * (digamma_y_a - digamma_a + digamma_delta); - const auto dphi_term2 = - (1.0 - p) * (digamma_n_y_b - digamma_b + digamma_delta); - - grad_phi += dphi_term1 + dphi_term2; - } - - const auto grad_theta = -grad_phi * phi; - gsl_vector_set(output, n_factors, grad_theta); -} - -[[nodiscard]] static double -neg_loglik(const gsl_vector *params, void *object) { - auto reg = static_cast(object); - return -log_likelihood(params, *reg); -} - -static void -neg_gradient(const gsl_vector *params, void *object, gsl_vector *output) { - auto reg = static_cast(object); - gradient(params, *reg, output); - gsl_vector_scale(output, -1.0); -} - -static void -neg_loglik_and_grad(const gsl_vector *params, void *object, double *loglik_val, - gsl_vector *d_loglik_val) { - *loglik_val = neg_loglik(params, object); - neg_gradient(params, object, d_loglik_val); -} - -void -fit_regression_model(Regression &r, std::vector &p_estimates, - double &dispersion_estimate) { - static constexpr auto init_dispersion_param = -2.5; - const auto stepsize = Regression::stepsize; - const auto max_iter = Regression::max_iter; - const auto n_groups = r.n_groups(); - r.cache.resize(n_groups); // make sure scratch space is allocated - - const std::size_t n_params = r.n_params(); - const auto tol = std::sqrt(n_params) * r.n_samples() * Regression::tolerance; - // clang-format off - auto loglik_bundle = gsl_multimin_function_fdf{ - &neg_loglik, // objective function - &neg_gradient, // gradient - &neg_loglik_and_grad, // combined obj and grad - n_params, // number of model params - static_cast(&r) // parameters for objective and gradient functions - }; - // clang-format on - - // set the parameters: zero for "p" parameters and the final one for - // dispersion using the constant - auto params = gsl_vector_alloc(n_params); - gsl_vector_set_all(params, 0.0); - gsl_vector_set(params, n_params - 1, init_dispersion_param); - - // Alternatives: - // - gsl_multimin_fdfminimizer_conjugate_pr - // - gsl_multimin_fdfminimizer_conjugate_fr - // - gsl_multimin_fdfminimizer_vector_bfgs2 - // - gsl_multimin_fdfminimizer_steepest_descent - const auto minimizer = gsl_multimin_fdfminimizer_conjugate_pr; - auto s = gsl_multimin_fdfminimizer_alloc(minimizer, n_params); - - gsl_multimin_fdfminimizer_set(s, &loglik_bundle, params, stepsize, tol); - - int status = 0; - std::size_t iter = 0; - do { - status = gsl_multimin_fdfminimizer_iterate(s); // one iter and get status - if (status) - break; - // check status from gradient - status = gsl_multimin_test_gradient(s->gradient, tol); - } while (status == GSL_CONTINUE && ++iter < max_iter); - - /// ADS: the condition below might not always pass even if we are doing - /// ok. It's not clear how to check for failure. - - // if (status != GSL_SUCCESS) - // throw std::runtime_error("failed to fit model parameters"); - - const auto param_estimates = gsl_multimin_fdfminimizer_x(s); - - const auto &groups = r.design.groups; - p_estimates.clear(); - for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) - p_estimates.push_back(get_p(groups[g_idx], param_estimates)); - const auto disp_param = gsl_vector_get(param_estimates, n_params - 1); - dispersion_estimate = 1.0 / std::exp(disp_param); - - r.max_loglik = log_likelihood(param_estimates, r); - - gsl_multimin_fdfminimizer_free(s); - gsl_vector_free(params); -} diff --git a/src/radmeth/radmeth_optimize.hpp b/src/radmeth/radmeth_optimize.hpp deleted file mode 100644 index d0209fb2..00000000 --- a/src/radmeth/radmeth_optimize.hpp +++ /dev/null @@ -1,30 +0,0 @@ -/* Copyright (C) 2013-2023 University of Southern California and - * Egor Dolzhenko - * Andrew D Smith - * Guilherme Sena - * - * Authors: Andrew D. Smith and Egor Dolzhenko and Guilherme Sena - * - * 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. - */ - -#ifndef RADMETH_OPTIMIZE_HPP -#define RADMETH_OPTIMIZE_HPP - -#include - -struct Regression; - -void -fit_regression_model(Regression &r, std::vector &p_estimates, - double &dispersion_estimate); - -#endif diff --git a/src/radmeth/radmeth_optimize_gamma.cpp b/src/radmeth/radmeth_optimize_gamma.cpp new file mode 100644 index 00000000..d4efc75b --- /dev/null +++ b/src/radmeth/radmeth_optimize_gamma.cpp @@ -0,0 +1,420 @@ +/* Copyright (C) 2013-2025 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 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 "radmeth_optimize_gamma.hpp" +#include "radmeth_model.hpp" +#include "radmeth_optimize_params.hpp" + +#include +#include + +#include +#include +#include +#include +#include + +/* Coefficients for the Chebyschev polynomial for the digamma function in the + range 0-1 */ +// clang-format off +static constexpr std::array psi1_cs { + -0.038057080835217922, // == -0.019028540417608961*2 + 0.491415393029387130, + -0.056815747821244730, + 0.008357821225914313, + -0.001333232857994342, + 0.000220313287069308, + -0.000037040238178456, + 0.000006283793654854, + -0.000001071263908506, + 0.000000183128394654, + -0.000000031353509361, + 0.000000005372808776, + -0.000000000921168141, + 0.000000000157981265, + -0.000000000027098646, + 0.000000000004648722, + -0.000000000000797527, + 0.000000000000136827, + -0.000000000000023475, + 0.000000000000004027, + -0.000000000000000691, + 0.000000000000000118, + -0.000000000000000020 +}; +// clang-format on + +/* Alternate set of coefficients for the Chebyschev polynomial for the digamma + function */ +// clang-format off +static constexpr std::array psi2_cs = { + -0.0204749044678185, // == -0.01023745223390925*2 + -0.0101801271534859, + 0.0000559718725387, + -0.0000012917176570, + 0.0000000572858606, + -0.0000000038213539, + 0.0000000003397434, + -0.0000000000374838, + 0.0000000000048990, + -0.0000000000007344, + 0.0000000000001233, + -0.0000000000000228, + 0.0000000000000045, + -0.0000000000000009, + 0.0000000000000002, + -0.0000000000000000 +}; +// clang-format on + +template +[[nodiscard]] static inline double +chebyschev(const std::array &coeffs, std::uint32_t order, + const double y) { + const auto y2 = 2.0 * y; + double d = 0.0; + double dd = 0.0; + for (auto j = order; j >= 1; j--) { + const auto temp = d; + d = y2 * d - dd + coeffs[j]; + dd = temp; + } + return y * d - dd + 0.5 * coeffs[0]; +} + +// digamma for x non-negative +static inline auto +digamma(const double y) -> double { + static constexpr auto psi1_order = 22; // max=22; + static constexpr auto psi2_order = 15; // max=15; + if (y >= 2.0) { + const auto t = 8.0 / (y * y) - 1.0; + return std::log(y) - 0.5 / y + chebyschev(psi2_cs, psi2_order, t); + } + if (y < 1.0) + return -1.0 / y + chebyschev(psi1_cs, psi1_order, 2.0 * y - 1.0); + return chebyschev(psi1_cs, psi1_order, 2.0 * (y - 1.0) - 1.0); +} + +[[nodiscard]] static inline auto +logistic(const double x) -> double { + return 1.0 / (1.0 / std::exp(x) + 1.0); +} + +template +[[nodiscard]] static auto +get_p(const std::vector &v, const gsl_vector *params) -> double { + const auto a = v.data(); + return logistic(std::inner_product(a, a + std::size(v), params->data, 0.0)); +} + +static auto +get_cache_lgamma(const std::vector &group, + const gsl_vector *params, const double phi, + vars_cache &cache) { + const auto p = get_p(group, params); + const auto a = p * phi; + const auto b = (1.0 - p) * phi; + cache.p = p; + cache.a = a; + cache.b = b; + cache.lgamma_a = std::lgamma(a); + cache.lgamma_b = std::lgamma(b); + cache.lgamma_a_b = std::lgamma(a + b); +} + +template +[[nodiscard]] static auto +log_likelihood(const gsl_vector *params, Regression ®) -> double { + const auto phi = 1.0 / std::exp(gsl_vector_get(params, reg.n_factors())); + + const auto n_groups = reg.n_groups(); + const auto &groups = reg.design.groups; + auto &cache = reg.cache; + for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) + get_cache_lgamma(groups[g_idx], params, phi, cache[g_idx]); + + const auto &group_id = reg.design.group_id; + const auto &mc = reg.mc; + double ll = 0.0; + const auto n_samples = reg.n_samples(); + for (auto i = 0u; i < n_samples; ++i) { + const auto y = mc[i].n_meth; + const auto n = mc[i].n_reads; + + const auto &c = cache[group_id[i]]; + const auto a = c.a; + const auto b = c.b; + + // clang-format off + ll += ((std::lgamma(y + a) - c.lgamma_a) + + (std::lgamma(n - y + b) - c.lgamma_b) + + (c.lgamma_a_b - std::lgamma(n + a + b))); + // clang-format on + } + + return ll; +} + +static void +get_cache_digamma(const std::vector &group, + const gsl_vector *params, const double phi, + vars_cache &cache) { + const auto p = get_p(group, params); + const auto a = p * phi; + const auto b = (1.0 - p) * phi; + cache.p = p; + cache.a = a; + cache.b = b; + cache.digamma_a = digamma(a); + cache.digamma_b = digamma(b); + cache.digamma_a_b = digamma(a + b); +} + +template +static void +gradient(const gsl_vector *params, Regression ®, gsl_vector *output) { + const auto n_factors = reg.n_factors(); + const auto phi = 1.0 / std::exp(gsl_vector_get(params, n_factors)); + + const auto &mc = reg.mc; + const auto &matrix = reg.design.matrix; + + const auto n_groups = reg.n_groups(); + const auto &groups = reg.design.groups; + auto &cache = reg.cache; // ADS: reusing scratch space + for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) + get_cache_digamma(groups[g_idx], params, phi, cache[g_idx]); + + gsl_vector_set_all(output, 0.0); // init output to zero for all factors + const auto &data = output->data; + + const auto &group_id = reg.design.group_id; + double grad_phi = 0.0; + const auto n_samples = reg.n_samples(); + for (auto i = 0u; i < n_samples; ++i) { + const double y = mc[i].n_meth; // ADS: n_meth maybe int or float + const auto n = mc[i].n_reads; + + const auto &c = cache[group_id[i]]; + const auto p = c.p; + const auto a = c.a; + const auto b = c.b; + + const auto digamma_a = c.digamma_a; + const auto digamma_b = c.digamma_b; + const auto digamma_y_a = digamma(y + a); + const auto digamma_n_y_b = digamma(n - y + b); + + // grad wrt p + const auto dlogl_dp = + phi * (digamma_y_a - digamma_a - digamma_n_y_b + digamma_b); + const auto dp_delta = p * (1.0 - p); // chain rule: d p / d param + const auto dlogl_delta = dlogl_dp * dp_delta; // grad wrt params + + auto matrix_itr = std::cbegin(matrix[i]); + const auto data_end = data + n_factors; + for (auto data_itr = data; data_itr != data_end; ++data_itr) + *data_itr += (*matrix_itr++) * dlogl_delta; + + const auto digamma_delta = c.digamma_a_b - digamma(n + a + b); + const auto dphi_term1 = p * (digamma_y_a - digamma_a + digamma_delta); + const auto dphi_term2 = + (1.0 - p) * (digamma_n_y_b - digamma_b + digamma_delta); + + grad_phi += dphi_term1 + dphi_term2; + } + + const auto grad_theta = -grad_phi * phi; + gsl_vector_set(output, n_factors, grad_theta); +} + +template +[[nodiscard]] static double +neg_loglik(const gsl_vector *params, void *object) { + auto reg = static_cast *>(object); + return -log_likelihood(params, *reg); +} + +template +static void +neg_gradient(const gsl_vector *params, void *object, gsl_vector *output) { + auto reg = static_cast *>(object); + gradient(params, *reg, output); + gsl_vector_scale(output, -1.0); +} + +template +static void +neg_loglik_and_grad(const gsl_vector *params, void *object, double *loglik_val, + gsl_vector *d_loglik_val) { + *loglik_val = neg_loglik(params, object); + neg_gradient(params, object, d_loglik_val); +} + +template +static void +fit_regression_model(Regression &r, std::vector &p_estimates, + double &dispersion_estimate) { + static constexpr auto init_dispersion_param = -2.5; + const auto stepsize = radmeth_optimize_params::stepsize; + const auto max_iter = radmeth_optimize_params::max_iter; + const auto n_groups = r.n_groups(); + r.cache.resize(n_groups); // make sure scratch space is allocated + + const std::size_t n_params = r.n_params(); + const auto tol = radmeth_optimize_params::tolerance; + + // set the parameters: zero for "p" parameters and the final one for + // dispersion using the constant + auto params = gsl_vector_alloc(n_params); + gsl_vector_set_all(params, 0.0); + gsl_vector_set(params, n_params - 1, init_dispersion_param); + + auto minex_func = gsl_multimin_function{ + &neg_loglik, + n_params, + static_cast(&r), + }; + + // set initial step size for all dims + auto step_sizes = gsl_vector_alloc(n_params); + gsl_vector_set_all(step_sizes, stepsize); + + const auto minimizer = gsl_multimin_fminimizer_nmsimplex2; + auto s = gsl_multimin_fminimizer_alloc(minimizer, n_params); + gsl_multimin_fminimizer_set(s, &minex_func, params, step_sizes); + + int status = 0; + std::size_t iter = 0; + + do { + status = gsl_multimin_fminimizer_iterate(s); // one iter and get status + if (status) + break; + + const auto size = gsl_multimin_fminimizer_size(s); + status = gsl_multimin_test_size(size, tol); + } while (status == GSL_CONTINUE && ++iter < max_iter); + // ADS: can't use (status != GSL_SUCCESS) + + const auto param_estimates = gsl_multimin_fminimizer_x(s); + + const auto &groups = r.design.groups; + p_estimates.clear(); + for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) + p_estimates.push_back(get_p(groups[g_idx], param_estimates)); + const auto disp_param = gsl_vector_get(param_estimates, n_params - 1); + dispersion_estimate = 1.0 / std::exp(disp_param); + + r.max_loglik = log_likelihood(param_estimates, r); + + gsl_multimin_fminimizer_free(s); + gsl_vector_free(step_sizes); + gsl_vector_free(params); +} + +template +static void +fit_regression_model_fdf(Regression &r, std::vector &p_estimates, + double &dispersion_estimate) { + static constexpr auto init_dispersion_param = -2.5; + const auto stepsize = radmeth_optimize_params::stepsize; + const auto max_iter = radmeth_optimize_params::max_iter; + const auto n_groups = r.n_groups(); + r.cache.resize(n_groups); // make sure scratch space is allocated + + const std::size_t n_params = r.n_params(); + const auto tol = + std::sqrt(n_params) * r.n_samples() * radmeth_optimize_params::tolerance; + // clang-format off + auto loglik_bundle = gsl_multimin_function_fdf{ + &neg_loglik, // objective function + &neg_gradient, // gradient + &neg_loglik_and_grad, // combined obj and grad + n_params, // number of model params + static_cast(&r) // parameters for objective and gradient functions + }; + // clang-format on + + // set the parameters: zero for "p" parameters and the final one for + // dispersion using the constant + auto params = gsl_vector_alloc(n_params); + gsl_vector_set_all(params, 0.0); + gsl_vector_set(params, n_params - 1, init_dispersion_param); + + // Alternatives: + // - gsl_multimin_fdfminimizer_conjugate_pr + // - gsl_multimin_fdfminimizer_conjugate_fr + // - gsl_multimin_fdfminimizer_vector_bfgs2 + // - gsl_multimin_fdfminimizer_steepest_descent + const auto minimizer = gsl_multimin_fdfminimizer_vector_bfgs2; + auto s = gsl_multimin_fdfminimizer_alloc(minimizer, n_params); + + gsl_multimin_fdfminimizer_set(s, &loglik_bundle, params, stepsize, tol); + + int status = 0; + std::size_t iter = 0; + + do { + status = gsl_multimin_fdfminimizer_iterate(s); // one iter and get status + if (status) + break; + // check status from gradient + status = gsl_multimin_test_gradient(s->gradient, tol); + } while (status == GSL_CONTINUE && ++iter < max_iter); + // ADS: can't use (status != GSL_SUCCESS) + + const auto param_estimates = gsl_multimin_fdfminimizer_x(s); + const auto &groups = r.design.groups; + p_estimates.clear(); + for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) + p_estimates.push_back(get_p(groups[g_idx], param_estimates)); + const auto disp_param = gsl_vector_get(param_estimates, n_params - 1); + dispersion_estimate = 1.0 / std::exp(disp_param); + + r.max_loglik = log_likelihood(param_estimates, r); + + gsl_multimin_fdfminimizer_free(s); + gsl_vector_free(params); +} + +void +fit_regression_model_gamma(Regression &r, + std::vector &p_estimates, + double &dispersion_estimate) { + fit_regression_model(r, p_estimates, dispersion_estimate); +} + +void +fit_regression_model_gamma(Regression &r, + std::vector &p_estimates, + double &dispersion_estimate) { + fit_regression_model(r, p_estimates, dispersion_estimate); +} + +void +fit_regression_model_gamma_fdf(Regression &r, + std::vector &p_estimates, + double &dispersion_estimate) { + fit_regression_model_fdf(r, p_estimates, dispersion_estimate); +} + +void +fit_regression_model_gamma_fdf(Regression &r, + std::vector &p_estimates, + double &dispersion_estimate) { + fit_regression_model_fdf(r, p_estimates, dispersion_estimate); +} diff --git a/src/radmeth/radmeth_optimize_gamma.hpp b/src/radmeth/radmeth_optimize_gamma.hpp new file mode 100644 index 00000000..1949b6fb --- /dev/null +++ b/src/radmeth/radmeth_optimize_gamma.hpp @@ -0,0 +1,44 @@ +/* Copyright (C) 2025 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 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. + */ + +#ifndef RADMETH_OPTIMIZE_GAMMA_HPP +#define RADMETH_OPTIMIZE_GAMMA_HPP + +#include +#include + +template struct Regression; + +void +fit_regression_model_gamma(Regression &r, + std::vector &p_estimates, + double &dispersion_estimate); + +void +fit_regression_model_gamma(Regression &r, + std::vector &p_estimates, + double &dispersion_estimate); + +void +fit_regression_model_gamma_fdf(Regression &r, + std::vector &p_estimates, + double &dispersion_estimate); + +void +fit_regression_model_gamma_fdf(Regression &r, + std::vector &p_estimates, + double &dispersion_estimate); + +#endif // RADMETH_OPTIMIZE_GAMMA_HPP diff --git a/src/radmeth/radmeth_optimize_params.hpp b/src/radmeth/radmeth_optimize_params.hpp new file mode 100644 index 00000000..a9cbd64a --- /dev/null +++ b/src/radmeth/radmeth_optimize_params.hpp @@ -0,0 +1,27 @@ +/* Copyright (C) 2025 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 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. + */ + +#ifndef RADMETH_OPTIMIZE_PARAMS_HPP +#define RADMETH_OPTIMIZE_PARAMS_HPP + +#include + +namespace radmeth_optimize_params { +inline double tolerance = 1e-4; +inline double stepsize = 0.01; +inline std::uint32_t max_iter = 250; +}; // namespace radmeth_optimize_params + +#endif // RADMETH_OPTIMIZE_PARAMS_HPP diff --git a/src/radmeth/radmeth_optimize_series.cpp b/src/radmeth/radmeth_optimize_series.cpp new file mode 100644 index 00000000..bd54f087 --- /dev/null +++ b/src/radmeth/radmeth_optimize_series.cpp @@ -0,0 +1,308 @@ +/* Copyright (C) 2025 Andrew D. Smith + * + * Author: Andrew D. Smith + * + * This program is free software: you can redistribute it and/or modify it + * under the terms 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. + */ + +#include "radmeth_optimize_series.hpp" + +#include "radmeth_model.hpp" +#include "radmeth_optimize_params.hpp" + +#include +#include + +#include +#include +#include +#include +#include + +[[nodiscard]] static inline auto +logistic(const double x) -> double { + return 1.0 / (1.0 / std::exp(x) + 1.0); +} + +template +[[nodiscard]] static inline auto +get_p(const std::vector &v, const gsl_vector *params) -> double { + const auto a = v.data(); + return logistic(std::inner_product(a, a + std::size(v), params->data, 0.0)); +} + +static inline auto +set_max_r_count(Regression ®) { + const auto &cumul = reg.cumul; + const auto max_itr = std::max_element( + std::cbegin(cumul), std::cend(cumul), [](const auto &a, const auto &b) { + return std::size(a.r_counts) < std::size(b.r_counts); + }); + reg.max_r_count = std::size(max_itr->r_counts); + // ADS: avoid the realloc that can happen even for resize(smaller_size) + if (reg.max_r_count > std::size(reg.cache_log1p_factors)) { + reg.cache_log1p_factors.resize(reg.max_r_count); + reg.cache_dispersion_effect.resize(reg.max_r_count); + } +} + +static inline auto +get_cached_log1p_factors(Regression ®, const double phi) { + const std::size_t max_k = reg.max_r_count; + auto &cache_log1p_factors = reg.cache_log1p_factors; + for (std::size_t k = 0; k < max_k; ++k) + cache_log1p_factors[k] = std::log1p(phi * (k - 1.0)); +} + +static inline auto +get_cached_dispersion_effect(Regression ®, const double phi) { + const std::size_t max_k = reg.max_r_count; + auto &cache_dispersion_effect = reg.cache_dispersion_effect; + double j = -1.0; + const auto lim = std::cbegin(cache_dispersion_effect) + max_k; + for (auto it = std::begin(cache_dispersion_effect); it != lim; ++it, ++j) + *it = j / (1.0 + phi * j); +} + +[[nodiscard]] static double +log_likelihood(const gsl_vector *params, Regression ®) { + const auto phi = logistic(gsl_vector_get(params, reg.design.n_factors())); + const auto one_minus_phi = 1.0 - phi; + + const auto n_groups = reg.n_groups(); + const auto &groups = reg.design.groups; + const auto &cumul = reg.cumul; + + // ADS: precompute log1p(phi * (k - 1.0)), reused for each group + get_cached_log1p_factors(reg, phi); + const auto &log1p_fact_v = reg.cache_log1p_factors; + + double log_lik = 0.0; + for (std::size_t g_idx = 0; g_idx < n_groups; ++g_idx) { + const auto p = get_p(groups[g_idx], params); + const auto one_minus_p = 1.0 - p; + + const auto term1 = one_minus_phi * p; + const auto &cumul_y = cumul[g_idx].m_counts; + std::size_t sz = std::size(cumul_y); + for (std::size_t k = 0; k < sz; ++k) + log_lik += cumul_y[k] * std::log(term1 + phi * k); + + const auto term2 = one_minus_phi * one_minus_p; + const auto &cumul_d = cumul[g_idx].d_counts; + sz = std::size(cumul_d); + for (std::size_t k = 0; k < sz; ++k) + log_lik += cumul_d[k] * std::log(term2 + phi * k); + + const auto &cumul_n = cumul[g_idx].r_counts; + sz = std::size(cumul_n); + for (std::size_t k = 0; k < sz; ++k) + log_lik -= cumul_n[k] * log1p_fact_v[k]; + } + return log_lik; +} + +static void +gradient(const gsl_vector *params, Regression ®, + gsl_vector *output) { + const auto n_factors = reg.design.n_factors(); + const auto phi = logistic(gsl_vector_get(params, n_factors)); + const auto one_minus_phi = 1.0 - phi; + + const auto n_groups = reg.n_groups(); + const auto &groups = reg.design.groups; + const auto &cumul = reg.cumul; + + auto &p_v = reg.p_v; + for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) + p_v[g_idx] = get_p(groups[g_idx], params); + + get_cached_dispersion_effect(reg, phi); // (k-1)/(1 + phi(k-1)) + const auto &dispersion_effect = reg.cache_dispersion_effect; + + // init output to zero for all factors + gsl_vector_set_all(output, 0.0); + const auto &data = output->data; + + double disp_deriv = 0.0; + for (std::size_t g_idx = 0; g_idx < n_groups; ++g_idx) { + const auto p = p_v[g_idx]; + const auto one_minus_p = 1.0 - p; + + double deriv = 0.0; + + const auto denom_term1 = one_minus_phi * p; + const auto &cumul_y = cumul[g_idx].m_counts; + const auto y_lim = std::size(cumul_y); + for (auto k = 0u; k < y_lim; ++k) { + const auto common_factor = cumul_y[k] / (denom_term1 + phi * k); + deriv += common_factor; + disp_deriv += (k - p) * common_factor; + } + + const auto denom_term2 = one_minus_phi * one_minus_p; + const auto &cumul_d = cumul[g_idx].d_counts; + const auto d_lim = std::size(cumul_d); + for (auto k = 0u; k < d_lim; ++k) { + const auto common_factor = cumul_d[k] / (denom_term2 + phi * k); + deriv -= common_factor; + disp_deriv += (k - one_minus_p) * common_factor; + } + + const auto &cumul_n = cumul[g_idx].r_counts; + const auto n_lim = std::size(cumul_n); + for (auto k = 0u; k < n_lim; ++k) + disp_deriv -= cumul_n[k] * dispersion_effect[k]; + + const auto &g = groups[g_idx]; + const auto denom_term1_one_minus_p = denom_term1 * one_minus_p; + for (auto fact_idx = 0u; fact_idx < n_factors; ++fact_idx) { + const auto level = g[fact_idx]; + if (level == 0) + continue; + data[fact_idx] += deriv * (denom_term1_one_minus_p * level); + } + } + gsl_vector_set(output, n_factors, disp_deriv * (phi * one_minus_phi)); +} + +[[nodiscard]] static double +neg_loglik(const gsl_vector *params, void *object) { + auto reg = static_cast *>(object); + return -log_likelihood(params, *reg); +} + +static void +neg_gradient(const gsl_vector *params, void *object, gsl_vector *output) { + auto reg = static_cast *>(object); + gradient(params, *reg, output); + gsl_vector_scale(output, -1.0); +} + +static void +neg_loglik_and_grad(const gsl_vector *params, void *object, double *loglik_val, + gsl_vector *d_loglik_val) { + *loglik_val = neg_loglik(params, object); + neg_gradient(params, object, d_loglik_val); +} + +static void +get_cumulative(const std::vector &group_id, + const std::uint32_t n_groups, + const std::vector> &mc, + std::vector &cumul) { + const auto n_cols = std::size(mc); + cumul.clear(); + cumul.resize(n_groups); + + const auto comp_cumul = [&](auto get_value, auto get_vector) { + // phase 1: determine max value for each group + for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) { + std::uint32_t max_v{}; + for (auto c_idx = 0u; c_idx < n_cols; ++c_idx) { + if (group_id[c_idx] == g_idx) { + const auto val = get_value(mc[c_idx]); + if (val > max_v) + max_v = val; + } + } + get_vector(cumul[g_idx]).resize(max_v, 0); + } + + // phase 2: fill cumulative counts + for (auto c_idx = 0u; c_idx < n_cols; ++c_idx) { + const auto g_idx = group_id[c_idx]; + const auto val = get_value(mc[c_idx]); + auto &vec = get_vector(cumul[g_idx]); + for (auto i = 0u; i < val; ++i) + ++vec[i]; + } + }; + // call 3 times: m_counts, r_counts, d_counts + comp_cumul( + [](const mcounts &m) { return m.n_meth; }, + [](cumul_counts &c) -> std::vector & { return c.m_counts; }); + + comp_cumul( + [](const mcounts &m) { return m.n_reads; }, + [](cumul_counts &c) -> std::vector & { return c.r_counts; }); + + comp_cumul( + [](const mcounts &m) { return m.n_reads - m.n_meth; }, + [](cumul_counts &c) -> std::vector & { return c.d_counts; }); +} + +void +fit_regression_model(Regression &r, + std::vector &p_estimates, + double &dispersion_estimate) { + static constexpr auto init_dispersion_param = -2.5; + const auto stepsize = radmeth_optimize_params::stepsize; + const auto max_iter = radmeth_optimize_params::max_iter; + + const auto n_groups = r.n_groups(); + get_cumulative(r.design.group_id, n_groups, r.mc, r.cumul); + set_max_r_count(r); + + r.p_v.resize(n_groups); + + const std::size_t n_params = r.n_params(); + const auto tol = + std::sqrt(n_params) * r.n_samples() * radmeth_optimize_params::tolerance; + // clang-format off + auto loglik_bundle = gsl_multimin_function_fdf{ + &neg_loglik, // objective function + &neg_gradient, // gradient + &neg_loglik_and_grad, // combined obj and grad + n_params, // number of model params + static_cast(&r) // params for objective and gradients + }; + // clang-format on + + // parameters: 0 for "p" params and final one is a constant for dispersion + auto params = gsl_vector_alloc(n_params); + gsl_vector_set_all(params, 0.0); + gsl_vector_set(params, n_params - 1, init_dispersion_param); + + // Alternatives: + // - gsl_multimin_fdfminimizer_conjugate_pr + // - gsl_multimin_fdfminimizer_conjugate_fr + // - gsl_multimin_fdfminimizer_vector_bfgs2 + // - gsl_multimin_fdfminimizer_steepest_descent + const auto minimizer = gsl_multimin_fdfminimizer_vector_bfgs2; + auto s = gsl_multimin_fdfminimizer_alloc(minimizer, n_params); + gsl_multimin_fdfminimizer_set(s, &loglik_bundle, params, stepsize, tol); + + int status = 0; + std::size_t iter = 0; + + do { + status = gsl_multimin_fdfminimizer_iterate(s); // one iter and get status + if (status) + break; + // check status from gradient + status = gsl_multimin_test_gradient(s->gradient, tol); + } while (status == GSL_CONTINUE && ++iter < max_iter); + // ADS: can't use (status != GSL_SUCCESS) + + const auto param_estimates = gsl_multimin_fdfminimizer_x(s); + const auto &groups = r.design.groups; + p_estimates.clear(); + for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) + p_estimates.push_back(get_p(groups[g_idx], param_estimates)); + const auto disp_param = gsl_vector_get(param_estimates, n_params - 1); + dispersion_estimate = 1.0 / std::exp(disp_param); + + r.max_loglik = log_likelihood(param_estimates, r); + + gsl_multimin_fdfminimizer_free(s); + gsl_vector_free(params); +} diff --git a/src/radmeth/radmeth_optimize_series.hpp b/src/radmeth/radmeth_optimize_series.hpp new file mode 100644 index 00000000..12b4b312 --- /dev/null +++ b/src/radmeth/radmeth_optimize_series.hpp @@ -0,0 +1,29 @@ +/* Copyright (C) 2025 Andrew D. + * + * 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 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. + */ + +#ifndef RADMETH_OPTIMIZE_SERIES_HPP +#define RADMETH_OPTIMIZE_SERIES_HPP + +#include +#include + +template struct Regression; + +void +fit_regression_model(Regression &r, + std::vector &p_estimates, + double &dispersion_estimate); + +#endif // RADMETH_OPTIMIZE_SERIES_HPP diff --git a/src/radmeth/radmeth_utils.cpp b/src/radmeth/radmeth_utils.cpp new file mode 100644 index 00000000..54a39c1c --- /dev/null +++ b/src/radmeth/radmeth_utils.cpp @@ -0,0 +1,144 @@ +/* Copyright (C) 2025 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 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 "radmeth_utils.hpp" + +#include +#include +#include +#include +#include +#include +#include + +[[nodiscard]] std::string +format_duration(const std::chrono::duration elapsed) { + static constexpr auto s_per_h = 3600; + static constexpr auto s_per_m = 60; + const double tot_s = elapsed.count(); + + // break down into hours, minutes, seconds + const std::uint32_t hours = tot_s / 3600; + const std::uint32_t minutes = (static_cast(tot_s) % s_per_h) / s_per_m; + const double seconds = tot_s - (hours * s_per_h) - (minutes * s_per_m); + + std::ostringstream oss; + oss << std::setfill('0') << std::setw(2) << hours << ":" << std::setfill('0') + << std::setw(2) << minutes << ":" << std::fixed << std::setprecision(2) + << std::setw(5) << seconds; + return oss.str(); +} + +file_progress::file_progress(const std::string &filename) : + one_thousand_over_filesize{1000.0 / std::filesystem::file_size(filename)} {} + +void +file_progress::operator()(std::ifstream &in) { + const std::size_t curr_offset = + in.eof() ? 1000 : in.tellg() * one_thousand_over_filesize; + if (curr_offset <= prev_offset) + return; + std::ios old_state(nullptr); + old_state.copyfmt(std::cerr); + std::cerr << "\r[progress: " << std::setw(5) << std::fixed + << std::setprecision(1) << (curr_offset / 10.0) + << (curr_offset == 1000 ? "%]\n" : "%]"); + std::cerr.copyfmt(old_state); + prev_offset = (curr_offset == 1000) ? std::numeric_limits::max() + : curr_offset; +} + +// Series representation for the lower incomplete gamma P(a,x) +[[nodiscard]] static double +gamma_p_series(const double a, const double x) { + constexpr auto eps = std::numeric_limits::epsilon(); + constexpr auto max_iter = 100; + double sum = 1.0 / a; + double term = sum; + for (auto n = 1; n < max_iter; ++n) { + term *= x / (a + n); + sum += term; + if (term < eps * sum) + break; + } + return sum * std::exp(-x + a * std::log(x) - std::lgamma(a)); +} + +[[nodiscard]] static inline double +safe_floor(const double x, const double floor_val) { + return std::abs(x) < floor_val ? floor_val : x; +} + +// Continued fraction representation for the upper incomplete gamma Q(a,x) +[[nodiscard]] static double +gamma_q_contfrac(const double a, const double x) { + constexpr auto epsilon = std::numeric_limits::epsilon(); + constexpr auto fpmin = std::numeric_limits::min() / epsilon; + constexpr auto max_iter = 100; + double b = x + 1 - a; + double c = 1.0 / fpmin; + double d = 1.0 / b; + double h = d; + for (auto i = 1; i < max_iter; ++i) { + const double an = -i * (i - a); + b += 2; + d = safe_floor(an * d + b, fpmin); + c = safe_floor(b + an / c, fpmin); + d = 1.0 / d; + const double delta = d * c; + h *= delta; + if (std::abs(delta - 1.0) < epsilon) + break; + } + return std::exp(-x + a * std::log(x) - std::lgamma(a)) * h; +} + +// Regularized lower incomplete gamma P(a,x) +[[nodiscard]] static double +gamma_p(const double a, const double x) { + if (x < 0 || a <= 0) + return 0.0; + if (x == 0) + return 0.0; + if (x < a + 1.0) + return gamma_p_series(a, x); + return 1.0 - gamma_q_contfrac(a, x); +} + +// chi-square CDF: P(k/2, x/2) +[[nodiscard]] static double +chi_square_cdf(const double x, const double k) { + return gamma_p(k * 0.5, x * 0.5); +} + +// Given the maximum likelihood estimates of the full and reduced models, the +// function outputs the p-value of the log-likelihood ratio. *Note* that it is +// assumed that the reduced model has one fewer factor than the reduced model. +[[nodiscard]] double +llr_test(const double null_loglik, const double full_loglik) { + // The log-likelihood ratio statistic. + const double log_lik_stat = -2 * (null_loglik - full_loglik); + + // It is assumed that null model has one fewer factor than the full + // model. Hence the number of degrees of freedom is 1. + const std::size_t degrees_of_freedom = 1; + + // Log-likelihood ratio statistic has a chi-sqare distribution. + const double chisq_p = chi_square_cdf(log_lik_stat, degrees_of_freedom); + + const double p_value = 1.0 - chisq_p; + + return p_value; +} diff --git a/src/radmeth/radmeth_utils.hpp b/src/radmeth/radmeth_utils.hpp new file mode 100644 index 00000000..cae962c1 --- /dev/null +++ b/src/radmeth/radmeth_utils.hpp @@ -0,0 +1,43 @@ +/* Copyright (C) 2025 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 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. + */ + +#ifndef RADMETH_UTILS_HPP +#define RADMETH_UTILS_HPP + +#include +#include +#include +#include + +[[nodiscard]] std::string +format_duration(const std::chrono::duration elapsed); + +struct file_progress { + double one_thousand_over_filesize{}; + std::size_t prev_offset{}; + explicit file_progress(const std::string &filename); + void + operator()(std::ifstream &in); +}; + +[[nodiscard]] double +llr_test(const double null_loglik, const double full_loglik); + +[[nodiscard]] inline double +overdispersion_factor(const std::uint32_t n_samples, const double dispersion) { + return (n_samples - 1) / (dispersion + 1); +} + +#endif // RADMETH_UTILS_HPP