使用贝叶斯方法进行共定位的Coloc软件

前言

这篇文章的书写，是受启发于《post-GWAS：使用coloc进行共定位分析(Colocalization)》，目的是想了解一下软件coloc是如何巧妙的运用贝叶斯方法的

共定位的基本原理

关于coloc的基本原理可以参考这四篇文章：

1. 《Eliciting priors and relaxing the single causal variant assumption in colocalisation analyses》
2. 《Bayesian Test for Colocalisation between Pairs of Genetic Association Studies Using Summary Statistics》
3. 《A Bayesian framework for multiple trait colocalization from summary association statistics》
4. 《Bayes Factors for Genome-Wide Association Studies: Comparison with P-values》

我们还是以GWAS数据和eQTL数据为例：

注：位于上方的区间理解为GWAS数据，位于下方的区间理解为eQTL数据；橙色表征该snp与 trait 1 显著相关，蓝色表征该snp与 trait 2 显著相关

对于某一段基因组区间来说，将会有下面四种假设：

1. H₀ ：No association with either trait
2. H₁ : Association with trait 1, not with trait 2
3. H₂ : Association with trait 2, not with trait 1
4. H₃ : Association with trait 1 and trait 2, two independent SNPs
5. H₄ : Association with trait 1 and trait 2, one shared SNP

那么怎么理解共定位呢？
首先，无论是eQTL的数据，还是GWAS的数据，都是基于混合线性模型将基因型和性状给联系起来的。但是无论是eQTL的数据，还是GWAS的数据，一般就只能关联一种性状；那么共定位的目的就是要将某些关联到多个性状的snp给找出来

而eQTL的优势是能够找到snp能够调控某基因的表达，从而影响某一种性状（eQTL往往是通过snp影响某个基因表达，而这个基因控制着某一种性状，从而将snp和性状联系起来）；而gwas的又是在全基因组的角度，寻找与某性状相关联的snp

如下图：

对于某一段的基因组区间来说，0代表与trait无关，1代表与trait有关，那么第三张图就表示该snp与两种性状都相关联，因此共定位就是想找出这一些与两个性状都相关的snp

贝叶斯因子

假设：

那么等式定义为贝叶斯因子

贝叶斯因子介于0—1之间，越接近1，代表越支持H1；越接近0，则越支持H0

由式子我们可以看出贝叶斯因子是连接先验概率和后验概率的桥梁，然而在该项目中，计算贝叶斯因子是一件比较困难的事情，因此作者采用了ABF(Approximate Bayes Factor)来近似代替贝叶斯因子。

ABF利用gwas，eQTL数据中的p_val与贝叶斯因子之间建立逻辑回归模型，做拟合，然后得到表达式

具体推导：
对于H0与H1两种情况的事件来说，贝叶斯因子可以写作是：

其中，i 代表每一个snp，p_ik表示在H1条件下发生的概率，1-p_ik表示在H0条件下发生的概率
而 z_i代表gwas，eQTL数据中的p_val，经过均值为0的正态分布转换后的似然值；事实上，p_val越小，代表该位点与性状的关联越显著（越能接受H1），因此p_val对应转换后的似然值越大，即z_i越大，进而推出 p_ik 越大，1-p_ik 越小，则H1发生的概率越大，即贝叶斯因子越大，越没有理由拒绝H1
因此，最终有

z_k代表gwas，eQTL数据中的p_val，经过正态分布转换后的似然值的总和
总结：这样做拟合，我们就只需要输入gwas，eQTL数据中的p_val，就可以得到对应的ABF了，这样的拟合连接了p_val和ABF，使计算更为方便了

对于四种假设，我们计算后验概率有：

而p0，p1，p2，p12分别代表H0，H1，H2，H4假设下的先验概率

代码分析

注明出处：《post-GWAS：使用coloc进行共定位分析(Colocalization)》，我用了这篇文章的示例数据

# install
library(remotes)
install_github("chr1swallace/coloc",build_vignettes=TRUE)
library("coloc")
library(dplyr)

# 读取文件
gwas <- read.table(file="GWAS.txt", header=T)
eqtl <- read.table(file="eQTL.txt", header=T)

# 运行
## 选取gwas数据和eQTL数据共有的snp(rs_id)
input <- merge(eqtl, gwas, by="rs_id", all=FALSE, suffixes=c("_eqtl","_gwas"))

# type="cc" 表示二元性状
# type="quant" 表示连续性状
result <- coloc.abf(dataset1=list(pvalues=input$pval_nominal_gwas, type="cc", s=0.33, N=50000), 
                    dataset2=list(pvalues=input$pval_nominal_eqtl, type="quant", N=10000), MAF=input$maf)

# 筛选后验概率大于0.95的位点
need_result=result$results %>% filter(SNP.PP.H4 > 0.95)

从源代码分析，对于函数coloc.abf()

# 读取数据
dataset1=list(pvalues=input$pval_nominal_gwas, type="cc", s=0.33, N=50000)
dataset2=list(pvalues=input$pval_nominal_eqtl, type="quant", N=10000)
MAF=input$maf

对于eQTL的数据：

1. id代表不同的snp
2. variant_id 包含该snp的位置信息和突变信息
3. gene_id 代表这些snp调控的基因，一般做此类分析我们focus在一个基因上就好，这一个基因通常控制着一个性状
4. tss_distance 代表某snp距离这个基因的距离
5. rna_count 代表该基因表达量
6. maf 代表次等位基因频率
7. slope 可以理解为某snp对该基因调控的效应
8. p值为显著性

对于gwas数据：

1. id代表不同的snp
2. variant_id 包含该snp的位置信息和突变信息
3. beta 可以理解为某snp对该基因调控的效应
4. p值为显著性

