技术讲座

用ggplot绘制火山图无基础教程

来源:小博发布时间:2016-08-10


什么是火山图?

火山图是一类用来展示组间差异数据的图像,因为在生物体发生变化时从全局角度而言大部分的基因表达没有或着发生了很小程度的变化,只有少部分基因的表达发生了显著的变化,我们根据计算结果在图上按照p value和我们自定义的fold change值来将这部分显著变化的基因用不同颜色标示出来以区分,这类图像往往呈现类似火山爆发的样子,于是就被叫做“火山图”(volcano plot)了。
 

 
而在本质上,火山图只是一类特别一点儿的散点图罢了,只不过需要在散点图上划分出几个区域来并对每个区域内的点涂上不同的颜色而已。
什么是ggplot2
ggplot2是R语言的一个包,也是一个绘图系统,它吸收了R语言基础绘图及lattice绘图的优点,并应用图像图形学的理论来进行图像的绘制。它的作者及主要维护者是Hadley Wickham。
ggplot2的官方网站是http://ggplot2.org/
ggplot2的学习资源:
http://docs.ggplot2.org/current/ 这是官方的文档站点,所有信息均是最新的,这是传统渠道所不具备的优势,比如随着版本的更新有些在下面要介绍的两本书里出现的语法已发生了改变,这些只能靠在线文档来获取信息了;
<R Graphics Cookbook>,Winston Chang著,他也是ggplot2的维护者,这本书的中文版名字叫做《R数据可视化手册》,特点是详尽到甚至啰嗦地讲解了ggplot2的每一点用法,适合当做案头参考书。
<ggplot2: Elegant Graphics for Data Analysis>,Hadley Wickham亲著,中文版名字叫做《ggplot2:数据分析与图形艺术》。这本书的评价褒贬不一,但是在相关图书比较缺乏的时代这本书还是值得一读的。
 
