RSS fine-mapping: pass finemappingOpts straight through to susie_rss

Per review: the fine-mapping parameters are not pipeline arguments. susieRssPipeline
and rssAnalysisPipeline drop L/L_greedy/min_abs_corr/median_abs_corr/R_finite/R_mismatch as
named args; finemappingOpts is a free-form list forwarded as-is into susie_rss, so unset keys
inherit susieR defaults and a run with none set matches a manual susie_rss call.

min_abs_corr (default 0.8) and median_abs_corr (default NULL) are credible-set purity options
applied at susie_get_cs, not the fit, so they are isolated from finemappingOpts and routed to
susie_get_cs across the reported coverages. No L_greedy clamp.

median_abs_corr requires a susieR that provides it (GitHub master); it is passed to susie_get_cs
only when the user sets it, so older susieR keeps working and the DESCRIPTION constraint stays
>= 0.16.2 (the pixi CI environment pins an older conda susieR). Full suite 3076 pass / 0 fail; on
chr21 AD_Bellenguez the unset pipeline matches a manual susie_rss run (max abs PIP diff 0), explicit
step-1 params reproduce the step-1 credible sets exactly, and median_abs_corr routing is confirmed.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

Author: Yining97

R/susie_wrapper.R

@@ -238,6 +238,7 @@ postprocessFinemappingFits <- function(fits, dataX, dataY = NULL,
                                        region = NULL,
                                        priorEffTol = 1e-9,
                                        minAbsCorr = 0.8,
+                                       medianAbsCorr = NULL,
                                        csInput = NULL) {
   fits <- fits[!vapply(fits, is.null, logical(1))]
   if (length(fits) == 0) stop("At least one fine-mapping fit must be supplied.")
@@ -257,6 +258,7 @@ postprocessFinemappingFits <- function(fits, dataX, dataY = NULL,
       signalCutoff = signalCutoff, otherQuantities = otherQuantities,
       region = region,
       priorEffTol = priorEffTol, minAbsCorr = minAbsCorr,
+      medianAbsCorr = medianAbsCorr,
       csInput = csInput
     )
   })
