#!/TJPROJ2/GB/PUBLIC/software/GB_TR/mRNA/miniconda3/envs/r_3.5.0/bin/Rscript

suppressMessages({
  library(stringr)
  library(ggplot2)
  library(dplyr)
  library(RColorBrewer)
  #install.packages("ggnewscale")
  library(ggnewscale)
  library(argparser)})

argv <- arg_parser('绘制多组差异火山图')
argv <- add_argument(argv,"--diffdir", help="the differential analysis results dir")
argv <- add_argument(argv,"--compares", help="the compare names")
argv <- add_argument(argv,"--fc", help="the foldchange vlaue")
argv <- add_argument(argv,"--pvalue", help="the p value")
argv <- add_argument(argv,"--padj", help="the p adjust value")
argv <- add_argument(argv,"--outdir", help="the output dir")
argv <- parse_args(argv)

diffdir <- argv$diffdir
compares <- argv$compares
fc <- argv$fc
p <- argv$pvalue
q <- argv$padj
outdir <- argv$outdir

if (!is.na(fc)){fc <-as.numeric(fc)}
if (!is.na(p)){p <-as.numeric(p)}
if (!is.na(q)){q <-as.numeric(q)}

compare_names <- unlist(str_split(compares, pattern = ','))
compare_names <- unlist(str_split(compares, pattern = ','))

volcano_data <- data.frame()
y_max <- c()
for (compare in compare_names){
  deg_file <- paste(diffdir,"/",compare,"/",compare,"_deg.xls", sep = '')
  deg <- read.table(deg_file, sep = '\t', header = T, quote = '',fill = T)
  tmp_volcano_data <- select(deg, gene_id, log2FoldChange, pvalue, padj)
  tmp_volcano_data$compare <- compare
  volcano_data <- rbind(volcano_data, tmp_volcano_data)
  y_max <- c(y_max, max(abs(deg$log2FoldChange)))
}

volcano_data$diff <- 'no_diff'
if (is.na(p)) {
  volcano_data$diff[abs(volcano_data$log2FoldChange) >=log2(fc) & volcano_data$padj <=0.05] <- 'diff'
}
if (is.na(q)){
  volcano_data$diff[abs(volcano_data$log2FoldChange) >=log2(fc) & volcano_data$pvalue <=0.05] <- 'diff'
}

data2 <- data.frame(x = compare_names,
                    y=0,
                    label = compare_names)

datbar<-data.frame(x=compare_names,
                   y=y_max)
if (length(compare_names) <= 12){
  color = brewer.pal(length(compare_names), "Set3")
} else {
  color = rainbow(length(compare_names))
}

P <- ggplot()+
  geom_col(data=datbar,aes(x=x,y=y),fill="grey",alpha=0.5)+
  geom_col(data=datbar,aes(x=x,y=-y),fill="grey",alpha=0.5)+
  geom_jitter(data=volcano_data,
              aes(x=compare,y=log2FoldChange,
                  color=diff))+
  scale_color_manual(name=NULL,
                     values = c("red","darkgrey"))+
  ggnewscale::new_scale_fill()+
  theme_minimal()+
  theme(axis.line.y = element_line(),
        axis.ticks.y = element_line(),
        panel.grid = element_blank(),
        legend.position = "top",
        legend.justification = c(1,0),
        legend.direction = "vertical",
        axis.text.x = element_blank())+
  labs(x="Compare",y="log2FoldChange")+
  geom_tile(data=data2,
            aes(x=x,y=y,fill=x),
            height=3,color="black",
            alpha=0.9,
            show.legend = F)+
  scale_fill_manual(values = color, "Paired")+
  geom_text(data=data2,aes(x=x,y=y,label=label))
  
pdf(file=paste(outdir,'compares_volcano.pdf', sep = '/'))
print(P)
dev.off()
svg(filename=paste(outdir,'compares_volcano.svg', sep = '/'))
print(P)
dev.off()
ggsave(file=paste(outdir,'compares_volcano.png', sep = '/'),type='cairo-png',plot=P,
       device ='png', width =7, height = 7,units = 'in',dpi = 300,
       bg='white')

新增部分：

测试路径： /TJPROJ6/RNA_SH/personal_dir/fengjie/Personal_analysis/volcano_multi

/TJPROJ2/GB/PUBLIC/software/GB_TR/mRNA/miniconda3/envs/r_3.5.0/bin/Rscript volcano_multi --das /TJPROJ7/GB_TR/PJ_GB/mRNA/med/2011/andong/X101SC23063193-Z01-F001.B1-16_Homo.20230825/Differential/1.deglist/HSYAvsSF/HSYAvsSF_deg.xls,/TJPROJ7/GB_TR/PJ_GB/mRNA/med/2011/andong/X101SC23063193-Z01-F001.B1-16_Homo.20230825/Differential/1.deglist/HFDvsSF/HFDvsSF_deg.xls,/TJPROJ7/GB_TR/PJ_GB/mRNA/med/2011/andong/X101SC23063193-Z01-F001.B1-16_Homo.20230825/Differential/1.deglist/HSYAvsHFD/HSYAvsHFD_deg.xls --fc 2 --pvalue 0.05 --outdir /TJPROJ6/RNA_SH/personal_dir/fengjie/Personal_analysis/volcano_multi