使用coloc 进行 QTL 共定位Colocalization_gwas eqtl coloc

GWAS找到显著信号位点后，需要解释显著信号位点如何影响表型。 常见的一个解释方法是共定位分析。

主流的共定位分析包括：

1）GWAS和eQTL共定位；

2）GWAS和sQTL共定位；

3）GWAS和meQTL共定位；

4）GWAS和pQTL共定位；

其中，GWAS和eQTL共定位应用最为广泛。

共定位分析旨在确定两个性状在给定基因组区域中可能共享的因果变异，本文中所说的共定位是基于贝叶斯推断的共定位分析，使用的软件是 R 语言中的 coloc ​包。

install.packages("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)

抛砖引玉-共定位的原理与算法

官方对于 coloc ​的介绍如下：

The coloc package can be used to perform genetic colocalisation analysis of two potentially related phenotypes, to ask whether they share common genetic causal variant(s) in a given region

因此，共定位的目的是为了检验输入的两种表型在给定的区域内是否共享同一个因果变异，本文中的共定位是基于贝叶斯推断的共定位，关于贝叶斯推断的原理，请参考以下的资料：

1. 文献 1：Bayesian Test for Colocalisation between Pairs of Genetic Association Studies Using Summary Statistics
2. 文献 2：Eliciting priors and relaxing the single causal variant assumption in colocalization analyses
3. coloc 官方文档

具体来说，当检测到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进行共定位分析。

共定位分析的准备工作

相信你在阅读官方文档时注意到了这段话：

A single genomic region does not correspond to the whole genome, nor to a whole chromosome. Coloc also does not automatically split chromosomes or a genome into regions. It is assumed the user can look at their data, identify a region with overlapping GWAS signals between two studies, and decide on the set of SNPs to include.

意思是软件不会为你自动选择共定位区域，而是把共定位区域的选择工作交给用户来完成，而且强调了共定位区域不能为整条染色体或者全基因组。因此，定义共定位区域是共定位分析中特别重要的一部分，而共定位区域是为了共定位数据对（也就是你要检验的两种共定位表型）而服务的，所以一般的流程是先确定要检验的表型，再确定要检验的区域。

待检验表型的选择

所谓待检验的表型，指的是最后用来计算共定位概率的配对数据。

在 GWAS 分析中，一般只针对一个性状进行关联分析，

而在 QTL 分析中，往往同时对很多表型进行关联，

        例如在 eQTL 分析中，每一个基因都代表一个表型，

        在 caQTL 中，每一个开放区域都代表一个表型。

此时，我们检验的表型就是每一个基因-开放区域配对数据，因此，我们需要首先确定所有的配对数据，然后为它们分别指定共定位区域。

在这里，我参考了 coloc ​官网的推荐文献，进行了简单的整理，感兴趣的学者可以阅读一下原文进行细节的探究：

这些文献中，很多都提到了一个概念独立位点，独立位点指的是一个区域中，没有其他 SNP 与此 SNP 的连锁程度超过给定阈值，这个阈值一般是 �2<0.2 ，也有一些研究使用了 �2<0.1 ，这部分工作可以使用 plink 完成。

共定位区域的选择

选择共定位区域是为了根据先验知识来计算上一步获得的所有数据对的共定位后验概率，计算的时候会根据共定位区域内所有的 SNP 的效应值、MAF 等统计量计算共定位的概率。因此共定位区域的选择往往依赖于共定位数据对的选择。上面的表格中也列出了其他研究中使用的共定位区域。coloc ​官方网站给出了两种方案，分别是基于 tagSNP 和基于遗传距离，这里的遗传距离可以自己根据实际情况指定。

准备 coloc 需要的数据格式

首先，我们先看一下 coloc ​需要的数据有哪些：

##  $ beta    : Named num [1:50] 0.288 0.3 0.334 0.444 0.494 ...
##   ..- attr(*, "names")= chr [1:50] "s1" "s2" "s3" "s4" ...
##  $ varbeta : Named num [1:50] 0.00681 0.0105 0.00733 0.00591 0.01514 ...
##   ..- attr(*, "names")= chr [1:50] "s1" "s2" "s3" "s4" ...
##  $ snp     : chr [1:50] "s1" "s2" "s3" "s4" ...
##  $ position: int [1:50] 1 2 3 4 5 6 7 8 9 10 ...
##  $ type    : chr "quant"
##  $ sdY     : num 1.12

总结一下，必要的信息有以下几类：

1. SNP 的基础信息，包括 SNP 的 ID（不一定是 rsID）和 SNP 位置
2. 关联分析的效应信息，包括 beta 值和效应方差方差 varbeta，如果没有这一项，就需要有 P、MAF 和 N
3. sdY，Y 的标准差，如果没有这一项，就需要下面两项：
   1. MAF，次等位基因频率
   2. N，样本量
4. type，分析的类型，有 quant ​和 cc ​两种，分别代表数量性状关联和 Case/Control 分析

不同的关联分析软件对于统计信息有着不同的命名，请根据你使用的分析软件的结果进行调整和计算，我使用的数据是用 MatrixEQTL 生成的，它的结果包含以下几列：

• p：未校正的 p 值
• fdr：校正后的 p 值
• beta：效应值，也就是线性回归的斜率
• t−statistic：T 检验的统计量

可以看出，默认结果中不包含 varbeta，因此需要我们自己进行计算，我们采用以下公式计算 varbeta：

�������=(������)2

最后，需要注意的是，coloc​ 接受的数据格式是列表，而不是数据框，而且 type​ 和 sdY​ 需要在转换成列表后再指定，也就是说，type​ 只需要指定一个值，具体格式参考上面的官方文件格式实例。

进行共定位分析

生成 coloc 需要的数据格式后，我们就可以进行共定位分析了，上面的部分已经说过了，共定位分析是以共定位数据对为单位的，也就是说，每一个共定位数据对都对应一个共定位区域，每一个共定位区域都对应两种表型在共定位区域内的所有 SNP。简单来说，有几个数据对，就有几个共定位区域，就要分别从两种性状的表型中提取几次 SNP，这个过程可以通过循环实现，也可以一次性把所有的数据提取到列表里，逐个进行分析。我采用的是后者。

共定位思想

我们将要使用的数据是 eQTL 和 meQTL 的共定位，我们的研究思路是：

1. 在 eQTL 中找出每一个基因的 lead-eSNP
2. 找出与这个 lead-eSNP 重合的所有 meSNP，这些 meSNP 可能对应着不同的甲基化探针
3. 找出这些甲基化探针对应的 lead-meSNP
4. 计算第三步中的 lead-meSNP 与第一步中的 lead-eSNP 的连锁强度（ �2 ），选择连锁强度最高的 lead-meSNP 对应的甲基化探针
5. 使用每一个基因的上下游各 1mb 范围作为共定位区域进行分析

下载练习数据

下面，我将使用公共数据库中的数据进行共定位，所使用的数据分别是来自 PancanQTL[1] 的 cis-eQTL 数据和来自 Pancan-meQTL[2] 的 cis-meQTL 数据，这里使用乳腺癌的数据，如果想根据这篇教程练习，请先下载相关数据，下载地址如下：

• eQTL 下载地址
• meQTL 下载地址

处理原始数据

推荐大家使用 data.table ​进行数据的读取，使用 tidyverse ​进行数据处理，代码简洁优雅，功能强大。假设大家已经把数据导入，并命名为 eQTL ​和 meQTL​，下面的分析都以此为基础进行分析。

suppressMessages(library(tidyverse))
suppressMessages(library(data.table))

数据预处理

此处的 MAF 是用 TCGA 的数据进行计算的，这里提供下载，下载地址

# 导入MAF信息
maf = fread("BRCA_MAF.txt")
meQTL = fread("BRCA_tumor.cis_meQTL.tsv") %>%
    rename(
        pvalue   = `p-value`,
        `t-stat` = status
    ) %>%
    mutate(varbeta = (beta / `t-stat`) ^ 2) %>%
    separate(
        col    = "alleles",
        into   = c("ref", "alt"),
        sep    = "/",
        fill   = "right",
        remove = TRUE
    ) %>%
    separate(
        col    = "snp_position",
        into   = c("chr", "position"),
        sep    = ":",
        fill   = "right",
        remove = TRUE
    ) %>%
    mutate_at(
        "position",
        as.numeric
    ) %>%
    left_join(
        y = maf,
        by = c("snp", "position", "ref", "alt")
    ) %>%
    select(snp, chr, position, ref, alt, maf, probes, beta, varbeta, pvalue)
eQTL = fread("BRCA_cis_eQTL_fdr005.tsv") %>%
    rename(pvalue  = `p-value`) %>%
    mutate(varbeta = (beta / `t-stat`) ^ 2) %>%
    separate(
        col    = "alleles",
        into   = c("ref", "alt"),
        sep    = "/",
        fill   = "right",
        remove = TRUE
    ) %>%
    left_join(
        y  = maf,
        by = c("snp", "position", "ref", "alt")
    ) %>%
    select(snp, chr, position, ref, alt, maf, gene, beta, varbeta, pvalue)

提取 leadQTL

lead_eSNP = eQTL %>%
    group_by(gene) %>%
    arrange(pvalue) %>%
    distinct(gene, .keep_all = TRUE) %>%
    rename_at(
        vars(c("beta", "varbeta", "pvalue")))
​
lead_meSNP = meQTL %>%
    group_by(probes) %>%
    arrange(pvalue) %>%
    distinct(probes, .keep_all = TRUE)

定义共定位数据对

coloc_pair_list = apply(
    lead_eSNP %>%
        mutate(gene2 = gene) %>%
        column_to_rownames("gene2") %>%
        mutate_at(
            .vars = vars(c("position", ends_with("_eqtl"))),
            .funs = as.numeric
        ),
    MARGIN = 1,
    FUN = function(x) {
        result                   = as.list(x)
        result[["maf"]]          = as.numeric(result[["maf"]])
        result[["position"]]     = as.numeric(result[["position"]])
        result[["beta_eqtl"]]    = as.numeric(result[["beta_eqtl"]])
        result[["varbeta_eqtl"]] = as.numeric(result[["varbeta_eqtl"]])
        result[["pvalue_eqtl"]]  = as.numeric(result[["pvalue_eqtl"]])
        return(result)
    },
    simplify = FALSE
)

找出 meQTL 中与每一个 lead-eSNP 重合的 meSNP

overlapped_meSNP = lapply(
    X = coloc_pair_list,
    FUN = function(x) {
        meQTL %>% filter(snp %in% x[["snp"]])
    }
)

计算连锁强度并确定最高连锁强度的甲基化探针

这一步明白原理就可以了，计算 LD 的方法有很多种，这里选择了 LDlinkR​ 软件包，这个工具的在线网页服务地址：LDlink。下面的代码会将连锁强度最高的甲基化探针以及对应的统计量信息输出。

lead_probe = mapply(
    FUN = function(coloc_list, meqtl_overlap) {
        # 没有重叠的情况直接输出空结果
        if(length(meqtl_overlap[["snp"]]) == 0) {
            return(list(
                probe        = NA,
                beta_meqtl   = NA,
                varbeta      = NA,
                pvalue_meqtl = NA
            ))
        }
        # 多个meSNP对应到同一个探针时，直接输出探针
        if(length(unique(meqtl_overlap[["probes"]])) == 1) {
            return(
                list(
                    probe         = meqtl_overlap[["probes"]][1],
                    beta_meqtl    = meqtl_overlap[["beta"]][1],
                    varbeta_meqtl = meqtl_overlap[["varbeta"]][1],
                    pvalue_meqtl  = meqtl_overlap[["pvalue"]][1]
                )
            )
        } else {
            lead_esnp = coloc_list[["snp"]]
            # 获得重合meSNP对应的探针的lead—meSNP用来计算LD
            lead_mesnps_query = unique(lead_meSNP$snp[lead_meSNP$probes %in% meqtl_overlap$probes])
            ld_with_lead_esnp = sapply(
                X = lead_mesnps_query,
                FUN = function(x) {
                    if (length(x) == 0) {
                        return(0)
                    } else {
                        tryCatch({
                            ld_matrix = LDlinkR::LDpair(
                                var1  = x,
                                var2  = lead_esnp,
                                pop   = "CEU",
                                # 这里的token需要自己申请
                                token = "自己申请"
                            )
                            return(ld_matrix[["r2"]][1])
                        },
                        error = function(error) {
                            print(paste(x, lead_esnp, error, sep = "\t"))
                            return(0)
                        })
                    }
                },
                simplify = "array"
            )
            max_ld = max(ld_with_lead_esnp)
            if(max_ld == 0 | is.na(max_ld)) {
                return(list(
                    probe        = NA,
                    beta_meqtl   = NA,
                    varbeta      = NA,
                    pvalue_meqtl = NA
                ))
            } else {
                max_ld_snp = names(ld_with_lead_esnp)[which(ld_with_lead_esnp == max_ld)]
                # 如果多个probe都有最大连锁，取最显著的
                max_ld_probe   = (lead_meSNP %>% filter(snp %in% max_ld_snp) %>% arrange(pvalue))$probes[1]
                max_ld_beta    = (lead_meSNP %>% filter(snp %in% max_ld_snp) %>% arrange(pvalue))$beta[1]
                max_ld_varbeta = (lead_meSNP %>% filter(snp %in% max_ld_snp) %>% arrange(pvalue))$varbeta[1]
                max_ld_pvalue  = (lead_meSNP %>% filter(snp %in% max_ld_snp) %>% arrange(pvalue))$pvalue[1]
                return(
                    list(
                        probe        = max_ld_probe,
                        beta_meqtl   = max_ld_beta,
                        varbeta      = max_ld_varbeta,
                        pvalue_meqtl = max_ld_pvalue
                    )
                )
            }
        }
    },
    coloc_pair_list,
    overlapped_meSNP,
    SIMPLIFY = FALSE
)

提取共定位区域内的所有 SNP

我们定义共定位区域为每一个基因的顺式关联窗口，也就是上下游各 1mb，由于我们使用的 eQTL 的关联距离也是 1mb，因此我们直接提取与每一个基因关联的所有 cis-eQTL，它们都位于我们的共定位区域中。

# 提取共定位区域内的所有eSNP
eSNP_by_gene = sapply(
    X = unique(eQTL[["gene"]]),
    FUN = function(egene) {
        eQTL %>% filter(gene == egene)
    },
    simplify = FALSE,
    USE.NAMES = TRUE
)
​
# 提取共定位区域内的所有meSNP
meSNP_by_gene = lapply(
    X = eSNP_by_gene,
    FUN = function(esnps) {
        meQTL %>% filter(snp %in% esnps[["snp"]])
    }
)

进行共定位分析

共定位的时候需要注意，有很多可能导致数据出问题的情况，比如没有与 eQTL 重叠的 meQTL，所使用的 SNP 不是单方向突变而无法计算 LD 等等，在写函数的时候需要设计一个异常捕获来处理这些情况。

coloc_result = mapply(
    FUN = function(df_1, df_2, n_1, n_2, type_1 = "quant", type_2 = "quant") {
    # 如果没有重叠的SNP，就跳过共定位
        if (nrow(df_1) == 0 | nrow(df_2) == 0) {
            return(
                list(
                    n_snps = 0,
                    PP3    = 0,
                    PP4    = 0
                )
            )
        }
        df_1 = df_1 %>%
            rename(MAF = maf, p = pvalue) %>%
            arrange(p) %>%
            # coloc要求不能有重复的SNP，所以只保留更显著的
            distinct(snp, .keep_all = TRUE) %>%
            select(snp, position, p, beta, varbeta, MAF) %>%
            filter(!is.na(MAF)) %>%
            as.list()
        df_1[["N"]] = n_1
        df_1[["type"]] = type_1
 
        df_2 = df_2 %>%
            rename(MAF = maf, p = pvalue) %>%
            arrange(p) %>%
            distinct(snp, .keep_all = TRUE) %>%
            select(snp, position, p, beta, varbeta, MAF) %>%
            filter(!is.na(MAF)) %>%
            as.list()
        df_2[["N"]] = n_2
        df_2[["type"]] = type_2
 
        if (length(df_1[["snp"]])== 0 | length(df_2[["snp"]]) == 0) {
            return(
                list(
                    n_snps = 0,
                    PP3    = 0,
                    PP4    = 0
                )
            )
        }
 
        if (is.null(coloc::check_dataset(df_1)) & is.null(coloc::check_dataset(df_2))) {
            invisible(capture.output({
                coloc_result = coloc::coloc.abf(
                    dataset1 = df_1,
                    dataset2 = df_2
                )
            }))
            return(
                list(
                    n_snps = coloc_result[["summary"]][["nsnps"]],
                    PP3    = coloc_result[["summary"]][["PP.H3.abf"]],
                    PP4    = coloc_result[["summary"]][["PP.H4.abf"]]
                )
            )
        } else {
            return(
                list(
                    n_snps = 0,
                    PP3    = 0,
                    PP4    = 0
                )
            )
        }
    },
    df_1      = eSNP_by_gene,
    df_2      = meSNP_by_gene,
    n_1       = 1092,
    n_2       = 664,
    SIMPLIFY  = FALSE,
    USE.NAMES = TRUE
)

合并所有数据

在提取 lead-eSNP 的过程中，我们获得了共定位的基因，然后我们通过计算连锁强度获得了共定位的甲基化探针，最后我们进行了共定位检验，现在我们将这些信息合并到一个列表中，最终生成一个数据表方便输出，现在我们有三个列表，分别保存着 SNP 的信息与 eQTL 的基因信息、甲基化探针信息、共定位结果

# 把所有的必要信息组合起来
final_coloc_result_list = mapply(
    FUN = function(coloc_pairs_eqtl, coloc_pairs_meqtl, coloc_result) {
        return(c(
            coloc_pairs_eqtl,
            coloc_pairs_meqtl,
            coloc_result
        ))
    },
    coloc_pairs_eqtl  = coloc_pair_list,
    coloc_pairs_meqtl = lead_probe,
    coloc_result      = coloc_result,
    SIMPLIFY          = FALSE
)
 
# 将列表合并成数据框
final_coloc_result_table = do.call(rbind, final_coloc_result_list) %>%
    as.data.frame() %>%
    # 转换后数据类型全部丢失，需要手动设置回来
    mutate_at(
        .vars = vars(c("snp", "chr", "ref", "alt", "gene", "probe")),
        .funs = as.character
    ) %>%
    mutate_at(
        .vars = vars(-c("snp", "chr", "ref", "alt", "gene", "probe")),
        .funs = as.numeric
    ) %>%
    mutate_all(
        .funs = function(x) {
            ifelse(is.na(x) | x == "NA", NA, x)
        }
    ) %>%
    # 如果探针缺失，则共定位信号为0
    mutate(PP4 = ifelse(is.na(probe), 0, PP4)) %>%
    arrange(desc(PP4))

合并后的部分共定位结果如下所示：

参考

1. ^Jing Gong, Shufang Mei, Chunjie Liu, Yu Xiang, Youqiong Ye, Zhao Zhang, Jing Feng, Renyan Liu, Lixia Diao, An-Yuan Guo, Xiaoping Miao, Leng Han, PancanQTL: systematic identification of cis-eQTLs and trans-eQTLs in 33 cancer types, Nucleic Acids Research, Volume 46, Issue D1, 4 January 2018, Pages D971–D976 https://doi.org/10.1093/nar/gkx861IF: 14.9 Q1
2. ^Jing Gong, Hao Wan, Shufang Mei, Hang Ruan, Zhao Zhang, Chunjie Liu, An-Yuan Guo, Lixia Diao, Xiaoping Miao, Leng Han, Pancan-meQTL: a database to systematically evaluate the effects of genetic variants on methylation in human cancer, Nucleic Acids Research, Volume 47, Issue D1, 08 January 2019, Pages D1066–D1072 https://doi.org/10.1093/nar/gky814IF: 14.9 Q1

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 位点共定位的工具）