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

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,

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

@@ -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)))
 })

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)
+})