---
title: "Colocalization Practical: From SuSiE to coloc.abf and coloc.susie"
author: "Masahiro Kanai"
date: 'Updated: `r Sys.Date()`'
output:
  rmarkdown::html_document:
    code_folding: show
    toc: true
    toc_float:
      collapsed: false
    df_print: paged
    number_sections: false
---

```{r init, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(eval = TRUE)
knitr::opts_chunk$set(include = TRUE)
knitr::opts_chunk$set(warning = FALSE)
knitr::opts_chunk$set(message = FALSE)
knitr::opts_chunk$set(fig.width = 8, fig.height = 3.5)
```

# Background

This practical follows directly from the fine-mapping practical. There, you simulated a
locus, ran a GWAS, and fine-mapped it with **SuSiE**, recovering credible sets of putative
causal variants. Here we take a second step. We ask whether a genetic signal at the same
locus is **shared** between two traits, a question known as **colocalization**.

The motivating biology is simple. A GWAS variant is most often non-coding, and the most
proximal molecular trait it can affect is gene expression. If a GWAS signal and an
expression QTL (eQTL) signal are driven by the *same* causal variant, the eQTL nominates a
candidate gene and direction of effect for the GWAS locus. If they merely sit nearby but are
driven by *distinct* variants in linkage disequilibrium (LD), the eQTL is a red herring.
Colocalization is the formal test that separates these two cases.

## The five hypotheses

`coloc` casts colocalization as a Bayesian model comparison over five mutually exclusive
hypotheses for a locus tested in two studies.

| Hypothesis | Meaning |
|:-:|:--|
| **H0** | No association with either trait |
| **H1** | Association with trait 1 only |
| **H2** | Association with trait 2 only |
| **H3** | Association with both traits, **distinct** causal variants |
| **H4** | Association with both traits, **one shared** causal variant |

Each hypothesis receives a posterior probability (`PP.H0` ... `PP.H4`) that sums to one.
A large **PP.H4** is evidence of colocalization. A large **PP.H3** means both traits are
associated but through different variants.

## Two flavors of coloc

- **`coloc.abf`** needs only marginal summary statistics and no LD matrix. That is its real
  advantage. Use it as the fallback **when you lack a reliable in-sample LD matrix**.
- **`coloc.abf` assumes at most one causal variant per locus.** Treat this as a limitation to
  work around, not a property you can usually assume holds. Choosing `coloc.abf` because you
  believe the locus has a single signal is the wrong reason to choose it.
- **`coloc.susie`** relaxes that assumption to **multiple causal variants per locus**. It
  takes the SuSiE fine-mapping output for each trait and tests every pair of credible sets,
  reporting a posterior probability of a shared causal variant (`PP.H4`) for each pair.
  Whenever you have a trustworthy in-sample LD matrix, this is the method to prefer.

**Note.** Real loci frequently harbor multiple causal variants, and their causal effects can
be positively or negatively correlated through LD. A second signal can then distort, or even
mask, the marginal association of the first, which breaks the single-causal-variant
assumption (Kanai et al. 2021; Zhang et al. 2023).

We will run both on the same locus and see why the distinction matters.

***

# Setup

```{r libs}
library(susieR)  # SuSiE fine-mapping
library(coloc)   # colocalization
```

***

# Recap: the fine-mapping locus and trait-1 SuSiE

We start from the locus built in the fine-mapping practical, loading its frozen genotype
matrix `ISG_2026_genotypes.rds`. We use a saved matrix because `MASS::mvrnorm` performs
differently across platforms even with the same seed, so regenerating it live would not be
reproducible. The commented chunk below shows how the file was generated once.

```{r recap-geno}
p <- 100   # variants in the locus
N <- 1000  # individuals

X <- readRDS("ISG_2026_genotypes.rds")  # frozen 1000 x 100 genotype matrix
in_sample_LD <- cov(scale(X))           # in-sample LD used by SuSiE
```

```{r make-geno, eval=FALSE}
# How ISG_2026_genotypes.rds was generated:
# library(MASS)
# LD <- matrix(0.3, nrow = p, ncol = p)
# diag(LD) <- 1.0
# set.seed(4)
# X <- mvrnorm(n = N, mu = rep(0, p), Sigma = LD)
# saveRDS(X, "ISG_2026_genotypes.rds")
```

