1 Objectives

In this practical, we will run SBayesRC, which allows us to fit multi-million SNPs simultaneously in practice and incorporate functional annotations. SBayesRC extends SBayesR by letting SNP effects share information through functional genomic annotations (Zheng et al., Nature Genetics 2024; GCTB tutorial).

Note: R code is shown with a light blue background, while terminal commands are shown with a light orange background.

2 Annotation

In real analyses, annotations often mark whether a SNP falls in a functional class (e.g. coding sequence, enhancer, promoter, or evolutionarily conserved region). Each annotation is usually a 0/1 indicator (in or out of that class). You can also use quantitative annotation values. GCTB expects the first two columns to be SNP ID and an intercept column of ones for all SNPs, followed by the annotation indicators. For example:

SNP Intercept Coding_UCSC Conserved_LindbladToh CTCF_Hoffman DGF_ENCODE
rs539032812 1 0 0 0 0
rs12238997 1 0 0 0 1
rs371890604 1 0 0 0 0
rs200531508 1 0 0 0 0
rs144155419 1 0 0 0 0
rs116587930 1 0 0 0 0
rs148120343 1 0 0 0 0
rs61770163 1 0 0 1 0

For teaching we supply annotations.txt: 11 indicator columns (each SNP is in or out of that category). For the first 10 columns, each lists 50 of the 100 simulated causal SNPs together with 50 additional SNPs drawn at random from the non-causal genome-wide list. For the last column, we randomly sampled 1 or 0 with equal probability for all SNPs. That design is artificial but keeps the file compact and gives SBayesRC a controlled analogue of “enriched” functional categories.

3 Run SBayesRC with genomic annotations

Use the same imputed GWAS file as SBayesR:gwas.imputed.ma. The annotation file must include all the SNPs in the LD reference (could have more SNPs). annotations.txt has 13 columns: SNP ID, a column of ones for all SNPs (intercept), then 11 annotation indicators (0/1).

Set --sbayes RC and pass the annotation file. The run usually takes a little longer than SBayesR.

gctb --ldm-eigen ldm_eigen \
--gwas-summary gwas.imputed.ma \
--annot annotations.txt \
--sbayes RC \
--chain-length 500 \
--burn-in 100 \
--thread 4 \
--out sbayesrc

Q1: while SBayesRC is running, think about what would be the potential benefit of incorporating these annotations into the prediction model?

Show answer Annotations let the model up-weight SNP effect sizes in categories that enrich for causal variants (here, Anno1–10), borrowing information across SNPs in the same functional class and potentially improving prediction when annotations truly capture biology.


Besides sbayesrc.snpRes and sbayesrc.parRes, GCTB writes sbayesrc.parSetRes (annotation-level parameters) and sbayesrc.enrich (heritability enrichment per annotation). Example enrichment output (ignore marginal vs joint enrichment which is not relevant to the theme of this practical):

Parameter Annotation Mean SD PostProbAboveOne
Marginal_Heritability_Enrichment Intercept 1.000000 0.000000 0.000000
Marginal_Heritability_Enrichment Anno1 1160.841064 115.182472 1.000000
Marginal_Heritability_Enrichment Anno2 815.068909 88.783798 1.000000
Marginal_Heritability_Enrichment Anno3 960.037109 82.534836 1.000000
Marginal_Heritability_Enrichment Anno4 1563.750366 98.788918 1.000000
Marginal_Heritability_Enrichment Anno5 878.377441 96.052391 1.000000
Marginal_Heritability_Enrichment Anno6 1348.684204 128.180038 1.000000
Marginal_Heritability_Enrichment Anno7 1190.434814 82.715630 1.000000
Marginal_Heritability_Enrichment Anno8 1037.594116 81.555962 1.000000
Marginal_Heritability_Enrichment Anno9 835.093140 92.254410 1.000000
Marginal_Heritability_Enrichment Anno10 1364.147583 85.787529 1.000000
Marginal_Heritability_Enrichment Anno11 0.930061 0.088228 0.220000

Q2: The result shows strong heritability enrichment in Anno1-10 but no enrichment in Anno11. Is this as expected?

Show answer Strong enrichment in Anno1–10 and none in Anno11 is expected (each of Anno1–10 was built to contain half the causal SNPs; Anno11 is random noise). SBayesRC produces more high-PIP SNPs than SBayesR because annotation sharing pulls signal toward annotated regions and helps identify the causal variants.


Q3: How many SNPs with high PIP (check sbayesrc.snpRes)? How does it compare to SBayesR? What does it indicate?

Show answer SBayesRC produces more high-PIP SNPs than SBayesR because annotation sharing pulls signal toward annotated regions and helps identify the causal variants.


5 Run SBayesRC with human height data

Here, we will use the data from UKB height and 96 functional annotations from S-LDSC model to demonstrate how to run SBayesRC using GCTB. Note that you can also run SBayesRC using a R package that we have developed: https://github.com/zhilizheng/SBayesRC.

There are 7.4M SNPs in the summary statistics for UKB height. The first step is to match and impute the missing SNPs in the LD reference. We will use the LD reference data generated above.

The default MCMC chain length is 3000 with the first 1000 iterations as burn-in. Here, for the interest of time, we set a shorter chain length because it takes more time to run model incorporating annotations. In your own real trait analysis, it is recommended to use the default setting or even longer chain.

