Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
66 changes: 64 additions & 2 deletions R/fineMappingPipeline.R
Original file line number Diff line number Diff line change
Expand Up @@ -265,6 +265,38 @@
#' passing (purity/coverage-filtered) credible set are widened. When
#' \code{TRUE}, every effect \code{L} is widened (including filtered-out ones,
#' labelled \code{L<k>} instead of \code{cs<k>}).
#' @param serFallback Logical (length 1, default \code{FALSE}).
#' \code{QtlSumStats} / \code{GwasSumStats} (summary-statistics SuSiE-RSS
#' paths) only. When \code{TRUE}, after each standard multi-effect SuSiE-RSS
#' fit (\code{susie} / \code{susieAsh}) the pipeline reads susieR's
#' finite-sample R diagnostics and, if
#' \code{fit$R_finite_diagnostics$R_reliability_flag} is \code{TRUE}, reports
#' the single-effect (\code{ser_model}) result for that region instead of the
#' multi-effect fit. Defaults \code{FALSE} so existing callers are unchanged.
#' See \code{keepFullFit} for retaining the multi-effect fit.
#' @param rFinite \code{QtlSumStats} / \code{GwasSumStats} only. Finite-sample
#' size for susieR's \code{R_finite} correction, forwarded to
#' \code{susieR::susie_rss()}. \code{NULL} (default) uses susieR's default,
#' except when a finite/EB mode is active (\code{serFallback=TRUE} or
#' \code{rMismatch != "none"}) and \code{rFinite} is \code{NULL}, in which
#' case it defaults to the LD-panel sample size \code{getNSamples(ldSketch)}.
#' @param rMismatch \code{QtlSumStats} / \code{GwasSumStats} only. LD-mismatch
#' mode forwarded to \code{susieR::susie_rss()} as \code{R_mismatch} (e.g.
#' \code{"eb"} for empirical Bayes). Default \code{"none"} (susieR's default).
#' @param rMismatchMethod \code{QtlSumStats} / \code{GwasSumStats} only. Optional
#' \code{R_mismatch_method} forwarded to \code{susieR::susie_rss()} when
#' non-\code{NULL}.
#' @param checkPrior \code{QtlSumStats} / \code{GwasSumStats} only. Optional
#' \code{check_prior} forwarded to \code{susieR::susie_rss()} when
#' non-\code{NULL}.
#' @param keepFullFit \code{QtlSumStats} / \code{GwasSumStats} only. Controls
#' retention of the pre-fallback multi-effect SuSiE-RSS fit when
#' \code{serFallback=TRUE}: \code{"fallback"} (default) keeps it only for
#' regions that fell back to SER; \code{"all"} keeps it for every region;
#' \code{"none"} keeps none. The retained fit and the decision are stored on
#' the entry's SuSiE fit and read via \code{getSusieFit(res)$multiEffectFit},
#' \code{getSusieFit(res)$R_reliability_flag}, and
#' \code{getSusieFit(res)$serFallbackUsed}.
#' @param ... Reserved for future per-method arguments.
#'
#' @return A \code{\link{FineMappingResult}} collection keyed by
Expand Down Expand Up @@ -1716,6 +1748,12 @@ setMethod("fineMappingPipeline", "QtlSumStats",
fullFit = FALSE,
fullFitAlphaOnly = TRUE,
includeAllCs = FALSE,
serFallback = FALSE,
rFinite = NULL,
rMismatch = "none",
rMismatchMethod = NULL,
checkPrior = NULL,
keepFullFit = "fallback",
...) {
.fmAssertQcd(data)
parsedJointSpec <- parseJointSpecification(jointSpecification, data)
Expand Down Expand Up @@ -1776,6 +1814,12 @@ setMethod("fineMappingPipeline", "QtlSumStats",
}

ldSketch <- getLdSketch(data)
# When a finite-sample R / EB LD-mismatch mode is active (serFallback on, or
# rMismatch other than "none") and rFinite is unset, default it to the
# LD-panel sample size (mirrors the GwasSumStats method).
rFiniteResolved <- if (is.null(rFinite) &&
(isTRUE(serFallback) || !identical(rMismatch, "none")))
getNSamples(ldSketch) else rFinite

rowStudy <- character(0)
rowContext <- character(0)
Expand Down Expand Up @@ -1826,7 +1870,10 @@ setMethod("fineMappingPipeline", "QtlSumStats",
signalCutoff, minAbsCorr, methodArgs, verbose,
label = sprintf("(study='%s', context='%s', trait='%s')", st, ctx, tr),
af = afByVar, fullFit = fullFit,
fullFitAlphaOnly = fullFitAlphaOnly, includeAllCs = includeAllCs)
fullFitAlphaOnly = fullFitAlphaOnly, includeAllCs = includeAllCs,
serFallback = serFallback, rFinite = rFiniteResolved,
rMismatch = rMismatch, rMismatchMethod = rMismatchMethod,
checkPrior = checkPrior, keepFullFit = keepFullFit)
# The method column carries the bare token, independent of the
# postprocess class.
for (tk in names(ents)) pushRow(st, ctx, tr, tk, ents[[tk]])
Expand Down Expand Up @@ -1903,6 +1950,12 @@ setMethod("fineMappingPipeline", "GwasSumStats",
fullFit = FALSE,
fullFitAlphaOnly = TRUE,
includeAllCs = FALSE,
serFallback = FALSE,
rFinite = NULL,
rMismatch = "none",
rMismatchMethod = NULL,
checkPrior = NULL,
keepFullFit = "fallback",
...) {
.fmAssertQcd(data)
norm <- .fmNormalizeMethods(methods, L = L, Lgreedy = Lgreedy)
Expand All @@ -1917,6 +1970,12 @@ setMethod("fineMappingPipeline", "GwasSumStats",
# across the entire entry's variant set; no in-pipeline LD-block
# partitioning.
ldSketch <- getLdSketch(data)
# When a finite-sample R / EB LD-mismatch mode is active (serFallback on, or
# rMismatch other than "none") and rFinite is unset, default it to the
# LD-panel sample size (the notebook's `B`).
rFiniteResolved <- if (is.null(rFinite) &&
(isTRUE(serFallback) || !identical(rMismatch, "none")))
getNSamples(ldSketch) else rFinite
studyCol <- as.character(data$study)

rowStudy <- character(0)
Expand Down Expand Up @@ -1977,7 +2036,10 @@ setMethod("fineMappingPipeline", "GwasSumStats",
signalCutoff, minAbsCorr, methodArgs, verbose,
label = sprintf("GWAS (study='%s', region='%s')", st, region_id),
af = afByVar, fullFit = fullFit,
fullFitAlphaOnly = fullFitAlphaOnly, includeAllCs = includeAllCs)
fullFitAlphaOnly = fullFitAlphaOnly, includeAllCs = includeAllCs,
serFallback = serFallback, rFinite = rFiniteResolved,
rMismatch = rMismatch, rMismatchMethod = rMismatchMethod,
checkPrior = checkPrior, keepFullFit = keepFullFit)
for (tk in names(ents)) pushRow(st, tk, region_id, ents[[tk]])
}

