Bayesian Fine-Mapping with BayesC

Authors

Affiliations

Palle Duun Rohde

Genomic Medicine, Department of Health Science and Technology, Aalborg University, Denmark

Peter Sørensen

Center for Quantitative Genetics and Genomics, Aarhus University, Denmark

Introduction

A genome-wide association study (GWAS) often identifies an associated region rather than a unique causal variant. Nearby SNPs are correlated through linkage disequilibrium (LD), so several SNPs can show strong marginal association even when only one SNP is truly causal.

Fine-mapping asks a more focused question:

Which SNPs in the associated region remain plausible causal explanations after accounting for LD?

In this tutorial we use BayesC as a simple Bayesian fine-mapping model. BayesC is a spike-and-slab model in which each SNP is either inactive with effect exactly zero, or active with an effect drawn from a normal distribution.

Companion presentation: Bayesian Fine-Mapping with BayesC — Slides

The workflow mirrors the accompanying slide deck:

load a small genotype region;
simulate one causal SNP;
perform marginal GWAS tests;
inspect LD among the candidate SNPs;
fit all SNPs jointly with BayesC;
calculate posterior inclusion probabilities (PIPs) and posterior mean effects;
rank candidate SNPs;
construct a 0.9 local credible set (LCS) by ranking SNPs in the local LD neighborhood and summing their raw PIPs.

Everything is implemented in base R. The only external data file required is:

../data/genotypes.txt

Load the Genotype Data

Read the genotype matrix and select the same 30-SNP region used in the slides.

Code

geno <- read.delim(
  "../data/genotypes.txt",
  header = TRUE,
  row.names = 1,
  check.names = FALSE
)

X_all <- as.matrix(geno)

region <- 114:143
X <- scale(
  X_all[, region, drop = FALSE],
  center = TRUE,
  scale = TRUE
)

n <- nrow(X)
m <- ncol(X)
snp <- colnames(X)

c(n = n, m = m)

  n   m 
503  30

The rows are individuals and the columns are SNPs. In this example the locus contains 30 candidate SNPs. Each genotype column is centered and standardized to unit sample variance, so the prior variance calculations and marker cross-products use a common scale.

Simulate One Causal SNP

We simulate a phenotype using SNP 11 as the causal variant.

Code

set.seed(123)

causal <- 11

beta_true <- rep(0, m)
beta_true[causal] <- 0.35

g <- drop(X %*% beta_true)

h2 <- 0.30
vare_true <- var(g) * (1 - h2) / h2

y <- g + rnorm(n, 0, sqrt(vare_true))

The causal SNP is:

Code

snp[causal]

[1] "rs200766023_G"

The simulation is deliberately simple: there is one true causal SNP in the region. This makes it possible to see directly whether the fine-mapping model concentrates posterior support on the correct marker or distributes support across LD-correlated alternatives.

Standard GWAS: One SNP at a Time

A standard GWAS tests each SNP separately. For SNP \(j\) we fit

\[ y_i = \mu + x_{ij}\beta_j + e_i. \]

The resulting effect estimate, standard error, and \(p\)-value describe marginal association evidence. They do not condition on the other SNPs in the region.

We calculate these quantities for all 30 SNPs.

Code

gwas_beta <- numeric(m)
gwas_se <- numeric(m)
gwas_p <- numeric(m)

for (j in seq_len(m)) {
  fit_j <- lm(y ~ X[, j])
  coef_j <- summary(fit_j)$coefficients

  gwas_beta[j] <- unname(coef_j[2, 1])
  gwas_se[j] <- unname(coef_j[2, 2])
  gwas_p[j] <- unname(coef_j[2, 4])
}

gwas_results <- data.frame(
  marker = seq_len(m),
  SNP = snp,
  beta_hat = gwas_beta,
  SE = gwas_se,
  p_value = gwas_p,
  causal = seq_len(m) == causal
)

gwas_results[order(gwas_results$p_value), ]

   marker           SNP    beta_hat         SE      p_value causal