```{r recap-pheno}
# Trait 1: a GWAS phenotype with L = 3 causal variants
L <- 3
h2g <- 0.1

set.seed(1)
causal_config <- rep(0, p)
causal_ind1 <- sample(1:p, L, replace = FALSE)   # the true causal variants
causal_config[causal_ind1] <- 1

per_snp_h2g <- h2g / L
effect_sizes <- rnorm(L, mean = 0, sd = sqrt(per_snp_h2g))
beta1 <- causal_config
beta1[causal_ind1] <- effect_sizes

genetic_effect <- X %*% beta1
var_g <- var(genetic_effect)
sigma_squared <- ifelse(1 - var_g > 0.05, 1 - var_g, 0.1)  # environmental noise
y1 <- as.numeric(genetic_effect + rnorm(N, mean = 0, sd = sqrt(sigma_squared)))

cat(sprintf("Trait 1 causal variants: %s\n", paste(causal_ind1, collapse = ", ")))
```

A small helper to run a single-variant GWAS, the same marginal regression as the fine-mapping
practical.

```{r gwas-fn}
run_gwas <- function(X, y) {
  p <- ncol(X)
  beta <- se <- numeric(p)
  for (i in 1:p) {
    fit <- summary(lm(y ~ X[, i] - 1))$coefficients
    beta[i] <- fit[1, "Estimate"]
    se[i]   <- fit[1, "Std. Error"]
  }
  data.frame(beta = beta, se = se, z = beta / se,
             pval = pchisq((beta / se)^2, df = 1, lower.tail = FALSE))
}

gwas1 <- run_gwas(X, y1)
```

Fine-map trait 1 with SuSiE, exactly as in the fine-mapping practical.

```{r susie1}
fm1 <- susie_rss(bhat = gwas1$beta, shat = gwas1$se, n = N,
                 R = in_sample_LD, var_y = var(y1),
                 L = L + 1, estimate_residual_variance = TRUE)

fm1$sets$cs  # credible sets for trait 1
```

Trait 1 has three causal variants (68, 39, and 1), but only variant 39 carries a strong enough
signal to form a credible set. SuSiE localizes the signals it can resolve confidently.

***

# Simulate a second trait: an eQTL

Now we create the QTL phenotype to colocalize against the GWAS. We reuse the **same
genotypes**, so both traits share the same LD, as a GWAS and an eQTL measured in the same
region would. We give the eQTL two causal variants:

1. one **shared** with the GWAS (variant 39, the strong trait-1 signal), and
2. one **private** to the eQTL.

This mirrors the realistic situation that motivates `coloc.susie`. A locus carries several
independent signals, only some of which are shared across traits.

```{r sim-eqtl}
set.seed(2026)

shared_snp  <- 39                 # also causal for the GWAS (trait 1)
private_snp <- 20                 # causal for the eQTL only
causal_ind2 <- c(shared_snp, private_snp)

L2  <- length(causal_ind2)
h2g2 <- 0.15                      # eQTLs often explain more variance than GWAS loci
per_snp_h2g2 <- h2g2 / L2

beta2 <- rep(0, p)
eff2 <- rnorm(L2, mean = 0, sd = sqrt(per_snp_h2g2))
beta2[causal_ind2] <- eff2

genetic_effect2 <- X %*% beta2
var_g2 <- var(genetic_effect2)
sigma_squared2 <- ifelse(1 - var_g2 > 0.05, 1 - var_g2, 0.1)
y2 <- as.numeric(genetic_effect2 + rnorm(N, mean = 0, sd = sqrt(sigma_squared2)))

cat(sprintf("Trait 2 (eQTL) causal variants: %s\n", paste(causal_ind2, collapse = ", ")))
```

Run the eQTL GWAS and fine-map it with SuSiE, the same workflow as trait 1.

```{r susie2}
gwas2 <- run_gwas(X, y2)

fm2 <- susie_rss(bhat = gwas2$beta, shat = gwas2$se, n = N,
                 R = in_sample_LD, var_y = var(y2),
                 L = L2 + 1, estimate_residual_variance = TRUE)

fm2$sets$cs  # credible sets for trait 2
```

***

# Visualize the two signals

Before testing colocalization, always look at the locus. First the marginal association,
plotted as $-\log_{10}(p)$ against variant position. True causal variants are red, and the
lead (most significant) variant of each trait is labeled.

```{r plot-assoc, fig.height=6}
par(mfrow = c(2, 1), mar = c(4, 4, 2, 1))

plot_assoc <- function(gwas, causal, trait) {
  logp <- -log10(gwas$pval)
  lead <- which.max(logp)
  plot(1:p, logp, pch = 16, col = "grey60",
       xlab = "Variant", ylab = expression(-log[10](p)),
       main = sprintf("%s marginal association", trait))
  points(causal, logp[causal], pch = 16, col = "red")
  text(lead, logp[lead], labels = paste0("lead: variant ", lead), pos = 2, cex = 0.9)
}

plot_assoc(gwas1, causal_ind1, "Trait 1 (GWAS)")
plot_assoc(gwas2, causal_ind2, "Trait 2 (eQTL)")
```

