pipeline fixes

DESCRIPTION

@@ -20,7 +20,6 @@ Imports:
     IRanges,
     MASS,
     Matrix,
-    MungeSumstats,
     Rsamtools,
     S4Vectors,
     SummarizedExperiment,

NAMESPACE

@@ -21,6 +21,7 @@ export(TwasWeightsEntry)
 export(adjustPips)
 export(alignVariantNames)
 export(alleleQc)
+export(assembleCtwasInputs)
 export(autoDecision)
 export(bLassoWeights)
 export(batchLoadTwasWeights)
@@ -61,6 +62,7 @@ export(enetWeights)
 export(enforceDesignFullRank)
 export(enlocPipeline)
 export(ensembleWeights)
+export(estCtwasParam)
 export(estimateH2)
 export(estimateSparsity)
 export(extractBlockGenotypes)
@@ -72,6 +74,7 @@ export(filterRelatedness)
 export(filterVariantsByLdReference)
 export(findOverlappingRegions)
 export(fineMappingPipeline)
+export(finemapCtwasRegions)
 export(fitFsusie)
 export(fitMashContrast)
 export(fitMvsusie)
@@ -225,6 +228,7 @@ export(rssAnalysisPipeline)
 export(sanitizeMashData)
 export(scadRssWeights)
 export(scadWeights)
+export(screenCtwasRegions)
 export(sdpr)
 export(sdprWeights)
 export(slalom)
@@ -352,7 +356,6 @@ importFrom(IRanges,DataFrameList)
 importFrom(IRanges,IRanges)
 importFrom(IRanges,findOverlaps)
 importFrom(MASS,ginv)
-importFrom(MungeSumstats,standardise_header)
 importFrom(Rsamtools,TabixFile)
 importFrom(Rsamtools,asBcf)
 importFrom(Rsamtools,headerTabix)

R/GwasFineMappingResult.R

@@ -1,9 +1,13 @@
 # =============================================================================
 # GwasFineMappingResult S4 class
 # -----------------------------------------------------------------------------
-# DFrame-subclass collection keyed by the identity tuple (study, method).
-# Each row holds a FineMappingEntry payload for one GWAS study at one
-# fine-mapping method, covering a single LD block. Class-level slots:
+# DFrame-subclass collection keyed by the identity tuple
+# (study, method, region_id). Each row holds a FineMappingEntry payload
+# for one GWAS study at one fine-mapping method over one LD block.
+# Multiple rows per (study, method) are allowed when they differ on
+# region_id — the genome-wide-across-blocks shape that
+# qtlEnrichmentPipeline / colocPipeline expect when sweeping a study
+# across LD blocks. Class-level slots:
 #   * ldSketch   GenotypeHandle for the LD reference; required for the
 #                LD-block-indexed susieRSS workflow.
 # Methods that take per-row selectors accept (study, method) and ignore