11     11 rs200766023_G  0.34990887 0.02322565 1.420561e-42   TRUE
8       8 rs199676333_C -0.30211540 0.02452817 1.223803e-30  FALSE
2       2 rs143764217_C  0.29621070 0.02467150 2.311365e-29  FALSE
14     14 rs149771619_A  0.27890547 0.02507087 7.645056e-26  FALSE
13     13 rs372428423_A  0.27785589 0.02509412 1.220905e-25  FALSE
17     17 rs199985493_T  0.27289436 0.02520256 1.077893e-24  FALSE
4       4  rs77714226_T -0.27012336 0.02526207 3.549144e-24  FALSE
1       1 rs144666531_G  0.26658638 0.02533696 1.583992e-23  FALSE
3       3 rs138139885_C  0.26573288 0.02535485 2.263028e-23  FALSE
7       7   rs2716416_C  0.26569930 0.02535555 2.294936e-23  FALSE
26     26  rs78913776_T  0.26081248 0.02545662 1.714950e-22  FALSE
27     27 rs371582893_A  0.25955620 0.02548224 2.851915e-22  FALSE
23     23 rs148271204_G  0.25894083 0.02549473 3.654194e-22  FALSE
20     20 rs201317463_T  0.25210072 0.02563120 5.441767e-21  FALSE
12     12 rs201631627_T -0.24130531 0.02583777 3.168011e-19  FALSE
29     29 rs114833654_T  0.23125685 0.02602054 1.128395e-17  FALSE
25     25 rs146675715_T  0.23106385 0.02602396 1.206224e-17  FALSE
10     10 rs201331129_A  0.22066305 0.02620352 3.953248e-16  FALSE
18     18 rs201146901_C -0.21368041 0.02631878 3.676519e-15  FALSE
21     21 rs199697997_T -0.20244187 0.02649552 1.107620e-13  FALSE
19     19 rs201944730_T -0.19903847 0.02654694 2.973910e-13  FALSE
28     28 rs370880731_C -0.15890098 0.02708166 8.050312e-09  FALSE
16     16 rs370414481_A  0.09326122 0.02768490 8.135171e-04  FALSE
6       6   rs2007911_C  0.08119491 0.02776068 3.602798e-03  FALSE
9       9 rs368834055_A  0.07647489 0.02778742 6.135385e-03  FALSE
5       5 rs368129211_G  0.05397037 0.02789266 5.356129e-02  FALSE
15     15  rs77287660_A -0.04378200 0.02792827 1.175922e-01  FALSE
22     22 rs370167421_C -0.02979486 0.02796502 2.871932e-01  FALSE
24     24   rs2212121_T -0.02574751 0.02797304 3.577857e-01  FALSE
30     30 rs201836474_C  0.00582386 0.02799548 8.352912e-01  FALSE

Local Manhattan Plot

A local Manhattan-style plot shows the marginal GWAS signal across the locus.

Code

plot(
  seq_len(m),
  -log10(gwas_p),
  type = "h",
  lwd = 4,
  xlab = "SNP position in the local region",
  ylab = expression(-log[10](p)),
  main = "Marginal GWAS signal in the locus"
)

points(
  causal,
  -log10(gwas_p[causal]),
  pch = 19,
  cex = 1.4,
  col = "red"
)

abline(v = causal, col = "red", lty = 2)

Several nearby SNPs can look associated because they tag the same underlying signal through LD. The GWAS therefore identifies an associated region, but it does not necessarily identify a unique causal variant.

Inspect Linkage Disequilibrium

The LD correlation matrix for the candidate SNPs is

\[ R = \operatorname{cor}(X). \]

Code

R <- cor(X)

image(
  R,
  axes = FALSE,
  main = "LD correlation matrix"
)

axis(
  1,
  at = seq(0, 1, length.out = m),
  labels = seq_len(m),
  cex.axis = 0.7
)

axis(
  2,
  at = seq(0, 1, length.out = m),
  labels = seq_len(m),
  cex.axis = 0.7
)

We can also inspect LD with the true causal SNP directly.

Code

plot(
  seq_len(m),
  R[causal, ]^2,
  type = "h",
  lwd = 3,
  xlab = "SNP",
  ylab = expression(r^2~"with causal SNP"),
  main = "LD around the causal SNP"
)

abline(v = causal, col = "red", lty = 2)

A non-causal SNP can have a strong marginal GWAS signal simply because it is strongly correlated with the causal SNP.

Fine-Mapping: Inputs and Outputs

The main inputs are:

marker-level association information, represented here by the phenotype and genotype data and summarized by \(\hat\beta_j\) and \(\mathrm{SE}(\hat\beta_j)\);
LD among candidate SNPs, represented by \(R\) or equivalently by the genotype cross-product structure;
a sparse joint model, here BayesC.

The main outputs are:

posterior inclusion probabilities (PIPs);
posterior mean SNP effects;
ranked candidate variants;
a local credible set summarizing the strongest signal and its LD-related alternatives.