准备工作
一台安装了R的电脑((废话);
预先安装了ggplot2包,安装命令是install.packages(“ggplot2”)。
加载包以供使用,命令是library(ggplot2)。注意R包不是安装了就万事大吉了,每次启动R后要使用某个包前必须要加载它才可以用。
准备一份demo数据用以画图,数据内要包含每一个变量的两类参数:p value和fold change;
这里预先准备了一份demo数据,感兴趣的可以发送标题内含有“火山图”三个字的任意邮件到rwzx_zyq@qq.com,会收到自动回复送来的本文相关所有内容。
绘图流程
简单地解释这个过程就是画了一个散点图,再将图分隔为几个区域,每个区域点涂上了不同颜色,最后给图添加文字说明等。


绘图代码

绘图结果
需要说明的是这里的demo数据是一份microRNA的芯片数据,microRNA相对基因来说数目较少,因此差异的microRNA(图上红色和蓝色的点)也就比较少了。大家用转录组表达谱等数据来画会有更密集的点,就好像本文第一幅图里那样。
代码详解
> setwd(“E:/RWorkbench/Volcano”)
设置工作目录,即绘图数据文件所在的目录,这里展示的是windows系统下绝对路径的设定,也可以用相对路径。比如我的默认工作目录就是E:\RWorkbench,则此时只需要用命令setwd(“./Volcano”)即可,注意目录路径的字母大小写不要有误。顺便说一下,查看当前所在目录的命令是getwd()。
> library(ggplot2)
加载ggplot2包准备进行绘图工作,注意R语言里使用一个包除了需要安装它之外每次使用前还需要加载这个包,加载某个包的命令就是library(),不加载就没法用这个包进行工作——哪怕你安装了它。
> demo <- read.table(file=”ESCC_VS_Control.txt”,header = TRUE, row.names =1)
这部分是数据读取部分,将我们存放在ESCC_VS_Control.txt内的文件读入R内,并给读入后的数据取了一个叫做demo的名字以便后面的操作。关于更多数据读入的内容可以查看一些参考资料及函数帮助。
> volcano <- ggplot(demo, aes(x= log2(foldchange), y= -1*log10(pvalues)))
生成了一个名为volcano的对象,这个数据对象内存储了做图的基本信息,这些信息包括:
做图数据来自哪里?——demo
X轴如何定义?——log2(foldchange),即demo数据里名为foldchange的那一列需要取以2为底的对数
Y轴如何定义?——-1*log10(pvalues),即对于demo数据里pvalues那一列的值取以10为底的对数后再乘以-1,目的是便于做图,譬如两个值分别是0.0000012和0.4,若是画图前不做处理则在坐标轴上两者距离较远,区分度也低,而经过-log10()转换后分别变为5.921和0.398,这样做图就比较容易了。
这一部分代码的基本语法是      ggplot(数据对象,aes(x轴定义,y轴定义)。这里的aes是aesthetics的缩写。
> volcano+geom_point(aes(color=significant))+ scale_color_manual(values =c(“red”,”grey”, “blue”))+ labs(title=”Volcanoplot”,x=”log2(fold change)”, y=”-log10(p value)”)+geom_hline(yintercept=1.30103,linetype=3)+geom_vline(xintercept=c(-1,1),linetype=3)+theme(legend.position=’none’)
这一部分是做图的主体了,运行完毕以后就可以预览成品的火山图了。这部分代码看着很长很吓人,但是只要详细解析下就会发现它是由几部分构成,各部分间用“+”连接。
geom是geometric的缩写,表明图像的几何特征,scale是标度,labs是labels。
第一部分:volcano对象,这是用于做图的数据对象,
第二部分:geom_point()
第三部分:scale_color_manual()
第四部分:labs()
第五部分:geom_hline()和geom_vline()
第六部分:theme()
 
下面对这个长代码进行详解。
第一部分是volcano对象这个不用多说,是前一个命令生成的,用于告诉ggplot2用什么东西来做图,两个坐标轴分别如何定义。
第二部分:geom_point(aes(color=significant))。这部分用到了一个geom_point()的命令,这个就是ggplot2用于绘制散点图的命令了,括号内的内容是指定做图系统对不同的变量使用不同的颜色,这里用到的变量分类依据是来自原始数据内的significant列,这一列包括三个内容:up、down和no。分别表示显著上调、下调和不显著表达的基因等。
第三部分:scale_color_manual(values =c(“red”,”grey”, “blue”)),前面一部分声明了要把图上所有点按照significant列的内容不同分为三类,此处scale_color_manual()的作用是给这三类点指定不同的颜色。此处命令是分别指定了红、灰、蓝三色对应up、no、down三类点。注意到了没有?这里三种点的分类是按照字母顺序进行的,颜色匹配也自动对应这一顺序。因此我们可以试着把顺序变换一下来查看效果,比如scale_color_manual(values =c(“grey”,”red”, “blue”))试试看。更多关于顺序的定义有兴趣大家可以去学习一下R语言中因子相关的知识,这里不涉及了。
此外,在颜色指定上除了常见的英文单词表示的颜色外我们还可以用颜色代码来指定,这里提供了一个简便易用的屏幕取色器供大家使用,大家可以用取色器在屏幕任意地方取色,然后把取到的颜色代码填写到scale_color_manual(values =c(“*”,”*”, “*”))三个可被取代的星号那里。
第四部分:labs(title=”Volcanoplot”,x=”log2(fold change)”, y=”-log10(p value)”)。这一部分就比较简单了,分别为图取一个标题以及图上x轴、y轴取一个名称。三部分内容对应于代码内三个引号里的内容。建议使用英文不要用中文。
第五部分:geom_hline(yintercept=1.30103,linetype=3)+geom_vline(xintercept=c(-1,1),linetype=3),其实这部分在ggplot2里合起来叫做geom_abline(),此处它的作用是在图上绘制出几条虚线来将图上坐标点分隔为几个区域。我们按照火山图的定义需要画3条线:平行于x轴的p value=0.05的那一条线以及垂直于x坐标轴的分别代表fold change>2及<0.5(也就是1/2)的那两条线。
我们知道,当p value=0.05时它的-log10值是1.30103,因此在geom_hline(yintercept=1.30103,linetype=3)这部分的参数就是如此设定,linetype=3代表这是一个细细的虚线条。
当fold change>2或<0.5时log2()值分别为1和-1,因此在geom_vline(xintercept=c(-1,1),linetype=3)部分我们设置参数为-1和1.
实际上,R也有简单的计算器功能,我们在修改参数是可以在R内直接计算这些值,比如
> -1*log10(0.05)[1] 1.30103
直接计算得到p value=0.05时-log10p()的值。
长代码的最后一部分是theme(legend.position=’none’),这部分的作用很简单:去掉图例。大家可以把这部分删除后看看运行效果就知道了。
 
至此,火山图就已经画完了,我们可以用R语言运行环境内的命令保存图片,也可以用ggplot2自带的保存命令保存图片,比如:
>ggsave(“ESCC_vs_Control_volcano.png”)
>Saving 7.62 x 5.84 in image
这样就把图片保存到本地了。