From be71415aa824feff7f6bc6e30648b07b80ee5e34 Mon Sep 17 00:00:00 2001 From: Yining97 Date: Thu, 11 Jun 2026 11:26:40 -0400 Subject: [PATCH 1/3] feat(univariate): make af the single source of truth (maf optional, derived from af) Step 1.2 follow-up to rss-top-loci-af. In univariate_analysis_pipeline(): - maf becomes optional (still 3rd positional arg; existing callers unaffected). - af is the single source of truth: when af is supplied, the maf_cutoff filtering MAF is derived from it (min(af, 1-af)); a supplied maf is a fallback only and, if it disagrees with the af-derived value, af wins (with a warning). - neither af nor maf supplied + maf_cutoff set -> clear error; maf is never exported. - roxygen regenerated; man page updated. Tests (test_univariate_pipeline.R): af-only -> derived; maf-only -> af NA; both-disagree -> warning; neither+maf_cutoff -> error; neither+no-cutoff -> OK; plus a behavioral test that af-derived MAF filters identically to the equivalent maf (sub-cutoff variant removed; af>0.5 ensures the minor-allele derivation). Co-Authored-By: Claude Opus 4.8 --- R/univariate_pipeline.R | 44 +++++++++--- man/univariate_analysis_pipeline.Rd | 11 ++- tests/testthat/test_univariate_pipeline.R | 88 +++++++++++++++++++++++ 3 files changed, 130 insertions(+), 13 deletions(-) diff --git a/R/univariate_pipeline.R b/R/univariate_pipeline.R index 9ba9c95f..c22ca40e 100644 --- a/R/univariate_pipeline.R +++ b/R/univariate_pipeline.R @@ -8,8 +8,13 @@ #' @param Y A vector of phenotype measurements. #' @param X_scalar A scalar or vector to rescale X to its original scale. #' @param Y_scalar A scalar to rescale Y to its original scale. -#' @param maf A vector of minor allele frequencies for each variant in X. Used -#' only for \code{maf_cutoff} filtering; never exported. +#' @param maf Optional vector of minor allele frequencies for each variant in X, +#' used ONLY for \code{maf_cutoff} filtering and never exported. \code{af} is +#' the single source of truth: when \code{af} is supplied the filtering MAF is +#' derived from it (\code{min(af, 1 - af)}) and a supplied \code{maf} is +#' ignored (with a warning if they disagree). Default NULL; if neither +#' \code{maf} nor \code{af} is supplied and \code{maf_cutoff} is set, the call +#' errors. #' @param af Optional vector of directional effect-allele frequencies (frequency #' of \code{a1}) aligned to the columns of X. When supplied it is exported as #' the \code{top_loci$af} column; when NULL, \code{af} is \code{NA_real_}. @@ -64,7 +69,7 @@ univariate_analysis_pipeline <- function( # input data X, Y, - maf, + maf = NULL, af = NULL, X_scalar = 1, Y_scalar = 1, @@ -99,15 +104,34 @@ univariate_analysis_pipeline <- function( if (!is.matrix(X) || !is.numeric(X)) stop("X must be a numeric matrix") if (!is.vector(Y) && !(is.matrix(Y) && ncol(Y) == 1) || !is.numeric(Y)) stop("Y must be a numeric vector or a single column matrix") if (nrow(X) != length(Y)) stop("X and Y must have the same number of rows/length") - if (!is.numeric(maf) || length(maf) != ncol(X)) stop("maf must be a numeric vector with length equal to the number of columns in X") - if (any(maf < 0 | maf > 1)) stop("maf values must be between 0 and 1") - # af (directional effect-allele frequency) is optional. When supplied it must - # align with X columns; it is exported as top_loci$af. maf stays directionless - # and is used only for maf_cutoff filtering. + # maf is optional (directionless, used ONLY for maf_cutoff filtering, never + # exported). af (directional effect-allele frequency) is optional and is the + # single source of truth: when supplied it is exported as top_loci$af and the + # filtering MAF is derived from it. + if (!is.null(maf)) { + if (!is.numeric(maf) || length(maf) != ncol(X)) stop("maf must be NULL or a numeric vector with length equal to the number of columns in X") + if (any(maf < 0 | maf > 1, na.rm = TRUE)) stop("maf values must be between 0 and 1") + } if (!is.null(af)) { if (!is.numeric(af) || length(af) != ncol(X)) stop("af must be NULL or a numeric vector with length equal to the number of columns in X") if (any(af < 0 | af > 1, na.rm = TRUE)) stop("af values must be between 0 and 1") } + # Single source of truth = af. When af is available, derive the filtering MAF + # from it (min(af, 1 - af)); a supplied directionless maf is only a fallback + # and, if it disagrees with the af-derived value, af wins (with a warning). + if (!is.null(af)) { + af_derived_maf <- pmin(af, 1 - af) + if (!is.null(maf) && any(abs(maf - af_derived_maf) > 1e-6, na.rm = TRUE)) { + warning("Both 'maf' and 'af' were supplied and disagree; using the ", + "af-derived MAF for filtering (af is the single source of truth).") + } + maf <- af_derived_maf + } + # If a MAF cutoff is requested, a frequency must be available to derive it. + if (is.null(maf) && !is.null(maf_cutoff) && is.numeric(maf_cutoff) && maf_cutoff > 0) { + stop("maf_cutoff is set but neither 'af' nor 'maf' was supplied; provide ", + "one so MAF can be derived for filtering.") + } if (!is.numeric(X_scalar) || (length(X_scalar) != 1 && length(X_scalar) != ncol(X))) stop("X_scalar must be a numeric scalar or vector with length equal to the number of columns in X") if (!is.numeric(Y_scalar) || length(Y_scalar) != 1) stop("Y_scalar must be a numeric scalar") if (!is.numeric(L) || L <= 0) stop("L must be a positive integer") @@ -165,7 +189,7 @@ univariate_analysis_pipeline <- function( if (!is.null(ld_reference_meta_file)) { variants_kept <- filter_variants_by_ld_reference(colnames(X), ld_reference_meta_file) X <- X[, variants_kept$data, drop = FALSE] - maf <- maf[variants_kept$idx] + if (!is.null(maf)) maf <- maf[variants_kept$idx] if (!is.null(af)) af <- af[variants_kept$idx] if (length(X_scalar) > 1) X_scalar <- X_scalar[variants_kept$idx] } @@ -174,7 +198,7 @@ univariate_analysis_pipeline <- function( if (!is.null(imiss_cutoff) || !is.null(maf_cutoff)) { X_filtered <- filter_X(X, imiss_cutoff, maf_cutoff, var_thresh = xvar_cutoff, maf = maf, X_variance = X_variance) kept_indices <- match(colnames(X_filtered), colnames(X)) - maf <- maf[kept_indices] + if (!is.null(maf)) maf <- maf[kept_indices] if (!is.null(af)) af <- af[kept_indices] if (length(X_scalar) > 1) X_scalar <- X_scalar[kept_indices] X <- X_filtered diff --git a/man/univariate_analysis_pipeline.Rd b/man/univariate_analysis_pipeline.Rd index 95b3943e..84a355b0 100644 --- a/man/univariate_analysis_pipeline.Rd +++ b/man/univariate_analysis_pipeline.Rd @@ -7,7 +7,7 @@ univariate_analysis_pipeline( X, Y, - maf, + maf = NULL, af = NULL, X_scalar = 1, Y_scalar = 1, @@ -41,8 +41,13 @@ univariate_analysis_pipeline( \item{Y}{A vector of phenotype measurements.} -\item{maf}{A vector of minor allele frequencies for each variant in X. Used -only for \code{maf_cutoff} filtering; never exported.} +\item{maf}{Optional vector of minor allele frequencies for each variant in X, +used ONLY for \code{maf_cutoff} filtering and never exported. \code{af} is +the single source of truth: when \code{af} is supplied the filtering MAF is +derived from it (\code{min(af, 1 - af)}) and a supplied \code{maf} is +ignored (with a warning if they disagree). Default NULL; if neither +\code{maf} nor \code{af} is supplied and \code{maf_cutoff} is set, the call +errors.} \item{af}{Optional vector of directional effect-allele frequencies (frequency of \code{a1}) aligned to the columns of X. When supplied it is exported as diff --git a/tests/testthat/test_univariate_pipeline.R b/tests/testthat/test_univariate_pipeline.R index 219d1a46..99fc755d 100644 --- a/tests/testthat/test_univariate_pipeline.R +++ b/tests/testthat/test_univariate_pipeline.R @@ -2414,3 +2414,91 @@ test_that("rss: mixture LD_data (list of X panels) preserves list shape into sus expect_true(all(sapply(susie_X_arg, is.matrix))) expect_true(any(grepl("susie_rss", names(result)))) }) + +# =========================================================================== +# rss-af-single-source (Step 1.2): af is the single source of truth +# =========================================================================== +.uap_af_xy <- function(n = 20, p = 5, seed = 7) { + set.seed(seed) + X <- matrix(rnorm(n * p), n, p) + colnames(X) <- paste0("chr1:", seq_len(p), ":A:G") + rownames(X) <- paste0("s", seq_len(n)) + Y <- as.numeric(X[, 1] * 2 + rnorm(n)) + names(Y) <- rownames(X) + list(X = X, Y = Y, p = p) +} + +test_that("af supplied (no maf): MAF derived from af, af exported, no maf column", { + skip_if_not_installed("susieR") + d <- .uap_af_xy() + af <- runif(d$p, 0.1, 0.5) + res <- univariate_analysis_pipeline(X = d$X, Y = d$Y, af = af, + maf_cutoff = 0.01, twas_weights = FALSE, L = 5, L_greedy = 5) + expect_true("top_loci" %in% names(res)) + expect_true("af" %in% colnames(res$top_loci)) + expect_false("maf" %in% colnames(res$top_loci)) + if (nrow(res$top_loci) > 0) expect_true(all(!is.na(res$top_loci$af))) +}) + +test_that("maf supplied (no af): maf used for filtering, top_loci$af is NA", { + skip_if_not_installed("susieR") + d <- .uap_af_xy() + maf <- runif(d$p, 0.1, 0.5) + res <- univariate_analysis_pipeline(X = d$X, Y = d$Y, maf = maf, + maf_cutoff = 0.01, twas_weights = FALSE, L = 5, L_greedy = 5) + expect_true("af" %in% colnames(res$top_loci)) + if (nrow(res$top_loci) > 0) expect_true(all(is.na(res$top_loci$af))) +}) + +test_that("both af and maf supplied and disagreeing emits a warning (af wins)", { + skip_if_not_installed("susieR") + d <- .uap_af_xy() + af <- rep(0.30, d$p) # min(af, 1-af) = 0.30 + maf <- rep(0.40, d$p) # disagrees with the af-derived 0.30 + expect_warning( + univariate_analysis_pipeline(X = d$X, Y = d$Y, af = af, maf = maf, + maf_cutoff = 0.01, twas_weights = FALSE, L = 5, L_greedy = 5), + "disagree" + ) +}) + +test_that("neither af nor maf with maf_cutoff set errors clearly", { + skip_if_not_installed("susieR") + d <- .uap_af_xy() + expect_error( + univariate_analysis_pipeline(X = d$X, Y = d$Y, + maf_cutoff = 0.01, twas_weights = FALSE, L = 5, L_greedy = 5), + "maf_cutoff is set but neither" + ) +}) + +test_that("neither af nor maf without maf_cutoff runs; top_loci$af is NA", { + skip_if_not_installed("susieR") + d <- .uap_af_xy() + res <- univariate_analysis_pipeline(X = d$X, Y = d$Y, + maf_cutoff = NULL, twas_weights = FALSE, L = 5, L_greedy = 5) + expect_true("top_loci" %in% names(res)) + expect_true("af" %in% colnames(res$top_loci)) + if (nrow(res$top_loci) > 0) expect_true(all(is.na(res$top_loci$af))) +}) + +test_that("af-derived MAF drives maf_cutoff filtering identically to the equivalent maf", { + # Behavioral test of the edit: supplying `af` must filter via min(af, 1-af), + # giving the SAME fit as supplying that directionless maf explicitly. Includes + # an af > 0.5 (variant 2) so a wrong derivation (using af directly, not the + # minor allele) would diverge; and a sub-cutoff variant (variant 1) so the + # cutoff actually removes something. + skip_if_not_installed("susieR") + d <- .uap_af_xy(p = 6) + af <- c(0.005, 0.85, 0.20, 0.45, 0.15, 0.40) # variant 1 below cutoff; variant 2 af>0.5 + res_af <- univariate_analysis_pipeline(X = d$X, Y = d$Y, af = af, + maf_cutoff = 0.01, twas_weights = FALSE, L = 5, L_greedy = 5) + res_maf <- univariate_analysis_pipeline(X = d$X, Y = d$Y, maf = pmin(af, 1 - af), + maf_cutoff = 0.01, twas_weights = FALSE, L = 5, L_greedy = 5) + # Identical on every column except af (res_maf's af is NA) => the af-derived + # MAF was used for filtering exactly like an explicit maf. + cols <- setdiff(names(res_af$top_loci), "af") + expect_equal(res_af$top_loci[cols], res_maf$top_loci[cols]) + # And the sub-cutoff variant (chr1:1:A:G, MAF 0.005 < 0.01) was filtered out: + expect_false("chr1:1:A:G" %in% res_af$top_loci$variant) +}) From d1237ebddec2f233a1dd9f8450fb98c5b58412eb Mon Sep 17 00:00:00 2001 From: Yining97 Date: Thu, 11 Jun 2026 12:56:06 -0400 Subject: [PATCH 2/3] fix: multivariate caller still passed removed maf arg to postprocess_finemapping_fits The maf->af rename in #513 missed multivariate_analysis_pipeline, which broke test_multivariate_pipeline.R. The multivariate maf is directionless and only used for filtering, so leave af unset (top_loci$af = NA) instead of relabeling it. Added a test. Pending follow-up: give multivariate_analysis_pipeline its own optional directional af input so it can export a real af, like univariate does. --- R/multivariate_pipeline.R | 5 ++++- tests/testthat/test_multivariate_pipeline.R | 6 ++++++ 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/R/multivariate_pipeline.R b/R/multivariate_pipeline.R index c75a9a18..0e313654 100644 --- a/R/multivariate_pipeline.R +++ b/R/multivariate_pipeline.R @@ -432,7 +432,10 @@ multivariate_analysis_pipeline <- function( data_y = NULL, X_scalar = 1, y_scalar = 1, - maf = maf, + # `maf` here is directionless (used only for filtering above); do NOT export + # it as `af` (effect-allele frequency). Leave `af` NULL -> top_loci$af = NA + # until the multivariate path carries a proven directional AF (follow-up, + # mirroring the univariate `af-single-source` work). coverage = coverage[1], secondary_coverage = sec_coverage, signal_cutoff = signal_cutoff, diff --git a/tests/testthat/test_multivariate_pipeline.R b/tests/testthat/test_multivariate_pipeline.R index 9ef1e33e..4a37e671 100644 --- a/tests/testthat/test_multivariate_pipeline.R +++ b/tests/testthat/test_multivariate_pipeline.R @@ -619,4 +619,10 @@ test_that("pipeline propagates outcome_names from mvsusie through post-processin trimmed <- getTrimmedFit(result$finemapping_result) expect_equal(trimmed$coef, fake_coef[-1, , drop = FALSE]) expect_equal(dim(trimmed$clfsr), c(L, p, r)) + # rss-top-loci-af regression: top_loci exports `af`, not the directionless + # `maf`. The multivariate path has no proven directional AF, so it must NOT + # relabel its filtering `maf` as `af` — `af` is present but NA, no `maf` column. + expect_true("af" %in% colnames(result$top_loci)) + expect_false("maf" %in% colnames(result$top_loci)) + expect_true(all(is.na(result$top_loci$af))) }) From b5026875be5063eb42aeee231939be3d0e52bf5f Mon Sep 17 00:00:00 2001 From: Yining97 Date: Thu, 11 Jun 2026 15:29:50 -0400 Subject: [PATCH 3/3] multivariate: make maf optional and add af, matching univariate Same af handling as univariate_analysis_pipeline: maf optional, optional directional af exported as top_loci$af, filtering MAF derived from af when supplied (maf is a fallback, never exported), error when neither is given with maf_cutoff. Updates the old maf-required validation tests and adds an af-supplied test. test_multivariate_pipeline.R: 0 fail / 53 pass. --- R/multivariate_pipeline.R | 50 ++++++++++++++---- man/multivariate_analysis_pipeline.Rd | 15 +++++- tests/testthat/test_multivariate_pipeline.R | 56 +++++++++++++++++++-- 3 files changed, 104 insertions(+), 17 deletions(-) diff --git a/R/multivariate_pipeline.R b/R/multivariate_pipeline.R index 0e313654..936180a3 100644 --- a/R/multivariate_pipeline.R +++ b/R/multivariate_pipeline.R @@ -137,7 +137,16 @@ region_data_to_mvsusie_rss_input <- function(sumstat_data, ld_name = NULL) { #' #' @param X A matrix of genotype data where rows represent samples and columns represent genetic variants. #' @param Y A matrix of phenotype measurements, representing samples and columns represent conditions. -#' @param maf A list of vectors for minor allele frequencies for each variant in X. +#' @param maf Optional vector of minor allele frequencies for each variant in X, +#' used ONLY for \code{maf_cutoff} filtering and never exported. When \code{af} +#' is supplied the filtering MAF is derived from it (\code{min(af, 1 - af)}) and +#' a supplied \code{maf} is ignored (with a warning if they disagree). Default +#' NULL; if neither \code{maf} nor \code{af} is supplied and \code{maf_cutoff} +#' is set, the call errors. +#' @param af Optional vector of directional effect-allele frequencies (frequency +#' of \code{a1}) aligned to the columns of X. When supplied it is exported as +#' the \code{top_loci$af} column; when NULL, \code{af} is \code{NA_real_}. +#' Default NULL. #' @param L Maximum number of components in mvSuSiE. Default is 30. #' @param L_greedy Initial greedy number of components in mvSuSiE. Default is 5. #' @param ld_reference_meta_file An optional path to a file containing linkage disequilibrium reference data. If provided, variants in X are filtered based on this reference. @@ -199,7 +208,8 @@ multivariate_analysis_pipeline <- function( # input data X, Y, - maf, + maf = NULL, + af = NULL, X_variance = NULL, other_quantities = list(), region = NULL, @@ -332,8 +342,29 @@ multivariate_analysis_pipeline <- function( if (!is.matrix(X) || !is.numeric(X)) stop("X must be a numeric matrix") if (!is.matrix(Y) || !is.numeric(Y)) stop("Y must be a numeric matrix") if (nrow(X) != nrow(Y)) stop("X and Y must have the same number of rows") - if (!is.numeric(maf) || length(maf) != ncol(X)) stop("maf must be a numeric vector with length equal to the number of columns in X") - if (any(maf < 0 | maf > 1)) stop("maf values must be between 0 and 1") + if (!is.null(maf)) { + if (!is.numeric(maf) || length(maf) != ncol(X)) stop("maf must be NULL or a numeric vector with length equal to the number of columns in X") + if (any(maf < 0 | maf > 1, na.rm = TRUE)) stop("maf values must be between 0 and 1") + } + if (!is.null(af)) { + if (!is.numeric(af) || length(af) != ncol(X)) stop("af must be NULL or a numeric vector with length equal to the number of columns in X") + if (any(af < 0 | af > 1, na.rm = TRUE)) stop("af values must be between 0 and 1") + } + # Single source of truth = af. When af is available, derive the filtering MAF + # from it (min(af, 1 - af)); a supplied directionless maf is only a fallback + # and, if it disagrees with the af-derived value, af wins (with a warning). + if (!is.null(af)) { + af_derived_maf <- pmin(af, 1 - af) + if (!is.null(maf) && any(abs(maf - af_derived_maf) > 1e-6, na.rm = TRUE)) { + warning("Both 'maf' and 'af' were supplied and disagree; using the ", + "af-derived MAF for filtering (af is the single source of truth).") + } + maf <- af_derived_maf + } + if (is.null(maf) && !is.null(maf_cutoff) && is.numeric(maf_cutoff) && maf_cutoff > 0) { + stop("maf_cutoff is set but neither 'af' nor 'maf' was supplied; provide ", + "one so MAF can be derived for filtering.") + } if (!is.numeric(L) || L <= 0) stop("L must be a positive integer") if (!is.null(L_greedy) && (!is.numeric(L_greedy) || L_greedy <= 0)) stop("L_greedy must be NULL or a positive integer") if (!is.null(L_greedy)) L_greedy <- min(L_greedy, L) @@ -356,13 +387,15 @@ multivariate_analysis_pipeline <- function( if (!is.null(ld_reference_meta_file)) { variants_kept <- filter_variants_by_ld_reference(colnames(X), ld_reference_meta_file) X <- X[, variants_kept$data, drop = FALSE] - maf <- maf[variants_kept$idx] + if (!is.null(maf)) maf <- maf[variants_kept$idx] + if (!is.null(af)) af <- af[variants_kept$idx] } # filter X based on Y subjects if (!is.null(imiss_cutoff) || !is.null(maf_cutoff)) { X <- filter_X_with_Y(X, Y, imiss_cutoff, maf_cutoff, var_thresh = xvar_cutoff, maf = maf, X_variance = X_variance) - maf <- maf[colnames(X)] + if (!is.null(maf)) maf <- maf[colnames(X)] + if (!is.null(af)) af <- af[colnames(X)] } # filter data driven prior matrices @@ -432,10 +465,7 @@ multivariate_analysis_pipeline <- function( data_y = NULL, X_scalar = 1, y_scalar = 1, - # `maf` here is directionless (used only for filtering above); do NOT export - # it as `af` (effect-allele frequency). Leave `af` NULL -> top_loci$af = NA - # until the multivariate path carries a proven directional AF (follow-up, - # mirroring the univariate `af-single-source` work). + af = af, coverage = coverage[1], secondary_coverage = sec_coverage, signal_cutoff = signal_cutoff, diff --git a/man/multivariate_analysis_pipeline.Rd b/man/multivariate_analysis_pipeline.Rd index e21e6280..51902c2a 100644 --- a/man/multivariate_analysis_pipeline.Rd +++ b/man/multivariate_analysis_pipeline.Rd @@ -7,7 +7,8 @@ multivariate_analysis_pipeline( X, Y, - maf, + maf = NULL, + af = NULL, X_variance = NULL, other_quantities = list(), region = NULL, @@ -41,7 +42,17 @@ multivariate_analysis_pipeline( \item{Y}{A matrix of phenotype measurements, representing samples and columns represent conditions.} -\item{maf}{A list of vectors for minor allele frequencies for each variant in X.} +\item{maf}{Optional vector of minor allele frequencies for each variant in X, +used ONLY for \code{maf_cutoff} filtering and never exported. When \code{af} +is supplied the filtering MAF is derived from it (\code{min(af, 1 - af)}) and +a supplied \code{maf} is ignored (with a warning if they disagree). Default +NULL; if neither \code{maf} nor \code{af} is supplied and \code{maf_cutoff} +is set, the call errors.} + +\item{af}{Optional vector of directional effect-allele frequencies (frequency +of \code{a1}) aligned to the columns of X. When supplied it is exported as +the \code{top_loci$af} column; when NULL, \code{af} is \code{NA_real_}. +Default NULL.} \item{ld_reference_meta_file}{An optional path to a file containing linkage disequilibrium reference data. If provided, variants in X are filtered based on this reference.} diff --git a/tests/testthat/test_multivariate_pipeline.R b/tests/testthat/test_multivariate_pipeline.R index 4a37e671..a8b35b8d 100644 --- a/tests/testthat/test_multivariate_pipeline.R +++ b/tests/testthat/test_multivariate_pipeline.R @@ -150,7 +150,7 @@ test_that("multivariate_analysis_pipeline errors when maf has wrong length", { multivariate_analysis_pipeline( X = d$X, Y = d$Y, maf = c(0.1, 0.2), pip_cutoff_to_skip = 0 ), - "maf must be a numeric vector with length equal to the number of columns in X" + "maf must be NULL or a numeric vector with length equal to the number of columns in X" ) }) @@ -162,7 +162,7 @@ test_that("multivariate_analysis_pipeline errors when maf is not numeric", { multivariate_analysis_pipeline( X = d$X, Y = d$Y, maf = rep("0.1", ncol(d$X)), pip_cutoff_to_skip = 0 ), - "maf must be a numeric vector" + "maf must be NULL or a numeric vector" ) }) @@ -372,15 +372,17 @@ test_that("pipeline warns and returns empty list when Y has single column", { expect_equal(length(result), 0) }) -test_that("multivariate_analysis_pipeline errors when maf is NULL", { +test_that("multivariate_analysis_pipeline errors when maf is NULL and maf_cutoff is set without af", { skip_if(!requireNamespace("mvsusieR", quietly = TRUE), "mvsusieR not installed, input validation unreachable") d <- make_mv_data() + # maf is now optional; it only errors if maf_cutoff needs a frequency and + # neither maf nor af is supplied. expect_error( multivariate_analysis_pipeline( - X = d$X, Y = d$Y, maf = NULL, pip_cutoff_to_skip = 0 + X = d$X, Y = d$Y, maf_cutoff = 0.01, pip_cutoff_to_skip = 0 ), - "maf must be a numeric vector" + "maf_cutoff is set but neither" ) }) @@ -626,3 +628,47 @@ test_that("pipeline propagates outcome_names from mvsusie through post-processin expect_false("maf" %in% colnames(result$top_loci)) expect_true(all(is.na(result$top_loci$af))) }) + +test_that("multivariate exports a supplied af and derives the filtering MAF from it", { + skip_if(!requireNamespace("mvsusieR", quietly = TRUE), "mvsusieR not installed") + skip_if(!requireNamespace("susieR", quietly = TRUE), "susieR not installed") + d <- make_mv_data() + r <- ncol(d$Y); p <- ncol(d$X); L <- 5 + vnames <- colnames(d$X); cnames <- colnames(d$Y) + fake_coef <- matrix(rnorm((p + 1) * r), nrow = p + 1, ncol = r) + prior_U <- list(udd_1 = matrix(0.5, r, r) + 0.5 * diag(r), udd_2 = diag(r)) + for (k in seq_along(prior_U)) rownames(prior_U[[k]]) <- colnames(prior_U[[k]]) <- cnames + prior_mats <- list(U = prior_U, w = c(udd_1 = 0.5, udd_2 = 0.5)) + fake_w0 <- c("udd_1_a" = 0.4, "udd_2_a" = 0.4, "null" = 0.2) + local_mocked_bindings( + mrmash_wrapper = function(X, Y, ...) list(V = diag(ncol(Y)), w0 = fake_w0, + w1 = matrix(0.1, nrow = ncol(X), ncol = 1)), + ) + local_mocked_bindings( + mvsusie = function(...) list( + pip = setNames(rep(0.5, p), vnames), # high PIP -> retained as top_loci rows + alpha = matrix(1 / p, nrow = L, ncol = p), + lbf_variable = matrix(0, nrow = L, ncol = p), + V = rep(1, L), sets = list(cs = NULL, requested_coverage = 0.95), + niter = 10, outcome_names = cnames, + conditional_lfsr = array(0.5, dim = c(L, p, r)) + ), + create_mixture_prior = function(...) list(matrices = prior_U, weights = c(0.5, 0.5)), + coef.mvsusie = function(...) fake_coef, + .package = "mvsusieR" + ) + af <- setNames(seq(0.1, 0.5, length.out = p), vnames) # directional, aligned to X cols + run <- function(...) multivariate_analysis_pipeline( + X = d$X, Y = d$Y, ..., pip_cutoff_to_skip = 0, signal_cutoff = 0.025, + data_driven_prior_matrices = prior_mats, twas_weights = FALSE) + res_af <- run(af = af) + res_maf <- run(maf = pmin(af, 1 - af)) + # af supplied -> exported with real values; no maf column + expect_true(nrow(res_af$top_loci) > 0) + expect_false("maf" %in% colnames(res_af$top_loci)) + expect_true(all(!is.na(res_af$top_loci$af))) + expect_true(all(res_af$top_loci$af %in% af)) + # af-derived MAF filters identically to supplying that maf (same fit, only af differs) + cols <- setdiff(names(res_af$top_loci), "af") + expect_equal(res_af$top_loci[cols], res_maf$top_loci[cols]) +})