Expand Down
62 changes: 57 additions & 5 deletions R/fineMappingWrappers.R
Original file line number Diff line number Diff line change
Expand Up @@ -1990,13 +1990,24 @@ mergeSusieCs <- function(fineMappingResult, coverage = 0.95) {
# the same unmappable_effects switch, chained init, and userArgs merge.
# @noRd
.fmFitSusieRss <- function(z, R, n, token, chainFromInf = NULL,
coverage = 0.95, userArgs = NULL) {
coverage = 0.95, userArgs = NULL,
rFinite = NULL, rMismatch = "none",
rMismatchMethod = NULL, checkPrior = NULL) {
info <- .fineMappingMethodCapabilities[[token]]
if (is.null(info) || identical(info$unmappableEffects, NA_character_)) {
stop(".fmFitSusieRss: token '", token, "' is not a SuSiE-family method.")
}
baseArgs <- list(z = z, R = R, n = n, coverage = coverage,
unmappable_effects = info$unmappableEffects)
# Finite-sample R / EB LD-mismatch knobs (see the GWAS SER-fallback path).
# rFinite = NULL removes the element -> susie_rss default; rMismatch = "none"
# is the susie_rss default. These sit in baseArgs so they survive the chained
# branch's modifyList and the non-chained branch's userArgs merge, while user
# methodArgs (folded in after) still override them.
baseArgs$R_finite <- rFinite
baseArgs$R_mismatch <- rMismatch
if (!is.null(rMismatchMethod)) baseArgs$R_mismatch_method <- rMismatchMethod
if (!is.null(checkPrior)) baseArgs$check_prior <- checkPrior
if (!is.null(chainFromInf) && token != "susieInf") {
# SuSiE-RSS(-ash) initialised from a SuSiE-inf fit; userArgs folded into the
# arg prep so L_greedy is clamped rather than passed through raw.
Expand Down Expand Up @@ -2102,17 +2113,26 @@ mergeSusieCs <- function(fineMappingResult, coverage = 0.95) {
secondaryCoverage, signalCutoff, minAbsCorr,
methodArgs, verbose, label, af = NULL,
fullFit = FALSE, fullFitAlphaOnly = TRUE,
includeAllCs = FALSE) {
includeAllCs = FALSE,
serFallback = FALSE, rFinite = NULL,
rMismatch = "none", rMismatchMethod = NULL,
checkPrior = NULL, keepFullFit = "fallback") {
chainLocal <- .fmResolveSusieChain(toRun, addSusieInf)
infFit <- NULL
if (chainLocal$runInf) {
if (verbose >= 1)
message(sprintf("Fitting susieInf (RSS) for %s ...", label))
infFit <- .fmFitSusieRss(z, R, n, "susieInf", coverage = coverage,
userArgs = methodArgs[["susieInf"]])
userArgs = methodArgs[["susieInf"]],
rFinite = rFinite, rMismatch = rMismatch,
rMismatchMethod = rMismatchMethod,
checkPrior = checkPrior)
}
out <- list()
for (tk in toRun) {
multiFit <- NULL
flag <- NA
isStd <- FALSE
if (tk == "susieInf") {
if (!chainLocal$keepInf) next
fit <- infFit
Expand All @@ -2124,20 +2144,52 @@ mergeSusieCs <- function(fineMappingResult, coverage = 0.95) {
fit <- .fmFitSusieSer(z, n, coverage = coverage,
userArgs = methodArgs[["ser"]])
} else {
# Standard multi-effect SuSiE-RSS (susie / susieAsh). This is the only
# branch that can carry susieR's finite-sample R diagnostics and honour
# the SER fallback.
isStd <- TRUE
chainFrom <- if ((tk == "susie" && chainLocal$chainSusie) ||
(tk == "susieAsh" && chainLocal$chainAsh))
infFit else NULL
if (verbose >= 1)
message(sprintf("Fitting %s (RSS) for %s ...", tk, label))
fit <- .fmFitSusieRss(z, R, n, tk, chainFromInf = chainFrom,
coverage = coverage, userArgs = methodArgs[[tk]])
coverage = coverage, userArgs = methodArgs[[tk]],
rFinite = rFinite, rMismatch = rMismatch,
rMismatchMethod = rMismatchMethod,
checkPrior = checkPrior)
rfd <- fit$R_finite_diagnostics
flag <- if (!is.null(rfd) && !is.null(rfd$R_reliability_flag))
isTRUE(rfd$R_reliability_flag) else NA
if (isTRUE(serFallback) && isTRUE(flag) && !is.null(rfd$ser_model)) {
# susieR flagged the multi-effect fit as unreliable: report the
# single-effect (ser_model) result instead, keeping the multi-effect
# fit aside for retention below.
multiFit <- fit
fit <- .setFinemappingFitClass(rfd$ser_model, "susieRss")
}
}
out[[tk]] <- .fmPostprocessOne(
ent <- .fmPostprocessOne(
fit = fit, method = "susieRss", dataX = R, dataY = list(z = z),
coverage = coverage, secondaryCoverage = secondaryCoverage,
signalCutoff = signalCutoff, minAbsCorr = minAbsCorr, af = af,
csInput = "Xcorr", fullFit = fullFit,
fullFitAlphaOnly = fullFitAlphaOnly, includeAllCs = includeAllCs)
if (isStd && isTRUE(serFallback)) {
# Record the reliability decision (and, per keepFullFit, the retained
# multi-effect fit) on the entry's susieFit list. Gated on serFallback so
# that with the feature off (the default on every sumstat path) the entry
# is byte-identical to before this change (no extra diagnostic fields).
sf <- ent@susieFit
sf$R_reliability_flag <- flag
sf$serFallbackUsed <- isTRUE(flag)
if (!is.null(multiFit) && keepFullFit %in% c("fallback", "all"))
sf$multiEffectFit <- multiFit
else if (identical(keepFullFit, "all"))
sf$multiEffectFit <- fit
ent@susieFit <- sf
}
out[[tk]] <- ent
}
out
}
50 changes: 50 additions & 0 deletions man/fineMappingPipeline.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

Loading
Loading