RSS QC summary logging + fine-mapping vignette on current API

R/allele_qc.R

@@ -32,7 +32,13 @@ NULL
 #' @return An \code{AlleleQcResult} S4 object. Use
 #'   \code{getHarmonizedData()} to recover the post-QC variant
 #'   data.frame and \code{getQcSummary()} to inspect the per-variant
-#'   merge/flip/strand diagnostics.
+#'   merge/flip/strand diagnostics. The object also carries a
+#'   \code{"qcCounts"} attribute (\code{attr(x, "qcCounts")}) — a list with
+#'   \code{considered}, \code{signFlip}, \code{strandFlip} (corrections applied
+#'   to retained variants), \code{kept}, \code{dropped}, and the drop breakdown
+#'   \code{droppedIndel} / \code{droppedAmbiguous} / \code{droppedOther} — used
+#'   for QC summary logging. Existing callers that read the data.frame are
+#'   unaffected.
 #' @importFrom magrittr %>%
 #' @importFrom dplyr mutate inner_join filter pull select everything row_number if_else any_of all_of rename
 #' @importFrom vctrs vec_duplicate_detect
@@ -98,7 +104,11 @@ matchRefPanel <- function(targetData, refVariants, colToFlip = NULL,
 
   if (nrow(matchResult) == 0) {
 	warning("No matching variants found between target data and reference variants.")
-	return(AlleleQcResult(harmonizedData = matchResult, qcSummary = matchResult))
+	emptyOut <- AlleleQcResult(harmonizedData = matchResult, qcSummary = matchResult)
+	attr(emptyOut, "qcCounts") <- list(considered = 0L, signFlip = 0L, strandFlip = 0L,
+	                                   kept = 0L, dropped = 0L, droppedIndel = 0L,
+	                                   droppedAmbiguous = 0L, droppedOther = 0L)
+	return(emptyOut)
   }
     # match target & ref by chrom and position
   matchResult = matchResult %>%
@@ -209,7 +219,31 @@ matchRefPanel <- function(targetData, refVariants, colToFlip = NULL,
 	stop("Duplicated variants with different values found. Please check the input data and determine which to keep.")
   }
 
-  return(AlleleQcResult(harmonizedData = result, qcSummary = matchResult))
+  # QC summary counts for logging. Corrections (sign/strand flip) are counted on
+  # retained variants; drops are split by reason. Computed from the per-variant
+  # flags still on matchResult (default removeUnmatched = TRUE path); attached as
+  # a `qcCounts` attribute so callers that read the data frame are unaffected.
+  qcCounts <- if (all(c("sign_flip", "strand_flip", "keep", "INDEL", "strand_unambiguous") %in% names(matchResult))) {
+    cnt <- list(
+      considered = nrow(matchResult),
+      signFlip = sum(matchResult$sign_flip & matchResult$keep, na.rm = TRUE),
+      strandFlip = sum(matchResult$strand_flip & matchResult$keep, na.rm = TRUE),
+      kept = sum(matchResult$keep, na.rm = TRUE),
+      dropped = sum(!matchResult$keep, na.rm = TRUE),
+      droppedIndel = sum(!matchResult$keep & matchResult$INDEL, na.rm = TRUE),
+      droppedAmbiguous = sum(!matchResult$keep & !matchResult$INDEL &
+                               matchResult$strand_flip & !matchResult$strand_unambiguous, na.rm = TRUE)
+    )
+    cnt$droppedOther <- cnt$dropped - cnt$droppedIndel - cnt$droppedAmbiguous
+    cnt
+  } else {
+    list(considered = nrow(matchResult), signFlip = 0L, strandFlip = 0L,
+         kept = nrow(matchResult), dropped = 0L, droppedIndel = 0L,
+         droppedAmbiguous = 0L, droppedOther = 0L)
+  }
+  out <- AlleleQcResult(harmonizedData = result, qcSummary = matchResult)
+  attr(out, "qcCounts") <- qcCounts
+  return(out)
 }
 
 #' @rdname matchRefPanel

R/sumstats_qc.R

@@ -58,6 +58,7 @@ rssBasicQc <- function(sumstats, ldData, skipRegion = NULL, keepIndel = TRUE,
     removeStrandAmbiguous = TRUE
   )
 
