fix allele filtering

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),

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)

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)
     }

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

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)

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", {