RNA-Seq分析中（一）——用R画带基因名标签的火山图

火山图

火山图（volcano plot）常用于显著差异基因表达的展示，包含显著和差异两个重要信息。显著性指P值小于0.05，差异性常用FoldChange值>=2作为筛选标准。

一、火山图简介

那么如何看懂一张火山图所包含的信息呢？首先需要知道，火山图的横坐标通常用log2（fold change）表示，差异越大的基因分布在两端，纵坐标用-log10（pvalue）表示，T检验显著性P值的负对数。由于P值越小表示越显著，所以我们进行-log10（P value）转化后，转化值越大表示差异越显著。通常差异倍数越大的基因T检验越显著，所以左上角和右上角的值往往是我们关注的。

二、火山图绘制

1. 数据准备（IL-1B_1_IL-1B_2_vs_Cu5_1_Cu5_2.all.xls）

第一列为基因ID，倒数第二列列为FDR，最后一列列为log2FC。

2. 我们使用R语言来绘制火山图，代码如下：

# 安装ggpubr和ggthemes包
#install.packages("ggpubr")
#install.packages("ggthemes")
# 加载ggpubr和ggthemes包
library(ggpubr)
library(ggthemes)

# 读取数据
deg.data <- read.table("IL-1B_1_IL-1B_2_vs_Cu5_1_Cu5_2.all.xls",header=TRUE,check.names =F, sep = "\t")
# 显示文件
head(deg.data)

# 对差异FDR进行log10转换
deg.data$logQ <- -log10(deg.data$FDR)

# 绘制基本热图
ggscatter(deg.data, x = "log2FC", y = "logQ") + theme_base()

# 新加一列Group
deg.data$Group = "normal"
# 将adj.P.Val小于0.05，logFC大于2
# 将adj.P.Val小于0.05，logFC小于2的基因设为显著下调基因
deg.data$Group[which( (deg.data$FDR < 0.01) & (deg.data$log2FC > 1) )] = "up"
deg.data$Group[which( (deg.data$FDR < 0.01) & (deg.data$log2FC < -1) )] = "down"
# 查看上调和下调基因数
table(deg.data$Group)

# 绘制新的火山图
ggscatter(deg.data, x = "log2FC", y = "logQ",
          color = "Group") + theme_base()

# 改变火山图颜色(palette)和点的大小(size)
ggscatter(deg.data, x = "log2FC", y = "logQ",
          color = "Group", 
          palette = c("green", "gray", "red"),
          size = 1) + theme_base()

# 为火山图添加logP分界线条(geom_hline)和logFC分界线条(geom_vline)
ggscatter(deg.data, x = "log2FC", y = "logQ",
          color = "Group", 
          palette = c("green", "gray", "red"),
          size = 1) + theme_base() + 
  geom_hline(yintercept = 1.30, linetype="dashed") +
  geom_vline(xintercept = c(-2,2), linetype="dashed")


# 新加一列Label
deg.data$Label = ""
# 对差异表达基因的log2FC值进行从大到小排序取top 20
deg.data <- deg.data[order(abs(deg.data$log2FC),decreasing = T), ]
log2FC.genes <- head(deg.data$ID, 20)
#head(deg.data)
# 对差异表达基因的p值进行从小到大排top 20
# deg.data <- deg.data[order(deg.data$logQ)]
# 高表达的基因中，选择FDR最小的20
deg.data <- deg.data[order(abs(deg.data$logQ),decreasing = T), ]
fdr.genes <- head(deg.data$ID, 20)

# 将log2FC.genes 和fdr.genes合并，并加入到Label
deg.top20.genes <- c(as.character(log2FC.genes), as.character(fdr.genes))
deg.data$Label[match(deg.top20.genes, deg.data$ID)] <- deg.top20.genes
#print (deg.data$Label)



# 改变火山图点的颜色和坐标轴标注，使图片更美观
ggscatter(deg.data, x = "log2FC", y = "logQ",
          color = "Group", 
          palette = c("#00BA38", "#BBBBBB", "#F8766D"),
          size = 2,
          label = deg.data$Label, 
          font.label = 8, 
          repel = T,
          xlab = "log2FC", 
          ylab = "-log10(FDR)") + 
  theme_base() + 
  geom_hline(yintercept = 1.30, linetype="dashed") +
  geom_vline(xintercept = c(-2,2), linetype="dashed")

三、火山图结果

火山图.png