首页 > 其他分享> 文章详细

post-GWAS：使用coloc进行共定位分析(Colocalization)

2021-07-23 20:31:37 阅读：364 来源： 互联网

标签：定位 GWAS rs Colocalization eQTL 表型 coloc 位点

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

感谢小可爱们多年来的陪伴，我与你们一起成长~

标签：定位,GWAS,rs,Colocalization,eQTL,表型,coloc,位点
来源： https://blog.csdn.net/wendy_milk/article/details/119043642

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

ICode9