From 2aa40e0860589f1b00b28cac1d68f4bbcb2cd591 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Fri, 12 Sep 2025 16:22:54 -0700 Subject: [PATCH 1/6] src/radmeth/radmeth_model.hpp: added a function to the Regression and Design classes to remove from the Design data structures any samples that do not appear in the big data matrix --- src/radmeth/radmeth_model.hpp | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/src/radmeth/radmeth_model.hpp b/src/radmeth/radmeth_model.hpp index 1c42cb17..fb4b5ac9 100644 --- a/src/radmeth/radmeth_model.hpp +++ b/src/radmeth/radmeth_model.hpp @@ -34,10 +34,10 @@ struct cumul_counts { struct Design { std::vector factor_names; std::vector sample_names; - std::vector> matrix; - std::vector> tmatrix; - std::vector> groups; - std::vector group_id; + 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 { @@ -56,6 +56,9 @@ struct Design { [[nodiscard]] Design drop_factor(const std::size_t factor_idx); + + void + order_samples(const std::vector &ordered_names); }; std::istream & @@ -121,6 +124,11 @@ struct Regression { n_samples() const { return design.n_samples(); } + + void + order_samples(const std::vector &ordered_names) { + design.order_samples(ordered_names); + } }; #endif From 3896cc949e6fd2e744becea051b618bc6b46531c Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Fri, 12 Sep 2025 16:23:50 -0700 Subject: [PATCH 2/6] src/radmeth/radmeth_model.cpp: added the implementation of Design::order_samples which ensures consistent order between the design and the big data matrix, while also removing any samples from the design that do not appear in the data matrix. Untested --- src/radmeth/radmeth_model.cpp | 39 +++++++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/src/radmeth/radmeth_model.cpp b/src/radmeth/radmeth_model.cpp index 5ca694a6..8c77b5f8 100644 --- a/src/radmeth/radmeth_model.cpp +++ b/src/radmeth/radmeth_model.cpp @@ -22,6 +22,7 @@ #include #include #include +#include #include double Regression::tolerance = 1e-4; @@ -167,6 +168,44 @@ Design::drop_factor(const std::size_t factor_idx) { return design; } +void +Design::order_samples(const std::vector &ordered_names) { + const auto lookup = [&] { + std::unordered_map tmp; + std::uint32_t idx{}; + for (const auto &name : ordered_names) + tmp.emplace(name, idx++); + return tmp; + }(); + + auto tmp_sample_names = sample_names; + auto tmp_matrix = matrix; + auto tmp_group_id = group_id; + // factor names should not change + // groups should not change + // tmatrix will be changed after matrix + + const auto n_names = std::size(ordered_names); + for (auto i = 0u; i < std::size(tmp_sample_names); ++i) { + const auto itr = lookup.find(tmp_sample_names[i]); + if (itr == std::cend(lookup)) + throw std::runtime_error("sample not found: " + tmp_sample_names[i]); + const std::uint32_t idx = itr->second; + sample_names[idx] = tmp_sample_names[i]; + matrix[idx] = tmp_matrix[i]; + group_id[idx] = tmp_group_id[i]; + } + + // make sure no extra entries for samples in the design but that don't have + // data in the table + sample_names.resize(n_names); + matrix.resize(n_names); + group_id.resize(n_names); + + // update tmatrix using matrix + transpose(matrix, tmatrix); +} + std::ostream & operator<<(std::ostream &out, const Design &design) { static constexpr std::uint32_t max_samples_to_report = 20; From f9650eb1a392f6ea7655ca051398562484132c22 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Fri, 12 Sep 2025 16:24:43 -0700 Subject: [PATCH 3/6] src/radmeth/radmeth.cpp: using the Regression::order_samples to ensure samples have the same order in the data matrix and the design matrix, rather than just throwing an error when recover should be possibe --- src/radmeth/radmeth.cpp | 24 +++++++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) diff --git a/src/radmeth/radmeth.cpp b/src/radmeth/radmeth.cpp index 9ace1fb4..da28acdb 100644 --- a/src/radmeth/radmeth.cpp +++ b/src/radmeth/radmeth.cpp @@ -78,7 +78,7 @@ struct file_progress { } }; -static bool +[[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)); @@ -90,6 +90,16 @@ consistent_sample_names(const Regression ®, const std::string &header) { 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; +} + // 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. @@ -195,6 +205,7 @@ radmeth(const bool show_progress, const bool more_na_info, 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, sample_names_header)) { // clang-format off const auto message = @@ -435,6 +446,17 @@ main_radmeth(int argc, char *argv[]) { if (verbose) std::cerr << "Null model:\n" << null_model.design << '\n'; + { + 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); + } + // clang-format off if (verbose) std::cerr << "Output columns:\n" From 93880052153066bd222612de0871ac03c7494f68 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Fri, 12 Sep 2025 17:00:11 -0700 Subject: [PATCH 4/6] src/radmeth/radmeth_model.cpp: fixed the subsetting re-ordering of samples which was inverted --- src/radmeth/radmeth_model.cpp | 49 +++++++++++++++++++---------------- 1 file changed, 26 insertions(+), 23 deletions(-) diff --git a/src/radmeth/radmeth_model.cpp b/src/radmeth/radmeth_model.cpp index 8c77b5f8..0df305cc 100644 --- a/src/radmeth/radmeth_model.cpp +++ b/src/radmeth/radmeth_model.cpp @@ -170,37 +170,40 @@ Design::drop_factor(const std::size_t factor_idx) { void Design::order_samples(const std::vector &ordered_names) { - const auto lookup = [&] { - std::unordered_map tmp; - std::uint32_t idx{}; - for (const auto &name : ordered_names) - tmp.emplace(name, idx++); - return tmp; + // Build lookup: sample name -> original index + const auto sample_to_index = [&] { + std::unordered_map sample_to_index; + std::uint32_t index = 0; + for (const auto &sample : sample_names) + sample_to_index.emplace(sample, index++); + return sample_to_index; }(); - auto tmp_sample_names = sample_names; - auto tmp_matrix = matrix; - auto tmp_group_id = group_id; + std::vector ord_sample_names; + std::vector> ord_matrix; + std::vector ord_group_id; // factor names should not change // groups should not change // tmatrix will be changed after matrix const auto n_names = std::size(ordered_names); - for (auto i = 0u; i < std::size(tmp_sample_names); ++i) { - const auto itr = lookup.find(tmp_sample_names[i]); - if (itr == std::cend(lookup)) - throw std::runtime_error("sample not found: " + tmp_sample_names[i]); - const std::uint32_t idx = itr->second; - sample_names[idx] = tmp_sample_names[i]; - matrix[idx] = tmp_matrix[i]; - group_id[idx] = tmp_group_id[i]; + ord_sample_names.reserve(n_names); + ord_matrix.reserve(n_names); + ord_group_id.reserve(n_names); + + for (const auto &name : ordered_names) { + const auto it = sample_to_index.find(name); + if (it == std::cend(sample_to_index)) + throw std::runtime_error("Sample not found: " + name); + + const auto original_index = it->second; + ord_sample_names.push_back(sample_names[original_index]); + ord_matrix.push_back(matrix[original_index]); + ord_group_id.push_back(group_id[original_index]); } - - // make sure no extra entries for samples in the design but that don't have - // data in the table - sample_names.resize(n_names); - matrix.resize(n_names); - group_id.resize(n_names); + sample_names = std::move(ord_sample_names); + matrix = std::move(ord_matrix); + group_id = std::move(ord_group_id); // update tmatrix using matrix transpose(matrix, tmatrix); From 9b22fb9a2bf93eccd6be3eb8f829506f47c3e70a Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sat, 13 Sep 2025 15:03:27 -0700 Subject: [PATCH 5/6] src/radmeth/radmeth.cpp: putting the code to ensure sample order into a function --- src/radmeth/radmeth.cpp | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/src/radmeth/radmeth.cpp b/src/radmeth/radmeth.cpp index da28acdb..511c9c46 100644 --- a/src/radmeth/radmeth.cpp +++ b/src/radmeth/radmeth.cpp @@ -361,6 +361,19 @@ 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 { @@ -446,16 +459,7 @@ main_radmeth(int argc, char *argv[]) { if (verbose) std::cerr << "Null model:\n" << null_model.design << '\n'; - { - 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); - } + ensure_sample_order(table_filename, alt_model, null_model); // clang-format off if (verbose) From cc4d2190ae1417812ce36a941f72672e482fffae Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 14 Sep 2025 18:36:28 -0700 Subject: [PATCH 6/6] src/radmeth/radmeth.cpp: using my own chisquared function instead of GSL; works very similarly --- src/radmeth/radmeth.cpp | 71 ++++++++++++++++++++++++++++++++++++++--- 1 file changed, 67 insertions(+), 4 deletions(-) diff --git a/src/radmeth/radmeth.cpp b/src/radmeth/radmeth.cpp index 511c9c46..789304c5 100644 --- a/src/radmeth/radmeth.cpp +++ b/src/radmeth/radmeth.cpp @@ -23,8 +23,6 @@ #include "smithlab_os.hpp" #include "smithlab_utils.hpp" -#include // GSL header for chisqrd distribution - #include #include #include @@ -32,6 +30,7 @@ #include #include #include +#include #include #include #include @@ -100,10 +99,73 @@ get_sample_names_from_header(const std::string &header) { 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. -static double +[[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); @@ -113,7 +175,8 @@ llr_test(const double null_loglik, const double full_loglik) { const std::size_t degrees_of_freedom = 1; // Log-likelihood ratio statistic has a chi-sqare distribution. - const double chisq_p = gsl_cdf_chisq_P(log_lik_stat, degrees_of_freedom); + const double chisq_p = chi_square_cdf(log_lik_stat, degrees_of_freedom); + const double p_value = 1.0 - chisq_p; return p_value;