在火山圖上標記基因

要玩圖勋篓，離不開哈德雷大神的《R數(shù)據(jù)科學》梢杭，第1章和21章是專門講圖的宜狐，我寫過對應的筆記：
http://www.reibang.com/p/4a154f6f0de7
http://www.reibang.com/p/bf0f12246865

關于火山圖加標簽的需求盒延，這里有幾種方法來實現(xiàn)节槐。

示例數(shù)據(jù)

方法一的示例數(shù)據(jù)是data.Rdata酒唉，方法二三的示例數(shù)據(jù)是test.Rdata旱眯。我將數(shù)據(jù)打包放在了“生信星球”公眾號后臺喂链，回復“火山圖”即可獲得铐姚。你解壓后雙擊文件夾里的volcano.Rproj策肝，復制粘貼運行本文代碼即可。

方法一：利用空字符串“”

原理：空字符串“”=nothing

關于空字符串谦屑，我曾寫過一篇文章來講他：http://www.reibang.com/p/aef98f3fc7d8

這種方法的參照是幫助文檔里的一段代碼：
(先準備好包)

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:篩選部分基因，用于顯示在圖上

想在圖上做修改坦康，一半是調參數(shù)竣付，一半是調數(shù)據(jù)。我們現(xiàn)在要做的就是調數(shù)據(jù)：要標記的滞欠，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)

但是我發(fā)現(xiàn)逸绎，這個只是適用于數(shù)據(jù)量比較小的時候惹恃，這個例子只有170個點，而一般來說火山圖數(shù)以萬計的行棺牧，用這個方法容易失敗巫糙。下午嘗試了幾次大的數(shù)據(jù)，結果Rstudio無一例外的嘎嘣了颊乘。

方法二：看R數(shù)據(jù)科學

代碼來源

以下代碼出自R數(shù)據(jù)科學筆記第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ù)框的內容是你想要標記的基因，這里根據(jù)logFC和Pvalue的大小來篩選识啦，可以自定義閾值來調整要顯示的基因的數(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ù)科學第一章和第21章颓哮。但21章是整本書的錯誤重災區(qū)家妆，請看我的筆記有改正后的代碼。

方法三：ggpubr的函數(shù)有現(xiàn)成的參數(shù)

這個函數(shù)叫ggscatter冕茅，還是用剛才的test數(shù)據(jù)來做伤极。

代碼來源

網易云課堂的GEO數(shù)據(jù)庫挖掘實戰(zhàn)演練視頻

https://study.163.com/course/introduction.htm?courseId=1005840004&share=2&shareId=400000000562025#/courseDetail?tab=1

由于ggpubr寫縱坐標時直接寫-log10（P.value）不識別，可采取迂回策略姨伤，改列名哨坪，完事再在圖上改縱軸標簽。

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") 
          )

然后加標簽乍楚，是現(xiàn)成的參數(shù)“l(fā)abel.select”当编。接受的參數(shù)數(shù)據(jù)結構應該是向量。

可以手動選一二三四個感興趣的基因

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]
          )