R-带指定基因标签的火山图

一、数据准备

二、R代码

# 清除当前环境中的所有对象
rm(list=ls())

# 获取当前工作目录
getwd()

# 设置工作路径
#setwd("") #设置工作路径

# 列出当前工作目录下的所有文件
list.files()

# 加载使用的包
library(openxlsx) # 用于读写Excel文件
library(ggplot2) # 用于数据可视化
library(ggrepel) # 用于在ggplot2中添加带排斥效果的文本标签

# 加载绘制火山图的数据表格
data<-read.xlsx("DEG_SYMBOL.xlsx") # 从Excel文件读取数据

# 将padj列转换为数值类型
data$padj<-as.numeric(data$padj)

# 创建一个新的列'threshold'，根据padj值和log2FoldChange值来判断基因是上调、下调还是无显著变化
data$threshold <- factor(ifelse(data$padj < 0.05 & abs(data$log2FoldChange) >= 1, 
                                ifelse(data$log2FoldChange >= 1, 'up', 'down'), 'no'),
                         levels=c('up', 'down', 'no'))

# 统计每个分类的数量
table(data$threshold)

# 计算上调基因的数量
up<-as.numeric(table(data$threshold)[1])
#up

# 计算下调基因的数量
down<-as.numeric(table(data$threshold)[2])
#down

# 计算无显著变化基因的数量
no.diff<-as.numeric(table(data$threshold)[3])
#no.diff

# 重新计算 'up'、'down' 和 'no' 的计数，这次使用sum函数
up <- sum(data$padj < 0.05 & data$log2FoldChange >= 1)
down <- sum(data$padj < 0.05 & data$log2FoldChange < -1)
no.diff <- sum(!(data$padj < 0.05) | (abs(data$log2FoldChange) < 1))

# 创建因子列 'threshold' 并设置标签，显示每个分类的数量
data$threshold <- factor(ifelse(data$padj < 0.05 & abs(data$log2FoldChange) >= 1,
                                ifelse(data$log2FoldChange >= 1, 'up', 'down'), 'no'),
                         levels = c('up', 'down', 'no'),
                         labels = c(paste("up:", up), paste("down:", down), paste("no:", no.diff)))

# 将带有标签的数据写入新的Excel文件
write.xlsx(data,"data_label.xlsx")

# 开始创建PDF文件，用于保存火山图
pdf("voc_plot.pdf")
#head(data)

# 设置火山图的最低和最高log2FoldChange值(阈值线位置)
left_line <- -1
right_line <- 1

# 查看数据的前几行
head(data)

# 绘制火山图
# 创建火山图对象P
plot <- ggplot(data, aes(x = log2FoldChange, y = -log10(padj), color = threshold)) +
  # 添加点图层，设置点的大小为1，透明度为0.5
  geom_point(size = 1, alpha = 0.5) + 
  # 使用白色背景主题
  theme_bw() + 
  # 自定义图表主题，包括图例标题为空，文本大小为14，坐标轴标题和文本大小为14和12，图例文本大小为12
  theme(legend.title = element_blank(),
        text = element_text(size = 14), 
        axis.title = element_text(size = 14), 
        axis.text = element_text(size = 12), 
        legend.text = element_text(size = 12)) + 
  # 设置y轴标签为'-log10 (p.adj)'
  ylab('-log10 (p.adj)') + 
  # 设置x轴标签为'log2 (FoldChange)'
  xlab('log2 (FoldChange)') + 
  # 添加垂直虚线，表示log2FoldChange的阈值，线型为3（虚线），颜色为黑色，线宽为0.5
  geom_vline(xintercept = c(left_line, right_line), lty = 3, col = "black", lwd = 0.5) + 
  # 添加水平虚线，表示p值调整后的阈值，线型为3（虚线），颜色为黑色，线宽为0.5
  geom_hline(yintercept = -log10(0.05), lty = 3, col = "black", lwd = 0.5) + 
  # 对于数据中show列为"YES"的点，添加排斥效果的文本标签，显示基因符号，文本大小为4，线段颜色为黑色
  geom_text_repel(data = data[data$show == "YES",], 
                  aes(label = Symbol), 
                  size = 4, 
                  segment.color = "black", 
                  color = "black", 
                  min.segment.length = 0, 
                  max.overlaps = 18000) + 
  # 手动设置颜色比例尺，为不同的阈值类别设置不同的颜色
  scale_color_manual(values = c("#0000FF", "#FF0000", "#696969")) + 
  # 设置坐标系，x轴限制为从负无穷到log2FoldChange的最大值
  coord_cartesian(xlim = c(NA, max(data$log2FoldChange)))

print(plot)

# 关闭PDF设备，完成火山图的保存
dev.off()

三、结果