@@ -17,7 +21,7 @@ setClass("GwasFineMappingResult",
   contains = "FineMappingResultBase",
   validity = function(object) {
     errors <- character()
-    required <- c("study", "method", "entry")
+    required <- c("study", "method", "region_id", "entry")
     missingCols <- setdiff(required, names(object))
     if (length(missingCols) > 0L)
       errors <- c(errors, paste("missing columns:",
@@ -32,10 +36,10 @@ setClass("GwasFineMappingResult",
       if (!all(entryTypes))
         errors <- c(errors,
           "every element of the `entry` column must be a FineMappingEntry")
-      keyDf <- as.data.frame(object[, c("study", "method")])
+      keyDf <- as.data.frame(object[, c("study", "method", "region_id")])
       if (anyDuplicated(keyDf))
         errors <- c(errors,
-          "(study, method) tuple uniqueness violated")
+          "(study, method, region_id) tuple uniqueness violated")
     }
     if (!is.null(object@ldSketch) &&
         !methods::is(object@ldSketch, "GenotypeHandle")) {
@@ -77,26 +81,41 @@ setClass("GwasFineMappingResult",
 
 #' @title Create a GwasFineMappingResult Collection
 #' @description Construct a \code{GwasFineMappingResult} DFrame-subclass
-#'   collection from per-(study, method) tuples and a list of
-#'   \code{FineMappingEntry} payloads. The collection represents one LD
-#'   block of GWAS fine-mapping fits; build a separate collection per
-#'   block when sweeping the genome.
+#'   collection from per-(study, method, region_id) tuples and a list of
+#'   \code{FineMappingEntry} payloads. The collection can represent
+#'   either a single LD block (one row per (study, method)) or a
+#'   genome-wide sweep across blocks (multiple rows per (study, method),
+#'   each tagged with its own region_id).
 #' @param study Character vector of study identifiers (per tuple).
 #' @param method Character vector of fine-mapping method names (per tuple).
 #' @param entry List / \code{SimpleList} of \code{FineMappingEntry} objects.
+#' @param region_id Character vector of LD-block identifiers (per
+#'   tuple). When omitted (\code{NULL}), defaults to a per-row synthetic
+#'   id (\code{"region_1"}, \code{"region_2"}, ...) so the
+#'   (\code{study}, \code{method}, \code{region_id}) triple is unique.
+#'   Supplying meaningful labels (e.g. \code{"chr22_10516173_17414263"})
+#'   is preferred for downstream consumers that join on region.
 #' @param ldSketch An optional \code{GenotypeHandle}.
 #' @return A \code{GwasFineMappingResult} object.
 #' @export
 GwasFineMappingResult <- function(study, method, entry,
+                                  region_id = NULL,
                                   ldSketch = NULL) {
   n <- length(study)
   if (length(method) != n || length(entry) != n) {
     stop("`study`, `method`, and `entry` must all have the same length.")
   }
+  if (is.null(region_id)) {
+    region_id <- paste0("region_", seq_len(n))
+  } else if (length(region_id) != n) {
+    stop("`region_id` must have the same length as `study` (got ",
+         length(region_id), " vs ", n, ").")
+  }
   cols <- list(
-    study  = as.character(study),
-    method = as.character(method),
-    entry  = S4Vectors::SimpleList(entry)
+    study     = as.character(study),
+    method    = as.character(method),
+    region_id = as.character(region_id),
+    entry     = S4Vectors::SimpleList(entry)
   )
   df <- do.call(S4Vectors::DataFrame,
                 c(cols, list(check.names = FALSE)))
@@ -121,9 +140,10 @@ setMethod("getContexts", "GwasFineMappingResult", function(x) NULL)
 #' @export
 setMethod("getTraits", "GwasFineMappingResult", function(x) NULL)
 
-# Per-tuple lookup for the GWAS variant (2-tuple instead of 4-tuple).
-# The generic accepts the full set of selectors; context/trait args are
-# ignored for GwasFineMappingResult.
+# Per-tuple lookup keyed by (study, method, region_id). The generic
+# accepts the full set of selectors; context/trait args are ignored for
+# GwasFineMappingResult. `region` (passed via `...`) is the per-block
+# disambiguator for multi-block genome-wide collections.
 #' @rdname getFineMappingResult
 #' @export
 setMethod("getFineMappingResult", "GwasFineMappingResult",
@@ -136,52 +156,58 @@ setMethod("getFineMappingResult", "GwasFineMappingResult",
 #' @export
 setMethod("getPip", "GwasFineMappingResult",
   function(x, study = NULL, context = NULL, trait = NULL, method = NULL,
-           returnList = FALSE, ...) {
-    entry <- getFineMappingResult(x, study = study, method = method)
-    getPip(entry)
+           region = NULL, returnList = FALSE, ...) {
+    idx <- .tupleSelectRowGwasFmr(x, study, method, region)
+    getPip(x$entry[[idx]])
   })
 
 #' @rdname getCs
 #' @export
 setMethod("getCs", "GwasFineMappingResult",
-  function(x, study = NULL, context = NULL, trait = NULL, method = NULL, ...) {
-    entry <- getFineMappingResult(x, study = study, method = method)
-    getCs(entry)
+  function(x, study = NULL, context = NULL, trait = NULL, method = NULL,
+           region = NULL, ...) {
+    idx <- .tupleSelectRowGwasFmr(x, study, method, region)
+    getCs(x$entry[[idx]])
   })
 
 #' @rdname getTopLoci
 #' @export
 setMethod("getTopLoci", "GwasFineMappingResult",
   function(x, type = c("data.frame", "GRanges"),
            signalCutoff = 0.025,
-           study = NULL, context = NULL, trait = NULL, method = NULL, ...) {
-    entry <- getFineMappingResult(x, study = study, method = method)
-    getTopLoci(entry, type = match.arg(type), signalCutoff = signalCutoff)
+           study = NULL, context = NULL, trait = NULL, method = NULL,
+           region = NULL, ...) {
+    idx <- .tupleSelectRowGwasFmr(x, study, method, region)
+    getTopLoci(x$entry[[idx]], type = match.arg(type),
+               signalCutoff = signalCutoff)
   })
 
 #' @rdname getMarginalEffects
 #' @export
 setMethod("getMarginalEffects", "GwasFineMappingResult",
   function(x, maxPval = NULL,
-           study = NULL, context = NULL, trait = NULL, method = NULL, ...) {
-    entry <- getFineMappingResult(x, study = study, method = method)
-    getMarginalEffects(entry, maxPval = maxPval)
+           study = NULL, context = NULL, trait = NULL, method = NULL,
+           region = NULL, ...) {
+    idx <- .tupleSelectRowGwasFmr(x, study, method, region)
+    getMarginalEffects(x$entry[[idx]], maxPval = maxPval)
   })
 
 #' @rdname getSusieFit
 #' @export
 setMethod("getSusieFit", "GwasFineMappingResult",
-  function(x, study = NULL, context = NULL, trait = NULL, method = NULL, ...) {
-    entry <- getFineMappingResult(x, study = study, method = method)
-    getSusieFit(entry)
+  function(x, study = NULL, context = NULL, trait = NULL, method = NULL,
+           region = NULL, ...) {
+    idx <- .tupleSelectRowGwasFmr(x, study, method, region)
+    getSusieFit(x$entry[[idx]])
   })
 
 #' @rdname getVariantIds
 #' @export
 setMethod("getVariantIds", "GwasFineMappingResult",
-  function(x, study = NULL, context = NULL, trait = NULL, method = NULL, ...) {
-    entry <- getFineMappingResult(x, study = study, method = method)
-    getVariantIds(entry)
+  function(x, study = NULL, context = NULL, trait = NULL, method = NULL,
+           region = NULL, ...) {
+    idx <- .tupleSelectRowGwasFmr(x, study, method, region)
+    getVariantIds(x$entry[[idx]])
   })
 
 

R/colocPipeline.R

@@ -81,9 +81,10 @@
 #' @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, qtlContext)
-#'   enrichment factors with columns \code{gwasStudy}, \code{qtlContext},
-#'   \code{enrichment}. Output of \code{\link{qtlEnrichmentPipeline}}.
+#' @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{p12} prior is scaled to
 #'   \code{min(p12 * (1 + enrichment), p12Max)} (the enrichment-informed
 #'   colocalization variant, "enloc"). Pairs without a matching
@@ -132,9 +133,9 @@ colocPipeline <- function(qtlFineMappingResult,
   if (useEnrichment) {
     if (!is.data.frame(enrichment))
       stop("`enrichment` must be a data.frame with at least gwasStudy, ",
-           "qtlContext, enrichment columns (output of ",
+           "qtlStudy, qtlContext, enrichment columns (output of ",
            "qtlEnrichmentPipeline).")
-    required <- c("gwasStudy", "qtlContext", "enrichment")
+    required <- c("gwasStudy", "qtlStudy", "qtlContext", "enrichment")
     missingCols <- setdiff(required, colnames(enrichment))
     if (length(missingCols) > 0L)
       stop("`enrichment` is missing column(s): ",
@@ -231,14 +232,16 @@ colocPipeline <- function(qtlFineMappingResult,
       qAligned <- aligned$qtl
       gAligned <- aligned$gwas
 
-      # Enrichment-informed p12: per-pair scaling capped at p12Max.
-      # Baseline p12 used when no enrichment table or no matching row.
+      # Enrichment-informed p12: per-(gwasStudy, qtlStudy, qtlContext)
+      # scaling capped at p12Max. Baseline p12 used when no enrichment
+      # table or no matching row.
       if (useEnrichment) {
-        enRow <- .colocLookupEnrichment(enrichment, gInfo$study, qContext)
+        enRow <- .colocLookupEnrichment(enrichment, gInfo$study,
+                                        qStudy, qContext)
         if (is.na(enRow)) {
           warning(sprintf(
-            "colocPipeline: no enrichment entry for (gwasStudy='%s', qtlContext='%s'); using baseline p12.",
-            gInfo$study, qContext))
+            "colocPipeline: no enrichment entry for (gwasStudy='%s', qtlStudy='%s', qtlContext='%s'); using baseline p12.",
+            gInfo$study, qStudy, qContext))
           enRow <- 0
         }
         p12Used <- min(p12 * (1 + enRow), p12Max)
@@ -504,13 +507,14 @@ colocPipeline <- function(qtlFineMappingResult,
   base
 }
 
-# Look up the enrichment factor for a (gwasStudy, qtlContext) pair in
-# the user-supplied enrichment table. Returns NA when the pair is not
-# present; the caller falls back to the baseline p12 and emits a
-# warning.
+# 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.
 # @noRd
-.colocLookupEnrichment <- function(enrichment, gwasStudy, qtlContext) {
+.colocLookupEnrichment <- function(enrichment, gwasStudy, qtlStudy, qtlContext) {
   idx <- which(as.character(enrichment$gwasStudy)  == gwasStudy &
+               as.character(enrichment$qtlStudy)   == qtlStudy &
                as.character(enrichment$qtlContext) == qtlContext)
   if (length(idx) == 0L) return(NA_real_)
   as.numeric(enrichment$enrichment[[idx[[1L]]]])