我们看看源代码：

# 定义H1，H2，H4的先验概率
## 默认概率如下
p1=1e-4
p2=1e-4
p12=1e-5

# 预处理data，计算gwas和eQTL的贝叶斯因子
## (1) process.dataset
df1 <- process.dataset(d=dataset1, suffix="df1")
df2 <- process.dataset(d=dataset2, suffix="df2")

merged.df <- merge(df1,df2)
# 新增一列gwas和eQTL的 logABF 的加和值
merged.df$internal.sum.lABF <- with(merged.df, lABF.df1 + lABF.df2)
## add SNP.PP.H4 - post prob that each SNP is THE causal variant for a shared signal
my.denom.log.abf <- logsum(merged.df$internal.sum.lABF)

## 计算 H4 的贝叶斯因子
merged.df$SNP.PP.H4 <- exp(merged.df$internal.sum.lABF - my.denom.log.abf)

pp.abf <- combine.abf(merged.df$lABF.df1, merged.df$lABF.df2, p1, p2, p12)  
common.snps <- nrow(merged.df)
results <- c(nsnps=common.snps, pp.abf)

output<-list(summary=results,
             results=merged.df,
             priors=c(p1=p1,p2=p2,p12=p12))

(1).对于函数process.dataset()

d$snp <- sprintf("SNP.%s",1:length(d$pvalues))
df <- data.frame(pvalues = d$pvalues,
                   MAF = d$MAF,
                   N=d$N,
                   snp=as.character(d$snp))    
snp.index <- which(colnames(df)=="snp")
colnames(df)[-snp.index] <- paste(colnames(df)[-snp.index], suffix, sep=".")
keep <- which(df$MAF>0 & df$pvalues > 0) # all p values and MAF > 0
df <- df[keep,]

## 计算 ABF(Approximate Bayes Factor)
abf <- approx.bf.p(p=df$pvalues, f=df$MAF, type=d$type, N=df$N, s=d$s, suffix=suffix)
df <- cbind(df, abf)

1.函数approx.bf.p()

p=df$pvalues
f=df$MAF
type=d$type
N=df$N
s=d$s
suffix=suffix
  
if(type=="quant") {
    sd.prior <- 0.15
    V <- Var.data(f, N)
} else {
    sd.prior <- 0.2
    V <- Var.data.cc(f, N, s)
}
# 将gwas和eQTL的p_val转换为正态分布的似然值
z <- qnorm(0.5 * p, lower.tail = FALSE)
## Shrinkage factor: ratio of the prior variance to the total variance
## 计算校正因子 r
r <- sd.prior^2 / (sd.prior^2 + V)
## Approximate BF  # I want ln scale to compare in log natural scale with LR diff
## 计算 ABF 的log值
lABF = 0.5 * (log(1-r) + (r * z^2))
ret <- data.frame(V,z,r,lABF)
if(!is.null(suffix))
   colnames(ret) <- paste(colnames(ret), suffix, sep=".")

(2).对于函数combine.abf()

combine.abf <- function(l1, l2, p1, p2, p12) {
  lsum <- l1 + l2
  lH0.abf <- 0
  ## 计算H1的后验概率
  lH1.abf <- log(p1) + logsum(l1)
  ## 计算H2的后验概率
  lH2.abf <- log(p2) + logsum(l2)
  ## 计算H3的后验概率
  lH3.abf <- log(p1) + log(p2) + logdiff(logsum(l1) + logsum(l2), logsum(lsum))
  ## 计算H4的后验概率
  lH4.abf <- log(p12) + logsum(lsum)
  
  all.abf <- c(lH0.abf, lH1.abf, lH2.abf, lH3.abf, lH4.abf)
  my.denom.log.abf <- logsum(all.abf)
  pp.abf <- exp(all.abf - my.denom.log.abf)
  names(pp.abf) <- paste("PP.H", (1:length(pp.abf)) - 1, ".abf", sep = "")
  return(pp.abf)
}

代码中注释的计算后验概率原理如下图：

那么我们找到第五项的后验概率高的话，代表该snp调控两种性状，并且将gwas和eQTL数据联系到了一起

结果解读

对于其中一个snp来说

1. snp 代表此snp的名字
2. pvalues.df1（df2）代表gwas（eQTL）数据的p_val
3. MAF.df1（df2）代表eQTL数据的MAF值，gwas数据的那一列MAF也由eQTL数据的MAF值代替
4. v.df1（df2）代表gwas（eQTL）数据的效应值
5. z.df1（df2）代表gwas（eQTL）数据的p_val，经过均值为0的正态分布转换后的似然值
6. r.df1（df2）代表gwas（eQTL）数据的矫正因子 r
7. lABF.df1（df2）代表gwas（eQTL）数据的H0，H1，H2，H3，H4
   log(贝叶斯因子)之和
8. internal.sum.lABF 代表 lABF.df1 + lABF.df2
9. SNP.PP.H4 代表H4的贝叶斯因子

# 新增一列gwas和eQTL的 logABF 的加和值
## lABF.df1 + lABF.df2 相当于某snp gwas和eQTL数据贝叶斯因子的乘积
merged.df$internal.sum.lABF <- with(merged.df, lABF.df1 + lABF.df2)
## add SNP.PP.H4 - post prob that each SNP is THE causal variant for a shared signal
my.denom.log.abf <- logsum(merged.df$internal.sum.lABF)

## 计算 H4 的贝叶斯因子
merged.df$SNP.PP.H4 <- exp(merged.df$internal.sum.lABF - my.denom.log.abf)

如上图所示，SNP.PP.H4 = 1代表有强烈的理由接受H4