【R实战】如何优雅智能的画出带标注的火山图
2021-12-16 本文已影响0人
xizzy
本来说这期是要更新LncRNA鉴定思路的，但是看到上一期人气低迷，就暂时晚点在弄了，这次就来跟大家分享一下火山图
什么是火山图？
顾名思义，火山图就是一个长得有点像火山的图形。其本质就是一个散点图，通常用来进行展示基因表达量。所以绘制一个火山图我们需要以下几个关键的值：
1.基因的表达量；
2.基因的显著值（Pvalue/FDR）；
差异peak表达量图.png
如上图所示展示的是一个exomePeak2进行差异分析后得到的结果，我们应该如何对这组数据绘制火山图呢。首先我们定位好两个关键的列
DiffModLog2FC就是表达量，pvalue/padj就是显著值。明白了这个我们就可以开始绘制火山图了。
我们直接上代码
library(ggplot2)
library(ggrepel)
library(optparse)
#作者：xizzy 公众号：生信咖啡
#为了方便我们以后一劳永逸的使用，把常用的一些作图参数用传参的方式进行配置
option_list = list(
make_option( c('-i','--expression'), type = 'character', help='Input your gene expression file'),
make_option( c('-n','--fcname'),type = 'character', help = 'Input your fc column name. Attetion () in R would be convert to .', default = 'log2.FC.' ),
make_option( c('-f','--foldchange'), type = 'double', default = 1, help = 'Input your foldchange choose level, default is 1' ),
make_option( c('-s','--signame'), type = 'character', help = 'Input your filter method. For example: Pvaule / FDR, it must be same with column name.' ),
make_option( c('-c','--cutoff'), type = 'double', help = 'Input your cutoff level, default is 0.05', default = 0.05 ),
make_option( c('-a','--annofile'), type = 'character', help = 'Input your annote file, if use this paramter volcano plot would mark the annote file point. This paramter is option.', default = 'no_this_file' ),
make_option( c('-o','--output'), type = 'character', help = 'Input the output prefix, default is ./volcano' , default = './volcano'),
make_option( '--ymax', type = 'integer', help = 'Y continuous max, default is 50', default = 50),
make_option( '--xmax', type = 'integer', help = 'X continuous max, default is 10', default = 10)
);
opt_parser = OptionParser(option_list=option_list);
opt = parse_args(opt_parser);
#读取我们的数据，如果你数据是csv就改分隔符为逗号
DE_data=read.table(file=opt$expression,
sep="\t",header=TRUE)
if( opt$annofile != 'no_this_file' ){
mark_df =read.table(opt$annofile,sep='\t',header=T) #标注特定点的文件，样式同DE_data
}
#配置一波显著水平，将数据分为上调、下调、不显著三种，并计算数量
DE_data$significant = 'No-DEGs'
DE_data$significant[DE_data[,opt$signame] < opt$cutoff & DE_data[, opt$fcname] < -opt$foldchange ] = 'down'
DE_data$significant[DE_data[,opt$signame] < opt$cutoff & DE_data[, opt$fcname] > opt$foldchange ] = 'up'
up_count = nrow(subset(DE_data, significant == 'up'))
down_count = nrow(subset(DE_data, significant == 'down'))
no_count = nrow(subset(DE_data, significant == 'No-DEGs'))
DE_data$significant[DE_data$significant == 'down'] = paste('Down:',down_count)
DE_data$significant[DE_data$significant == 'up'] = paste('Up:',up_count)
DE_data$significant[DE_data$significant == 'No-DEGs'] = paste('No-DEGs:',no_count)
#开始绘图
DE_data$logsig = -log10(DE_data[,opt$signame])
r = ggplot(DE_data,aes_string(x = opt$fcname , y = 'logsig' )) + theme_set(theme_bw(base_size = 8/.pt)) + theme_gray(base_family = "arial") #首先是映射数据，x轴画表达量，y画显著水平，并修改了一下字体
rp2 = r + geom_hline(yintercept = -log10(opt$cutoff),linetype = 4)+
geom_vline(xintercept = c(-foldchange,foldchange),linetype=4) #添加显著水平的线
rp3 = rp2 + geom_point(aes(color=significant),alpha = 0.4,size = exp(1)^(-log10(nrow(DE_data)))*20)+ #根据数量的多少自动调整点的大小，这里用了一个单调递减的函数去调控
theme_bw()+
scale_color_manual(values = c('#3288BD','gray','#D53E4F'))+
theme(panel.grid.major =element_blank(),panel.grid.minor=element_blank()) #去背景颜色，设置点的颜色
rp4 = rp3 +
guides(colour = guide_legend(override.aes = list(size=2))) +
theme(axis.ticks.y = element_line(color = c(NA, NA, NA, NA))) #固定图例中点的大小
rp4 = rp4 + theme(panel.border = element_blank()) +
theme(axis.line = element_line(colour = 'black')) #去边框，并为加上x和y轴的线条，这点在很多好的期刊都是这样要求的
rp4 = rp4 + xlab(bquote(log[2]~(Fold~ ~Change))) + ylab(bquote(-log[10]~(Significance))) + scale_y_continuous(limits=c(0,opt$ymax)) + scale_x_continuous(limits=c( -opt$xmax, opt$xmax)) #添加x和y轴的名字，并配置坐标范围
if( opt$annofile != 'no_this_file' ){ #如果有特别需要标注的点的数据，在图中标识出来
mark_df$logsig = -log10(mark_df[,opt$signame])
rp5 = rp4 + geom_label_repel(data = mark_df ,aes_string( label = opt$annoname),
color="black",size=6/.pt,
box.padding=unit(0.6, "lines"),
point.padding=unit(0.6, "lines"),
segment.colour = "grey50")
}
#保存图片
ggsave( paste0(opt$output,".png"), width = 5, height = 4)
ggsave( paste0(opt$output,".pdf"), width = 5, height = 4)
结果展示
让我们实际运行一下看看，使用命令行
Rscript Volcano.r -i test.xls -n DiffModLog2FC -s pvalue --ymax 10
输入文件为test.xls, 表达量的列名叫做DiffModLog2FC, 显著值的列名用pvalue，y的最大值限制到10
火山图.png
这样一张可以发表水平的火山图就完成了，下面在演示一下如何标注关键基因。
这里为了演示，我就把其中的一个name改为TP53然后让它标识出来，文件如下：
注释文件展示
使用代码：
Rscript Volcano.r -i test.xls -n DiffModLog2FC -s pvalue --ymax 10 --annofile mark.xls --annoname name
带标注的火山图
至此，优雅的自动化发表级，可带标注的火山图绘制完成！