The key difference from GWAS is that fine-mapping evaluates the SNPs jointly.

BayesC Fine-Mapping Model

At the individual-data level we use

\[ y = \mu + Xb + e, \qquad e \sim \mathcal{N}(0,v_e I_N). \]

Here:

\(X\) is the genotype matrix;
\(b=(b_1,\ldots,b_M)^\top\) is the vector of SNP effects;
\(v_e\) is the residual variance.

Because all SNPs are fitted together, LD enters the joint likelihood through the genotype cross-product structure \(X^\top X\).

This is the important conceptual transition from GWAS to fine-mapping:

\[ \text{marginal association} \quad\longrightarrow\quad \text{joint posterior evidence}. \]

Correlated SNPs can now compete to explain the same association signal.

BayesC Prior

BayesC assumes that each SNP is either inactive or active. Introduce an indicator

\[ d_j \in \{0,1\}. \]

The prior is

\[ d_j \sim \operatorname{Bernoulli}(\pi), \]

with

\[ b_j = 0 \qquad\text{if }d_j=0, \]

and

\[ b_j\mid d_j=1 \sim \mathcal{N}(0,v_b). \]

Equivalently,

\[ b_j \sim (1-\pi)\delta_0 + \pi \mathcal{N}(0,v_b). \]

We use

Code

pi_active <- 0.03

c(
  prior_P_active = pi_active,
  prior_expected_active_SNPs = m * pi_active
)

            prior_P_active prior_expected_active_SNPs 
                      0.03                       0.90

Thus each SNP has prior probability 0.03 of being active, corresponding to a prior expectation of 0.9 active SNPs among the 30 markers.

For this tutorial we keep \(\pi\) fixed. The goal is to study the local fine-mapping problem, not to estimate genome-wide architecture parameters from only 30 SNPs.

What the Gibbs Sampler Does

At each MCMC iteration the sampler updates every SNP conditional on the current values of the other SNP effects.

For SNP \(j\) it considers two possibilities:

\[ d_j=0 \quad\Rightarrow\quad b_j=0, \]

\[ d_j=1 \quad\Rightarrow\quad b_j\neq0. \]

Because the residual signal changes whenever another SNP is included or excluded, the inclusion probability for SNP \(j\) automatically reflects competition with LD-correlated SNPs.

The sampler also updates the non-zero effect variance and residual variance.

Base-R BayesC Gibbs Sampler

The following implementation is intentionally compact and teaching-oriented. It uses only base R.

Code

