欠小伙伴一个火山图，这次给你们补上
学习生信和绘图技能，就点上方蓝字关注我们
分析背景
火山图常用于转录组研究中，也能应用于基因组，蛋白质组，代谢组等统计数据。它归根到底是散点图的一种，它能将统计学中的统计显著性量度（如p value，矫正后的P值）和变化幅度相结合，从而能够帮助人们识别那些变化幅度较大且具有统计学意义的数据（如基因，蛋白等），如下图所示。
图形分解
图中每一点代表所研究的基因，黑色的点代表这些基因在实验组和对照组的表达情况不具有差异性，绿色和红色的点代表上调基因和下调基因。
横坐标log2FC代表差异倍数，离坐标原点越远，说明差异倍数越大
纵坐标-log10(PValue)代表差异结果的显著性，越往上的点，说明差异结果越显著。
图1B与图1A向比较，图1B将实验组和对照组中差异倍数最大的前二十个基因的名称在图中展示出来，包括十个上调基因和十个下调基因，当然在后面小编将教小伙伴们如何将感兴趣的基因（如：抑癌基因，管家基因，促癌基因）展示在图上中。
函数式编程的优势
与零碎的代码相比，反复调用，参数简单易学，对初学者更加友好。
与网页绘图小工具相比，不要要注册，不需要付费，方面经济。
对有经验的程序员，函数式编程能够提高程序员的编程水平，减少代码的重复使用，高效快捷
敲黑板，划重点
小编今天把小伙伴们转录组分析常用的火山图封闭成函数，并使用ggplot2，小伙伴们拿到函数以后，只需要将数据整理成固定的格式，代入函数中，即可绘制成高质量，高颜值的火山图，用AI小修一下，即可插入到自己的文章中。
分析方法
# 安装R包if (!requireNamespace("ggplot2", quietly = TRUE))install.packages("ggplot2",repos = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/")if (!requireNamespace("RColorBrewer", quietly = TRUE))install.packages("RColorBrewer",repos = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/")if (!requireNamespace("data.table", quietly = TRUE))install.packages("data.table",repos = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/")# 自定义函数## 快速读入数据readFlie=function(input,type,row=T,header=T){# input 为读入文件的路径，type为读入文件的类型，格式为‘.txt’或‘.csv’,row=T,将文件的第一列设置为列名library(data.table,quietly = TRUE)if(type=='txt'){dat = fread(input,header = header,sep='\t',stringsAsFactors = F,check.names = F)if(row){dat = as.data.frame(dat,stringsAsFactors = F)rownames(dat) = dat[,1]dat = dat[,-1]}else{dat = as.data.frame(dat,stringsAsFactors = F)}}else{dat = fread(input,header = header,sep=',',stringsAsFactors = F,check.names = F)if(row){dat = as.data.frame(dat,stringsAsFactors = F)rownames(dat) = dat[,1]dat = dat[,-1]}else{dat = as.data.frame(dat,stringsAsFactors = F)}}return(dat)}## 绘制火山图wn_volcano = function(Symbol,logFC,Pvalue,cutFC=1,cutPvalue=0.05,colors=NULL,ylab='-log10(FDR)',xlab='log2(FoldChange)'){require(RColorBrewer,quietly = T,warn.conflicts =F)require(ggplot2,quietly = T,warn.conflicts =F)# 整理数据state=rep('None',length(Symbol))state[which(logFC>cutFC&Pvalue<cutPvalue)]='Up'state[which(logFC< -cutFC&Pvalue<cutPvalue)]='Down'df = data.frame('logFC'=logFC,'FDR'=Pvalue,'state'=state,stringsAsFactors = F,row.names = Symbol )if (is.null(colors)){colors = c(brewer.pal(9,'Set1'))[c(3,9,1)]}else{colors = colors}# 绘图p1=ggplot(data = df, aes(x = logFC, y = -log10(FDR)))p1=p1+geom_point(alpha=0.5, size=3, aes(color=state))p1=p1+scale_color_manual(values=colors,limits = c("Down",'None', "Up"),name='State')p1=p1+geom_vline(xintercept=c(-cutFC,cutFC),lty=3,col="azure4",lwd=1)p1=p1+geom_hline(yintercept = -log10(cutPvalue),lty=3,col="azure4",lwd=1)p1=p1+ylab(ylab)+xlab(xlab)p1=p1+theme_bw()return(p1)}## 保存图片,只支持ggplot对象savePlots=function(path,plot,type=c('pdf','png','tiff')[1],width=10,height=8,dpi=300){# path表示保存图片路径,需要加上相应的文件扩展名称library(ggplot2)if(type=='pdf'){ggsave(filename = path,plot = plot,width = width,height = height,device = 'pdf')}else if(type=='png'){ggsave(filename = path,plot = plot,width = width,height = height,device = 'png',dpi = dpi)}else{ggsave(filename = path,plot = plot,width = width,height = height,device = 'tiff',dpi = dpi)}}
实战演练
扫描关注微信公众号，后台回复”火山图“，获取源代码和测试数据
# 读入数据df = readFlie('./volcano_file.txt',type = 'txt')fg = wn_volcano(Symbol = rownames(df),logFC=df$logFC,Pvalue=df$FDR,cutFC=1,cutPvalue=0.05,xlab='log2(FoldChange)',colors = c('#34e610','#495057','#ff1616'))# 展示图片fg# 保存图片savePlots(path = './fg.pdf',plot = fg,type = 'pdf',width = 10,height = 8)savePlots(path = './fg.png',plot = fg,type = 'png',width = 10,height = 8,dpi = 300)savePlots(path = './fg.tiff',plot = fg,type = 'tiff',width = 10,height = 8,dpi = 600)
高级绘图