🤩 CMplot | 完美复刻Nature上的曼哈顿图（一）
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)
牛奶最后祝大家早日不卷!~
举报
腾讯云开发者