@@ -324,6 +326,7 @@ postprocessFinemappingFit.susiF <- function(fit, method = "fsusie", csInput = NU
                                              region = NULL,
                                              priorEffTol = 1e-9,
                                              minAbsCorr = 0.8,
+                                             medianAbsCorr = NULL,
                                              csInput = c("X", "Xcorr", "fsusie")) {
   csInput <- match.arg(csInput)
   variantNames <- extractVariantNames(fit)
@@ -332,7 +335,7 @@ postprocessFinemappingFit.susiF <- function(fit, method = "fsusie", csInput = NU
   csTables <- computeCsTables(
     fit, dataX = dataX, coverage = coverage,
     secondaryCoverage = secondaryCoverage, method = method,
-    csInput = csInput, minAbsCorr = minAbsCorr
+    csInput = csInput, minAbsCorr = minAbsCorr, medianAbsCorr = medianAbsCorr
   )
   topLoci <- buildTopLoci(
     fit, csTables, variantNames = variantNames, sumstats = sumstats,
@@ -432,7 +435,7 @@ selectEffects <- function(fit, priorEffTol = 1e-9) {
 computeCsTables <- function(fit, dataX, coverage = NULL,
                             secondaryCoverage = c(0.7, 0.5),
                             method = "susie", csInput = c("X", "Xcorr", "fsusie"),
-                            minAbsCorr = 0.8) {
+                            minAbsCorr = 0.8, medianAbsCorr = NULL) {
   csInput <- match.arg(csInput)
   primaryCoverage <- coverage
   if (is.null(primaryCoverage)) primaryCoverage <- fit$sets$requested_coverage
@@ -441,15 +444,16 @@ computeCsTables <- function(fit, dataX, coverage = NULL,
   coverages <- coverages[!is.na(coverages)]
 
   tables <- lapply(coverages, function(cov) {
-    computeCsTable(fit, dataX, coverage = cov, csInput = csInput, minAbsCorr = minAbsCorr)
+    computeCsTable(fit, dataX, coverage = cov, csInput = csInput,
+                   minAbsCorr = minAbsCorr, medianAbsCorr = medianAbsCorr)
   })
   names(tables) <- vapply(coverages, formatCsColumn, character(1), method = method)
   attr(tables, "coverage") <- coverages
   tables
 }
 
 computeCsTable <- function(fit, dataX, coverage, csInput = c("X", "Xcorr", "fsusie"),
-                           minAbsCorr = 0.8) {
+                           minAbsCorr = 0.8, medianAbsCorr = NULL) {
   csInput <- match.arg(csInput)
   if (csInput == "fsusie") {
     sets <- tryCatch(
@@ -470,12 +474,18 @@ computeCsTable <- function(fit, dataX, coverage, csInput = c("X", "Xcorr", "fsus
     return(list(sets = sets, cs_corr = csCorr, pip = fit$pip))
   }
 
+  # Purity thresholds for credible-set extraction. min_abs_corr / median_abs_corr
+  # are isolated from finemappingOpts upstream and routed here; pass each only
+  # when set. `fit` is passed positionally as `res`.
+  csArgs <- list(coverage = coverage)
+  if (!is.null(minAbsCorr)) csArgs$min_abs_corr <- minAbsCorr
+  if (!is.null(medianAbsCorr)) csArgs$median_abs_corr <- medianAbsCorr
   if (csInput == "X") {
-    sets <- susie_get_cs(fit, X = dataX, coverage = coverage, min_abs_corr = minAbsCorr)
+    sets <- do.call(susie_get_cs, c(list(fit), csArgs, list(X = dataX)))
     out <- list(sets = sets, pip = fit$pip)
     out$cs_corr <- get_cs_correlation(out, X = dataX)
   } else {
-    sets <- susie_get_cs(fit, Xcorr = dataX, coverage = coverage, min_abs_corr = minAbsCorr)
+    sets <- do.call(susie_get_cs, c(list(fit), csArgs, list(Xcorr = dataX)))
     out <- list(sets = sets, pip = fit$pip)
     out$cs_corr <- get_cs_correlation(out, Xcorr = dataX)
   }
@@ -896,8 +906,13 @@ adjustSusieWeights <- function(twasWeightsResults, keepVariants, runAlleleQc = T
 #' @param ldMat LD correlation matrix. Mutually exclusive with xMat.
 #' @param xMat Genotype matrix (samples x variants). Mutually exclusive with ldMat.
 #' @param n Sample size.
-#' @param L Maximum number of causal configurations (default: 30).
-#' @param lGreedy Initial greedy number of causal configurations (default: 5).
+#' @param finemappingOpts Free-form list of fine-mapping options forwarded as-is
+#'   into \code{susieR::susie_rss()} (e.g. \code{L}, \code{L_greedy},
+#'   \code{R_finite}, \code{R_mismatch}). A key supplied is passed through; a key
+#'   omitted inherits \code{susie_rss}'s own default. Two purity keys are special:
+#'   \code{min_abs_corr} (default \code{0.8}) and \code{median_abs_corr} (default
+#'   \code{NULL}) are isolated from this list and routed to
+#'   \code{susieR::susie_get_cs()} for credible-set extraction, not the fit.
 #' @param analysisMethod Iteration mode for the \code{"susie_rss"} fit:
 #'   \code{"susie_rss"} (default, normal IBSS), \code{"single_effect"} (L=1,
 #'   single iteration), or \code{"bayesian_conditional_regression"}
@@ -919,12 +934,6 @@ adjustSusieWeights <- function(twasWeightsResults, keepVariants, runAlleleQc = T
 #' @param coverage Coverage level (default: 0.95).
 #' @param secondaryCoverage Secondary coverage levels (default: c(0.7, 0.5)).
 #' @param signalCutoff PIP cutoff for selecting top loci (default: 0.1).
-#' @param minAbsCorr Minimum absolute correlation for CS purity (default: 0.8).
-#' @param rFinite Controls variance inflation to account for estimating
-#'   the R matrix from a finite reference panel. NULL (default): no
-#'   variance inflation. Passed directly to susie_rss.
-#' @param rMismatch LD mismatch correction method passed directly to susie_rss.
-#'   Default NULL disables mismatch correction.
 #' @param ... Additional parameters passed to susie_rss. Supplying
 #'   \code{var_y} here, together with \code{beta} and \code{se} columns in
 #'   \code{sumstats}, selects the \code{bhat/shat/var_y} sufficient-statistic
@@ -940,19 +949,24 @@ adjustSusieWeights <- function(twasWeightsResults, keepVariants, runAlleleQc = T
 #' @importFrom dplyr arrange select
 #' @export
 susieRssPipeline <- function(sumstats, ldMat = NULL, xMat = NULL, n = NULL,
-                             L = 30, lGreedy = 5,
                              analysisMethod = c("susie_rss", "single_effect", "bayesian_conditional_regression"),
                              methods = NULL,
                              addSusieInf = TRUE,
                              coverage = 0.95,
                              secondaryCoverage = c(0.7, 0.5),
                              signalCutoff = 0.1,
-                             minAbsCorr = 0.8,
-                             rFinite = NULL, rMismatch = NULL, ...) {
+                             finemappingOpts = list(), ...) {
   analysisMethod <- match.arg(analysisMethod)
   if (is.null(ldMat) && is.null(xMat)) stop("Either ldMat or xMat must be provided.")
   if (!is.null(ldMat) && !is.null(xMat)) stop("Only one of ldMat or xMat should be provided, not both.")
-  if (!is.null(lGreedy)) lGreedy <- min(lGreedy, L)
+  # Fine-mapping options are a free-form passthrough to susie_rss. Isolate the
+  # two credible-set purity options (they act at susie_get_cs extraction, not the
+  # fit) and route them separately; everything else flows into the susie_rss call.
+  if (!is.list(finemappingOpts)) stop("finemappingOpts must be a list.")
+  minAbsCorr <- finemappingOpts$min_abs_corr %||% 0.8
+  medianAbsCorr <- finemappingOpts$median_abs_corr   # NULL when absent (susie_get_cs default)
+  finemappingOpts$min_abs_corr <- NULL
+  finemappingOpts$median_abs_corr <- NULL
 
   # Resolve effective methods. NULL => legacy single-method via analysisMethod.
   validRssMethods <- c("susie_rss", "susie_inf_rss", "susie_ash_rss")
@@ -990,10 +1004,12 @@ susieRssPipeline <- function(sumstats, ldMat = NULL, xMat = NULL, n = NULL,
     names(z) <- rownames(sumstats)
   }
 
-  dots <- list(...)
-  varY <- dots$varY
-  dots$varY <- NULL
-  if (!is.null(dots$bhat) || !is.null(dots$shat)) {
+  # Free-form passthrough to susie_rss: finemappingOpts (purity keys already
+  # isolated above) plus any extra named args. Forwarded verbatim into the fit.
+  fitOpts <- c(finemappingOpts, list(...))
+  varY <- fitOpts$varY
+  fitOpts$varY <- NULL
+  if (!is.null(fitOpts$bhat) || !is.null(fitOpts$shat)) {
     stop("Pass summary effects as 'beta' and 'se' columns in sumstats; ",
          "susieRssPipeline constructs bhat and shat internally.")
   }
@@ -1017,33 +1033,34 @@ susieRssPipeline <- function(sumstats, ldMat = NULL, xMat = NULL, n = NULL,
     shat <- sumstats$se
     names(bhat) <- names(shat) <- names(z)
     common <- c(list(bhat = bhat, shat = shat, var_y = varY, n = n,
-                     coverage = coverage, R_finite = rFinite,
-                     R_mismatch = rMismatch), dots)
+                     coverage = coverage), fitOpts)
   } else {
-    common <- c(list(z = z, n = n, coverage = coverage,
-                     R_finite = rFinite, R_mismatch = rMismatch), dots)
+    common <- c(list(z = z, n = n, coverage = coverage), fitOpts)
   }
   if (!is.null(xMat)) common$X <- xMat else common$R <- ldMat
 
+  # Method-specific overrides are applied with modifyList so they win over any
+  # passthrough key of the same name (e.g. single_effect forces L = 1) without a
+  # duplicate-argument error.
   fitOneSusieRss <- function() {
     if (analysisMethod == "single_effect") {
-      do.call(susie_rss, c(common, list(L = 1, L_greedy = NULL, max_iter = 1)))
+      do.call(susie_rss, modifyList(common, list(L = 1, L_greedy = NULL, max_iter = 1)))
     } else if (analysisMethod == "bayesian_conditional_regression") {
-      do.call(susie_rss, c(common, list(L = L, L_greedy = lGreedy, max_iter = 1)))
+      do.call(susie_rss, modifyList(common, list(max_iter = 1)))
     } else {
-      do.call(susie_rss, c(common, list(L = L, L_greedy = lGreedy)))
+      do.call(susie_rss, common)
     }
   }
   fitOneSusieInfRss <- function() {
-    do.call(susie_rss, c(common, list(L = L, L_greedy = lGreedy,
-                                       unmappable_effects = "inf",
-                                       convergence_method = "pip",
-                                       refine = FALSE, model_init = NULL)))
+    do.call(susie_rss, modifyList(common,
+                         list(unmappable_effects = "inf",
+                              convergence_method = "pip",
+                              refine = FALSE, model_init = NULL)))
   }
   fitOneSusieAshRss <- function() {
-    do.call(susie_rss, c(common, list(L = L, L_greedy = lGreedy,
-                                       unmappable_effects = "ash",
-                                       convergence_method = "pip")))
+    do.call(susie_rss, modifyList(common,
+                         list(unmappable_effects = "ash",
+                              convergence_method = "pip")))
   }
 
   fittedModels <- list()
@@ -1056,11 +1073,11 @@ susieRssPipeline <- function(sumstats, ldMat = NULL, xMat = NULL, n = NULL,
       identical(fitMethods, "bayesian_conditional_regression")) {
     if (chainInfToSusieRss) {
       chainedArgs <- prepareSusieFromInfArgs(
-        list(L = L, L_greedy = lGreedy),
+        list(L = common$L, L_greedy = common$L_greedy),
         fittedModels[["susie_inf_rss"]], refineDefault = TRUE,
         unmappableEffects = "none"
       )
-      rssFit <- do.call(susie_rss, c(common, chainedArgs))
+      rssFit <- do.call(susie_rss, modifyList(common, chainedArgs))
     } else {
       rssFit <- fitOneSusieRss()
     }
@@ -1071,11 +1088,11 @@ susieRssPipeline <- function(sumstats, ldMat = NULL, xMat = NULL, n = NULL,
   if ("susie_ash_rss" %in% fitMethods) {
     if (chainInfToSusieAshRss) {
       chainedArgs <- prepareSusieFromInfArgs(
-        list(L = L, L_greedy = lGreedy),
+        list(L = common$L, L_greedy = common$L_greedy),
         fittedModels[["susie_inf_rss"]], refineDefault = NULL,
         unmappableEffects = "ash"
       )
-      ashFit <- do.call(susie_rss, c(common, chainedArgs))
+      ashFit <- do.call(susie_rss, modifyList(common, chainedArgs))
     } else {
       ashFit <- fitOneSusieAshRss()
     }
@@ -1113,6 +1130,7 @@ susieRssPipeline <- function(sumstats, ldMat = NULL, xMat = NULL, n = NULL,
     secondaryCoverage = secondaryCoverage,
     signalCutoff = signalCutoff,
     minAbsCorr = minAbsCorr,
+    medianAbsCorr = medianAbsCorr,
     csInput = ppCsInput
   )
   # Primary method preference: "susie_rss" > other names > first fit

R/univariate_pipeline.R

@@ -519,15 +519,17 @@ regionDataToSusieRssInput <- function(rssInput, ldData) {
 #' @param addSusieInf Logical controlling chained init when
 #'   \code{"susie_inf_rss"} is in \code{methods} alongside
 #'   \code{"susie_rss"} and/or \code{"susie_ash_rss"}. Default \code{TRUE}.
-#' @param finemappingOpts List of fine-mapping options (L, L_greedy, coverage,
-#'   signal_cutoff, min_abs_corr).
+#' @param finemappingOpts Free-form list of fine-mapping options. \code{coverage}
+#'   and \code{signal_cutoff} are pipeline-reporting choices kept here; everything
+#'   else (e.g. \code{L}, \code{L_greedy}, \code{R_finite}, \code{R_mismatch}) is
+#'   forwarded as-is into \code{susieR::susie_rss()} — supplied keys pass through,
+#'   omitted keys inherit susieR defaults (a run with the fit params unset matches
+#'   a manual \code{susie_rss()} call). The purity keys \code{min_abs_corr}
+#'   (default \code{0.8}) and \code{median_abs_corr} (default \code{NULL}) are
+#'   isolated and routed to \code{susieR::susie_get_cs()} instead of the fit.
 #' @param impute Whether to impute missing variants via RAISS (default TRUE).
 #' @param imputeOpts List of imputation options (rcond, R2_threshold, minimum_ld, lamb).
 #' @param pipCutoffToSkip PIP threshold for early stopping (default 0, no skip).
-#' @param rFinite Controls variance inflation to account for finite reference LD.
-#'   Passed to \code{susieR::susie_rss()}.
-#' @param rMismatch LD mismatch correction method passed to \code{susieR::susie_rss()}.
-#'   Default NULL disables mismatch correction.
 #' @param keepIndel Whether to keep indel variants (default TRUE).
 #' @param commentString Comment character for sumstat file (default "#").
 #' @param diagnostics Whether to include diagnostic info (default FALSE).
@@ -547,18 +549,20 @@ rssAnalysisPipeline <- function(
     methods = NULL,
     addSusieInf = TRUE,
     finemappingOpts = list(
-      L = 20, L_greedy = 5,
-      coverage = c(0.95, 0.7, 0.5), signal_cutoff = 0.025,
-      min_abs_corr = 0.8
+      coverage = c(0.95, 0.7, 0.5), signal_cutoff = 0.025
     ),
     impute = TRUE, imputeOpts = list(rcond = 0.01, R2_threshold = 0.6, minimum_ld = 5, lamb = 0.01),
-    pipCutoffToSkip = 0, rFinite = NULL, rMismatch = NULL,
+    pipCutoffToSkip = 0,
     keepIndel = TRUE, commentString = "#", diagnostics = FALSE,
     binaryTraitModel = c("rss", "ols")) {
   binaryTraitModel <- match.arg(binaryTraitModel)
   if (!is(ldData, "LdData")) {
     stop("ldData must be an LdData object")
   }
+  # R_finite / R_mismatch are susie_rss fit options supplied via finemappingOpts;
+  # source them once here for the QC stage (which also accepts them).
+  rFinite <- finemappingOpts$R_finite %||% finemappingOpts$rFinite
+  rMismatch <- finemappingOpts$R_mismatch %||% finemappingOpts$rMismatch
   res <- list()
   rssInput <- loadRssData(
     sumstatPath = sumstatPath, columnFilePath = columnFilePath,
@@ -648,20 +652,22 @@ rssAnalysisPipeline <- function(
     priCoverage <- finemappingOpts$coverage[1]
     secCoverage <- if (length(finemappingOpts$coverage) > 1) finemappingOpts$coverage[-1] else NULL
 
-    finemappingOptsL_greedy <- finemappingOpts$L_greedy %||% finemappingOpts$lGreedy
-    finemappingOptsSignalCutoff <- finemappingOpts$signal_cutoff %||% finemappingOpts$signalCutoff
-    finemappingOptsMinAbsCorr <- finemappingOpts$min_abs_corr %||% finemappingOpts$minAbsCorr
+    finemappingOptsSignalCutoff <- finemappingOpts$signal_cutoff %||% finemappingOpts$signalCutoff %||% 0.025
+    # The fit/purity passthrough = finemappingOpts minus the pipeline-reporting keys
+    # (coverage / signal_cutoff), which susieRssPipeline takes as separate arguments.
+    # susieRssPipeline isolates min_abs_corr/median_abs_corr and forwards the rest to susie_rss.
+    finemappingOptsForFit <- finemappingOpts
+    finemappingOptsForFit$coverage <- NULL
+    finemappingOptsForFit$signal_cutoff <- NULL
+    finemappingOptsForFit$signalCutoff <- NULL
     res <- do.call(susieRssPipeline, c(susieReady, list(
-      L = finemappingOpts$L, lGreedy = finemappingOptsL_greedy,
       analysisMethod = finemappingMethod,
       methods = methods,
       addSusieInf = addSusieInf,
       coverage = priCoverage,
       secondaryCoverage = secCoverage,
       signalCutoff = finemappingOptsSignalCutoff,
-      minAbsCorr = finemappingOptsMinAbsCorr,
-      rFinite = rFinite,
-      rMismatch = rMismatch
+      finemappingOpts = finemappingOptsForFit
     )))
     if (!identical(zMismatchQc, "none") || isTRUE(alleleFlipKriging)) {
       res$outlier_number <- qcResults$outlier_number
@@ -686,15 +692,14 @@ rssAnalysisPipeline <- function(
       list(sumstats = sumstats, n = n, var_y = varY),
       ldData
     )$susie_rss_input
+    fmFit <- finemappingOpts
+    fmFit$coverage <- NULL; fmFit$signal_cutoff <- NULL; fmFit$signalCutoff <- NULL
     do.call(susieRssPipeline, c(reanalysisInput, list(
-      L = finemappingOpts$L, lGreedy = finemappingOpts$L_greedy %||% finemappingOpts$lGreedy,
       analysisMethod = method,
       coverage = priCoverage,
       secondaryCoverage = secCoverage,
       signalCutoff = finemappingOptsSignalCutoff,
-      minAbsCorr = finemappingOptsMinAbsCorr,
-      rFinite = rFinite,
-      rMismatch = rMismatch
+      finemappingOpts = fmFit
     )))
   }