前言
这篇文章的书写，是受启发于《post-GWAS：使用coloc进行共定位分析(Colocalization)》，目的是想了解一下软件coloc是如何巧妙的运用贝叶斯方法的
共定位的基本原理
关于coloc的基本原理可以参考这四篇文章：
- 《Eliciting priors and relaxing the single causal variant assumption in colocalisation analyses》
- 《Bayesian Test for Colocalisation between Pairs of Genetic Association Studies Using Summary Statistics》
- 《A Bayesian framework for multiple trait colocalization from summary association statistics》
- 《Bayes Factors for Genome-Wide Association Studies: Comparison with P-values》
我们还是以GWAS数据和eQTL数据为例：
注：位于上方的区间理解为GWAS数据，位于下方的区间理解为eQTL数据；橙色表征该snp与 trait 1 显著相关，蓝色表征该snp与 trait 2 显著相关
对于某一段基因组区间来说，将会有下面四种假设：
- H0 ：No association with either trait
- H1 : Association with trait 1, not with trait 2
- H2 : Association with trait 2, not with trait 1
- H3 : Association with trait 1 and trait 2, two independent SNPs
- H4 : 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，pik表示在H1条件下发生的概率，1-pik表示在H0条件下发生的概率
而 zi代表gwas，eQTL数据中的p_val，经过均值为0的正态分布转换后的似然值；事实上，p_val越小，代表该位点与性状的关联越显著（越能接受H1），因此p_val对应转换后的似然值越大，即zi越大，进而推出 pik 越大，1-pik 越小，则H1发生的概率越大，即贝叶斯因子越大，越没有理由拒绝H1
因此，最终有zk代表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的数据：
eQTL
- id代表不同的snp
- variant_id 包含该snp的位置信息和突变信息
- gene_id 代表这些snp调控的基因，一般做此类分析我们focus在一个基因上就好，这一个基因通常控制着一个性状
- tss_distance 代表某snp距离这个基因的距离
- rna_count 代表该基因表达量
- maf 代表次等位基因频率
- slope 可以理解为某snp对该基因调控的效应
- p值为显著性
对于gwas数据：
gwas
- id代表不同的snp
- variant_id 包含该snp的位置信息和突变信息
- beta 可以理解为某snp对该基因调控的效应
- 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来说
- snp 代表此snp的名字
- pvalues.df1（df2）代表gwas（eQTL）数据的p_val
- MAF.df1（df2）代表eQTL数据的MAF值，gwas数据的那一列MAF也由eQTL数据的MAF值代替
- v.df1（df2）代表gwas（eQTL）数据的效应值
- z.df1（df2）代表gwas（eQTL）数据的p_val，经过均值为0的正态分布转换后的似然值
- r.df1（df2）代表gwas（eQTL）数据的矫正因子 r
- lABF.df1（df2）代表gwas（eQTL）数据的H0，H1，H2，H3，H4
log(贝叶斯因子)之和- internal.sum.lABF 代表 lABF.df1 + lABF.df2
- 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
