前言
這篇文章的書寫锡溯，是受啟發(fā)于《post-GWAS：使用coloc進行共定位分析(Colocalization)》，目的是想了解一下軟件coloc是如何巧妙的運用貝葉斯方法的
共定位的基本原理
關(guān)于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數(shù)據(jù)和eQTL數(shù)據(jù)為例：
注：位于上方的區(qū)間理解為GWAS數(shù)據(jù)氧猬，位于下方的區(qū)間理解為eQTL數(shù)據(jù)垛叨；橙色表征該snp與 trait 1 顯著相關(guān)，藍色表征該snp與 trait 2 顯著相關(guān)
對于某一段基因組區(qū)間來說儡遮，將會有下面四種假設(shè)：
- 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的數(shù)據(jù)，還是GWAS的數(shù)據(jù)荐类，都是基于混合線性模型將基因型和性狀給聯(lián)系起來的怖现。但是無論是eQTL的數(shù)據(jù)，還是GWAS的數(shù)據(jù)玉罐，一般就只能關(guān)聯(lián)一種性狀屈嗤；那么共定位的目的就是要將某些關(guān)聯(lián)到多個性狀的snp給找出來
而eQTL的優(yōu)勢是能夠找到snp能夠調(diào)控某基因的表達潘拨，從而影響某一種性狀（eQTL往往是通過snp影響某個基因表達，而這個基因控制著某一種性狀恢共，從而將snp和性狀聯(lián)系起來）；而gwas的又是在全基因組的角度璧亚，尋找與某性狀相關(guān)聯(lián)的snp
如下圖：
對于某一段的基因組區(qū)間來說讨韭，0代表與trait無關(guān)，1代表與trait有關(guān)癣蟋，那么第三張圖就表示該snp與兩種性狀都相關(guān)聯(lián)透硝，因此共定位就是想找出這一些與兩個性狀都相關(guān)的snp
貝葉斯因子
假設(shè)：
那么等式
貝葉斯因子的大小
貝葉斯因子介于0—1之間，越接近1疯搅，代表越支持H1濒生；越接近0，則越支持H0
由式子我們可以看出貝葉斯因子是連接先驗概率和后驗概率的橋梁幔欧，然而在該項目中罪治，計算貝葉斯因子是一件比較困難的事情，因此作者采用了ABF(Approximate Bayes Factor)來近似代替貝葉斯因子礁蔗。
ABF利用gwas觉义，eQTL數(shù)據(jù)中的p_val與貝葉斯因子之間建立邏輯回歸模型，做擬合浴井，然后得到表達式
具體推導：
對于H0與H1兩種情況的事件來說晒骇，貝葉斯因子可以寫作是：其中，i 代表每一個snp磺浙，pik表示在H1條件下發(fā)生的概率洪囤，1-pik表示在H0條件下發(fā)生的概率
而 zi代表gwas，eQTL數(shù)據(jù)中的p_val撕氧，經(jīng)過均值為0的正態(tài)分布轉(zhuǎn)換后的似然值瘤缩；事實上，p_val越小伦泥，代表該位點與性狀的關(guān)聯(lián)越顯著（越能接受H1）款咖，因此p_val對應(yīng)轉(zhuǎn)換后的似然值越大，即zi越大奄喂，進而推出 pik 越大铐殃，1-pik 越小，則H1發(fā)生的概率越大跨新，即貝葉斯因子越大富腊，越?jīng)]有理由拒絕H1
因此，最終有zk代表gwas域帐，eQTL數(shù)據(jù)中的p_val赘被，經(jīng)過正態(tài)分布轉(zhuǎn)換后的似然值的總和
總結(jié)：這樣做擬合是整，我們就只需要輸入gwas，eQTL數(shù)據(jù)中的p_val民假，就可以得到對應(yīng)的ABF了浮入，這樣的擬合連接了p_val和ABF，使計算更為方便了
對于四種假設(shè)羊异，我們計算后驗概率有：
而p0事秀，p1，p2野舶，p12分別代表H0易迹，H1，H2平道，H4假設(shè)下的先驗概率
代碼分析
注明出處：《post-GWAS：使用coloc進行共定位分析(Colocalization)》睹欲，我用了這篇文章的示例數(shù)據(jù)
# 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數(shù)據(jù)和eQTL數(shù)據(jù)共有的snp(rs_id)
input <- merge(eqtl, gwas, by="rs_id", all=FALSE, suffixes=c("_eqtl","_gwas"))
# type="cc" 表示二元性狀
# type="quant" 表示連續(xù)性狀
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)
從源代碼分析，對于函數(shù)coloc.abf()
# 讀取數(shù)據(jù)
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的數(shù)據(jù)：
eQTL
- id代表不同的snp
- variant_id 包含該snp的位置信息和突變信息
- gene_id 代表這些snp調(diào)控的基因一屋，一般做此類分析我們focus在一個基因上就好窘疮，這一個基因通常控制著一個性狀
- tss_distance 代表某snp距離這個基因的距離
- rna_count 代表該基因表達量
- maf 代表次等位基因頻率
- slope 可以理解為某snp對該基因調(diào)控的效應(yīng)
- p值為顯著性
對于gwas數(shù)據(jù)：
gwas
- id代表不同的snp
- variant_id 包含該snp的位置信息和突變信息
- beta 可以理解為某snp對該基因調(diào)控的效應(yīng)
- p值為顯著性
我們看看源代碼：
# 定義H1冀墨，H2考余，H4的先驗概率
## 默認概率如下
p1=1e-4
p2=1e-4
p12=1e-5
# 預(yù)處理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)
# 新增一列g(shù)was和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).對于函數(shù)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.函數(shù)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轉(zhuǎn)換為正態(tài)分布的似然值
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).對于函數(shù)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調(diào)控兩種性狀楚堤，并且將gwas和eQTL數(shù)據(jù)聯(lián)系到了一起
結(jié)果解讀
對于其中一個snp來說
- snp 代表此snp的名字
- pvalues.df1（df2）代表gwas（eQTL）數(shù)據(jù)的p_val
- MAF.df1（df2）代表eQTL數(shù)據(jù)的MAF值，gwas數(shù)據(jù)的那一列MAF也由eQTL數(shù)據(jù)的MAF值代替
- v.df1（df2）代表gwas（eQTL）數(shù)據(jù)的效應(yīng)值
- z.df1（df2）代表gwas（eQTL）數(shù)據(jù)的p_val含懊，經(jīng)過均值為0的正態(tài)分布轉(zhuǎn)換后的似然值
- r.df1（df2）代表gwas（eQTL）數(shù)據(jù)的矯正因子 r
- lABF.df1（df2）代表gwas（eQTL）數(shù)據(jù)的H0身冬，H1，H2岔乔，H3酥筝，H4
log(貝葉斯因子)之和- internal.sum.lABF 代表 lABF.df1 + lABF.df2
- SNP.PP.H4 代表H4的貝葉斯因子
# 新增一列g(shù)was和eQTL的 logABF 的加和值
## lABF.df1 + lABF.df2 相當于某snp gwas和eQTL數(shù)據(jù)貝葉斯因子的乘積
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