fit_bayesc <- function(
  X,
  y,
  niter = 2000,
  pi = 0.03,
  start_h2 = 0.30
) {

  n <- nrow(X)
  m <- ncol(X)

  stopifnot(
    pi > 0,
    pi < 1,
    start_h2 > 0,
    start_h2 < 1
  )

  # ------------------------------------------------------------
  # Starting values
  # ------------------------------------------------------------
  vary <- var(y)
  varg <- vary * start_h2
  vare <- vary * (1 - start_h2)

  # With one non-zero component, E(number active) = m * pi.
  varb <- varg / (m * pi)

  # Scaled inverse-chi-square priors for varb and vare.
  nub <- 4
  nue <- 4

  scale_b <- (nub - 2) / nub * varb
  scale_e <- (nue - 2) / nue * vare

  # ------------------------------------------------------------
  # State variables
  # ------------------------------------------------------------
  beta <- numeric(m)
  delta <- integer(m)
  mu <- mean(y)

  # ycorr is the current residual y - mu - X beta.
  ycorr <- y - mu
  xpx <- colSums(X^2)

  beta_mcmc <- matrix(
    0,
    nrow = niter,
    ncol = m
  )

  delta_mcmc <- matrix(
    0L,
    nrow = niter,
    ncol = m
  )

  par_mcmc <- matrix(
    NA_real_,
    nrow = niter,
    ncol = 5
  )

  colnames(par_mcmc) <- c(
    "Nactive",
    "Varb",
    "Vare",
    "Varg",
    "h2"
  )

  # ------------------------------------------------------------
  # Gibbs sampler
  # ------------------------------------------------------------
  for (iter in seq_len(niter)) {

    # 1. Update the intercept.
    ycorr <- ycorr + mu

    mu_hat <- mean(ycorr)
    mu <- rnorm(
      1,
      mean = mu_hat,
      sd = sqrt(vare / n)
    )

    ycorr <- ycorr - mu

    # 2. Update SNPs one at a time.
    for (j in seq_len(m)) {

      xj <- X[, j]
      old_beta <- beta[j]

      # Add the current SNP effect back before updating it.
      residual_without_j <- ycorr + xj * old_beta

      rhs <- sum(xj * residual_without_j)

      # Posterior distribution if the SNP is active.
      Cj <- xpx[j] + vare / varb
      post_mean <- rhs / Cj
      post_var <- vare / Cj

      # Integrated evidence for the active component.
      log_prob_zero <- log(1 - pi)

      log_prob_active <-
        log(pi) +
        0.5 * (
          log(vare) -
          log(varb) -
          log(Cj) +
          rhs^2 / (vare * Cj)
        )

      max_log <- max(
        log_prob_zero,
        log_prob_active
      )

      prob <- exp(
        c(
          log_prob_zero,
          log_prob_active
        ) - max_log
      )

      prob <- prob / sum(prob)

      delta[j] <- sample(
        c(0L, 1L),
        size = 1,
        prob = prob
      )

      if (delta[j] == 0L) {
        beta[j] <- 0
      } else {
        beta[j] <- rnorm(
          1,
          mean = post_mean,
          sd = sqrt(post_var)
        )
      }

      # Update residual after the new SNP effect is sampled.
      ycorr <- residual_without_j - xj * beta[j]
    }

    # 3. Update the variance of active SNP effects.
    active <- which(delta == 1L)
    nactive <- length(active)

    shape_b <- 0.5 * (nub + nactive)
    rate_b <- 0.5 * nub * scale_b

    if (nactive > 0) {
      rate_b <-
        rate_b +
        0.5 * sum(beta[active]^2)
    }

    varb <- 1 / rgamma(
      1,
      shape = shape_b,
      rate = rate_b
    )

    # 4. Update residual variance.
    resid <- y - mu - drop(X %*% beta)

    shape_e <- 0.5 * (nue + n)
    rate_e <-
      0.5 * (
        nue * scale_e +
        sum(resid^2)
      )

    vare <- 1 / rgamma(
      1,
      shape = shape_e,
      rate = rate_e
    )

    # 5. Derived quantities and storage.
    varg_iter <- var(drop(X %*% beta))
    h2_iter <- varg_iter / (varg_iter + vare)

    beta_mcmc[iter, ] <- beta
    delta_mcmc[iter, ] <- delta

    par_mcmc[iter, ] <- c(
      nactive,
      varb,
      vare,
      varg_iter,
      h2_iter
    )
  }

  colnames(beta_mcmc) <- colnames(X)
  colnames(delta_mcmc) <- colnames(X)

  list(
    beta = beta_mcmc,
    delta = delta_mcmc,
    par = par_mcmc,
    input = list(
      niter = niter,
      pi = pi,
      start_h2 = start_h2
    )
  )
}

The important part for fine-mapping is the sampled inclusion indicator \(d_j\). A SNP that is repeatedly required to explain the signal will be active in a large fraction of posterior samples.

Fit BayesC

Use the same prior probability as above and keep the first 500 iterations as burn-in.

Code

set.seed(1)

niter <- 2000
burn_in <- 500

fit <- fit_bayesc(
  X = X,
  y = y,
  niter = niter,
  pi = pi_active,
  start_h2 = 0.30
)

keep <- (burn_in + 1):niter

Posterior Inclusion Probabilities

The posterior inclusion probability for SNP \(j\) is

\[ \operatorname{PIP}_j = P(d_j=1\mid\mathcal D). \]

From the MCMC samples we estimate it as

\[ \widehat{\operatorname{PIP}}_j = \frac{1}{T} \sum_{t=1}^T d_j^{(t)}. \]

Code

PIP <- colMeans(
  fit$delta[keep, , drop = FALSE]
)

The PIP comes from the joint BayesC model. SNP \(j\) is updated conditional on the current values of the other SNPs, and the final PIP averages over uncertainty in those other posterior states.

This is sometimes called a marginal PIP because the other model states have been marginalized out. It should not be confused with a marginal GWAS analysis: the underlying fit is joint.

Posterior Mean Effects

The posterior mean SNP effect is

\[ E(b_j\mid\mathcal D), \]

estimated by averaging the sampled SNP effects over all retained iterations, including the iterations in which the SNP is inactive.

Code

beta_post <- colMeans(
  fit$beta[keep, , drop = FALSE]
)

PIP and posterior mean effect answer different questions:

PIP: how strongly does the posterior support inclusion of this SNP?
posterior mean effect: what is the average effect after model averaging over inclusion and exclusion?

