Tolerate flip-twin collapse in .ldFromSketch instead of erroring

R/fineMappingPipeline.R

@@ -715,6 +715,22 @@ combineFineMappingResults <- function(..., ldSketch = NULL) {
   .ldFromSketch(ldSketch, variantIds, label = ".fmLdFromSketch")
 }
 
+# A flip-twin region makes .ldFromSketch return an LD matrix over the
+# deduplicated panel variants (a variant and its ref/alt-swapped twin collapse to
+# one), fewer than requested, with the survivors in attr(,"keptVariantIds").
+# Realign the association vectors (z / n / af) to those survivors so the LD matrix
+# and z/N stay in lockstep before the RSS fit. A no-op when nothing collapsed.
+# @noRd
+.fmAlignToKeptLd <- function(ldMat, variantIds, z, n, af = NULL) {
+  kept <- attr(ldMat, "keptVariantIds")
+  if (is.null(kept) || length(kept) == length(variantIds))
+    return(list(ldMat = ldMat, variantIds = variantIds, z = z, n = n, af = af))
+  ki <- match(kept, variantIds)
+  list(ldMat = ldMat, variantIds = kept, z = z[ki],
+       n  = if (length(n)  == length(variantIds)) n[ki]  else n,
+       af = if (!is.null(af) && length(af) == length(variantIds)) af[ki] else af)
+}
+
 
 # Wrap one finemapping fit into a FineMappingEntry via the surviving
 # post-processing helpers (postprocessFinemappingFits +
