曼哈顿图的本质是一种散点图,通过其独特的视觉呈现方式,能够清晰地展示大量数据点中的显著特征,其在全基因组关联研究(GWAS)中发挥着重要作用。曼哈顿图通过直观展示每个SNP的关联强度(以P值的负对数为指标),帮助研究人员快速定位到潜在的关联位点。
R语言作为专门为统计分析和数据处理开发的编程语言,功能化十分强大,有很多的专业的包可以下载并使用。其中CMplot包就是专门为绘制曼哈顿图开发的包,下载并加载CMplot包的代码如下:
install.packages('CMplot')
library(CMplot)
为了方便处理数据我们还需要加载以下的包,并设置绘图时,坐标轴数字不使用科学计数法:
library(openxlsx)
library(tidyverse)
# 设置当数字的绝对值大于10^200的时候使用科学计数法
options(scipen=200)
为了方便说明CMplot包中的CMplot函数的使用方法,我们首先制作一个假数据,代码如下,结果如图1:
# 制作假数据,输出,注意,以下的程序必须满足满足假数据的形式,否则不发成功运行
# 制作两个假数据
set.seed(123)
Pos = sample(10000:50000,1000,replace = FALSE)
data1 = data.frame(SNP = paste0('gen1_',Pos),
Chromosome = sample(1:5,1000,replace = TRUE),
Position = Pos,
P.value = runif(1000,min = 0.431231,max = 0.99999))
set.seed(321)
Pos = sample(10000:50000,1000,replace = FALSE)
data2 = data.frame(SNP = paste0('gen1_',Pos),
Chromosome = sample(1:5,1000,replace = TRUE),
Position = Pos,
P.value = runif(1000,min = 0.501,max = 0.99999))
# 一定一定把数据做成这个样子再开始下面的程序
write.xlsx(data1,'fake_data1.xlsx')
write.xlsx(data2,'fake_data2.xlsx')
fake_data1.xlsx的一部分数据
fake_data2.xlsx的一部分数据
现在假设我们有两组数据,两组数据的变量从左到右都是按照SNP标记、染色体、位置和P值排列的。当我们准备好类似上面的数据后可以参考下面的代码将数据合并起来,并计算相对位置,代码如下:
data1.1 = data1%>%
mutate(var = 'gene1')
data2.1 = data2%>%
mutate(var = 'gene2')
data3 = rbind(data1.1,data2.1)
colnames(data3)[2] = 'Chr'
chr_len = data3%>%
group_by(Chr)%>%
summarise(Chr_len = max(Position))
chr_pos = chr_len%>%
mutate(total = cumsum(Chr_len)-Chr_len)%>%
select(-Chr_len)
data4 = data3%>%
left_join(chr_pos,by ='Chr')%>%
arrange(Chr,Position)%>%
mutate(Position_cum = Position+total)
data5 = data4%>%
select(SNP,Chr,Position,P.value,var)%>%
pivot_wider(id_cols = c('SNP','Chr','Position'),names_from = 'var',values_from = 'P.value')
write.xlsx(data5,'fake_data5.xlsx')
最终的data5就可以被用来画曼哈顿图了,data5的一部分示例如下:
data5的一部分
下面我们介绍一些CMplot函数里面常用的参数:
col 设置颜色
cex/pch 设置点的大小/形状
bin.size 设置SNP密度图中的窗口大小
cex.axis 设置坐标轴字体和标签字体的大小
plot.type 设置不同的绘图类型,可以设定为 "d", "c", "m", "q" or "b",其中d是SNP密度图,c是环形曼哈顿图,m是曼哈顿图,q是qq图,b是同时画环形曼哈顿图、曼哈顿图和qq图
LOG10 = TRUE 数值轴进行LOG10转换
threshold:阈值并添加阈值线 / threshold.col:阈值线的颜色 / threshold.lwd:阈值线的宽度/ threshold.lty:阈值线的类型
signal.cex:显著点的大小 /signal.pch:显著点的形状 /signal.col:显著点的颜色
cir.legend:设置是否显示图例 /cir.legend.cex:图例字体大小 /cir.legend.col :图例颜色
使用我们的假数据data1、data、data5来绘制曼哈顿图,有如下的几种情况:
1. 如果只想绘制一个数据对应的曼哈顿图,例如使用data1来画,代码如下,运行下面的代码,程序会自动把图片输出到rmd文件所在的文件夹里面,效果如图:
CMplot(data1,plot.type="m",LOG10=TRUE,
amplify=TRUE,chr.den.col=c("darkgreen","yellow","red"),bin.size=1e6,
signal.col=c("red","green"),signal.cex=c(1,1),signal.pch=c(19,19),file.name = 'gene1',file = 'jpg')
Rect_Manhtn.gene1.jpg
2. 如果多个数据画在一个图里面,但是分多层,例如使用data5来画,代码如下,效果如图:
CMplot(data5,plot.type="m",LOG10=TRUE,multracks = TRUE,
amplify=TRUE,chr.den.col=c("darkgreen","yellow","red"),bin.size=1e6,
signal.col=c("red","green"),signal.cex=c(1,1),signal.pch=c(19,19),file.name = 'all',file = 'jpg')
Multi-tracks_Manhtn.all.jpg
3. 如果多个数据画在一个图里面,点混在一起,代码如下,效果如图:
Multi-traits_Manhtn.all-in-one.jpg
4. 如果想绘制QQ图,可以设置plot.type = 'q',代码如下,效果如图:
CMplot(data5,plot.type = 'q',
col = c('yellow','purple'),
amplify = TRUE,
multraits = TRUE,
file.output = TRUE,
file.name = 'gene1&2',file = 'jpg')
Multi-traits_QQplot.gene1&2.jpg
5. 环形曼哈顿图,设置plot.type = 'c',代码如下,效果如图:
CMplot(data1,plot.type='c',file.name = 'gene1',file = 'jpg')
Cir_Manhtn.gene1.jpg
效果与正常的曼哈顿图是差不多的,只是形状变了。
6. 最后,如果需要SNP密度图,设置plot.type = 'd',代码如下,效果如图:
CMplot(data5,plot.type = 'd',
col = c('yellow','purple'),
amplify = TRUE,
multraits = TRUE,
file.output = TRUE,
file.name = 'gene1&2',file = 'jpg')
Marker_Density.gene1&2.jpg
最后,本文中的数据都是程序生成的假数据,只作为示例使用,没有任何实际意义。
本文所有的代码、数据集、图片等文件可以通过下面的链接获得:
链接:https://pan.baidu.com/s/1dqrqB2dXbQD0-ASrpUUogA?pwd=73t1
提取码:73t1
转发
0
评论
0