From b7d2234db4d09e87794c2c76a3a9a97e02a6f597 Mon Sep 17 00:00:00 2001 From: Daniel Nachun Date: Fri, 4 Sep 2026 17:23:27 -0700 Subject: [PATCH 1/2] fix other susie args --- R/fineMappingWrappers.R | 24 ++++++++++++++-------- R/twasWeights.R | 4 ++-- man/fineMappingPipeline.Rd | 15 +++++++------- man/mvsusieRssWeights.Rd | 11 ++++++---- man/mvsusieWeights.Rd | 11 ++++++---- tests/testthat/test_fineMappingWrappers.R | 25 ++++++++++++++++++++++- tests/testthat/test_twasWeights.R | 6 ++++-- 7 files changed, 68 insertions(+), 28 deletions(-) diff --git a/R/fineMappingWrappers.R b/R/fineMappingWrappers.R index bed515f6..fde45fa3 100644 --- a/R/fineMappingWrappers.R +++ b/R/fineMappingWrappers.R @@ -2708,8 +2708,11 @@ susieAshRssWeights <- function( #' @param Y Phenotype matrix. Required when `mvsusieFit` is NULL. #' @param priorVariance Optional mvSuSiE prior variance list. #' @param residualVariance Optional residual variance matrix. -#' @param L Maximum number of components. -#' @param LGreedy Initial greedy number of components. +#' @param L Maximum number of components. Default \code{10}, matching +#' \code{mvsusieR::mvsusie}. +#' @param LGreedy Integer or \code{NULL}. Number of greedily-added components. +#' \code{NULL} (default) disables the greedy loop and fits \code{L} +#' directly. #' @param verbose If TRUE, prints mvSuSiE fitting progress. #' @param ... Additional arguments passed to `fitMvsusie()` when fitting. #' @return Matrix of variant weights. @@ -2729,8 +2732,8 @@ mvsusieWeights <- function( Y = NULL, priorVariance = NULL, residualVariance = NULL, - L = 30, - LGreedy = 5, + L = 10, + LGreedy = NULL, verbose = FALSE, ... ) { @@ -3021,8 +3024,11 @@ fsusieWeights <- function( #' NULL, \code{mvsusieR::create_mixture_prior()} is used with \code{R = #' ncol(stat$z)}. #' @param residualVariance Optional residual covariance matrix. -#' @param L Maximum number of single effects (default 30). -#' @param LGreedy Initial greedy effect count (default 5). +#' @param L Maximum number of single effects. Default \code{10}, matching +#' \code{mvsusieR::mvsusie}. +#' @param LGreedy Integer or \code{NULL}. Number of greedily-added effects. +#' \code{NULL} (default) disables the greedy loop and fits \code{L} +#' directly. #' @param retainFit If TRUE, attaches the fitted object as an attribute. #' @param ... Additional arguments forwarded to \code{mvsusieR::mvsusie_rss}. #' @@ -3051,8 +3057,8 @@ mvsusieRssWeights <- function( mvsusieRssFit = NULL, priorVariance = NULL, residualVariance = NULL, - L = 30, - LGreedy = 5, + L = 10, + LGreedy = NULL, retainFit = FALSE, ... ) { @@ -3115,6 +3121,8 @@ mvsusieRssWeights <- function( Z = Z, R = LD, N = nScalar, + L = L, + L_greedy = LGreedy, prior_variance = priorVariance, residual_variance = residualVariance, ... diff --git a/R/twasWeights.R b/R/twasWeights.R index f19958b1..1e3699bd 100644 --- a/R/twasWeights.R +++ b/R/twasWeights.R @@ -576,7 +576,7 @@ setMethod("show", "TwasWeights", function(object) { susie = list( fn = "susie_weights", impl = "susieWeights", - args = list(refine = FALSE, L = 20, L_greedy = 5) + args = list(refine = FALSE, L = 10) ), susieAsh = list( fn = "susie_ash_weights", @@ -659,7 +659,7 @@ setMethod("show", "TwasWeights", function(object) { mvsusie = list( fn = "mvsusie_weights", impl = "mvsusieWeights", - args = list(L = 30, L_greedy = 5) + args = list(L = 10) ), mrmash = list( fn = "mrmash_weights", diff --git a/man/fineMappingPipeline.Rd b/man/fineMappingPipeline.Rd index b9b6b454..bbc21be5 100644 --- a/man/fineMappingPipeline.Rd +++ b/man/fineMappingPipeline.Rd @@ -28,8 +28,8 @@ fineMappingPipeline(data, ...) jointRegions = FALSE, jointSpecification = NULL, addSusieInf = TRUE, - L = 20L, - Lgreedy = 5L, + L = 10L, + Lgreedy = NULL, coverage = 0.95, secondaryCoverage = c(0.7, 0.5), signalCutoff = 0.025, @@ -135,8 +135,8 @@ fineMappingPipeline(data, ...) data, methods, addSusieInf = TRUE, - L = 20L, - Lgreedy = 5L, + L = 10L, + Lgreedy = NULL, coverage = 0.95, secondaryCoverage = c(0.7, 0.5), signalCutoff = 0.025, @@ -251,10 +251,11 @@ alongside \code{susie} and/or \code{susieAsh}, controls whether the SuSiE-inf fit initialises the chained downstream method(s). Default \code{TRUE}.} -\item{L}{Integer. Maximum number of SuSiE single effects. Default \code{20}.} +\item{L}{Integer. Maximum number of SuSiE single effects. Default \code{10}.} -\item{Lgreedy}{Integer. Number of greedily-added effects in the SuSiE-inf -refinement. Default \code{5}.} +\item{Lgreedy}{Integer or \code{NULL}. Number of greedily-added effects in the +SuSiE-inf refinement (the greedy-L loop). \code{NULL} (default) disables the +greedy loop and fits \code{L} directly.} \item{coverage}{Primary credible-set coverage (numeric, length 1). Default \code{0.95}.} diff --git a/man/mvsusieRssWeights.Rd b/man/mvsusieRssWeights.Rd index 65439691..b3505080 100644 --- a/man/mvsusieRssWeights.Rd +++ b/man/mvsusieRssWeights.Rd @@ -10,8 +10,8 @@ mvsusieRssWeights( mvsusieRssFit = NULL, priorVariance = NULL, residualVariance = NULL, - L = 30, - LGreedy = 5, + L = 10, + LGreedy = NULL, retainFit = FALSE, ... ) @@ -30,9 +30,12 @@ ncol(stat$z)}.} \item{residualVariance}{Optional residual covariance matrix.} -\item{L}{Maximum number of single effects (default 30).} +\item{L}{Maximum number of single effects. Default \code{10}, matching +\code{mvsusieR::mvsusie}.} -\item{LGreedy}{Initial greedy effect count (default 5).} +\item{LGreedy}{Integer or \code{NULL}. Number of greedily-added effects. +\code{NULL} (default) disables the greedy loop and fits \code{L} +directly.} \item{retainFit}{If TRUE, attaches the fitted object as an attribute.} diff --git a/man/mvsusieWeights.Rd b/man/mvsusieWeights.Rd index 1c1cd4a8..b52b2c1c 100644 --- a/man/mvsusieWeights.Rd +++ b/man/mvsusieWeights.Rd @@ -10,8 +10,8 @@ mvsusieWeights( Y = NULL, priorVariance = NULL, residualVariance = NULL, - L = 30, - LGreedy = 5, + L = 10, + LGreedy = NULL, verbose = FALSE, ... ) @@ -27,9 +27,12 @@ mvsusieWeights( \item{residualVariance}{Optional residual variance matrix.} -\item{L}{Maximum number of components.} +\item{L}{Maximum number of components. Default \code{10}, matching +\code{mvsusieR::mvsusie}.} -\item{LGreedy}{Initial greedy number of components.} +\item{LGreedy}{Integer or \code{NULL}. Number of greedily-added components. +\code{NULL} (default) disables the greedy loop and fits \code{L} +directly.} \item{verbose}{If TRUE, prints mvSuSiE fitting progress.} diff --git a/tests/testthat/test_fineMappingWrappers.R b/tests/testthat/test_fineMappingWrappers.R index 33be971e..09cbe6aa 100644 --- a/tests/testthat/test_fineMappingWrappers.R +++ b/tests/testthat/test_fineMappingWrappers.R @@ -2720,11 +2720,34 @@ test_that("susieAshRssWeights returns weights of length p", { test_that("mvsusieRssWeights fits mvsusie_rss and returns p x K weights", { skip_if_not_installed("mvsusieR") m <- .rrwMulti(n = 80, p = 8, K = 2) - w <- mvsusieRssWeights(m$stat, m$LD, L = 5, LGreedy = 2) + # LGreedy stays at its NULL default: mvsusieR's greedy-L loop errors on + # small data ("Values and their weights should have equal length"), which + # is why the greedy loop is off unless a caller asks for it. + w <- mvsusieRssWeights(m$stat, m$LD, L = 5) expect_equal(dim(w), c(m$p, m$K)) expect_true(all(is.finite(w))) }) +test_that("mvsusieRssWeights forwards L / clamped LGreedy to mvsusie_rss", { + skip_if_not_installed("mvsusieR") + m <- .rrwMulti(n = 80, p = 8, K = 2) + fakeCoef <- matrix(rnorm((m$p + 1) * m$K), nrow = m$p + 1, ncol = m$K) + captured <- list() + local_mocked_bindings( + create_mixture_prior = function(...) list(), + mvsusie_rss = function(...) { + captured <<- list(...) + "mock_fit" + }, + coef.mvsusie = function(...) fakeCoef, + .package = "mvsusieR" + ) + w <- mvsusieRssWeights(m$stat, m$LD, L = 3, LGreedy = 7) + expect_equal(dim(w), c(m$p, m$K)) + expect_equal(captured$L, 3) + expect_equal(captured$L_greedy, 3) # clamped to min(LGreedy, L) +}) + test_that("mvsusieRssWeights errors on single-context stat$z", { skip_if_not_installed("mvsusieR") f <- .rrwStatLd() diff --git a/tests/testthat/test_twasWeights.R b/tests/testthat/test_twasWeights.R index f184e041..cd50e12a 100644 --- a/tests/testthat/test_twasWeights.R +++ b/tests/testthat/test_twasWeights.R @@ -141,8 +141,10 @@ test_that(".twas_method_lookup: unknown method produces error", { test_that(".twas_method_lookup: default args are set for susie and mrash", { result <- pecotmr:::.twasMethodLookup("fast_default") expect_equal(result$susie_weights$refine, FALSE) - expect_equal(result$susie_weights$L, 20) - expect_equal(result$susie_weights$L_greedy, 5) + # Matches susieR::susie's own defaults (L = min(10, p), greedy loop off), + # as fineMappingPipeline does. + expect_equal(result$susie_weights$L, 10) + expect_null(result$susie_weights$L_greedy) expect_equal(result$mrash_weights$initPriorSd, TRUE) expect_equal(result$mrash_weights$max.iter, 100) }) From 08ab4a2388016369fe5ceb67ef1ac6ca67d3fee8 Mon Sep 17 00:00:00 2001 From: Daniel Nachun Date: Mon, 7 Sep 2026 18:36:52 -0700 Subject: [PATCH 2/2] make colocPipeline more flexible --- R/ColocResult.R | 40 +- R/colocPipeline.R | 361 ++++++++++------ R/qtlEnrichmentPipeline.R | 455 +++++++++++++------- R/tupleSelectors.R | 45 ++ man/ColocResult.Rd | 14 +- man/as.data.frame-ColocResult-method.Rd | 4 +- man/colocPipeline.Rd | 103 +++-- man/qtlEnrichment.Rd | 10 +- man/qtlEnrichmentPipeline.Rd | 70 ++- src/qtl_enrichment.h | 26 +- tests/testthat/test_ColocResult.R | 23 +- tests/testthat/test_colocPipeline.R | 199 ++++++++- tests/testthat/test_qtlEnrichmentPipeline.R | 337 +++++++++++---- 13 files changed, 1219 insertions(+), 468 deletions(-) diff --git a/R/ColocResult.R b/R/ColocResult.R index 35c9a62c..60b3fd18 100644 --- a/R/ColocResult.R +++ b/R/ColocResult.R @@ -97,6 +97,8 @@ methods::setValidity("ColocResult", function(object) { "trait", "method", "gwasStudy", + "gwasContext", + "gwasTrait", "gwasMethod", "blockId", "qtlCs", @@ -169,8 +171,8 @@ setMethod("show", "ColocResult", function(object) { return(invisible(NULL)) } md <- mcols(object, use.names = FALSE) - cat(" QTL studies :", str_flatten(unique(md$study), ", "), "\n") - cat(" GWAS studies:", str_flatten(unique(md$gwasStudy), ", "), "\n") + cat(" studies :", str_flatten(unique(md$study), ", "), "\n") + cat(" paired with :", str_flatten(unique(md$gwasStudy), ", "), "\n") cat(" variants :", sum(lengths(object)), "across all pairs\n") cat( " max PP.H4 :", @@ -187,10 +189,12 @@ setMethod("show", "ColocResult", function(object) { #' table and the per-pair variant tables that go with it. Callers normally #' get one from \code{\link{colocPipeline}} rather than building it directly. #' @param pairs A data frame with one row per tested pair, carrying at least -#' the identity columns (\code{study}, \code{context}, \code{trait}, -#' \code{method}, \code{gwasStudy}, \code{gwasMethod}), \code{blockId}, -#' \code{qtlCs}, \code{gwasCs}, \code{nSnps} and \code{PP.H0.abf} through -#' \code{PP.H4.abf}. +#' the identity columns of both sides (\code{study}, \code{context}, +#' \code{trait}, \code{method} and \code{gwasStudy}, \code{gwasContext}, +#' \code{gwasTrait}, \code{gwasMethod}), \code{blockId}, \code{qtlCs}, +#' \code{gwasCs}, \code{nSnps} and \code{PP.H0.abf} through +#' \code{PP.H4.abf}. A GWAS side has no context or trait axis, so those two +#' columns are \code{NA} for it. #' @param variants A list, parallel to \code{pairs}' rows, of per-pair data #' frames with a \code{variant_id} column and a \code{SNP.PP.H4} column. #' @param ldSketch Optional genotype panel (see \code{\link{readGenotypes}}) @@ -200,7 +204,9 @@ setMethod("show", "ColocResult", function(object) { #' @examples #' pairs <- data.frame( #' study = "s1", context = "c1", trait = "g1", method = "susie", -#' gwasStudy = "G1", gwasMethod = "susie", blockId = "chr1_1_1000", +#' gwasStudy = "G1", gwasContext = NA_character_, +#' gwasTrait = NA_character_, gwasMethod = "susie", +#' blockId = "chr1_1_1000", #' qtlCs = 1L, gwasCs = 1L, nSnps = 2L, #' PP.H0.abf = 0.1, PP.H1.abf = 0.1, PP.H2.abf = 0.1, #' PP.H3.abf = 0.1, PP.H4.abf = 0.6 @@ -301,10 +307,22 @@ ColocResult <- function(pairs, variants, ldSketch = NULL) { # ---- views ------------------------------------------------------------------ # The columns that identify a gene-level unit: everything that is fixed within -# one QTL molecular trait tested against one GWAS study. +# one first-side trait tested against one second-side trait. Both sides carry +# the full (study, context, trait, method) tuple, since either may be a QTL +# collection whose rows differ only on context or trait -- grouping on the +# study alone would pool two distinct molecular phenotypes into one unit. # @noRd .crGeneCols <- function() { - c("study", "context", "trait", "method", "gwasStudy", "gwasMethod") + c( + "study", + "context", + "trait", + "method", + "gwasStudy", + "gwasContext", + "gwasTrait", + "gwasMethod" + ) } #' @rdname colocViews @@ -560,7 +578,9 @@ setMethod( #' @examples #' pairs <- data.frame( #' study = "s1", context = "c1", trait = "g1", method = "susie", -#' gwasStudy = "G1", gwasMethod = "susie", blockId = "chr1_1_1000", +#' gwasStudy = "G1", gwasContext = NA_character_, +#' gwasTrait = NA_character_, gwasMethod = "susie", +#' blockId = "chr1_1_1000", #' qtlCs = 1L, gwasCs = 1L, nSnps = 2L, #' PP.H0.abf = 0.1, PP.H1.abf = 0.1, PP.H2.abf = 0.1, #' PP.H3.abf = 0.1, PP.H4.abf = 0.6 diff --git a/R/colocPipeline.R b/R/colocPipeline.R index 7dc9e2fe..79c9855b 100644 --- a/R/colocPipeline.R +++ b/R/colocPipeline.R @@ -1,10 +1,23 @@ -#' @title Colocalization Pipeline (coloc.bf_bf over QTL + GWAS LBF matrices) -#' @description Per-region pipeline that pairs a QTL -#' \code{\link{QtlFineMappingResult}} with a GWAS fine-mapping result (either -#' supplied directly as a \code{\link{GwasFineMappingResult}} or computed -#' inline from a \code{\link{GwasSumStats}}) and runs -#' \code{coloc::coloc.bf_bf} per (QTL tuple, GWAS tuple) pair to produce -#' per-pair colocalization posterior probabilities PP.H0-PP.H4. +#' @title Colocalization Pipeline (coloc.bf_bf over paired LBF matrices) +#' @description Per-region pipeline that pairs two fine-mapping result +#' collections and runs \code{coloc::coloc.bf_bf} per (first-side tuple, +#' second-side tuple) pair to produce per-pair colocalization posterior +#' probabilities PP.H0-PP.H4. +#' +#' Either side may be a \code{\link{QtlFineMappingResult}} or a +#' \code{\link{GwasFineMappingResult}}: QTL-GWAS, QTL-QTL (two molecular +#' phenotypes) and GWAS-GWAS (two diseases) all run through one code path, +#' since nothing below the identity tuple depends on which flavour a side +#' is. The second side may also be handed in as summary statistics +#' (\code{\link{QtlSumStats}} or \code{\link{GwasSumStats}}), which are +#' fine-mapped inline. +#' +#' The argument names keep the QTL / GWAS wording of the common case. What +#' they mean generally is: \code{qtlFineMappingResult} is the side whose +#' identity is reported in the unprefixed \code{study} / \code{context} / +#' \code{trait} / \code{method} columns, and \code{gwasInput} is the side +#' reported in the \code{gwas}-prefixed ones. A GWAS side has no context or +#' trait axis, so those two columns are \code{NA} for it. #' #' @section Why \code{coloc.bf_bf} and not \code{coloc.susie}: #' The prior \code{colocWrapper} (now stubbed) used @@ -20,32 +33,37 @@ #' keeps effects at a secondary coverage; otherwise the default #' filter drops effects whose prior variance is below #' \code{priorTol}. -#' \item \strong{Multiple-GWAS batching}: when several GWAS -#' fine-mapping rows fall in the same region they are merged into -#' one combined LBF matrix per QTL pair (one \code{coloc.bf_bf} -#' call covers them all). +#' \item \strong{Per-tuple LBF reuse}: each second-side tuple's LBF +#' matrix is extracted once and scored against every first-side +#' tuple, so the filtering above is applied once per tuple rather +#' than once per pair. #' } #' This pipeline preserves all three. #' -#' GWAS input dispatch: +#' Second-side input dispatch: #' \itemize{ -#' \item \code{gwasInput} is a \code{\link{GwasSumStats}}: GWAS -#' fine-mapping is performed inline by +#' \item \code{gwasInput} is a \code{\link{QtlSumStats}} or a +#' \code{\link{GwasSumStats}}: it is fine-mapped inline by #' \code{\link{fineMappingPipeline}} with the supplied #' \code{finemappingMethods} (default \code{"susie"}). -#' \item \code{gwasInput} is a \code{\link{GwasFineMappingResult}}: -#' used directly; no inline fine-mapping. +#' \item \code{gwasInput} is a fine-mapping result: used directly; no +#' inline fine-mapping. #' } #' #' @section LD-sketch identity check: If #' \code{getLdSketch(qtlFineMappingResult)} is non-\code{NULL}, it must match -#' the LD sketch on \code{gwasInput}. Mismatch is a hard error. When the QTL -#' FMR's \code{ldSketch} is \code{NULL} (individual-level fit), the validation -#' is skipped on the QTL side. +#' the LD sketch on \code{gwasInput}. Mismatch is a hard error. When the +#' first side's \code{ldSketch} is \code{NULL} (individual-level fit), the +#' validation is skipped on that side and the second side's panel is what the +#' result carries forward. #' -#' @param qtlFineMappingResult A \code{\link{QtlFineMappingResult}} (required). -#' @param gwasInput Either a \code{\link{GwasSumStats}} or a -#' \code{\link{GwasFineMappingResult}}. +#' @param qtlFineMappingResult The first side: a +#' \code{\link{QtlFineMappingResult}} or a +#' \code{\link{GwasFineMappingResult}} (required). +#' @param gwasInput The second side: a \code{\link{QtlFineMappingResult}}, a +#' \code{\link{GwasFineMappingResult}}, or the summary statistics to +#' fine-map inline (\code{\link{QtlSumStats}} or +#' \code{\link{GwasSumStats}}). #' @param filterLbfCs Logical. When \code{TRUE} (and \code{filterLbfCsSecondary} #' is \code{NULL}), keep only effects that produced a credible set #' (\code{trimmedFit$sets$cs_index}). Default \code{FALSE}. @@ -68,15 +86,17 @@ #' @param p12 Prior probability of shared signal per variant. Default #' \code{5e-6}. #' @param finemappingMethods Character vector forwarded to -#' \code{\link{fineMappingPipeline}} when \code{gwasInput} is a -#' \code{GwasSumStats}. Default \code{"susie"}. -#' @param returnGwasFineMapping Logical. When \code{TRUE}, attach the computed -#' \code{GwasFineMappingResult} on the returned data frame as attribute -#' \code{"gwasFineMapping"}. Default \code{FALSE}. -#' @param enrichment Optional data.frame of per-(gwasStudy, qtlStudy, -#' qtlContext) enrichment factors with columns \code{gwasStudy}, -#' \code{qtlStudy}, \code{qtlContext}, \code{enrichment}. Output of -#' \code{\link{qtlEnrichmentPipeline}}. When non-\code{NULL}, each pair's +#' \code{\link{fineMappingPipeline}} when \code{gwasInput} is summary +#' statistics rather than a fine-mapping result. Default \code{"susie"}. +#' @param returnGwasFineMapping Logical. When \code{TRUE}, attach the +#' fine-mapping result computed from \code{gwasInput} on the returned object +#' as attribute \code{"gwasFineMapping"}. Default \code{FALSE}. +#' @param enrichment Optional data.frame of per-pair enrichment factors with +#' columns \code{gwasStudy}, \code{qtlStudy}, \code{qtlContext}, +#' \code{enrichment}, and optionally \code{gwasContext} / \code{gwasTrait} +#' (which the join uses when present, and which +#' \code{\link{qtlEnrichmentPipeline}} emits for a QTL outcome side). Output +#' of \code{\link{qtlEnrichmentPipeline}}. When non-\code{NULL}, each pair's #' \code{p12} prior is scaled to \code{min(p12 * (1 + enrichment), p12Max)} #' (the enrichment-informed colocalization variant, "enloc"). Pairs without a #' matching enrichment row fall back to the baseline \code{p12} with a @@ -105,11 +125,12 @@ #' is coding-invariant, so no sign change is needed); when FALSE, match on #' exact alleles only, so a ref/alt swap is treated as a distinct variant. #' @param ... Additional arguments forwarded to \code{coloc::coloc.bf_bf}. -#' @return A \code{\linkS4class{ColocResult}}: one element per tested (QTL -#' credible set, GWAS credible set, block) pair, holding that pair's aligned -#' variants with their \code{SNP.PP.H4}. Pair-level metadata carries the -#' identity columns (\code{study}, \code{context}, \code{trait}, -#' \code{method}, \code{gwasStudy}, \code{gwasMethod}), the block and stable +#' @return A \code{\linkS4class{ColocResult}}: one element per tested +#' (first-side credible set, second-side credible set, block) pair, holding +#' that pair's aligned variants with their \code{SNP.PP.H4}. Pair-level +#' metadata carries the identity columns (\code{study}, \code{context}, +#' \code{trait}, \code{method}, \code{gwasStudy}, \code{gwasContext}, +#' \code{gwasTrait}, \code{gwasMethod}), the block and stable #' credible-set ids (\code{blockId}, \code{qtlCs}, \code{gwasCs}), the #' standard coloc fields (\code{idx1}, \code{idx2}, \code{nSnps}, #' \code{hit1}, \code{hit2}, \code{PP.H0.abf} \ldots \code{PP.H4.abf}) and @@ -128,6 +149,12 @@ #' data(gwasFineMappingLbfExample) #' colocPipeline(qtlFineMappingLbfExample, #' gwasInput = gwasFineMappingLbfExample) +#' # Either side may be a QTL result. Pairing this collection against itself +#' # colocalizes its two contexts, and reports the second side's context and +#' # trait in gwasContext / gwasTrait. +#' res <- colocPipeline(qtlFineMappingLbfExample, +#' gwasInput = qtlFineMappingLbfExample) +#' unique(getColocPairs(res)[, c("context", "gwasContext")]) #' @export colocPipeline <- function( qtlFineMappingResult, @@ -194,19 +221,22 @@ colocPipeline <- function( abort(msg) # nocov end } - if (!methods::is(p$qtlFineMappingResult, "QtlFineMappingResult")) { + if (!methods::is(p$qtlFineMappingResult, "FineMappingResultBase")) { msg <- glue( - "`qtlFineMappingResult` must be a QtlFineMappingResult ", + "`qtlFineMappingResult` must be a QtlFineMappingResult or a ", + "GwasFineMappingResult ", "(got class '{class(p$qtlFineMappingResult)[[1L]]}')." ) abort(msg) } if ( - !methods::is(p$gwasInput, "GwasSumStats") && - !methods::is(p$gwasInput, "GwasFineMappingResult") + !methods::is(p$gwasInput, "SumStatsBase") && + !methods::is(p$gwasInput, "FineMappingResultBase") ) { msg <- glue( - "`gwasInput` must be a GwasSumStats or a GwasFineMappingResult ", + "`gwasInput` must be a fine-mapping result ", + "(QtlFineMappingResult / GwasFineMappingResult) or summary ", + "statistics (QtlSumStats / GwasSumStats) ", "(got class '{class(p$gwasInput)[[1L]]}')." ) abort(msg) @@ -235,20 +265,41 @@ colocPipeline <- function( ) abort(msg) } + .colocValidateEnrichmentKeys(enrichment) invisible(NULL) } -# Resolve the GWAS side to a GwasFineMappingResult (fine-map QC'd sumstats when -# a GwasSumStats is passed). +# One factor per pair: rows that repeat the identity the lookup joins on would +# make the applied enrichment depend on row order, so they are refused here +# rather than resolved by taking the first. +# @noRd +.colocValidateEnrichmentKeys <- function(enrichment) { + idCols <- intersect( + c("gwasStudy", "gwasContext", "gwasTrait", "qtlStudy", "qtlContext"), + colnames(enrichment) + ) + ids <- select(as_tibble(enrichment), all_of(idCols)) + if (nrow(distinct(ids)) == nrow(ids)) { + return(invisible(NULL)) + } + msg <- glue( + "`enrichment` has repeated ({str_flatten(idCols, ', ')}) rows; ", + "each pair needs exactly one enrichment factor." + ) + abort(msg) +} + +# Resolve the second side to a fine-mapping collection (fine-map QC'd sumstats +# when summary statistics are passed, whichever flavour they are). # @noRd .colocResolveGwasFmr <- function(gwasInput, finemappingMethods) { - if (methods::is(gwasInput, "GwasFineMappingResult")) { + if (methods::is(gwasInput, "FineMappingResultBase")) { return(gwasInput) } if (length(getQcInfo(gwasInput)) == 0L) { msg <- glue( - "colocPipeline: gwasInput (GwasSumStats) has no QC record. ", - "Call summaryStatsQc() first." + "colocPipeline: gwasInput ({class(gwasInput)[[1L]]}) has no QC ", + "record. Call summaryStatsQc() first." ) abort(msg) } @@ -278,11 +329,13 @@ colocPipeline <- function( # The LD reference the result carries forward, so getColocCredibleSets() can # recompute purity (section 3.7) without being handed a sketch separately. The -# QTL and GWAS sketches are already required to match by -# .colocRequireMatchingLdSketches, so either one identifies the panel. +# two sides' sketches are already required to match by +# .colocRequireMatchingLdSketches, so either one identifies the panel -- but a +# first side fit on individual-level data carries none, and then the second +# side's panel is the only one there is. # @noRd .colocLdSketch <- function(p) { - getLdSketch(p$qtlFineMappingResult) + getLdSketch(p$qtlFineMappingResult) %||% getLdSketch(p$gwasFmr) } # Empty-result early return (attaching the GWAS fine-mapping when requested). @@ -292,7 +345,7 @@ colocPipeline <- function( enriched = p$useEnrichment, ldSketch = .colocLdSketch(p) ) - if (p$returnGwasFineMapping && methods::is(p$gwasInput, "GwasSumStats")) { + if (p$returnGwasFineMapping && methods::is(p$gwasInput, "SumStatsBase")) { attr(out, "gwasFineMapping") <- p$gwasFmr } out @@ -316,7 +369,7 @@ colocPipeline <- function( q$retainedMass <- qLbfInfo$retainedMass q$effect <- qLbfInfo$effect compact(map( - names(p$gwasLbfByPair), + p$gwasLbfByPair, .colocScorePairAt, qLbfInfo = qLbfInfo, p = p, @@ -324,27 +377,32 @@ colocPipeline <- function( )) } -# Identity + row payload + log label for a QTL tuple. +# Identity + row payload + log label for one first-side tuple. # @noRd .colocQtlTupleInfo <- function(qi, p) { fmr <- p$qtlFineMappingResult - study <- as.character(fmr$study)[[qi]] - context <- as.character(fmr$context)[[qi]] - trait <- as.character(fmr$trait)[[qi]] - method <- as.character(fmr$method)[[qi]] - list( - study = study, - context = context, - trait = trait, - method = method, - parts = .fmrRowParts(fmr, qi), - label = glue( - "QTL (study='{study}', context='{context}', ", - "trait='{trait}', method='{method}')" + ident <- .colocTupleIdentity(fmr, qi) + c( + ident, + list( + parts = .fmrRowParts(fmr, qi), + label = .fmrTupleLabel(.fmrSideName(fmr), ident) ) ) } +# The identity tuple of one fine-mapping row, whichever flavour the collection +# is (see .fmrIdentityAt for the absent-axis rule). +# @noRd +.colocTupleIdentity <- function(fmr, ri) { + list( + study = .fmrIdentityAt(fmr, "study", ri), + context = .fmrIdentityAt(fmr, "context", ri), + trait = .fmrIdentityAt(fmr, "trait", ri), + method = .fmrIdentityAt(fmr, "method", ri) + ) +} + # Score one (QTL, GWAS) pair via coloc.bf_bf -> a summary row, or NULL when the # variants don't align or coloc fails / returns no summary. # @noRd @@ -355,7 +413,7 @@ colocPipeline <- function( if (is.null(aligned)) { return(NULL) } - p12Info <- .colocResolveP12(p, gInfo$study, q$study, q$context) + p12Info <- .colocResolveP12(p, gInfo, q) pairRes <- .colocRunPair(aligned, p, p12Info$p12Used, q, gInfo) if (is.null(pairRes) || is.null(pairRes$summary)) { return(NULL) @@ -373,16 +431,16 @@ colocPipeline <- function( # Enrichment-informed p12 (per-(gwasStudy, qtlStudy, qtlContext) scaling capped # at p12Max; baseline p12 with no enrichment table / no matching row). # @noRd -.colocResolveP12 <- function(p, gwasStudy, qStudy, qContext) { +.colocResolveP12 <- function(p, gInfo, q) { if (!p$useEnrichment) { return(list(enRow = NA_real_, p12Used = p$p12)) } - enRow <- .colocLookupEnrichment(p$enrichment, gwasStudy, qStudy, qContext) + enRow <- .colocLookupEnrichment(p$enrichment, gInfo, q) if (is.na(enRow)) { msg <- glue( "colocPipeline: no enrichment entry for ", - "(gwasStudy='{gwasStudy}', qtlStudy='{qStudy}', ", - "qtlContext='{qContext}'); using baseline p12." + "(gwasStudy='{gInfo$study}', qtlStudy='{q$study}', ", + "qtlContext='{q$context}'); using baseline p12." ) warn(msg) enRow <- 0 @@ -407,10 +465,8 @@ colocPipeline <- function( exec(coloc::coloc.bf_bf, !!!colocArgs), error = function(e) { msg <- glue( - "colocPipeline: coloc.bf_bf failed for QTL ", - "(study='{q$study}', context='{q$context}', ", - "trait='{q$trait}', method='{q$method}') x GWAS ", - "(study='{gInfo$study}', method='{gInfo$method}'): ", + "colocPipeline: coloc.bf_bf failed for ", + "{q$label} x {gInfo$label}: ", "{conditionMessage(e)}" ) warn(msg) @@ -429,6 +485,8 @@ colocPipeline <- function( sm$trait <- q$trait sm$method <- q$method sm$gwasStudy <- gInfo$study + sm$gwasContext <- gInfo$context + sm$gwasTrait <- gInfo$trait sm$gwasMethod <- gInfo$method # idx1 / idx2 index the LBF rows handed to coloc.bf_bf, which is exactly # what retainedMass runs parallel to -- so the mass reported here is the @@ -471,7 +529,7 @@ colocPipeline <- function( # @noRd .colocFinalize <- function(results, p) { out <- .colocAssemble(results, p$useEnrichment, .colocLdSketch(p)) - if (p$returnGwasFineMapping && methods::is(p$gwasInput, "GwasSumStats")) { + if (p$returnGwasFineMapping && methods::is(p$gwasInput, "SumStatsBase")) { attr(out, "gwasFineMapping") <- p$gwasFmr } out @@ -497,6 +555,8 @@ colocPipeline <- function( "trait", "method", "gwasStudy", + "gwasContext", + "gwasTrait", "gwasMethod", "blockId", "qtlCs", @@ -710,10 +770,14 @@ colocPipeline <- function( lbfMatrix[, !is.na(colnames(lbfMatrix)), drop = FALSE] } -# Build a per-GWAS-tuple LBF matrix list, keyed by "study|method". -# Within each key we stack multiple FMR rows row-wise (the legacy -# "combined GWAS LBF" pattern), drop NA columns, and replace NAs with -# 0 so a fresh QTL pairing always lands on the same coordinate frame. +# The second side's LBF matrices, one record per row of the collection, each +# carrying the row's identity so the pair it scores can be named. +# +# One record per ROW rather than per identity key: a QTL second side has many +# rows sharing (study, method, block) and differing only on context / trait, so +# keying on the GWAS 2-tuple would let one trait's matrix silently replace +# another's. Nothing downstream indexes this list by name, so positional +# records make the collision impossible instead of merely unlikely. # @noRd .colocPreextractGwasLbf <- function( gwasFmr, @@ -722,39 +786,59 @@ colocPipeline <- function( filterLbfCsConcentration, priorTol ) { - keys <- str_c( - as.character(gwasFmr$study), - as.character(gwasFmr$method), - as.character(.colocGwasBlockIds(gwasFmr)), - sep = "||" + if (nrow(gwasFmr) == 0L) { + return(list()) + } + compact(map( + seq_len(nrow(gwasFmr)), + .colocGwasLbfAt, + gwasFmr = gwasFmr, + blockIds = .colocGwasBlockIds(gwasFmr), + side = .fmrSideName(gwasFmr), + filterLbfCs = filterLbfCs, + filterLbfCsSecondary = filterLbfCsSecondary, + filterLbfCsConcentration = filterLbfCsConcentration, + priorTol = priorTol + )) +} + +# One second-side row's LBF matrix plus its identity, or NULL when the row has +# no usable LBF. +# @noRd +.colocGwasLbfAt <- function( + ri, + gwasFmr, + blockIds, + side, + filterLbfCs, + filterLbfCsSecondary, + filterLbfCsConcentration, + priorTol +) { + ident <- .colocTupleIdentity(gwasFmr, ri) + blockId <- blockIds[[ri]] + label <- .fmrTupleLabel(side, ident, block = blockId) + info <- .colocExtractLbfFromEntry( + .fmrRowParts(gwasFmr, ri), + filterLbfCs, + filterLbfCsSecondary, + filterLbfCsConcentration, + priorTol, + label = label ) - out <- list() - for (ri in seq_len(nrow(gwasFmr))) { - parts <- str_split(keys[[ri]], "\\|\\|")[[1L]] - info <- .colocExtractLbfFromEntry( - .fmrRowParts(gwasFmr, ri), - filterLbfCs, - filterLbfCsSecondary, - filterLbfCsConcentration, - priorTol, - label = glue( - "GWAS (study='{parts[[1L]]}', method='{parts[[2L]]}', ", - "block='{parts[[3L]]}')" - ) - ) - if (is.null(info)) { - next - } - out[[keys[[ri]]]] <- list( + if (is.null(info)) { + return(NULL) + } + c( + ident, + list( lbf = info$lbf, retainedMass = info$retainedMass, effect = info$effect, - study = parts[[1L]], - method = parts[[2L]], - blockId = parts[[3L]] + blockId = blockId, + label = label ) - } - out + ) } # The LD block each GWAS fine-mapping row was computed on. The element's own @@ -807,6 +891,8 @@ colocPipeline <- function( trait = character(0), method = character(0), gwasStudy = character(0), + gwasContext = character(0), + gwasTrait = character(0), gwasMethod = character(0), blockId = character(0), qtlCs = integer(0), @@ -829,28 +915,54 @@ colocPipeline <- function( ColocResult(base, list(), ldSketch = ldSketch) } -# Look up the enrichment factor for a (gwasStudy, qtlStudy, qtlContext) -# triple in the user-supplied enrichment table. Returns NA when the -# triple is not present; the caller falls back to the baseline p12 and -# emits a warning. +# Look up this pair's enrichment factor in the user-supplied enrichment table. +# +# The join uses whichever identity columns the table carries: the +# (gwasStudy, qtlStudy, qtlContext) triple always, plus gwasContext / +# gwasTrait when qtlEnrichmentPipeline ran with a QTL outcome side. Without +# those two, one study's molecular traits are indistinguishable in the table +# and every one of them would take the same row. Returns NA when the pair is +# not present; the caller falls back to the baseline p12 and warns. # @noRd -.colocLookupEnrichment <- function( - enrichment, - gwasStudy, - qtlStudy, - qtlContext -) { - idx <- which( - as.character(enrichment$gwasStudy) == gwasStudy & - as.character(enrichment$qtlStudy) == qtlStudy & - as.character(enrichment$qtlContext) == qtlContext +.colocLookupEnrichment <- function(enrichment, gInfo, q) { + wanted <- .colocEnrichmentKey(enrichment, gInfo, q) + hits <- map( + names(wanted), + .colocEnrichmentColumnMatches, + enrichment = enrichment, + wanted = wanted ) + idx <- which(reduce(hits, `&`)) if (length(idx) == 0L) { return(NA_real_) } as.numeric(enrichment$enrichment[[idx[[1L]]]]) } +# The identity a pair is looked up by, narrowed to the columns the table has. +# @noRd +.colocEnrichmentKey <- function(enrichment, gInfo, q) { + wanted <- list( + gwasStudy = gInfo$study, + gwasContext = gInfo$context, + gwasTrait = gInfo$trait, + qtlStudy = q$study, + qtlContext = q$context + ) + wanted[is_in(names(wanted), colnames(enrichment))] +} + +# One key column's row match. An axis neither side has is NA on both, and +# matches -- `==` would evaluate to NA there and drop every row. +# @noRd +.colocEnrichmentColumnMatches <- function(column, enrichment, wanted) { + values <- as.character(enrichment[[column]]) + if (is.na(wanted[[column]])) { + return(is.na(values)) + } + !is.na(values) & values == wanted[[column]] +} + # Ensure each row data.frame from coloc.bf_bf carries the standard PP # columns even when the underlying call produced a slightly different # shape. @@ -901,10 +1013,11 @@ colocPipeline <- function( # The variant ids of one fine-mapping entry (S4 slot; not pluckable by name). # @noRd -# Score the QTL LBF against GWAS pair `gKey` -> a summary row (or NULL). +# Score the first side's LBF against one second-side record -> a summary row +# (or NULL). # @noRd -.colocScorePairAt <- function(gKey, qLbfInfo, p, q) { - .colocScorePair(qLbfInfo$lbf, p$gwasLbfByPair[[gKey]], q, p) +.colocScorePairAt <- function(gInfo, qLbfInfo, p, q) { + .colocScorePair(qLbfInfo$lbf, gInfo, q, p) } # TRUE when a credible set has fewer than `maxSize` variants. diff --git a/R/qtlEnrichmentPipeline.R b/R/qtlEnrichmentPipeline.R index a2460675..bcaf5a7a 100644 --- a/R/qtlEnrichmentPipeline.R +++ b/R/qtlEnrichmentPipeline.R @@ -1,32 +1,43 @@ #' @title QTL Enrichment Pipeline (Genome-Wide) -#' @description Genome-wide pipeline that computes per-pair (GWAS study, QTL -#' context) enrichment estimates by passing the GWAS PIP vector and the QTL +#' @description Genome-wide pipeline that computes per-pair enrichment +#' estimates by passing an outcome PIP vector and a set of annotation #' credible-set posteriors to \code{\link{qtlEnrichment}}. The returned table #' feeds \code{\link{colocPipeline}} via its \code{enrichment} argument. #' +#' Either side may be a \code{\link{QtlFineMappingResult}} or a +#' \code{\link{GwasFineMappingResult}}, so QTL-in-GWAS, QTL-in-QTL and +#' GWAS-in-GWAS enrichment all run through one code path. The argument names +#' keep the QTL / GWAS wording of the common case; what they mean generally +#' is that \code{gwasFineMappingResult} is the \strong{outcome} whose PIPs +#' are scanned and \code{qtlFineMappingResult} is the \strong{annotation} +#' whose region fits are tested for enrichment within them. +#' #' \strong{Not gene-parallelisable}: the enrichment estimator runs over the full -#' genome of GWAS PIPs and the full collection of QTL fits at once. +#' genome of outcome PIPs and the full collection of annotation fits at once. #' #' @section Inputs: #' \itemize{ -#' \item \code{gwasFineMappingResult}: a genome-wide -#' \code{\link{GwasFineMappingResult}} (one row per (study, LD -#' block) tuple). Each entry's \code{FineMappingRow$trimmedFit} -#' must carry a \code{pip} vector. -#' \item \code{qtlFineMappingResult}: the genome-wide -#' \code{\link{QtlFineMappingResult}}. Each entry's +#' \item \code{gwasFineMappingResult}: the genome-wide outcome collection. +#' Each entry's \code{FineMappingRow$trimmedFit} must carry a \code{pip} +#' vector. One PIP vector is built per outcome trait -- keyed by +#' \code{study} for a GWAS collection (one trait per study) and by +#' (\code{study}, \code{context}, \code{trait}) for a QTL one, since a +#' variant's PIP differs between molecular traits and pooling them would +#' collide. +#' \item \code{qtlFineMappingResult}: the annotation collection. Each entry's #' \code{trimmedFit} must carry \code{alpha}, \code{pip}, and -#' prior-variance fields (\code{V}). +#' prior-variance fields (\code{V}). Its region fits are pooled per +#' (\code{study}, \code{context}) -- per \code{study} alone for a GWAS +#' collection, which has no context axis. #' } #' -#' @section LD-sketch identity check: The GWAS \code{FineMappingResultBase} must -#' have a non-NULL \code{ldSketch} (RSS-derived). If the QTL FMR also has a -#' non-NULL \code{ldSketch}, the two must match exactly. When the QTL FMR's -#' \code{ldSketch} is NULL (individual-level QTL fit), validation is skipped -#' on the QTL side. +#' @section LD-sketch identity check: A GWAS outcome collection must have a +#' non-NULL \code{ldSketch} (it should be RSS-derived). Where both sides +#' carry one, the two must match exactly; a \code{NULL} on either side (an +#' individual-level fit) skips the check. #' -#' @param gwasFineMappingResult See above. -#' @param qtlFineMappingResult See above. +#' @param gwasFineMappingResult The outcome side; see above. +#' @param qtlFineMappingResult The annotation side; see above. #' @param numGwas Number of GWAS variants used to estimate \code{piGwas}. When #' \code{NULL} (default) it is estimated from the data -- bias warning applies #' if the input PIP vector is not genome-wide. @@ -42,11 +53,26 @@ #' \code{\link{qtlEnrichment}} for reproducible multiple imputation. #' \code{NULL} (default) draws a nondeterministic seed. #' @param ... Additional arguments forwarded to \code{\link{qtlEnrichment}}. -#' @return A tibble with one row per (gwasStudy, qtlStudy, qtlContext) -#' triple and columns \code{gwasStudy}, \code{qtlStudy}, \code{qtlContext}, -#' \code{enrichment}, \code{enrichmentSe}, \code{enrichmentLogOdds}, plus any -#' extras the underlying estimator emits. Suitable as the \code{enrichment} -#' argument to \code{\link{colocPipeline}} (which joins on the same triple). +#' @return A tibble with one row per (outcome trait, annotation unit) pair. +#' The identity columns are \code{gwasStudy}, \code{gwasContext}, +#' \code{gwasTrait}, \code{qtlStudy}, \code{qtlContext}; the axes a side does +#' not have are \code{NA} (\code{gwasContext} / \code{gwasTrait} for a GWAS +#' outcome, \code{qtlContext} for a GWAS annotation). Suitable as the +#' \code{enrichment} argument to \code{\link{colocPipeline}}, which joins on +#' those columns. +#' +#' The estimates are \code{enrichmentLogOdds}, the enrichment parameter +#' \eqn{a_1} on the log-odds scale, with its standard error +#' \code{enrichmentSe}; \code{enrichment} is the same quantity as a +#' multiplicative factor, \eqn{e^{a_1} - 1}, which is what +#' \code{colocPipeline} scales \code{p12} by (so \eqn{a_1 = 0} leaves the +#' prior untouched). \code{enrichmentLogOddsNoShrinkage} and +#' \code{enrichmentSeNoShrinkage} are the same estimate before shrinkage, +#' \code{intercept} / \code{interceptSe} are \eqn{a_0}, and \code{colocP1}, +#' \code{colocP2}, \code{colocP12} are the enrichment-informed coloc priors +#' the estimator derives from \eqn{(a_0, a_1)} -- an alternative to scaling +#' a baseline \code{p12}. \code{effectiveMiRounds} is how many +#' multiple-imputation rounds survived outlier filtering. #' @examples #' data(gwasFineMappingExample) #' data(qtlFineMappingExample) @@ -70,56 +96,79 @@ qtlEnrichmentPipeline <- function( p <- as.list(environment()) p$dots <- list(...) p <- .enrPrepare(p) - results <- list_flatten(map(p$gwasStudies, .enrScoreGwasStudy, p = p)) + results <- list_flatten(map( + seq_len(nrow(p$gwasTuples)), + .enrScoreOutcomeTuple, + p = p + )) .enrAssemble(results) } # Validate the input classes + LD-sketch presence / identity. # @noRd .enrValidateInputs <- function(gwasFineMappingResult, qtlFineMappingResult) { - if (!methods::is(gwasFineMappingResult, "GwasFineMappingResult")) { - abort("`gwasFineMappingResult` must be a GwasFineMappingResult.") + if (!methods::is(gwasFineMappingResult, "FineMappingResultBase")) { + msg <- glue( + "`gwasFineMappingResult` must be a GwasFineMappingResult or a ", + "QtlFineMappingResult ", + "(got class '{class(gwasFineMappingResult)[[1L]]}')." + ) + abort(msg) } - if (!methods::is(qtlFineMappingResult, "QtlFineMappingResult")) { - abort("`qtlFineMappingResult` must be a QtlFineMappingResult.") + if (!methods::is(qtlFineMappingResult, "FineMappingResultBase")) { + msg <- glue( + "`qtlFineMappingResult` must be a QtlFineMappingResult or a ", + "GwasFineMappingResult ", + "(got class '{class(qtlFineMappingResult)[[1L]]}')." + ) + abort(msg) } - gwasLd <- getLdSketch(gwasFineMappingResult) - if (is.null(gwasLd)) { + outcomeLd <- getLdSketch(gwasFineMappingResult) + if ( + is.null(outcomeLd) && + methods::is(gwasFineMappingResult, "GwasFineMappingResult") + ) { msg <- glue( "qtlEnrichmentPipeline: the GWAS FineMappingResult must have a ", "non-NULL ldSketch (it should be RSS-derived)." ) abort(msg) } - .colocRequireMatchingLdSketches(getLdSketch(qtlFineMappingResult), gwasLd) + # Lenient rather than qtl-required: a QTL outcome side may be an + # individual-level fit carrying no panel, and the requirement above already + # covers the RSS-derived GWAS case. + .requireMatchingLdSketches( + getLdSketch(qtlFineMappingResult), + outcomeLd, + pipelineName = "qtlEnrichmentPipeline", + nullPolicy = "lenient" + ) invisible(NULL) } -# Hoist the GWAS-study-independent work out of the double loop: per-study GWAS -# PIP vectors, the union variant-name panel, per-tuple QTL regions, and each -# tuple's one-time alignment to the union panel (errors captured as values). +# Hoist the outcome-independent work out of the double loop: per-trait outcome +# PIP vectors, the union variant-name panel, per-tuple annotation regions, and +# each tuple's one-time alignment to the union panel (errors captured as +# values). # @noRd .enrPrepare <- function(p) { - p$gwasStudies <- unique(as.character(p$gwasFineMappingResult$study)) - # Iterate the joint (study, context) key: context alone would silently merge - # two studies sharing a context label, giving wrong enrichment estimates. - p$qtlTuples <- distinct(tibble( - qtlStudy = as.character(p$qtlFineMappingResult$study), - qtlContext = as.character(p$qtlFineMappingResult$context) - )) - if (length(p$gwasStudies) == 0L || nrow(p$qtlTuples) == 0L) { + p$gwasTuples <- .enrOutcomeTuples(p$gwasFineMappingResult) + p$qtlTuples <- .enrAnnotationTuples(p$qtlFineMappingResult) + if (nrow(p$gwasTuples) == 0L || nrow(p$qtlTuples) == 0L) { msg <- glue( - "qtlEnrichmentPipeline: no (gwasStudy, qtlStudy, qtlContext) ", - "triples to compute (one of the inputs has zero rows)." + "qtlEnrichmentPipeline: no (outcome, annotation) pairs to ", + "compute (one of the inputs has zero rows)." ) abort(msg) } - p$gwasPipByStudy <- .enrGwasPipByStudy( - p$gwasFineMappingResult, - p$gwasStudies + p$gwasPipByTuple <- map( + seq_len(nrow(p$gwasTuples)), + .enrGwasPipForRow, + gwasTuples = p$gwasTuples, + fmr = p$gwasFineMappingResult ) unionGwasNames <- unique(unlist( - map(p$gwasPipByStudy, names), + map(p$gwasPipByTuple, names), use.names = FALSE )) p$qtlRegionsByTuple <- map( @@ -136,28 +185,77 @@ qtlEnrichmentPipeline <- function( p } -# Per-study genome-wide GWAS PIP vectors (named by variant id). A loop rather -# than map because the varying study is `.enrBuildGwasPipVector`'s second arg. +# The outcome side's per-trait keys: one PIP vector is built per key. A GWAS +# collection carries one trait per study, so its QTL-only axes read NA and the +# key collapses to (study). # @noRd -.enrGwasPipByStudy <- function(gwasFineMappingResult, gwasStudies) { - out <- vector("list", length(gwasStudies)) - for (i in seq_along(gwasStudies)) { - out[[i]] <- .enrBuildGwasPipVector( - gwasFineMappingResult, - gwasStudies[[i]] - ) +.enrOutcomeTuples <- function(fmr) { + distinct(tibble( + gwasStudy = .fmrIdentityColumn(fmr, "study"), + gwasContext = .fmrIdentityColumn(fmr, "context"), + gwasTrait = .fmrIdentityColumn(fmr, "trait") + )) +} + +# The annotation side's keys: (study, context), the unit whose region fits are +# pooled into one enrichment estimate. The joint key matters -- context alone +# would silently merge two studies sharing a context label -- and a GWAS +# annotation pools its blocks under one NA-context key. +# @noRd +.enrAnnotationTuples <- function(fmr) { + distinct(tibble( + qtlStudy = .fmrIdentityColumn(fmr, "study"), + qtlContext = .fmrIdentityColumn(fmr, "context") + )) +} + +# One tuple-table row as an identity list addressed by the collection's OWN +# column names; the table prefixes them so the result can name both sides. +# @noRd +.enrOutcomeIdent <- function(gwasTuples, k) { + list( + study = gwasTuples$gwasStudy[[k]], + context = gwasTuples$gwasContext[[k]], + trait = gwasTuples$gwasTrait[[k]] + ) +} + +# @noRd +.enrAnnotationIdent <- function(qtlTuples, k) { + list( + study = qtlTuples$qtlStudy[[k]], + context = qtlTuples$qtlContext[[k]] + ) +} + +# Row indices of `fmr` matching an identity tuple. An axis the collection does +# not have is NA on both sides and matches, rather than excluding every row. +# @noRd +.enrMatchRows <- function(fmr, ident) { + hits <- map(names(ident), .enrColumnMatches, fmr = fmr, ident = ident) + which(reduce(hits, `&`, .init = rep(TRUE, nrow(fmr)))) +} + +# @noRd +.enrColumnMatches <- function(column, fmr, ident) { + values <- .fmrIdentityColumn(fmr, column) + wanted <- ident[[column]] + if (is.na(wanted)) { + return(is.na(values)) } - set_names(out, gwasStudies) + !is.na(values) & values == wanted } -# QTL SuSiE region list for the k-th (study, context) tuple. +# The outcome PIP vector for the k-th outcome tuple. +# @noRd +.enrGwasPipForRow <- function(k, gwasTuples, fmr) { + .enrBuildGwasPipVector(fmr, .enrOutcomeIdent(gwasTuples, k)) +} + +# Annotation SuSiE region list for the k-th (study, context) tuple. # @noRd .enrQtlRegionsForRow <- function(k, qtlTuples, fmr) { - .enrBuildQtlRegionsList( - fmr, - qtlTuples$qtlStudy[[k]], - qtlTuples$qtlContext[[k]] - ) + .enrBuildQtlRegionsList(fmr, .enrAnnotationIdent(qtlTuples, k)) } # Align one tuple's regions to the union GWAS panel, capturing any error as a @@ -170,15 +268,15 @@ qtlEnrichmentPipeline <- function( ) } -# Score one GWAS study against every QTL tuple -> enrichment records (empty when -# the study has no usable PIPs). +# Score one outcome trait against every annotation tuple -> enrichment records +# (empty when the outcome has no usable PIPs). # @noRd -.enrScoreGwasStudy <- function(gStudy, p) { - gwasPip <- p$gwasPipByStudy[[gStudy]] +.enrScoreOutcomeTuple <- function(gi, p) { + gwasPip <- p$gwasPipByTuple[[gi]] if (length(gwasPip) == 0L) { msg <- glue( "qtlEnrichmentPipeline: no usable PIPs for ", - "gwasStudy='{gStudy}'; skipping." + "{.enrOutcomeLabel(p, gi)}; skipping." ) warn(msg) return(list()) @@ -186,41 +284,57 @@ qtlEnrichmentPipeline <- function( compact(map( seq_len(nrow(p$qtlTuples)), .enrScoreTuple, - gStudy = gStudy, + gi = gi, gwasPip = gwasPip, p = p )) } -# Score one (gwasStudy, qtl tuple) pair -> an enrichment record, or NULL when -# the tuple has no regions or qtlEnrichment fails. +# Score one (outcome trait, annotation tuple) pair -> an enrichment record, or +# NULL when the tuple has no regions or qtlEnrichment fails. # @noRd -.enrScoreTuple <- function(k, gStudy, gwasPip, p) { - qStudy <- p$qtlTuples$qtlStudy[[k]] - qContext <- p$qtlTuples$qtlContext[[k]] +.enrScoreTuple <- function(k, gi, gwasPip, p) { if (length(p$qtlRegionsByTuple[[k]]) == 0L) { msg <- glue( - "qtlEnrichmentPipeline: no usable QTL regions for ", - "(qtlStudy='{qStudy}', qtlContext='{qContext}'); skipping." + "qtlEnrichmentPipeline: no usable regions for ", + "{.enrAnnotationLabel(p, k)}; skipping." ) warn(msg) return(NULL) } - enr <- .enrRunEnrichment(gStudy, gwasPip, k, qStudy, qContext, p) + enr <- .enrRunEnrichment(gi, gwasPip, k, p) if (is.null(enr)) { return(NULL) } - row <- .enrFlattenEnrichment(enr) - row$gwasStudy <- gStudy - row$qtlStudy <- qStudy - row$qtlContext <- qContext - row + c( + .enrFlattenEnrichment(enr), + as.list(p$gwasTuples[gi, , drop = FALSE]), + as.list(p$qtlTuples[k, , drop = FALSE]) + ) +} + +# Human-readable identities for the warnings above, naming each side by its own +# flavour and only the axes it has. +# @noRd +.enrOutcomeLabel <- function(p, gi) { + .fmrTupleLabel( + .fmrSideName(p$gwasFineMappingResult), + .enrOutcomeIdent(p$gwasTuples, gi) + ) +} + +# @noRd +.enrAnnotationLabel <- function(p, k) { + .fmrTupleLabel( + .fmrSideName(p$qtlFineMappingResult), + .enrAnnotationIdent(p$qtlTuples, k) + ) } # Run qtlEnrichment for a pair (with the pre-aligned regions), warning + NULL on # failure. alignNames = FALSE reuses the shared per-tuple alignment. # @noRd -.enrRunEnrichment <- function(gStudy, gwasPip, k, qStudy, qContext, p) { +.enrRunEnrichment <- function(gi, gwasPip, k, p) { aligned <- p$alignedByTuple[[k]] tryCatch( { @@ -247,8 +361,8 @@ qtlEnrichmentPipeline <- function( eMsg <- conditionMessage(e) msg <- glue( "qtlEnrichmentPipeline: qtlEnrichment failed for ", - "(gwasStudy='{gStudy}', qtlStudy='{qStudy}', ", - "qtlContext='{qContext}'): {eMsg}" + "{.enrOutcomeLabel(p, gi)} x ", + "{.enrAnnotationLabel(p, k)}: {eMsg}" ) warn(msg) NULL @@ -263,21 +377,29 @@ qtlEnrichmentPipeline <- function( return(.enrEmptyResult()) } out <- bind_rows(results) - idCols <- c("gwasStudy", "qtlStudy", "qtlContext") - select(out, all_of(idCols), everything()) + select(out, all_of(.enrIdCols()), everything()) +} + +# The identity columns of a result row: the outcome trait's tuple, then the +# annotation unit's. +# @noRd +.enrIdCols <- function() { + c("gwasStudy", "gwasContext", "gwasTrait", "qtlStudy", "qtlContext") } -# The empty enrichment result table. +# The empty enrichment result table, with the same columns a populated one has. # @noRd .enrEmptyResult <- function() { - tibble( - gwasStudy = character(0), - qtlStudy = character(0), - qtlContext = character(0), - enrichment = numeric(0), - enrichmentSe = numeric(0), - enrichmentLogOdds = numeric(0) + idCols <- set_names( + rep(list(character(0)), length(.enrIdCols())), + .enrIdCols() ) + valueNames <- names(.enrNaEnrichment()) + valueCols <- set_names( + rep(list(numeric(0)), length(valueNames)), + valueNames + ) + as_tibble(c(idCols, valueCols)) } # ============================================================================= @@ -293,15 +415,16 @@ qtlEnrichmentPipeline <- function( map(regions, .enrAlignRegion, unionGwasNames = unionGwasNames) } -# Build a named GWAS PIP vector for one study. Walks every row of the -# GwasFineMappingResult tagged with that study, extracts the per-row -# pip from each FineMappingRow, and concatenates with variant-id -# names. Errors if any single variant appears with conflicting PIP -# values across rows. +# Build a named outcome PIP vector for one trait. Walks every row of the +# collection carrying that identity, extracts the per-row pip from each +# FineMappingRow, and concatenates with variant-id names. Errors if any single +# variant appears with conflicting PIP values across rows -- which is why the +# identity is the full trait tuple: pooling two molecular traits of one study +# would collide on every variant they share. #' @importFrom dplyr add_count #' @noRd -.enrBuildGwasPipVector <- function(gwasFmr, gStudy) { - idx <- which(as.character(gwasFmr$study) == gStudy) +.enrBuildGwasPipVector <- function(gwasFmr, ident) { + idx <- .enrMatchRows(gwasFmr, ident) if (length(idx) == 0L) { return(numeric(0)) } @@ -357,17 +480,14 @@ qtlEnrichmentPipeline <- function( set_names(byId$pip, byId$id) } -# Build the per-(qtlStudy, qtlContext) list of region fits in the shape -# that qtlEnrichment expects: list(d) where each d carries -# alpha, pip, prior_variance (V). Filters on BOTH study and context so -# entries from different studies that happen to share a context label -# are not pooled into one enrichment estimate. +# Build the per-(study, context) list of region fits in the shape that +# qtlEnrichment expects: list(d) where each d carries alpha, pip, +# prior_variance (V). Filters on BOTH study and context so entries from +# different studies that happen to share a context label are not pooled into +# one enrichment estimate. # @noRd -.enrBuildQtlRegionsList <- function(qtlFmr, qStudy, qContext) { - idx <- which( - as.character(qtlFmr$study) == qStudy & - as.character(qtlFmr$context) == qContext - ) +.enrBuildQtlRegionsList <- function(qtlFmr, ident) { + idx <- .enrMatchRows(qtlFmr, ident) if (length(idx) == 0L) { return(list()) } @@ -412,56 +532,67 @@ qtlEnrichmentPipeline <- function( } } -# Coerce qtlEnrichment's variable-shape output into a single-row -# named list with the canonical columns the caller documents. The -# underlying estimator returns either a list with named numeric scalars -# (enrichment, enrichmentSe, enrichmentLogOdds, ...) or a matrix/df -- -# this helper handles both. +# Project qtlEnrichment()'s output onto the columns this pipeline publishes. +# +# The field names are the estimator's own, written in src/qtl_enrichment.h and +# shared verbatim with upstream fastenloc's enloc.enrich.out, so they are +# matched literally rather than guessed at. The shrinkage estimates are the +# ones reported, matching upstream, which likewise switches to the shrunk a1 +# before deriving the coloc priors. +# +# `enrichment` is expm1 of the log-odds a1 rather than a1 itself, because +# colocPipeline consumes it as `p12 * (1 + enrichment)`, which is then exactly +# the enloc-adjusted prior `p12 * exp(a1)` -- and leaves p12 untouched for an +# unenriched annotation (a1 = 0). expm1 is bounded below by -1, so a depleted +# annotation shrinks p12 towards 0 rather than past it. # @noRd .enrFlattenEnrichment <- function(enr) { - if (is.list(enr) && is.null(dim(enr))) { - list( - enrichment = .enrPickScalar("enrichment", enr), - enrichmentSe = .enrPickScalar("enrichmentSe", enr), - enrichmentLogOdds = .enrPickScalar("enrichmentLogOdds", enr) - ) - } else if (is.matrix(enr) || is.data.frame(enr)) { - df <- as_tibble(enr, .name_repair = "minimal") - if (nrow(df) == 0L) { - list( - enrichment = NA_real_, - enrichmentSe = NA_real_, - enrichmentLogOdds = NA_real_ - ) - } else { - list( - enrichment = .enrPickColumn(df, c("enrichment", "Enrichment")), - enrichmentSe = .enrPickColumn( - df, - c("enrichmentSe", "se", "stderr") - ), - enrichmentLogOdds = .enrPickColumn( - df, - c("enrichmentLogOdds", "logOdds", "log_odds") - ) - ) - } - } else { - list( - enrichment = NA_real_, - enrichmentSe = NA_real_, - enrichmentLogOdds = NA_real_ + if (!is.list(enr) || !is_in("Enrichment (w/ shrinkage)", names(enr))) { + msg <- glue( + "qtlEnrichmentPipeline: the enrichment estimator returned no ", + "'Enrichment (w/ shrinkage)' field, so every estimate is ", + "reported as NA. Expected the field names written by ", + "qtl_enrichment.h." ) + warn(msg) + return(.enrNaEnrichment()) } + logOdds <- .enrPickScalar("Enrichment (w/ shrinkage)", enr) + list( + enrichment = expm1(logOdds), + enrichmentSe = .enrPickScalar("sd (w/ shrinkage)", enr), + enrichmentLogOdds = logOdds, + enrichmentLogOddsNoShrinkage = .enrPickScalar( + "Enrichment (no shrinkage)", + enr + ), + enrichmentSeNoShrinkage = .enrPickScalar("sd (no shrinkage)", enr), + intercept = .enrPickScalar("Intercept", enr), + interceptSe = .enrPickScalar("sd (intercept)", enr), + colocP1 = .enrPickScalar("Alternative (coloc) p1", enr), + colocP2 = .enrPickScalar("Alternative (coloc) p2", enr), + colocP12 = .enrPickScalar("Alternative (coloc) p12", enr), + effectiveMiRounds = .enrPickScalar("Effective MI rounds", enr) + ) } +# The value columns, all unmeasured. Also the single source of the value-column +# schema, so the empty result cannot drift from the populated one. # @noRd -.enrPickColumn <- function(df, candidates) { - hit <- intersect(candidates, colnames(df)) - if (length(hit) == 0L) { - return(NA_real_) - } - as.numeric(df[[hit[[1L]]]][[1L]]) +.enrNaEnrichment <- function() { + list( + enrichment = NA_real_, + enrichmentSe = NA_real_, + enrichmentLogOdds = NA_real_, + enrichmentLogOddsNoShrinkage = NA_real_, + enrichmentSeNoShrinkage = NA_real_, + intercept = NA_real_, + interceptSe = NA_real_, + colocP1 = NA_real_, + colocP2 = NA_real_, + colocP12 = NA_real_, + effectiveMiRounds = NA_real_ + ) } @@ -520,7 +651,15 @@ qtlEnrichmentPipeline <- function( #' imputation sampler; each imputation round derives its own seed from it, so #' a fixed \code{seed} gives reproducible results. \code{NULL} (default) draws #' a nondeterministic seed. -#' @return A list of enrichment parameter estimates +#' @return A named list of enrichment parameter estimates, carrying the +#' fields the C++ estimator writes -- \code{Intercept}, +#' \code{sd (intercept)}, \code{Enrichment (no shrinkage)}, +#' \code{Enrichment (w/ shrinkage)}, \code{sd (no shrinkage)}, +#' \code{sd (w/ shrinkage)}, \code{Alternative (coloc) p1} / \code{p2} / +#' \code{p12} and \code{Effective MI rounds} -- plus +#' \code{unused_xqtl_variants}, the QTL variants of each region that no GWAS +#' variant matched. The names are upstream fastenloc's; +#' \code{\link{qtlEnrichmentPipeline}} is what renames them to a tidy table. #' #' @examples #' @@ -582,7 +721,7 @@ qtlEnrichment <- function( unmatchedVariants <- map(aligned, "unmatched_variants") susieQtlRegions <- map(aligned, .enrStripUnmatched) # cpp11 requires exact integer types for int parameters. - en <- list(qtlEnrichmentRcpp( + en <- qtlEnrichmentRcpp( rGwasPip = gwasPip, rQtlSusieFit = susieQtlRegions, piGwas = piGwas, @@ -593,7 +732,7 @@ qtlEnrichment <- function( besselCorrection = besselCorrection, numThreads = as.integer(numThreads), seed = if (is.null(seed)) NULL else as.integer(seed) - )) + ) en$unused_xqtl_variants <- unmatchedVariants en } diff --git a/R/tupleSelectors.R b/R/tupleSelectors.R index 436699cc..1bb8194c 100644 --- a/R/tupleSelectors.R +++ b/R/tupleSelectors.R @@ -33,6 +33,51 @@ mcols(x)[[k]] } +# ---- flavour-agnostic identity reads ---------------------------------------- +# `context` and `trait` are QTL-only axes: a GwasFineMappingResult is keyed by +# (study, method) alone. Pipelines that pair the two flavours -- colocPipeline, +# qtlEnrichmentPipeline -- need "this side has no context" as a VALUE, so these +# report an absent axis as NA instead of erroring on a missing column. + +# One identity column of a collection, as character, NA-filled when the +# collection has no such axis (length always matches the collection). +# @noRd +.fmrIdentityColumn <- function(x, column) { + values <- .tupleColumn(x, column) + if (is.null(values)) { + return(rep(NA_character_, nrow(x))) + } + as.character(values) +} + +# The same, for a single row. Indexed before coercion so reading one row does +# not convert the whole column. +# @noRd +.fmrIdentityAt <- function(x, column, i) { + values <- .tupleColumn(x, column) + if (is.null(values)) { + return(NA_character_) + } + as.character(values[[i]]) +} + +# Which flavour of fine-mapping collection a side is, for messages. +# @noRd +.fmrSideName <- function(x) { + if (methods::is(x, "QtlFineMappingResult")) "QTL" else "GWAS" +} + +# A human-readable identity for warnings: the side's flavour plus the identity +# fields it actually carries (the axes a GWAS side does not have are dropped +# rather than reported as NA). +# @noRd +.fmrTupleLabel <- function(side, ident, block = NULL) { + fields <- compact(c(ident, list(block = block))) + fields <- fields[!map_lgl(fields, is.na)] + shown <- str_c(names(fields), "='", unlist(fields), "'") + glue("{side} ({str_flatten(shown, ', ')})") +} + # One collection ELEMENT, whichever shape `x` has. On a RangedTupleList the # elements are the container itself; the DFrame-backed collections still keep # them in an `entry` column. Both shapes coexist until every collection has diff --git a/man/ColocResult.Rd b/man/ColocResult.Rd index 40f92d88..11b78314 100644 --- a/man/ColocResult.Rd +++ b/man/ColocResult.Rd @@ -8,10 +8,12 @@ ColocResult(pairs, variants, ldSketch = NULL) } \arguments{ \item{pairs}{A data frame with one row per tested pair, carrying at least -the identity columns (\code{study}, \code{context}, \code{trait}, -\code{method}, \code{gwasStudy}, \code{gwasMethod}), \code{blockId}, -\code{qtlCs}, \code{gwasCs}, \code{nSnps} and \code{PP.H0.abf} through -\code{PP.H4.abf}.} +the identity columns of both sides (\code{study}, \code{context}, +\code{trait}, \code{method} and \code{gwasStudy}, \code{gwasContext}, +\code{gwasTrait}, \code{gwasMethod}), \code{blockId}, \code{qtlCs}, +\code{gwasCs}, \code{nSnps} and \code{PP.H0.abf} through +\code{PP.H4.abf}. A GWAS side has no context or trait axis, so those two +columns are \code{NA} for it.} \item{variants}{A list, parallel to \code{pairs}' rows, of per-pair data frames with a \code{variant_id} column and a \code{SNP.PP.H4} column.} @@ -31,7 +33,9 @@ Assemble a \code{\linkS4class{ColocResult}} from a pair-level \examples{ pairs <- data.frame( study = "s1", context = "c1", trait = "g1", method = "susie", - gwasStudy = "G1", gwasMethod = "susie", blockId = "chr1_1_1000", + gwasStudy = "G1", gwasContext = NA_character_, + gwasTrait = NA_character_, gwasMethod = "susie", + blockId = "chr1_1_1000", qtlCs = 1L, gwasCs = 1L, nSnps = 2L, PP.H0.abf = 0.1, PP.H1.abf = 0.1, PP.H2.abf = 0.1, PP.H3.abf = 0.1, PP.H4.abf = 0.6 diff --git a/man/as.data.frame-ColocResult-method.Rd b/man/as.data.frame-ColocResult-method.Rd index 76c6dd39..86f8242a 100644 --- a/man/as.data.frame-ColocResult-method.Rd +++ b/man/as.data.frame-ColocResult-method.Rd @@ -25,7 +25,9 @@ Returns the pair-level view, which is the flat table \examples{ pairs <- data.frame( study = "s1", context = "c1", trait = "g1", method = "susie", - gwasStudy = "G1", gwasMethod = "susie", blockId = "chr1_1_1000", + gwasStudy = "G1", gwasContext = NA_character_, + gwasTrait = NA_character_, gwasMethod = "susie", + blockId = "chr1_1_1000", qtlCs = 1L, gwasCs = 1L, nSnps = 2L, PP.H0.abf = 0.1, PP.H1.abf = 0.1, PP.H2.abf = 0.1, PP.H3.abf = 0.1, PP.H4.abf = 0.6 diff --git a/man/colocPipeline.Rd b/man/colocPipeline.Rd index e2753f83..b1a5e1b1 100644 --- a/man/colocPipeline.Rd +++ b/man/colocPipeline.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/colocPipeline.R \name{colocPipeline} \alias{colocPipeline} -\title{Colocalization Pipeline (coloc.bf_bf over QTL + GWAS LBF matrices)} +\title{Colocalization Pipeline (coloc.bf_bf over paired LBF matrices)} \usage{ colocPipeline( qtlFineMappingResult, @@ -24,10 +24,14 @@ colocPipeline( ) } \arguments{ -\item{qtlFineMappingResult}{A \code{\link{QtlFineMappingResult}} (required).} +\item{qtlFineMappingResult}{The first side: a +\code{\link{QtlFineMappingResult}} or a +\code{\link{GwasFineMappingResult}} (required).} -\item{gwasInput}{Either a \code{\link{GwasSumStats}} or a -\code{\link{GwasFineMappingResult}}.} +\item{gwasInput}{The second side: a \code{\link{QtlFineMappingResult}}, a +\code{\link{GwasFineMappingResult}}, or the summary statistics to +fine-map inline (\code{\link{QtlSumStats}} or +\code{\link{GwasSumStats}}).} \item{filterLbfCs}{Logical. When \code{TRUE} (and \code{filterLbfCsSecondary} is \code{NULL}), keep only effects that produced a credible set @@ -58,17 +62,19 @@ consulted when \code{filterLbfCsSecondary} is non-NULL. Default \code{0.5}.} \code{5e-6}.} \item{finemappingMethods}{Character vector forwarded to -\code{\link{fineMappingPipeline}} when \code{gwasInput} is a -\code{GwasSumStats}. Default \code{"susie"}.} - -\item{returnGwasFineMapping}{Logical. When \code{TRUE}, attach the computed -\code{GwasFineMappingResult} on the returned data frame as attribute -\code{"gwasFineMapping"}. Default \code{FALSE}.} - -\item{enrichment}{Optional data.frame of per-(gwasStudy, qtlStudy, -qtlContext) enrichment factors with columns \code{gwasStudy}, -\code{qtlStudy}, \code{qtlContext}, \code{enrichment}. Output of -\code{\link{qtlEnrichmentPipeline}}. When non-\code{NULL}, each pair's +\code{\link{fineMappingPipeline}} when \code{gwasInput} is summary +statistics rather than a fine-mapping result. Default \code{"susie"}.} + +\item{returnGwasFineMapping}{Logical. When \code{TRUE}, attach the +fine-mapping result computed from \code{gwasInput} on the returned object +as attribute \code{"gwasFineMapping"}. Default \code{FALSE}.} + +\item{enrichment}{Optional data.frame of per-pair enrichment factors with +columns \code{gwasStudy}, \code{qtlStudy}, \code{qtlContext}, +\code{enrichment}, and optionally \code{gwasContext} / \code{gwasTrait} +(which the join uses when present, and which +\code{\link{qtlEnrichmentPipeline}} emits for a QTL outcome side). Output +of \code{\link{qtlEnrichmentPipeline}}. When non-\code{NULL}, each pair's \code{p12} prior is scaled to \code{min(p12 * (1 + enrichment), p12Max)} (the enrichment-informed colocalization variant, "enloc"). Pairs without a matching enrichment row fall back to the baseline \code{p12} with a @@ -103,11 +109,12 @@ exact alleles only, so a ref/alt swap is treated as a distinct variant.} \item{...}{Additional arguments forwarded to \code{coloc::coloc.bf_bf}.} } \value{ -A \code{\linkS4class{ColocResult}}: one element per tested (QTL - credible set, GWAS credible set, block) pair, holding that pair's aligned - variants with their \code{SNP.PP.H4}. Pair-level metadata carries the - identity columns (\code{study}, \code{context}, \code{trait}, - \code{method}, \code{gwasStudy}, \code{gwasMethod}), the block and stable +A \code{\linkS4class{ColocResult}}: one element per tested + (first-side credible set, second-side credible set, block) pair, holding + that pair's aligned variants with their \code{SNP.PP.H4}. Pair-level + metadata carries the identity columns (\code{study}, \code{context}, + \code{trait}, \code{method}, \code{gwasStudy}, \code{gwasContext}, + \code{gwasTrait}, \code{gwasMethod}), the block and stable credible-set ids (\code{blockId}, \code{qtlCs}, \code{gwasCs}), the standard coloc fields (\code{idx1}, \code{idx2}, \code{nSnps}, \code{hit1}, \code{hit2}, \code{PP.H0.abf} \ldots \code{PP.H4.abf}) and @@ -123,12 +130,25 @@ A \code{\linkS4class{ColocResult}}: one element per tested (QTL \code{\link{getColocGenes}}. } \description{ -Per-region pipeline that pairs a QTL - \code{\link{QtlFineMappingResult}} with a GWAS fine-mapping result (either - supplied directly as a \code{\link{GwasFineMappingResult}} or computed - inline from a \code{\link{GwasSumStats}}) and runs - \code{coloc::coloc.bf_bf} per (QTL tuple, GWAS tuple) pair to produce - per-pair colocalization posterior probabilities PP.H0-PP.H4. +Per-region pipeline that pairs two fine-mapping result + collections and runs \code{coloc::coloc.bf_bf} per (first-side tuple, + second-side tuple) pair to produce per-pair colocalization posterior + probabilities PP.H0-PP.H4. + + Either side may be a \code{\link{QtlFineMappingResult}} or a + \code{\link{GwasFineMappingResult}}: QTL-GWAS, QTL-QTL (two molecular + phenotypes) and GWAS-GWAS (two diseases) all run through one code path, + since nothing below the identity tuple depends on which flavour a side + is. The second side may also be handed in as summary statistics + (\code{\link{QtlSumStats}} or \code{\link{GwasSumStats}}), which are + fine-mapped inline. + + The argument names keep the QTL / GWAS wording of the common case. What + they mean generally is: \code{qtlFineMappingResult} is the side whose + identity is reported in the unprefixed \code{study} / \code{context} / + \code{trait} / \code{method} columns, and \code{gwasInput} is the side + reported in the \code{gwas}-prefixed ones. A GWAS side has no context or + trait axis, so those two columns are \code{NA} for it. } \section{Why \code{coloc.bf_bf} and not \code{coloc.susie}}{ @@ -145,30 +165,31 @@ directly. That choice carries three behaviours that keeps effects at a secondary coverage; otherwise the default filter drops effects whose prior variance is below \code{priorTol}. - \item \strong{Multiple-GWAS batching}: when several GWAS - fine-mapping rows fall in the same region they are merged into - one combined LBF matrix per QTL pair (one \code{coloc.bf_bf} - call covers them all). + \item \strong{Per-tuple LBF reuse}: each second-side tuple's LBF + matrix is extracted once and scored against every first-side + tuple, so the filtering above is applied once per tuple rather + than once per pair. } This pipeline preserves all three. - GWAS input dispatch: + Second-side input dispatch: \itemize{ - \item \code{gwasInput} is a \code{\link{GwasSumStats}}: GWAS - fine-mapping is performed inline by + \item \code{gwasInput} is a \code{\link{QtlSumStats}} or a + \code{\link{GwasSumStats}}: it is fine-mapped inline by \code{\link{fineMappingPipeline}} with the supplied \code{finemappingMethods} (default \code{"susie"}). - \item \code{gwasInput} is a \code{\link{GwasFineMappingResult}}: - used directly; no inline fine-mapping. + \item \code{gwasInput} is a fine-mapping result: used directly; no + inline fine-mapping. } } \section{LD-sketch identity check}{ If \code{getLdSketch(qtlFineMappingResult)} is non-\code{NULL}, it must match - the LD sketch on \code{gwasInput}. Mismatch is a hard error. When the QTL - FMR's \code{ldSketch} is \code{NULL} (individual-level fit), the validation - is skipped on the QTL side. + the LD sketch on \code{gwasInput}. Mismatch is a hard error. When the + first side's \code{ldSketch} is \code{NULL} (individual-level fit), the + validation is skipped on that side and the second side's panel is what the + result carries forward. } \examples{ @@ -176,4 +197,10 @@ data(qtlFineMappingLbfExample) data(gwasFineMappingLbfExample) colocPipeline(qtlFineMappingLbfExample, gwasInput = gwasFineMappingLbfExample) +# Either side may be a QTL result. Pairing this collection against itself +# colocalizes its two contexts, and reports the second side's context and +# trait in gwasContext / gwasTrait. +res <- colocPipeline(qtlFineMappingLbfExample, + gwasInput = qtlFineMappingLbfExample) +unique(getColocPairs(res)[, c("context", "gwasContext")]) } diff --git a/man/qtlEnrichment.Rd b/man/qtlEnrichment.Rd index e1216d0c..8acc477c 100644 --- a/man/qtlEnrichment.Rd +++ b/man/qtlEnrichment.Rd @@ -66,7 +66,15 @@ a fixed \code{seed} gives reproducible results. \code{NULL} (default) draws a nondeterministic seed.} } \value{ -A list of enrichment parameter estimates +A named list of enrichment parameter estimates, carrying the + fields the C++ estimator writes -- \code{Intercept}, + \code{sd (intercept)}, \code{Enrichment (no shrinkage)}, + \code{Enrichment (w/ shrinkage)}, \code{sd (no shrinkage)}, + \code{sd (w/ shrinkage)}, \code{Alternative (coloc) p1} / \code{p2} / + \code{p12} and \code{Effective MI rounds} -- plus + \code{unused_xqtl_variants}, the QTL variants of each region that no GWAS + variant matched. The names are upstream fastenloc's; + \code{\link{qtlEnrichmentPipeline}} is what renames them to a tidy table. } \description{ Largely follows from fastenloc diff --git a/man/qtlEnrichmentPipeline.Rd b/man/qtlEnrichmentPipeline.Rd index cd5558fa..e132d638 100644 --- a/man/qtlEnrichmentPipeline.Rd +++ b/man/qtlEnrichmentPipeline.Rd @@ -17,9 +17,9 @@ qtlEnrichmentPipeline( ) } \arguments{ -\item{gwasFineMappingResult}{See above.} +\item{gwasFineMappingResult}{The outcome side; see above.} -\item{qtlFineMappingResult}{See above.} +\item{qtlFineMappingResult}{The annotation side; see above.} \item{numGwas}{Number of GWAS variants used to estimate \code{piGwas}. When \code{NULL} (default) it is estimated from the data -- bias warning applies @@ -44,41 +44,67 @@ if the input PIP vector is not genome-wide.} \item{...}{Additional arguments forwarded to \code{\link{qtlEnrichment}}.} } \value{ -A tibble with one row per (gwasStudy, qtlStudy, qtlContext) - triple and columns \code{gwasStudy}, \code{qtlStudy}, \code{qtlContext}, - \code{enrichment}, \code{enrichmentSe}, \code{enrichmentLogOdds}, plus any - extras the underlying estimator emits. Suitable as the \code{enrichment} - argument to \code{\link{colocPipeline}} (which joins on the same triple). +A tibble with one row per (outcome trait, annotation unit) pair. + The identity columns are \code{gwasStudy}, \code{gwasContext}, + \code{gwasTrait}, \code{qtlStudy}, \code{qtlContext}; the axes a side does + not have are \code{NA} (\code{gwasContext} / \code{gwasTrait} for a GWAS + outcome, \code{qtlContext} for a GWAS annotation). Suitable as the + \code{enrichment} argument to \code{\link{colocPipeline}}, which joins on + those columns. + + The estimates are \code{enrichmentLogOdds}, the enrichment parameter + \eqn{a_1} on the log-odds scale, with its standard error + \code{enrichmentSe}; \code{enrichment} is the same quantity as a + multiplicative factor, \eqn{e^{a_1} - 1}, which is what + \code{colocPipeline} scales \code{p12} by (so \eqn{a_1 = 0} leaves the + prior untouched). \code{enrichmentLogOddsNoShrinkage} and + \code{enrichmentSeNoShrinkage} are the same estimate before shrinkage, + \code{intercept} / \code{interceptSe} are \eqn{a_0}, and \code{colocP1}, + \code{colocP2}, \code{colocP12} are the enrichment-informed coloc priors + the estimator derives from \eqn{(a_0, a_1)} -- an alternative to scaling + a baseline \code{p12}. \code{effectiveMiRounds} is how many + multiple-imputation rounds survived outlier filtering. } \description{ -Genome-wide pipeline that computes per-pair (GWAS study, QTL - context) enrichment estimates by passing the GWAS PIP vector and the QTL +Genome-wide pipeline that computes per-pair enrichment + estimates by passing an outcome PIP vector and a set of annotation credible-set posteriors to \code{\link{qtlEnrichment}}. The returned table feeds \code{\link{colocPipeline}} via its \code{enrichment} argument. + Either side may be a \code{\link{QtlFineMappingResult}} or a + \code{\link{GwasFineMappingResult}}, so QTL-in-GWAS, QTL-in-QTL and + GWAS-in-GWAS enrichment all run through one code path. The argument names + keep the QTL / GWAS wording of the common case; what they mean generally + is that \code{gwasFineMappingResult} is the \strong{outcome} whose PIPs + are scanned and \code{qtlFineMappingResult} is the \strong{annotation} + whose region fits are tested for enrichment within them. + \strong{Not gene-parallelisable}: the enrichment estimator runs over the full -genome of GWAS PIPs and the full collection of QTL fits at once. +genome of outcome PIPs and the full collection of annotation fits at once. } \section{Inputs}{ \itemize{ - \item \code{gwasFineMappingResult}: a genome-wide - \code{\link{GwasFineMappingResult}} (one row per (study, LD - block) tuple). Each entry's \code{FineMappingRow$trimmedFit} - must carry a \code{pip} vector. - \item \code{qtlFineMappingResult}: the genome-wide - \code{\link{QtlFineMappingResult}}. Each entry's + \item \code{gwasFineMappingResult}: the genome-wide outcome collection. + Each entry's \code{FineMappingRow$trimmedFit} must carry a \code{pip} + vector. One PIP vector is built per outcome trait -- keyed by + \code{study} for a GWAS collection (one trait per study) and by + (\code{study}, \code{context}, \code{trait}) for a QTL one, since a + variant's PIP differs between molecular traits and pooling them would + collide. + \item \code{qtlFineMappingResult}: the annotation collection. Each entry's \code{trimmedFit} must carry \code{alpha}, \code{pip}, and - prior-variance fields (\code{V}). + prior-variance fields (\code{V}). Its region fits are pooled per + (\code{study}, \code{context}) -- per \code{study} alone for a GWAS + collection, which has no context axis. } } \section{LD-sketch identity check}{ - The GWAS \code{FineMappingResultBase} must - have a non-NULL \code{ldSketch} (RSS-derived). If the QTL FMR also has a - non-NULL \code{ldSketch}, the two must match exactly. When the QTL FMR's - \code{ldSketch} is NULL (individual-level QTL fit), validation is skipped - on the QTL side. + A GWAS outcome collection must have a + non-NULL \code{ldSketch} (it should be RSS-derived). Where both sides + carry one, the two must match exactly; a \code{NULL} on either side (an + individual-level fit) skips the check. } \examples{ diff --git a/src/qtl_enrichment.h b/src/qtl_enrichment.h index c9735e59..88b4e8cb 100644 --- a/src/qtl_enrichment.h +++ b/src/qtl_enrichment.h @@ -10,6 +10,7 @@ #include #include #include +#include using namespace cpp11; using namespace arma; @@ -116,6 +117,17 @@ std::vector filter_outlier_indices( } } + // A filter that removes EVERY round is not a filter: it means the + // dispersion estimate was degenerate rather than that every round was an + // outlier. With a single MI round and the Bessel correction the sd is + // 0/0, so no |x - mean| <= 3 * NaN comparison holds and every round is + // dropped -- which then divides by zero throughout the MI combination and + // reports NaN for every estimate. Keep them all instead. + if(kept.empty()) { + kept.resize(n); + std::iota(kept.begin(), kept.end(), 0); + } + return kept; } @@ -313,8 +325,16 @@ std::map qtl_enrichment_workhorse( bv0 += pow(a0_vec[k] - a0_est, 2.0); bv1 += pow(a1_vec[k] - a1_est, 2.0); } - bv0 /= (m - 1); - bv1 /= (m - 1); + // Rubin's rules need two imputations to estimate a between-round variance; + // with one surviving round only the within-round variance is measurable, + // and dividing by (m - 1) would make every downstream estimate NaN. + if (m > 1) { + bv0 /= (m - 1); + bv1 /= (m - 1); + } else { + bv0 = 0; + bv1 = 0; + } var0 /= m; var1 /= m; @@ -348,7 +368,7 @@ std::map qtl_enrichment_workhorse( for (size_t i = 0; i < m; i++) { bv1_shrink += pow(a1_shrink_vec[i] - a1_shrink_est, 2.0); } - bv1_shrink /= (m - 1); + bv1_shrink = (m > 1) ? bv1_shrink / (m - 1) : 0.0; double sd1_shrink = sqrt(var1_shrink + bv1_shrink * (m + 1.0) / m); // 2. Shrink the combined estimate again diff --git a/tests/testthat/test_ColocResult.R b/tests/testthat/test_ColocResult.R index 215021fd..6c1127f9 100644 --- a/tests/testthat/test_ColocResult.R +++ b/tests/testthat/test_ColocResult.R @@ -9,14 +9,16 @@ gwasCs = 1L, blockId = "chr1_1_1000", pp4 = 0.6, - trait = "g1" + trait = "g1", + gwasTrait = NA_character_ ) { n <- max( length(qtlCs), length(gwasCs), length(blockId), length(pp4), - length(trait) + length(trait), + length(gwasTrait) ) data.frame( study = "s1", @@ -24,6 +26,8 @@ trait = trait, method = "susie", gwasStudy = "G1", + gwasContext = NA_character_, + gwasTrait = gwasTrait, gwasMethod = "susie", blockId = blockId, qtlCs = as.integer(qtlCs), @@ -139,6 +143,17 @@ test_that("getColocGenes: keeps distinct genes apart", { expect_setequal(genes$trait, c("g1", "g2")) }) +test_that("getColocGenes: keeps distinct second-side traits apart", { + # A QTL-QTL colocalization has a trait axis on BOTH sides. Grouping on the + # second side's study alone would pool two molecular phenotypes into one + # unit and noisy-OR their posteriors together. + pairs <- .cr_pairs(gwasTrait = c("h1", "h2"), pp4 = c(0.5, 0.25)) + x <- ColocResult(pairs, .cr_variants(2L)) + genes <- getColocGenes(x) + expect_equal(nrow(genes), 2L) + expect_setequal(genes$gwasTrait, c("h1", "h2")) +}) + test_that("getColocVariants(pooled): pools by the same rule as genes", { pairs <- .cr_pairs( qtlCs = c(1L, 2L), @@ -244,8 +259,8 @@ test_that(".crFirstAtLeast: absorbs floating-point error, not real shortfall", { test_that("show summarizes the pairs, studies and best PP.H4", { x <- ColocResult(.cr_pairs(), .cr_variants()) expect_output(show(x), "ColocResult with 1 colocalized pair\\(s\\)") - expect_output(show(x), "QTL studies") - expect_output(show(x), "GWAS studies") + expect_output(show(x), "studies") + expect_output(show(x), "paired with") expect_output(show(x), "variants") expect_output(show(x), "max PP.H4") }) diff --git a/tests/testthat/test_colocPipeline.R b/tests/testthat/test_colocPipeline.R index 5a02eff4..d0d42ce4 100644 --- a/tests/testthat/test_colocPipeline.R +++ b/tests/testthat/test_colocPipeline.R @@ -137,6 +137,38 @@ context("colocPipeline") ) } +.cp_makeQtlSumstats <- function(study = "Q1", qc = TRUE) { + gr <- GenomicRanges::GRanges( + seqnames = "chr1", + ranges = IRanges::IRanges( + start = seq(100L, by = 100L, length.out = 5L), + width = 1L + ) + ) + S4Vectors::mcols(gr) <- S4Vectors::DataFrame( + SNP = sprintf("chr1:%d:A:G", 100L * (1:5)), + A1 = rep("A", 5), + A2 = rep("G", 5), + Z = rnorm(5), + N = rep(1000L, 5) + ) + QtlSumStats( + study = study, + context = "c1", + trait = "t1", + entry = list(gr), + genome = "hg19", + ldSketch = .cp_makeHandle(), + qcInfo = if (qc) list(step1 = "ok") else list() + ) +} + +# The enrichment lookup takes each side's identity list, as the scoring loop +# builds it. +.cp_side <- function(study, context = NA_character_, trait = NA_character_) { + list(study = study, context = context, trait = trait) +} + .cp_mockColocBfBf <- function() { function(qLbf, gLbf, p1, p2, p12, ...) { list( @@ -159,20 +191,20 @@ context("colocPipeline") # Input-type validation # =========================================================================== -test_that("colocPipeline: rejects non-QtlFineMappingResult qtlFmr", { +test_that("colocPipeline: rejects a non-fine-mapping first side", { expect_error( colocPipeline( qtlFineMappingResult = "no", gwasInput = .cp_makeGwasFmr() ), - "must be a QtlFineMappingResult" + "must be a QtlFineMappingResult or a GwasFineMappingResult" ) }) -test_that("colocPipeline: rejects gwasInput that is neither GwasSumStats nor GwasFineMappingResult", { +test_that("colocPipeline: rejects an unusable second side", { expect_error( colocPipeline(qtlFineMappingResult = .cp_makeQtlFmr(), gwasInput = 42L), - "must be a GwasSumStats or a GwasFineMappingResult" + "must be a fine-mapping result" ) }) @@ -247,6 +279,94 @@ test_that("colocPipeline: returns one row per (QTL tuple, GWAS tuple) pair", { expect_setequal(out$gwasStudy, c("G1", "G2")) }) +# =========================================================================== +# Either side may be a QTL or a GWAS fine-mapping result +# =========================================================================== + +test_that("colocPipeline: pairs two QTL fine-mapping results", { + qfmr <- .cp_makeQtlFmr() + # The two second-side rows share (study, method, block) and differ only on + # trait -- keyed on the GWAS 2-tuple, one would silently replace the other. + other <- .cp_makeQtlFmr( + tuples = list( + c("Q2", "c1", "t1", "susie"), + c("Q2", "c1", "t2", "susie") + ) + ) + local_mocked_bindings(coloc.bf_bf = .cp_mockColocBfBf(), .package = "coloc") + out <- suppressWarnings(colocPipeline( + qtlFineMappingResult = qfmr, + gwasInput = other + )) + expect_equal(nrow(out), 2L) + expect_setequal(out$gwasStudy, "Q2") + expect_setequal(out$gwasContext, "c1") + expect_setequal(out$gwasTrait, c("t1", "t2")) +}) + +test_that("colocPipeline: pairs two GWAS fine-mapping results", { + first <- .cp_makeGwasFmr(tuples = list(c("G1", "susie"))) + second <- .cp_makeGwasFmr(tuples = list(c("G2", "susie"))) + local_mocked_bindings(coloc.bf_bf = .cp_mockColocBfBf(), .package = "coloc") + out <- suppressWarnings(colocPipeline( + qtlFineMappingResult = first, + gwasInput = second + )) + expect_equal(nrow(out), 1L) + expect_equal(as.character(out$study), "G1") + expect_equal(as.character(out$gwasStudy), "G2") + # Neither side has a context or trait axis, so all four are reported as NA + # rather than invented. + expect_true(all(is.na(c( + out$context, + out$trait, + out$gwasContext, + out$gwasTrait + )))) +}) + +test_that("colocPipeline: a GWAS first side names itself in warnings", { + first <- .cp_makeGwasFmr() + second <- .cp_makeGwasFmr(tuples = list(c("G2", "susie"))) + local_mocked_bindings( + coloc.bf_bf = function(...) stop("synthetic test failure"), + .package = "coloc" + ) + expect_warning( + colocPipeline(qtlFineMappingResult = first, gwasInput = second), + "GWAS \\(study='G1', method='susie'\\)" + ) +}) + +test_that("colocPipeline: resolves QtlSumStats by fine-mapping it", { + first <- .cp_makeGwasFmr() + qss <- .cp_makeQtlSumstats() + resolved <- .cp_makeQtlFmr() + local_mocked_bindings(coloc.bf_bf = .cp_mockColocBfBf(), .package = "coloc") + local_mocked_bindings( + fineMappingPipeline = function(data, methods, ...) resolved, + .package = "pecotmr" + ) + out <- suppressWarnings(colocPipeline( + qtlFineMappingResult = first, + gwasInput = qss, + returnGwasFineMapping = TRUE, + adjustPips = FALSE + )) + expect_equal(as.character(out$gwasTrait), "t1") + expect_identical(attr(out, "gwasFineMapping"), resolved) +}) + +test_that("colocPipeline: rejects an un-QCd QtlSumStats second side", { + expect_error( + colocPipeline( + qtlFineMappingResult = .cp_makeQtlFmr(), + gwasInput = .cp_makeQtlSumstats(qc = FALSE) + ), + "has no QC record" + ) +}) + test_that("colocPipeline: resolves GwasSumStats via fineMappingPipeline (mocked)", { qfmr <- .cp_makeQtlFmr() gss <- .cp_makeGwasSumstats() @@ -351,6 +471,8 @@ test_that("colocPipeline: empty result has the documented schema", { "trait", "method", "gwasStudy", + "gwasContext", + "gwasTrait", "gwasMethod", "blockId", "qtlCs", @@ -504,7 +626,14 @@ test_that(".colocLookupEnrichment: returns the value for a (gwasStudy, qtlStudy, enrichment = c(2.0, 3.5), stringsAsFactors = FALSE ) - expect_equal(pecotmr:::.colocLookupEnrichment(enr, "G2", "Q1", "c1"), 3.5) + expect_equal( + pecotmr:::.colocLookupEnrichment( + enr, + .cp_side("G2"), + .cp_side("Q1", "c1") + ), + 3.5 + ) }) test_that(".colocLookupEnrichment: returns NA when no row matches", { @@ -517,16 +646,14 @@ test_that(".colocLookupEnrichment: returns NA when no row matches", { ) expect_true(is.na(pecotmr:::.colocLookupEnrichment( enr, - "ghost", - "Q1", - "c1" + .cp_side("ghost"), + .cp_side("Q1", "c1") ))) # qtlStudy mismatch also a miss. expect_true(is.na(pecotmr:::.colocLookupEnrichment( enr, - "G1", - "Qghost", - "c1" + .cp_side("G1"), + .cp_side("Qghost", "c1") ))) }) @@ -615,6 +742,56 @@ test_that("colocPipeline: enrichment hit scales p12 and emits enrichment/p12Used expect_equal(unique(out$p12Used), min(5e-6 * 3, 1e-3)) }) +test_that("colocPipeline: enrichment joins on the second side's trait", { + # A QTL second side puts two traits under one study. Joining on the study + # alone would give both the first row's factor. + qfmr <- .cp_makeQtlFmr() + other <- .cp_makeQtlFmr( + tuples = list( + c("Q2", "c1", "t1", "susie"), + c("Q2", "c1", "t2", "susie") + ) + ) + enr <- data.frame( + gwasStudy = c("Q2", "Q2"), + gwasContext = c("c1", "c1"), + gwasTrait = c("t1", "t2"), + qtlStudy = c("Q1", "Q1"), + qtlContext = c("c1", "c1"), + enrichment = c(1.0, 3.0), + stringsAsFactors = FALSE + ) + local_mocked_bindings(coloc.bf_bf = .cp_mockColocBfBf(), .package = "coloc") + out <- suppressWarnings(colocPipeline( + qtlFineMappingResult = qfmr, + gwasInput = other, + enrichment = enr + )) + pairs <- getColocPairs(out) + expect_equal( + pairs$enrichment[order(pairs$gwasTrait)], + c(1.0, 3.0) + ) +}) + +test_that("colocPipeline: rejects an enrichment table with repeated keys", { + enr <- data.frame( + gwasStudy = c("G1", "G1"), + qtlStudy = c("Q1", "Q1"), + qtlContext = c("c1", "c1"), + enrichment = c(2.0, 3.0), + stringsAsFactors = FALSE + ) + expect_error( + colocPipeline( + qtlFineMappingResult = .cp_makeQtlFmr(), + gwasInput = .cp_makeGwasFmr(), + enrichment = enr + ), + "repeated" + ) +}) + test_that("colocPipeline: enrichment miss warns and falls back to baseline p12", { qfmr <- .cp_makeQtlFmr() gfmr <- .cp_makeGwasFmr() diff --git a/tests/testthat/test_qtlEnrichmentPipeline.R b/tests/testthat/test_qtlEnrichmentPipeline.R index 87cbf317..5e9e7484 100644 --- a/tests/testthat/test_qtlEnrichmentPipeline.R +++ b/tests/testthat/test_qtlEnrichmentPipeline.R @@ -5,6 +5,20 @@ context("qtlEnrichmentPipeline") # small fixture, but the heavy mixture-of-enrichment estimator never fires. # =========================================================================== +# The internal builders take an identity list keyed by the collection's own +# columns; the axes a GWAS collection lacks are NA. +.qep_gwasIdent <- function( + study, + context = NA_character_, + trait = NA_character_ +) { + list(study = study, context = context, trait = trait) +} + +.qep_qtlIdent <- function(study, context) { + list(study = study, context = context) +} + .qep_makeHandle <- function( snp_n = 6L, n_samples = 30L, @@ -101,16 +115,162 @@ context("qtlEnrichmentPipeline") } # Mock that returns a plausible enrichment list. +# The shape qtlEnrichment() really returns: the C++ estimator's own field +# names (src/qtl_enrichment.h), not the pipeline's output columns. A mock that +# invented the output names is what let a total field-name mismatch -- +# every estimate NA on real data -- sit undetected. +# +# `value` is the enrichment the pipeline should REPORT, so it is encoded here +# as the log-odds the estimator would have produced for it. .qep_mockEnrichment <- function(value = 1.5) { + logOdds <- log1p(value) function(gwasPip, susieQtlRegions, ...) { list( - enrichment = value, - enrichmentSe = 0.1, - enrichmentLogOdds = log(value) + "Intercept" = -6.5, + "sd (intercept)" = 0.5, + "Enrichment (no shrinkage)" = logOdds, + "Enrichment (w/ shrinkage)" = logOdds, + "sd (no shrinkage)" = 0.2, + "sd (w/ shrinkage)" = 0.1, + "Alternative (coloc) p1" = 1e-4, + "Alternative (coloc) p2" = 1e-4, + "Alternative (coloc) p12" = 5e-6, + "Effective MI rounds" = 25, + unused_xqtl_variants = list() ) } } +# =========================================================================== +# End-to-end against the real estimator (no mock) +# =========================================================================== + +test_that("qtlEnrichmentPipeline: the real estimator fills the value columns", { + # Unmocked on purpose. The C++ kernel's field names are the contract, and + # a mismatch between them and this pipeline's output columns is invisible + # to a mocked run -- which is how every estimate came back NA on real data + # while the mocked tests stayed green. + ids <- sprintf("chr1:%d:A:G", 100L * seq_len(50L)) + gwasPip <- rep(0.01, 50L) + gwasPip[c(5L, 20L, 35L)] <- c(0.8, 0.6, 0.9) + alpha <- matrix(0.001, nrow = 2L, ncol = 50L) + alpha[1L, 5L] <- 0.95 + alpha[2L, 20L] <- 0.95 + alpha <- alpha / rowSums(alpha) + sketch <- .qep_makeHandle() + gfmr <- GwasFineMappingResult( + study = "G1", + method = "susie", + entry = list(fineMappingRow( + variantIds = ids, + susieFit = list( + alpha = matrix(1 / 50, nrow = 1L, ncol = 50L), + pip = setNames(gwasPip, ids), + V = 0.5 + ), + topLoci = data.frame(variant_id = ids, pip = gwasPip) + )), + ldSketch = sketch + ) + qfmr <- QtlFineMappingResult( + study = "Q1", + context = "c1", + trait = "t1", + method = "susie", + entry = list(fineMappingRow( + variantIds = ids, + susieFit = list( + alpha = alpha, + pip = setNames(colSums(alpha), ids), + V = c(0.5, 0.3) + ), + topLoci = data.frame(variant_id = ids, pip = colSums(alpha)) + )), + ldSketch = sketch + ) + suppressWarnings(capture.output( + out <- qtlEnrichmentPipeline( + gwasFineMappingResult = gfmr, + qtlFineMappingResult = qfmr, + impN = 5, + seed = 1L + ) + )) + expect_equal(nrow(out), 1L) + valueCols <- names(pecotmr:::.enrNaEnrichment()) + expect_true(all(valueCols %in% colnames(out))) + expect_true(all(is.finite(unlist(out[, valueCols])))) + # The multiplicative factor and the log-odds must agree. + expect_equal(out$enrichment, expm1(out$enrichmentLogOdds)) +}) + +# =========================================================================== +# Either side may be a QTL or a GWAS fine-mapping result +# =========================================================================== + +test_that("qtlEnrichmentPipeline: keys a QTL outcome side by its trait", { + # Two molecular traits over the same variants with DIFFERENT pips. Keyed by + # study alone they collide -- the pipeline would abort on conflicting PIPs + # -- and one estimate would stand in for both traits. + outcome <- QtlFineMappingResult( + study = c("Q1", "Q1"), + context = c("c1", "c1"), + trait = c("t1", "t2"), + method = c("susie", "susie"), + entry = list( + .qep_makeFmEntry(pip = c(0.9, 0.05, 0.02, 0.02, 0.01)), + .qep_makeFmEntry(pip = c(0.1, 0.7, 0.1, 0.05, 0.05)) + ), + ldSketch = .qep_makeHandle() + ) + local_mocked_bindings( + qtlEnrichment = .qep_mockEnrichment(2.0), + .package = "pecotmr" + ) + out <- qtlEnrichmentPipeline( + gwasFineMappingResult = outcome, + qtlFineMappingResult = .qep_makeQtlFmr() + ) + expect_equal(nrow(out), 2L) + expect_setequal(out$gwasStudy, "Q1") + expect_setequal(out$gwasContext, "c1") + expect_setequal(out$gwasTrait, c("t1", "t2")) +}) + +test_that("qtlEnrichmentPipeline: pairs two GWAS collections", { + local_mocked_bindings( + qtlEnrichment = .qep_mockEnrichment(2.0), + .package = "pecotmr" + ) + out <- qtlEnrichmentPipeline( + gwasFineMappingResult = .qep_makeGwasFmr(studies = "G1"), + qtlFineMappingResult = .qep_makeGwasFmr(studies = "G2") + ) + expect_equal(nrow(out), 1L) + expect_equal(out$gwasStudy, "G1") + expect_equal(out$qtlStudy, "G2") + # Neither side has a context or trait axis, so all three read NA. + expect_true(all(is.na(c( + out$gwasContext, + out$gwasTrait, + out$qtlContext + )))) +}) + +test_that("qtlEnrichmentPipeline: a QTL outcome side may carry no ldSketch", { + # The RSS-derived requirement is a GWAS-side contract; an individual-level + # QTL outcome has no panel to check. + local_mocked_bindings( + qtlEnrichment = .qep_mockEnrichment(2.0), + .package = "pecotmr" + ) + out <- qtlEnrichmentPipeline( + gwasFineMappingResult = .qep_makeQtlFmr(with_sketch = FALSE), + qtlFineMappingResult = .qep_makeQtlFmr() + ) + expect_equal(nrow(out), 1L) +}) + # =========================================================================== # Input-type validation # =========================================================================== @@ -214,11 +374,7 @@ test_that("qtlEnrichmentPipeline: distinguishes two QTL studies that share a con local_mocked_bindings( qtlEnrichment = function(gwasPip, susieQtlRegions, ...) { capturedRegions[[length(capturedRegions) + 1L]] <<- susieQtlRegions - list( - enrichment = 2.0, - enrichmentSe = 0.1, - enrichmentLogOdds = log(2) - ) + .qep_mockEnrichment(2.0)() }, .package = "pecotmr" ) @@ -282,11 +438,11 @@ test_that("qtlEnrichmentPipeline: empty input collections yield the empty schema colnames(out), c( "gwasStudy", + "gwasContext", + "gwasTrait", "qtlStudy", "qtlContext", - "enrichment", - "enrichmentSe", - "enrichmentLogOdds" + names(pecotmr:::.enrNaEnrichment()) ) ) }) @@ -297,7 +453,7 @@ test_that("qtlEnrichmentPipeline: empty input collections yield the empty schema test_that(".enrBuildGwasPipVector: extracts pip per study", { gfmr <- .qep_makeGwasFmr() - out <- pecotmr:::.enrBuildGwasPipVector(gfmr, "G1") + out <- pecotmr:::.enrBuildGwasPipVector(gfmr, .qep_gwasIdent("G1")) expect_equal(length(out), 3L) expect_setequal(names(out), sprintf("chr1:%d:A:G", 100L * (1:3))) }) @@ -321,7 +477,7 @@ test_that(".enrBuildGwasPipVector: deduplicates identical PIPs across blocks", { entry = list(e1, e2), ldSketch = .qep_makeHandle() ) - out <- pecotmr:::.enrBuildGwasPipVector(g, "G1") + out <- pecotmr:::.enrBuildGwasPipVector(g, .qep_gwasIdent("G1")) expect_setequal( names(out), c("chr1:100:A:G", "chr1:200:A:G", "chr1:300:A:G") @@ -351,14 +507,14 @@ test_that(".enrBuildGwasPipVector: conflicting PIPs across blocks errors", { ldSketch = .qep_makeHandle() ) expect_error( - pecotmr:::.enrBuildGwasPipVector(g, "G1"), + pecotmr:::.enrBuildGwasPipVector(g, .qep_gwasIdent("G1")), "conflicting PIPs" ) }) test_that(".enrBuildQtlRegionsList: returns per-entry fit shapes for a (study, context) hit", { qfmr <- .qep_makeQtlFmr(contexts = c("c1", "c2")) - out <- pecotmr:::.enrBuildQtlRegionsList(qfmr, "Q1", "c1") + out <- pecotmr:::.enrBuildQtlRegionsList(qfmr, .qep_qtlIdent("Q1", "c1")) expect_equal(length(out), 1L) expect_true(!is.null(out[[1L]]$alpha)) expect_true(!is.null(out[[1L]]$pip)) @@ -368,12 +524,18 @@ test_that(".enrBuildQtlRegionsList: returns empty list when the (study, context) qfmr <- .qep_makeQtlFmr(contexts = c("c1", "c2")) # Correct context but wrong study -> no hit, even though context exists. expect_equal( - length(pecotmr:::.enrBuildQtlRegionsList(qfmr, "Q_ghost", "c1")), + length(pecotmr:::.enrBuildQtlRegionsList( + qfmr, + .qep_qtlIdent("Q_ghost", "c1") + )), 0L ) # Correct study but wrong context. expect_equal( - length(pecotmr:::.enrBuildQtlRegionsList(qfmr, "Q1", "c_ghost")), + length(pecotmr:::.enrBuildQtlRegionsList( + qfmr, + .qep_qtlIdent("Q1", "c_ghost") + )), 0L ) }) @@ -435,7 +597,9 @@ test_that("qtlEnrichment: real C++ kernel returns the expected keys (numGwas + p verbose = FALSE ) expect_type(res, "list") - en <- res[[1L]] + # Flat, and addressable by name: the estimates used to sit inside an + # unnamed first element, where every by-name read of them returned NULL. + en <- res expectedKeys <- c( "Intercept", "Enrichment (no shrinkage)", @@ -450,6 +614,30 @@ test_that("qtlEnrichment: real C++ kernel returns the expected keys (numGwas + p expect_true(all(is.finite(unlist(en[expectedKeys])))) }) +test_that("qtlEnrichment: a single MI round is estimable, not NaN", { + # With one round the Bessel-corrected dispersion is 0/0, so no + # |x - mean| <= 3 * NaN comparison held, the outlier filter dropped the + # only round, and dividing by (m - 1) = 0 made every estimate NaN. + fx <- .qep_makeRealKernelInputs() + res <- qtlEnrichment( + gwasPip = fx$gwasPip, + susieQtlRegions = fx$susieQtlRegions, + numGwas = 5000, + piQtl = 0.5, + impN = 1, + numThreads = 1, + verbose = FALSE, + seed = 1L + ) + expect_equal(res[["Effective MI rounds"]], 1) + expect_true(all(is.finite(c( + res[["Enrichment (w/ shrinkage)"]], + res[["sd (w/ shrinkage)"]], + res[["Intercept"]], + res[["Alternative (coloc) p12"]] + )))) +}) + test_that("qtlEnrichment: numGwas omitted -> estimates piGwas from data + warns", { fx <- .qep_makeRealKernelInputs(nSnps = 30, causalIdx = c(5, 15)) expect_warning( @@ -620,80 +808,38 @@ test_that("qtlEnrichment: alignNames=FALSE recomputes only the unmatched set", { # .enrFlattenEnrichment(): shape coercion + NA fallbacks # =========================================================================== -test_that(".enrFlattenEnrichment: numeric scalar falls back to all-NA", { - out <- pecotmr:::.enrFlattenEnrichment(1.5) - expect_equal( - out, - list( - enrichment = NA_real_, - enrichmentSe = NA_real_, - enrichmentLogOdds = NA_real_ - ) - ) -}) - -test_that(".enrFlattenEnrichment: plain list picks named scalar fields", { - out <- pecotmr:::.enrFlattenEnrichment( - list(enrichment = 2.0, enrichmentSe = 0.1, enrichmentLogOdds = log(2)) - ) - expect_equal(out$enrichment, 2.0) - expect_equal(out$enrichmentSe, 0.1) - expect_equal(out$enrichmentLogOdds, log(2)) -}) - -test_that(".enrFlattenEnrichment: non-empty matrix picks columns by name", { - m <- matrix(c(2.0, 0.1, log(2)), nrow = 1) - colnames(m) <- c("enrichment", "enrichmentSe", "enrichmentLogOdds") - out <- pecotmr:::.enrFlattenEnrichment(m) +test_that(".enrFlattenEnrichment: reads the estimator's own field names", { + out <- pecotmr:::.enrFlattenEnrichment(.qep_mockEnrichment(2.0)()) + # enrichment is expm1 of the log-odds, so colocPipeline's + # p12 * (1 + enrichment) is the enloc prior p12 * exp(a1). expect_equal(out$enrichment, 2.0) + expect_equal(out$enrichmentLogOdds, log(3)) expect_equal(out$enrichmentSe, 0.1) - expect_equal(out$enrichmentLogOdds, log(2)) -}) - -test_that(".enrFlattenEnrichment: empty data.frame falls back to all-NA", { - df <- data.frame( - enrichment = numeric(0), - enrichmentSe = numeric(0), - enrichmentLogOdds = numeric(0) - ) - out <- pecotmr:::.enrFlattenEnrichment(df) - expect_equal( - out, - list( - enrichment = NA_real_, - enrichmentSe = NA_real_, - enrichmentLogOdds = NA_real_ - ) - ) + expect_equal(out$enrichmentSeNoShrinkage, 0.2) + expect_equal(out$intercept, -6.5) + expect_equal(out$colocP12, 5e-6) + expect_equal(out$effectiveMiRounds, 25) }) -test_that(".enrFlattenEnrichment: non-empty data.frame resolves alternate column names", { - df <- data.frame(Enrichment = 3.0, se = 0.2, log_odds = 1.1) - out <- pecotmr:::.enrFlattenEnrichment(df) - expect_equal(out$enrichment, 3.0) - expect_equal(out$enrichmentSe, 0.2) - expect_equal(out$enrichmentLogOdds, 1.1) +test_that(".enrFlattenEnrichment: a zero log-odds leaves p12 unscaled", { + mock <- .qep_mockEnrichment(2.0)() + mock[["Enrichment (w/ shrinkage)"]] <- 0 + expect_equal(pecotmr:::.enrFlattenEnrichment(mock)$enrichment, 0) }) - -# =========================================================================== -# .enrPickColumn(): candidate column resolution -# =========================================================================== - -test_that(".enrPickColumn: returns the first matching column's value", { - df <- data.frame(se = 0.5, enrichment = 2.5) - expect_equal( - pecotmr:::.enrPickColumn(df, c("enrichment", "Enrichment")), - 2.5 +test_that(".enrFlattenEnrichment: an unrecognized shape warns, not NAs", { + # Silent all-NA here is exactly the failure that hid a field-name + # mismatch between this pipeline and the C++ estimator. + expect_warning( + out <- pecotmr:::.enrFlattenEnrichment(1.5), + "returned no 'Enrichment \\(w/ shrinkage\\)' field" ) -}) - -test_that(".enrPickColumn: returns NA_real_ when no candidate matches", { - df <- data.frame(foo = 1, bar = 2) - expect_identical( - pecotmr:::.enrPickColumn(df, c("enrichment", "Enrichment")), - NA_real_ + expect_equal(out, pecotmr:::.enrNaEnrichment()) + expect_warning( + out <- pecotmr:::.enrFlattenEnrichment(list(enrichment = 2.0)), + "returned no 'Enrichment" ) + expect_true(is.na(out$enrichment)) }) @@ -704,7 +850,7 @@ test_that(".enrPickColumn: returns NA_real_ when no candidate matches", { test_that(".enrBuildGwasPipVector: unknown study yields numeric(0)", { gfmr <- .qep_makeGwasFmr() expect_identical( - pecotmr:::.enrBuildGwasPipVector(gfmr, "GHOST"), + pecotmr:::.enrBuildGwasPipVector(gfmr, .qep_gwasIdent("GHOST")), numeric(0) ) }) @@ -728,7 +874,10 @@ test_that(".enrBuildGwasPipVector: skips a block whose pip length disagrees with entry = list(badEntry), ldSketch = .qep_makeHandle() ) - expect_identical(pecotmr:::.enrBuildGwasPipVector(g, "G1"), numeric(0)) + expect_identical( + pecotmr:::.enrBuildGwasPipVector(g, .qep_gwasIdent("G1")), + numeric(0) + ) }) @@ -755,7 +904,10 @@ test_that(".enrBuildQtlRegionsList: skips an entry whose fit lacks alpha/pip", { ldSketch = NULL ) expect_equal( - length(pecotmr:::.enrBuildQtlRegionsList(qfmr, "Q1", "c1")), + length(pecotmr:::.enrBuildQtlRegionsList( + qfmr, + .qep_qtlIdent("Q1", "c1") + )), 0L ) }) @@ -783,7 +935,10 @@ test_that(".enrBuildQtlRegionsList: skips a fit with no V and no prior_variance" ldSketch = NULL ) expect_equal( - length(pecotmr:::.enrBuildQtlRegionsList(qfmr, "Q1", "c1")), + length(pecotmr:::.enrBuildQtlRegionsList( + qfmr, + .qep_qtlIdent("Q1", "c1") + )), 0L ) }) @@ -811,7 +966,7 @@ test_that(".enrBuildQtlRegionsList: names an unnamed pip from the entry's varian entry = list(entry), ldSketch = NULL ) - out <- pecotmr:::.enrBuildQtlRegionsList(qfmr, "Q1", "c1") + out <- pecotmr:::.enrBuildQtlRegionsList(qfmr, .qep_qtlIdent("Q1", "c1")) expect_equal(length(out), 1L) expect_equal(names(out[[1L]]$pip), c("chr1:100:A:G", "chr1:200:A:G")) expect_equal(out[[1L]]$prior_variance, 0.1) @@ -838,7 +993,7 @@ test_that("qtlEnrichmentPipeline: empty QTL collection errors with the no-triple gwasFineMappingResult = gfmr, qtlFineMappingResult = qfmrEmpty ), - "triples to compute" + "pairs to compute" ) }) @@ -897,7 +1052,7 @@ test_that("qtlEnrichmentPipeline: a tuple with no usable QTL regions warns and i gwasFineMappingResult = gfmr, qtlFineMappingResult = qfmr ), - "no usable QTL regions" + "no usable regions" ) expect_equal(nrow(out), 0L) })