@@ -1999,6 +2015,9 @@ setMethod("fineMappingPipeline", "QtlSumStats",
         afByVar <- if ("MAF" %in% colnames(.qmc))
           setNames(as.numeric(.qmc$MAF), as.character(.qmc$SNP))[variantIds] else NULL
         ldMat <- .fmLdFromSketch(ldSketch, variantIds)
+        .al <- .fmAlignToKeptLd(ldMat, variantIds, z, n, afByVar)
+        ldMat <- .al$ldMat; variantIds <- .al$variantIds
+        z <- .al$z; n <- .al$n; afByVar <- .al$af
         names(z) <- variantIds
 
         ents <- .fmFitRssBlock(
@@ -2182,6 +2201,9 @@ setMethod("fineMappingPipeline", "GwasSumStats",
       if (length(toRun) == 0L) next
 
       ldMat <- .fmLdFromSketch(ldSketch, variantIds)
+      .al <- .fmAlignToKeptLd(ldMat, variantIds, z, n, afByVar)
+      ldMat <- .al$ldMat; variantIds <- .al$variantIds
+      z <- .al$z; n <- .al$n; afByVar <- .al$af
       names(z) <- variantIds
       ents <- .fmFitRssBlock(
         z, ldMat, n, toRun, addSusieInf, coverage, secondaryCoverage,

R/ld.R

@@ -562,10 +562,27 @@ loadLdFromGenotype <- function(genotypePath, region,
   # original ids are kept as the returned labels, in the requested order.
   m <- matchVariants(variantIds, as.character(snpInfo$SNP),
                      removeStrandAmbiguous = FALSE)
-  nMissing <- length(variantIds) - length(m$idxA)
-  if (nMissing > 0L && onMissing == "error") {
-    stop(sprintf("%s: %d variant id(s) not present in the LD sketch panel.",
-                 label, nMissing))
+  # An unmatched requested id falls into one of two very different cases:
+  #   (a) genuine absence  -- its (chrom, pos) is not on the panel at all;
+  #   (b) flip-twin collapse -- its (chrom, pos) IS on the panel, but
+  #       matchVariants(allowFlip) treats a variant and its ref/alt-swapped twin
+  #       as the same variant and keeps the sibling, so the pair merged to one
+  #       match. This happens when the requested set carries both orientations of
+  #       one variant (e.g. a panel that stores both, imputed by RAISS as two).
+  # Case (a) is a real data problem; case (b) is benign -- the variant is present,
+  # just deduplicated. Only (a) counts toward the onMissing="error" stop. The LD
+  # matrix is built on the matched (deduplicated) set either way, and the survivors
+  # are exposed via keptVariantIds so callers can realign their z/N vectors.
+  unmatched <- setdiff(seq_along(variantIds), m$idxA)
+  if (length(unmatched) > 0L) {
+    reqP     <- parseVariantId(variantIds[unmatched])
+    reqKey   <- paste(canonChrom(as.character(reqP$chrom)), reqP$pos)
+    panelKey <- paste(canonChrom(as.character(snpInfo$CHR)), snpInfo$BP)
+    nAbsent  <- sum(!(reqKey %in% panelKey))   # (chrom,pos) truly not on panel
+    if (nAbsent > 0L && onMissing == "error") {
+      stop(sprintf("%s: %d variant id(s) not present in the LD sketch panel.",
+                   label, nAbsent))
+    }
   }
   if (length(m$idxA) == 0L) return(NULL)
   o       <- order(m$idxA)            # restore the caller's requested order
@@ -575,9 +592,9 @@ loadLdFromGenotype <- function(genotypePath, region,
   colnames(geno) <- keptIds
   ldMat <- computeLd(geno, method = "sample")
   dimnames(ldMat) <- list(keptIds, keptIds)
-  if (onMissing == "drop") {
-    attr(ldMat, "keptVariantIds") <- keptIds
-  }
+  # Always expose the retained ids: a flip-twin collapse drops entries even in
+  # "error" mode, so every caller needs to be able to realign to the survivors.
+  attr(ldMat, "keptVariantIds") <- keptIds
   ldMat
 }
 

R/sumstatsQc.R

@@ -2594,6 +2594,13 @@ krigingOutlierQc <- function(zScore, R, n, variantIds = NULL,
   # match with chr-prefix tolerance, strand-ambiguous variants kept; errors if
   # any variant is absent from the panel).
   R <- .ldFromSketch(ldSketch, variantIds, label = "summaryStatsQc: zMismatchQc")
+  # A flip-twin collapse leaves R over the deduplicated survivors; drop the
+  # collapsed duplicate rows from df so df$Z and R align for ldMismatchQc.
+  kept <- attr(R, "keptVariantIds")
+  if (!is.null(kept) && length(kept) != length(variantIds)) {
+    df <- df[match(kept, df$SNP), , drop = FALSE]
+    variantIds <- kept
+  }
   qc <- ldMismatchQc(zScore = df$Z, R = R, nSample = getNSamples(ldSketch),
                      method = method)
   # slalom / dentist can leave NA in the outlier column when their
@@ -2941,6 +2948,12 @@ krigingOutlierQc <- function(zScore, R, n, variantIds = NULL,
     nKrIn <- nrow(df)
     R <- .ldFromSketch(ldSketch, df$SNP,
                        label = "summaryStatsQc: kriging prefilter")
+    # A flip-twin (a variant and its ref/alt swap both present) collapses in the
+    # panel match, so R covers the deduplicated survivors; drop the collapsed
+    # duplicate rows from df so df and R stay aligned for krigingOutlierQc.
+    .kept <- attr(R, "keptVariantIds")
+    if (!is.null(.kept) && length(.kept) != nrow(df))
+      df <- df[match(.kept, df$SNP), , drop = FALSE]
     nKrig <- if (!is.null(opts$nForPip) && is.finite(opts$nForPip)) opts$nForPip
              else stats::median(as.numeric(df$N), na.rm = TRUE)
     kr <- krigingOutlierQc(df$Z, R, n = nKrig, variantIds = df$SNP)

tests/testthat/test_fineMappingPipeline.R

@@ -1547,6 +1547,61 @@ test_that("fineMappingPipeline(GwasSumStats): runs end-to-end with mocked RSS fi
   expect_setequal(getMethodNames(res), "susie")
 })
 
+test_that(".fmAlignToKeptLd realigns z/n/af to keptVariantIds (no-op when nothing collapsed)", {
+  vids <- c("1_100_A_G", "1_100_G_A", "1_200_A_G")   # middle is a flip-twin of the first
+  z <- c(1.1, -1.1, 2.2); n <- c(1000L, 1000L, 2000L); af <- c(0.1, 0.9, 0.2)
+  # no collapse: keptVariantIds == vids -> everything passes through unchanged.
+  m0 <- matrix(0, 3, 3); attr(m0, "keptVariantIds") <- vids
+  a0 <- pecotmr:::.fmAlignToKeptLd(m0, vids, z, n, af)
+  expect_identical(a0$variantIds, vids); expect_identical(a0$z, z)
+  expect_identical(a0$n, n); expect_identical(a0$af, af)
+  # collapse: the twin dropped -> z/n/af subset to the survivors, in kept order.
+  kept <- c("1_100_A_G", "1_200_A_G"); m1 <- matrix(0, 2, 2)
+  attr(m1, "keptVariantIds") <- kept
+  a1 <- pecotmr:::.fmAlignToKeptLd(m1, vids, z, n, af)
+  expect_identical(a1$variantIds, kept)
+  expect_identical(a1$z, c(1.1, 2.2))
+  expect_identical(a1$n, c(1000L, 2000L))
+  expect_identical(a1$af, c(0.1, 0.2))
+  # a scalar n (per-study) is passed through unchanged.
+  a2 <- pecotmr:::.fmAlignToKeptLd(m1, vids, z, n = 5000L, af = NULL)
+  expect_identical(a2$n, 5000L); expect_null(a2$af)
+})
+
+test_that("fineMappingPipeline(GwasSumStats): a flip-twin region fits on the deduplicated set", {
+  # panel carries one orientation per position (parseable chr_pos_ref_alt ids).
+  pan <- new("GenotypeHandle", path = "/tmp/ft.gds", format = "gds",
+    snpInfo = data.frame(
+      SNP = c("1_100_A_G", "1_200_A_G", "1_300_A_G"),
+      CHR = rep("1", 3), BP = c(100L, 200L, 300L),
+      A1 = rep("G", 3), A2 = rep("A", 3), stringsAsFactors = FALSE),
+    nSamples = 40L, sampleIds = paste0("s", seq_len(40L)), pgenPtr = NULL)
+  # entry carries 1_100_A_G AND its ref/alt swap 1_100_G_A -> a flip-twin.
+  gr <- GenomicRanges::GRanges("chr1",
+    IRanges::IRanges(c(100L, 100L, 200L, 300L), width = 1L))
+  S4Vectors::mcols(gr) <- S4Vectors::DataFrame(
+    SNP = c("1_100_A_G", "1_100_G_A", "1_200_A_G", "1_300_A_G"),
+    A1 = c("G", "A", "G", "G"), A2 = c("A", "G", "A", "A"),
+    Z = c(1.5, -1.5, 2.0, 0.5), N = rep(1000L, 4))
+  gss <- GwasSumStats(study = "G1", entry = list(gr), genome = "hg19",
+                      ldSketch = pan, qcInfo = list(step1 = "ok"))
+  cap <- new.env()
+  local_mocked_bindings(
+    extractBlockGenotypes = .fmp_mockExtractor(),
+    .fmFitSusieRss = function(z, R, n, token, ...) {
+      cap$zlen <- length(z); cap$rdim <- dim(R)
+      list(token = token, n_variants = length(z)) },
+    .fmPostprocessOne = .fmp_mockPostprocess(),
+    .package = "pecotmr")
+  # must NOT stop with "not present in the LD sketch panel"; the twin collapses.
+  res <- suppressMessages(
+    fineMappingPipeline(gss, methods = "susie", addSusieInf = FALSE))
+  expect_s4_class(res, "GwasFineMappingResult")
+  # 4 requested variants -> 3 after the flip-twin dedup; z and R stay consistent.
+  expect_equal(cap$zlen, 3L)
+  expect_equal(cap$rdim, c(3L, 3L))
+})
+
 # ---- PIP-screen graceful skip: a screened region -> empty result, not error ----
 
 # A GwasSumStats whose (single) entry was emptied by summaryStatsQc's PIP screen:

tests/testthat/test_ld.R

@@ -2848,6 +2848,49 @@ test_that(".ldFromSketch drops variants absent from the panel when onMissing='dr
   expect_equal(unname(diag(m)), c(1, 1))
 })
 
+# =============================================================================
+# Additional coverage: .ldFromSketch tolerates flip-twin collapse (not fatal)
+# =============================================================================
+
+test_that(".ldFromSketch tolerates a flip-twin in the request under onMissing='error'", {
+  skip_if_not_installed("pgenlibr")
+  h  <- readGenotypes(file.path(geno_test_data_dir, "test_variants"), format = "plink2")
+  si <- getSnpInfo(h)
+  swap <- function(id) { p <- strsplit(id, "_")[[1]]; k <- length(p)
+                         paste(c(p[seq_len(k - 2)], p[k], p[k - 1]), collapse = "_") }
+  v <- as.character(si$SNP[1:3])
+  # request carries v[1] AND its ref/alt swap: matchVariants(allowFlip) collapses
+  # the pair. Both resolve to the same panel position, so default onMissing="error"
+  # must NOT stop -- the twin is deduplicated, not absent.
+  m <- pecotmr:::.ldFromSketch(h, c(v[1], swap(v[1]), v[2], v[3]), label = "flip-twin")
+  expect_equal(dim(m), c(3L, 3L))
+  expect_equal(attr(m, "keptVariantIds"), v)             # swapped twin dropped
+  expect_false(swap(v[1]) %in% attr(m, "keptVariantIds"))
+})
+
+test_that(".ldFromSketch still errors on a genuine absence under onMissing='error'", {
+  skip_if_not_installed("pgenlibr")
+  h  <- readGenotypes(file.path(geno_test_data_dir, "test_variants"), format = "plink2")
+  si <- getSnpInfo(h)
+  # a variant whose (chrom, pos) is not on the panel at all is a real absence.
+  absent <- paste("chr21", "999999999", "A", "G", sep = "_")
+  expect_error(
+    pecotmr:::.ldFromSketch(h, c(as.character(si$SNP[1]), absent), label = "absent"),
+    "not present in the LD sketch panel")
+})
+
+test_that("matchVariants collapses only ref/alt swaps, not distinct ALTs at one position", {
+  a   <- "chr21_17513228_C_G"
+  p   <- strsplit(a, "_")[[1]]
+  swp <- paste(c(p[1:2], p[4], p[3]), collapse = "_")     # ref/alt swap of a
+  alt <- paste(c(p[1:3], "T"), collapse = "_")            # different ALT, same pos
+  # a + its swap collapse to one; a + a-different-alt stay two (distinct variants).
+  expect_equal(length(pecotmr:::matchVariants(c(a, swp), c(a, swp),
+                 removeStrandAmbiguous = FALSE)$idxA), 1L)
+  expect_equal(length(pecotmr:::matchVariants(c(a, alt), c(a, alt),
+                 removeStrandAmbiguous = FALSE)$idxA), 2L)
+})
+
 # =============================================================================
 # Additional coverage: .ldFromSketch reconciles chr-prefix conventions
 # =============================================================================