方法学|共定位分析(Colocalization analysis)

昵称69125444 2023-09-26 发布于广西

展开全文

https://chr1swallace./coloc/index.html

01.共定位分析的目的

共定位分析通常用于鉴定两个表型是否由某一区域的同一个causal variant驱动，从而加强两个表型之前的关联证据。这些表型可以是分子表型(如蛋白质，mRNA等)，也可以是常见的连续型性状(如血压，BMI等)或二分类疾病(如糖尿病，AD等)。

02.共定位分析的假设

在给定区域中，共定位分析的前提假定是两个性状中的每一个性状在该区域最多有1个真正的causal variant，由此产生了五个互斥的模型假定(H0-H4)，这五个模型假定是前提假定下的所有可能的关联情况

共定位分析的过程中，会为上面每个模型产生后验概率(PP.H0-PP.H4)，五个模型后验概率总和为1，当某一模型的后验概率越高，说明对应的模型假定在给定数据的情况下更有可能成立，当然在一般的分析中，我们更希望H4假定成立，因为H4模型假定表示两个性状由同一个causal variant驱动。一般当PP.H4>0.75时(这个值可视情况调整)，我们认为H4模型假定成立。

H4模型假定成立，说明该区域内存在一个variant是两个性状共享的，但具体是哪个variant并不清楚，所以共定位分析除了给一个区域的五个模型假定计算后验概率外，它也对区域内的每个SNP计算SNP的后验概率(SNP.PP.H4)，以评估哪个variant是最有可能的causal variant。

03.共定位分析的数据要求

For summary data（eQTL，pQTL，GWAS数据等）:

二分类表型需要SNP，CHR，BP，A1，A2，BETA，VAR

连续型表型需要SNP，CHR，BP，A1，A2，BETA，VAR，MAF，N

注意：

1.共定位分析是在一个基因组区域进行，所以至少一个数据里有CHR，BP，方便提取区域的SNP

2.共定位需要的是方差，而不是标准差，方差等于SE的平方，VAR = SE^2

3.共定位分析要求不能有重复SNP，并且SNP在两个数据中不能有缺失值

4.数据最好对齐效应等位基因

关于区域的选择：一般这个区域在1MB左右即可，这样保证区域内有几千个SNP可以用来计算，这个区域的选择一般是包含该区域在两个性状中的显著信号即可，如果是QTL数据和GWAS数据进行共定位，一般这个区域由QTL的cis-region决定，如果是对两个性状的GWAS数据进行共定位，一般是根据对GWAS数据clump后的某一区域的top SNP加减500kb决定。

04.共定位分析的代码演示

##################两个二分类表型#注意，此代码仅为示例#################rm(list=ls())library(dplyr)library(coloc)setwd('E:/project/')#读取数据data1 = data.table::fread('./1.DATA/AD',sep = '\t',header=T)data2 = data.table::fread('./1.DATA/PD',sep = '\t',header=T)

#提取一个区域data1 = data1 %>% filter(CHR==9,BP >= 4981602-500000,BP <= 4981602 500000)data2 = data2 %>% filter(CHR==9,BP >= 4981602-500000,BP <= 4981602 500000)

#合并并去重data = merge(data1,data2,by='SNP')data = data[!duplicated(data$SNP),]

#对齐效应等位基因data = data %>% filter((A1.x==A1.y&A2.x==A2.y)|(A1.x==A2.y&A2.x==A1.y)) data = data %>% mutate(BETA.y = ifelse(A1.x==A1.y,BETA.y,-BETA.y))

#计算方差VARdata$VAR.x = data$SE.x^2data$VAR.y = data$SE.y^2data = data[data$VAR.x!=0 & data$VAR.y!=0 ,]

#拆分整理data1 = data[,c('BETA.x','VAR.x','SNP')]data2 = data[,c('BETA.y','VAR.y','SNP')]colnames(data1)=c('beta','varbeta','snp')colnames(data2)=c('beta','varbeta','snp')data1 = as.list(data1)data2 = as.list(data2)

#声明表型类型,二分类表型'cc',连续型表型'quant'data1$type = 'cc'data2$type = 'cc'

#coloc分析，p1,p2,p12为先验概率参数，下面的值是默认参数res = coloc.abf(data1,data2,p1=1e-4,p2=1e-4,p12=1e-5)

##################一个连续一个二分类#注意，此代码仅为示例#################rm(list=ls())library(dplyr)library(coloc)setwd('E:/project/')#读取数据data1 = data.table::fread('./1.DATA/AD',sep = '\t',header=T)data2 = data.table::fread('./1.DATA/BMI',sep = '\t',header=T)
#提取一个区域data1 = data1 %>% filter(CHR==9,BP >= 4981602-500000,BP <= 4981602 500000)data2 = data2 %>% filter(CHR==9,BP >= 4981602-500000,BP <= 4981602 500000)
#合并并去重data = merge(data1,data2,by='SNP')data = data[!duplicated(data$SNP),]
#对齐效应等位基因data = data %>% filter((A1.x==A1.y&A2.x==A2.y)|(A1.x==A2.y&A2.x==A1.y)) data = data %>% mutate(BETA.y = ifelse(A1.x==A1.y,BETA.y,-BETA.y))
#计算方差VARdata$VAR.x = data$SE.x^2data$VAR.y = data$SE.y^2data = data[data$VAR.x!=0 & data$VAR.y!=0 ,]
#拆分整理data1 = data[,c('BETA.x','VAR.x','SNP')]data2 = data[,c('BETA.y','VAR.y','SNP','MAF','N')]colnames(data1)=c('beta','varbeta','snp')colnames(data2)=c('beta','varbeta','snp','MAF','N')data1 = as.list(data1)data2 = as.list(data2)
#声明表型类型,二分类表型'cc',连续型表型'quant'data1$type = 'cc'data2$type = 'quant'
#coloc分析，p1,p2,p12为先验概率参数，下面的值是默认参数res = coloc.abf(data1,data2,p1=1e-4,p2=1e-4,p12=1e-5)

04.共定位分析的结果

上面代码中的res变量保存了共定位分析的结果

res$summary保存了该区域的五个模型的后验概率，PP.H4.abf>0.75，认为两个性状在该区域存在共定位的证据，即该区域存在共享的causal variant。

res$results保存了该区域每个SNP的SNP.PP.H4，区域中SNP.PP.H4最大的SNP是最有可能共享的那个causal variant。我们可以看到rs4603的SNP.PP.H4=0.906，说明rs4603是最有可能的causal variant。

05.该算法的缺陷及扩展

首先，coloc假定区域里一个性状的causal variant最多有一个，这并不在任何情况下都成立，不过，coloc最近加入的SuSiE算法放宽了这一假定，SuSiE适用于区域中存在多个causal variants的情况。
其次，coloc只能针对两个性状，而不是多个性状，不过，想同时对多个性状做共定位分析，可以考虑其他算法，如HyprColoc,moloc等。

SuSiE:https://chr1swallace./coloc/articles/a06_SuSiE.html

HyprColoc:https://github.com/cnfoley/hyprcoloc

moloc:https://github.com/clagiamba/moloc