转录组-差异基因火山图

##增强版火山图1##

getwd()
gene_info<-read.csv("gene_info.csv")#gene信息表
gene_exp<-read.csv("gene_exp.csv")#行名=sample_info列名；gene表达量数据表
de<-read.csv("de.csv")#差异基因表

library(tidyverse)
de_result<-left_join(de,gene_info,by=c("id"="GID"))#关联基因信息表
de_result2<-left_join(de_result,gene_exp,by=c("id"="GID"))%>%#关联基因表达表
  arrange(desc(abs(log2FoldChange)))#按abs(log2FoldChange)从大到小/小到大？排序


#查看差异基因总数
group_by(de_result2,direction)%>%
  summarise(count=n())
#install.packages("EnhancedVolcano")
#BiocManager::install("EnhancedVolcano")
library(EnhancedVolcano)
library(ggplot2)


#BiocManager::install("ggrepel")
library(ggrepel)
EnhancedVolcano(de_result2,
                lab=de_result2$id,
                x="log2FoldChange",
                y="padj",
                FCcutoff=2,#改变FC阈值，默认为1
                pCutoff=0.01)#改变pvalue阈值，默认为0.05
#保存到csv
write.csv(de_result2,
          file="de_result2.csv",
          row.names = F,quote = F)
save(de_result2,file="de_result3.rdate")

##火山图2##

getwd()
de1<-read.csv("diff_detail.csv")#全部基因信息表
library(tidyverse)
#添加调控信息，添加一列为regulate
de1 %>% mutate(regulate = case_when(log2FoldChange>1&padj<0.05 ~ "up",
                                    log2FoldChange<(-1)&padj<0.05 ~ "down",
                                    TRUE ~ "ns")) ->de1

#统计个数
table(de1$regulate)
#ggplot2绘图
library(ggplot2)
ggplot(de1,aes(log2FoldChange,-log10(padj),color = regulate))+
  geom_point()+
  scale_color_manual(values = c("blue", "gray", "red"))+
  xlim(-10,10)