要玩圖，離不開哈德雷大神的《R數(shù)據(jù)科學(xué)》，第1章和21章是專門講圖的，我寫過對應(yīng)的筆記：
http://m.itdecent.cn/p/4a154f6f0de7
http://m.itdecent.cn/p/bf0f12246865
關(guān)于火山圖加標(biāo)簽的需求，這里有幾種方法來實現(xiàn)。
示例數(shù)據(jù)
方法一的示例數(shù)據(jù)是data.Rdata，方法二三的示例數(shù)據(jù)是test.Rdata。我將數(shù)據(jù)打包放在了“生信星球”公眾號后臺，回復(fù)“火山圖”即可獲得。你解壓后雙擊文件夾里的volcano.Rproj，復(fù)制粘貼運行本文代碼即可。
方法一：利用空字符串“”
原理：空字符串“”=nothing
關(guān)于空字符串，我曾寫過一篇文章來講他：http://m.itdecent.cn/p/aef98f3fc7d8
這種方法的參照是幫助文檔里的一段代碼：
(先準(zhǔn)備好包)
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))
做出的圖是這樣：
可以看到，一部分點有標(biāo)簽， 一部分沒有，思路就是把不要標(biāo)簽的部分變成空字符串“”。
學(xué)以致用
火山圖的本質(zhì)就是點圖，那么在火山圖上標(biāo)記部分基因，就是在點圖上標(biāo)記部分點。
參考這個思路為火山圖加標(biāo)簽：
（美圖預(yù)警）
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:篩選部分基因，用于顯示在圖上
想在圖上做修改，一半是調(diào)參數(shù)，一半是調(diào)數(shù)據(jù)。我們現(xiàn)在要做的就是調(diào)數(shù)據(jù)：要標(biāo)記的，label=基因，無需標(biāo)記的，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)
但是我發(fā)現(xiàn)，這個只是適用于數(shù)據(jù)量比較小的時候，這個例子只有170個點，而一般來說火山圖數(shù)以萬計的行，用這個方法容易失敗。下午嘗試了幾次大的數(shù)據(jù)，結(jié)果Rstudio無一例外的嘎嘣了。
方法二：看R數(shù)據(jù)科學(xué)
代碼來源
以下代碼出自R數(shù)據(jù)科學(xué)筆記第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
)
這個方法適用于較大的數(shù)據(jù)。
端詳代碼找思路
1.從原來數(shù)據(jù)中挑選了一部分，生成新數(shù)據(jù)
2.用新數(shù)據(jù)作圖，向原數(shù)據(jù)做的點圖上疊加兩個圖層，一個空心點圖，一個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：生成用于添加圖層的新數(shù)據(jù)
?重點在這里
新數(shù)據(jù)框的內(nèi)容是你想要標(biāo)記的基因，這里根據(jù)logFC和Pvalue的大小來篩選，可以自定義閾值來調(diào)整要顯示的基因的數(shù)量：
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"
)
加號連接兩句代碼就實現(xiàn)了圖層的疊加，如果對ggplot2不了解，請看R數(shù)據(jù)科學(xué)第一章和第21章。但21章是整本書的錯誤重災(zāi)區(qū)，請看我的筆記有改正后的代碼。
方法三：ggpubr的函數(shù)有現(xiàn)成的參數(shù)
這個函數(shù)叫g(shù)gscatter，還是用剛才的test數(shù)據(jù)來做。
代碼來源
網(wǎng)易云課堂的GEO數(shù)據(jù)庫挖掘?qū)崙?zhàn)演練視頻
由于ggpubr寫縱坐標(biāo)時直接寫-log10（P.value）不識別，可采取迂回策略，改列名，完事再在圖上改縱軸標(biāo)簽。
load("test.Rdata")
if(!require(ggpubr))install.packages("ggpubr")
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")
)
然后加標(biāo)簽，是現(xiàn)成的參數(shù)“l(fā)abel.select”。接受的參數(shù)數(shù)據(jù)結(jié)構(gòu)應(yīng)該是向量。
可以手動選一二三四個感興趣的基因
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]
)