Fine-Mapping Results

Combine the marginal GWAS evidence, LD, true simulation information, and BayesC posterior summaries.

Code

results <- data.frame(
  marker = seq_len(m),
  SNP = snp,
  marginal_beta = gwas_beta,
  marginal_p = gwas_p,
  r2_with_causal = R[causal, ]^2,
  true_effect = beta_true,
  PIP = PIP,
  posterior_mean_effect = beta_post,
  causal = seq_len(m) == causal
)

ranked_results <-
  results[
    order(-results$PIP),
    ,
    drop = FALSE
  ]

head(ranked_results, 10)

              marker           SNP marginal_beta   marginal_p r2_with_causal
rs200766023_G     11 rs200766023_G    0.34990887 1.420561e-42     1.00000000
rs77714226_T       4  rs77714226_T   -0.27012336 3.549144e-24     0.74015325
rs2716416_C        7   rs2716416_C    0.26569930 2.294936e-23     0.69876803
rs143764217_C      2 rs143764217_C    0.29621070 2.311365e-29     0.78917002
rs146675715_T     25 rs146675715_T    0.23106385 1.206224e-17     0.49484714
rs149771619_A     14 rs149771619_A    0.27890547 7.645056e-26     0.61846167
rs199985493_T     17 rs199985493_T    0.27289436 1.077893e-24     0.54557553
rs368129211_G      5 rs368129211_G    0.05397037 5.356129e-02     0.05815418
rs370880731_C     28 rs370880731_C   -0.15890098 8.050312e-09     0.13961262
rs148271204_G     23 rs148271204_G    0.25894083 3.654194e-22     0.53566334
              true_effect         PIP posterior_mean_effect causal
rs200766023_G        0.35 1.000000000          0.3607370502   TRUE
rs77714226_T         0.00 0.116666667          0.0130069629  FALSE
rs2716416_C          0.00 0.024666667         -0.0019739998  FALSE
rs143764217_C        0.00 0.016000000         -0.0004875785  FALSE
rs146675715_T        0.00 0.011333333         -0.0006590190  FALSE
rs149771619_A        0.00 0.010000000          0.0006546902  FALSE
rs199985493_T        0.00 0.007333333          0.0002872568  FALSE
rs368129211_G        0.00 0.006666667         -0.0001970794  FALSE
rs370880731_C        0.00 0.006666667         -0.0001951910  FALSE
rs148271204_G        0.00 0.006000000          0.0001140686  FALSE

The important comparison is between the marginal GWAS ranking and the posterior fine-mapping ranking. Several SNPs may show strong marginal association, while the joint BayesC analysis can concentrate posterior support more strongly on the SNPs that best explain the signal after accounting for LD.

Plot the PIPs

Code

plot(
  seq_len(m),
  PIP,
  type = "h",
  lwd = 5,
  ylim = c(0, 1),
  xlab = "SNP index in the 30-marker region",
  ylab = "Posterior inclusion probability",
  main = "BayesC fine-mapping result"
)

points(
  causal,
  PIP[causal],
  pch = 19,
  cex = 1.4,
  col = "red"
)

abline(v = causal, col = "red", lty = 2)

legend(
  "topright",
  legend = c("BayesC PIP", "causal SNP"),
  col = c("black", "red"),
  lwd = c(5, 1),
  lty = c(1, 2),
  pch = c(NA, 19),
  bty = "n"
)

The causal SNP should receive strong posterior support. LD proxies can still retain non-zero PIP because they provide alternative explanations of the same association signal.

Quantify Fine-Mapping Uncertainty

Because this is a simulation, we know which SNP is causal and can calculate a few useful summaries.

Code

noncausal <- setdiff(
  seq_len(m),
  causal
)

best_noncausal <-
  noncausal[
    which.max(PIP[noncausal])
  ]

causal_rank <-
  rank(
    -PIP,
    ties.method = "min"
  )[causal]

finemap_summary <- data.frame(
  quantity = c(
    "Causal SNP PIP",
    "Causal SNP rank",
    "Best non-causal PIP",
    "PIP gap",
    "Posterior expected active SNPs",
    "Best competitor marker",
    "Best competitor SNP"
  ),
  value = c(
    round(PIP[causal], 3),
    causal_rank,
    round(PIP[best_noncausal], 3),
    round(
      PIP[causal] -
      PIP[best_noncausal],
      3
    ),
    round(sum(PIP), 3),
    best_noncausal,
    snp[best_noncausal]
  )
)

finemap_summary

                        quantity        value
