fix LD sketch mismatch

R/causalInferencePipeline.R

@@ -16,11 +16,16 @@
 #'           MR computed from \code{fineMappingResult}.
 #'   }
 #'
-#' @section LD-sketch identity check: If a QTL input (TwasWeights or
+#' @section LD-sketch compatibility check: If a QTL input (TwasWeights or
 #'   QtlFineMappingResult) carries a non-\code{NULL} \code{ldSketch}, it must
-#'   match the \code{ldSketch} on \code{gwasSumStats}. Mismatch is a hard error.
-#'   A QTL input with \code{ldSketch = NULL} (the fit was learned from
-#'   individual-level data) skips the validation for that input.
+#'   come from the same reference panel as the \code{ldSketch} on
+#'   \code{gwasSumStats} --- same samples, same allele orientation on the
+#'   shared variants. The two need NOT carry the same variants: QC-ing the two
+#'   sides separately trims each sketch to its own surviving variants, and a
+#'   partial overlap only warns. A different sample set, a swapped A1/A2 on a
+#'   shared variant, or no shared variant at all is a hard error. A QTL input
+#'   with \code{ldSketch = NULL} (the fit was learned from individual-level
+#'   data) skips the validation for that input.
 #'
 #' @section Output shape:
 #' A long-format \code{GRanges} with one row per
