共定位分析是一种用贝叶斯检验来分析两种可能相关的表型是否在给定区域内共享共同的遗传因果变异的方法，本文中，我们使用 coloc 包来进行共定位分析。
数据的准备
MAF的计算
我们的共定位的目的是为了找出caQTL、eQTL与meQTL的共定位，主角是我们的caQTL数据。coloc需要SNP的次等位基因频率（MAF）以及样本量信息，这两部分信息需要从基因型文件中提取，因此，我们的第一步就是从基因型数据中提取caQTL中的caSNP的次等位基因频率和样本量信息，由于基因型矩阵分别使用0、1、2对AA、Aa和aa进行了编码，因此对于较小频率的突变等位基因a，其等位基因频率为：
\[AF_{a} = \frac{0N_{AA} + 1N_{Aa} + 2N_{aa}}{2N_{AA} + 2N_{Aa} + 2N_{aa}}，\ 其中N_{AA}、N_{Aa}、N_{aa}分别代表三种基因型的样本数 \]
上面的这个公式实质上等价于基因型矩阵每一行的均值的一半，因此我们可以使用 mean(row) / 2 来计算每一个SNP的alt等位基因频率，然后再判断此频率是否在5%到95%之间。
library("readr")
library("tidyverse")
# 导入原始的caQTL文件
caqtl_raw = read_tsv("caqtl_raw.tsv") %>% distinct(snps, .keep_all = TRUE)
# 导入eQTL文件
eqtl = read_tsv("eqtl.tsv", col_types = cols())
eqtl_unique = eqtl %>% arrange(`p-value`) %>% distinct(snp, .keep_all = TRUE)
# 计算等位基因频率的函数
calculate_AF = function(row) {
x = as.numeric(row)
AF = round(mean(x[!is.na(x)]) / 2, 3)
return(min(AF, 1-AF))
}
# 导入基因型文件并计算MAF
genotype = read_tsv("BRCA_genotype")
genotype$maf = apply(genotype, 1, calculate_AF)
# 筛选caQTL和eQTL中出现的SNP
genotype_caqtl = genotype %>% filter(ID %in% caqtl_raw$snps)
genotype_caqtl = genotype_caqtl[match(caqtl_raw$snps, genotype_caqtl$ID),]
genotype_caqtl$maf = apply(genotype_caqtl, 1, calculate_AF)
genotype_eqtl = genotype %>%
mutate(ID = gsub(":.*", "", ID)) %>%
filter(ID %in% intersect(ID, eqtl_unique$snp)) %>%
distinct(ID, .keep_all = TRUE)
eqtl_unique = eqtl_unique %>% filter(snp %in% intersect(snp, genotype_eqtl$ID))
genotype_eqtl = genotype_eqtl[match(eqtl_unique$snp, genotype_eqtl$ID),]
genotype_eqtl$maf = apply(genotype_eqtl, 1, calculate_AF)
处理两种QTL的数据为coloc需要的格式
coloc不接受数据表格式，因此需要将数据包格式的数据转换为列表的格式，对于定量数据，需要以列表的形式指定额外的属性：
type: "quant"：数据的类型MAF：SNP的次等位基因频率N：数据的样本量
数据的处理方式如下：
# 处理caQTL数据为coloc需要的格式
caqtl = read_tsv("caQTL.tsv", col_types = cols())
caqtl_clean = caqtl %>%
select(caSNP_ID,caSNP_chr,caSNP_pos,std_error,beta,pvalue) %>%
mutate(varbeta = std_error^2) %>%
select(-std_error) %>%
rename("snp" = caSNP_ID, "p" = pvalue) %>%
unite(col = "position", caSNP_chr,caSNP_pos,sep = ":") %>%
arrange(p) %>%
distinct(snp, .keep_all = TRUE) %>%
as.list()
# 处理eQTL数据为coloc需要的格式
eqtl_clean = eqtl_unique %>%
rename("tstat" = "t-stat", "p" = "p-value") %>%
mutate(varbeta = (beta / tstat) ^ 2) %>%
select(snp, chr, position, beta, varbeta, p) %>%
unite(col = "position", chr, position, sep = ":") %>%
arrange(p) %>%
distinct(snp, .keep_all = TRUE) %>%
as.list()
# 转换为列表格式并加入上一部分的结果
caqtl_clean$type = "quant"
caqtl_clean$MAF = genotype_caqtl$maf
caqtl_clean$N = 74
eqtl_clean$type = "quant"
eqtl_clean$MAF = genotype_eqtl$maf
eqtl_clean$N = 1093
# 检查数据表是否有问题
check_dataset(caqtl_clean) # NULL
check_dataset(eqtl_clean) # NULL
coloc分析以及结果解读
如果数据的检查函数返回 NULL 说明数据没有问题，接着就可以进行共定位分析了，使用的函数是 coloc.abf
coloc_results = coloc.abf(eqtl_clean, caqtl_clean)
根据共定位算法，结果中的后验概率类型包含以下5种，分别是：
H0：两种性状在区域内都没有遗传学关联H1：只有性状1在区域内具有遗传学关联H2：只有性状2在区域内具有遗传学关联H3：两种性状在区域内具有不同的遗传学关联H4：两种性状在区域内具有同一的遗传学关联（共定位证据）
我们选取后验概率比较高的点进行查看
coloc_results$results %>% filter(SNP.PP.H4 > 0.75) %>% select(snp)
发现只有 rs2031947 一个结果，其余的SNP都不显著，我们在caQTL和eQTL数据库中进行寻找，结果如下：
可以看到caQTL和eQTL的趋势比较一致，而与meQTL相反，而且这个SNP对应的基因 GSTM1 曾被预测与多种癌症相关。