The GWAS lead is variant 39, whereas the eQTL lead is variant 20. The shared variant 39 is
present in the eQTL but it is *not* the eQTL's strongest signal. Keep this in mind when we run
`coloc.abf`, which only ever compares the single lead per trait.

Next, the SuSiE fine-mapping output. `susieR` ships its own plotting function, `susie_plot`.
With `y = "PIP"` it draws the posterior inclusion probabilities and colors variants by
credible-set membership. Passing the true effect vector to `b` overlays the simulated causal
variants, which is handy for a teaching example. (`susieR` also provides
`susie_plot(model, y = "z")` for z-scores, plus `susie_plot_iteration` and
`susie_plot_changepoint` for diagnostics.)

```{r susie-plot, fig.height=6}
par(mfrow = c(2, 1), mar = c(4, 4, 2, 1))

b1_true <- numeric(p); b1_true[causal_ind1] <- 1
susie_plot(fm1, y = "PIP", b = b1_true, main = "Trait 1 (GWAS) SuSiE")

b2_true <- numeric(p); b2_true[causal_ind2] <- 1
susie_plot(fm2, y = "PIP", b = b2_true, main = "Trait 2 (eQTL) SuSiE")
```

Each credible set is drawn in its own color, and the true causal variants are circled. Note
that trait 1's two weaker causal variants (68 and 1) are not captured in any credible set.
SuSiE reports only the signals it can localize confidently.

**Q1: Which variant falls in a credible set for both traits, and which is private to the
eQTL?**

<details>
<summary>Answer</summary>
Variant 39 is in a credible set for both the GWAS and the eQTL, so it is the shared signal.
Variant 20 is in a credible set for the eQTL only, so it is private to the eQTL. Importantly,
variant 20 is the eQTL's *strongest* signal, so the two traits' lead variants differ. A
single-causal-variant method like `coloc.abf` will therefore lean toward H3 (distinct
variants), as we confirm below.
</details>

***

# coloc.abf: the single-causal-variant test

`coloc.abf` needs a small list per trait. For a quantitative trait we supply the marginal
effect size `beta`, its variance `varbeta = se^2`, the variant names `snp`, the sample size
`N`, and `sdY`, the standard deviation of the phenotype. We attach simple variant names so
both traits refer to the same SNPs.

```{r build-datasets}
snp_names <- paste0("SNP", 1:p)

D1 <- list(beta = gwas1$beta, varbeta = gwas1$se^2, snp = snp_names,
           type = "quant", N = N, sdY = sd(y1))
D2 <- list(beta = gwas2$beta, varbeta = gwas2$se^2, snp = snp_names,
           type = "quant", N = N, sdY = sd(y2))

# coloc checks that each dataset is well formed
check_dataset(D1)
check_dataset(D2)
```

```{r run-abf}
res.abf <- coloc.abf(dataset1 = D1, dataset2 = D2)
round(res.abf$summary, 4)
```

**Q2: Which hypothesis does `coloc.abf` favor, and what is it concluding about the locus?**

<details>
<summary>Answer</summary>
`PP.H3` is close to 1, so `coloc.abf` concludes that both traits are associated but through
**distinct** causal variants. Under its one-causal-variant assumption it compares only the
single strongest signal in each trait. The GWAS lead is variant 39, whereas the eQTL lead is
the private variant 20, so the single signals look different and the verdict is H3.
</details>

Notice the trap. The locus genuinely *does* share variant 39, but `coloc.abf` cannot see it
because the eQTL's strongest signal sits elsewhere. The single-causal-variant assumption has
thrown away the shared signal. This is exactly the situation `coloc.susie` was built for.

***

# coloc.susie: allowing multiple causal variants

`coloc.susie` takes the SuSiE objects directly. It needs the variant names on the fitted
object so it can match SNPs across traits, so we attach them. (`runsusie()` does this for
you automatically; here we add them to the objects produced above.)

```{r name-susie}
name_susie <- function(S, snps) {
  colnames(S$alpha)        <- snps
  colnames(S$lbf_variable) <- snps
  names(S$pip)             <- snps
  S
}

fm1 <- name_susie(fm1, snp_names)
fm2 <- name_susie(fm2, snp_names)
```

```{r run-susie-coloc}
res.susie <- coloc.susie(fm1, fm2)
# coloc.susie returns a data.table, which does not auto-print inside knitr; coerce to print it
as.data.frame(res.susie$summary)
```

`coloc.susie` returns **one row per pair of credible sets**, with `hit1` and `hit2` naming
the lead variant of each set. Read the table by looking for rows with high `PP.H4`.

