随着测序技术的提升及成本的不断降低，基因表达谱数据也更精确及更容易获得，测几个样本做个差异表达分析，接下来就是GO功能及KEGG信号通路分析，如果觉得工作量不够的话，再做个蛋白相互作用（PPI）分析，不得不说这些已经成为套路，不，标准流程了。今天主要介绍这些套路的可视化及批量化。
能够做GO、KEGG分析的工具太多了，但其实在做的时候有一两个好用的够用就行，我个人喜欢的两个在线的工具如webgenstal(http://www.webgestalt.org)，DAVID (https://david-d.ncifcrf.gov/)。前者05年发布在NAR上，13年更新，被引超过1000；后者09年发表在NAR上，被引超过3000，虽然DAVID的界面及其丑陋，但不能掩盖其实用的功能。这两工具使用起来都是导入目的基因，背景基因列表（如果有），然后选择要分析的功能，选择统计学阈值，导出结果文件即可。
接下来就是对结果的可视化，使用恰当美观的图来展示数据是生信分析的一个重要步骤，相信没谁愿意面对黑压压的一堆数字。Webgestal在分析后会自动生成图片，如下面的图：
但这样就会有个问题，比如我只想在GO的BP结果中展示与肿瘤相关的或只想展示与代谢相关的p<0.05的通路，其他结果直接去掉，这就需要先把需要的通路挑出来，然后再手动作图，而不是自动生成的图。在可视化工具的选择上面， R语言ggplot2包将数据统计和作图傻瓜式的结合了起来，自动配色，高清图像导出，批量作图。我一般习惯从DAVID获得GO、KEGG结果列表，然后挑选要展示的通路，保存为tsv格式文件，如下：
然后导入R作图，代码及图如下：
rm(list=ls())
library(ggplot2)
library(Cairo)
setwd("/home/ ")
GO_BP <- read.table("./enh_statistics/A549_GO_BP_spe.tsv",header = T,sep="\t")
png_path="./figure/GO_BP.png"
CairoPNG(png_path, width = 12, height = 7, units='in', dpi=600)
ggplot(data=GO_BP)+
geom_bar(aes(x=reorder(Term,Count),y=Count, fill=-log10(PValue)), stat='identity') +
coord_flip() +
scale_fill_gradient(expression(-log["10"](P.value)),low="red", high = "blue") +
xlab("") +
ylab("Gene count") +
scale_y_continuous(expand=c(0, 0))+
theme(
axis.text.x=element_text(color="black",size=rel(1.5)),
axis.text.y=element_text(color="black", size=rel(1.6)),
axis.title.x = element_text(color="black", size=rel(1.6)),
legend.text=element_text(color="black",size=rel(1.0)),
legend.title = element_text(color="black",size=rel(1.1))
# legend.position=c(0,1),legend.justification=c(-1,0)
# legend.position="top",
)
dev.off()
Kegg输入文本格式如下：
Kegg气泡图代码如下：
rm(list=ls())
library(Cairo)
library(stringr)
setwd("/home/")
pathway = read.table("./enh_statistics/A549_KEGG.tsv",header=T,sep="\t")
pathway$Term<-str_split_fixed(pathway$Term,":",2)[,2]
png_path="./figure/KEGG.png"
CairoPNG(png_path, width = 5.9, height = 3, units='in', dpi=600)
ggplot(pathway,aes(x=Fold.Enrichment,y=Term)) +
geom_point(aes(size=Count,color=-1*log10(PValue)))+
scale_colour_gradient(low="green",high="red")+
labs(
color=expression(-log[10](P.value)),
size="Gene number",
x="Fold enrichment"
# y="Pathway name",
# title="Pathway enrichment")
)+
theme_bw()+
theme(
axis.text.y = element_text(size = rel(1.3)),
axis.title.x = element_text(size=rel(1.3)),
axis.title.y = element_blank()
)
dev.off()
Kegg的气泡图结果如下：
如果实在用不来代码，也可以使用在线工具WeGO(http://wego.genomics.org.cn/cgi-bin/wego/index.pl) 进行绘制，结果是这样的:
除了可以对基因进行GO分析外，还可以对基因组原件如沉默子、绝缘子、增强子等顺势调控元件进行GO注释，另外，也可以对Chip-seq结果的候选保守性功能元件sequence motif进行GO分析。如果需要对N个组织或样本的数据进行GO分析，那最好使用bioconductor的R包批量实现，这些将在后续陆续整理发布。
写在最后，关于可视化作图，其实matlab的可视化功能比R要更加强大，但matlab并非免费，于是开源、功能较为全面的R就备受推崇。但不得不说的是相比于前面提到的作图工具，Javascript的图表库是最漂亮的，有兴趣的同学可以搜下ECHARTS，是百度的良心作。
下面这张图是ggplot2的作者Hadley Wickham，让我们记住这个美国男人。
网友评论