差异分析后画火山图
> library(ggplot2)
> library(ggrepel)
> res = na.omit(res)#去NA
> res$type = ifelse(res$log2FoldChange >= 1.5 & res$padj <=0.001,"up_regulated", ifelse(res$log2FoldChange<= -1.5 & res$padj <=0.001,"down_regulated","none siginificant"))
> table(res$type)#列出上调、下调基因的数目
> title = paste0('The number of up gene is ', nrow(res[res$group == 'up_regulated',]),'\nThe number of down gene is ', nrow(res[res$group == 'down_regulated',]))
> res = as.data.frame(res)#不转换会出现报错:`data` must be a data frame
> ggplot(res, aes(x = log2FoldChange, y = -log10(padj)))+geom_point(aes(color = group)) + scale_color_manual(values = c("purple","grey","navy blue"),limits = c('up_regulated','none siginificant'," down_regulated")) + theme_bw(base_size = 20) + theme(plot.title = element_text(size=15,hjust = 0.5) + theme_classic())+ggtitle(title)
上述代码可画出最简单的火山图,如果想加点花里胡哨的东西↓
> p = ggplot(res, aes(log2FoldChange, -log10(padj),color = type))
> p +geom_point()
> x_lim = max(res$log2FoldChange, - res$log2FoldChange)#设置x轴坐标
> gg = p + geom_point(size =1) + xlim(-x_lim, x_lim) +labs(x= "log2FoldChange", y = "-log10(padj)") +scale_color_manual(values = c("#A52A2A","grey","#f8766d")) + geom_hline(aes(yintercept=-1*log10(0.05)),colour="black", linetype="dashed") + geom_vline(xintercept=c(-2,2),colour="black", linetype="dashed")
> print(gg)