+  bqCounts <- attr(alleleFlip, "qcCounts")
   qced <- getHarmonizedData(alleleFlip)
   if (!is.null(skipRegion)) {
     skipTable <- tibble(region = skipRegion) %>%
@@ -81,13 +82,15 @@ rssBasicQc <- function(sumstats, ldData, skipRegion = NULL, keepIndel = TRUE,
   nRef <- ldData@nRef
 
   if (!isTRUE(returnLdMat)) {
-    return(QcResult(
+    resNoLd <- QcResult(
       ldData = NULL,
       rssInput = list(sumstats = sumstatsProcessed, n = NA_real_, varY = NA_real_),
       preprocess = list(),
       outlierNumber = 0L,
       skipped = FALSE
-    ))
+    )
+    attr(resNoLd, "qcCounts") <- bqCounts
+    return(resNoLd)
   }
 
   # Align and subset LD by mapping core IDs (strip trailing build suffix) to exact LD IDs
@@ -113,14 +116,16 @@ rssBasicQc <- function(sumstats, ldData, skipRegion = NULL, keepIndel = TRUE,
 
   ldMatProcessed <- ldMatrix[sumstatsProcessed$variant_id, sumstatsProcessed$variant_id, drop = FALSE]
 
-  QcResult(
+  resFinal <- QcResult(
     ldData = .qcLdDataFromMatrix(ldMatProcessed, sumstatsProcessed$variant_id,
                                  hasGenotype = FALSE, nRef = nRef),
     rssInput = list(sumstats = sumstatsProcessed, n = NA_real_, varY = NA_real_),
     preprocess = list(),
     outlierNumber = 0L,
     skipped = FALSE
   )
+  attr(resFinal, "qcCounts") <- bqCounts
+  resFinal
 }
 
 
@@ -537,6 +542,10 @@ summaryStatsQc <- function(sumstats, ldData, n = NULL,
                       returnLdMat = !hasGenotype)
   sumstats <- getRssInput(basic)$sumstats
   basicLd <- getLdData(basic)
+  # QC logging: counts surfaced by matchRefPanel (via rssBasicQc) plus the study
+  # input size, used for "N of M" per-step messages and the per-study rollup.
+  hCounts <- attr(basic, "qcCounts")
+  nStudyIn <- nrow(rssInput$sumstats)
   R_mat <- if (is.null(basicLd)) NULL else getCorrelation(basicLd)
   n <- rssInput$n
   varY <- rssInput$varY
@@ -582,8 +591,15 @@ summaryStatsQc <- function(sumstats, ldData, n = NULL,
       nRef = ldDataForQc@nRef
     )
   }
-  message("QC track: basic harmonization retained ", nrow(sumstats),
-          " variants for summary-stat study ", study, ".")
+  harmKept <- nrow(sumstats)
+  harmDropped <- nStudyIn - harmKept
+  message("QC track: basic harmonization kept ", harmKept, " of ", nStudyIn,
+          " variant(s) for summary-stat study ", study,
+          if (!is.null(hCounts)) {
+            paste0(" (corrected: sign-flipped ", hCounts$signFlip,
+                   ", strand-flipped ", hCounts$strandFlip,
+                   "; dropped ", harmDropped, ").")
+          } else ".")
   if (nrow(sumstats) < 2) {
     return(skippedResult(sumstats, referenceForVariants(sumstats$variant_id),
                          "fewer than two variants remain after basic QC"))
@@ -607,10 +623,15 @@ summaryStatsQc <- function(sumstats, ldData, n = NULL,
   }
 
   outlierNumber <- 0L
+  krRemoved <- 0L
+  mismatchRemoved <- 0L
+  imputedBefore <- NA_integer_
+  imputedAfter <- NA_integer_
   # Optional kriging LD-consistency prefilter (opt-in). Runs before the heavier
   # SLALOM/DENTIST QC, or standalone when zMismatchQc == "none".
   if (isTRUE(alleleFlipKriging) && nrow(sumstats) >= 2) {
     message("QC track: running kriging LD-consistency prefilter for summary-stat study ", study, ".")
+    nKrIn <- nrow(sumstats)
     krLd <- ldDataWithLocalR(sumstats)
     rKr <- getCorrelation(krLd)
     rKr <- rKr[sumstats$variant_id, sumstats$variant_id, drop = FALSE]
@@ -621,11 +642,13 @@ summaryStatsQc <- function(sumstats, ldData, n = NULL,
       if (!hasGenotype) R_mat <- R_mat[sumstats$variant_id, sumstats$variant_id, drop = FALSE]
       outlierNumber <- outlierNumber + nKr
     }
-    message("QC track: kriging prefilter removed ", nKr,
+    krRemoved <- nKr
+    message("QC track: kriging prefilter removed ", nKr, " of ", nKrIn,
             " LD-inconsistent variant(s) for summary-stat study ", study, ".")
   }
   if (!is.null(zMismatchQc) && !identical(zMismatchQc, "none")) {
     message("QC track: running ", zMismatchQc, " z-score/LD-mismatch QC for summary-stat study ", study, ".")
+    nMmIn <- nrow(sumstats)
     qc <- summaryStatsQc(sumstats = sumstats, ldData = ldDataWithLocalR(sumstats),
                          n = n, method = zMismatchQc)
     qcRss <- getRssInput(qc)
@@ -634,11 +657,13 @@ summaryStatsQc <- function(sumstats, ldData, n = NULL,
     R_mat <- if (is.null(qcLd)) NULL else getCorrelation(qcLd)
     ldMismatchOutliers <- getOutlierNumber(qc)
     outlierNumber <- outlierNumber + ldMismatchOutliers
-    message("QC track: removed ", ldMismatchOutliers,
+    mismatchRemoved <- ldMismatchOutliers
+    message("QC track: removed ", ldMismatchOutliers, " of ", nMmIn,
             " LD-mismatch outlier(s) for summary-stat study ", study, ".")
   }
   if (isTRUE(impute)) {
     message("QC track: running imputation for summary-stat study ", study, ".")
+    imputedBefore <- nrow(sumstats)
     imputed <- if (hasGenotype) {
       X_ref_scaled <- scale(X_ref)
       X_ref_scaled[is.na(X_ref_scaled)] <- 0
@@ -658,7 +683,28 @@ summaryStatsQc <- function(sumstats, ldData, n = NULL,
     }
     sumstats <- imputed$resultFilter
     if (!is.null(imputed$ldMat)) R_mat <- imputed$ldMat
+    imputedAfter <- nrow(sumstats)
+    message("QC track: imputation ", imputedBefore, " -> ", imputedAfter,
+            " variant(s) (net ", sprintf("%+d", imputedAfter - imputedBefore),
+            ") for summary-stat study ", study, ".")
   }
+  # Per-study QC rollup: corrected (sign/strand flip, retained), removed (drops at
+  # each stage), imputed (added). Kept as distinct categories because imputation
+  # adds variants back, so "in -> out" is not monotonic. Skipped stages omitted.
+  correctedSeg <- if (!is.null(hCounts)) {
+    paste0("sign-flip ", hCounts$signFlip, ", strand-flip ", hCounts$strandFlip)
+  } else "n/a"
+  removedSegs <- character(0)
+  removedSegs <- c(removedSegs, paste0("harmonization ", harmDropped))
+  if (isTRUE(alleleFlipKriging)) removedSegs <- c(removedSegs, paste0("kriging ", krRemoved))
+  if (!identical(zMismatchQc, "none")) removedSegs <- c(removedSegs, paste0("mismatch ", mismatchRemoved))
+  impSeg <- if (isTRUE(impute) && !is.na(imputedAfter)) {
+    paste0(" | imputed ", sprintf("%+d", imputedAfter - imputedBefore))
+  } else ""
+  message("QC summary [", study, "]: ", nStudyIn, " in -> ", nrow(sumstats),
+          " out | corrected: ", correctedSeg,
+          " | removed: ", paste(removedSegs, collapse = ", "), impSeg)
+
   finalVars <- sumstats$variant_id
   finalLdMat <- if (hasGenotype) referenceForVariants(finalVars) else R_mat
   preprocessSumstats <- preprocess$sumstats

tests/testthat/test_allele_qc.R

@@ -472,3 +472,35 @@ test_that("colToComplement errors on a missing column name", {
     "not found in targetData"
   )
 })
+
+# ===========================================================================
+# QC summary counts (rss-vignette-and-qc-logging)
+# ===========================================================================
+
+test_that("matchRefPanel surfaces sign/strand/dropped counts via qcCounts attribute", {
+  # Deterministic fixture at four shared positions:
+  #   100 exact match | 200 sign flip | 300 strand flip (A/G, unambiguous) | 400 allele mismatch (dropped)
+  target <- data.frame(
+    chrom = c(1, 1, 1, 1), pos = c(100, 200, 300, 400),
+    A2 = c("A", "A", "A", "A"), A1 = c("G", "G", "G", "G"),
+    z = c(1, 2, 3, 4), stringsAsFactors = FALSE
+  )
+  ref <- data.frame(
+    chrom = c(1, 1, 1, 1), pos = c(100, 200, 300, 400),
+    A2 = c("A", "G", "T", "C"), A1 = c("G", "A", "C", "A"), stringsAsFactors = FALSE
+  )
+  res <- matchRefPanel(target, ref, colToFlip = "z", matchMinProp = 0)
+
+  # Default return shape unchanged: kept variants only in harmonizedData.
+  expect_s4_class(res, "AlleleQcResult")
+  expect_equal(nrow(getHarmonizedData(res)), 3L)  # exact + sign + strand kept; mismatch dropped
+
+  # Counts surfaced as an attribute (additive; non-RSS callers ignore it).
+  cnt <- attr(res, "qcCounts")
+  expect_false(is.null(cnt))
+  expect_equal(cnt$considered, 4L)
+  expect_equal(cnt$signFlip, 1L)
+  expect_equal(cnt$strandFlip, 1L)
+  expect_equal(cnt$kept, 3L)
+  expect_equal(cnt$dropped, 1L)
+})

tests/testthat/test_sumstats_qc.R

@@ -415,7 +415,7 @@ test_that("summaryStatsQc basic genotype-backed path does not compute LD", {
   expect_message(
     result <- summaryStatsQc(rssInput = rss_input, ldData = LD_data_geno,
                                zMismatchQc = "none", impute = FALSE),
-    "basic harmonization retained"
+    "basic harmonization kept"
   )
   result_ld <- getLdData(result)
   result_geno <- getGenotypes(result_ld)
@@ -526,7 +526,7 @@ test_that("summaryStatsQc treats NULL qc_method as basic-only none", {
       zMismatchQc = NULL,
       impute = FALSE
     ),
-    "basic harmonization retained"
+    "basic harmonization kept"
   )
   expect_equal(nrow(getRssInput(result)$sumstats), nrow(td$sumstats))
 })
@@ -673,7 +673,7 @@ test_that("kriging is not run by default (alleleFlipKriging = FALSE)", {
       rssInput = rssInput, ldData = td$ldData,
       zMismatchQc = "none", impute = FALSE
     ),
-    "basic harmonization retained"
+    "basic harmonization kept"
   )
   expect_equal(nrow(getRssInput(result)$sumstats), nrow(td$sumstats))
 })
@@ -701,3 +701,46 @@ test_that("alleleFlipKriging drops a planted outlier and composes with zMismatch
   expect_equal(seenN, nrow(td$sumstats) - 1)   # kriging removed one before dentist
   expect_equal(getOutlierNumber(result), 1)
 })
+
+# ===========================================================================
+# QC summary logging (rss-vignette-and-qc-logging)
+# ===========================================================================
+
+test_that("QC track reports kept-of-considered and a per-study rollup", {
+  td <- make_test_sumstats_ld(6)
+  msgs <- capture_messages(
+    summaryStatsQc(rssInput = list(sumstats = td$sumstats, n = 1000, varY = 1),
+                   ldData = td$ldData, zMismatchQc = "none", impute = FALSE,
+                   returnOnSkip = "preprocess", study = "studyA")
+  )
+  joined <- paste(msgs, collapse = "")
+  expect_match(joined, "harmoni", ignore.case = TRUE)
+  expect_match(joined, "of [0-9]+")        # "N of M" denominator framing
+  expect_match(joined, "QC summary")       # per-study rollup line
+})
+
+test_that("kriging outliers are reported as removed-of-N (denominator), never flipped", {
+  td <- make_test_sumstats_ld(6)
+  msgs <- capture_messages(
+    summaryStatsQc(rssInput = list(sumstats = td$sumstats, n = 1000, varY = 1),
+                   ldData = td$ldData, zMismatchQc = "none",
+                   alleleFlipKriging = TRUE, impute = FALSE,
+                   returnOnSkip = "preprocess", study = "studyA")
+  )
+  joined <- paste(msgs, collapse = "")
+  expect_match(joined, "kriging", ignore.case = TRUE)
+  expect_match(joined, "removed [0-9]+ of [0-9]+", ignore.case = TRUE)
+})
+
+test_that("skipped QC steps are not reported as actions", {
+  td <- make_test_sumstats_ld(6)
+  msgs <- capture_messages(
+    summaryStatsQc(rssInput = list(sumstats = td$sumstats, n = 1000, varY = 1),
+                   ldData = td$ldData, zMismatchQc = "none",
+                   alleleFlipKriging = FALSE, impute = FALSE,
+                   returnOnSkip = "preprocess", study = "studyA")
+  )
+  joined <- paste(msgs, collapse = "")
+  expect_match(joined, "QC summary")
+  expect_false(grepl("imputed \\+[1-9]", joined))   # imputation disabled -> no positive imputed count
+})