健明

要玩圖，離不開哈德雷大神的ggplot2，《R數(shù)據(jù)科學(xué)》第1章和21章是專門講圖的，我寫過對(duì)應(yīng)的筆記：

• R數(shù)據(jù)科學(xué)--詳解ggplot2

• R數(shù)據(jù)科學(xué)--第21章 圖形化溝通

• 完整的R數(shù)據(jù)科學(xué)筆記目錄（點(diǎn)擊這里）

關(guān)于火山圖加標(biāo)簽的需求，這里有幾種方法來實(shí)現(xiàn)。

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

方法一的示例數(shù)據(jù)是data.Rdata，方法二三的示例數(shù)據(jù)是test.Rdata。我將數(shù)據(jù)打包放在了“生信技能樹”公眾號(hào)后臺(tái)，回復(fù)“火山圖”即可獲得。你解壓后雙擊文件夾里的volcano.Rproj，復(fù)制粘貼運(yùn)行本文代碼即可。

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

原理：空字符串“”=nothing

關(guān)于空字符串，我曾寫過一篇文章來講他：R數(shù)據(jù)框里的空格子不是NA是什么

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

代碼來源

下面代碼來源于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))

做出的圖是這樣：

學(xué)以致用

火山圖的本質(zhì)就是點(diǎn)圖，那么在火山圖上標(biāo)記部分基因，就是在點(diǎn)圖上標(biāo)記部分點(diǎn)。

參考這個(gè)思路為火山圖加標(biāo)簽：

（美圖預(yù)警）

step1:先把圖畫出來

step2:篩選部分基因，用于顯示在圖上

想在圖上做修改，一半是調(diào)參數(shù)，一半是調(diào)數(shù)據(jù)。我們現(xiàn)在要做的就是調(diào)數(shù)據(jù)：要標(biāo)記的，label=基因，無需標(biāo)記的，label=“”。

?重點(diǎn)就在這里：

data$label=ifelse(data$p.value < 0.000001 & abs(log2(data$FC)) >= 1,data$symbol,'')

step3：將文字圖層疊加上去

但是我發(fā)現(xiàn)，這個(gè)只是適用于數(shù)據(jù)量比較小的時(shí)候，這個(gè)例子只有170個(gè)點(diǎn)，而一般來說火山圖數(shù)以萬計(jì)的行，用這個(gè)方法容易失敗。下午嘗試了幾次大的數(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
  )

端詳代碼找思路

1.從原來數(shù)據(jù)中挑選了一部分，生成新數(shù)據(jù)
2.用新數(shù)據(jù)作圖，向原數(shù)據(jù)做的點(diǎn)圖上疊加兩個(gè)圖層，一個(gè)空心點(diǎn)圖，一個(gè)geom_label_repel。

step1：先把火山圖畫出

step2：生成用于添加圖層的新數(shù)據(jù)

?重點(diǎn)在這里

新數(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:新圖層疊加到原圖上去

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

這個(gè)函數(shù)叫g(shù)gscatter，還是用剛才的test數(shù)據(jù)來做。

代碼來源

當(dāng)然是群主在GitHub的的800M的GEO數(shù)據(jù)挖掘代碼啦，還有配套視頻：

由于ggpubr寫縱坐標(biāo)時(shí)直接寫-log10（P.value）不識(shí)別，可采取迂回策略，改列名，完事再在圖上改縱軸標(biāo)簽。

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

可以手動(dòng)選一二三四個(gè)感興趣的基因

也可以用向量取子集的方法來選很多個(gè)

比如差異基因前30個(gè)

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