Skip to content

Skip PIP-screened regions gracefully instead of erroring #565

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
danielnachun merged 1 commit into from
Aug 4, 2026
Merged

Skip PIP-screened regions gracefully instead of erroring #565

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
72 changes: 66 additions & 6 deletions R/fineMappingPipeline.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -107,6 +107,11 @@
#' carries the bare token (\code{"susie"},
#' \code{"susieInf"}, \code{"mvsusie"}, ...) only. QC
#' provenance is recorded on the sumstats' \code{qcInfo}.
#' \item An entry that \code{summaryStatsQc(pipCutoffToSkip = ...)} screened
#' out (recorded as \code{qcInfo$entryAudit[[i]]$pipScreenSkipped}, and
#' emptied to 0 variants) is \strong{skipped}, not fit: it produces no
#' row and a message with the screen reason. An all-screened collection
#' yields a valid empty result rather than an error.
#' }
#'
#' @param data A \code{QtlDataset}, \code{MultiStudyQtlDataset},
Expand Down Expand Up @@ -586,8 +591,9 @@ setGeneric("fineMappingPipeline",
jointTraits = NULL,
region = NULL,
traitPos = NULL,
ldSketch = NULL) {
if (length(entries) == 0L) {
ldSketch = NULL,
allowEmpty = FALSE) {
if (length(entries) == 0L && !allowEmpty) {
stop("fineMappingPipeline: no (study, context, trait, method) tuples ",
"produced a fine-mapping result.")
}
Expand All @@ -612,8 +618,9 @@ setGeneric("fineMappingPipeline",
# pass them explicitly so downstream consumers can join on region.
# @noRd
.fmBuildGwasResult <- function(studies, methods, entries,
region_ids = NULL, ldSketch = NULL) {
if (length(entries) == 0L) {
region_ids = NULL, ldSketch = NULL,
allowEmpty = FALSE) {
if (length(entries) == 0L && !allowEmpty) {
stop("fineMappingPipeline: no (study, method, region_id) tuples produced a ",
"fine-mapping result.")
}
Expand Down Expand Up @@ -1120,6 +1127,25 @@ combineFineMappingResults <- function(..., ldSketch = NULL) {
n = stats::median(df$N, na.rm = TRUE))
}

# Whether entry `i` of a QC'd SumStats was deliberately screened out (and why),
# so fineMappingPipeline can skip it gracefully instead of erroring on a
# 0-variant entry. `summaryStatsQc(pipCutoffToSkip = ...)` empties a no-signal
# region and records qcInfo$entryAudit[[i]]$pipScreenSkipped (+ pipScreenReason);
# an entry may also be empty for other reasons. Returns list(skipped, reason).
.fmEntrySkipInfo <- function(data, i) {
ea <- tryCatch(getQcInfo(data)$entryAudit[[i]], error = function(e) NULL)
screened <- isTRUE(ea$pipScreenSkipped)
entry <- data$entry[[i]]
empty <- is.null(entry) || length(entry) == 0L
reason <-
if (!is.null(ea$pipScreenReason) && nzchar(as.character(ea$pipScreenReason)))
as.character(ea$pipScreenReason)
else if (screened) "no signal above the PIP pre-screen cutoff"
else if (empty) "empty entry (no variants)"
else NA_character_
list(skipped = isTRUE(screened || empty), reason = reason)
}


# =============================================================================
# Per-fold cross-validation of fine-mapping methods
Expand Down Expand Up @@ -1825,6 +1851,7 @@ setMethod("fineMappingPipeline", "QtlSumStats",
rowTrait <- character(0)
rowMethod <- character(0)
rowEntries <- list()
nSkipped <- 0L
pushRow <- function(st, ctx, tr, mt, ent) {
rowStudy <<- c(rowStudy, st)
rowContext <<- c(rowContext, ctx)
Expand All @@ -1850,6 +1877,17 @@ setMethod("fineMappingPipeline", "QtlSumStats",
}
if (length(toRun) == 0L) next

# A trait screened out by summaryStatsQc(pipCutoffToSkip) is empty here;
# skip it gracefully (no row, no error) rather than tripping .fmExtractZN.
skip <- .fmEntrySkipInfo(data, i)
if (isTRUE(skip$skipped)) {
nSkipped <- nSkipped + 1L
if (verbose >= 1)
message(sprintf(
"fineMappingPipeline(QtlSumStats): entry %d (study='%s', context='%s', trait='%s') skipped: %s",
i, st, ctx, tr, skip$reason))
next
}
entry <- data$entry[[i]]
zn <- .fmExtractZN(entry,
sprintf("fineMappingPipeline(QtlSumStats): entry %d (study='%s', context='%s', trait='%s')", i, st, ctx, tr))
Expand Down Expand Up @@ -1917,8 +1955,15 @@ setMethod("fineMappingPipeline", "QtlSumStats",
ldSketch = ldSketch)
else NULL
if (is.null(jointResult)) {
if (is.null(perTupleResult))
if (is.null(perTupleResult)) {
# All traits screened out by summaryStatsQc(pipCutoffToSkip) -> a valid
# empty result, not an error.
if (nSkipped > 0L)
return(.fmBuildQtlResult(character(0), character(0), character(0),
character(0), list(), ldSketch = ldSketch,
allowEmpty = TRUE))
stop("fineMappingPipeline(QtlSumStats): no entries produced a result.")
}
return(perTupleResult)
}
if (is.null(perTupleResult)) return(jointResult)
Expand Down Expand Up @@ -1980,6 +2025,7 @@ setMethod("fineMappingPipeline", "GwasSumStats",
rowMethod <- character(0)
rowRegion <- character(0)
rowEntries <- list()
nSkipped <- 0L
pushRow <- function(st, mt, rg, ent) {
rowStudy <<- c(rowStudy, st)
rowMethod <<- c(rowMethod, mt)
Expand All @@ -1990,6 +2036,17 @@ setMethod("fineMappingPipeline", "GwasSumStats",
for (i in seq_len(nrow(data))) {
st <- studyCol[[i]]
gr <- data$entry[[i]]
# A region screened out by summaryStatsQc(pipCutoffToSkip) is empty here;
# skip it gracefully (no row, no error) rather than tripping .fmExtractZN.
skip <- .fmEntrySkipInfo(data, i)
if (isTRUE(skip$skipped)) {
nSkipped <- nSkipped + 1L
if (verbose >= 1)
message(sprintf(
"fineMappingPipeline(GwasSumStats): study='%s' region skipped: %s",
st, skip$reason))
next
}
zn <- .fmExtractZN(gr,
sprintf("fineMappingPipeline(GwasSumStats): study='%s'", st))
variantIds <- zn$variantIds
Expand Down Expand Up @@ -2041,9 +2098,12 @@ setMethod("fineMappingPipeline", "GwasSumStats",
for (tk in names(ents)) pushRow(st, tk, region_id, ents[[tk]])
}

# An all-screened (or empty-input) collection legitimately yields a 0-row
# result -- allow it instead of erroring "no ... tuples produced a result".
.fmBuildGwasResult(rowStudy, rowMethod, rowEntries,
region_ids = rowRegion,
ldSketch = ldSketch)
ldSketch = ldSketch,
allowEmpty = (nSkipped > 0L || nrow(data) == 0L))
})


Expand Down
64 changes: 64 additions & 0 deletions tests/testthat/test_fineMappingPipeline.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1515,6 +1515,70 @@ test_that("fineMappingPipeline(GwasSumStats): runs end-to-end with mocked RSS fi
expect_setequal(getMethodNames(res), "susie")
})

# ---- PIP-screen graceful skip: a screened region -> empty result, not error ----

# A GwasSumStats whose (single) entry was emptied by summaryStatsQc's PIP screen:
# 0-variant entry + qcInfo$entryAudit[[1]]$pipScreenSkipped = TRUE (+ reason).
.fmp_makeScreenedGwas <- function(study = "G1",
reason = "no signals above PIP threshold 0.025") {
GwasSumStats(
study = study,
entry = list(GenomicRanges::GRanges()),
genome = "hg19",
ldSketch = .fmp_makeHandle(),
qcInfo = list(entryAudit = list(list(pipScreenSkipped = TRUE,
pipScreenReason = reason))))
}

test_that("fineMappingPipeline(GwasSumStats): a PIP-screened region yields an empty result, not an error", {
gss <- .fmp_makeScreenedGwas()
expect_message(
res <- fineMappingPipeline(gss, methods = "susie", addSusieInf = FALSE),
"region skipped")
expect_s4_class(res, "GwasFineMappingResult")
expect_equal(nrow(res), 0L) # graceful skip, not stop("no ... tuples")
})

test_that("fineMappingPipeline(GwasSumStats): mixed screened + real keeps only the real region", {
gss <- GwasSumStats(
study = c("Gskip", "Greal"),
entry = list(GenomicRanges::GRanges(), .fmp_makeSumstatsGr()),
genome = "hg19",
ldSketch = .fmp_makeHandle(),
qcInfo = list(entryAudit = list(
list(pipScreenSkipped = TRUE, pipScreenReason = "no signal"), list())))
local_mocked_bindings(
extractBlockGenotypes = .fmp_mockExtractor(),
.fmFitSusieRss = .fmp_mockFitRss(),
.fmPostprocessOne = .fmp_mockPostprocess(),
.package = "pecotmr")
res <- suppressMessages(
fineMappingPipeline(gss, methods = "susie", addSusieInf = FALSE))
expect_equal(unique(as.character(res$study)), "Greal") # screened study absent
expect_gt(nrow(res), 0L)
})

test_that("fineMappingPipeline(GwasSumStats): a 0-variant entry is skipped even without the flag", {
gss <- GwasSumStats(study = "G1", entry = list(GenomicRanges::GRanges()),
genome = "hg19", ldSketch = .fmp_makeHandle(), qcInfo = list(step1 = "ok"))
res <- suppressMessages(
fineMappingPipeline(gss, methods = "susie", addSusieInf = FALSE))
expect_s4_class(res, "GwasFineMappingResult")
expect_equal(nrow(res), 0L)
})

test_that("fineMappingPipeline(QtlSumStats): a screened trait is skipped gracefully", {
qss <- QtlSumStats(study = "Q1", context = "c1", trait = "t1",
entry = list(GenomicRanges::GRanges()), genome = "hg19",
ldSketch = .fmp_makeHandle(),
qcInfo = list(entryAudit = list(list(pipScreenSkipped = TRUE,
pipScreenReason = "no signal"))))
res <- suppressMessages(
fineMappingPipeline(qss, methods = "susie", addSusieInf = FALSE))
expect_s4_class(res, "QtlFineMappingResult")
expect_equal(nrow(res), 0L)
})

# ---- SER fallback (shared sumstat SuSiE-RSS path: GwasSumStats + QtlSumStats) ----

test_that("fineMappingPipeline(GwasSumStats): serFallback + reliable R keeps multi-effect", {
Expand Down
Loading