1                 Causal SNP PIP            1
2                Causal SNP rank            1
3            Best non-causal PIP        0.117
4                        PIP gap        0.883
5 Posterior expected active SNPs        1.277
6         Best competitor marker            4
7            Best competitor SNP rs77714226_T

The summaries have direct interpretations:

causal SNP PIP: posterior support for the true causal SNP;
causal rank: whether the true causal SNP is the leading candidate;
best non-causal PIP: support retained by the strongest alternative;
PIP gap: separation between the causal SNP and its strongest competitor;
sum of PIPs: posterior expected number of active SNPs in the region.

The last identity follows because

\[ E(N_{\mathrm{active}}\mid\mathcal D) = E\left(\sum_j d_j\mid\mathcal D\right) = \sum_j \operatorname{PIP}_j. \]

PIP and LD

A useful diagnostic is to compare PIP with LD to the causal SNP.

Code

plot(
  R[causal, ]^2,
  PIP,
  pch = 19,
  xlab = expression(r^2~"with causal SNP"),
  ylab = "BayesC PIP",
  main = "LD and posterior inclusion probability"
)

points(
  R[causal, causal]^2,
  PIP[causal],
  pch = 19,
  col = "red",
  cex = 1.5
)

Strong LD does not automatically imply high PIP. Instead, it tells us which SNPs contain similar information. The joint posterior determines how much support each of those correlated alternatives retains.

PIP Is Not the Same as Effect Size

Code

plot(
  PIP,
  abs(beta_post),
  pch = 19,
  xlab = "Posterior inclusion probability",
  ylab = "Absolute posterior mean effect",
  main = "PIP versus posterior mean effect"
)

points(
  PIP[causal],
  abs(beta_post[causal]),
  pch = 19,
  col = "red",
  cex = 1.5
)

A SNP can have appreciable PIP but a modest posterior mean effect if it is active only in part of the posterior or if its active effects are small. Conversely, a SNP selected less frequently can have a relatively large effect when it is selected.

For fine-mapping, PIP is the main marker-level measure of support for inclusion.

From PIPs to a Local Credible Set

We now use the same local credible-set construction as in the slide deck.

The procedure is:

identify the SNP with the highest PIP;
define the local LD neighborhood around that focal SNP, here using \(r^2\ge0.5\);
rank SNPs in that neighborhood by decreasing raw PIP;
add PIPs until the cumulative sum reaches 0.90.

Thus the 0.9 local credible set \(\mathcal C\) satisfies

\[ \sum_{j\in\mathcal C} \operatorname{PIP}_j \ge 0.90. \]

We do not normalize the PIPs by their total across the locus.

Because BayesC permits several SNPs to be active, this cumulative-raw-PIP LCS is a practical localization summary. Without a single-causal-variant assumption or a configuration-based calculation, it is not a formal \(90\%\) posterior coverage set.

Construct the 0.9 Local Credible Set

Code

make_local_credible_set <- function(
  PIP,
  R,
  SNP,
  coverage = 0.90,
  r2_threshold = 0.50
) {

  # SNP with highest posterior inclusion probability.
  focal <- which.max(PIP)

  # Local LD neighborhood of the focal SNP.
  local <- which(
    R[focal, ]^2 >= r2_threshold
  )

  # Rank local candidates by raw PIP.
  ord <-
    local[
      order(
        PIP[local],
        decreasing = TRUE
      )
    ]

  cumulative_PIP <- cumsum(PIP[ord])

  # Smallest number of ranked SNPs reaching the target.
  if (any(cumulative_PIP >= coverage)) {
    n_keep <- which(
      cumulative_PIP >= coverage
    )[1]
  } else {
    n_keep <- length(ord)
    warning(
      "Cumulative PIP in the local LD neighborhood does not reach the requested coverage."
    )
  }

  in_set <- seq_along(ord) <= n_keep

  table <- data.frame(
    rank = seq_along(ord),
    marker = ord,
    SNP = SNP[ord],
    PIP = PIP[ord],
    r2_with_focal = R[focal, ord]^2,
    cumulative_PIP = cumulative_PIP,
    in_LCS = in_set
  )

  list(
    focal = focal,
    local_candidates = local,
    LCS = ord[in_set],
    table = table,
    coverage = coverage,
    r2_threshold = r2_threshold
  )
}

lcs <- make_local_credible_set(
  PIP = PIP,
  R = R,
  SNP = snp,
  coverage = 0.90,
  r2_threshold = 0.50
)

