--- title: "Modeling and Accounting for LD Reference Mismatch in Summary Statistics Fine-mapping" author: "Gao Wang" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Modeling and Accounting for LD Reference Mismatch in Summary Statistics Fine-mapping} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set(collapse = TRUE,comment = "#",fig.width = 5, fig.height = 3,fig.align = "center", fig.cap = " ",dpi = 120) ``` This vignette gives a practical workflow for using `susie_rss` when the GWAS summary statistics and the supplied LD from external reference samples mismatch, inducing high false positives in fine-mapping. This vignette complements the [summary statistics vignette](finemapping_summary_statistics.html) and the older [RSS diagnostic vignette](susierss_diagnostic.html). Because an effect variant induces marginal associations at correlated variants due to LD, SuSiE-RSS uses the input LD matrix to subtract the expected contribution of the fitted effect across the region at each single effect update. If that LD matrix mismatches, the subtraction is imperfect. The leftover residual can then surface as additional sparse association signals, creating spurious additional association signals. The model introduced here targets a specific mechanism by which LD mismatch affects SuSiE-RSS: when the LD matrix is estimated from a finite genotype reference panel, or from a reference panel whose population differs from the target GWAS population. The approach here is to detect this LD mismatch pattern, quantify the uncertainty it introduces, and attenuate signals that depend strongly on uncertain LD subtraction. In the current `susie_rss` implementation, `R_finite` accounts for finite LD reference panel uncertainty, and `R_mismatch = "eb"` estimates additional LD mismatch component --- mainly motivated as population drift --- from the data. A separate diagnostic flags the reliability of fine-mapped CS as well as the region by measuring signal attenuation and residual artifacts to caution users to inspect the results. ```{r} library(susieR) count_cs <- function(fit) { if (is.null(fit$sets) || is.null(fit$sets$cs)) return(0L) length(fit$sets$cs) } ``` ## Recap: a "preventable" allele coding artifact The older [RSS diagnostic vignette](susierss_diagnostic.html) discusses allele coding artifacts as source of LD mismatch --- specifically the "allele flip". We briefly revisit the same idea here because it provides a useful contrast with the our finite reference and population-drift LD mismatch model developed below. ```{r} data(SummaryConsistency, package = "susieR") if (!exists("SummaryConsistency")) load(file.path("..", "data", "SummaryConsistency.RData")) zflip <- SummaryConsistency$z Rtoy <- SummaryConsistency$ldref flip_id <- SummaryConsistency$flip_id signal_id <- SummaryConsistency$signal_id ntoy <- 10000 btoy <- numeric(length(zflip)) btoy[signal_id] <- zflip[signal_id] ``` With the allele coding problem present, ordinary `susie_rss` reports the true signal and an additional signal at the mismatched variant. ```{r} fit_toy_bad <- susie_rss(zflip, Rtoy, n = ntoy) susie_plot(fit_toy_bad, y = "PIP", b = btoy, add_legend = TRUE, main = "with allele coding artifact") points(SummaryConsistency$flip_id, fit_toy_bad$pip[flip_id], col = 7, pch = 16) ``` The kriging diagnostic for allele flip highlights the discordant variant. ```{r} kr_toy <- kriging_rss(zflip, Rtoy, n = ntoy) kr_toy$plot ``` After correcting that allele encoding, the extra credible set disappears. ```{r} zfix <- zflip zfix[flip_id] <- -zfix[flip_id] fit_toy_fixed <- susie_rss(zfix, Rtoy, n = ntoy) susie_plot(fit_toy_fixed, y = "PIP", b = btoy, add_legend = TRUE, main = "after allele correction") ``` In practice, it is useful to run this basic check before modeling LD mismatch. When there is an obvious low level allele coding artifact, it is the best practice to fix it directly through **proper bioinformatics procedures**. This can be done using many bioinformatics tools, for example `pecotmr::allele_qc` from [this package](https://github.com/StatFunGen/pecotmr). In this toy example, `susie_rss` reports `r count_cs(fit_toy_bad)` credible sets before the correction and `r count_cs(fit_toy_fixed)` after the correction, and the kriging rule identifies the deliberately misaligned variant. ## Example using real GWAS data and a 1000 Genomes European LD reference panel The example below from real data below shows a different situation where the kriging outlier check does not find isolated variants to remove, but model-based LD uncertainty still changes the fine-mapping result. The summary statistics came from a real GWAS study of a genomic region. The supplied LD matrix `R` was computed from a 1000 Genomes European LD reference panel with roughly 500 samples. The full region used in the original analysis contains more than 5,000 variants. In the full region, ordinary `susie_rss` reported 29 credible sets, whereas setting `R_finite = 500` reduced this to 3 credible sets, and adding `R_mismatch = "eb"` reduced it to 1 credible set. To keep the vignette small, the example included here uses a subset of 500 variants. The rest of this vignette walks through this smaller example step by step, showing how to use `R_finite`, `R_mismatch`, and the accompanying diagnostics to evaluate LD mismatch in practice. ```{r} data(rss_mismatch_example, package = "susieR") if (!exists("rss_mismatch_example")) load(file.path("..", "data", "rss_mismatch_example.RData")) z <- rss_mismatch_example$z bhat <- rss_mismatch_example$bhat shat <- rss_mismatch_example$shat n <- rss_mismatch_example$n R <- rss_mismatch_example$R str(rss_mismatch_example) ``` ```{r} length(z) dim(R) ``` The marginal z-scores show the association pattern before any fine-mapping model is fitted. ```{r} susie_plot(z, y = "z") ``` The marginal z-scores show one apparent signal spread over many variants. In this subset, the smallest marginal p-value is approximately $6.0 \times 10^{-4}$. We define a few small helpers for the comparisons below. ```{r} diag_scalar <- function(fit, name) { d <- fit$R_finite_diagnostics if (is.null(d) || is.null(d[[name]])) return(NA_real_) as.numeric(d[[name]][1]) } diag_flag <- function(fit, name) { d <- fit$R_finite_diagnostics if (is.null(d) || is.null(d[[name]])) return(NA) as.logical(d[[name]][1]) } fit_table <- function(fits) { data.frame( method = names(fits), n_cs = vapply(fits, count_cs, integer(1)), niter = vapply(fits, function(x) x$niter, integer(1)), converged = vapply(fits, function(x) isTRUE(x$converged), logical(1)), Q_art = vapply(fits, diag_scalar, numeric(1), name = "Q_art"), R_sensitivity_flag = vapply(fits, diag_flag, logical(1), name = "R_sensitivity_flag"), R_reliability_flag = vapply(fits, diag_flag, logical(1), name = "R_reliability_flag") ) } ``` ## Ordinary SuSiE-RSS First run `susie_rss` without any LD mismatch correction. ```{r} fit_rss <- susie_rss(bhat = bhat, shat = shat, R = R, n = n) count_cs(fit_rss) ``` ```{r} susie_plot(fit_rss, y = "PIP", add_legend = TRUE, main = "susie_rss") ``` ```{r} fit_rss$sets$purity cs_corr <- get_cs_correlation(fit_rss, Xcorr = R) round(cs_corr, 3) ``` In this subset, the ordinary fit reports `r count_cs(fit_rss)` credible sets and does not converge within `r fit_rss$niter` iterations. The reported credible sets are tight, with minimum purity at least `r round(min(fit_rss$sets$purity$min.abs.corr), 3)`, but several of them are nearly perfectly correlated with each other. This combination, many high-purity credible sets that are mutually correlated, is a typical sign of LD mismatch. ## Kriging diagnostic The older RSS diagnostic checks whether individual z-scores are surprising conditional on the other z-scores and the supplied LD matrix. This is useful for finding isolated allele coding problems. ```{r} kr <- kriging_rss(z, R, n = n) kr$plot ``` The default rule labels possible allele switches when `logLR > 2` and `abs(z) > 2`. ```{r} outlier <- which(kr$conditional_dist$logLR > 2 & abs(kr$conditional_dist$z) > 2) length(outlier) ``` For this example, the kriging diagnostic identifies `r length(outlier)` isolated allele switch outliers. This illustrates that kriging, designed for checking allele flips, does not work on this example which is not of that type of artifacts. It is not designed to model the broader finite LD reference error and mismatch between the LD reference panel and the target GWAS population that can affect an entire region. ## Finite LD reference panel correction The LD matrix here was computed from an LD reference panel of about 500 individuals. Setting `R_finite = 500` tells `susie_rss` to account for this finite LD reference panel uncertainty in each SER update. ```{r} fit_B500 <- susie_rss(bhat = bhat, shat = shat, R = R, n = n, R_finite = 500) count_cs(fit_B500) ``` ```{r} susie_plot(fit_B500, y = "PIP", add_legend = TRUE, main = "R_finite = 500") ``` ```{r} vapply(fit_B500$sets$cs, length, integer(1)) fit_B500$sets$purity cs_corr_B500 <- get_cs_correlation(fit_B500, Xcorr = R) round(cs_corr_B500, 3) ``` In this subset of 500 variants, the finite LD reference correction reduces the result from `r count_cs(fit_rss)` credible sets to `r count_cs(fit_B500)`. In the original larger region, using the same correction reduced the ordinary `susie_rss` result from 29 credible sets to 3. The credible sets also become less artificially sharp. Their sizes change from 2, 9, 9, 11 and 9 in the ordinary fit to 20, 3, 46 and 46, and the lowest purity drops from 0.9831 to 0.9496. The CS correlations are less uniformly extreme, but some remain high: one pair is still 1.000 and another is 0.978. Thus `R_finite = 500` helps, but it does not fully resolve the mismatch pattern. ## Empirical Bayes mismatch correction Finite LD reference uncertainty only accounts for the fact that the LD matrix was estimated from a finite panel. It does not account for systematic mismatch due to drift between the LD reference population and the GWAS population. The option `R_mismatch = "eb"` estimates an additional variance component shared across the region from the fitted residuals. ```{r} fit_B500_eb <- susie_rss(bhat = bhat, shat = shat, R = R, n = n, R_finite = 500, R_mismatch = "eb") count_cs(fit_B500_eb) ``` ```{r} susie_plot(fit_B500_eb, y = "PIP", add_legend = TRUE, main = "R_finite = 500, R_mismatch = eb") ``` Here, adding `R_mismatch = "eb"` reduces the result to `r count_cs(fit_B500_eb)` credible set and the fit converges. One can also run the mismatch model without specifying a finite LD reference panel size. This treats the continuous mismatch variance as the only LD uncertainty component. ```{r} fit_Binf_eb <- susie_rss(bhat = bhat, shat = shat, R = R, n = n, R_mismatch = "eb") count_cs(fit_Binf_eb) ``` ```{r} susie_plot(fit_Binf_eb, y = "PIP", add_legend = TRUE, main = "R_mismatch = eb") ``` This fit also reports `r count_cs(fit_Binf_eb)` credible set and converges quickly. ```{r} vapply(fit_B500_eb$sets$cs, length, integer(1)) fit_B500_eb$sets$purity max(fit_B500_eb$pip) ``` In this subset, the two EB mismatch fits report the same 95% credible set. It contains 34 variants, has minimum purity 0.8609, and has coverage 0.9543; the largest PIP is 0.1142. Compared with the earlier fits, the EB adjustment keeps one signal but represents it as a broad credible set rather than several small, highly correlated credible sets. This is more consistent with the weak marginal association in this region. In practice we recommend setting both `R_finite` and `R_mismatch` unless your reference LD is based on large samples such as current UK Biobank (300,000 to 500,000) in which case setting e.g. `R_finite = 300000` will have no practical impact to the results. ## Diagnostic metrics The warnings above should not be ignored. Even after the EB adjustment, `susie_rss` reports diagnostic warnings, indicating that the fine mapping result remain uncertain. The fit object returns several diagnostic metrics to summarize this uncertainty. ```{r} fits <- list(susie_rss = fit_rss, R_finite_500 = fit_B500, R_finite_500_eb = fit_B500_eb, R_mismatch_eb = fit_Binf_eb) knitr::kable(fit_table(fits), digits = 3) ``` There are two relevant reliability diagnostics. `Q_art` measures whether the final residual points into directions where the supplied LD matrix has little support. `R_sensitivity_flag` measures whether reported credible sets depend strongly on the LD uncertainty penalty. If either diagnostic is triggered, `R_reliability_flag` is `TRUE`, and the fit should be interpreted cautiously. In this subset, `R_finite = 500` is flagged by the credible set sensitivity diagnostic. The two EB mismatch fits are flagged by the residual space diagnostic, even though they report only one credible set. This is why the warning points users to `fit$R_finite_diagnostics` rather than treating the corrected fit as automatically reliable. ## Diagnostic with one effect When the reliability flag is raised, the fit object provides a diagnostic model using single-effect regression (SER): ```{r} ser <- fit_Binf_eb$R_finite_diagnostics$ser_model names(ser) ``` This is a credible set model with one effect as in Maller et al. (2012). For this model, the single row of `alpha` is the PIP. ```{r} identical(as.numeric(ser$alpha[1, ]), as.numeric(ser$pip)) count_cs(ser) vapply(ser$sets$cs, length, integer(1)) ser$sets$purity max(ser$pip) identical(ser$sets$cs, fit_Binf_eb$sets$cs) isTRUE(all.equal(as.numeric(ser$pip), as.numeric(fit_Binf_eb$pip))) ``` ```{r} susie_plot(ser, y = "PIP", add_legend = TRUE, main = "diagnostic with one effect") ``` In this example, the one-effect diagnostic gives the same result as the EB mismatch fit above: one 34-variant credible set, minimum purity 0.8609, and maximum PIP 0.1142. This exact match is merely a coincidence, not a general property of the method. When an EB mismatch fit produces a reliability warning, compare the EB result with this one-effect diagnostic rather than treating the EB fit as automatically reliable. ## Summary For routine analyses with GWAS summary statistics and an external LD reference panel: - First run basic allele and variant alignment quality control and checks. The kriging diagnostic can be useful for such obvious, localized allele coding artifacts. - Fit ordinary `susie_rss` and look for signs of instability, such as nonconvergence, multiple highly correlated credible sets, or many small credible sets in a weak association region. - Use `R_finite` when the LD matrix was estimated from a finite LD reference panel; add `R_mismatch = "eb"` when broader LD mismatch remains after accounting for finite reference uncertainty. - Inspect `fit$R_finite_diagnostics`, including `Q_art`, `R_sensitivity_flag`, and `R_reliability_flag`. - If a reliability warning is raised, compare the EB result with the SER diagnostic model, and interpret the fine-mapping result conservatively with the SER model as a safer fall-back. ## Caution: allele coding artifacts are outside this model The old RSS vignette example is based on the out-of-sample LD reference example from Zou et al. (2022). Because its LD matrix was computed from a finite reference panel of 500 individuals, it is appropriate to set `R_finite = 500` when using the new diagnostics. This setting can raise a reliability warning, but that warning is not an allele flip correction and it does not replace allele alignment QC. ```{r, message=FALSE, warning=FALSE} fit_toy_bad_B500_eb <- susie_rss(zflip, Rtoy, n = ntoy, R_finite = 500, R_mismatch = "eb", track_fit = TRUE) d_toy <- fit_toy_bad_B500_eb$R_finite_diagnostics knitr::kable(data.frame( n_cs = count_cs(fit_toy_bad_B500_eb), flipped_variant_pip = unname(fit_toy_bad_B500_eb$pip[SummaryConsistency$flip_id]), Q_art = d_toy$Q_art, residual_artifact_flag = d_toy$artifact_flag, CS_sensitivity_flag = d_toy$R_sensitivity_flag, overall_reliability_flag = d_toy$R_reliability_flag ), digits = 3, row.names = FALSE) ``` ```{r} susie_plot(fit_toy_bad_B500_eb, y = "PIP", b = btoy, add_legend = TRUE, main = "allele artifact with R_finite = 500, R_mismatch = eb") points(SummaryConsistency$flip_id, fit_toy_bad_B500_eb$pip[SummaryConsistency$flip_id], col = 7, pch = 16) ``` This fit still reports 2 credible sets, and the flipped variant still has PIP close to 1. The residual artifact diagnostic does not flag this fit: `Q_art` is small and `residual_artifact_flag` is `FALSE`. However, the CS-level sensitivity diagnostic is `TRUE`, so the overall reliability flag is also `TRUE`. The warning therefore comes from CS sensitivity, not from `Q_art`. This behavior is expected. The model in this vignette is designed for finite reference error and population drift in the LD matrix. It does not include an allele flip parameter, and it does not try to decide whether an individual z-score should have its sign reversed. That is deliberate: allele alignment is a data-processing problem that should be handled before fitting the LD mismatch model. The reason is visible from the z-scores. Given the true signal and the LD matrix, the flipped variant should have a positive marginal association. Instead its observed z-score has the opposite sign. ```{r} flip_expected <- Rtoy[SummaryConsistency$flip_id, SummaryConsistency$signal_id] * zflip[SummaryConsistency$signal_id] knitr::kable(data.frame( flipped_variant_z = zflip[SummaryConsistency$flip_id], LD_implied_z_from_signal = flip_expected, residual_after_signal = zflip[SummaryConsistency$flip_id] - flip_expected ), digits = 3) ``` This is a deterministic sign error at one variant, not the stochastic LD uncertainty modeled by `R_finite` or `R_mismatch = "eb"`. The compact track history can then be used to follow this specific artifact. We do not need to inspect the full object. For this toy example, the relevant rows are the single effects whose top variable is either the true signal or the known flipped variant. After the first IBSS iteration, one effect is already assigned almost entirely to the flipped variant. ```{r} track_toy <- fit_toy_bad_B500_eb$trace track_effect <- subset(track_toy$effect, iteration %in% 1:3 & top_variable %in% c(SummaryConsistency$signal_id, SummaryConsistency$flip_id)) track_story <- data.frame( iteration = track_effect$iteration, effect = track_effect$effect, variant = paste0(track_effect$top_variable_name, " (", ifelse( track_effect$top_variable == SummaryConsistency$signal_id, "true signal", "flipped variant" ), ")"), alpha_at_variant = track_effect$top_alpha, max_variable_lbf = track_effect$top_lbf, single_effect_lbf = track_effect$effect_lbf ) knitr::kable(track_story, digits = 3) ``` Because this is a toy example, we know which variant is the true signal and which variant was deliberately flipped. The track table itself only shows the fitted single effect summaries: one effect is already highly concentrated on the flipped variant after the first iteration. The CS-level diagnostic also shows why this happens. The fitted log BF and the log BF attenuation are different quantities: the former describes the fitted signal strength, whereas the latter measures how much that signal depends on uncertain LD. ```{r} cs_story <- track_toy$cs_sensitivity final_effect <- track_toy$effect[ track_toy$effect$iteration == max(track_toy$effect$iteration), c("effect", "top_lbf") ] cs_story <- merge(cs_story, final_effect, by = "effect", all.x = TRUE, sort = FALSE) cs_story <- cs_story[cs_story$is_sensitive, c("cs", "effect", "cs_size", "top_lbf", "cs_max_bf_attenuation", "threshold", "cs_sensitivity_label", "top_sensitive_variable_name")] names(cs_story) <- c("cs", "effect", "cs_size", "max_log_BF", "log_BF_attenuation", "threshold", "sensitivity", "sensitive_variant") knitr::kable(cs_story, digits = 3) ``` For both credible sets, the attenuation is above the threshold `log(20) = 3.0`, so both are labeled `highly_sensitive`. This is why `CS_sensitivity_flag` is `TRUE`, and why the overall reliability flag is `TRUE`, even though `Q_art = 0.012` and the residual artifact flag is `FALSE`. Thus, an allele flip can create a convincing sparse effect before the LD mismatch model has anything useful to adjust. In this example, `R_finite = 500` still raises a reliability warning through CS-level sensitivity, but that warning should be interpreted as a cue to inspect alignment, not as an automatic correction of the allele coding problem. In practice this warning is not guaranteed: for example, fitting this same toy example with `R_mismatch = "eb"` but without `R_finite = 500` does not raise the reliability warning. ```{r cleanup, include=FALSE} invisible(NULL) ``` ## Session information ```{r session-info} sessionInfo() ``` ## Reference Maller JB, McVean G, Byrnes J, Vukcevic D, Palin K, Su Z, et al. Bayesian refinement of association signals for 14 loci in 3 common diseases. Nature Genetics 44, 1294-1301 (2012).