volcano_plot-详细注释

东北林业大学_LXG_曾可豪

已于 2023-07-06 20:44:24 修改

阅读量159

点赞数

文章标签： r语言

于 2023-07-06 16:49:12 首次发布

本文链接：https://blog.csdn.net/buchongmingle/article/details/131524190

版权

该文介绍了一种使用R语言和DESeq2包处理基因表达数据的方法，包括准备数据、设定阈值、基础绘图以及为关键基因添加标签。通过ggplot2库创建火山图，用于展示基因的log2FoldChange和padj值，同时突出显示感兴趣的基因。

摘要由CSDN通过智能技术生成

一、如果你有改进的意见，方便的话希望能评论告诉我，我十分想要将我的代码改进到傻瓜式操作。如果你也在学习R和linux中，欢迎私信，共同学习进步。

二、准备数据

需要准备的文件有：1、DESeq2的结果文件“DEseq2_allresults.csv”。2、感兴趣或通路中关键的基因列表“labelgene_list.txt”。

其中，DEseq2_allresults.csv文件行名为基因id，必须具备的列名是：X（基因id的列名，可以改为“geneid”）、log2FoldChange、padj和change（up或down），如何获得这个文件请参考我的关于DESeq2的另一篇文章。

“labelgene_list.txt”文件内容十分简单，就是基因ID（当然基因名也可以）。

三、整理基本的绘图数据格式

要做的就是提取需要的数据列，使用which会方便一点（对我来说）。

同时也需要设置“显著”和“差异”的阈值。

data1 <- read.table("DEseq2_allresults.csv", header = TRUE, sep = ",")
Target_IndexValue <- c(1,
                       which((colnames(data1) %in% "log2FoldChange") == 1),
                       which((colnames(data1) %in% "padj") == 1),
                       which((colnames(data1) %in% "change") == 1))#检索目的列的索引
data2 <- data1[,Target_IndexValue,drop = FALSE]#取出目的列
#data3 <- data2[complete.cases(data2), ]#删除具有缺失值NA的行（DESeq2代码中已删除）

#设置padj和log2FC的阈值
#0.05、0.01、0.001
cut_off_padj1= 0.05
cut_off_padj2= 0.01
cut_off_padj3= 0.001
cut_off_padj4= 0.0001#统计显著性
cut_off_log2FC = 1#差异倍数值

四、基本的绘图

其实注释很清楚了。

library(ggplot2)
p <- ggplot(
  #数据来源
  #xy值和颜色的映射关系，颜色取决于change列
  data2, aes(x = log2FoldChange, y = -log10(padj), colour=change)) +
  geom_point(alpha=0.4, size=2) +#点的透明度和大小
  scale_color_manual(values=c("#546de5", "#d2dae2","#ff4757"))+#点的颜色
  #边界线，数字0-6代表线的类型
  geom_vline(xintercept=c(-1,1),lty=4,col="black",lwd=0.8) +
  geom_hline(yintercept = -log10(cut_off_padj1),lty=4,col="black",lwd=0.8) +
  #坐标轴标签
  labs(x="log2(FoldChange)",
       y="-log10(p.adj)")+
  #设置图片背景
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5), #标题居中，hjust = 0.5表示水平居中对齐
        legend.position="right", #图例相对位置
        legend.title = element_blank()) +#图例无标题
  #调整坐标轴标题和坐标轴刻度标签的字号和字体样式
  theme(axis.title.x =element_text(size=14,face = "bold"), 
      axis.title.y=element_text(size=14,face = "bold"), 
      axis.text = element_text(size = 14,face = "bold"))

p

五、为基因添加标签

可以为显著且差异倍数大的基因加标签，但更推荐为关键基因添加标签。

library(ggrepel)
target_gene_list <- as.character(unlist(read.table("labelgene_list.txt")))
#显著大差异的基因绘制标签
#data2$label <- ifelse(data2$padj < cut_off_padj4 & abs(data2$log2FoldChange) >= 10, data2$X, "")
#感兴趣的基因绘制标签
data2$label <- ifelse(data2$X %in% target_gene_list, data2$X, "")#注意！"%in%"需要被比较者为字符或数值向量