**Q3: Which pair of signals colocalizes, and what is its PP.H4?**

<details>
<summary>Answer</summary>
The pair `hit1 = SNP39`, `hit2 = SNP39` has `PP.H4` close to 1. The shared variant 39 is
recovered as a colocalizing signal. Every pairing that involves the eQTL's private variant
20 returns high `PP.H3` instead, correctly flagging those as distinct.
</details>

**Q4: `coloc.abf` said H3 and `coloc.susie` found a shared signal. Are they contradicting
each other?**

<details>
<summary>Answer</summary>
No. They answer different questions. `coloc.abf` asks whether the **single** strongest
signal is shared, and at this locus it is not (39 versus 20). `coloc.susie` decomposes the
locus into independent signals and asks the question per signal, so it can find that one of
several signals, variant 39, is shared while another is private. When a locus may have more
than one causal variant, prefer `coloc.susie`.
</details>

***

# Interpretation and practical guidance

- **Always look at the locus first.** Marginal Manhattan plots and SuSiE PIPs tell you how
  many signals are present before you trust any single number.
- **Let LD availability pick the method, not your guess about signal count.** Prefer
  `coloc.susie` whenever you have a reliable in-sample LD matrix,
  because you usually cannot know in advance that a locus carries a single causal variant.
  Fall back to `coloc.abf` when no trustworthy LD is available, accepting its
  single-causal-variant assumption as the price of needing only summary statistics.
- **LD must match the summary statistics.** SuSiE and `coloc.susie` are sensitive to LD
  misspecification. Use an in-sample LD matrix, or a reference panel that matches the GWAS
  ancestry, exactly as in the fine-mapping practical.
- **Priors matter.** `coloc.abf` uses prior probabilities `p1`, `p2`, and `p12`. The default
  `p12 = 1e-5` is reasonable but worth a sensitivity check, for example with
  `sensitivity(res.abf, "H4 > 0.5")`.
- **PP.H4 is necessary, not sufficient.** A high PP.H4 says the data are consistent with one
  shared variant. It does not by itself prove causality, and it can be inflated by shared LD
  structure or poorly matched reference panels.

**Q5: You fine-map a GWAS locus and an eQTL and each resolves a single credible set, both
centered on the same variant. Which method would you reach for, and what result do you
expect?**

<details>
<summary>Answer</summary>
You already have the in-sample LD that fine-mapping required, so reach for `coloc.susie`. It
reduces to the same shared-variant conclusion here, with high `PP.H4`, but it stays correct
if a weak second signal is present that the single credible set did not resolve. Reserve
`coloc.abf` for the case where you have only marginal summary statistics and no reliable LD
matrix. Its single-causal-variant assumption is a real limitation, not a safe default,
because correlated or anti-correlated causal effects at the same locus are common
(Kanai et al. 2021; Zhang et al. 2023).
</details>

***

# Resources

- Giambartolomei et al. 2014, the original `coloc.abf` method.
  [PLoS Genet](https://doi.org/10.1371/journal.pgen.1004383)
- Wallace 2021, `coloc.susie` for multiple causal variants.
  [PLoS Genet](https://doi.org/10.1371/journal.pgen.1009440)
- Kanai et al. 2021, fine-mapping across diverse populations, with examples of negatively
  correlated causal variants at a locus (e.g. *FGFR4*, *CETP*).
  [medRxiv](https://doi.org/10.1101/2021.09.03.21262975)
- Zhang et al. 2023, pervasive correlations between causal effects of proximal SNPs, from
  linkage masking under stabilizing selection.
  [medRxiv](https://doi.org/10.1101/2023.12.04.23299391)
- `coloc` package vignettes:
  [chr1glh24.github.io/coloc](https://chr1swallace.github.io/coloc/)

***

# Appendix: using `runsusie` instead

When you fine-map specifically with colocalization downstream, `coloc::runsusie()` is a
cleaner entry point than calling `susie_rss` directly. It wraps `susie_rss`, carries the LD
and SNP names through automatically, and returns an object ready for `coloc.susie`. It expects
a coloc-style dataset that also includes the LD matrix.

```{r runsusie-demo, eval=FALSE}
# runsusie checks that the LD matrix is labelled with the same SNP names as the dataset
rownames(in_sample_LD) <- snp_names
colnames(in_sample_LD) <- snp_names

D1_ld <- c(D1, list(LD = in_sample_LD, position = 1:p))
D2_ld <- c(D2, list(LD = in_sample_LD, position = 1:p))

s1 <- runsusie(D1_ld)
s2 <- runsusie(D2_ld)

coloc.susie(s1, s2)$summary
```
