「R」数据可视化1: 火山图

2020-07-03 11:58:52 浏览数 (1)

在生物领域我们常常使用R语言对数据可视化。在对数据可视化的时候,我们需要明确想要展示的信息,从而选择最为合适的图突出该信息。本系列文章将介绍多种基于不同R包的作图方法,希望能够帮助到各位读者。

什么是火山图

火山图,是形如火山喷发的一种图形展示方法,常被用于展示差异,比如差异基因、差异微生物等等。

上图就是一张典型的火山图,描述了差异基因的情况。该火山图的y轴是-log10(Qvalue),即qvalue(pvalue校正后的值)取-log10,因此数值越高说明qvalue越小即越显著。横坐标是Log2 fold change,即对fold change取log2,所以越靠两侧的点(每个点代表一个基因),其基因表达量上调或者下调幅度越大。

一般来说在差异基因分析过程中,筛选标准通常认为qvalue小于0.05且foldchange的绝对值大于2为差异基因。

当然在差异基因数量过多的时候,我们可以调整筛选标准以获得相对适合数量的差异基因再进行下游的富集分析,比如要求qvalue小于0.001,或者要求foldchange的绝对值大于5等等。当差异基因数量过少的时候,我们可以考虑将foldchange的绝对值变为1.5,或者考虑选择pvalue小于0.05。

而图中的虚线就是根据自己的筛选标准确定添加。其中两条竖线(x=-2和x=2)说明该筛选标准是要求foldchange的绝对值大于4。横线(大胆猜测是在y=2处),说明要求qvalue小于0.01。当Foldchange没有那么明显的时候,我们的横轴也可以选择展示Foldchange,不取log2。

火山图要怎么画

(1) 需要什么格式的数据

从网上找了一个别人的RNA-seq数据 ,让我们一起来瞧瞧。

数据部分

第一列是基因名,第二列是计算好的log2FoldChange(通常我们用DESeq等软件分析得到的结果里包含这一个值),第三列是pvalue,第四列是校正后的p即padj。如果你已经确定要画pvalue或者padj那只要三列即可。

(2) 如何使用ggpot2做火山图

能够做火山图的方法有很多,有一些RNA-seq分析的包中自带了画火山图的函数。利用R自带的基础画图函数也可以画,但是鉴于之后我们都几乎都选择ggplot2包进行作图,所以只展示如何用ggplot2包画图。

代码语言:javascript复制
#加载包
library(ggplot2)
#读取数据
Dat<-read.table('./results.txt',header = T,stringsAsFactors = F,check.names = F,sep=' ')
#作图
ggplot(Dat,aes(x=log2FoldChange,y=-log10(padj))) 
  geom_point()

通过简单的几行代码我们就能够得到一个最基本(简陋)的火山图。

然后我们来美化这张火山图

代码语言:javascript复制
#加载包
library(ggplot2)
library(ggrepel)
#读取数据
Dat<-read.table('./results.txt',header = T,stringsAsFactors = F,check.names = F,sep=' ')
#确定是上调还是下调,用于给图中点上色)
Dat$threshold = factor(ifelse(Dat$padj < 0.05 & abs(Dat$log2FoldChange) >= 1, ifelse(Dat$log2FoldChange>= 1 ,'Up','Down'),'NoSignifi'),levels=c('Up','Down','NoSignifi'))


ggplot(Dat,aes(x=log2FoldChange,y=-log10(padj),color=threshold)) 
  geom_point() 
  scale_color_manual(values=c("#DC143C","#00008B","#808080")) #确定点的颜色
  geom_text_repel(
    data = Dat[Dat$padj<0.05&abs(Dat$log2FoldChange)>1,],
    aes(label = Gene),
    size = 3,
    segment.color = "black", show.legend = FALSE ) #添加关注的点的基因名
  theme_bw() #修改图片背景
  theme(
    legend.title = element_blank()#不显示图例标题
  ) 
  ylab('-log10 (p-adj)') #修改y轴名称
  xlab('log2 (FoldChange)') #修改x轴名称
  geom_vline(xintercept=c(-1,1),lty=3,col="black",lwd=0.5)  #添加横线|FoldChange|>2
  geom_hline(yintercept = -log10(0.05),lty=3,col="black",lwd=0.5)#添加竖线padj<0.05

Volcano

一张精致的火山图就做好啦。

0 人点赞