====== 多比较组火山图 ======
#!/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')
{{:个性化条目:多比较组和火山图.png?400|}}
新增部分:
测试路径:
/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
{{:compares_volcano.png?400|}}