多种方法在火山图上标记感兴趣基因（差异基因，或者通路）
健明
全国巡讲课程结束后的一个月持续答疑环节，被问的最多的问题居然是如何在差异分析后的火山图上面标记出来感兴趣的基因，这里有必要派我们杰出能干的小洁老师出马！
要玩图，离不开哈德雷大神的ggplot2，《R数据科学》第1章和21章是专门讲图的，我写过对应的笔记：
关于火山图加标签的需求，这里有几种方法来实现。
示例数据
方法一的示例数据是data.Rdata，方法二三的示例数据是test.Rdata。我将数据打包放在了“生信技能树”公众号后台，回复“火山图”即可获得。你解压后双击文件夹里的volcano.Rproj，复制粘贴运行本文代码即可。
方法一：利用空字符串“”
原理：空字符串“”=nothing
关于空字符串，我曾写过一篇文章来讲他：R数据框里的空格子不是NA是什么
这种方法的参照是帮助文档里的一段代码： (先准备好包)
if(!require(ggplot2)) install.packages("ggplot2")
if(!require(ggrepel)) install.packages("ggrepel")
if(!require(dplyr))install.packages("dplyr")
library(ggplot2)
library(ggrepel)
library(dplyr)
代码来源
下面代码来源于geom_text_repel的帮助文档
p <- ggplot(mtcars,
aes(wt, mpg, label = rownames(mtcars), colour = factor(cyl))) +
geom_point()
# Hide some of the labels, but repel from all data points
mtcars$label <- rownames(mtcars)
mtcars$label[1:15] <- ""
p + geom_text_repel(data = mtcars, aes(wt, mpg, label = label))
做出的图是这样：
可以看到，一部分点有标签， 一部分没有，思路就是把不要标签的部分变成空字符串“”。
学以致用
火山图的本质就是点图，那么在火山图上标记部分基因，就是在点图上标记部分点。
参考这个思路为火山图加标签：
（美图预警）
step1:先把图画出来
load("data.Rdata")
head(data)
# symbol p.value FC change
#1 PCMTD2 1.53544e-11 1.3548360 Stable
#2 KIAA0087 6.71382e-13 0.7314603 Stable
#3 AFAP1L1 4.24611e-12 0.6284560 Stable
#4 CHMP1A 3.76821e-09 1.6035994 Stable
#5 TRERF1 1.80652e-08 0.6875469 Stable
#6 C8B 7.88047e-04 1.2374303 Stable
data$change = ifelse(data$p.value < 0.000001 & abs(log2(data$FC)) >= 1,
ifelse(log2(data$FC)> 1 ,'Up','Down'),
'Stable')
p <- ggplot(data = data,
aes(x = log2(data$FC),
y = -log10(data$p.value),
colour=change,
label = data$symbol)) +
geom_point(alpha=0.4, size=3.5) +
scale_color_manual(values=c("blue", "grey","red"))+
xlim(c(-4.5, 4.5)) +
geom_vline(xintercept=c(-1,1),lty=4,col="black",lwd=0.8) +
geom_hline(yintercept = -log10(0.000001),lty=4,col="black",lwd=0.8) +
labs(x="log2(fold change)",
y="-log10 (p-value)",
title="Differential metabolites") +
theme_bw()+
theme(plot.title = element_text(hjust = 0.5),
legend.position="right",
legend.title = element_blank())
p
step2:筛选部分基因，用于显示在图上
想在图上做修改，一半是调参数，一半是调数据。我们现在要做的就是调数据：要标记的，label=基因，无需标记的，label=“”。
⭐重点就在这里：
data$label=ifelse(data$p.value < 0.000001 & abs(log2(data$FC)) >= 1,data$symbol,"")
step3：将文字图层叠加上去
p+geom_text_repel(data = data, aes(x = log2(data$FC),
y = -log10(data$p.value),
label = label),
size = 3,box.padding = unit(0.5, "lines"),
point.padding = unit(0.8, "lines"),
segment.color = "black",
show.legend = FALSE)
但是我发现，这个只是适用于数据量比较小的时候，这个例子只有170个点，而一般来说火山图数以万计的行，用这个方法容易失败。下午尝试了几次大的数据，结果Rstudio无一例外的嘎嘣了。
方法二：看R数据科学
代码来源
以下代码出自R数据科学笔记第21章，原书第312页：
best_in_class <- mpg %>%
group_by(class) %>%
filter(row_number(desc(hwy)) == 1)
ggplot(mpg, aes(displ, hwy)) +
geom_point(aes(color = class)) +
geom_point(size = 3, shape = 1, data = best_in_class) +
ggrepel::geom_label_repel(
aes(label = model),
data = best_in_class
)
这个方法适用于较大的数据。
端详代码找思路
1.从原来数据中挑选了一部分，生成新数据 2.用新数据作图，向原数据做的点图上叠加两个图层，一个空心点图，一个geom_label_repel。
step1：先把火山图画出
load("test.Rdata")
p <- ggplot(data = test,
aes(x = logFC,
y = `-log10(P.value)`)) +
geom_point(alpha=0.4, size=3.5,
aes(color=change)) +
scale_color_manual(values=c("blue", "grey","red"))+
geom_vline(xintercept=c(-1,1),lty=4,col="black",lwd=0.8) +
geom_hline(yintercept = -log10(0.01),lty=4,col="black",lwd=0.8) +
theme_bw()
p
step2：生成用于添加图层的新数据
⭐重点在这里
新数据框的内容是你想要标记的基因，这里根据logFC和Pvalue的大小来筛选，可以自定义阈值来调整要显示的基因的数量：
for_label <- test %>%
filter(abs(logFC) >4& `-log10(P.value)`> -log10(0.000001))
step3:新图层叠加到原图上去
p +
geom_point(size = 3, shape = 1, data = for_label) +
ggrepel::geom_label_repel(
aes(label = symbol),
data = for_label,
color="black"
)
加号连接两句代码就实现了图层的叠加，如果对ggplot2不了解，请看R数据科学第1章和第21章。但21章是整本书的错误重灾区，请看我的笔记有改正后的代码。
方法三：ggpubr的函数有现成的参数
这个函数叫ggscatter，还是用刚才的test数据来做。
代码来源
当然是群主在GitHub的的800M的GEO数据挖掘代码啦，还有配套视频：
由于ggpubr写纵坐标时直接写-log10（P.value）不识别，可采取迂回策略，改列名，完事再在图上改纵轴标签。
load("test.Rdata")
if(!require(ggpubr))install.packages("ggplubr")
library(ggpubr)
colnames(test)[4] <- "v"
ggscatter(test,
x = "logFC",
y ="v",
ylab="-log10(P.value)",
size=0.5,
color = "change",
palette = c("#00AFBB", "#999999", "#FC4E07")
)
然后加标签，是现成的参数“label.select”。接受的参数数据结构应该是向量。
可以手动选一二三四个感兴趣的基因
ggscatter(test,
x = "logFC",
y = "v",
ylab="-log10(P.value)",
color = "change",
size = 0.5,
label = "symbol",
repel = T,
palette = c("#00AFBB", "#999999", "#FC4E07") ,
#label.select = dat$symbol[1:30] ,
label.select = c("CD36", "DUSP6", "DCT", "SPRY2", "MOXD1", "ETV4" )
)
也可以用向量取子集的方法来选很多个
比如差异基因前30个
ggscatter(test,
x = "logFC",
y = "v",
ylab="-log10(P.value)",
color = "change",
size = 0.5,
label = "symbol",
repel = T,
palette = c("#00AFBB", "#999999", "#FC4E07") ,
label.select = test$symbol[1:30]
)
