Derive Wald z from beta/se when z is absent

R/manifestLoaders.R

@@ -390,8 +390,9 @@ NULL
 # per-variant N_CASE + N_CONTROL are both present (from which summaryStatsQc
 # derives the effective sample size), OR `allowNoN = TRUE` because the caller
 # has a study-level scalar (nCase/nControl or nSample) that summaryStatsQc will
-# fill N from. `mapping` (named std -> source) overrides the default aliases
-# per key.
+# fill N from. Z is required UNLESS BETA + SE are both present, in which case
+# the Wald z (Z = BETA / SE) is derived; a supplied Z always takes precedence.
+# `mapping` (named std -> source) overrides the default aliases per key.
 .resolveSumstatCols <- function(df, columnMapping, label, allowNoN = FALSE) {
   # Strip a leading '#' from the first column name so a '#CHR'-style header
   # (common in GWAS TSVs) resolves the same way on the plain-text and tabix
@@ -415,14 +416,28 @@ NULL
     }
     intersect(.sumstatColumnAliases[[key]], names(df))[1L]
   }
-  required <- c("chrom", "pos", "variant_id", "A1", "A2", "Z")
+  required <- c("chrom", "pos", "variant_id", "A1", "A2")
   resolved <- setNames(lapply(required, resolveKey), required)
   missingKeys <- required[vapply(resolved, function(x) is.na(x) || is.null(x), logical(1))]
   if (length(missingKeys) > 0L) {
     stop(label, ": sumstats file is missing required field(s): ",
          paste(missingKeys, collapse = ", "),
          " (supply a columnMapping if the source names differ).")
   }
+  # z-score: a Z column, OR BETA + SE (from which the Wald z = beta/se is
+  # derived, matching susie_rss()'s z_method="wald"). At least one is required;
+  # when both a Z column and BETA/SE are present, the Z column takes precedence.
+  zSrc      <- resolveKey("Z")
+  betaSrc   <- resolveKey("BETA")
+  seSrc     <- resolveKey("SE")
+  hasZ      <- !is.na(zSrc)    && !is.null(zSrc)
+  hasBetaSe <- !is.na(betaSrc) && !is.null(betaSrc) &&
+               !is.na(seSrc)   && !is.null(seSrc)
+  if (!hasZ && !hasBetaSe) {
+    stop(label, ": sumstats file needs a z field (z/Z), or beta + se to derive ",
+         "the Wald z-score beta/se (supply a columnMapping if the source names ",
+         "differ).")
+  }
   # Sample size: an N column, OR per-variant N_CASE + N_CONTROL counts (from
   # which summaryStatsQc(effectiveN=) derives the effective sample size). At
   # least one is required; both may be present (N_CASE/N_CONTROL take priority
@@ -445,8 +460,10 @@ NULL
     SNP   = as.character(df[[resolved$variant_id]]),
     A1    = as.character(df[[resolved$A1]]),
     A2    = as.character(df[[resolved$A2]]),
-    Z     = as.numeric(df[[resolved$Z]]),
     stringsAsFactors = FALSE)
+  # Z takes precedence when present; otherwise derive the Wald z from beta/se.
+  out$Z <- if (hasZ) as.numeric(df[[zSrc]]) else
+             as.numeric(df[[betaSrc]]) / as.numeric(df[[seSrc]])
   if (hasN) out$N <- as.numeric(df[[nSrc]])
   # Per-variant case/control counts -> N_CASE / N_CONTROL (kept by
   # .dfToEntryGranges; summaryStatsQc(effectiveN=) consumes them).
