From 6d48d6137c15654c3289add19a8a0960130aa833 Mon Sep 17 00:00:00 2001 From: Daniel Nachun Date: Fri, 4 Sep 2026 11:09:46 -0700 Subject: [PATCH] fix variant normalization --- R/sumstatsQc.R | 17 ++++++++++++++--- tests/testthat/test_sumstatsQc.R | 25 +++++++++++++++++++++++++ 2 files changed, 39 insertions(+), 3 deletions(-) diff --git a/R/sumstatsQc.R b/R/sumstatsQc.R index 1e907463..3ad17b83 100644 --- a/R/sumstatsQc.R +++ b/R/sumstatsQc.R @@ -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) } diff --git a/tests/testthat/test_sumstatsQc.R b/tests/testthat/test_sumstatsQc.R index f06345e1..51bb06ed 100644 --- a/tests/testthat/test_sumstatsQc.R +++ b/tests/testthat/test_sumstatsQc.R @@ -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