more sumstats fixes

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)
 })
 

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)
 }
 

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)
 })
 

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)
   }
 }

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)
 

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
 # ---------------------------------------------------------------------------

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"),

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),