Inspect the focal SNP and the local candidate set.

Code

c(
  focal_marker = lcs$focal,
  focal_SNP = snp[lcs$focal],
  number_local_candidates = length(lcs$local_candidates),
  LCS_size = length(lcs$LCS),
  causal_in_LCS = causal %in% lcs$LCS
)

focal_marker.rs200766023_G                  focal_SNP 
                      "11"            "rs200766023_G" 
   number_local_candidates                   LCS_size 
                      "13"                        "1" 
             causal_in_LCS 
                    "TRUE"

Inspect the ranked local candidates.

Code

lcs$table

              rank marker           SNP         PIP r2_with_focal
rs200766023_G    1     11 rs200766023_G 1.000000000     1.0000000
rs77714226_T     2      4  rs77714226_T 0.116666667     0.7401532
rs2716416_C      3      7   rs2716416_C 0.024666667     0.6987680
rs143764217_C    4      2 rs143764217_C 0.016000000     0.7891700
rs149771619_A    5     14 rs149771619_A 0.010000000     0.6184617
rs199985493_T    6     17 rs199985493_T 0.007333333     0.5455755
rs148271204_G    7     23 rs148271204_G 0.006000000     0.5356633
rs371582893_A    8     27 rs371582893_A 0.005333333     0.5007914
rs199676333_C    9      8 rs199676333_C 0.004666667     0.8156835
rs138139885_C   10      3 rs138139885_C 0.004000000     0.6066870
rs78913776_T    11     26  rs78913776_T 0.004000000     0.5434041
rs144666531_G   12      1 rs144666531_G 0.002666667     0.5671530
rs372428423_A   13     13 rs372428423_A 0.002000000     0.6135135
              cumulative_PIP in_LCS
rs200766023_G       1.000000   TRUE
rs77714226_T        1.116667  FALSE
rs2716416_C         1.141333  FALSE
rs143764217_C       1.157333  FALSE
rs149771619_A       1.167333  FALSE
rs199985493_T       1.174667  FALSE
rs148271204_G       1.180667  FALSE
rs371582893_A       1.186000  FALSE
rs199676333_C       1.190667  FALSE
rs138139885_C       1.194667  FALSE
rs78913776_T        1.198667  FALSE
rs144666531_G       1.201333  FALSE
rs372428423_A       1.203333  FALSE

The rows with in_LCS = TRUE form the 0.9 local credible set.

If the leading SNP itself has PIP above 0.90, the 0.9 LCS contains only that SNP. If posterior support is distributed among several LD-correlated variants, several SNPs will be needed before the cumulative PIP reaches 0.90.

Visualize Construction of the LCS

Code

plot(
  lcs$table$rank,
  lcs$table$cumulative_PIP,
  type = "b",
  pch = 19,
  xlab = "Rank within local LD neighborhood",
  ylab = "Cumulative PIP",
  ylim = c(0, 1),
  main = "Building the 0.9 local credible set"
)

abline(
  h = 0.90,
  lty = 2,
  col = "red"
)

abline(
  v = length(lcs$LCS),
  lty = 2,
  col = "blue"
)

The LCS converts the marker-by-marker PIP output into an interpretable local candidate set. A small LCS indicates that posterior support is strongly concentrated, whereas a larger LCS reflects greater uncertainty among LD-related alternatives.

This is a local, signal-oriented summary. BayesC itself still allows more than one SNP to be active in a posterior sample.

Compare GWAS and Fine-Mapping Rankings

It is useful to compare the top SNPs under marginal GWAS evidence and BayesC PIP.

Code

top_gwas <- order(gwas_p)[1:5]
top_pip <- order(PIP, decreasing = TRUE)[1:5]

comparison <- data.frame(
  rank = 1:5,
  GWAS_marker = top_gwas,
  GWAS_SNP = snp[top_gwas],
  GWAS_p = gwas_p[top_gwas],
  BayesC_marker = top_pip,
  BayesC_SNP = snp[top_pip],
  BayesC_PIP = PIP[top_pip]
)

comparison

              rank GWAS_marker      GWAS_SNP       GWAS_p BayesC_marker
rs200766023_G    1          11 rs200766023_G 1.420561e-42            11
rs77714226_T     2           8 rs199676333_C 1.223803e-30             4
rs2716416_C      3           2 rs143764217_C 2.311365e-29             7
rs143764217_C    4          14 rs149771619_A 7.645056e-26             2
rs146675715_T    5          13 rs372428423_A 1.220905e-25            25
                 BayesC_SNP BayesC_PIP
