一文做会漂亮的火山图

一、通过limma包对输入数据进行处理

1、归一化处理

在利用limma包进行差异分析处理之前，要对数据进行归一化处理：

输入文件1

在使用limma包之前首先要对数据进行标准化处理

rm(list = ls())
setwd('/lab412C/LSM/RNA-SEQ/Volcant')
data1 <- read.csv(file = "data2.csv",header = T,sep = ",")
rownames(data1)<-data1[,1] 
data2<-data1[,-1]
qx <- as.numeric(quantile(data2, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <- (qx[5] > 100) ||
  (qx[6]-qx[1] > 50 && qx[2] > 0) ||
  (qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)
LogC
data <- log2(data2[,]+1)
data <- data[which(rowSums(data) > 0),]

如果在这一步你没有对数据进行归一化等处理，那么极大值就会掩盖极小值且火山图无法显示出喷射状态， 就是中间的分叉分不开，非常难看！
示例——经过归一化处理之后的数据：

归一化后的数据

2、limma包的差异分析处理

group_list <- c('c','c','c','E','E','E')
library(limma)
design=model.matrix(~factor(group_list))
fit=lmFit(data,design)
fit=eBayes(fit)
deg=topTable(fit,coef=2,number = Inf)
write.csv(x=deg,file='/lab412C/LSM/RNA-SEQ/Volcant/deg1.csv')

数据展示

注意：
这里的x名称应该为gene,所以可以先导出来csv文件编辑一下 在输入，也可以自己用R再设置一下就可以了
；因为后续我们需要对deg$gene这一列进行处理，不建议用X这种名称直接进行处理哦！

deg1 <- read.csv(file = "deg1.csv",header = T,sep = ",")
deg <- deg1
logFC_t=1 #不同的阈值，筛选到的差异基因数量就不一样，后面的超几何分布检验结果就大相径庭。
change=ifelse(deg$P.Value>0.05,'stable', 
              ifelse( deg$logFC >logFC_t,'up', 
                      ifelse( deg$logFC < -logFC_t,'down','stable') )
)

接下来的数据应该是这样的：

二、作图

这一部分主要包括两部分，首先是对adj.P.Val取对数，另外需要根据logFC的标准定义gene的上下调

logFC_t=1 #不同的阈值，筛选到的差异基因数量就不一样，后面的超几何分布检验结果就大相径庭。
change=ifelse(deg$P.Value>0.05,'stable', 
              ifelse( deg$logFC >logFC_t,'up', 
                      ifelse( deg$logFC < -logFC_t,'down','stable') )
)
deg$logP <- -log10(deg$adj.P.Val)
library(ggpubr)
library(ggthemes)
ggscatter(deg,x='logFC',y='logP')+theme_base()

deg <- mutate(deg,change)
table(deg$change)
ggscatter(deg, x = "logFC", y = "logP",color = "change",palette = c("#9999FF", "gray" , "#FF9999"),size=1 )+ theme_base()
##加分界线：
ggscatter(deg, x = "logFC", y = "logP",color = "change",palette = c("#9999FF", "gray" , "#FF9999"),size=1 )+ theme_base()+
  geom_hline(yintercept = 0.43 , linetype ="dashed")+ 
  geom_vline(xintercept = c(-1,1), linetype= "dashed")

##加gene_name
deg$label= ""
deg <- deg[order(deg$adj.P.Val), ]
up.gene <- head(deg$gene[which(deg$change=="up")],10)
down.gene <- head(deg$gene[which(deg$change=="down")],10)
deg.top10.genes <- c(as.character(up.gene),as.character(down.gene))
deg$label[match(deg.top10.genes,deg$gene)] <- deg.top10.genes

ggscatter(deg, x = "logFC", y = "logP",color = "change",palette = c("#9999FF", "gray" , "#FF9999"),
          size=1,
          lable = deg$label,
          font.label = 8,
          repel = T ,
          xlab = "log2FoldChange",
          ylab = "-log10(Adjust P-value)",)+ theme_base()+
  geom_hline(yintercept = 0.43 , linetype ="dashed")+ 
  geom_vline(xintercept = c(-1,1), linetype= "dashed")

结果：

记得点赞分享哦！

以上参考一篇分有意思的推文:
https://cloud.tencent.com/developer/article/1512442