p + geom_label_repel(
    data = data2[data2$X %in% target_gene_list,],
    aes(label = label),
    size = 3,
    segment.color = "black",
    show.legend = FALSE,
    box.padding = unit(0.5, "lines"),
    point.padding = unit(0.8, "lines"),
    max.overlaps = 20,
    )

六、整体的代码

#重来
rm(list = ls())
setwd("D:/great_work/R_wd/volcano_plot/")

#整理数据
data1 <- read.table("DEseq2_allresults.csv", header = TRUE, sep = ",")
Target_IndexValue <- c(1,
                       which((colnames(data1) %in% "log2FoldChange") == 1),
                       which((colnames(data1) %in% "padj") == 1),
                       which((colnames(data1) %in% "change") == 1))#检索目的列的索引
data2 <- data1[,Target_IndexValue,drop = FALSE]#取出目的列
#data3 <- data2[complete.cases(data2), ]#删除具有缺失值NA的行（DESeq2代码中已删除）

#设置padj和log2FC的阈值
#0.05、0.01、0.001
cut_off_padj1= 0.05
cut_off_padj2= 0.01
cut_off_padj3= 0.001
cut_off_padj4= 0.0001#统计显著性
cut_off_log2FC = 1#差异倍数值

#根据阈值参数，上调基因设置为‘up’，下调基因设置为‘Down’，无差异设置为‘Stable’，并保存到change列中
# data3$change = ifelse(data3$padj < cut_off_padj1 & abs(data2$log2FoldChange) >= cut_off_log2FC, 
#                         ifelse(data2$log2FoldChange > cut_off_log2FC ,'up','down'),
#                         'none')
# head(data3)

#绘图函数
library(ggplot2)
p <- ggplot(
  #数据来源
  #xy值和颜色的映射关系，颜色取决于change列
  data2, aes(x = log2FoldChange, y = -log10(padj), colour=change)) +
  geom_point(alpha=0.4, size=2) +#点的透明度和大小
  scale_color_manual(values=c("#546de5", "#d2dae2","#ff4757"))+#点的颜色
  #边界线，数字0-6代表线的类型
  geom_vline(xintercept=c(-1,1),lty=4,col="black",lwd=0.8) +
  geom_hline(yintercept = -log10(cut_off_padj1),lty=4,col="black",lwd=0.8) +
  #坐标轴标签
  labs(x="log2(FoldChange)",
       y="-log10(p.adj)")+
  #设置图片背景
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5), #标题居中，hjust = 0.5表示水平居中对齐
        legend.position="right", #图例相对位置
        legend.title = element_blank()) +#图例无标题
  #调整坐标轴标题和坐标轴刻度标签的字号和字体样式
  theme(axis.title.x =element_text(size=14,face = "bold"), 
      axis.title.y=element_text(size=14,face = "bold"), 
      axis.text = element_text(size = 14,face = "bold"))

p

#标记特别关注的基因和显著大差异的基因
#将满足条件的基因放置在label列
library(ggrepel)
target_gene_list <- as.character(unlist(read.table("labelgene_list.txt")))
#显著大差异的基因绘制标签
#data2$label <- ifelse(data2$padj < cut_off_padj4 & abs(data2$log2FoldChange) >= 10, data2$X, "")
#感兴趣的基因绘制标签
data2$label <- ifelse(data2$X %in% target_gene_list, data2$X, "")#注意！"%in%"需要被比较者为字符或数值向量



p + geom_label_repel(
    data = data2[data2$X %in% target_gene_list,],
    aes(label = label),
    size = 3,
    segment.color = "black",
    show.legend = FALSE,
    box.padding = unit(0.5, "lines"),
    point.padding = unit(0.8, "lines"),
    max.overlaps = 20,
    )

欢迎交流，共同进步