diff --git a/R/manifestLoaders.R b/R/manifestLoaders.R index 80e83c93..380d9905 100644 --- a/R/manifestLoaders.R +++ b/R/manifestLoaders.R @@ -330,6 +330,8 @@ NULL A2 = c("A2", "a2"), Z = c("z", "Z"), N = c("n_sample", "N", "n"), + N_CASE = c("N_CASE", "n_case", "ncase"), + N_CONTROL = c("N_CONTROL", "n_control", "ncontrol"), BETA = c("beta", "BETA"), SE = c("se", "SE"), P = c("p", "P", "pvalue", "pval"), @@ -347,6 +349,8 @@ NULL A2 = c("A2", "a2"), Z = c("z", "Z"), N = c("n_sample", "N", "n"), + N_CASE = c("n_case", "N_CASE", "ncase"), + N_CONTROL = c("n_control", "N_CONTROL", "ncontrol"), BETA = c("beta", "BETA"), SE = c("se", "SE"), P = c("p", "P", "pvalue"), @@ -382,7 +386,9 @@ NULL # Standardise the columns of a raw sumstats data.frame into the canonical # schema expected by .dfToEntryGranges: chrom, pos, SNP, A1, A2, Z, N (+ -# optional BETA, SE, P, MAF, INFO). `mapping` (named std -> source) overrides +# optional N_CASE, N_CONTROL, BETA, SE, P, MAF, INFO). N is required UNLESS +# per-variant N_CASE + N_CONTROL are both present, from which summaryStatsQc +# derives the effective sample size. `mapping` (named std -> source) overrides # the default aliases per key. .resolveSumstatCols <- function(df, columnMapping, label) { # Strip a leading '#' from the first column name so a '#CHR'-style header @@ -407,7 +413,7 @@ NULL } intersect(.sumstatColumnAliases[[key]], names(df))[1L] } - required <- c("chrom", "pos", "variant_id", "A1", "A2", "Z", "N") + required <- c("chrom", "pos", "variant_id", "A1", "A2", "Z") resolved <- setNames(lapply(required, resolveKey), required) missingKeys <- required[vapply(resolved, function(x) is.na(x) || is.null(x), logical(1))] if (length(missingKeys) > 0L) { @@ -415,6 +421,21 @@ NULL paste(missingKeys, collapse = ", "), " (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 + # downstream when effectiveN is on). + nSrc <- resolveKey("N") + ncaseSrc <- resolveKey("N_CASE") + ncontrolSrc <- resolveKey("N_CONTROL") + hasN <- !is.na(nSrc) && !is.null(nSrc) + hasCounts <- !is.na(ncaseSrc) && !is.null(ncaseSrc) && + !is.na(ncontrolSrc) && !is.null(ncontrolSrc) + if (!hasN && !hasCounts) { + stop(label, ": sumstats file needs an N field (n_sample/N/n), or ", + "N_CASE + N_CONTROL columns for the effective sample size ", + "(supply a columnMapping if the source names differ).") + } out <- data.frame( chrom = as.character(df[[resolved$chrom]]), pos = as.integer(df[[resolved$pos]]), @@ -422,8 +443,14 @@ NULL A1 = as.character(df[[resolved$A1]]), A2 = as.character(df[[resolved$A2]]), Z = as.numeric(df[[resolved$Z]]), - N = as.numeric(df[[resolved$N]]), stringsAsFactors = FALSE) + if (hasN) out$N <- as.numeric(df[[nSrc]]) + # Per-variant case/control counts -> N_CASE / N_CONTROL (kept by + # .dfToEntryGranges; summaryStatsQc(effectiveN=) consumes them). + if (hasCounts) { + out$N_CASE <- as.numeric(df[[ncaseSrc]]) + out$N_CONTROL <- as.numeric(df[[ncontrolSrc]]) + } for (key in c("BETA", "SE", "P", "MAF", "INFO")) { src <- resolveKey(key) if (!is.na(src) && !is.null(src)) out[[key]] <- as.numeric(df[[src]]) diff --git a/tests/testthat/test_manifestLoaders.R b/tests/testthat/test_manifestLoaders.R index 65647861..93d50ec2 100644 --- a/tests/testthat/test_manifestLoaders.R +++ b/tests/testthat/test_manifestLoaders.R @@ -216,6 +216,35 @@ test_that("loadGwasSumStatsFromManifest resolves a YAML column mapping", { c("G", "C", "G", "C", "A")) }) +# --------------------------------------------------------------------------- +# Case/control counts -> effective sample size +# --------------------------------------------------------------------------- + +test_that("loadGwasSumStatsFromManifest reads per-variant N_CASE/N_CONTROL (no N)", { + tmp <- withr::local_tempdir() + df <- .toyGwasDf(5) + df$n_sample <- NULL # no N column ... + df$n_case <- rep(1000L, 5) # ... only per-variant case/control counts + df$n_control <- rep(2000L, 5) + ssPath <- .writeSumstatsTsv(df, file.path(tmp, "cc.tsv")) + manifest <- data.frame(study = "cc", sumStatsPath = ssPath, + stringsAsFactors = FALSE) + obj <- loadGwasSumStatsFromManifest(manifest, genome = "hg38", + ldSketch = .toyLdSketch()) + expect_s4_class(obj, "GwasSumStats") + mc <- S4Vectors::mcols(obj$entry[[1L]]) + expect_true(all(c("N_CASE", "N_CONTROL") %in% colnames(mc))) + expect_false("N" %in% colnames(mc)) # effective N is derived later by summaryStatsQc + expect_equal(as.integer(mc$N_CASE), rep(1000L, 5)) +}) + +test_that(".resolveSumstatCols requires an N field or N_CASE + N_CONTROL", { + df <- .toyGwasDf(3) + df$n_sample <- NULL # drop the only sample-size source + expect_error(pecotmr:::.resolveSumstatCols(df, NULL, "lbl"), + "N_CASE + N_CONTROL", fixed = TRUE) +}) + test_that(".readColumnMapping accepts a named list and a YAML file alike", { tmp <- withr::local_tempdir() mp <- file.path(tmp, "m.yaml")