本来说这期是要更新LncRNA鉴定思路的，但是看到上一期人气低迷，就暂时晚点在弄了，这次就来跟大家分享一下火山图

什么是火山图？

顾名思义，火山图就是一个长得有点像火山的图形。其本质就是一个散点图，通常用来进行展示基因表达量。所以绘制一个火山图我们需要以下几个关键的值：
1.基因的表达量；
2.基因的显著值（Pvalue/FDR）；

差异peak表达量图.png

我们直接上代码

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

注释文件展示

带标注的火山图

至此，优雅的自动化发表级，可带标注的火山图绘制完成！

还不赶快关注我，以及我的公众号生信咖啡吗？