Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
17 changes: 14 additions & 3 deletions R/sumstatsQc.R
Original file line number Diff line number Diff line change
Expand Up @@ -4412,10 +4412,21 @@ krigingOutlierQc <- function(
if (length(ids) == 0L) {
return(ldSketch)
}
keep <- is_in(
normalizeVariantId(.ldSketchMatchIds(ldSketch)),
normalizeVariantId(unique(ids))
# Match on the (chrom, pos, allele) tuple rather than the raw id string.
# The entry ids have just been canonicalized to chr:pos:A2:A1 by the QC
# above, while panel ids are passed through verbatim by
# `.repairVariantIds()` whenever their allele fields are already valid DNA
# -- so a panel keyed chr:pos:A1:A2 (which is how a PLINK .bim writes ids;
# a .pvar writes REF:ALT and is therefore already canonical) matches
# nothing under string equality and the panel is silently emptied. This is
# the matcher `.ldFromSketchMatch()` and every other LD lookup already use.
panelIds <- .ldSketchMatchIds(ldSketch)
matched <- matchVariants(
panelIds,
unique(ids),
removeStrandAmbiguous = FALSE
)
keep <- is_in(seq_along(panelIds), matched$idxA)
.ldSketchSubset(ldSketch, keep)
}

Expand Down
25 changes: 25 additions & 0 deletions tests/testthat/test_sumstatsQc.R
Original file line number Diff line number Diff line change
Expand Up @@ -3324,6 +3324,31 @@ test_that("a tag-named panel entry survives the end-to-end sketch subset", {
)
})

test_that("the sketch survives a panel keyed chr:pos:A1:A2 (PLINK .bim order)", {
# Pre-fix regression: harmonization canonicalizes the entry ids to
# chr:pos:A2:A1, while panel ids are passed through verbatim whenever their
# allele fields are already valid DNA. A panel keyed the other way round --
# which is how a PLINK .bim writes ids, whereas a .pvar writes REF:ALT and
# is already canonical -- matched nothing under string equality, so
# .subsetSketchToIds() silently emptied the sketch and RSS fine-mapping /
# TWAS weights aborted later with "not present in the LD sketch panel".
h <- .ssQ_makeHandleVid()
si <- getSnpInfo(h)
h@snpInfo$SNP <- paste0("chr1:", si$BP, ":", si$A1, ":", si$A2)
ss <- GwasSumStats(
study = "g1",
entry = list(.ssQ_makeEntryGr(
snp_ids = paste0("chr1:", c(100L, 200L, 300L, 400L), ":G:A")
)),
genome = "hg19",
ldSketch = h
)
out <- summaryStatsQc(ss, pipCutoffToSkip = 0, nCutoff = 0)
kept <- pecotmr:::.ldSketchMatchIds(getLdSketch(out))
expect_length(kept, 4L)
expect_setequal(pecotmr:::parseVariantId(kept)$pos, c(100L, 200L, 300L, 400L))
})

test_that("summaryStatsQc: slalom z-mismatch resolves sign-flipped variants against the panel", {
# Pre-fix regression: a sign-flipped variant kept its input-orientation SNP,
# which is absent from the panel, so .applyLdMismatchQcToEntry errored with
Expand Down