Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 9 additions & 1 deletion R/gwasSumStats.R
Original file line number Diff line number Diff line change
Expand Up @@ -99,12 +99,19 @@ NULL
#' consumers (e.g. \code{\link{colocboostPipeline}}) use the effective
#' sample size \code{4 / (1/nCase + 1/nControl)} in place of the
#' per-variant \code{N}.
#' @param nSample Optional per-study total sample size (numeric; default
#' \code{NULL}). Attached only when supplied (length 1 or length(study)).
#' Used as the study-level fallback for the per-variant \code{N} when a study
#' has no per-variant \code{N} column and no case/control counts. Named
#' \code{nSample} to avoid clashing with \code{getNSamples()} (the LD-panel
#' sample size).
#' @param ... Additional per-study columns to attach to the collection.
#' @return A \code{GwasSumStats} object.
#' @export
GwasSumStats <- function(study, entry, genome, ldSketch = NULL,
varY = NA_real_, nCase = NULL,
nControl = NULL, qcInfo = list(), ...) {
nControl = NULL, nSample = NULL,
qcInfo = list(), ...) {
if (missing(study) || missing(entry) || missing(genome)) {
stop("`study`, `entry`, and `genome` are all required.")
}
Expand Down Expand Up @@ -139,6 +146,7 @@ GwasSumStats <- function(study, entry, genome, ldSketch = NULL,
# non-case/control studies in a mixed collection.
if (!is.null(nCase)) cols$nCase <- recycle(nCase, "nCase")
if (!is.null(nControl)) cols$nControl <- recycle(nControl, "nControl")
if (!is.null(nSample)) cols$nSample <- recycle(nSample, "nSample")
extras <- list(...)
for (nm in names(extras)) cols[[nm]] <- extras[[nm]]
df <- do.call(S4Vectors::DataFrame, c(cols, list(check.names = FALSE)))
Expand Down
62 changes: 44 additions & 18 deletions R/manifestLoaders.R
Original file line number Diff line number Diff line change
Expand Up @@ -387,10 +387,12 @@ 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 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) {
# 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.
.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
# read paths, and so explicit columnMapping values match the bare name.
Expand Down Expand Up @@ -431,9 +433,10 @@ NULL
hasN <- !is.na(nSrc) && !is.null(nSrc)
hasCounts <- !is.na(ncaseSrc) && !is.null(ncaseSrc) &&
!is.na(ncontrolSrc) && !is.null(ncontrolSrc)
if (!hasN && !hasCounts) {
if (!hasN && !hasCounts && !allowNoN) {
stop(label, ": sumstats file needs an N field (n_sample/N/n), or ",
"N_CASE + N_CONTROL columns for the effective sample size ",
"N_CASE + N_CONTROL columns for the effective sample size, or a ",
"study-level scalar (nCase/nControl or nSample) in the manifest ",
"(supply a columnMapping if the source names differ).")
}
out <- data.frame(
Expand Down Expand Up @@ -554,7 +557,8 @@ NULL

# Read one delimited-text sumstats file. Region reads only when a <path>.tbi
# sidecar exists; otherwise the whole file is read and any region is ignored.
.readSumStatsText <- function(path, region, columnMapping, label) {
.readSumStatsText <- function(path, region, columnMapping, label,
allowNoN = FALSE) {
hasTbi <- file.exists(paste0(path, ".tbi"))
raw <- if (!is.null(region) && hasTbi) {
.readTabixRegion(path, region)
Expand All @@ -570,29 +574,31 @@ NULL
col_types = readr::cols(.default = readr::col_character())),
check.names = FALSE)
}
.resolveSumstatCols(raw, columnMapping, label)
.resolveSumstatCols(raw, columnMapping, label, allowNoN = allowNoN)
}

# Dispatch a sumstats file to the text or GWAS-VCF reader. BCF is rejected.
.readSumStatsFile <- function(path, region, columnMapping, sampleSelect,
formatMapping, label) {
formatMapping, label, allowNoN = FALSE) {
lower <- tolower(path)
if (endsWith(lower, ".bcf")) {
stop(label, ": BCF sumstats are not supported; convert '", path,
"' to a bgzipped, tabix-indexed VCF.")
}
if (grepl("\\.vcf(\\.b?gz)?$", lower)) {
# GWAS-VCF always carries an SS (sample-size) FORMAT field, so the study
# scalar never has to stand in for a per-variant N here.
.readSumStatsVcf(path, region, sampleSelect, formatMapping, label)
} else {
.readSumStatsText(path, region, columnMapping, label)
.readSumStatsText(path, region, columnMapping, label, allowNoN = allowNoN)
}
}

# Read one sumstats source into an entry GRanges (canonical mcols).
.loadSumStatsEntry <- function(path, region, columnMapping, sampleSelect,
formatMapping, label) {
formatMapping, label, allowNoN = FALSE) {
df <- .readSumStatsFile(path, region, columnMapping, sampleSelect,
formatMapping, label)
formatMapping, label, allowNoN = allowNoN)
.dfToEntryGranges(df)
}

Expand Down Expand Up @@ -802,6 +808,7 @@ loadQtlDatasetFromManifest <- function(manifest, study = NULL,
"columnMappingFile"),
nCase = c("nCase", "n_case"),
nControl = c("nControl", "n_control"),
nSample = c("nSample", "n_sample"),
varY = c("varY", "var_y"),
genome = c("genome"),
ldSketchPath = c("ldSketchPath", "ld_sketch_path", "ld_sketch"))
Expand All @@ -824,8 +831,11 @@ loadQtlDatasetFromManifest <- function(manifest, study = NULL,
#' @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}
#' (optional), \code{varY} (optional), and the single-valued \code{genome} /
#' \code{ldSketchPath}.
#' (optional), \code{nSample} (optional study-level total N), \code{varY}
#' (optional), and the single-valued \code{genome} / \code{ldSketchPath}.
#' When a row supplies a study-level scalar (\code{nCase} + \code{nControl},
#' or \code{nSample}), its sumstats file need not carry a per-variant \code{N}
#' column; \code{\link{summaryStatsQc}} fills \code{N} from the scalar.
#' @param genome Genome build; reconciled with a \code{genome} column.
#' @param ldSketch A \code{\link{GenotypeHandle}} or a spec
#' (path/prefix/genoMeta); reconciled with an \code{ldSketchPath} column.
Expand Down Expand Up @@ -859,6 +869,20 @@ loadGwasSumStatsFromManifest <- function(manifest, genome = NULL,
genome <- .reconcileScalar(df$genome, genome, "genome")
ldSketch <- .resolveLdSketchInput(df, ldSketch, base)

# Study-level sample-size scalars (from the manifest). When a row carries a
# usable scalar (n_case + n_control, or n_sample), the sumstats file need not
# supply a per-variant N: summaryStatsQc fills N from the scalar. Compute the
# per-row scalar presence up front so it can gate the entry reader's N check.
nCaseCol <- if ("nCase" %in% names(df)) as.numeric(df$nCase) else NULL
nControlCol <- if ("nControl" %in% names(df)) as.numeric(df$nControl) else NULL
nSampleCol <- if ("nSample" %in% names(df)) as.numeric(df$nSample) else NULL
hasStudyScalar <- function(i) {
ccOk <- !is.null(nCaseCol) && !is.null(nControlCol) &&
is.finite(nCaseCol[[i]]) && is.finite(nControlCol[[i]])
nOk <- !is.null(nSampleCol) && is.finite(nSampleCol[[i]])
isTRUE(ccOk) || isTRUE(nOk)
}

entries <- lapply(seq_len(nrow(df)), function(i) {
label <- paste0("GwasSumStats[study=", df$study[[i]], "]")
mapping <- if ("columnMapping" %in% names(df) &&
Expand All @@ -869,16 +893,18 @@ loadGwasSumStatsFromManifest <- function(manifest, genome = NULL,
columnMapping
}
gr <- .loadSumStatsEntry(.resolveRel(as.character(df$sumStatsPath[[i]]), base),
region, mapping, sampleSelect, formatMapping, label)
region, mapping, sampleSelect, formatMapping, label,
allowNoN = hasStudyScalar(i))
.checkLdContainment(ldSketch, gr, minLdOverlapWarn, label)
gr
})

args <- list(study = as.character(df$study), entry = entries,
genome = genome, ldSketch = ldSketch)
if ("nCase" %in% names(df)) args$nCase <- as.numeric(df$nCase)
if ("nControl" %in% names(df)) args$nControl <- as.numeric(df$nControl)
if ("varY" %in% names(df)) args$varY <- as.numeric(df$varY)
if (!is.null(nCaseCol)) args$nCase <- nCaseCol
if (!is.null(nControlCol)) args$nControl <- nControlCol
if (!is.null(nSampleCol)) args$nSample <- nSampleCol
if ("varY" %in% names(df)) args$varY <- as.numeric(df$varY)
do.call(GwasSumStats, args)
}

Expand Down
72 changes: 42 additions & 30 deletions R/sumstatsQc.R
Original file line number Diff line number Diff line change
Expand Up @@ -2591,56 +2591,63 @@ krigingOutlierQc <- function(zScore, R, n, variantIds = NULL,
list(df = df, skipped = FALSE)
}

# Internal: canonicalize the working per-variant `N` to the effective sample
# size for case/control input. Counts are resolved by priority --- per-variant
# N_CASE / N_CONTROL columns first, else study-level nCase / nControl scalars.
# Returns list(df=, nSource=). nSource is one of "effective" (N set from
# counts), "column" (existing N used as-is), "total" (raw total from counts,
# escape hatch), or NA_character_ (no counts and no N). `emit` logs the
# counts-win override.
# Internal: canonicalize the working per-variant `N`. The N source is resolved
# by a four-level priority: (1) per-variant N_CASE / N_CONTROL columns, (2) study
# -level nCase / nControl scalars, (3) a per-variant N column, (4) a study-level
# nSample scalar (total N). Levels 1-2 give the effective sample size (default)
# or the raw total (escape hatch). Returns list(df=, nSource=), where nSource is
# "effective" | "column" | "total" | "study-n" | NA_character_ (no source).
# `emit` logs the counts-win override. A study-level scalar fills `df$N` even
# when the entry has no per-variant N column.
.resolveEffectiveN <- function(df, opts, emit) {
hasCols <- all(c("N_CASE", "N_CONTROL") %in% colnames(df))
hasScalar <- !is.null(opts$nCase) && !is.null(opts$nControl) &&
length(opts$nCase) == 1L && length(opts$nControl) == 1L &&
is.finite(opts$nCase) && is.finite(opts$nControl) &&
opts$nCase > 0 && opts$nControl > 0
hasNSample <- !is.null(opts$nSample) && length(opts$nSample) == 1L &&
is.finite(opts$nSample) && opts$nSample > 0
hasN <- "N" %in% colnames(df)
nRow <- nrow(df)

if (!isTRUE(opts$effectiveN)) {
# Escape hatch: use the N column as-is; when there is no N but counts are
# present, fall back to the raw total (n_case + n_control), no override.
if (!hasN && (hasCols || hasScalar)) {
if (hasCols) {
df$N <- as.numeric(df$N_CASE) + as.numeric(df$N_CONTROL)
} else {
df$N <- rep(opts$nCase + opts$nControl, nRow)
}
# Escape hatch: prefer the raw N column; else raw total from counts; else the
# study-level total N. No effective-N conversion, no override.
if (hasN) return(list(df = df, nSource = "column"))
if (hasCols) {
df$N <- as.numeric(df$N_CASE) + as.numeric(df$N_CONTROL)
return(list(df = df, nSource = "total"))
}
if (hasScalar) {
df$N <- rep(opts$nCase + opts$nControl, nRow)
return(list(df = df, nSource = "total"))
}
return(list(df = df, nSource = if (hasN) "column" else NA_character_))
if (hasNSample) {
df$N <- rep(opts$nSample, nRow)
return(list(df = df, nSource = "study-n"))
}
return(list(df = df, nSource = NA_character_))
}

# Default on: derive the effective sample size from the counts when present.
# Default: per-variant c/c -> study c/c -> per-variant N -> study nSample.
if (hasCols) {
neff <- effectiveN(df$N_CASE, df$N_CONTROL)
if (hasN) {
emit("QC track: N overridden by effective N from per-variant ",
"n_case/n_control.")
}
df$N <- neff
if (hasN) emit("QC track: N overridden by effective N from per-variant ",
"n_case/n_control.")
df$N <- effectiveN(df$N_CASE, df$N_CONTROL)
return(list(df = df, nSource = "effective"))
}
if (hasScalar) {
neff <- rep(effectiveN(opts$nCase, opts$nControl), nRow)
if (hasN) {
emit("QC track: N overridden by effective N from study ",
"nCase/nControl.")
}
df$N <- neff
if (hasN) emit("QC track: N overridden by effective N from study ",
"nCase/nControl.")
df$N <- rep(effectiveN(opts$nCase, opts$nControl), nRow)
return(list(df = df, nSource = "effective"))
}
list(df = df, nSource = if (hasN) "column" else NA_character_)
if (hasN) return(list(df = df, nSource = "column"))
if (hasNSample) {
df$N <- rep(opts$nSample, nRow)
return(list(df = df, nSource = "study-n"))
}
list(df = df, nSource = NA_character_)
}

# Internal: per-entry pipeline. Returns the cleaned GRanges and an audit list.
Expand Down Expand Up @@ -3158,6 +3165,9 @@ summaryStatsQc <- function(sumstats,
# no per-variant N_CASE / N_CONTROL columns). NULL for quantitative traits.
opts$nCase <- if ("nCase" %in% names(sumstats)) as.numeric(sumstats$nCase)[[i]] else NULL
opts$nControl <- if ("nControl" %in% names(sumstats)) as.numeric(sumstats$nControl)[[i]] else NULL
# Study-level total N: the level-4 fallback used when the entry has no
# per-variant N_CASE/N_CONTROL, no study nCase/nControl, and no N column.
opts$nSample <- if ("nSample" %in% names(sumstats)) as.numeric(sumstats$nSample)[[i]] else NULL
# Per-entry label woven into QC log messages and the rollup. For
# QtlSumStats it's (study/context/trait); for GwasSumStats it's the
# study identifier.
Expand Down Expand Up @@ -3214,6 +3224,8 @@ summaryStatsQc <- function(sumstats,
as.numeric(sumstats$nCase) else NULL,
nControl = if ("nControl" %in% names(sumstats))
as.numeric(sumstats$nControl) else NULL,
nSample = if ("nSample" %in% names(sumstats))
as.numeric(sumstats$nSample) else NULL,
qcInfo = qcInfo)
} else {
QtlSumStats(
Expand Down
7 changes: 7 additions & 0 deletions man/GwasSumStats.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

7 changes: 5 additions & 2 deletions man/loadGwasSumStatsFromManifest.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

49 changes: 49 additions & 0 deletions tests/testthat/test_manifestLoaders.R
Original file line number Diff line number Diff line change
Expand Up @@ -245,6 +245,55 @@ test_that(".resolveSumstatCols requires an N field or N_CASE + N_CONTROL", {
"N_CASE + N_CONTROL", fixed = TRUE)
})

test_that(".resolveSumstatCols allows no N when allowNoN = TRUE (study scalar)", {
df <- .toyGwasDf(3)
df$n_sample <- NULL # no per-variant N, no counts
out <- pecotmr:::.resolveSumstatCols(df, NULL, "lbl", allowNoN = TRUE)
expect_false("N" %in% colnames(out)) # built without N ...
expect_true(all(c("SNP", "A1", "A2", "Z") %in% colnames(out))) # ... core kept
})

test_that("loadGwasSumStatsFromManifest builds from a study nSample scalar (no per-variant N)", {
tmp <- withr::local_tempdir()
df <- .toyGwasDf(5)
df$n_sample <- NULL # sumstats has no per-variant N and no counts
ssPath <- .writeSumstatsTsv(df, file.path(tmp, "studyn.tsv"))
manifest <- data.frame(study = "sn", sumStatsPath = ssPath,
n_sample = 487511, stringsAsFactors = FALSE)
obj <- loadGwasSumStatsFromManifest(manifest, genome = "hg38",
ldSketch = .toyLdSketch())
expect_s4_class(obj, "GwasSumStats")
expect_false("N" %in% colnames(S4Vectors::mcols(obj$entry[[1L]])))
expect_equal(as.numeric(obj$nSample), 487511) # forwarded to the slot
})

test_that("loadGwasSumStatsFromManifest builds from study case/control scalars (no per-variant N)", {
tmp <- withr::local_tempdir()
df <- .toyGwasDf(5)
df$n_sample <- NULL
ssPath <- .writeSumstatsTsv(df, file.path(tmp, "studycc.tsv"))
manifest <- data.frame(study = "cc", sumStatsPath = ssPath,
n_case = 5000, n_control = 15000, stringsAsFactors = FALSE)
obj <- loadGwasSumStatsFromManifest(manifest, genome = "hg38",
ldSketch = .toyLdSketch())
expect_s4_class(obj, "GwasSumStats")
expect_false("N" %in% colnames(S4Vectors::mcols(obj$entry[[1L]])))
expect_equal(as.numeric(obj$nCase), 5000)
expect_equal(as.numeric(obj$nControl), 15000)
})

test_that("loadGwasSumStatsFromManifest still errors when a study has no N source at all", {
tmp <- withr::local_tempdir()
df <- .toyGwasDf(5)
df$n_sample <- NULL # no per-variant N, no counts ...
ssPath <- .writeSumstatsTsv(df, file.path(tmp, "non.tsv"))
manifest <- data.frame(study = "x", sumStatsPath = ssPath,
stringsAsFactors = FALSE) # ... and no study scalar either
expect_error(loadGwasSumStatsFromManifest(manifest, genome = "hg38",
ldSketch = .toyLdSketch()),
"needs an N field")
})

test_that(".readColumnMapping accepts a named list and a YAML file alike", {
tmp <- withr::local_tempdir()
mp <- file.path(tmp, "m.yaml")
Expand Down
Loading
Loading