diff --git a/R/colocboostPipeline.R b/R/colocboostPipeline.R index eb9c599e..af01227f 100644 --- a/R/colocboostPipeline.R +++ b/R/colocboostPipeline.R @@ -11,6 +11,9 @@ #' shared LD reference (\code{ldSketch}). Must already have #' been passed through \code{\link{summaryStatsQc}} (the #' pipeline rejects inputs whose \code{getQcInfo()} is empty). +#' Every analysis variant is available on this input, including +#' \code{xqtlColoc}, which colocalizes the QTL studies against +#' each other with no GWAS involved. #' \item \code{MultiStudyQtlDataset} -- a mixture of one or more #' individual-level \code{QtlDataset} studies and an optional #' \code{QtlSumStats} collection. @@ -50,7 +53,12 @@ #' @section Analysis variants: #' \itemize{ #' \item \code{xqtlColoc} (default \code{TRUE}): run a colocboost -#' model over the QTL contexts only (individual-level inputs). +#' model over the QTL outcomes only, excluding +#' \code{gwasSumStats}. Works for either data form -- the +#' individual-level contexts, the summary-level QTL studies, or +#' both together when a \code{MultiStudyQtlDataset} carries a +#' mixture. This is the only variant that honors +#' \code{focalTrait}. #' \item \code{jointGwas} (default \code{FALSE}): run a non-focal #' colocboost model that combines all QTL contexts/studies #' with the supplied \code{gwasSumStats} studies. @@ -80,7 +88,9 @@ #' analysis to; \code{NULL} (default) uses all samples. #' @param focalTrait Optional trait name; when supplied and present in the #' assembled outcome list, the colocboost xQTL-only run uses it as the focal -#' outcome. +#' outcome. Only \code{xqtlColoc} reads it: \code{jointGwas} is +#' non-focal by construction and \code{separateGwas} always makes the +#' GWAS study focal. #' @param xqtlColoc,jointGwas,separateGwas Logical flags selecting which #' colocboost variants to run. #' @param pipCutoffToSkip Individual-level pre-filter (ports the legacy @@ -780,9 +790,25 @@ setGeneric("colocboostPipeline", function(qtlData, gwasSumStats = NULL, ...) { ) } +# A requested analysis with nothing to run on would otherwise be a silent +# no-op: the caller gets an empty ColocBoostResult whose getComputingTime() +# entries are all NULL, with no indication of why. Say so. +# @noRd +.cbWarnNoData <- function(flag, needed) { + msg <- glue( + "colocboostPipeline: {flag} = TRUE was requested, but there is no ", + "{needed} to run it on. Skipping it -- the returned ", + "ColocBoostResult will hold no confidence sets from this analysis ", + "and its getComputingTime() entry will be NULL." + ) + warn(msg) +} + # Shared dispatch: accepts a fully-prepared individual bundle (possibly # NULL) plus a sumstat bundle (possibly empty) and runs the three -# colocboost variants the user requested. +# colocboost variants the user requested. `qtlSumstatBundle` is the QTL-side +# subset of `sumstatBundle` (see .cbDriver); only the xQTL-only run uses it, +# so that a GWAS study never becomes an xQTL-only outcome. .cbRunVariants <- function( individualBundle, sumstatBundle, @@ -791,11 +817,14 @@ setGeneric("colocboostPipeline", function(qtlData, gwasSumStats = NULL, ...) { separateGwas, focalTrait, dotArgs, - qtlLdSketch = NULL + qtlLdSketch = NULL, + qtlSumstatBundle = NULL ) { results <- .cbEmptyResult() hasInd <- !is.null(individualBundle) hasSs <- length(sumstatBundle$sumstat) > 0L + qtlSumstatBundle <- qtlSumstatBundle %||% .cbMergeSumstatBundles(list()) + hasQtlSs <- length(qtlSumstatBundle$sumstat) > 0L if (!hasInd && !hasSs) { msg <- glue( "colocboostPipeline: no QTL inputs remain after selection. ", @@ -804,25 +833,48 @@ setGeneric("colocboostPipeline", function(qtlData, gwasSumStats = NULL, ...) { inform(msg) return(.cbEmptyResultObject()) } - if (isTRUE(xqtlColoc) && hasInd) { - run <- .cbRunXqtlOnly(individualBundle, focalTrait, dotArgs) - results$xqtl_coloc <- run$result - results$computing_time$Analysis$xqtl_coloc <- run$time - } - if (isTRUE(jointGwas) && hasSs) { - run <- .cbRunJointGwas(individualBundle, sumstatBundle, hasInd, dotArgs) - results$joint_gwas <- run$result - results$computing_time$Analysis$joint_gwas <- run$time - } - if (isTRUE(separateGwas) && hasSs) { - run <- .cbRunSeparateGwas( - individualBundle, - sumstatBundle, - hasInd, - dotArgs - ) - results$separate_gwas <- run$result - results$computing_time$Analysis$separate_gwas <- run$time + if (isTRUE(xqtlColoc)) { + if (hasInd || hasQtlSs) { + run <- .cbRunXqtlOnly( + individualBundle, + qtlSumstatBundle, + hasInd, + focalTrait, + dotArgs + ) + results$xqtl_coloc <- run$result + results$computing_time$Analysis$xqtl_coloc <- run$time + } else { + .cbWarnNoData("xqtlColoc", "QTL data") + } + } + if (isTRUE(jointGwas)) { + if (hasSs) { + run <- .cbRunJointGwas( + individualBundle, + sumstatBundle, + hasInd, + dotArgs + ) + results$joint_gwas <- run$result + results$computing_time$Analysis$joint_gwas <- run$time + } else { + .cbWarnNoData("jointGwas", "summary-statistic data") + } + } + if (isTRUE(separateGwas)) { + if (hasSs) { + run <- .cbRunSeparateGwas( + individualBundle, + sumstatBundle, + hasInd, + dotArgs + ) + results$separate_gwas <- run$result + results$computing_time$Analysis$separate_gwas <- run$time + } else { + .cbWarnNoData("separateGwas", "summary-statistic data") + } } .cbToResultObject( results, @@ -886,28 +938,59 @@ setGeneric("colocboostPipeline", function(qtlData, gwasSumStats = NULL, ...) { } # xQTL-only ColocBoost run -> list(result, time). +# +# Either side may be absent: `individualBundle` is NULL for a summary-level +# QTL input, and `sumstatBundle` is empty for an individual-level one. It holds +# the QTL-side sumstats ONLY -- a GWAS study must never be pulled into the +# xQTL-only analysis, which is why this does not take the merged bundle the +# joint / separate runs use. # @noRd -.cbRunXqtlOnly <- function(individualBundle, focalTrait, dotArgs) { - traits <- individualBundle$outcomeNames +.cbRunXqtlOnly <- function( + individualBundle, + sumstatBundle, + hasInd, + focalTrait, + dotArgs +) { + traits <- c( + if (hasInd) individualBundle$outcomeNames else character(), + names(sumstatBundle$sumstat) + ) focalIdx <- if (!is.null(focalTrait) && is_in(focalTrait, traits)) { which(traits == focalTrait) } else { NULL } - nCtx <- length(individualBundle$Y) - msg <- glue( - "====== Performing xQTL-only ColocBoost on {nCtx} contexts. =====" - ) + nCtx <- if (hasInd) length(individualBundle$Y) else 0L + nSs <- length(sumstatBundle$sumstat) + msg <- if (nSs > 0L) { + glue( + "====== Performing xQTL-only ColocBoost on {nCtx} contexts ", + "and {nSs} summary-statistic studies. =====" + ) + } else { + glue( + "====== Performing xQTL-only ColocBoost on {nCtx} contexts. =====" + ) + } inform(msg) + ldArgs <- if (nSs > 0L) .cbBuildLdArgs(sumstatBundle$LD) else list() args <- c( list( - X = individualBundle$X, - Y = individualBundle$Y, - dict_YX = individualBundle$dict_YX, + X = if (hasInd) individualBundle$X else NULL, + Y = if (hasInd) individualBundle$Y else NULL, + dict_YX = if (hasInd) individualBundle$dict_YX else NULL, + sumstat = if (nSs > 0L) sumstatBundle$sumstat else NULL, + dict_sumstatLD = if (nSs > 0L) { + sumstatBundle$dict_sumstatLD + } else { + NULL + }, outcome_names = traits, focal_outcome_idx = focalIdx, output_level = 2 ), + ldArgs, dotArgs ) run <- .cbRun("xQTL-only ColocBoost", args) @@ -1190,6 +1273,12 @@ setGeneric("colocboostPipeline", function(qtlData, gwasSumStats = NULL, ...) { combinedPairs <- harmonized$pairs } sumstatBundle <- .cbMergeSumstatBundles(combinedPairs) + # The xQTL-only run gets its own bundle over just the QTL-side pairs, so + # a GWAS study is never treated as an xQTL outcome. Rebuilding it through + # .cbMergeSumstatBundles (rather than subsetting the merged one) keeps the + # deduplicated LD list and dict_sumstatLD consistent for the subset. + qtlKeys <- intersect(names(combinedPairs), names(qtlPairs)) + qtlSumstatBundle <- .cbMergeSumstatBundles(combinedPairs[qtlKeys]) .cbRunVariants( individualBundle, sumstatBundle, @@ -1198,7 +1287,8 @@ setGeneric("colocboostPipeline", function(qtlData, gwasSumStats = NULL, ...) { separateGwas, focalTrait, dotArgs, - qtlLdSketch = qtlLdSketch + qtlLdSketch = qtlLdSketch, + qtlSumstatBundle = qtlSumstatBundle ) } diff --git a/man/colocboostPipeline.Rd b/man/colocboostPipeline.Rd index 470e1adc..f1ca6ec3 100644 --- a/man/colocboostPipeline.Rd +++ b/man/colocboostPipeline.Rd @@ -106,7 +106,9 @@ accessors).} \item{focalTrait}{Optional trait name; when supplied and present in the assembled outcome list, the colocboost xQTL-only run uses it as the focal -outcome.} +outcome. Only \code{xqtlColoc} reads it: \code{jointGwas} is +non-focal by construction and \code{separateGwas} always makes the +GWAS study focal.} \item{xqtlColoc, jointGwas, separateGwas}{Logical flags selecting which colocboost variants to run.} @@ -172,6 +174,9 @@ Protocol-level multi-trait colocalization analysis using shared LD reference (\code{ldSketch}). Must already have been passed through \code{\link{summaryStatsQc}} (the pipeline rejects inputs whose \code{getQcInfo()} is empty). + Every analysis variant is available on this input, including + \code{xqtlColoc}, which colocalizes the QTL studies against + each other with no GWAS involved. \item \code{MultiStudyQtlDataset} -- a mixture of one or more individual-level \code{QtlDataset} studies and an optional \code{QtlSumStats} collection. @@ -214,7 +219,12 @@ for either side; colocboost has its own variable-selection algorithm. \itemize{ \item \code{xqtlColoc} (default \code{TRUE}): run a colocboost - model over the QTL contexts only (individual-level inputs). + model over the QTL outcomes only, excluding + \code{gwasSumStats}. Works for either data form -- the + individual-level contexts, the summary-level QTL studies, or + both together when a \code{MultiStudyQtlDataset} carries a + mixture. This is the only variant that honors + \code{focalTrait}. \item \code{jointGwas} (default \code{FALSE}): run a non-focal colocboost model that combines all QTL contexts/studies with the supplied \code{gwasSumStats} studies. diff --git a/tests/testthat/test_colocboostPipeline.R b/tests/testthat/test_colocboostPipeline.R index d54815fe..f6b68352 100644 --- a/tests/testthat/test_colocboostPipeline.R +++ b/tests/testthat/test_colocboostPipeline.R @@ -1188,7 +1188,12 @@ test_that(".cbResidualizedX reports why genotypes were unavailable", { # The underlying message is carried through so the skip is diagnosable. expect_message( res <- pecotmr:::.cbResidualizedX( - NULL, "c1", NULL, NULL, NULL, NULL + NULL, + "c1", + NULL, + NULL, + NULL, + NULL ), "residualized genotypes unavailable: kaboom" ) @@ -1233,17 +1238,162 @@ test_that(".cbRunXqtlOnly passes a focal outcome through as an index", { X = list(), dict_YX = NULL ) - run <- suppressMessages(pecotmr:::.cbRunXqtlOnly(bundle, "tB", list())) + empty <- pecotmr:::.cbMergeSumstatBundles(list()) + run <- suppressMessages( + pecotmr:::.cbRunXqtlOnly(bundle, empty, TRUE, "tB", list()) + ) # colocboost wants a position, not a name. expect_equal(run$result$focal, 2L) absent <- suppressMessages( - pecotmr:::.cbRunXqtlOnly(bundle, "nope", list()) + pecotmr:::.cbRunXqtlOnly(bundle, empty, TRUE, "nope", list()) ) expect_null(absent$result$focal) - none <- suppressMessages(pecotmr:::.cbRunXqtlOnly(bundle, NULL, list())) + none <- suppressMessages( + pecotmr:::.cbRunXqtlOnly(bundle, empty, TRUE, NULL, list()) + ) expect_null(none$result$focal) }) +# A summary-level QTL side used to make xqtlColoc a silent no-op: the run was +# gated on having an individual bundle, so the default flags returned an empty +# ColocBoostResult with all-NULL timings and no message. +test_that(".cbRunXqtlOnly runs on a summary-level QTL side alone", { + captured <- NULL + local_mocked_bindings( + .cbRun = function(label, args) { + captured <<- args + list(result = list(ran = TRUE), time = 0) + }, + .package = "pecotmr" + ) + ssBundle <- pecotmr:::.cbMergeSumstatBundles(list( + qtlA = list(sumstat = list(z = 1), LD = diag(2)), + qtlB = list(sumstat = list(z = 2), LD = diag(2)) + )) + run <- suppressMessages( + pecotmr:::.cbRunXqtlOnly(NULL, ssBundle, FALSE, "qtlB", list()) + ) + expect_equal(run$result$ran, TRUE) + expect_null(captured$X) + expect_null(captured$Y) + expect_null(captured$dict_YX) + expect_equal(names(captured$sumstat), c("qtlA", "qtlB")) + expect_equal(captured$outcome_names, c("qtlA", "qtlB")) + # focalTrait is honored on the summary-level side too. + expect_equal(captured$focal_outcome_idx, 2L) + # Identical LD matrices dedupe to one, so the dict points both at it. + expect_equal(unname(captured$dict_sumstatLD[, "LD"]), c(1L, 1L)) +}) + +test_that(".cbRunXqtlOnly combines individual and summary-level QTL sides", { + captured <- NULL + local_mocked_bindings( + .cbRun = function(label, args) { + captured <<- args + list(result = NULL, time = 0) + }, + .package = "pecotmr" + ) + bundle <- list( + outcomeNames = c("tA", "tB"), + Y = list(1, 2), + X = list(), + dict_YX = NULL + ) + ssBundle <- pecotmr:::.cbMergeSumstatBundles(list( + qtlC = list(sumstat = list(z = 3), LD = diag(2)) + )) + suppressMessages( + pecotmr:::.cbRunXqtlOnly(bundle, ssBundle, TRUE, "qtlC", list()) + ) + expect_equal(captured$outcome_names, c("tA", "tB", "qtlC")) + expect_equal(captured$focal_outcome_idx, 3L) +}) + +test_that(".cbRunVariants: xqtlColoc runs on a QTL-only sumstat bundle", { + called <- character(0) + local_mocked_bindings( + .cbRunXqtlOnly = function(...) { + called <<- c(called, "xqtl") + list(result = list(stub = TRUE), time = 1) + }, + .cbOutcomeInfo = function(...) pecotmr:::.cbEmptyOutcomeInfo(), + .package = "pecotmr" + ) + merged <- pecotmr:::.cbMergeSumstatBundles(list( + qtlA = list(sumstat = list(z = 1), LD = diag(2)), + gwasG = list(sumstat = list(z = 2), LD = diag(2)) + )) + qtlOnly <- pecotmr:::.cbMergeSumstatBundles(list( + qtlA = list(sumstat = list(z = 1), LD = diag(2)) + )) + out <- suppressMessages(pecotmr:::.cbRunVariants( + NULL, + merged, + xqtlColoc = TRUE, + jointGwas = FALSE, + separateGwas = FALSE, + focalTrait = NULL, + dotArgs = list(), + qtlSumstatBundle = qtlOnly + )) + expect_equal(called, "xqtl") + expect_false(is.null(getComputingTime(out)$Analysis$xqtl_coloc)) +}) + +test_that(".cbRunVariants warns instead of silently skipping an analysis", { + local_mocked_bindings( + .cbOutcomeInfo = function(...) pecotmr:::.cbEmptyOutcomeInfo(), + .package = "pecotmr" + ) + # Individual-level QTL side, no sumstats anywhere: the two GWAS variants + # cannot run, and the caller is told so rather than getting a bare empty. + ind <- list( + outcomeNames = "tA", + Y = list(1), + X = list(), + dict_YX = NULL + ) + warnings <- capture_warnings( + suppressMessages(pecotmr:::.cbRunVariants( + ind, + pecotmr:::.cbMergeSumstatBundles(list()), + xqtlColoc = FALSE, + jointGwas = TRUE, + separateGwas = TRUE, + focalTrait = NULL, + dotArgs = list() + )) + ) + expect_length(warnings, 2L) + expect_match(warnings[[1L]], "jointGwas = TRUE was requested") + expect_match(warnings[[2L]], "separateGwas = TRUE was requested") +}) + +test_that(".cbRunVariants warns when xqtlColoc has only GWAS sumstats", { + local_mocked_bindings( + .cbOutcomeInfo = function(...) pecotmr:::.cbEmptyOutcomeInfo(), + .package = "pecotmr" + ) + gwasOnly <- pecotmr:::.cbMergeSumstatBundles(list( + gwasG = list(sumstat = list(z = 1), LD = diag(2)) + )) + out <- expect_warning( + suppressMessages(pecotmr:::.cbRunVariants( + NULL, + gwasOnly, + xqtlColoc = TRUE, + jointGwas = FALSE, + separateGwas = FALSE, + focalTrait = NULL, + dotArgs = list(), + qtlSumstatBundle = pecotmr:::.cbMergeSumstatBundles(list()) + )), + "xqtlColoc = TRUE was requested" + ) + expect_null(getComputingTime(out)$Analysis$xqtl_coloc) +}) + test_that(".cbAppendGwasPairs disambiguates a colliding study key", { local_mocked_bindings( .cbRequireSumStatsQc = function(...) invisible(NULL),