TCGA数据分析系列之火山图 / 开普饭

前面我们做了TCGA的差异分析，并且用ggplot2验证了差异分析的准确性，TCGA差异分析及ggplot作图验证，而差异分析后一般会又热图和火山图，热图我们之前也有说过热图系列1，R语言学习系列之“多变的热图”，今天我们来了解一下火山图的画法

加载数据

加载差异分析的结果

rm(list = ls())load(file='limma差异分析结果.Rda')library(ggplot2)class(DEG_limma_voom)volcano<-DEG_limma_voom

数据处理

我们先把差异分析结果增加一列，以区分上调下调和无差异的基因。我们以adj.P.Val < 0.05，logFC>2或者logFC<-2作为有明显差异的基因

volcano$type[(volcano$adj.P.Val > 0.05|volcano$adj.P.Val=="NA")|(volcano$logFC < 2)& volcano$logFC > -2] <- "none significant"

volcano$type[volcano$adj.P.Val <= 0.05 & volcano$logFC >= 2] <- "up-regulated"volcano$type[volcano$adj.P.Val <= 0.05 & volcano$logFC <= -2] <- "down-regulated"

纵坐标一般取-log10，绘图看一下

p = ggplot(volcano,aes(logFC,-1*log10(adj.P.Val),color=type)) p + geom_point()

给横坐标设定界限，自定义颜色，加上临界值的虚线

x_lim <- max(volcano$logFC,-volcano$logFC) gg=p + geom_point( size=1) + xlim(-x_lim,x_lim) +labs(x="log2(FC)",y="-log10(adjP)")+ scale_color_manual(values =c("#A52A2A","grey","#f8766d"))+ geom_hline(aes(yintercept=-1*log10(0.05)),colour="black", linetype="dashed") + geom_vline(xintercept=c(-2,2),colour="black", linetype="dashed")print(gg)

看起来还不错

加上基因名

如何把差异最明显的基因名标在火山图上呢？这里提供了几种方法可供选择

方法一

volcano1<-volcanovolcano1$gene_name<-row.names(volcano1)

library(dplyr)top <- volcano1 %>% group_by(type) %>% top_n(n = 10, wt = -1*log10(adj.P.Val)) %>% filter(type !='none significant') # 取上下调中显著差异表达基因前10个

gg + geom_text(data=top,aes(x=logFC, y=-1*log10(adj.P.Val),label=as.character(gene_name)),size=3)

基因名重叠的可以下载PFD格式，在AI里面修改

如果想把背景的灰色去掉们也有办法

library(ggthemes)gg + geom_text(data=top,aes(x=logFC, y=-1*log10(adj.P.Val),label=as.character(gene_name)),size=3)+theme_few()

方法二

library(ggrepel)up <- subset(volcano1, type == 'up-regulated')up <- up[order(up$adj.P.Val), ][1:10, ]down <- subset(volcano1, type == 'down-regulated')down <- down[order(down$adj.P.Val), ][1:10, ]p1 <- gg + theme(legend.position = 'right') + geom_text_repel(data = rbind(up, down), aes(x = logFC, y = -log10(adj.P.Val), label = gene_name), size = 3,box.padding = unit(0.5, 'lines'), segment.color = 'black', show.legend = FALSE) print(p1)

这样基因名就不会重叠在一起，还有箭头指示

方法三

library(ggrepel)up <- subset(volcano1, type == 'up-regulated')up <- up[order(up$adj.P.Val), ][1:10, ]down <- subset(volcano1, type == 'down-regulated')down <- down[order(down$adj.P.Val), ][1:10, ]p2 <- gg + theme(legend.position = 'right') + geom_label_repel(data = rbind(up, down), aes(label = gene_name), show.legend = FALSE)print(p1)

比方法二多了基因的框，看起来也不错。

最后再介绍一个画热图的在线网站：http://www.ehbio.com/ImageGP/index.php/Home/Index/Volcanoplot.html,只要准备好数据输入进去就行了，不过要注意，数据的行名和列名有要求，不能包含特殊字符，比如空格，逗号等，最好就只用字母表示就行。比如我们把数据整理成下面格式