diff --git a/R/gwasSumStats.R b/R/gwasSumStats.R index f9a03ecd..acbb84c7 100644 --- a/R/gwasSumStats.R +++ b/R/gwasSumStats.R @@ -233,6 +233,11 @@ setMethod("subsetChr", "GwasSumStats", function(x, chr) { genome = x@genome, ldSketch = x@ldSketch, varY = as.numeric(x$varY), + # Preserve the optional per-study case/control counts + total N through the + # chromosome subset (they are study-level scalars, not per-variant). + nCase = if ("nCase" %in% names(x)) as.numeric(x$nCase) else NULL, + nControl = if ("nControl" %in% names(x)) as.numeric(x$nControl) else NULL, + nSample = if ("nSample" %in% names(x)) as.numeric(x$nSample) else NULL, qcInfo = x@qcInfo) }) diff --git a/R/manifestLoaders.R b/R/manifestLoaders.R index a9565625..e8192158 100644 --- a/R/manifestLoaders.R +++ b/R/manifestLoaders.R @@ -838,6 +838,7 @@ loadQtlDatasetFromManifest <- function(manifest, study = NULL, "filePath"), columnMapping = c("columnMapping", "column_mapping", "column_mapping_file", "columnMappingFile"), + nSample = c("nSample", "n_sample"), varY = c("varY", "var_y"), genome = c("genome"), ldSketchPath = c("ldSketchPath", "ld_sketch_path", "ld_sketch")) @@ -928,10 +929,14 @@ loadGwasSumStatsFromManifest <- function(manifest, genome = NULL, do.call(GwasSumStats, args) } -# Build the entry list + tuple vectors for a QtlSumStats manifest. +# Build the entry list + tuple vectors for a QtlSumStats manifest. `allowNoN` +# is an optional per-row logical vector: when TRUE for a row, its sumstats file +# need not carry a per-variant N (a study-level nSample scalar fills it later in +# summaryStatsQc). NULL (default) requires a per-variant N on every row. .loadQtlSumStatsEntries <- function(df, base, region, ldSketch, minLdOverlapWarn, columnMapping, - sampleSelect, formatMapping) { + sampleSelect, formatMapping, + allowNoN = NULL) { lapply(seq_len(nrow(df)), function(i) { label <- paste0("QtlSumStats[", df$study[[i]], "/", df$context[[i]], "/", df$trait[[i]], "]") @@ -943,7 +948,8 @@ 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 = !is.null(allowNoN) && isTRUE(allowNoN[[i]])) if (!is.null(ldSketch)) { .checkLdContainment(ldSketch, gr, minLdOverlapWarn, label) } @@ -953,11 +959,19 @@ loadGwasSumStatsFromManifest <- function(manifest, genome = NULL, #' @title Load a QtlSumStats collection from a manifest #' @description Build a \code{\link{QtlSumStats}} from a manifest with one row -#' per \code{(study, context, trait)} tuple. No QC is run. +#' per \code{(study, context, trait)} tuple. No QC is run. 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}, \code{context}, \code{trait} (required), \code{sumStatsPath} -#' (required), \code{columnMapping} (optional), \code{varY} (optional), and -#' the single-valued \code{genome} / \code{ldSketchPath}. +#' (required), \code{columnMapping} (optional), \code{nSample} (optional +#' tuple-level total N), \code{varY} (optional), and the single-valued +#' \code{genome} / \code{ldSketchPath}. When a row supplies \code{nSample}, +#' its sumstats file need not carry a per-variant \code{N} column; +#' \code{\link{summaryStatsQc}} fills \code{N} from the scalar. (Unlike the +#' GWAS loader, there are no \code{nCase}/\code{nControl} columns: molecular +#' QTL traits are quantitative.) #' @param genome Genome build; reconciled with a \code{genome} column. #' @param ldSketch A \code{\link{GenotypeHandle}} or spec; reconciled with an #' \code{ldSketchPath} column. @@ -979,14 +993,23 @@ loadQtlSumStatsFromManifest <- function(manifest, genome = NULL, label = "QtlSumStats") genome <- .reconcileScalar(df$genome, genome, "genome") ldSketch <- .resolveLdSketchInput(df, ldSketch, base) + + # Tuple-level total-N scalar (from the manifest). When a row carries a usable + # nSample the sumstats file need not supply a per-variant N: summaryStatsQc + # fills N from the scalar. Gate the entry reader's N check per row on it. + nSampleCol <- if ("nSample" %in% names(df)) as.numeric(df$nSample) else NULL + allowNoN <- if (!is.null(nSampleCol)) is.finite(nSampleCol) else NULL + entries <- .loadQtlSumStatsEntries(df, base, region, ldSketch, minLdOverlapWarn, columnMapping, - sampleSelect, formatMapping) + sampleSelect, formatMapping, + allowNoN = allowNoN) args <- list(study = as.character(df$study), context = as.character(df$context), trait = as.character(df$trait), entry = entries, genome = genome, ldSketch = ldSketch) - if ("varY" %in% names(df)) args$varY <- as.numeric(df$varY) + if (!is.null(nSampleCol)) args$nSample <- nSampleCol + if ("varY" %in% names(df)) args$varY <- as.numeric(df$varY) do.call(QtlSumStats, args) } diff --git a/R/qtlSumStats.R b/R/qtlSumStats.R index 45b300ec..87a326f7 100644 --- a/R/qtlSumStats.R +++ b/R/qtlSumStats.R @@ -89,11 +89,19 @@ NULL #' @param ldSketch A \code{GenotypeHandle} carrying the LD reference. #' @param varY Optional numeric vector of per-tuple phenotype variances #' (\code{NA_real_} entries allowed). +#' @param nSample Optional per-tuple 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 tuple +#' has no per-variant \code{N} column. Named \code{nSample} to avoid clashing +#' with \code{getNSamples()} (the LD-panel sample size). Unlike GWAS, QTL +#' collections carry no case/control counts (molecular traits are +#' quantitative), so only this total-N fallback is exposed. #' @param ... Additional per-tuple columns to attach to the collection. #' @return A \code{QtlSumStats} object. #' @export QtlSumStats <- function(study, context, trait, entry, genome, ldSketch = NULL, - varY = NA_real_, qcInfo = list(), traitPos = NULL, ...) { + varY = NA_real_, nSample = NULL, qcInfo = list(), + traitPos = NULL, ...) { if (missing(study) || missing(context) || missing(trait) || missing(entry) || missing(genome)) { stop("`study`, `context`, `trait`, `entry`, and `genome` ", @@ -131,6 +139,17 @@ QtlSumStats <- function(study, context, trait, entry, genome, ldSketch = NULL, entry = S4Vectors::SimpleList(entry), varY = as.numeric(varY) ) + # nSample is OPTIONAL (per-tuple study-level total sample size): the column is + # attached only when supplied (default NULL), so QTL collections without it + # keep the original schema. summaryStatsQc() fills a missing per-variant N + # from this scalar. Provide length 1 or length(study). + if (!is.null(nSample)) { + if (length(nSample) == 1L && n > 1L) nSample <- rep(nSample, n) + if (length(nSample) != n) { + stop("`nSample` must have length 1 or length(study).") + } + cols$nSample <- as.numeric(nSample) + } # Optional trait-position provenance (one GRanges per trait; TSS = start()), # carried forward into QtlFineMappingResult / TwasWeights since the true trait # position cannot be inferred from summary statistics alone. @@ -273,6 +292,7 @@ setMethod("subsetChr", "QtlSumStats", function(x, chr) { genome = x@genome, ldSketch = x@ldSketch, varY = as.numeric(x$varY), + nSample = if ("nSample" %in% names(x)) as.numeric(x$nSample) else NULL, qcInfo = x@qcInfo) }) diff --git a/R/sumstatsQc.R b/R/sumstatsQc.R index 04260d26..45550542 100644 --- a/R/sumstatsQc.R +++ b/R/sumstatsQc.R @@ -3236,6 +3236,9 @@ summaryStatsQc <- function(sumstats, genome = getGenome(sumstats), ldSketch = getLdSketch(sumstats), varY = as.numeric(sumstats$varY), + # Preserve the optional per-tuple study-level total N through QC. + nSample = if ("nSample" %in% names(sumstats)) + as.numeric(sumstats$nSample) else NULL, qcInfo = qcInfo) } } diff --git a/man/GwasSumStats.Rd b/man/GwasSumStats.Rd index 4fa9c7b3..a719165a 100644 --- a/man/GwasSumStats.Rd +++ b/man/GwasSumStats.Rd @@ -43,10 +43,11 @@ sample size \code{4 / (1/nCase + 1/nControl)} in place of the per-variant \code{N}.} \item{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).} +\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).} \item{...}{Additional per-study columns to attach to the collection.} } diff --git a/man/QtlSumStats.Rd b/man/QtlSumStats.Rd index 898ab7eb..69360714 100644 --- a/man/QtlSumStats.Rd +++ b/man/QtlSumStats.Rd @@ -12,6 +12,7 @@ QtlSumStats( genome, ldSketch = NULL, varY = NA_real_, + nSample = NULL, qcInfo = list(), traitPos = NULL, ... @@ -36,6 +37,14 @@ because all entries share the same LD sketch.} \item{varY}{Optional numeric vector of per-tuple phenotype variances (\code{NA_real_} entries allowed).} +\item{nSample}{Optional per-tuple 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 tuple +has no per-variant \code{N} column. Named \code{nSample} to avoid clashing +with \code{getNSamples()} (the LD-panel sample size). Unlike GWAS, QTL +collections carry no case/control counts (molecular traits are +quantitative), so only this total-N fallback is exposed.} + \item{...}{Additional per-tuple columns to attach to the collection.} } \value{ diff --git a/man/loadQtlSumStatsFromManifest.Rd b/man/loadQtlSumStatsFromManifest.Rd index 259b8a7b..428d0c95 100644 --- a/man/loadQtlSumStatsFromManifest.Rd +++ b/man/loadQtlSumStatsFromManifest.Rd @@ -18,8 +18,13 @@ loadQtlSumStatsFromManifest( \arguments{ \item{manifest}{A data.frame or path. Columns (snake_case aliases accepted): \code{study}, \code{context}, \code{trait} (required), \code{sumStatsPath} -(required), \code{columnMapping} (optional), \code{varY} (optional), and -the single-valued \code{genome} / \code{ldSketchPath}.} +(required), \code{columnMapping} (optional), \code{nSample} (optional +tuple-level total N), \code{varY} (optional), and the single-valued +\code{genome} / \code{ldSketchPath}. When a row supplies \code{nSample}, +its sumstats file need not carry a per-variant \code{N} column; +\code{\link{summaryStatsQc}} fills \code{N} from the scalar. (Unlike the +GWAS loader, there are no \code{nCase}/\code{nControl} columns: molecular +QTL traits are quantitative.)} \item{genome}{Genome build; reconciled with a \code{genome} column.} @@ -34,5 +39,8 @@ A \code{QtlSumStats} object. } \description{ Build a \code{\link{QtlSumStats}} from a manifest with one row - per \code{(study, context, trait)} tuple. No QC is run. + per \code{(study, context, trait)} tuple. No QC is run. 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_GwasSumStats.R b/tests/testthat/test_GwasSumStats.R index a3058353..66966bd7 100644 --- a/tests/testthat/test_GwasSumStats.R +++ b/tests/testthat/test_GwasSumStats.R @@ -127,6 +127,28 @@ test_that("subsetchr() filters correctly", { }) +test_that("subsetChr() preserves study-level nCase/nControl/nSample scalars", { + # Regression: the chromosome subset rebuilds the GwasSumStats and previously + # dropped the optional per-study case/control counts + total N. + gr <- GenomicRanges::GRanges( + c("chr1", "chr1", "chr2"), + IRanges::IRanges(start = c(100L, 200L, 300L), width = 1L)) + S4Vectors::mcols(gr) <- S4Vectors::DataFrame( + SNP = paste0("rs", 1:3), A1 = rep("A", 3), A2 = rep("G", 3), + Z = c(1.0, -0.5, 2.0), N = rep(100L, 3)) + obj <- GwasSumStats( + study = "g1", entry = list(gr), genome = "hg19", + ldSketch = .sh_makeGenotypeHandle(), + nCase = 5000, nControl = 15000, nSample = 20000) + + sub <- subsetChr(obj, "chr1") + expect_equal(nSnps(sub), 2) + expect_equal(as.numeric(sub$nCase), 5000) + expect_equal(as.numeric(sub$nControl), 15000) + expect_equal(as.numeric(sub$nSample), 20000) +}) + + test_that("getvary() returns var_y and NULL cases", { df <- make_test_sumstats_df(5) diff --git a/tests/testthat/test_manifestLoaders.R b/tests/testthat/test_manifestLoaders.R index e14e5f36..98abd9e3 100644 --- a/tests/testthat/test_manifestLoaders.R +++ b/tests/testthat/test_manifestLoaders.R @@ -364,6 +364,52 @@ test_that("loadQtlSumStatsFromManifest builds a per-tuple collection", { expect_equal(as.character(obj$context), "Whole_Blood") }) +test_that("loadQtlSumStatsFromManifest derives the Wald z from beta/se (no z column)", { + tmp <- withr::local_tempdir() + df <- .toyGwasDf(5) + z_expected <- df$z + df$beta <- df$z * 2; df$se <- rep(2, 5) # beta/se reproduces z + df$z <- NULL # ... but no z column + ssPath <- .writeSumstatsTsv(df, file.path(tmp, "geneA.tsv")) + manifest <- data.frame(study = "eqtl", context = "Whole_Blood", + trait = "geneA", sumStatsPath = ssPath, + stringsAsFactors = FALSE) + obj <- loadQtlSumStatsFromManifest(manifest, genome = "hg38", + ldSketch = .toyLdSketch()) + expect_s4_class(obj, "QtlSumStats") + mc <- S4Vectors::mcols(obj$entry[[1L]]) + expect_equal(as.numeric(mc$Z), z_expected) # derived Wald z + expect_true(all(c("BETA", "SE") %in% colnames(mc))) # beta/se still attached +}) + +test_that("loadQtlSumStatsFromManifest builds from a tuple nSample scalar (no per-variant N)", { + tmp <- withr::local_tempdir() + df <- .toyGwasDf(5) + df$n_sample <- NULL # sumstats has no per-variant N + ssPath <- .writeSumstatsTsv(df, file.path(tmp, "geneA.tsv")) + manifest <- data.frame(study = "eqtl", context = "Whole_Blood", + trait = "geneA", sumStatsPath = ssPath, + n_sample = 838, stringsAsFactors = FALSE) + obj <- loadQtlSumStatsFromManifest(manifest, genome = "hg38", + ldSketch = .toyLdSketch()) + expect_s4_class(obj, "QtlSumStats") + expect_false("N" %in% colnames(S4Vectors::mcols(obj$entry[[1L]]))) + expect_equal(as.numeric(obj$nSample), 838) # forwarded to the slot +}) + +test_that("loadQtlSumStatsFromManifest still errors when a tuple has no N source at all", { + tmp <- withr::local_tempdir() + df <- .toyGwasDf(5) + df$n_sample <- NULL # no per-variant N ... + ssPath <- .writeSumstatsTsv(df, file.path(tmp, "geneA.tsv")) + manifest <- data.frame(study = "eqtl", context = "Whole_Blood", + trait = "geneA", sumStatsPath = ssPath, + stringsAsFactors = FALSE) # ... and no nSample scalar + expect_error(loadQtlSumStatsFromManifest(manifest, genome = "hg38", + ldSketch = .toyLdSketch()), + "needs an N field") +}) + # --------------------------------------------------------------------------- # QtlDataset loader # --------------------------------------------------------------------------- diff --git a/tests/testthat/test_qtlSumStats.R b/tests/testthat/test_qtlSumStats.R index 6cce50a1..7de3aeff 100644 --- a/tests/testthat/test_qtlSumStats.R +++ b/tests/testthat/test_qtlSumStats.R @@ -190,6 +190,45 @@ test_that("QtlSumStats: errors when varY length is neither 1 nor n", { ) }) +test_that("QtlSumStats: optional nSample column is absent by default, attached and recycled when supplied", { + bare <- .qtlMakeOne() + expect_false("nSample" %in% names(bare)) # schema unchanged by default + obj <- QtlSumStats( + study = c("s1", "s2"), + context = c("c1", "c1"), + trait = c("t1", "t1"), + entry = list(.qtlMakeEntryGr(2), .qtlMakeEntryGr(2)), + genome = "hg19", + ldSketch = .qtlMakeGenotypeHandle(), + nSample = 838) # scalar recycled to both tuples + expect_equal(as.numeric(obj$nSample), c(838, 838)) + expect_true(methods::validObject(obj)) +}) + +test_that("QtlSumStats: errors when nSample length is neither 1 nor n", { + expect_error( + QtlSumStats( + study = c("s1", "s2"), + context = c("c1", "c1"), + trait = c("t1", "t1"), + entry = list(.qtlMakeEntryGr(2), .qtlMakeEntryGr(2)), + genome = "hg19", + ldSketch = .qtlMakeGenotypeHandle(), + nSample = c(1, 2, 3)), + "length 1 or length\\(study\\)" + ) +}) + +test_that("QtlSumStats: subsetChr preserves the nSample column", { + gr <- .qtlMakeEntryGr(4, chr = "chr1") + obj <- QtlSumStats( + study = "s1", context = "c1", trait = "t1", + entry = list(gr), genome = "hg19", + ldSketch = .qtlMakeGenotypeHandle(), nSample = 838) + sub <- subsetChr(obj, "chr1") + expect_equal(as.numeric(sub$nSample), 838) +}) + test_that("QtlSumStats: accepts and stores extra per-tuple columns via ...", { obj <- QtlSumStats( study = c("s1", "s2"), diff --git a/tests/testthat/test_sumstatsQc.R b/tests/testthat/test_sumstatsQc.R index 2d155db3..947828c0 100644 --- a/tests/testthat/test_sumstatsQc.R +++ b/tests/testthat/test_sumstatsQc.R @@ -2912,6 +2912,24 @@ test_that("summaryStatsQc: study nSample is the level-4 fallback (no counts, no expect_identical(getQcInfo(res)$entryAudit[[1L]]$nSource, "study-n") }) +test_that("summaryStatsQc: QtlSumStats tuple nSample is the level-4 fallback too", { + # Parity with the GWAS study-n fallback: a QtlSumStats carrying only a + # tuple-level nSample (no per-variant N) fills N from the scalar and is + # preserved on the QC'd object. + gr <- .ssQ_makeEntryGr() + mc <- S4Vectors::mcols(gr) + mc$N <- NULL # remove per-variant N + S4Vectors::mcols(gr) <- mc + ss <- QtlSumStats(study = "s1", context = "c1", trait = "t1", + entry = list(gr), genome = "hg19", + ldSketch = .ssQ_makeHandle(), nSample = 838) + res <- summaryStatsQc(ss) + expect_s4_class(res, "QtlSumStats") + expect_true(all(.ssQ_entryNByPos(res$entry[[1L]]) == 838)) + expect_identical(getQcInfo(res)$entryAudit[[1L]]$nSource, "study-n") + expect_equal(as.numeric(res$nSample), 838) # slot preserved through QC +}) + test_that("summaryStatsQc: level precedence -- per-variant counts beat nSample; N column beats nSample", { # per-variant counts present alongside nSample -> counts win (effective). gr1 <- .ssQ_makeCCEntry(nCase = c(100, 200, 150, 250),