@@ -827,7 +844,10 @@ loadQtlDatasetFromManifest <- function(manifest, study = NULL,
 
 #' @title Load a GwasSumStats collection from a manifest
 #' @description Build a \code{\link{GwasSumStats}} from a manifest with one row
-#'   per study. No QC is run (the result carries \code{qcInfo = list()}).
+#'   per study. No QC is run (the result carries \code{qcInfo = list()}). Each
+#'   sumstats file needs a \code{z} column, or \code{beta}+\code{se} from which
+#'   the Wald z (\code{z = beta/se}) is derived when \code{z} is absent (a
+#'   supplied \code{z} takes precedence).
 #' @param manifest A data.frame or path. Columns (snake_case aliases accepted):
 #'   \code{study} (required, unique), \code{sumStatsPath} (required),
 #'   \code{columnMapping} (optional), \code{nCase} / \code{nControl}

tests/testthat/test_fineMappingWrappers.R

@@ -1783,6 +1783,10 @@ test_that(".fmFitSusieRss rejects a non-named-list rssControl", {
 })
 
 test_that(".fmFitSusieRss forwards rssControl to susie_rss as control", {
+  # susie_rss_control() only exists in susieR >= 0.16.6; on older susieR the
+  # mock binding can't resolve. rssControl is a 0.16.6-only feature anyway.
+  skip_if_not("susie_rss_control" %in% getNamespaceExports("susieR"),
+              "susieR < 0.16.6 has no susie_rss_control")
   cap <- new.env(parent = emptyenv())
   local_mocked_bindings(
     susie_rss_control = function(...) list(.tag = "ctrl", ...),

tests/testthat/test_manifestLoaders.R

@@ -245,6 +245,47 @@ test_that(".resolveSumstatCols requires an N field or N_CASE + N_CONTROL", {
                "N_CASE + N_CONTROL", fixed = TRUE)
 })
 
+test_that(".resolveSumstatCols derives the Wald z from beta/se when z is absent", {
+  df <- .toyGwasDf(4)
+  df$beta <- df$z * 0.5                 # se = 0.5 -> beta/se reproduces z
+  df$se   <- rep(0.5, 4)
+  df$z    <- NULL                       # no z column
+  out <- pecotmr:::.resolveSumstatCols(df, NULL, "lbl")
+  expect_equal(out$Z, .toyGwasDf(4)$z)  # Z == beta/se
+  expect_true(all(c("BETA", "SE") %in% colnames(out)))  # beta/se still attached
+})
+
+test_that(".resolveSumstatCols keeps a supplied z (beta/se do not overwrite it)", {
+  df <- .toyGwasDf(4)
+  df$beta <- rep(99, 4); df$se <- rep(1, 4)   # beta/se = 99, must be ignored
+  out <- pecotmr:::.resolveSumstatCols(df, NULL, "lbl")
+  expect_equal(out$Z, .toyGwasDf(4)$z)        # z column wins, not 99
+})
+
+test_that(".resolveSumstatCols errors when neither z nor beta+se is present", {
+  df <- .toyGwasDf(3)
+  df$z <- NULL                          # no z, no beta/se
+  expect_error(pecotmr:::.resolveSumstatCols(df, NULL, "lbl"),
+               "needs a z field", fixed = TRUE)
+})
+
+test_that("loadGwasSumStatsFromManifest builds from beta/se with no z column", {
+  tmp <- withr::local_tempdir()
+  df <- .toyGwasDf(5)
+  df$beta <- df$z * 2; df$se <- rep(2, 5)     # beta/se reproduces z
+  z_expected <- df$z
+  df$z <- NULL
+  ssPath <- .writeSumstatsTsv(df, file.path(tmp, "betase.tsv"))
+  manifest <- data.frame(study = "bs", sumStatsPath = ssPath,
+                         stringsAsFactors = FALSE)
+  obj <- loadGwasSumStatsFromManifest(manifest, genome = "hg38",
+                                      ldSketch = .toyLdSketch())
+  expect_s4_class(obj, "GwasSumStats")
+  mc <- S4Vectors::mcols(obj$entry[[1L]])
+  expect_true("Z" %in% colnames(mc))
+  expect_equal(as.numeric(mc$Z), z_expected)  # derived Wald z
+})
+
 test_that(".resolveSumstatCols allows no N when allowNoN = TRUE (study scalar)", {
   df <- .toyGwasDf(3)
   df$n_sample <- NULL                  # no per-variant N, no counts