2024-01-17小众物种的go分析
我本身研究的植物品种都比较小众，比如Brachypodium distachyon和Setaria italica，还有一些genome都不是很好的品种，在进行GO analysis时就比较困难。一般情况下有两种策略，一是直接使用该物种现有的GO annotation进行分析，另一种是根据蛋白质序列找到相近annotation比较好物种的homologs再进行分析。接下来就来详细讲一讲以上两种策略的实现方法。
一是使用该物种现有的GO annotation。对于植物来说，annotation比较齐全的有以下两个database：
agriGOv2.0 http://systemsbiology.cau.edu.cn/agriGOv2/index.php
PlantGSEA http://structuralbiology.cau.edu.cn/PlantGSEA/index.php
通常来讲，可以直接在网站在线进行分析，只需要上传gene list和设定相关参数就可以。需要注意的是，这两个网站都有gene ID格式的要求，要根据相应要求找到对应版本的gene ID。一般情况下都支持Phytozome最新版本的gene ID。如果不相符，可以重新mapping到支持版本的genome上。想要下载相应的数据，只需找到downloads即可。
agriGO的功能和物种都相对齐全，它整合了一些常用的下游网站。我一般常用以下两个功能：SEA和REVIGO。SEA是普通的GO enrichment analysis（Fig. 1），结果以列表的形式给出，同时也可以使用网站自带的可视化功能。REVIGO可以将缩减很长的GO结果，同时给出semantic similarity-based scatterplots，结果以图片的形式给出，同时也可以保存相应的R script在R中进行调整。REVIGO也可以使用单独的网站 http://revigo.irb.hr/ 。
图片.png
Fig. 1 agriGO v2.0 部分页面
PlantGSEA用法和agriGO相似（Fig. 2），右边表格里的物种都可以进行基本的GO analysis。除此之外，它还集合了一些其他的database，表格中显示为绿色的就是存在的data。
图片.png
Fig. 2 PlantGSEA页面
除了在线使用，下载的GO annotation也可以在R package里使用，比如clusterProfiler。
library("clusterProfiler")
data <- read.table("GO_annotation.txt",header = T,sep="\t")
gene <- read.table("gene_list.txt")
go2gene <- data[, c(2, 1)]
go2name <- data[, c(2, 3)]
go <- enricher(c(gene)$GeneId,TERM2GENE=go2gene,TERM2NAME=go2name)
go
图片.png
gene_list.txt格式如下
图片.png
如果使用本物种GO annotation结果不太满意，还可以blastp到相近annotation比较好的物种，比如拟南芥和水稻。通过本地blast实现比较快。如果有需要可以找我，我们找一期讲一讲。
除此之外，PANTHER网站也有比较丰富的植物GO annotation，它的优点是可视化比较漂亮。如果对图片要求不是很高，可以直接使用网站生成的图片。
http://www.pantherdb.org/geneListAnalysis.do
图片.png
Fig. 3 PANTHER页面
参考内容
-
Xin Yi, Zhou Du, Zhen Su. (2013) PlantGSEA: a Gene Set Enrichment Analysis toolkit for plant community. ***Nucleic Acids Research. ***doi:10.1093/nar/gkt281.
-
agriGO v2.0: Tian Tian, Yue Liu, Hengyu Yan, Qi You, Xin Yi, Zhou Du, Wenying Xu, Zhen Su; agriGO v2.0: a GO analysis toolkit for the agricultural community, 2017 update. Nucleic Acids Res 2017 gkx382. doi: 10.1093/nar/gkx382
-
使用clusterProfiler进行GO富集分析 https://www.jianshu.com/p/eac402c603d9
参考：https://mp.weixin.qq.com/s/8SQO0VqwvTXje9rGXv7DcQ
以上提到了如何做GO分析，这一篇来说一下怎么画简单的GO图。我们先来看看通过agriGO分析得出的结果，在图中这些内容都可以表达出来。对于FDR，建议在画图前进行简单处理，一般是-log10（FDR）。
图片.png
如果你不太熟悉R，或是想更直观出图，可以参考这篇文章中的在线工具。在线画图工具 除此之外，还有一个比较实用简便的本地工具TBtools，可以参考这篇文章。在这个工具里，除了可以使用已有结果画图，也可以直接进行GO富集分析。我的挣扎 与 TBtools 的开发
如果以上的方法都没法满足你的需求，可以试试使用R。对于简单的条形图（Fig. 1），甚至可以选择直接使用excel或者AI，方便调整颜色和坐标轴。如果想还原下图也非常简单，可以参考以下代码。输入的数据也可以参考下边的表格（Table 1）。
library(ggplot2)
t <- read.csv("GO.csv", header=T)
t$Term=factor(t$Term,levels=t$Term)
p=ggplot(t, aes(x=Term,y=log10_p_value, fill=Ontology))
#下一步是确定图表类型bar，翻转x和y坐标轴生成横向的条形图，确保是逆时针90°翻转
pbar=p+geom_bar(stat="identity")+coord_flip()+scale_x_discrete(limits=rev(levels(t$Term)))
#下一步是根据term分类来填充颜色，这三个颜色都是默认的
pr=pbar+scale_fill_discrete(name="Ontology",breaks=c("Biological process","Molecular function","Cellular component"))
#下一步去掉了边框和背景，调整了字体大小，确定了y轴坐标尺
p1= pr+theme_bw()+ scale_y_continuous(breaks=NULL)+ylim(c(0,20))+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(),panel.border = element_blank())+theme(axis.text.y = element_text(size = 15))+theme(axis.text.x = element_text(size = 15))+theme(axis.title.x = element_text(size = 15))+theme(axis.title.y = element_text(size = 15))+theme(legend.text = element_text(size = 15), legend.title = element_text(size = 15))
p1
图片.png
Fig. 1 我做过的其中的GO简单条形图。
Table. 1 GO.csv格式
图片.png
如果想画一个表达内容比较多的图片（Fig.2 ，Table. 2），可以试试下边的内容。Rich factor是overlapped基因数和此term所有基因数的比值。
library(ggplot2)
#row names是Class，GO_Name，GO_ID，Rich_factor，Number_of_Genes，q_value
t <- read.table("GO.csv",header = T,sep = ",")
x=t$Rich_Score
y=factor(t$GO_Name,levels = t$GO_Name)
p = ggplot(t,aes(x,y))
p1 = p + geom_point(aes(size=Number_of_Genes,color=q_value,shape=Class,))+scale_color_gradient(low = "olivedrab3", high = "DeepPink")+labs(x="Rich Factor",y="Pathway Name")+theme_bw()+theme(panel.grid.minor = element_blank(),panel.background = element_blank())
图片.png
Fig. 2 表达内容稍多的GO图
Table 2 GO.csv
图片.png
参考内容：
http://blog.sina.com.cn/s/blog_1704ff73a0102wtx4.html
https://www.jianshu.com/p/59428e69da05
上次说到怎么做GO分析，在做完GO分析后往往要专注某个或某些基因。如果你有了感兴趣的GO term，可以看一下里边包含了哪些基因。想要看这些基因的具体annotation，找到在这个pathway里重要的marker基因，有下边几个网站可以一定程度上减轻你的工作量。
1. PLAZA
https://bioinformatics.psb.ugent.be/plaza/versions/plaza_v4_5_monocots/
这个网站整合了多种单子叶和双子叶植物的基因组数据，还包含了一些物种不同测序品种的基因组。在这个网站可以搜索基因，GO term，binding sites和gene family。如果直接搜索某个基因名称，可以搜索出不同物种里该基因的情况。这个网站最大的优势是资源整合齐全，不用为了找某个基因连续换网站。而且网站更新较快，现在已经更新到4.5版本，界面也非常友好。
图片.png
图片.png
图片.png
**2. RIKEN FLcDNA database **
二穗短柄草http://brachy.bmep.riken.jp/ver.1/index.pl
大豆 http://spectra.psc.riken.jp/menta.cgi/rsoy/index
这系列的网站本质和Phytozome或者其他网站没有区别，但是网站做得很有意思，可以在搜索对话框里直接输入非常多基因ID，只要基因ID之间有空格就可以。结果也可以整批输出，不用一个一个输入还是很方便的。而且可以同时看到在几个主要database里这个基因的annotation，也能加快筛选速度。
图片.png
图片.png
除了以前介绍过的植物专门的GO分析网站，经常用来分析蛋白互作的网站STRING也可以进行GO分析。
https://string-db.org/
在Search中可以选择Multiple sequences或Multiple protein，直接复制基因ID，蛋白名称或蛋白序列都可以。但是每次最多输入2000个基因。
图片.png
结果一般以network的形式呈现，交互式的结果方便拖拽观察蛋白之间的关系，也可以直接点击单个蛋白查看其基本信息。
图片.png
如果想提高蛋白之间关系的可信度，可以在Settings里进行选择。可以选择edge之间是如何如何连接的，比如是来源于实验或者数据库。
图片.png
比较重要的筛选数值是minimum required interaction score，这个数值来源于KEGG database，数值越大表示可信度越高，同时保留的相互作用越少。默认选择- 0.4. 如果想提高可信性，可以选择-0.7，通常情况下选择-0.7保留下的edge会减少很多。
图片.png
下方的工具栏中选择Analysis可以看到这个nework的基本信息，也包括在几个database里的GO富集结果，比如KEGG Pathways和Uniprot Keywords。
图片.png
图片.png
如果想详细查看具体富集信息，下方可以直接下载保存，格式是tab分割的，阅读友好。
图片.png
在Export里可以看到每个相互作用中两个node的annotation，各种格式都可以下载。一般情况下对于不是特别常见的模式植物用处不是很大。
图片.png
因为STRING有Cluster database，还可以对你的结果进行简单的cluster，cluster的结果可以使用Excel处理或Cytospce来进行可视化。
图片.png
除了网站，STRING还有R package——stingDB，但是每次输出只能使用400个基因的network，而且引用的database是2016年更新v10，网站已经更新的v11，所以不太建议使用。在Cytospace里有stringApp，下载和使用都很方便，也可以同样使用GO分析，同时可以联合MCODE app进行cluster分析，这个教程搜索关键词就可以轻松找到，几乎傻瓜式使用，就是速度不是很快。
