摘要#
共定位分析旨在确定两个性状在给定基因组区域中可能共享的因果变异，本文中所说的共定位是基于贝叶斯推断的共定位分析，使用的软件是R语言中的coloc包。
抛砖引玉-共定位的原理与算法#
官方对于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
因此，共定位的目的是为了检验输入的两种表型在给定的区域内是否共享同一个因果变异，本文中的共定位是基于贝叶斯推断的共定位，关于贝叶斯推断的原理，请参考以下的资料：
共定位分析的准备工作#
相信你在阅读官方文档时注意到了这段话：
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官网的推荐文献，进行了简单的整理，感兴趣的学者可以阅读一下原文进行细节的探究：
| 共定位类型 | 文章地址 | 共定位区域 |
|---|---|---|
| eQTL-meQTL | Nature communications 2018 | leadSNP 上下游 250kb |
| eQTL-GWAS | Nature 2017 | 独立 GWAS 上下 1mb |
| eQTL-pQTL | Nature 2018 | 按照遗传距离划分区域 |
| meQTL-GWAS | AJRCCM 2018 | 甲基化位点上下游 500kb |
| pQTL-eQTL | Nature Communications 2018 | lead-pSNP 上下 1mb |
| caQTL-GWAS/eQTL | Nature Communications 2018 | 其他研究确定的区域 |
| GWAS-GWAS | Inflammatory Bowel Diseases 2018 | 每一个 SNP 的上下 50kb |
这些文献中，很多都提到了一个概念独立位点，独立位点指的是一个区域中，没有其他SNP与此SNP的连锁程度超过给定阈值，这个阈值一般是r2<0.2r2<0.2），选择连锁强度最高的lead-meSNP对应的甲基化探针
- 使用每一个基因的上下游各1mb范围作为共定位区域进行分析
下载练习数据#
下面，我将使用公共数据库中的数据进行共定位，所使用的数据分别是来自PancanQTL[1]的cis-eQTL数据和来自Pancan-meQTL[2]的cis-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软件包，这个工具的在线网页服务地址见我的友情链接。下面的代码会将连锁强度最高的甲基化探针以及对应的统计量信息输出。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))合并后的部分共定位结果如下所示：
snp chr position ref alt gene beta_eqtl varbeta_eqtl maf pvalue_eqtl probe beta_meqtl varbeta_meqtl pvalue_meqtl n_snps PP3 PP4 rs2239961 chr22 21363960 A G FLJ39582 0.67 0.0015315885 0.2202381 3.21e-58 cg17353431 -0.82 0.0010622003280752 1.4e-97 4 0 1 rs9896243 chr17 44826056 C G LRRC37A2 0.7 0.0017976367 0.15888278 1.01e-54 cg01570182 0.86 0.00381148899043469 1.01e-38 2 0 1 rs4982912 chr14 24903284 C T CBLN3 -0.5 0.001231148 0.31730769 2.57e-42 cg13105904 0.64 0.00173385275543321 1.42e-45 3 0 1 rs9897978 chr17 13900256 G T CDRT15P 0.49 0.0012662641 0.34798535 7.77e-40 cg11395062 -0.49 0.00200979534574591 1.25e-25 1 0 1 rs11868472 chr17 74701165 A C MXRA7 0.34 0.0007779477 0.43543956 4.19e-32 cg27546012 -0.42 0.0014711953462188 1.1e-25 1 0 1 rs4751321 chr10 132897429 A C TCERG1L -0.45 0.0015664051 0.24679487 2.39e-28 cg09858951 -0.45 0.0021041999831664 2.85e-21 1 0 1 rs16831518 chr1 160146645 C T ATP1A4 0.38 0.0016033737 0.19413919 1.45e-20 cg19308497 -0.51 0.00195311913350895 3.91e-28 3 4.70379291075399e-15 1 rs9896243 chr17 44826056 C G LRRC37A 0.38 0.0018860408 0.15888278 8.53e-18 cg01570182 0.86 0.00381148899043469 1.01e-38 2 1.47792889038308e-15 0.999999999999929 rs7132019 chr12 6992122 A G SPSB2 -0.49 0.0006479339 0.37957875 4.09e-71 cg26269324 -0.67 0.00191263489423201 2.26e-45 50 1.21303855863666e-13 0.999999999999886 rs2992756 chr1 18807339 T C KLHDC7A 0.39 0.0009356401 0.47149533 9.11e-35 cg05040210 -0.56 0.00171816693065586 9.17e-37 19 9.57129487693459e-13 0.999999999999034 参考资料#
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/gkx861 ↩︎
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/gky814 ↩︎
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
