diff --git a/R/causalInferencePipeline.R b/R/causalInferencePipeline.R index 5cbb02af..c15228ca 100644 --- a/R/causalInferencePipeline.R +++ b/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( diff --git a/R/colocPipeline.R b/R/colocPipeline.R index cf8565f5..8a21f9ec 100644 --- a/R/colocPipeline.R +++ b/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 diff --git a/R/colocboostPipeline.R b/R/colocboostPipeline.R index b7965153..eb9c599e 100644 --- a/R/colocboostPipeline.R +++ b/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 diff --git a/R/ctwasPipeline.R b/R/ctwasPipeline.R index 1c46ff42..b0d2e3c1 100644 --- a/R/ctwasPipeline.R +++ b/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") diff --git a/R/ld.R b/R/ld.R index aa868842..6f294c6a 100644 --- a/R/ld.R +++ b/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( diff --git a/R/qtlEnrichmentPipeline.R b/R/qtlEnrichmentPipeline.R index bcaf5a7a..9519f616 100644 --- a/R/qtlEnrichmentPipeline.R +++ b/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. diff --git a/man/causalInferencePipeline.Rd b/man/causalInferencePipeline.Rd index 2a55a945..e192f20b 100644 --- a/man/causalInferencePipeline.Rd +++ b/man/causalInferencePipeline.Rd @@ -117,12 +117,17 @@ Per-region pipeline that pairs QTL-derived weight vectors MR computed from \code{fineMappingResult}. } } -\section{LD-sketch identity check}{ +\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}{ diff --git a/man/colocPipeline.Rd b/man/colocPipeline.Rd index b1a5e1b1..7ffb7062 100644 --- a/man/colocPipeline.Rd +++ b/man/colocPipeline.Rd @@ -183,13 +183,19 @@ This pipeline preserves all three. } } -\section{LD-sketch identity check}{ +\section{LD-sketch compatibility 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. + \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. } \examples{ diff --git a/man/ctwasPipeline.Rd b/man/ctwasPipeline.Rd index 0a0eec01..bacb53fb 100644 --- a/man/ctwasPipeline.Rd +++ b/man/ctwasPipeline.Rd @@ -164,10 +164,14 @@ Pipeline that hands a per-block set of 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. } \examples{ diff --git a/man/qtlEnrichmentPipeline.Rd b/man/qtlEnrichmentPipeline.Rd index e132d638..c257aa45 100644 --- a/man/qtlEnrichmentPipeline.Rd +++ b/man/qtlEnrichmentPipeline.Rd @@ -100,11 +100,14 @@ genome of outcome PIPs and the full collection of annotation fits at once. } } -\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. } \examples{ diff --git a/tests/testthat/test_causalInferencePipeline.R b/tests/testthat/test_causalInferencePipeline.R index 05044a2d..6fb43e6f 100644 --- a/tests/testthat/test_causalInferencePipeline.R +++ b/tests/testthat/test_causalInferencePipeline.R @@ -224,9 +224,9 @@ test_that(".cipRequireMatchingLdSketches: NULL twas-side ldSketch is allowed", { expect_s4_class(out, "GRanges") }) -test_that(".cipRequireMatchingLdSketches: panel size mismatch errors", { - bigSketch <- .cip_makeHandle(snp_n = 7L) - twBig <- TwasWeights( +test_that(".cipRequireMatchingLdSketches: sample set mismatch errors", { + otherPanel <- .cip_makeHandle(sample_prefix = "other") + twOther <- TwasWeights( study = "Q1", context = "c1", trait = "t1", @@ -235,14 +235,28 @@ test_that(".cipRequireMatchingLdSketches: panel size mismatch errors", { variantIds = sprintf("chr1:%d:A:G", 100L * (1:5)), weights = rep(0.1, 5) )), - ldSketch = bigSketch + ldSketch = otherPanel ) expect_error( causalInferencePipeline( gwasSumStats = .cip_makeGwasSumstats(), - twasWeights = twBig + twasWeights = twOther ), - "differ in size" + "different sample sets" + ) +}) + +test_that(".cipRequireMatchingLdSketches: differently trimmed panels pass", { + # Independent QC of the QTL and GWAS sides leaves one LD reference as two + # overlapping-but-unequal panels, which the check reports rather than + # refuses. + expect_warning( + expect_null(pecotmr:::.cipRequireMatchingLdSketches( + .cip_makeHandle(snp_n = 6L), + .cip_makeHandle(snp_n = 7L), + label = "twasWeights" + )), + "share 6 variant" ) }) diff --git a/tests/testthat/test_colocPipeline.R b/tests/testthat/test_colocPipeline.R index dac04fda..ee2c318f 100644 --- a/tests/testthat/test_colocPipeline.R +++ b/tests/testthat/test_colocPipeline.R @@ -241,18 +241,30 @@ test_that(".colocRequireMatchingLdSketches: non-NULL qtl + NULL gwas errors", { ) }) -test_that(".colocRequireMatchingLdSketches: panel size mismatch errors", { +test_that(".colocRequireMatchingLdSketches: sample set mismatch errors", { qfmr <- .cp_makeQtlFmr() - bigSketch <- .cp_makeHandle(snp_n = 7L) + otherPanel <- .cp_makeHandle(sample_prefix = "other") gfmr <- GwasFineMappingResult( study = "G1", method = "susie", entry = list(.cp_makeFmEntry()), - ldSketch = bigSketch + ldSketch = otherPanel ) expect_error( colocPipeline(qtlFineMappingResult = qfmr, gwasInput = gfmr), - "ldSketch panels differ in size" + "different sample sets" + ) +}) + +test_that(".colocRequireMatchingLdSketches: differently trimmed panels pass", { + # The two sides QC'd separately keep different subsets of one LD + # reference; the panels overlap, so the check warns rather than aborting. + expect_warning( + expect_null(pecotmr:::.colocRequireMatchingLdSketches( + .cp_makeHandle(snp_n = 6L), + .cp_makeHandle(snp_n = 7L) + )), + "share 6 variant" ) }) diff --git a/tests/testthat/test_colocboostPipeline.R b/tests/testthat/test_colocboostPipeline.R index e7e04aaa..d54815fe 100644 --- a/tests/testthat/test_colocboostPipeline.R +++ b/tests/testthat/test_colocboostPipeline.R @@ -326,13 +326,25 @@ test_that(".cbRequireMatchingLdSketches: NULL sides are allowed", { )) }) -test_that(".cbRequireMatchingLdSketches: variant-count mismatch errors", { - expect_error( - pecotmr:::.cbRequireMatchingLdSketches( +test_that(".cbRequireMatchingLdSketches: accepts differently trimmed panels", { + # Separate QC of the two sides trims one shared LD reference to two + # overlapping-but-unequal variant sets; the check reports, not refuses. + expect_warning( + expect_null(pecotmr:::.cbRequireMatchingLdSketches( .cbp_makeHandle(snp_n = 4L), .cbp_makeHandle(snp_n = 5L) + )), + "share 4 variant" + ) +}) + +test_that(".cbRequireMatchingLdSketches: sample-set mismatch errors", { + expect_error( + pecotmr:::.cbRequireMatchingLdSketches( + .cbp_makeHandle(sample_prefix = "a"), + .cbp_makeHandle(sample_prefix = "b") ), - "differ in size" + "different sample sets" ) }) diff --git a/tests/testthat/test_ctwasPipeline.R b/tests/testthat/test_ctwasPipeline.R index 33ec6f74..d981f6e8 100644 --- a/tests/testthat/test_ctwasPipeline.R +++ b/tests/testthat/test_ctwasPipeline.R @@ -7,8 +7,8 @@ context("ctwasPipeline") # =========================================================================== # 12 = 6 variants x the 2 blocks the default GwasSumStats fixture builds. The -# LD-sketch identity check compares panels between the GWAS and weight sides, -# so every fixture has to draw on the same one. +# LD-sketch compatibility check compares panels between the GWAS and weight +# sides, so every fixture has to draw on the same one. .ctp_makeHandle <- function(snp_n = 12L, n_samples = 30L) { # Use a per-process tempfile so .ctwasLdPanelKey's file.exists check # succeeds against the fixture handle (real LD-sketch payloads exist @@ -427,12 +427,27 @@ test_that(".ctwasRequireMatchingLdSketches: NULL twas-side handle is allowed", { )) }) -test_that(".ctwasRequireMatchingLdSketches: panel-size mismatch errors", { +test_that(".ctwasRequireMatchingLdSketches: differently trimmed panels ok", { + # One LD reference QC'd separately on the two sides leaves each sketch + # trimmed to its own surviving variants, so the panels overlap without + # being identical. That is the normal case, not an error. twLd <- .ctp_makeHandle(snp_n = 5L) gwasLd <- .ctp_makeHandle(snp_n = 6L) + expect_warning( + expect_null(pecotmr:::.ctwasRequireMatchingLdSketches(twLd, gwasLd)), + "share 5 variant" + ) +}) + +test_that(".ctwasRequireMatchingLdSketches: disjoint panels error", { + twLd <- .ctp_makeHandle(snp_n = 3L) + gwasLd <- .ctp_makeHandle(snp_n = 3L) + si <- getSnpInfo(gwasLd) + si$BP <- si$BP + 1e6L + gwasLd@snpInfo <- si expect_error( pecotmr:::.ctwasRequireMatchingLdSketches(twLd, gwasLd), - "ldSketch panels differ in size" + "share no variant" ) }) diff --git a/tests/testthat/test_ld.R b/tests/testthat/test_ld.R index 5490c89b..7084e066 100644 --- a/tests/testthat/test_ld.R +++ b/tests/testthat/test_ld.R @@ -3295,7 +3295,7 @@ test_that(".requireMatchingLdSketches errors when slots are not GenotypeHandle", ) }) -test_that(".requireMatchingLdSketches errors when panels differ in a column", { +test_that(".requireMatchingLdSketches errors on swapped alleles", { skip_if_not_installed("pgenlibr") h <- readGenotypeHandle( file.path(geno_test_data_dir, "test_variants"), @@ -3303,7 +3303,8 @@ test_that(".requireMatchingLdSketches errors when panels differ in a column", { ) si <- getSnpInfo(h) si2 <- si - si2$A1[1] <- if (identical(si2$A1[1], "A")) "C" else "A" # mutate one allele + si2$A1[1] <- si$A2[1] # swap the coding of one shared variant + si2$A2[1] <- si$A1[1] h2 <- new( "GenotypeHandle", path = getPath(h), @@ -3316,7 +3317,107 @@ test_that(".requireMatchingLdSketches errors when panels differ in a column", { ) expect_error( pecotmr:::.requireMatchingLdSketches(h, h2, "testPipeline"), - "differ in column" + "swapped A1/A2" + ) +}) + +test_that(".requireMatchingLdSketches accepts panels trimmed differently", { + # The reported bug: one LD sketch shared by a QtlSumStats and a + # GwasSumStats, each trimmed by its own summaryStatsQc to a different + # surviving variant set. The panels overlap but are not identical. + si <- data.frame( + SNP = c("1:100:A:G", "1:200:C:T", "1:300:G:A", "1:400:T:C"), + CHR = rep("1", 4L), + BP = c(100L, 200L, 300L, 400L), + A1 = c("A", "C", "G", "T"), + A2 = c("G", "T", "A", "C"), + stringsAsFactors = FALSE + ) + mk <- function(rows) { + new( + "GenotypeHandle", + path = "/tmp/x", + format = "gds", + snpInfo = si[rows, , drop = FALSE], + nSamples = 3L, + sampleIds = str_c("s", 1:3), + pgenPtr = NULL, + chromPaths = character(0) + ) + } + expect_warning( + expect_null( + pecotmr:::.requireMatchingLdSketches( + mk(1:3), + mk(2:4), + "trimPipelineA" + ) + ), + "share 2 variant" + ) +}) + +test_that(".requireMatchingLdSketches errors when panels share no variant", { + mkSi <- function(bp) { + data.frame( + SNP = str_c("1:", bp, ":A:G"), + CHR = rep("1", length(bp)), + BP = as.integer(bp), + A1 = rep("A", length(bp)), + A2 = rep("G", length(bp)), + stringsAsFactors = FALSE + ) + } + mk <- function(bp) { + new( + "GenotypeHandle", + path = "/tmp/x", + format = "gds", + snpInfo = mkSi(bp), + nSamples = 3L, + sampleIds = str_c("s", 1:3), + pgenPtr = NULL, + chromPaths = character(0) + ) + } + expect_error( + pecotmr:::.requireMatchingLdSketches( + mk(c(100, 200)), + mk(c(900, 950)), + "disjointPipeline" + ), + "share no variant" + ) +}) + +test_that(".requireMatchingLdSketches errors on a different sample set", { + si <- data.frame( + SNP = "1:100:A:G", + CHR = "1", + BP = 100L, + A1 = "A", + A2 = "G", + stringsAsFactors = FALSE + ) + mk <- function(ids) { + new( + "GenotypeHandle", + path = "/tmp/x", + format = "gds", + snpInfo = si, + nSamples = length(ids), + sampleIds = ids, + pgenPtr = NULL, + chromPaths = character(0) + ) + } + expect_error( + pecotmr:::.requireMatchingLdSketches( + mk(str_c("s", 1:3)), + mk(str_c("t", 1:3)), + "samplePipeline" + ), + "different sample sets" ) }) @@ -4727,6 +4828,8 @@ test_that(".ldSketchNullGuard names the label in its strict error", { }) test_that(".ldSketchCheckContent rejects panels on different chromosomes", { + # Disjoint chromosomes leave the two panels with nothing in common, which + # is a different LD reference rather than a trimming difference. mkSe <- function(chr, pos) { g <- GenomicRanges::GRanges(chr, IRanges::IRanges(pos, width = 1)) S4Vectors::mcols(g)$SNP <- str_c(chr, ":", pos, ":A:G") @@ -4742,7 +4845,7 @@ test_that(".ldSketchCheckContent rejects panels on different chromosomes", { "myPipe", " between X and Y" ), - "differ in column CHR between X and Y" + "share no variant between X and Y" ) })