用R画个区域图(LocusZoomlikePlot)

位点区域图(Regional plot)

“

一些简单的刊误与优化：1. 我们使用 Plink 计算位点 LD 关系时，需要指定窗口大小，而在后续的区域图绘制过程中，我们的图形范围也将处于这个窗口内；2. 我更新了代码之后，你可以直接使用下述代码 + 标准的 GWAS(linear) summary（assoc）结果与 ld 计算结果导入并一键输出所有显著位点的区域图，但我还是建议阅读代码并自行根据需要撰写与优化

“

A regional association plot is essentially a zoomed-in Manhattan plot, allowing the researcher to look at associations in a small, pre-defined area of the genome.^[1]

区域图示例

一、 介绍

区域图联合 P 值与变异位点间的 LD 关系，可以帮助我们评价某些显著位点是否是稳定的，可信的。示例图中的区域图是在 LocusZoom^[2] 网站上绘制的。它能够根据研究者上传位点数据和一些公共基因数据库(如 1KG )的 LD 关系情况绘制该图，同时绘制区域位点的密度图和 mapping 这些位点所对应的基因。

美中不足的是，由于有时需要上传的关联数据较大，而且单图模式和其旧网站中无法保存历史上传数据的特性，目前绘制这样的区域图时间成本较高，往往仅绘制较少的位点区域图。其新的网站存在一系列的数据保密条款，此外对于上传数据的格式有着严格要求，而本地运行的 .js 和基于 linux 的命令行模式学习成本较高(Linux 命令行模式已停止支持)。因此目前这个网站对于新手来说还不太适合大量绘制。

区域图在评价显著位点时很直观明了，如果能够快速绘制，对于我们前期评价和筛选位点是很有帮助的。好在当样本量大于 500 时，数据本身的 LD 关系已经足够支持分析。此时就不再需要使用公共基因数据库的LD关系，也就免除了 SNPID(Chr:Pos:Ref:Alt) 转换 rsID 的痛苦。使用 *Plink --r2*^[3] 可以很方便的计算指定位点上下游区域的 LD 关系。同时联合我们的 GWAS summary 数据就能在其他软件中绘制出简易的区域图(无位点密度图与基因 mapping 结果)。

二、 LD 与 GWAS Summary

我们首先从 GWAS summary 数据中提取感兴趣位点的 SNPID list，该 SNPID 需要与基因数据(此处以 Pfile 数据为例)如 .pvar 中的位点ID格式保持一致。然后利用该 list 在基因数据中计算 LD 关系。实现代码如下。

plink --bfile ../LTL001/allchr \
      --r2 \
      --ld-snp-list sig.SNP_ID.txt \ # SNPID list
      --ld-window-kb 1000 \ # 滑动窗口大小
      --ld-window 99999 \ # 窗口内位点数量最大值
      --ld-window-r2 0 \ # 设定最终纳入 LD 结果的位点的 R2 限制
      --out 5e-06 # 结果生成 5e-06.ld 文件

观察 LocusZoom 绘制的区域图可以看出其横轴为染色体上的物理距离，纵轴为位点的关联 P 值，散点颜色映射为 LD r2 值的大小。这对于 ggplot2 来说并不难实现，所以我们需要的数据包括Chr、Pos、ID、P、R2。P 值来自 GWAS summary 数据，R2 来自 plink LD 结果。在 R 中加以整理并绘图就可以实现类似 LocusZoom plot的效果。

三、 R 绘制区域图

R 绘制区域图示例

关于区域图位点评价，我的理解是：1，显著位点因为与区域内周围基因存在一定的 LD 关系，因此与我们选定的显著为位点越接近的位点，其往往也应该是显著或者较为显著的(P 值相接近)；2，由于物理距离增加，位点间的相关性减少，R2 降低。因此，稳定可靠的位点从上到下，应该是存在位点聚集成“柱子的“，这提示高 LD 关系的位点 P 值相近；从中间向两边，其 R2 映射的颜色应该是“渐变”的，如上图所示，从上到下颜色渐变(R2 逐渐降低)，且形成“柱子”，这便是一个稳定可信的位点。

由于需要评价的位点可能较多，因此我使用循环一次性输出所有在 SNPID list 中位点的区域图。需要注意的是，1，绘制时记得对 P 值进行 -log10 转换；2，绘制时按需要可以纳入指定位点与上下游区域内所有位点的 LD 关系。代码如下。

library(data.table)
library(tidyverse)
library(ggsci)

assoc <- fread("../GWAS/LTL001_sep/assocall_control")
ld <- fread("control5e-06.ld")
assoc <- assoc[, c(1, 2, 3, 14)]
colnames(assoc) <- c("Chr", "Pos", "SNPID", "P")
ld <- ld[which(ld$SNP_B %in% assoc$SNPID), ]
dat.sig <- assoc[which(assoc$SNPID %in% ld$SNP_A), ]
dat.sig$R2.group <- as.factor("5")
ld$sigID <- as.factor(ld$SNP_A)

for (i in 1:length(levels(ld$sigID))) {
    dat <- subset(ld[, c(4, 5, 6, 7, 3)], ld$sigID == levels(ld$sigID)[i])
    colnames(dat) <- c("Chr", "Pos", "SNPID", "R2", "sigID")
    datP <- assoc[which(assoc$SNPID %in% dat$SNPID), c(3, 4)]
    dat <- merge(dat, datP, by = "SNPID")
    dat$R2.group <- as.factor(ifelse(dat$R2 < 0.2, "1", ifelse(dat$R2 < 0.4, "2", ifelse(dat$R2 < 0.6, "3", ifelse(dat$R2 < 0.8, "4", "5")))))
    dir <- c()
    dir <- paste("5e-06LZplots/", dat$SNPID[1], "_LZ_plot", ".pdf", sep = "")
    ggplot(dat, aes(x = Pos, y = -log(x = P, base = 10), fill = R2.group)) +
        geom_point(shape = 21, color = "black", size = 4.5, stroke = 0.75, alpha = 0.8) +
        geom_point(data = dat.sig[which(dat.sig$SNPID %in% dat$sigID), ], color = "black", size = 2, alpha = 1, stroke = 0.75, shape = 4) +
        theme_classic() +
        geom_hline(yintercept = 5.30103, color = "#8b0000", lty = 2) +
        geom_hline(yintercept = 6.30103, color = "#8b0000", lty = 2) +
        scale_fill_manual(
            values = c("1" = "#00007B", "2" = "#97CCF6", "3" = "#75FB4C", "4" = "#F2A93B", "5" = "#EA3323")
        )
    ggsave(file = dir)
}

为了画区域图和 LocusZoom 网站死磕一两天，最后还是因为时间原因没办法学会其技术文档，仅能绘制单图。但了解了其绘制原理后，发现如果不需要位点密度图和基因 mapping 情况的话，我们是可以自己绘制的。而且使用较为熟悉的 R 总比使用新工具要来的更容易。以上就是 R 中的简易区域图的绘制方法啦，全文完，感谢浏览。

参考资料

[1]

Regional Plot: https://rdrr.io/cran/QCGWAS/man/plot_regional.html#:~:text=A%20regional%20association%20plot%20is,defined%20area%20of%20the%20genome.

[2]

LocusZoom: http://locuszoom.org/

[3]

Plink1.9 -r: https://www.cog-genomics.org/plink/1.9/ld#r