@@ -466,7 +471,7 @@ causalInferencePipeline <- function(
 # Internal helpers
 # =============================================================================
 
-# Compare two GenotypeHandles for LD-sketch identity. Thin wrapper over
+# Check two LD sketches for reference-panel compatibility. Thin wrapper over
 # the shared `.requireMatchingLdSketches` helper (R/ld.R).
 .cipRequireMatchingLdSketches <- function(qtlLd, gwasLd, label) {
     .requireMatchingLdSketches(

R/colocPipeline.R

@@ -50,12 +50,18 @@
 #'           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
-#'   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.
+#' @section LD-sketch compatibility check: If
+#'   \code{getLdSketch(qtlFineMappingResult)} is non-\code{NULL}, it must come
+#'   from the same reference panel as the LD sketch on \code{gwasInput}: the
+#'   same samples, and the same allele orientation on the variants the two
+#'   carry in common. The two need NOT carry the same variants --- running
+#'   \code{\link{summaryStatsQc}} on the two sides separately normally leaves
+#'   each sketch trimmed to its own surviving variants, and LD is looked up per
+#'   variant, so a partial overlap only warns (once per session). No shared
+#'   variant at all, a different sample set, or a swapped A1/A2 on a shared
+#'   variant 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 The first side: a
 #'   \code{\link{QtlFineMappingResult}} or a
@@ -570,7 +576,7 @@ colocPipeline <- function(
 # Internal helpers
 # =============================================================================
 
-# LD-sketch identity check. Thin wrapper over the shared
+# LD-sketch compatibility check. Thin wrapper over the shared
 # `.requireMatchingLdSketches` helper (R/ld.R). Shared with
 # qtlEnrichmentPipeline.
 # @noRd

R/colocboostPipeline.R

@@ -718,7 +718,7 @@ setGeneric("colocboostPipeline", function(qtlData, gwasSumStats = NULL, ...) {
     bundle
 }
 
-# Compare two GenotypeHandles for the LD-sketch equality contract. Thin
+# Check two LD sketches for the reference-panel compatibility contract. Thin
 # wrapper over the shared `.requireMatchingLdSketches` helper (R/ld.R)
 # using the "lenient" null policy: a NULL on either side skips the check
 # (only colocboostPipeline allows that, since some bundles only have a

R/ctwasPipeline.R

@@ -17,9 +17,13 @@
 #'   cannot converge on a single region, so callers must supply at least two
 #'   blocks.
 #'
-#' @section LD-sketch identity check: Per block: \code{getLdSketch(twasWeights)}
-#'   (when non-NULL) must match \code{getLdSketch(gwasSumStats)}. Mismatch is a
-#'   hard error.
+#' @section LD-sketch compatibility check: Per block:
+#'   \code{getLdSketch(twasWeights)} (when non-NULL) must come from the same
+#'   reference panel as \code{getLdSketch(gwasSumStats)} --- same samples, same
+#'   allele orientation on the shared variants. The two need NOT carry the same
+#'   variants; a partial overlap, which is what QC-ing the two sides separately
+#'   produces, only warns. A different sample set, a swapped A1/A2 on a shared
+#'   variant, or no shared variant at all is a hard error.
 #'
 #' @param gwasSumStats A \code{\link{GwasSumStats}} whose elements are LD
 #'   blocks (at least two), keyed by its \code{blockId} column, with
@@ -1233,7 +1237,7 @@ mergeCtwasBoundaryRegions <- function(
 # Internal helpers
 # =============================================================================
 
-# LD-sketch identity check. Thin wrapper over the shared
+# LD-sketch compatibility check. Thin wrapper over the shared
 # `.requireMatchingLdSketches` helper (R/ld.R).
 .ctwasRequireMatchingLdSketches <- function(twLd, gwasLd) {
     .requireMatchingLdSketches(twLd, gwasLd, pipelineName = "ctwasPipeline")

R/ld.R

@@ -1581,16 +1581,21 @@ loadLdFromGenotype <- function(
     }
     list(mafCutoff = maf, macCutoff = mac, imissCutoff = imiss)
 }
-# ---------- LD sketch: cross-pipeline LD-panel equality check ----------
-
-# Internal: assert that two `GenotypeHandle` LD sketches describe the same
-# reference panel: same variant identity (chr-agnostic CHR via canonChrom,
-# exact BP/A1/A2, in the same order) and the same sampleIds. The SNP label is
-# not compared, so a pure chr-prefix difference does not fail; an allele swap
-# (different A1/A2) still does, since it means a different LD coding. Shared by
-# causalInferencePipeline, colocPipeline,
-# qtlEnrichmentPipeline, ctwasPipeline, and
-# colocboostPipeline.
+# ---------- LD sketch: cross-pipeline LD-panel compatibility check ----------
+
+# Internal: assert that two LD sketches describe the same reference panel --
+# the same samples, and the same allele coding (A1/A2 in the same orientation)
+# on the variants they share. The SNP label is not compared, so a pure
+# chr-prefix difference does not fail; an allele swap still does, since it
+# means the two panels code their LD in opposite directions.
+#
+# The two panels' variant SETS need not agree, and after independent QC they
+# normally do not: each object's sketch is trimmed to its own span and its own
+# surviving variants, so one LD reference used for a QTL and a GWAS collection
+# yields two differently trimmed sketches. Downstream LD lookups are all
+# id-matched, so a partial overlap is reported (once per session) rather than
+# refused; no overlap at all is an error. Shared by causalInferencePipeline,
+# colocPipeline, qtlEnrichmentPipeline, ctwasPipeline, and colocboostPipeline.
 #
 # NULL handling:
 #   nullPolicy = "qtl-required" (default): a NULL qtlLd skips the check; a
@@ -1628,7 +1633,17 @@ loadLdFromGenotype <- function(
     FALSE
 }
 
-# Both must be GenotypeHandles with matching panel size.
+# Both must be genotype panels. Panel SIZE is deliberately NOT compared: two
+# objects that share one LD sketch stop carrying identical variant sets as soon
+# as they are QC'd separately, because each object's sketch is trimmed to its
+# own position span at load (`.subsetSketchToRange`) and to its own surviving
+# variants at the end of QC (`.subsetSketchToIds`). A QTL and a GWAS collection
+# put through `summaryStatsQc` independently therefore diverge by construction,
+# which is the normal case rather than a mistake. Every downstream LD lookup
+# matches by variant id (`.ldFromSketch`), never by position in the panel, so
+# divergent variant sets cost only the variants one side lacks. What would not
+# be harmless -- two genuinely different reference panels -- is what
+# `.ldSketchCheckContent` catches, by sample set and by allele coding.
 .ldSketchCheckShape <- function(qtlLd, gwasLd, pipelineName, between) {
     if (!.ldIsPanel(qtlLd) || !.ldIsPanel(gwasLd)) {
         msg <- glue(
@@ -1637,40 +1652,137 @@ loadLdFromGenotype <- function(
         )
         abort(msg)
     }
-    nQ <- length(.ldSketchRanges(qtlLd))
-    nG <- length(.ldSketchRanges(gwasLd))
-    if (nQ != nG) {
+}
+
+# The (chrom, position, allele-pair) key the two panels are compared on. The
+# allele pair is order-insensitive, so an A1/A2 swap keys to the SAME variant
+# and is then reported as an allele-coding difference rather than looking like
+# two unrelated variants that happen not to overlap. A panel carrying no allele
+# columns falls back to chrom:position, which is all it can be keyed on.
+# @noRd
+.ldSketchVariantKeys <- function(x) {
+    gr <- .ldSketchRanges(x)
+    mc <- S4Vectors::mcols(gr)
+    stem <- str_c(
+        canonChrom(as.character(GenomicRanges::seqnames(gr))),
+        ":",
+        as.character(GenomicRanges::start(gr))
+    )
+    if (is.null(mc$A1) || is.null(mc$A2)) {
+        return(stem)
+    }
+    # str_c propagates NA, which would make a variant with an unknown allele
+    # match nothing on either side; an empty field keeps it comparable.
+    a1 <- as.character(mc$A1)
+    a2 <- as.character(mc$A2)
+    a1 <- if_else(is.na(a1), "", a1)
+    a2 <- if_else(is.na(a2), "", a2)
+    str_c(stem, ":", if_else(a1 < a2, a1, a2), ":", if_else(a1 < a2, a2, a1))
+}
+
+# The panel's A1 (effect / counted) allele, which fixes the sign of every
+# correlation the variant takes part in.
+# @noRd
+.ldSketchA1 <- function(x) {
+    as.character(S4Vectors::mcols(.ldSketchRanges(x))$A1)
+}
+
+# A shared variant whose A1/A2 are swapped between the panels is the same
+# variant coded in opposite directions, so the two panels' LD matrices disagree
+# in sign wherever it appears. That is a different reference, not a trimming
+# difference, and it is an error.
+# @noRd
+.ldSketchCheckAlleleCoding <- function(
+    qtlLd,
+    gwasLd,
+    qIdx,
+    gIdx,
+    pipelineName,
+    between
+) {
+    qA1 <- .ldSketchA1(qtlLd)[qIdx]
+    gA1 <- .ldSketchA1(gwasLd)[gIdx]
+    if (length(qA1) == 0L || length(gA1) == 0L) {
+        return(invisible(NULL))
+    }
+    nSwapped <- sum(qA1 != gA1, na.rm = TRUE)
+    if (nSwapped > 0L) {
         msg <- glue(
-            "{pipelineName}: ldSketch panels differ in size ({nQ} vs ",
-            "{nG} variants){between}; the two ldSketch panels ",
-            "must match exactly."
+            "{pipelineName}: {nSwapped} variant(s) present in both ldSketch ",
+            "panels carry swapped A1/A2 alleles{between}, so the two panels ",
+            "code their LD in opposite directions; use the same ldSketch on ",
+            "both."
         )
         abort(msg)
     }
+    invisible(NULL)
 }
 
-# Panels must agree on CHR/BP/A1/A2 columns and on the sample set.
-.ldSketchCheckContent <- function(qtlLd, gwasLd, pipelineName, between) {
-    qGr <- .ldSketchRanges(qtlLd)
-    gGr <- .ldSketchRanges(gwasLd)
-    if (!identical(.ldSketchChrom(qtlLd), .ldSketchChrom(gwasLd))) {
+# A partial overlap is the expected result of QC-ing the two sides separately,
+# so it is reported once per session rather than on every call: ctwasPipeline
+# and colocboostPipeline run this check once per region / per bundle.
+# @noRd
+.ldSketchReportOverlap <- function(nQ, nG, nShared, pipelineName, between) {
+    if (nShared == nQ && nShared == nG) {
+        return(invisible(NULL))
+    }
+    msg <- glue(
+        "{pipelineName}: the two ldSketch panels share {nShared} variant(s) ",
+        "of {nQ} (QTL side) and {nG} (GWAS side){between}. LD is looked up ",
+        "per variant, so only the shared ones contribute. Differing variant ",
+        "sets are expected when the two sides were QC'd separately."
+    )
+    warn(
+        msg,
+        .frequency = "once",
+        .frequency_id = str_c("pecotmrLdSketchOverlap-", pipelineName)
+    )
+    invisible(NULL)
+}
+
+# Compare the two panels variant by variant. A partial overlap passes; what
+# fails is no overlap at all (two unrelated panels) or a shared variant whose
+# alleles are swapped.
+# @noRd
+.ldSketchCheckOverlap <- function(qtlLd, gwasLd, pipelineName, between) {
+    qKey <- .ldSketchVariantKeys(qtlLd)
+    gKey <- .ldSketchVariantKeys(gwasLd)
+    # An emptied panel -- what a zero-variant object carries, since it
+    # references no LD -- contradicts nothing, so there is nothing to compare.
+    if (length(qKey) == 0L || length(gKey) == 0L) {
+        return(invisible(NULL))
+    }
+    idx <- match(qKey, gKey)
+    shared <- which(!is.na(idx))
+    if (length(shared) == 0L) {
         msg <- glue(
-            "{pipelineName}: ldSketch panels differ in column CHR",
-            "{between}; use the same ldSketch on both."
+            "{pipelineName}: the two ldSketch panels share no variant",
+            "{between} ({length(qKey)} vs {length(gKey)} variants); they ",
+            "describe different LD references."
         )
         abort(msg)
     }
-    qCol <- .ldRangesColumns(qGr)
-    gCol <- .ldRangesColumns(gGr)
-    for (col in c("BP", "A1", "A2")) {
-        if (!identical(qCol[[col]], gCol[[col]])) {
-            msg <- glue(
-                "{pipelineName}: ldSketch panels differ in column ",
-                "{col}{between}; use the same ldSketch on both."
-            )
-            abort(msg)
-        }
-    }
+    .ldSketchCheckAlleleCoding(
+        qtlLd,
+        gwasLd,
+        shared,
+        idx[shared],
+        pipelineName,
+        between
+    )
+    .ldSketchReportOverlap(
+        length(qKey),
+        length(gKey),
+        length(shared),
+        pipelineName,
+        between
+    )
+}
+
+# Panels must describe the same reference panel: the same samples, and the same
+# allele coding on the variants they share. Their variant SETS need not agree
+# -- see `.ldSketchCheckShape` for why.
+.ldSketchCheckContent <- function(qtlLd, gwasLd, pipelineName, between) {
     qIds <- .ldSketchSampleIds(qtlLd)
     gIds <- .ldSketchSampleIds(gwasLd)
     if (!identical(qIds, gIds)) {
@@ -1680,6 +1792,7 @@ loadLdFromGenotype <- function(
         )
         abort(msg)
     }
+    .ldSketchCheckOverlap(qtlLd, gwasLd, pipelineName, between)
 }
 
 .requireMatchingLdSketches <- function(

R/qtlEnrichmentPipeline.R

@@ -31,10 +31,13 @@
 #'     collection, which has no context axis.
 #' }
 #'
-#' @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.
+#' @section LD-sketch compatibility 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 come from the same reference panel --- same
+#'   samples, same allele orientation on the variants they share --- but need
+#'   NOT carry the same variants: QC-ing the two sides separately trims each
+#'   sketch to its own surviving variants, and a partial overlap only warns. A
+#'   \code{NULL} on either side (an individual-level fit) skips the check.
 #'
 #' @param gwasFineMappingResult The outcome side; see above.
 #' @param qtlFineMappingResult The annotation side; see above.