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

高级的火山图是这样的

火山图

本文涉及到的测试数据及代码已经给大家打包好，如下：

言归正传，火山图是散点图的一种展现形式，一般来说，x轴为实验组基因表达量比上对照组基因表达量的倍数差异（logFC），而y轴则为实验组比对照组之后的p值或者校正后的p值【-log10（adj.p.value）】，图中一个点代表一个基因，而颜色则代表他们是显著上调还是显著下调。

那么，可想而至，绘制火山图，需要三列数据，即logFC、adj.p.value和Symbol基因。这些数据正好是我们差异分析得到的。所以，火山图只是用来可视化那些测序数据差异分析结果而已。

在RNAseq测序中，使用较多的计算差异基因的软件为DESeq2和limma。其差异分析结果文件我们存储在DEGdata.txt文件中，如下：

DEGdata.txt

代码重要部分已被删除，完整代码+数据需要私信要！！！

代码重要部分已被删除，完整代码+数据需要私信要！！！

代码重要部分已被删除，完整代码+数据需要私信要！！！

options(stringsAsFactors = F)
options(warn = -1)

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

# 读取数据
deg.data <- read.table("DEGdata.txt", header = T, sep = "\t")
# 显示文件前6行
head(deg.data)

# 对差异p(adj.P.Val一列进行log10转换
deg.data$logP <- -log10(deg.data$adj.P.Val)

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

# 新加一列Group
# 将adj.P.Val小于0.05，logFC大于2
# 将adj.P.Val小于0.05，logFC小于2的基因设为显著下调基因
# 查看上调和下调基因数目
table(deg.data$Group)

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

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

# 为火山图添加logP分界线(geom_hline)和logFC分界线(geom_vline)


# 新加一列Label
deg.data$Label = ""
# 对差异表达基因的p值进行从小到大排序
deg.data <- deg.data[order(deg.data$adj.P.Val), ]
# 高表达的基因中，选择adj.P.Val最小的10个
up.genes <- head(deg.data$Symbol[which(deg.data$Group == "up-regulated")], 10)
# 低表达的adj.P.Val最小的10个
down.genes <- head(deg.data$Symbol[which(deg.data$Group == "down-regulated")], 10)
# 将up.genes和down.genes合并，并加入到Label中
deg.top10.genes <- c(as.character(up.genes), as.character(down.genes))
deg.data$Label[match(deg.top10.genes, deg.data$Symbol)] <- deg.top10.genes

# 为火山图添加显著差异表达前十名的基因




# 改变火山图点的颜色和坐标轴标注，使图片更美观
ggscatter(deg.data, x = "logFC", y = "logP",
          color = "Group", 
          palette = c("#2f5688", "#BBBBBB", "#CC0000"),
          size = 1,
          label = deg.data$Label, 
          font.label = 8, 
          repel = T,
          xlab = "log2FoldChange", 
          ylab = "-log10(Adjust P-value)") + 
  theme_base() + 
  geom_hline(yintercept = 1.30, linetype="dashed") +
  geom_vline(xintercept = c(-2,2), linetype="dashed")