diff --git a/NAMESPACE b/NAMESPACE index c0f1ef68..736321cc 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -67,6 +67,7 @@ export(dprAdaptiveGibbsWeights) export(dprGibbsWeights) export(dprVbWeights) export(dprWeights) +export(effectiveN) export(enetWeights) export(enforceDesignFullRank) export(enlocPipeline) diff --git a/R/colocboostPipeline.R b/R/colocboostPipeline.R index 30f70a6e..fb34a93b 100644 --- a/R/colocboostPipeline.R +++ b/R/colocboostPipeline.R @@ -349,7 +349,7 @@ setGeneric("colocboostPipeline", okCC <- !is.null(nCase) && !is.null(nControl) && !is.na(nCase) && !is.na(nControl) && nCase > 0 && nControl > 0 - nVal <- if (okCC) 4 / (1 / nCase + 1 / nControl) + nVal <- if (okCC) effectiveN(nCase, nControl) else if (!is.null(df$N)) df$N else NA_real_ ss <- data.frame( z = df$z, diff --git a/R/sumstatsQc.R b/R/sumstatsQc.R index f71f72dc..f86725f0 100644 --- a/R/sumstatsQc.R +++ b/R/sumstatsQc.R @@ -1996,6 +1996,33 @@ ldMismatchQc <- function(zScore, R = NULL, X = NULL, nSample = NULL, match.arg(zMismatchQc, c("none", "slalom", "dentist")) } +#' Effective sample size for a case/control study +#' +#' Computes the effective sample size +#' \code{N_eff = 4 / (1/nCase + 1/nControl) = 4 * nCase * nControl / +#' (nCase + nControl)} for case/control GWAS. Balanced studies +#' (\code{nCase == nControl}) recover the total \code{nCase + nControl}; +#' imbalanced studies give a smaller value, which is the statistically +#' correct sample size for the RSS likelihood, residual-variance estimation, +#' kriging, and the N-cutoff filter. Vectorized over \code{nCase} / +#' \code{nControl}; entries where either count is \code{NA} or \code{<= 0} +#' return \code{NA_real_}. +#' +#' @param nCase Numeric vector of case counts. +#' @param nControl Numeric vector of control counts. +#' @return Numeric vector of effective sample sizes (\code{NA_real_} where a +#' count is missing or non-positive). +#' @references Privé et al., "Identifying and correcting for misspecifications +#' in GWAS summary statistics and polygenic scores", HGG Advances 2022. +#' @export +effectiveN <- function(nCase, nControl) { + nCase <- as.numeric(nCase) + nControl <- as.numeric(nControl) + out <- 4 / (1 / nCase + 1 / nControl) + out[is.na(nCase) | is.na(nControl) | nCase <= 0 | nControl <= 0] <- NA_real_ + out +} + #' Kriging-style LD-consistency outlier QC #' #' Flags variants whose observed z-score is inconsistent with the value @@ -2114,7 +2141,7 @@ krigingOutlierQc <- function(zScore, R, n, variantIds = NULL, ranges = IRanges::IRanges(start = as.integer(df$pos), width = 1L)) if (!is.null(df$variant_id) && is.null(df$SNP)) df$SNP <- df$variant_id baseCols <- c("SNP", "A1", "A2", "Z", "N") - optCols <- c("MAF", "INFO", "BETA", "SE", "P") + optCols <- c("MAF", "INFO", "BETA", "SE", "P", "N_CASE", "N_CONTROL") use <- intersect(c(baseCols, optCols), colnames(df)) S4Vectors::mcols(gr) <- S4Vectors::DataFrame(df[, use, drop = FALSE]) gr @@ -2564,6 +2591,58 @@ 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. +.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 + 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) + } + return(list(df = df, nSource = "total")) + } + return(list(df = df, nSource = if (hasN) "column" else NA_character_)) + } + + # Default on: derive the effective sample size from the counts when present. + 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 + 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 + return(list(df = df, nSource = "effective")) + } + list(df = df, nSource = if (hasN) "column" else NA_character_) +} + # Internal: per-entry pipeline. Returns the cleaned GRanges and an audit list. .runEntrySummaryStatsQc <- function(gr, ldSketch, refGenome, opts, entryLabel = NULL) { @@ -2608,6 +2687,18 @@ krigingOutlierQc <- function(zScore, R, n, variantIds = NULL, " variant(s).") } + # Canonicalize N to the effective sample size for case/control input, + # BEFORE the N-cutoff filter so the filter, kriging (median df$N), and the + # rebuilt entry all consume N_eff. No-op for quantitative traits (no counts). + nRes <- .resolveEffectiveN(df, opts, emit) + df <- nRes$df + entryAudit$nSource <- nRes$nSource + # Keep the PIP-screen / kriging sample size consistent with the canonicalized + # N: when effective (or total) N was applied, nForPip (set upstream from the + # entry's raw N) must follow the same N the N-filter/fit use. + if (isTRUE(nRes$nSource %in% c("effective", "total")) && "N" %in% colnames(df)) + opts$nForPip <- stats::median(as.numeric(df$N), na.rm = TRUE) + # 2. Variant-content filters (MAF / INFO / N). Pure column-numeric # filters; the indel / strand-ambiguous variant-allele filtering happens # inside .matchAgainstSketch via harmonizeAlleles against the LD panel. @@ -2935,6 +3026,20 @@ krigingOutlierQc <- function(zScore, R, n, variantIds = NULL, #' @param alleleFlipKriging Logical (length 1). Opt-in kriging #' LD-consistency prefilter run before SLALOM/DENTIST. Default #' \code{FALSE}. +#' @param effectiveN Logical (length 1). When \code{TRUE} (default) and the +#' input carries case/control counts --- per-variant \code{N_CASE} / +#' \code{N_CONTROL} mcols, else the study-level \code{nCase} / +#' \code{nControl} scalars --- the working per-variant \code{N} is set to +#' the effective sample size \code{effectiveN(nCase, nControl)} BEFORE the +#' N-cutoff filter, so the filter, kriging, and the downstream fit all use +#' \code{N_eff}. When both counts and an \code{N} column are present the +#' counts win: \code{N} is overridden and the override is logged. Inputs +#' with no counts (quantitative traits) are unchanged. The escape hatch +#' \code{effectiveN = FALSE} restores the raw \code{N} column (or, when +#' there is no \code{N}, the raw total \code{nCase + nControl}) with no +#' override. \code{qcInfo$options$effectiveN} records the setting and each +#' entry's \code{nSource} is one of \code{"effective"}, \code{"column"}, +#' \code{"total"}, or \code{NA}. #' @param impute Logical (length 1). Run RAISS imputation against the #' \code{ldSketch}. Default \code{FALSE}. (Note: RAISS against the #' sketch is not yet fully wired for the new path; the option is @@ -2980,6 +3085,7 @@ summaryStatsQc <- function(sumstats, zMismatchQc = c("none", "slalom", "dentist"), alleleFlipKriging = FALSE, + effectiveN = TRUE, impute = FALSE, imputeOpts = list(rcond = 0.01, r2Threshold = 0.6, @@ -3024,6 +3130,9 @@ summaryStatsQc <- function(sumstats, pipCutoffToSkip = pipCutoffToSkip, zMismatchQc = zMismatchQc, alleleFlipKriging = alleleFlipKriging, + effectiveN = effectiveN, + nCase = NULL, + nControl = NULL, impute = impute, imputeOpts = imputeOpts, matchMinProp = matchMinProp, @@ -3045,6 +3154,10 @@ summaryStatsQc <- function(sumstats, opts$nForPip <- if ("N" %in% colnames(S4Vectors::mcols(sumstats$entry[[i]]))) stats::median(S4Vectors::mcols(sumstats$entry[[i]])$N, na.rm = TRUE) else NULL + # Per-entry study-level case/control counts (used only when the entry has + # 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 # 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. @@ -3075,6 +3188,7 @@ summaryStatsQc <- function(sumstats, nCutoff = nCutoff, zMismatchQc = zMismatchQc, alleleFlipKriging = alleleFlipKriging, + effectiveN = effectiveN, impute = impute, coerceNumeric = coerceNumeric, normalizeChr = normalizeChr, diff --git a/man/effectiveN.Rd b/man/effectiveN.Rd new file mode 100644 index 00000000..b76ba2d6 --- /dev/null +++ b/man/effectiveN.Rd @@ -0,0 +1,32 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/sumstatsQc.R +\name{effectiveN} +\alias{effectiveN} +\title{Effective sample size for a case/control study} +\usage{ +effectiveN(nCase, nControl) +} +\arguments{ +\item{nCase}{Numeric vector of case counts.} + +\item{nControl}{Numeric vector of control counts.} +} +\value{ +Numeric vector of effective sample sizes (\code{NA_real_} where a + count is missing or non-positive). +} +\description{ +Computes the effective sample size +\code{N_eff = 4 / (1/nCase + 1/nControl) = 4 * nCase * nControl / +(nCase + nControl)} for case/control GWAS. Balanced studies +(\code{nCase == nControl}) recover the total \code{nCase + nControl}; +imbalanced studies give a smaller value, which is the statistically +correct sample size for the RSS likelihood, residual-variance estimation, +kriging, and the N-cutoff filter. Vectorized over \code{nCase} / +\code{nControl}; entries where either count is \code{NA} or \code{<= 0} +return \code{NA_real_}. +} +\references{ +Privé et al., "Identifying and correcting for misspecifications + in GWAS summary statistics and polygenic scores", HGG Advances 2022. +} diff --git a/man/summaryStatsQc.Rd b/man/summaryStatsQc.Rd index 07af7c4d..cfc7d654 100644 --- a/man/summaryStatsQc.Rd +++ b/man/summaryStatsQc.Rd @@ -16,6 +16,7 @@ summaryStatsQc( pipCutoffToSkip = 0, zMismatchQc = c("none", "slalom", "dentist"), alleleFlipKriging = FALSE, + effectiveN = TRUE, impute = FALSE, imputeOpts = list(rcond = 0.01, r2Threshold = 0.6, minimumLd = 5, lamb = 0.01), matchMinProp = 0, @@ -70,6 +71,21 @@ Default 0 (no screen).} LD-consistency prefilter run before SLALOM/DENTIST. Default \code{FALSE}.} +\item{effectiveN}{Logical (length 1). When \code{TRUE} (default) and the +input carries case/control counts --- per-variant \code{N_CASE} / +\code{N_CONTROL} mcols, else the study-level \code{nCase} / +\code{nControl} scalars --- the working per-variant \code{N} is set to +the effective sample size \code{effectiveN(nCase, nControl)} BEFORE the +N-cutoff filter, so the filter, kriging, and the downstream fit all use +\code{N_eff}. When both counts and an \code{N} column are present the +counts win: \code{N} is overridden and the override is logged. Inputs +with no counts (quantitative traits) are unchanged. The escape hatch +\code{effectiveN = FALSE} restores the raw \code{N} column (or, when +there is no \code{N}, the raw total \code{nCase + nControl}) with no +override. \code{qcInfo$options$effectiveN} records the setting and each +entry's \code{nSource} is one of \code{"effective"}, \code{"column"}, +\code{"total"}, or \code{NA}.} + \item{impute}{Logical (length 1). Run RAISS imputation against the \code{ldSketch}. Default \code{FALSE}. (Note: RAISS against the sketch is not yet fully wired for the new path; the option is diff --git a/tests/testthat/test_sumstatsQc.R b/tests/testthat/test_sumstatsQc.R index 232b9754..f03bd738 100644 --- a/tests/testthat/test_sumstatsQc.R +++ b/tests/testthat/test_sumstatsQc.R @@ -2785,6 +2785,138 @@ test_that("summaryStatsQc: round-trips QtlSumStats inputs", { expect_equal(length(getQcInfo(res)$entryAudit), 1L) }) +# =========================================================================== +# effectiveN helper + case/control N canonicalization in summaryStatsQc +# =========================================================================== + +# Build a GwasSumStats entry carrying per-variant N_CASE / N_CONTROL mcols +# (optionally an N column too). Panel-matched (A1=A / A2=G at BP 100..) so the +# variants survive harmonization against .ssQ_makeHandle(). +.ssQ_makeCCEntry <- function(nCase, nControl, includeN = FALSE, + snp_ids = paste0("rs", 1:4), + positions = c(100L, 200L, 300L, 400L)) { + m <- length(snp_ids) + gr <- GenomicRanges::GRanges( + seqnames = rep("chr1", m), + ranges = IRanges::IRanges(start = positions, width = 1L)) + mc <- S4Vectors::DataFrame( + SNP = snp_ids, + A1 = rep("A", m), + A2 = rep("G", m), + Z = seq(1.0, by = 0.5, length.out = m), + N_CASE = as.numeric(rep_len(nCase, m)), + N_CONTROL = as.numeric(rep_len(nControl, m))) + if (includeN) mc$N <- rep(1000L, m) + S4Vectors::mcols(gr) <- mc + gr +} + +# Entry N ordered by genomic position (harmonization can reorder rows). +.ssQ_entryNByPos <- function(entry) { + o <- order(GenomicRanges::start(entry)) + as.numeric(S4Vectors::mcols(entry)$N[o]) +} + +test_that("effectiveN: balanced equals total, imbalanced is smaller, guards NA/<=0", { + # Balanced: 4/(1/500 + 1/500) == 1000 == total. + expect_equal(effectiveN(500, 500), 1000) + # Imbalanced: 4*100*900/1000 = 360 < 1000. + expect_equal(effectiveN(100, 900), 360) + expect_lt(effectiveN(100, 900), 100 + 900) + # NA / non-positive counts -> NA_real_. + expect_true(is.na(effectiveN(NA, 500))) + expect_true(is.na(effectiveN(500, NA))) + expect_true(is.na(effectiveN(0, 500))) + expect_true(is.na(effectiveN(-5, 500))) + # Vectorized, element-wise. + expect_equal(effectiveN(c(500, 100, 0), c(500, 900, 500)), + c(1000, 360, NA_real_)) +}) + +test_that("summaryStatsQc(effectiveN=TRUE): per-variant counts, no N -> N == N_eff", { + gr <- .ssQ_makeCCEntry(nCase = c(100, 200, 150, 250), + nControl = c(900, 800, 850, 750)) + ss <- GwasSumStats(study = "g1", entry = list(gr), genome = "hg19", + ldSketch = .ssQ_makeHandle()) + res <- summaryStatsQc(ss) + # 4*case*control/(case+control) per variant. + expect_equal(.ssQ_entryNByPos(res$entry[[1L]]), c(360, 640, 510, 750)) + expect_identical(getQcInfo(res)$entryAudit[[1L]]$nSource, "effective") +}) + +test_that("summaryStatsQc(effectiveN=TRUE): counts + N -> counts win, override logged", { + gr <- .ssQ_makeCCEntry(nCase = c(100, 200, 150, 250), + nControl = c(900, 800, 850, 750), + includeN = TRUE) + ss <- GwasSumStats(study = "g1", entry = list(gr), genome = "hg19", + ldSketch = .ssQ_makeHandle()) + expect_message(summaryStatsQc(ss), "overridden by effective N") + res <- summaryStatsQc(ss) + # N (was 1000) is replaced by the per-variant N_eff. + expect_equal(.ssQ_entryNByPos(res$entry[[1L]]), c(360, 640, 510, 750)) + expect_identical(getQcInfo(res)$entryAudit[[1L]]$nSource, "effective") +}) + +test_that("summaryStatsQc(effectiveN=TRUE): N only, no counts -> used as-is", { + ss <- .ssQ_makeGwasSumStats() # entry carries N = 1000, no counts + res <- summaryStatsQc(ss) + expect_true(all(.ssQ_entryNByPos(res$entry[[1L]]) == 1000)) + expect_identical(getQcInfo(res)$entryAudit[[1L]]$nSource, "column") +}) + +test_that("summaryStatsQc(effectiveN=FALSE): counts + N -> raw N, no override", { + gr <- .ssQ_makeCCEntry(nCase = c(100, 200, 150, 250), + nControl = c(900, 800, 850, 750), + includeN = TRUE) + ss <- GwasSumStats(study = "g1", entry = list(gr), genome = "hg19", + ldSketch = .ssQ_makeHandle()) + res <- summaryStatsQc(ss, effectiveN = FALSE) + expect_true(all(.ssQ_entryNByPos(res$entry[[1L]]) == 1000)) + expect_identical(getQcInfo(res)$entryAudit[[1L]]$nSource, "column") +}) + +test_that("summaryStatsQc(effectiveN=FALSE): counts only -> raw total, nSource='total'", { + gr <- .ssQ_makeCCEntry(nCase = c(100, 200, 150, 250), + nControl = c(900, 800, 850, 750)) + ss <- GwasSumStats(study = "g1", entry = list(gr), genome = "hg19", + ldSketch = .ssQ_makeHandle()) + res <- summaryStatsQc(ss, effectiveN = FALSE) + # Raw total n_case + n_control per variant. + expect_equal(.ssQ_entryNByPos(res$entry[[1L]]), c(1000, 1000, 1000, 1000)) + expect_identical(getQcInfo(res)$entryAudit[[1L]]$nSource, "total") +}) + +test_that("summaryStatsQc(effectiveN=TRUE): study-level scalars applied to all variants", { + # Entry has an N column but NO per-variant N_CASE/N_CONTROL; the scalars win. + gr <- .ssQ_makeEntryGr() + ss <- GwasSumStats(study = "g1", entry = list(gr), genome = "hg19", + ldSketch = .ssQ_makeHandle(), + nCase = 100, nControl = 900) + expect_message(summaryStatsQc(ss), "from study") + res <- summaryStatsQc(ss) + # 4*100*900/1000 = 360 for every variant. + expect_true(all(.ssQ_entryNByPos(res$entry[[1L]]) == 360)) + expect_identical(getQcInfo(res)$entryAudit[[1L]]$nSource, "effective") +}) + +test_that("summaryStatsQc: effectiveN recorded in qcInfo options", { + ss <- .ssQ_makeGwasSumStats() + expect_true(getQcInfo(summaryStatsQc(ss))$options$effectiveN) + expect_false(getQcInfo(summaryStatsQc(ss, effectiveN = FALSE))$options$effectiveN) +}) + +test_that("summaryStatsQc: quantitative QtlSumStats is a no-op for effective N", { + # No counts anywhere; entry carries an N column -> used as-is (nSource + # "column"), N untouched. + gr <- .ssQ_makeEntryGr() + ss <- QtlSumStats(study = "s1", context = "c1", trait = "t1", + entry = list(gr), genome = "hg19", + ldSketch = .ssQ_makeHandle()) + res <- summaryStatsQc(ss) + expect_true(all(.ssQ_entryNByPos(res$entry[[1L]]) == 1000)) + expect_identical(getQcInfo(res)$entryAudit[[1L]]$nSource, "column") +}) + # =========================================================================== # summaryStatsQc with LD-mismatch QC enabled (mocked extractor) # ===========================================================================