From d88d5f0760f9584058cd512fe1ef3ffce311711f Mon Sep 17 00:00:00 2001 From: Daniel Nachun Date: Mon, 31 Aug 2026 13:03:29 -0700 Subject: [PATCH] fix allele filtering --- R/ctwasPipeline.R | 9 +++- R/ld.R | 27 ++++++++++-- R/sumstatsQc.R | 40 ++++------------- R/variantId.R | 66 ++++++++++++++++++---------- tests/testthat/test_ctwasPipeline.R | 19 ++++++++ tests/testthat/test_sumstatsQc.R | 68 ++++++++++++++++++----------- tests/testthat/test_variantId.R | 67 ++++++++++++++++++++++++++++ 7 files changed, 213 insertions(+), 83 deletions(-) diff --git a/R/ctwasPipeline.R b/R/ctwasPipeline.R index 9e59d934..ae1314a3 100644 --- a/R/ctwasPipeline.R +++ b/R/ctwasPipeline.R @@ -2170,12 +2170,19 @@ asCtwasResult <- function(finemapResult, keepSnps = FALSE) { .ctwasSnpInfoForBlock <- function(gwasLd) { gr <- .ldSketchRanges(gwasLd) mc <- S4Vectors::mcols(gr) + # `.ldSketchMatchIds()`, not the raw SNP label: this id becomes the R + # dimnames, the `variance` names, `panelSnps` and the weight-harmonization + # reference id all at once, and every one of those is compared against a + # harmonized (reference-allele) id from the GWAS or the weights. A panel + # entry spelling a tag where an allele belongs would fail all four + # comparisons and silently drop the variant from the analysis. + # # Base data.frame (not tibble): ctwas indexes snp_map positionally # (df[, "pos"] -> vector for region-bound min/max); a tibble column is a # 1-col list and errors min(). data.frame( chrom = as.integer(.ldSketchChrom(gwasLd)), - id = as.character(mc$SNP), + id = .ldSketchMatchIds(gwasLd), pos = as.integer(GenomicRanges::start(gr)), alt = as.character(mc$A1), ref = as.character(mc$A2), diff --git a/R/ld.R b/R/ld.R index c2e54952..e3cbdad3 100644 --- a/R/ld.R +++ b/R/ld.R @@ -917,13 +917,34 @@ loadLdFromGenotype <- function( } # The panel's variant ids in the panel's OWN labelling. This is the mcols -# column, not `names()`: the names are normalized ids, while matching against -# summary statistics has to see the panel's raw labels. +# column, not `names()`: the names are normalized ids, while anything keyed on +# the panel FILE -- the `.afreq` sidecar, the ctwas snp_map -- has to see the +# panel's raw labels. # @noRd .ldSketchVariantIds <- function(x) { as.character(S4Vectors::mcols(.ldSketchRanges(x))$SNP) } +# The panel's variant ids as MATCHING keys: the raw labels, with any id whose +# allele slot holds a tag rather than an allele (chr21:13988152:INS:T) re- +# rendered from the panel's own A1 / A2 columns, which are correct. Use this +# wherever a panel id is matched against summary statistics or another id set, +# and `.ldSketchVariantIds()` only where the file's literal label is wanted. +# +# The repair is what lets an untidily-named panel entry be MATCHED rather than +# discarded: harmonization re-keys a surviving variant to the reference-allele +# id, so the panel side has to speak the same form or every later id lookup +# misses the variant it just harmonized. +# @noRd +.ldSketchMatchIds <- function(x) { + mc <- S4Vectors::mcols(.ldSketchRanges(x)) + ids <- as.character(mc$SNP) + if (is.null(mc$A1) || is.null(mc$A2)) { + return(ids) + } + .repairVariantIds(ids, as.character(mc$A2), as.character(mc$A1)) +} + # Per-variant chromosome, canonical (no chr prefix). # @noRd .ldSketchChrom <- function(x) { @@ -998,7 +1019,7 @@ loadLdFromGenotype <- function( # rsID panels; the caller's original ids and order are preserved. m <- matchVariants( variantIds, - .ldSketchVariantIds(ldSketch), + .ldSketchMatchIds(ldSketch), removeStrandAmbiguous = FALSE ) nMissing <- length(variantIds) - length(m$idxA) diff --git a/R/sumstatsQc.R b/R/sumstatsQc.R index f26a5812..5bd26c6c 100644 --- a/R/sumstatsQc.R +++ b/R/sumstatsQc.R @@ -1722,12 +1722,7 @@ autoDecision <- function(df, highCorrCols) { # correspondence between `knownZscores$z` and the LD rows indexed by `knowns`. # @noRd .raissUnsafeToImpute <- function(refPanelIds, knownIds) { - # unusable-id: the id does not encode its alleles (e.g. an R5 INS/DEL tag). - # It is dropped from harmonization for the same reason -- the sumstats cannot - # say which orientation was measured -- so imputing it would put back a - # variant the drop exists to exclude. - !.ldSketchIdUsable(refPanelIds) | - .raissFlipPairMask(refPanelIds) | + .raissFlipPairMask(refPanelIds) | .raissFlipOfKnownMask(refPanelIds, knownIds) } @@ -2853,6 +2848,9 @@ krigingOutlierQc <- function( # Build a refVariants data.frame (chrom, pos, A1, A2, variant_id) from the # panel's own variant ranges so harmonizeAlleles can join by (chrom, pos). +# The id comes from `.ldSketchMatchIds()`, not the raw SNP label: a panel entry +# that spells a tag where an allele belongs still harmonizes off its A1/A2 +# columns, and the repaired id is what the LD lookups downstream can find. .refVariantsFromSketch <- function(ldSketch) { gr <- .ldSketchRanges(ldSketch) mc <- S4Vectors::mcols(gr) @@ -2861,7 +2859,7 @@ krigingOutlierQc <- function( pos = as.integer(GenomicRanges::start(gr)), A1 = as.character(mc$A1), A2 = as.character(mc$A2), - variant_id = as.character(mc$SNP), + variant_id = .ldSketchMatchIds(ldSketch), stringsAsFactors = FALSE ) } @@ -3156,26 +3154,6 @@ krigingOutlierQc <- function( # harmonizeAlleles against the ldSketch's variant info. Threads the # variant-level filters (indels, strand-ambiguous, duplicates) through # so the LD-panel-anchored pass handles them in a single sweep. -# The harmonization reference set: the panel's variants, minus the entries -# whose id does not encode its alleles (e.g. R5 INS/DEL tags). Their A1/A2 -# columns may be correct, but the summary statistics carry no way to say which -# orientation was measured at that position, so a GWAS variant there must be -# dropped rather than aligned by a guess that could sign-flip a real signal -- -# the same "exclude what you can't interpret" rule the RAISS flip-pair / -# flip-of-known masks apply. Removing the reference row leaves the GWAS variant -# unmatched, so removeUnmatched removes it. -# -# Filtered HERE rather than inside .refVariantsFromSketch because the RAISS -# path indexes into that frame and would misalign against its dosage if it -# were filtered globally. -# @noRd -.sketchRefVariants <- function(ldSketch) { - filter( - .refVariantsFromSketch(ldSketch), - .ldSketchIdUsable(.data$variant_id) - ) -} - .matchAgainstSketch <- function( df, ldSketch, @@ -3184,7 +3162,7 @@ krigingOutlierQc <- function( removeStrandAmbiguous = TRUE, removeDups = TRUE ) { - refVariants <- .sketchRefVariants(ldSketch) + refVariants <- .refVariantsFromSketch(ldSketch) flipCandidates <- c("Z", "BETA") colToFlip <- intersect(flipCandidates, colnames(df)) if (length(colToFlip) == 0L) { @@ -3929,7 +3907,7 @@ krigingOutlierQc <- function( # TRUE did. dosage <- .ldSketchDosage(ldSketch, windowIdx, meanImpute = FALSE) colnames(dosage) <- normalizeVariantId( - .ldSketchVariantIds(ldSketch)[windowIdx] + .ldSketchMatchIds(ldSketch)[windowIdx] ) dosage <- dosage[, refPanel$variant_id, drop = FALSE] keep <- .qcRaissTargetMask(refPanel, knownZ, dosage, opts) @@ -4642,7 +4620,7 @@ krigingOutlierQc <- function( return(ldSketch) } keep <- is_in( - normalizeVariantId(.ldSketchVariantIds(ldSketch)), + normalizeVariantId(.ldSketchMatchIds(ldSketch)), normalizeVariantId(unique(ids)) ) .ldSketchSubset(ldSketch, keep) @@ -4792,7 +4770,7 @@ krigingOutlierQc <- function( if (is.null(ldSketch) || is.null(cutoffs)) { return(ldSketch) } - ids <- .ldSketchVariantIds(ldSketch) + ids <- .ldSketchMatchIds(ldSketch) if (length(ids) == 0L) { return(ldSketch) } diff --git a/R/variantId.R b/R/variantId.R index 84e07eab..56de6bc2 100644 --- a/R/variantId.R +++ b/R/variantId.R @@ -98,33 +98,55 @@ isSnpAlleles <- function(a1, a2) { replace_na(isSnp, FALSE) } -#' Test whether variant ids encode usable DNA alleles (not a tag like INS/DEL). +#' Re-render variant ids whose allele slot holds a tag instead of an allele. +#' +#' Some LD panels (e.g. ADSP R5) name a variant with a tag where an allele +#' belongs -- \code{chr21:13988152:INS:T} whose real alleles are +#' \code{A}/\code{AT} -- so the id string stops encoding its alleles even +#' though the panel's own A1/A2 columns are correct. The id-string matchers +#' (\code{.ldFromSketchMatch}, \code{.subsetSketchToIds}) key on the id and +#' cannot interpret such a record. The allele COLUMNS are authoritative, so the +#' id is rebuilt from them; the variant is never dropped for how its id happens +#' to spell an allele. +#' +#' Only an id that parses to \code{chr:pos:allele:allele} but whose allele +#' slots are not pure DNA is rewritten, and only where both supplied alleles are +#' present. An id that does not parse to the allele form -- an rsID, say -- is a +#' different, valid convention that \code{matchVariants} handles by exact +#' string, so it is left alone; so is a clean allele-encoded id, which keeps the +#' panel's own spelling. The DNA test is case-insensitive, so a lower-case id +#' (\code{chr1:100:a:g}) is clean, not a tag. Rewritten ids follow the +#' convention (chr prefix, field separator) detected on \code{ids}, so a +#' repaired id looks like its neighbours. #' -#' Some LD panels (e.g. ADSP R5) name a variant with a tag where an allele belongs -#' -- \code{chr21:13988152:INS:T} whose real alleles are \code{A}/\code{AT} -- so the -#' id no longer encodes its alleles, even though the panel's A1/A2 columns are correct. -#' The id-string matchers (\code{.ldFromSketchMatch}, \code{.subsetSketchToIds}) key on -#' the id and cannot interpret such a record. Returns TRUE only where both parsed -#' alleles are pure DNA; tests for DNA rather than the literal \code{INS}/\code{DEL}, -#' so any tagging convention is caught. A legitimate multi-base indel allele (e.g. -#' \code{TAATGG}) stays usable. -#' An id that does not parse to \code{chr:pos:allele:allele} at all -- e.g. an -#' rsID -- is a different, valid convention (its A1/A2 columns are clean and -#' \code{matchVariants} handles it by string), so it is NOT flagged; only an id -#' that parses to the allele form but whose allele slot holds a non-DNA tag is -#' unusable. #' @param ids Character vector of variant ids (\code{chr:pos:A2:A1}). -#' @return Logical vector, TRUE unless the id parses to a non-DNA allele. +#' @param A2,A1 Character vectors of authoritative reference / alternate +#' alleles, parallel to \code{ids}. +#' @return \code{ids} with the tag-allele entries re-rendered from +#' \code{A2}/\code{A1}. #' @noRd -.ldSketchIdUsable <- function(ids) { +.repairVariantIds <- function(ids, A2, A1) { + ids <- as.character(ids) + if (length(A2) != length(ids) || length(A1) != length(ids)) { + return(ids) + } p <- parseVariantId(ids) - parsed <- !is.na(p$A1) & !is.na(p$A2) - dnaOk <- replace_na( - str_detect(p$A1, "^[ACGT]+$") & str_detect(p$A2, "^[ACGT]+$"), - FALSE + dna <- regex("^[ACGT]+$", ignore_case = TRUE) + tagged <- !is.na(p$A1) & + !is.na(p$A2) & + !replace_na(str_detect(p$A1, dna) & str_detect(p$A2, dna), FALSE) + repair <- tagged & !is.na(A2) & !is.na(A1) & A2 != "" & A1 != "" + if (!any(repair)) { + return(ids) + } + ids[repair] <- formatVariantId( + p$chrom[repair], + p$pos[repair], + A2[repair], + A1[repair], + convention = attr(p, "convention") ) - # usable = not an allele-encoded id (rsID etc.) OR parsed and pure DNA. - !parsed | dnaOk + ids } # Backwards-compat alias diff --git a/tests/testthat/test_ctwasPipeline.R b/tests/testthat/test_ctwasPipeline.R index a7fc24be..e7d4ebff 100644 --- a/tests/testthat/test_ctwasPipeline.R +++ b/tests/testthat/test_ctwasPipeline.R @@ -526,6 +526,25 @@ test_that(".ctwasSnpInfoForBlock: returns ctwas-required columns chrom/id/pos/al expect_setequal(colnames(df), c("chrom", "id", "pos", "alt", "ref")) }) +test_that(".ctwasSnpInfoForBlock: repairs a tag-allele panel id from A1/A2", { + # The id becomes the R dimnames and `panelSnps`, both compared against + # harmonized (reference-allele) ids, so a tag id would drop the variant. + h <- .ctp_makeHandle(snp_n = 3L) + h@snpInfo$SNP[2L] <- "chr1:200:INS:A" + df <- pecotmr:::.ctwasSnpInfoForBlock(h) + expect_equal(df$id[2L], "chr1:200:G:A") + expect_equal(df$ref[2L], "G") + expect_equal(df$alt[2L], "A") + # the LD matrix and variance vector are keyed the same way + local_mocked_bindings( + extractBlockGenotypes = .ctp_mockExtractor(), + .package = "pecotmr" + ) + panel <- pecotmr:::.ctwasComputeFullPanelLd(h) + expect_true("chr1:200:G:A" %in% rownames(panel$R)) + expect_true("chr1:200:G:A" %in% names(panel$variance)) +}) + test_that(".ctwasLdPanelKey: returns the on-disk path for an existing GDS sketch", { handle <- .ctp_makeHandle() key <- pecotmr:::.ctwasLdPanelKey(handle) diff --git a/tests/testthat/test_sumstatsQc.R b/tests/testthat/test_sumstatsQc.R index 7833df40..85f5bea0 100644 --- a/tests/testthat/test_sumstatsQc.R +++ b/tests/testthat/test_sumstatsQc.R @@ -3250,30 +3250,27 @@ test_that("summaryStatsQc: harmonization re-keys SNP to the panel id and sign-fl }) # --------------------------------------------------------------------------- -# Uninterpretable panel ids (e.g. ADSP R5 INS/DEL tags): drop, don't guess +# Tag-allele panel ids (e.g. ADSP R5 INS/DEL): repair from the allele columns, +# never drop the variant # --------------------------------------------------------------------------- -test_that(".ldSketchIdUsable flags tag-allele ids but not SNPs, real indels, or rsIDs", { - ids <- c( - "chr21:13988031:T:C", - "chr1:788757:T:TAATGG", - "chr21:16298:C:T", - "chr21:13988152:INS:T", - "chr21:13988153:DEL:T", - "rs12345", - NA - ) - expect_equal( - pecotmr:::.ldSketchIdUsable(ids), - c(TRUE, TRUE, TRUE, FALSE, FALSE, TRUE, TRUE) - ) -}) - -test_that(".raissUnsafeToImpute now excludes tag-allele ids (third clause)", { - # A tag id has no flip twin, so it slipped through before; it must now be - # unsafe, while a normal SNP stays safe and the flip-pair behaviour is intact. - expect_true(pecotmr:::.raissUnsafeToImpute( - "chr21:13988152:INS:T", +test_that(".ldSketchMatchIds repairs tag ids and leaves clean labels alone", { + h <- .ssQ_makeHandleVid() # ids chr1:100:G:A .. chr1:800:G:A, A1=A A2=G + h@snpInfo$SNP[h@snpInfo$BP == 200L] <- "chr1:200:INS:A" + ids <- pecotmr:::.ldSketchMatchIds(h) + # the tag id is re-rendered from the panel's own A2/A1 columns ... + expect_equal(ids[2L], "chr1:200:G:A") + # ... and every other label is untouched. + expect_equal(ids[-2L], pecotmr:::.ldSketchVariantIds(h)[-2L]) + # the raw accessor still answers the panel FILE's label (.afreq keys on it) + expect_equal(pecotmr:::.ldSketchVariantIds(h)[2L], "chr1:200:INS:A") +}) + +test_that(".raissUnsafeToImpute does not refuse a tag-allele id", { + # A tag id is repaired at the panel boundary, so RAISS sees a normal indel + # id; only the flip-pair / flip-of-known ambiguities are unsafe. + expect_false(pecotmr:::.raissUnsafeToImpute( + "chr21:13988152:A:AT", character(0) )) expect_false(pecotmr:::.raissUnsafeToImpute( @@ -3285,7 +3282,7 @@ test_that(".raissUnsafeToImpute now excludes tag-allele ids (third clause)", { expect_true(all(pecotmr:::.raissUnsafeToImpute(fp, character(0)))) }) -test_that("harmonization drops a variant on a tag-named panel position, keeps normal ones", { +test_that("harmonization keeps a variant on a tag-named panel position", { skip_if_not_installed("pgenlibr") h <- .ssQ_makeHandleVid() # ids chr1:100:G:A .. chr1:800:G:A, A1=A A2=G df <- tibble( @@ -3303,9 +3300,28 @@ test_that("harmonization drops a variant on a tag-named panel position, keeps no # tag the pos-200 panel id (A1/A2 columns stay the real alleles) h@snpInfo$SNP[h@snpInfo$BP == 200L] <- "chr1:200:INS:A" outT <- pecotmr:::.matchAgainstSketch(df, h, matchMinProp = 0) - # pos-200 is dropped (its only reference entry is uninterpretable); pos-100 kept - expect_true(100L %in% pecotmr:::parseVariantId(outT$SNP)$pos) - expect_false(200L %in% pecotmr:::parseVariantId(outT$SNP)$pos) + # the tag costs nothing: both variants survive, with unchanged Z + expect_setequal(pecotmr:::parseVariantId(outT$SNP)$pos, c(100L, 200L)) + expect_equal(outT$SNP, out0$SNP) + expect_equal(outT$Z, out0$Z) +}) + +test_that("a tag-named panel entry survives the end-to-end sketch subset", { + h <- .ssQ_makeHandleVid() + h@snpInfo$SNP[h@snpInfo$BP == 200L] <- "chr1:200:INS:A" + ss <- GwasSumStats( + study = "g1", + entry = list(.ssQ_makeEntryGr()), + genome = "hg19", + ldSketch = h + ) + out <- summaryStatsQc(ss, pipCutoffToSkip = 0, nCutoff = 0) + snp <- as.character(S4Vectors::mcols(out[[1L]])$SNP) + expect_true("chr1:200:G:A" %in% snp) + # the retained sketch keeps the panel row the QC'd entry still refers to + expect_true( + "chr1:200:G:A" %in% pecotmr:::.ldSketchMatchIds(getLdSketch(out)) + ) }) test_that("summaryStatsQc: slalom z-mismatch resolves sign-flipped variants against the panel", { diff --git a/tests/testthat/test_variantId.R b/tests/testthat/test_variantId.R index 14994588..0ad545f9 100644 --- a/tests/testthat/test_variantId.R +++ b/tests/testthat/test_variantId.R @@ -115,6 +115,73 @@ test_that("normalizeVariantId: strips build suffix and re-emits canonical", { expect_equal(res, "chr1:100:A:G") }) +# =========================================================================== +# .repairVariantIds +# =========================================================================== + +test_that(".repairVariantIds: rebuilds a tag-allele id from the allele columns", { + ids <- c("chr21:13988152:INS:T", "chr21:13988153:DEL:T") + expect_equal( + pecotmr:::.repairVariantIds(ids, A2 = c("A", "GT"), A1 = c("AT", "G")), + c("chr21:13988152:A:AT", "chr21:13988153:GT:G") + ) +}) + +test_that(".repairVariantIds: leaves clean, lower-case and rsID labels alone", { + ids <- c( + "chr1:100:A:G", # plain SNP + "chr1:788757:T:TAATGG", # real insertion, spelled out + "chr1:200:a:g", # lower case is DNA, not a tag + "chr1:300_A_G", # underscore convention + "rs12345", # does not parse to the allele form + NA + ) + out <- pecotmr:::.repairVariantIds( + ids, + A2 = c("A", "T", "A", "A", "C", "A"), + A1 = c("G", "TAATGG", "G", "G", "T", "G") + ) + expect_equal(out, as.character(ids)) +}) + +test_that(".repairVariantIds: preserves the ids' own naming convention", { + ids <- c("1:100_A_G", "1:200_INS_G") + expect_equal( + pecotmr:::.repairVariantIds(ids, A2 = c("A", "A"), A1 = c("G", "AG")), + c("1:100_A_G", "1:200_A_AG") + ) +}) + +test_that(".repairVariantIds: keeps the id when the alleles cannot replace it", { + ids <- c("chr1:100:INS:T", "chr1:200:INS:T", "chr1:300:INS:T") + expect_equal( + pecotmr:::.repairVariantIds( + ids, + A2 = c(NA, "", "A"), + A1 = c("AT", "AT", "AT") + ), + c("chr1:100:INS:T", "chr1:200:INS:T", "chr1:300:A:AT") + ) + # A mis-shaped allele vector is a no-op rather than a recycling error. + expect_equal(pecotmr:::.repairVariantIds(ids, "A", "AT"), ids) + expect_equal( + pecotmr:::.repairVariantIds(character(0), character(0), character(0)), + character(0) + ) +}) + +test_that(".repairVariantIds: non-DNA alleles on both sides are still kept", { + # Nothing better can be said about these, but the point of the repair is + # that a variant is never DROPPED for how its id spells an allele. + ids <- c("chr1:100:A:*", "chr1:200:N:G", "chr1:300:A:") + out <- pecotmr:::.repairVariantIds( + ids, + A2 = c("A", "N", "A"), + A1 = c("*", "G", "") + ) + expect_equal(out, ids) +}) + # =========================================================================== # parseRegion # ===========================================================================