gwas="HT.ma"
annot="baselineLD.txt"

## Match and impute SNPs
gctb --ldm-eigen ldm_eigen --gwas-summary $gwas --impute-summary --thread 4 --out HT

## Run SBayesRC
gctb --ldm-eigen ldm_eigen \
     --gwas-summary HT.imputed.ma \
     --annot $annot \
     --sbayes RC \
     --thread 4 \
     --chain-length 500 \
     --burn-in 100 \
     --out HT_sbayesrc

You can check the SNP-based heritability enrichment across functional annotations in R code:

enrich = read.table("HT_sbayesrc.enrich", header=T)
enrich = subset(enrich, Parameter == "Marginal_Heritability_Enrichment")

# plot the per-SNP heritability enrichment for all annotations
library(ggplot2)
enrich$Annotation = factor(enrich$Annotation, levels=enrich$Annotation[order(enrich$Mean, decreasing = T)])
p1 = ggplot(enrich, aes(x=Annotation, y=Mean, fill=Annotation)) + geom_bar(stat="identity", position="dodge") + 
  geom_errorbar(aes(ymin=Mean-SD, ymax=Mean+SD, color=Annotation)) + geom_hline(yintercept=1, linetype="dashed") +
  xlab("Functional annotation") + ylab("Per-SNP heritability enrichment") +
  guides(fill="none", color="none") +
  theme(axis.text = element_text(size=5, angle = 90, vjust = 0.5, hjust=1))
ggsave("height_h2_enrich_all.pdf", p1, width=12, height=6)
p1

Highlight some important functional annotations:

# highlight some annotations
selected = c("Ancient_Sequence_Age_Human_Promoter", "Human_Promoter_Villar_ExAC", "Intron_UCSC", "Coding_UCSC", "Conserved_LindbladToh", "TSS_Hoffman", "UTR_5_UCSC", "Enhancer_Andersson", "UTR_3_UCSC", "Transcr_Hoffman", "Repressed_Hoffman")
p2 = ggplot(subset(enrich, Annotation %in% selected), aes(x=Annotation, y=Mean, fill=Annotation)) + geom_bar(stat="identity", position="dodge") + 
  geom_errorbar(aes(ymin=Mean-SD, ymax=Mean+SD, color=Annotation)) + geom_hline(yintercept=1, linetype="dashed") +
  xlab("Functional annotation") + ylab("Per-SNP heritability enrichment") +
  guides(fill="none", color="none") +
  theme(axis.text = element_text(angle = 90, vjust = 0.5, hjust=1))
ggsave("height_h2_enrich_highlight.pdf", p2, width=12, height=6)
p2

With SBayesRC result, we can gain some insights into the functional architecture of the trait:

library(stringr)
library(dplyr)
library(tidyr)

# Read results
res = read.table("HT_sbayesrc.parSetRes", header=T)
names(res) = c("Par","Anno","Est","SE","PP")  # the columns are parameter, annotation, estimate, standard error, and posterior probability of greater than zero

# Extract Pi parameter estimates
pi = res[str_detect(res$Par, "AnnoJointProb"), ]
pi = pi %>% separate(Par, into = c("Pi", "Component"), sep = "_pi")

# Identify the `Est` value for the `Intercept` row within each component
intercept_values <- pi %>% 
  filter(Anno == "Intercept") %>% 
  select(Component, Est) %>% 
  rename(Intercept_Est = Est)

# Merge these intercept values back to the original data frame
pi_with_intercept <- pi %>% left_join(intercept_values, by = "Component")

# Subtract the intercept value from the `Est` values
pi_with_intercept <- pi_with_intercept %>% mutate(Adjusted_Est = Est - Intercept_Est)

# Create a new data frame with the adjusted values
adjusted_pi <- pi_with_intercept %>% select(Component, Anno, Adjusted_Est, SE, PP)

# plot the result focusing on the selected annotation and omit component 5 because there is no signal in there
facet_labels <- c("1" = "Zero", "2" = "Small effect", "3" = "Moderate effect", "4" = "Large effect")
p3 = ggplot(subset(adjusted_pi, Anno %in% selected & Component != 5), aes(x=Anno, y=Adjusted_Est, fill=Anno)) + geom_bar(stat="identity", position="dodge") + 
  geom_errorbar(aes(ymin=Adjusted_Est-SE, ymax=Adjusted_Est+SE, color=Anno)) + geom_hline(yintercept=0, linetype="dashed") +
  facet_grid(.~Component, scales="free", labeller = labeller(Component = facet_labels)) +
  xlab("Functional annotation") + ylab("Deviation of Pi to average") +
  coord_flip() +
  guides(fill="none", color="none") +
  theme(axis.text = element_text(angle = 0, vjust = 0.5, hjust=1))
ggsave("height_pi_enrich_highlight.pdf", p3, width=12, height=6)
p3

Despite large standard errors, evolutionary conserved regions, coding sequences, promoters and transcription start sites (TSS) tend to be enriched in causal variants with large effect sizes, whereas repressed regions and introns tend to be depleted in causal variants. Such information learned from the data will, in turn, help with better estimate the SNP effects and therefore improve prediction. Since we cannot provide validation with real genotypes and phenotypes, we are not doing prediction for this UKB height data in this practical exercise.