From d082195d32ff8fe434934ee20922e0decffa96b1 Mon Sep 17 00:00:00 2001 From: Yining97 Date: Mon, 7 Sep 2026 00:00:00 -0400 Subject: [PATCH] Drop panel-unsupported variants at QC LD lookups instead of aborting summaryStatsQc aborts when a harmonized GWAS variant has no entry in the pruned LD panel (kriging prefilter: N variant id(s) not present in the LD sketch panel). The panel filter prunes once before harmonization, and harmonization drops variants the pruned panel can't cover -- but co-located indels survive that: the R5 panel puts an insertion and a deletion on the same POS column with tag ids, harmonization joins by (chrom, pos) so the deletion survives, then the kriging and z-mismatch lookups match by (chrom, pos, allele) tuple and can't resolve it once its partner row is pruned. .panelVariantFilter already documents that the drop-or-error call for a panel-absent variant belongs to .ldFromSketch's onMissing; the QC lookups were the only callers still using the default "error". Make .qcKrigingFlip and .applyLdMismatchQcToEntry pass onMissing="drop", align the entry to the panel-supported subset, and record the dropped count in the qcInfo audit. Leave fine-mapping's onMissing="error" as the post-QC invariant. --- R/sumstatsQc.R | 67 +++++++++++++++++++++--- tests/testthat/test_sumstatsQc.R | 88 +++++++++++++++++++++++++++++--- 2 files changed, 141 insertions(+), 14 deletions(-) diff --git a/R/sumstatsQc.R b/R/sumstatsQc.R index 3ad17b83..889b66b6 100644 --- a/R/sumstatsQc.R +++ b/R/sumstatsQc.R @@ -3333,13 +3333,32 @@ krigingOutlierQc <- function( abort("summaryStatsQc: ldMismatchQc requires SNP column on the entry.") } # Panel LD for the entry variants via the shared LD-from-sketch helper - # (tuple match with chr-prefix tolerance, strand-ambiguous variants kept; - # errors if any variant is absent from the panel). + # (tuple match with chr-prefix tolerance, strand-ambiguous variants kept). + # A variant can survive QC while its panel partner is removed by the panel + # MAF/MAC/missingness filter; .panelVariantFilter passes it through and + # leaves the drop to onMissing here, so drop it (aligning df with the + # returned LD) rather than aborting. + nMmIn <- nrow(df) R <- .ldFromSketch( ldSketch, variantIds, - label = "summaryStatsQc: zMismatchQc" + label = "summaryStatsQc: zMismatchQc", + onMissing = "drop" ) + if (is.null(R)) { + return(list( + df = df[0L, , drop = FALSE], outliers = 0L, diagnostics = NULL, + panelUnsupportedDropped = nMmIn + )) + } + keptIds <- attr(R, "keptVariantIds") + attr(R, "keptVariantIds") <- NULL + nPanelDrop <- 0L + if (!is.null(keptIds) && length(keptIds) < nrow(df)) { + nPanelDrop <- nrow(df) - length(keptIds) + df <- filter(df, is_in(.data$SNP, keptIds)) + variantIds <- df$SNP + } qc <- ldMismatchQc( zScore = df$Z, R = R, @@ -3367,7 +3386,8 @@ krigingOutlierQc <- function( list( df = filter(df, !outlierFlags), outliers = sum(outlierFlags), - diagnostics = diagnostics + diagnostics = diagnostics, + panelUnsupportedDropped = nPanelDrop ) } @@ -4028,8 +4048,31 @@ krigingOutlierQc <- function( R <- .ldFromSketch( ldSketch, df$SNP, - label = "summaryStatsQc: kriging prefilter" + label = "summaryStatsQc: kriging prefilter", + onMissing = "drop" ) + # A GWAS variant can survive QC while its LD-panel partner is removed by the + # panel MAF/MAC/missingness filter (common in the study, rare in the panel -- + # deletions especially). .panelVariantFilter passes such variants through and + # leaves the drop to onMissing here, so drop them (aligning df with the + # returned LD) rather than aborting the run. + if (is.null(R)) { + return(list( + df = df[0L, , drop = FALSE], count = 0L, + audit = list(krigingFlipped = 0L, panelUnsupportedDropped = nKrIn) + )) + } + keptIds <- attr(R, "keptVariantIds") + attr(R, "keptVariantIds") <- NULL + nPanelDrop <- 0L + if (!is.null(keptIds) && length(keptIds) < nrow(df)) { + nPanelDrop <- nrow(df) - length(keptIds) + df <- filter(df, is_in(.data$SNP, keptIds)) + .qcEmit( + lbl, "QC track: dropped ", nPanelDrop, " of ", nKrIn, + " variant(s) with no LD-panel entry after panel filtering." + ) + } nKrig <- if (!is.null(opts$nForPip) && is.finite(opts$nForPip)) { opts$nForPip } else { @@ -4054,7 +4097,11 @@ krigingOutlierQc <- function( list( df = df, count = nKr, - audit = list(krigingFlipped = nKr, krigingDiagnostics = kr$diagnostics) + audit = list( + krigingFlipped = nKr, + krigingDiagnostics = kr$diagnostics, + panelUnsupportedDropped = nPanelDrop + ) ) } @@ -4070,6 +4117,14 @@ krigingOutlierQc <- function( ldMismatchOutliersDropped = ldQc$outliers, ldMismatchMethod = opts$zMismatchQc ) + nPanelDrop <- ldQc$panelUnsupportedDropped %||% 0L + if (nPanelDrop > 0L) { + audit$panelUnsupportedDropped <- nPanelDrop + .qcEmit( + lbl, "QC track: dropped ", nPanelDrop, " of ", nMmIn, + " variant(s) with no LD-panel entry after panel filtering." + ) + } if (!is.null(ldQc$diagnostics)) { audit$ldMismatchDiagnostics <- ldQc$diagnostics } diff --git a/tests/testthat/test_sumstatsQc.R b/tests/testthat/test_sumstatsQc.R index 51bb06ed..bff06477 100644 --- a/tests/testthat/test_sumstatsQc.R +++ b/tests/testthat/test_sumstatsQc.R @@ -4498,18 +4498,90 @@ test_that(".applyLdMismatchQcToEntry: errors when SNP column is missing", { ) }) -test_that(".applyLdMismatchQcToEntry: errors on variants absent from the sketch", { +test_that(".applyLdMismatchQcToEntry: drops variants absent from the sketch", { + # A variant can survive study-side QC while its LD-panel partner is removed + # by the panel MAF/MAC/missingness filter, so it is absent from the sketch + # the QC lookup sees. .panelVariantFilter passes such a variant through and + # leaves the drop to .ldFromSketch's onMissing, so the z-mismatch lookup + # drops it (recording the count) instead of aborting the run. + local_mocked_bindings( + extractBlockGenotypes = .ssh_mockExtractor(), + .package = "pecotmr" + ) df <- data.frame( - SNP = c("rs1", "ghost"), - Z = c(1, 2), - N = c(1000, 1000), + SNP = c(paste0("rs", 1:6), "ghost"), + Z = c(seq(1.0, by = 0.5, length.out = 6L), 2.0), + N = rep(1000L, 7L), + stringsAsFactors = FALSE + ) + out <- pecotmr:::.applyLdMismatchQcToEntry( + df, + .ssh_makeHandle(snp_n = 6L), + method = "slalom" + ) + expect_false("ghost" %in% out$df$SNP) + expect_equal(out$panelUnsupportedDropped, 1L) + expect_true(all(out$df$SNP %in% paste0("rs", 1:6))) +}) + +test_that(".qcKrigingFlip: drops a panel-unsupported variant instead of aborting", { + # Same contract at the kriging prefilter: the orphan (no sketch entry) is + # dropped and counted, and the run continues on the panel-supported subset. + local_mocked_bindings( + extractBlockGenotypes = .ssh_mockExtractor(), + .package = "pecotmr" + ) + df <- data.frame( + SNP = c(paste0("rs", 1:6), "ghost"), + Z = c(seq(1.0, by = 0.5, length.out = 6L), 2.0), + N = rep(1000L, 7L), stringsAsFactors = FALSE ) + out <- pecotmr:::.qcKrigingFlip( + df, + .ssh_makeHandle(snp_n = 6L), + opts = list(alleleFlipKriging = TRUE, nForPip = NULL), + lbl = NA_character_ + ) + expect_false("ghost" %in% out$df$SNP) + expect_equal(out$audit$panelUnsupportedDropped, 1L) + expect_equal(nrow(out$df), 6L) +}) + +test_that(".qcKrigingFlip: no-orphan entry is unchanged (drop is a no-op)", { + # When every variant has a sketch entry, nothing is dropped and the result + # matches the pre-change behaviour (panelUnsupportedDropped == 0). + local_mocked_bindings( + extractBlockGenotypes = .ssh_mockExtractor(), + .package = "pecotmr" + ) + df <- data.frame( + SNP = paste0("rs", 1:6), + Z = seq(1.0, by = 0.5, length.out = 6L), + N = rep(1000L, 6L), + stringsAsFactors = FALSE + ) + out <- pecotmr:::.qcKrigingFlip( + df, + .ssh_makeHandle(snp_n = 6L), + opts = list(alleleFlipKriging = TRUE, nForPip = NULL), + lbl = NA_character_ + ) + expect_equal(out$audit$panelUnsupportedDropped, 0L) + expect_equal(nrow(out$df), 6L) + expect_setequal(out$df$SNP, paste0("rs", 1:6)) +}) + +test_that("fine-mapping keeps onMissing='error': a genuinely absent variant still aborts", { + # The post-QC invariant: QC now drops panel-unsupported variants, so every + # variant reaching fine-mapping has a panel entry. .ldFromSketch's default + # onMissing='error' (used by the fine-mapping LD build) MUST still abort on + # a genuinely absent variant, guaranteeing that invariant. expect_error( - pecotmr:::.applyLdMismatchQcToEntry( - df, - .ssh_makeHandle(), - method = "dentist" + pecotmr:::.ldFromSketch( + .ssh_makeHandle(snp_n = 6L), + c("rs1", "ghost"), + label = "fine-mapping" ), "not present in the LD sketch panel" )