diff --git a/R/manifestLoaders.R b/R/manifestLoaders.R index b1439f33..a9565625 100644 --- a/R/manifestLoaders.R +++ b/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,7 +416,7 @@ 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) { @@ -423,6 +424,20 @@ NULL 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} diff --git a/man/loadGwasSumStatsFromManifest.Rd b/man/loadGwasSumStatsFromManifest.Rd index 700fedae..5852e6d8 100644 --- a/man/loadGwasSumStatsFromManifest.Rd +++ b/man/loadGwasSumStatsFromManifest.Rd @@ -51,5 +51,8 @@ A \code{GwasSumStats} object. } \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). } diff --git a/tests/testthat/test_fineMappingWrappers.R b/tests/testthat/test_fineMappingWrappers.R index 278e1200..8529f565 100644 --- a/tests/testthat/test_fineMappingWrappers.R +++ b/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", ...), diff --git a/tests/testthat/test_manifestLoaders.R b/tests/testthat/test_manifestLoaders.R index a2cad765..e14e5f36 100644 --- a/tests/testthat/test_manifestLoaders.R +++ b/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