fix ctwas bug

R/GenotypeHandle.R

@@ -130,6 +130,16 @@ setMethod("show", "GenotypeHandle", function(object) {
 #'   \code{#chr,path} meta file or a named character vector
 #'   (names = chromosomes, values = payload paths/prefixes). Optionally pass
 #'   \code{format} via \code{...} to force a single backend for every shard.
+#' @param chroms Optional character vector of chromosomes. With \code{genoMeta}
+#'   (a sharded, one-file-per-chromosome panel) only the shards whose
+#'   chromosome is listed are read, skipping the rest -- an I/O optimisation
+#'   when the panel is genome-wide but only a few chromosomes are needed.
+#'   Chromosome labels are compared canonically (\code{"chr1"}/\code{"1"} match,
+#'   \code{23}/\code{X} etc.). If none of the requested chromosomes are present
+#'   in the panel, every shard is read (so the caller's own absence check can
+#'   report the mismatch). Only meaningful with \code{genoMeta}; supplying it
+#'   with any other source is an error (a single-file panel has no per-chromosome
+#'   shards to skip).
 #' @param ... Additional arguments forwarded to the format-specific reader.
 #' @return A \code{GenotypeHandle} object.
 #' @export
@@ -138,7 +148,7 @@ GenotypeHandle <- function(path = NULL,
                            bed = NULL, bim = NULL, fam = NULL,
                            pgen = NULL, pvar = NULL, psam = NULL,
                            ldMeta = NULL, region = NULL,
-                           genoMeta = NULL,
+                           genoMeta = NULL, chroms = NULL,
                            ...) {
   bedTrioGiven <- !is.null(bed) || !is.null(bim) || !is.null(fam)
   bedTrioComplete <- !is.null(bed) && !is.null(bim) && !is.null(fam)
@@ -173,6 +183,10 @@ GenotypeHandle <- function(path = NULL,
          "bed/bim/fam triplet, the pgen/pvar/psam triplet, `ldMeta`, or ",
          "`genoMeta` must be specified (got ", nSources, ").")
   }
+  if (!is.null(chroms) && !sources[["genoMeta"]]) {
+    stop("`chroms` restricts which per-chromosome shards are read and is only ",
+         "supported with `genoMeta` (a single-file panel has no shards to skip).")
+  }
 
   if (sources[["path"]]) {
     return(readGenotypes(path, ...))
@@ -193,7 +207,7 @@ GenotypeHandle <- function(path = NULL,
     return(.genotypeHandleFromLdMeta(ldMeta, region, ...))
   }
   if (sources[["genoMeta"]]) {
-    return(.genotypeHandleFromChromMeta(genoMeta, ...))
+    return(.genotypeHandleFromChromMeta(genoMeta, chroms = chroms, ...))
   }
 }
 
@@ -393,14 +407,28 @@ GenotypeHandle <- function(path = NULL,
 # metadata via the existing single-file readers, validates a single shared
 # format and identical sample IDs (same order, required for cross-shard
 # cbind), and row-binds the per-shard snpInfo into one global index space.
+# `chroms` (optional) restricts the read to the shards for those chromosomes:
+# the other per-chromosome files are never opened, which is the I/O win when a
+# genome-wide panel backs summary statistics on only a few chromosomes.
 #' @keywords internal
-.genotypeHandleFromChromMeta <- function(genoMeta, ...) {
+.genotypeHandleFromChromMeta <- function(genoMeta, chroms = NULL, ...) {
   dots   <- list(...)
   format <- dots$format
   parsed <- .parseChromMeta(genoMeta)
   if (nrow(parsed) == 0L)
     stop("GenotypeHandle(genoMeta): no chromosomes found in the meta input.")
 
+  # Keep only the shards whose (meta-declared) chromosome was requested, so the
+  # rest are never read. If nothing matches -- the requested chromosomes are
+  # absent from the panel -- fall back to reading every shard, leaving the
+  # caller's own containment/absence check to produce the usual diagnostic
+  # instead of a confusing empty-handle error here.
+  if (!is.null(chroms)) {
+    keep <- canonChrom(as.character(parsed$chrom)) %in%
+            canonChrom(as.character(chroms))
+    if (any(keep)) parsed <- parsed[keep, , drop = FALSE]
+  }
+
   shards <- lapply(parsed$path, .resolveGenotypeShard, format = format)
 
   formats <- vapply(shards, function(h) h@format, character(1))

R/colocboostPipeline.R

@@ -81,6 +81,14 @@
 #'   skipped. \code{0} (default) disables it; a negative value uses
 #'   \code{3 / n_variants}. (Summary-statistic skipping is handled upstream by
 #'   \code{\link{summaryStatsQc}}'s own \code{pipCutoffToSkip}.)
+#' @param absZCutoffToSkip,bfCutoffToSkip,logBfCutoffToSkip Alternative
+#'   individual-level pre-filter metrics used in place of
+#'   \code{pipCutoffToSkip}: drop an outcome unless its maximum marginal
+#'   \code{|z|} (\code{absZCutoffToSkip}), or its maximum per-variant
+#'   single-effect Bayes factor (\code{bfCutoffToSkip}) / log Bayes factor
+#'   (\code{logBfCutoffToSkip}) from the \code{L = 1} fit, exceeds the cutoff.
+#'   Scalars, each defaulting to 0 (off). Exactly one screening metric may be
+#'   enabled: setting any of these requires \code{pipCutoffToSkip = 0}.
 #' @param alleleFlip Logical, default \code{TRUE}. When TRUE, harmonize variants
 #'   across the individual X, sumstats, and LD by (chrom, pos) with ref/alt
 #'   swaps recognized (flipping z / residualized dosage / LD to a shared
@@ -160,29 +168,49 @@ setGeneric("colocboostPipeline",
   invisible(NULL)
 }
 
-# Resolve the per-context pipCutoffToSkip from either a scalar (applies to
-# every context) or a named vector keyed by context. Default 0 (no skip).
+# Resolve the per-context screen spec. The `pipCutoffToSkip` channel carries
+# EITHER a resolved screen object (list(metric, cutoff) for absZ/bf/logBf/pip --
+# applies uniformly to every context) OR the legacy PIP cutoff as a scalar (all
+# contexts) or a named vector keyed by context. Default 0 (no screen).
 .cbResolveCutoff <- function(pipCutoffToSkip, ctx) {
   if (is.null(pipCutoffToSkip) || length(pipCutoffToSkip) == 0L) return(0)
+  if (is.list(pipCutoffToSkip)) return(pipCutoffToSkip)   # uniform screen object
   if (!is.null(names(pipCutoffToSkip))) {
     if (ctx %in% names(pipCutoffToSkip)) return(pipCutoffToSkip[[ctx]])
     return(0)
   }
   pipCutoffToSkip[[1L]]
 }
 
+# Combine the four colocboost screen cutoffs into a single spec to thread
+# through the pipCutoffToSkip channel: a resolved screen object when a new
+# metric (absZ / bf / logBf) is set, otherwise the (possibly context-named)
+# legacy pipCutoffToSkip. Enforces one screening metric at a time.
+.cbScreenSpec <- function(pipCutoffToSkip, absZCutoffToSkip,
+                          bfCutoffToSkip, logBfCutoffToSkip) {
+  newScreen <- .resolveScreenMetric(0, absZCutoffToSkip,
+                                    bfCutoffToSkip, logBfCutoffToSkip)
+  pipOn <- !is.null(pipCutoffToSkip) && length(pipCutoffToSkip) > 0L &&
+           any(as.numeric(pipCutoffToSkip) != 0, na.rm = TRUE)
+  if (!is.null(newScreen) && pipOn)
+    stop("colocboostPipeline: only one signal screen may be enabled at a time; ",
+         "unset pipCutoffToSkip to use absZCutoffToSkip / bfCutoffToSkip / ",
+         "logBfCutoffToSkip.")
+  if (!is.null(newScreen)) newScreen else pipCutoffToSkip
+}
+
 # Per-outcome single-trait skip (ports the legacy qc_individual_data
 # pip_cutoff_to_skip): for each outcome column of Y, fit a single-effect
 # SuSiE (L = 1, max_iter = 100) on (X, Y[, j]) and keep the outcome only if
 # any variant's PIP exceeds the cutoff. A cutoff < 0 means 3 / n_variants.
 # Returns the retained Y (NULL when no outcome clears the threshold).
-.cbPipSkipOutcomes <- function(X, Y, cutoff) {
-  if (is.null(cutoff) || is.na(cutoff) || cutoff == 0) return(Y)
+.cbPipSkipOutcomes <- function(X, Y, spec) {
+  if (is.null(.asScreen(spec))) return(Y)
   # Single-effect screen per outcome, sharing the L = 1 SuSiE pre-screen
   # (.fmSerScreen) with the fine-mapping pipeline. fallback = FALSE: an outcome
   # that cannot be screened (too few samples / fit failure) is dropped.
   keep <- vapply(seq_len(ncol(Y)),
-                 function(j) .fmSerScreen(X, Y[, j], cutoff, fallback = FALSE),
+                 function(j) .fmSerScreen(X, Y[, j], spec, fallback = FALSE),
                  logical(1L))
   if (!any(keep)) return(NULL)
   Y[, keep, drop = FALSE]
@@ -251,11 +279,11 @@ setGeneric("colocboostPipeline",
     X <- X[common, , drop = FALSE]
     Y <- Y[common, , drop = FALSE]
     cutoffCtx <- .cbResolveCutoff(pipCutoffToSkip, ctx)
-    if (!is.null(cutoffCtx) && !is.na(cutoffCtx) && cutoffCtx != 0) {
+    if (!is.null(.asScreen(cutoffCtx))) {
       Y <- .cbPipSkipOutcomes(X, Y, cutoffCtx)
       if (is.null(Y) || ncol(Y) == 0L) {
         message("colocboostPipeline: skipping context '", ctx,
-                "' (no outcome cleared pipCutoffToSkip = ", cutoffCtx, ").")
+                "' (no outcome cleared the signal screen).")
         next
       }
     }
@@ -699,17 +727,22 @@ setMethod("colocboostPipeline", "QtlDataset",
            separateGwas = FALSE,
            samples = NULL,
            pipCutoffToSkip = 0,
+           absZCutoffToSkip = 0,
+           bfCutoffToSkip = 0,
+           logBfCutoffToSkip = 0,
            alleleFlip = TRUE,
            ...) {
     dotArgs <- list(...)
+    screenSpec <- .cbScreenSpec(pipCutoffToSkip, absZCutoffToSkip,
+                                bfCutoffToSkip, logBfCutoffToSkip)
     indBundle <- .cbIndividualBundle(
       qtlData,
       contexts     = contexts,
       traitId      = traitId,
       region       = region,
       cisWindow    = cisWindow,
       samples      = samples,
-      pipCutoffToSkip = pipCutoffToSkip)
+      pipCutoffToSkip = screenSpec)
     .cbDriver(indBundle, qtlPairs = list(), gwasSumStats,
               xqtlColoc, jointGwas, separateGwas,
               focalTrait, dotArgs, alleleFlip = alleleFlip)
@@ -755,9 +788,14 @@ setMethod("colocboostPipeline", "MultiStudyQtlDataset",
            separateGwas = FALSE,
            samples = NULL,
            pipCutoffToSkip = 0,
+           absZCutoffToSkip = 0,
+           bfCutoffToSkip = 0,
+           logBfCutoffToSkip = 0,
            alleleFlip = TRUE,
            ...) {
     dotArgs <- list(...)
+    screenSpec <- .cbScreenSpec(pipCutoffToSkip, absZCutoffToSkip,
+                                bfCutoffToSkip, logBfCutoffToSkip)
 
     # Aggregate the individual-level bundles across all QtlDataset
     # members. Per-study trait names are prefixed with "{study}:" so
@@ -777,7 +815,7 @@ setMethod("colocboostPipeline", "MultiStudyQtlDataset",
         region       = region,
         cisWindow    = cisWindow,
         samples      = samples,
-        pipCutoffToSkip = pipCutoffToSkip)
+        pipCutoffToSkip = screenSpec)
       if (is.null(sub)) next
       xOffset <- length(combinedX)
       yOffset <- length(combinedY)

R/ctwasPipeline.R

@@ -747,8 +747,26 @@ mergeCtwasBoundaryRegions <- function(finemapResult,
                               groupPriorVarStructure, thin, ncore,
                               extra = list()) {
   fitEm <- getFromNamespace("fit_EM", "ctwas")
+  # Mirror ctwas::est_param's degenerate-region skip before fitting. est_param
+  # drops regions with fewer than `min_var` total variables or fewer than
+  # `min_gene` genes -- whose per-region `sid` is unset -- and runs its prefit
+  # fit_EM only on the survivors. This prefit-only fallback bypasses est_param's
+  # p(single effect) selection gate, but it must apply the SAME filter, or
+  # ctwas::fit_EM errors inside extract_region_data ("regiondata$sid ... target
+  # is NULL") on a skipped region. Honor min_var / min_gene forwarded via `...`.
+  minVar  <- if (!is.null(extra$min_var))  as.integer(extra$min_var)  else 2L
+  minGene <- if (!is.null(extra$min_gene)) as.integer(extra$min_gene) else 1L
+  allRegionIds <- names(region_data)
+  nGid <- vapply(region_data, function(x) length(x$gid), integer(1L))
+  nSid <- vapply(region_data, function(x) length(x$sid), integer(1L))
+  keep <- rep(TRUE, length(region_data))
+  if (minVar  > 0L) keep <- keep & (nSid + nGid) >= minVar
+  if (minGene > 0L) keep <- keep & nGid >= minGene
+  fitRegionData <- region_data[keep]
+  if (length(fitRegionData) == 0L)
+    stop("No regions selected!")
   fitArgs <- list(
-    region_data               = region_data,
+    region_data               = fitRegionData,
     niter                     = as.integer(niterPrefit),
     group_prior_var_structure = groupPriorVarStructure,
     ncore                     = as.integer(ncore))
@@ -771,7 +789,7 @@ mergeCtwasBoundaryRegions <- function(finemapResult,
     group_prior_var_structure = groupPriorVarStructure,
     group_size                = groupSize,
     p_single_effect           = data.frame(
-      region_id        = names(region_data),
+      region_id        = allRegionIds,
       p_single_effect  = NA_real_,
       stringsAsFactors = FALSE))
 }