jamesbang

V1

2022/09/02阅读:56主题:雁栖湖

🤣 CMplot | 完美复刻Nature上的曼哈顿图(一)

11. 写在前面

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

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

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

22. 用到的包

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

33. 示例数据

data(pig60K)
data(cattle50K)


Note! 示例数据中的前三列分别是SNP的名称、染色体和位置,
其余各列是GWASp值或traitsGS/GP

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

44. 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 **

55. 修改细节

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)

66. 更进一步

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

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)

牛奶
牛奶
最后祝大家早日不卷!~

点个在看吧各位~ ✐.ɴɪᴄᴇ ᴅᴀʏ 〰

分类:

后端

标签:

后端

作者介绍

jamesbang
V1