GWAS找到显著信号位点后,需要解释显著信号位点如何影响表型。
常见的一个解释方法是共定位分析。
主流的共定位分析包括:
- 1)GWAS和eQTL共定位;
- 2)GWAS和sQTL共定位;
- 3)GWAS和meQTL共定位;
- 4)GWAS和pQTL共定位;
其中,GWAS和eQTL共定位应用最为广泛。
具体来说,当检测到GWAS信号和eQTL共定位时,我们会认为GWAS信号上的位点可能通过改变基因表达的生物学过程影响表型。
共定位分析有四种设想:
第一种设想 H0: 表型1(GWAS)和 表型2 (以eQTL为例)与某个基因组区域的所有SNP位点无显著相关;
第二种设想 H1/H2: 表型1(GWAS)或表型2(以eQTL为例)与某个基因组区域的SNP位点显著相关;
第三种设想 H3: 表型1(GWAS)和 表型2 (以eQTL为例)与某个基因组区域的SNP位点显著相关,但由不同的因果变异位点驱动;
第四种设想 H4: 表型1(GWAS)和 表型2 (以eQTL为例)与某个基因组区域的SNP位点显著相关,且由同一个因果变异位点驱动;
基于以上四种设想,我们希望第四种设想 H4 在统计学上概率更高,这样就能解释显著信号位点如何影响表型;
所以共定位分析,本质上是在检验第四种的后验概率;
讲完共定位分析的相关概念,接下来我们以“GWAS和eQTL共定位”为例讲一下如何使用coloc进行共定位分析。
1 下载、安装coloc
if(!require("remotes"))
install.packages("remotes")
install.packages("dplyr")
library(remotes)
install_github("chr1swallace/coloc",build_vignettes=TRUE)
library("coloc")
library(dplyr)
2 下载测试数据
测试数据请在公众号“bio生物信息”后台回复"coloc"获取。
3 分析
3.1 导入表型1(GWAS)数据:
gwas <- read.table(file="E:/path_to_GWAS/GWAS.txt", header=T);
GWAS数据包括:rs编号rs_id
,P值pval_nominal
等:
注意事项:如果表型是二分类变量(case和control),输入文件二选一:
1)rs编号
rs_id
、P值pval_nominal
、SNP的效应值beta
、效应值方差varbeta
;2)rs编号
rs_id
、P值pval_nominal
、case在所有样本中的比例s
3.2 导入表型2(eQTL)数据:
eqtl <- read.table(file="E:/path_to_eqtl/eQTL.txt", header=T);
eQTL数据包括:rs编号rs_id
,基因gene_id
,次等位基因频率maf
、P值pval_nominal
等:
注意事项:如果表型是连续型变量,输入文件三选一:
1)rs编号
rs_id
、P值pval_nominal
、表型的标准差sdY
;2)rs编号
rs_id
、P值pval_nominal
、效应值beta
,效应值方差varbeta
, 样本量N
,次等位基因频率MAF
;3)rs编号
rs_id
、P值pval_nominal
、次等位基因频率MAF
;
3.3 合并GWAS和eQTL数据:
input <- merge(eqtl, gwas, by="rs_id", all=FALSE, suffixes=c("_eqtl","_gwas"))
head(input)
3.4 共定位分析
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)
dataset1的type="cc"指的是GWAS的表型是二分类(case和control);
dataset2的type="quant"指的是eQTL的表型(基因表达量)是连续型
N指样本量;
3.5 筛选共定位的位点
通常情况下,很多文献认为PPA > 0.95的位点是共定位位点,也有一些文献会放松要求到0.75。
这里假定后验概率大于0.95为共定位位点:
library(dplyr)
need_result=result$results %>% filter(SNP.PP.H4 > 0.95)
结果如下:
从图上可以看出,SNP.4811位点的后验概率为1。怎么找到这个位点呢?可以通过对应的P值(1.81e-42)找到这个位点的rs号;
4 结果解读
对于输出结果,我们只需要关注最后一列信息SNP.PP.H4
,对应推文前面提到的第四种设想 H4。
SNP.PP.H4
表示的是GWAS显著信号和eQTL位点为同一个位点的后验概率,范围在0-1之间,0表示概率为0%,1表示概率为100%。后验概率越高越好。
5 注意事项
1)由于共定位分析是基于某个基因组区域进行计算,所以请不要把全基因组的信息都丢进去(偷懒该打),这么做一个是没意义,另外一个特别耗时,给计算机增加工作负担;
2)虽然我们没必要把基因组的全部信息丢进去,但也不意味着只放一个变异位点信息就行。
3)因此,正确的做法是,先提取GWAS中显著变异位点上下游区域(这个区域多大自己定,没有金标准)的GWAS summary数据,理想情况下,提取后显著变异位点所在基因组区域的SNP数量在1,000-10,000之间。
4)该方法的设想是只有一个causal 位点,所以如果表型1(GWAS)和 表型2 (以eQTL为例)在某个区域有多个显著位点的话,用该方法是定位不出结果的,这是该方法的局限,所以如果某个基因组区域存在多个显著位点,请考虑其他工具(允许多个causal 位点共定位的工具)
特别鸣谢:
https://chr1swallace.github.io/coloc/FAQ.html
https://hanruizhang.github.io/GWAS-eQTL-Colocalization/
https://rpubs.com/Charleen_D_Adams/743547
致谢橙子牛奶糖(陈文燕),请用参考模版:We thank the blogger (orange_milk_sugar, Wenyan Chen) for XXX
感谢小可爱们多年来的陪伴, 我与你们一起成长~