Step 1.2: make maf optional (derived from af) in univariate and multivariate pipelines

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,7 +465,7 @@ multivariate_analysis_pipeline <- function(
     data_y = NULL,
     X_scalar = 1,
     y_scalar = 1,
-    maf = maf,
+    af = af,
     coverage = coverage[1],
     secondary_coverage = sec_coverage,
     signal_cutoff = signal_cutoff,

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

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"
   )
 })
 
@@ -619,4 +621,54 @@ 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)))
+})
+
+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])
 })