🤩 CMplot | 完美复刻Nature上的曼哈顿图（一）

生信漫卷

发布于 2022-10-31 17:05:27

1.5K0

发布于 2022-10-31 17:05:27

文章被收录于专栏：R语言及实用科研软件R语言及实用科研软件

1. 写在前面

原文地址：https://www.nature.com/articles/s41562-022-01359-x

今天要复刻一下Nature human behaviour上的一张Manhattan plot

曼哈顿图的名字来源是因为其形如曼哈顿的天际线, 高耸在较低高度的“建筑物”上方的摩天大楼的轮廓。主要用于GWAS结果的展示。

2. 用到的包

rm(list = ls())
library(CMplot)

3. 示例数据

data(pig60K)
data(cattle50K)

Note! 示例数据中的前三列分别是SNP的名称、染色体和位置，其余各列是GWAS的p值或traits的GS/GP。

有时候你可能只是想画SNP的密度图，这个时候前3列就够了。

4. GWAS结果展示

CMplot(pig60K,
       type="p",
       plot.type="m",
       LOG10=TRUE,
       threshold=NULL,
       file="jpg",
       memo="",
       dpi=300,
       file.output=F,
       verbose=F,
       width=14,height=6,
       chr.labels.angle=0)

这里需要说明一下plot.type可选的有：

✅ "d" → SNP density plot ✅ "c" → circle-Manhattan plot ✅ "m" → Manhattan plot ✅ "q" → Q-Q plot ✅ "b" → both circle-Manhattan, Manhattan and Q-Q plots

5. 修改细节

5.1 更改颜色

CMplot(pig60K,
       col = c("#3E0A52", "#423D77","#3F678B",
               "#468C8D", "#5FB47F", "#9FD55C","#F9E956"),
       type="p",
       plot.type="m",
       LOG10=TRUE,
       threshold=NULL,
       file="jpg",
       memo="",
       dpi=300,
       file.output=F,
       verbose=F,
       width=14,height=6,
       chr.labels.angle=0)

5.2 标注基因

我们在这里假设存在几个基因名gene1,gene2,gene3, gene4.... balabala......

SNPs <- pig60K[pig60K[,5] < (0.05 / nrow(pig60K)), 1]
genes <- paste("GENE", 1:length(SNPs), sep="_")

CMplot(pig60K[,c(1:3,5)],
       plot.type="m",
       LOG10=TRUE,
       col= c("#3E0A52", "#423D77","#3F678B",
               "#468C8D", "#5FB47F", "#9FD55C","#F9E956"),
       highlight = SNPs,
       highlight.col = NULL,
       highlight.cex = 1,
       highlight.pch = c(15:17), 
       highlight.text = genes,      
       highlight.text.col = "black",
       threshold = 0.05/nrow(pig60K),
       amplify = FALSE,
       file = "jpg",
       memo = "",
       dpi = 300,
       file.output = F,
       verbose = F,
       width = 14,height = 6)

5.3 更改阈值线的颜色和类型

CMplot(pig60K[,c(1:3,5)],
       plot.type="m",
       LOG10=TRUE,
       col= c("#3E0A52", "#423D77","#3F678B",
               "#468C8D", "#5FB47F", "#9FD55C","#F9E956"),
       highlight = SNPs,
       highlight.col = NULL,
       highlight.cex = 1,
       highlight.pch = c(15:17), 
       highlight.text = genes,      
       highlight.text.col = "black",
       threshold = 0.05/nrow(pig60K),
       threshold.lty = 2, 
       threshold.col = "black",
       threshold.lwd = 3,
       amplify = FALSE,
       file = "jpg",
       memo = "",
       dpi = 300,
       file.output = F,
       verbose = F,
       width = 14,height = 6)

6. 更进一步

这里我们再加上染色体的图，锦上添花！~🤩

CMplot(pig60K[,c(1:3,5)],
       plot.type="m",
       LOG10=TRUE,
       col= c("#3E0A52", "#423D77","#3F678B",
               "#468C8D", "#5FB47F", "#9FD55C","#F9E956"),
       chr.den.col=c("darkgreen", "yellow", "red"),
       highlight = SNPs,
       highlight.col = NULL,
       highlight.cex = 1,
       highlight.pch = c(15:17), 
       highlight.text = genes,      
       highlight.text.col = "black",
       threshold = 0.05/nrow(pig60K),
       threshold.lty = 2, 
       threshold.col = "black",
       threshold.lwd = 3,
       amplify = FALSE,
       file = "jpg",
       memo = "",
       dpi = 300,
       file.output = F,
       verbose = F,
       width = 14,height = 6)