在火山图上标记基因(方法一)

   大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~

   就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~

   这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！

要玩图，离不开哈德雷大神的《R数据科学》，第1章和21章是专门讲图的，我写过对应的笔记：来自小洁的《R数据科学》目录

关于火山图加标签的需求，有几种方法来实现，今天介绍方法一。

示例数据

方法一的示例数据是data.Rdata，方法二三的示例数据是test.Rdata。我将数据打包放在了“生信星球”公众号后台，回复“火山图”即可获得。你解压后双击文件夹里的volcano.Rproj，复制粘贴运行本文代码即可。

原理：空字符串“”=nothing

关于空字符串，我曾写过一篇文章来讲他：R数据框里的空格子不是NA是什么

这种方法的参照是帮助文档里的一段代码：
(先准备好包)

if(!require(ggplot2)) install.packages('ggplot2')
if(!require(ggrepel)) install.packages('ggrepel')
library(ggplot2)
library(ggrepel)


代码来源

下面代码来源于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,'')

如果看不懂ifelse的语法，请看这一篇：从前提到条件语句只知道if else

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无一例外的嘎嘣了。

另外两种方法我将在后续推文中继续更新。