rs200766023_G rs200766023_G 1.00000000
rs77714226_T   rs77714226_T 0.11666667
rs2716416_C     rs2716416_C 0.02466667
rs143764217_C rs143764217_C 0.01600000
rs146675715_T rs146675715_T 0.01133333

The GWAS ranking is based on one-SNP-at-a-time tests. The BayesC ranking is based on a joint posterior model in which correlated SNPs compete to explain the phenotype.

MCMC Diagnostics

Fine-mapping summaries should be based on a sampler that has reached a stable posterior region. The following simple trace plots are useful teaching diagnostics.

Number of Active SNPs

Code

plot(
  fit$par[, "Nactive"],
  type = "l",
  xlab = "Iteration",
  ylab = "Number of active SNPs",
  main = "Trace plot: active SNPs"
)

abline(
  v = burn_in,
  lty = 2,
  col = "red"
)

The sampler can move among models containing different numbers of active SNPs. The posterior average can be compared with the sum of PIPs.

Code

c(
  mean_sampled_active =
    mean(
      fit$par[keep, "Nactive"]
    ),
  sum_PIP = sum(PIP)
)

mean_sampled_active             sum_PIP 
           1.276667            1.276667

These two quantities should agree up to Monte Carlo error.

Active-Effect Variance

Code

plot(
  fit$par[, "Varb"],
  type = "l",
  xlab = "Iteration",
  ylab = expression(v[b]),
  main = "Trace plot: active-effect variance"
)

abline(
  v = burn_in,
  lty = 2,
  col = "red"
)

Heritability

Code

plot(
  fit$par[, "h2"],
  type = "l",
  xlab = "Iteration",
  ylab = expression(h^2),
  main = "Trace plot: sampled heritability"
)

abline(
  h = h2,
  lty = 2,
  col = "blue"
)

abline(
  v = burn_in,
  lty = 2,
  col = "red"
)

Code

mean(fit$par[keep, "h2"])

[1] 0.3100284

The simulated value is \(h^2=0.30\). The posterior trace should fluctuate around a stable region rather than showing a persistent trend.

What Have We Learned?

This tutorial follows the complete fine-mapping workflow from association evidence to posterior candidate variants.

The main points are:

standard GWAS tests SNPs one at a time and therefore identifies associated markers, not necessarily causal variants;
LD causes several SNPs to carry similar association information;
BayesC fits the SNPs jointly under a sparse spike-and-slab prior;
each SNP is either inactive or active in each posterior sample;
the posterior inclusion probability

\[ \operatorname{PIP}_j = P(d_j=1\mid\mathcal D) \]

summarizes how strongly the joint posterior supports inclusion of SNP \(j\);

the posterior mean effect and PIP answer different questions;
correlated SNPs can retain non-zero PIP because they provide alternative explanations of the same signal;
\(\sum_j\operatorname{PIP}_j\) is the posterior expected number of active SNPs;
candidate SNPs can be ranked directly by PIP;
a 0.9 local credible set can be constructed by defining an LD neighborhood around the leading SNP, ranking those SNPs by raw PIP, and summing PIPs until the cumulative value reaches 0.90;
no normalization of the PIPs is used in this LCS construction.

The overall workflow is therefore

\[ \begin{aligned} \text{GWAS evidence + LD} &\longrightarrow \text{BayesC joint model} \\ &\longrightarrow \text{posterior samples} \\ &\longrightarrow \text{PIPs, effects, and local credible set}. \end{aligned} \]

This provides the foundation for the next tutorial, where prior inclusion probabilities can be informed by functional annotations.

Introduction

Load the Genotype Data

Simulate One Causal SNP

Standard GWAS: One SNP at a Time

Local Manhattan Plot

Inspect Linkage Disequilibrium

Fine-Mapping: Inputs and Outputs

BayesC Fine-Mapping Model

BayesC Prior

What the Gibbs Sampler Does

Base-R BayesC Gibbs Sampler

Fit BayesC

Posterior Inclusion Probabilities

Posterior Mean Effects

Fine-Mapping Results

Plot the PIPs

Quantify Fine-Mapping Uncertainty

PIP and LD

PIP Is Not the Same as Effect Size

From PIPs to a Local Credible Set

Construct the 0.9 Local Credible Set

Visualize Construction of the LCS

Compare GWAS and Fine-Mapping Rankings

MCMC Diagnostics

Number of Active SNPs

Active-Effect Variance

Heritability

What Have We Learned?

Further Reading