火山图绘制

一些作图笔记，防止自己忘记，以下得到的图结果是添加辅助线条，标记基因

args <- commandArgs(trailingOnly = TRUE)
diffdir <- args[1]
outdir <- args[2]
library(ggplot2)
library(ggrepel)
theme_set(theme_bw())
qvalue_limit = as.numeric(0.05)
## 示例路径
#1.diffpath:差异比较组，组别之间通过VS连接，比如A-vs-B
#2.outdir:Volcano_path/火山图的输出路径
diffpath <- diffdir
filepath <- list.files(path = diffpath, pattern = "-vs-", full.names = TRUE)

for (i in filepath) {
  infile <- paste0(i, "/DiffExpr.SYMBOL.csv")
  vsname <- basename(i)
  data <- read.csv(file=infile, header=T, sep=",", check.names=F)
  names(data)[1] <- "SYMBOL"
  data$stat <- factor(data$stat, levels=c("down-regulated", "up-regulated", "non-significant"))
  ## get up_10 ordered by Pvalue
  up <- data[which(data$stat == "up-regulated"),1][1:10]
  up_10genes <- data[which(data$SYMBOL %in% up), ]
  ## get down_10 ordered by Pvalue
  down <- data[which(data$stat == "down-regulated"),1][1:10]
  down_10genes <- data[which(data$SYMBOL %in% down), ]
  ## get deg numubers
  deg <- data[which(data$stat != "non-significant"), 1]
  down_num <- length(which(data$stat=="down-regulated"))
  up_num <- length(which(data$stat=="up-regulated"))
  deg_num <- down_num + up_num
  down_lab <- paste("down-regulated (",down_num,")", sep="")
  up_lab <- paste("up-regulated (",up_num,")", sep="")
  label <- paste0(vsname,"\nTotal Genes: ", deg_num,  sep="")
  p <- ggplot(data, aes(x=logFC, y=-log10(PValue)))
  p <- p + geom_point(aes(colour=stat)) + labs(title = label)
  p <- p+scale_colour_manual(values=c("blue","#FF4040","grey"),
                             breaks=c("down-regulated","up-regulated","non-significant"),
                             labels=c(down_lab,up_lab,"Not significant"))
  p <- p + labs(x="log2(Fold Change)", y="-log10(PValue)")
  p <- p + theme(plot.title=element_text(size=15, margin=margin(b=10), hjust=0.5, face="bold"))
  p <- p + theme(axis.title.x=element_text(size=12, face="bold"))
  p <- p + theme(axis.title.y=element_text(size=12, face="bold"))
  p <- p + theme(legend.title=element_text(size=8, face="bold"))
  p <- p + theme(legend.text=element_text(size=6))
  p<-p+geom_hline(yintercept=-log10(qvalue_limit),linetype="dotdash",size=0.4)
  p<-p+geom_vline(xintercept=1, linetype = 'dotdash', size = 0.4, alpha = 0.5)
  p<-p+geom_vline(xintercept=-1, linetype = 'dotdash', size = 0.4, alpha = 0.5)
  p<-p+scale_x_continuous(limits = c(-max(abs(data$logFC))-0.1, max(abs(data$logFC))+0.1))
  p<-p+theme(legend.position=c(0.8,0.9), legend.key = element_rect(fill = NULL), legend.background = element_blank(),legend.text = element_text(size=10))
  p<-p+theme_test()+theme(plot.margin = margin(0.5,0.5,0.5,0.5, unit = "cm"))
  p<-p+theme(plot.title=element_text(size=12,hjust=0.5),plot.subtitle=element_text(size=10,hjust=0.5))
#标记上调基因
  p <- p + geom_text_repel(data = up_10genes,label = up_10genes$SYMBOL ,aes(x=logFC, y=-log10(PValue)),size = 3.0, color = 'black',
                           show.legend = FALSE,max.overlaps = Inf,nudge_x = 0.7,  min.segment.length = 0.5)
标记下调基因
  p <- p + geom_text_repel(data = down_10genes,label = down_10genes$SYMBOL ,aes(x=logFC, y=-log10(PValue)),size = 3.0, color = 'black',
                           show.legend = FALSE,max.overlaps = Inf,nudge_x = -0.7, min.segment.length = 0.5)
  filepng <- paste0(outdir, vsname, ".volcano.png")
  filepdf <- paste0(outdir, vsname, ".volcano.pdf")
  ggsave(filename=filepng, units="px", width=2220, height=1800, dpi=300)
  ggsave(filename=filepdf,  width=10, height=7, device="pdf")
}

以上是在多组情况下进行循环作图，如果只有一个比较组，可以不用循环，结果图如下，如果继续绘制其他形式 的图，会及时补充