在非模式生物中,我们想要获取相应物种的GO或KEGG注释信息,或者描述信息的渠道很少。数据分布不集中,需要我们花费很长时间去收集,整合,才能得到相应的信息去做功能富集分析。下面以小麦为例,我们如何去收集功能富集前的准备工作。
一、功能注释
对一些未知功能的序列做初步的功能注释,通常的方法就是先下载一些功能数据库(NR、NT、uniprot、GO、KEGG),然后将序列进行比对,从而预测功能。但是这些数据库通常比较庞大,下载慢,处理起来也比较复杂。因此,我们一般采用在线做功能注释。有些网站例如eggNOG-MAPPER、PANNZER、g:Profiler等。我认为在线功能注释网站中比较快的是eggNOG-MAPPER,PANNZER专注做功能注释约有10年了,也是比较专业的,大家也可以多在几个网站中做,然后比较一下结果。取自己想要的,合适的结果。
下面就以eggNOG-MAPPER网站为例,做以下序列的功能注释:
示例数据取自小麦中国春v2.1版本的高可信蛋白序列。一般我们以蛋白序列作为功能注释的输入文件,示例文件名:Example_pep.fa
image.png
image.png
上面为填写后的网页,我一般采用默认参数,然后提交。在邮箱中会看到下面界面。
image.png
点击后会转到下面界面
image.png
开始后只需等待就可以啦
image.png
当运行完成之后会给你邮箱发送一个任务完成邮件。
任务完成啦,然后下载结果。
image.png
一般下载 out.emapper.annotations和out.emapper.annotations.xlsx两个结果就够了,其中out.emapper.annotations为后面处理数据的输入文件。
image.png
我们来看一下结果文件里面有什么内容
image.png
文件内容这么乱,想得到我们需要的输入文件,还需要进行进一步整理,下面就使用TBtools软件进行整理,确实很方便。使用GO&KEGG模块里面的eggNOG-mappper Helper功能进行处理。
image.png
输入文件为out.emapper文件
image.png
在目录下就会产生每个功能注释库所对应的蛋白质ID
image.png
out.emapper.annotations.GO文件里,有GO与蛋白质ID对应关系,这样我们就初步得到我们所需要的文件啦。
image.png
image.png
image.png
对于out.emapper.annotations.KEGG_Pathway文件中,我们可以用以下命令简单处理:
grep map out.emapper.annotations.KEGG_Pathway.txt > out.emapper.annotations.KEGG_map.tx
grep ko out.emapper.annotations.KEGG_Pathway.txt > out.emapper.annotations.KEGG_ko.txt
这样我们就能提取出只有map号和KO号的文件啦。
image.png
image.png
到此我们前景基因的功能注释就完成啦。
二、获取该物种所对应的GO、KEGG注释
获取物种与GO号所对应的关系的方法有以下几种方式:
1、通过GO官网进行下载http://current.geneontology.org/annotations/index.html或者http://current.geneontology.org/products/pages/downloads.html下载对应物种的注释信息。很可惜,官网里面并没有小麦的注释信息。
2、 从COA项目中下载ftp://ftp.ebi.ac.uk/pub/databases/GO/goa/,但是也没有小麦的注释信息。
3、从NCBI中下载数据ftp.ncbi.nih.gov
image.png
下载gene2go.gz文件,打开后根据Tax id进行筛选,小麦的Tax id是4565,很可惜也没有。
4、从Bioconductor 获取,虽然AnnotationHub包提供了常见物种的注释信息,却没有小麦的注释信息。
5、将小麦参考基因组蛋白质序列与uniprot蛋白质库进行blast,然后得到GO号。由于uniprot库下载速度慢,并且比较大,耗时耗力。我采用的是下面的方法得到背景基因。
6、将参考蛋白质序列在eggNOG-mapper在线网站进行比对,然后对比对结果进行处理,这样就能得到GO号以及KEGG的K num、map号、ko号等等。
image.png
三、GO、KEGG描述信息
GO描述信息可以在http://current.geneontology.org/ontology/go-basic.obo得到,是所有物种的描述信息,不影响做GO富集分析。进入下面页面后右击点击另存为,保存文件。然后提取文件中信息即可得到简化的描述信息。
使用下面python进行处理一下:
with open("go-basic.obo","r") as file:
lib={}
for line in file:
line=line.strip()
col_name=line.split(":")[0]
if col_name == "id":
id=line.split(" ",maxsplit=1)[1]
lib[id]=""
if col_name == "name":
name=line.split(" ",maxsplit=1)[1]
lib[id]=lib[id]+"@"+name
if col_name == "namespace":
namespace=line.split(" ",maxsplit=1)[1]
lib[id]=lib[id]+"@"+namespace
out=open("GO_basic_Description.txt","a+")
out.write("Class"+"\t"+"GO_IDs"+"\t"+"Description"+"\n")
for key in lib.keys():
go_id=key
go_name=lib[key].split("@")[1]
go_namespace=lib[key].split("@")[2]
if go_namespace == "molecular_function":
go_namespace="MF"
out.write(go_namespace+"\t"+go_id+"\t"+go_name+"\n")
if go_namespace == "biological_process":
go_namespace="BP"
out.write(go_namespace+"\t"+go_id+"\t"+go_name+"\n")
if go_namespace == "cellular_component":
go_namespace="CC"
out.write(go_namespace+"\t"+go_id+"\t"+go_name+"\n")
得到下面的文件 GO_basic_Description.txt:
KEGG描述信息同样也可以在KEGG官网中得到。参考网址https://www.kegg.jp/kegg/catalog/org_list.html。
image.png
然后找到小麦(Triticumaestivum)缩写为taes。点击raes进入小麦专栏
image.png
点击Brite hierarchy
image.png
然后在点击KEGG Orthology
image.png
然后在download栏中选择一种格式进行下载
image.png
json格式
image.png
KEG格式
image.png
发现上面两个格式不是我们想要的,还需要进行处理,显得有点麻烦。
下面推荐一种较为简便的方式,使用R包“KEGGREST”,这个包就比较简便,人性化。
下面代码就是提取相应的信息。只需要修改输出文件路径即可。
#获取KEGG数据库信息
#加载包
library(KEGGREST)
##查看KEGG数据库包含的数据
listDatabases()#
##获取pathway(所有物种)数据集中的数据
pathway<- keggList("pathway")
head(pathway)
#转换数据集,导出数据集
pathway_data<-as.data.frame(pathway)
write.table(pathway_data,"<path>/KEGG_pathway_allspacies_database.txt",row.names = T,col.names = F,sep = "\t")
##对单个(小麦)数据库进行物种的选择
taes_pathway <-keggList("pathway","taes")#taes是小麦的缩写
taes_pathway_data<-as.data.frame(taes_pathway)
write.table(taes_pathway_data,"<path>/KEGG_pathway_wheat_database.txt",row.names = T,col.names = F,sep = "\t")
##获取KO(所有基因)数据集中的数据
Ko<-keggList("ko")
Ko_data<-as.data.frame(Ko)
write.table(Ko_data,"<path>/KEGG_KO_allspacies_database.txt",row.names = T,col.names = F,sep = "\t")
##获取KO(小麦)数据集中的数据
taes_Ko<-keggList("ko","taes")
taes_Ko_data<-as.data.frame(taes_Ko)
write.table(taes_Ko_data,"<path>/KEGG_KO_wheat_database.txt",row.names = T,col.names = F,sep = "\t")
得到文件如下:
KEGG_pathway_allspacies_database.txt
KEGG_pathway_wheat_database.txt
KEGG_KO_allspacies_database.txt
KEGG_KO_wheat_database.txt
这个包不仅可以提取map号、ko号还可以提取其他的东西。下图就是所有可提的部分,只要修改以下keggList("pathway","taes")参数即可。
image.png
得到KEGG的描述信息后需要用一个python脚本简单处理一下,输入文件为KEGG_pathway_allspacies_database.txt、KEGG_KO_allspacies_database.txt
#将得到的KO、pathway文件进行处理,得到可以富集的描述信息文件
with open("KEGG_KO_allspacies_database.txt","r") as file:
out=open("KEGG_KO_allspacies_description.txt","a+")
for line in file:
line=line.strip()
if "[EC:" in line.split("\t",maxsplit=1)[1] :
desc=line.split("\t",maxsplit=1)[1].rsplit("[EC:",maxsplit=1)[0].split(";")[1].strip('"')
id=line.split("\t",maxsplit=1)[0].split(":")[1][:-1]
out.write(id+"\t"+desc+"\n")
else:
if ";" in line.split("\t",maxsplit=1)[1]:
out.write(line.split("\t",maxsplit=1)[0][:-1].split(":")[1]+"\t"+line.split("\t",maxsplit=1)[1].split(";")[1].strip('"')+"\n")
else:
out.write(line.split("\t",maxsplit=1)[0][:-1].split(":")[1]+"\t"+line.split("\t",maxsplit=1)[1].strip('"')+"\n")
with open("KEGG_pathway_allspacies_database.txt") as file1:
out1=open("KEGG_pathway__allspacies_description.txt","a+")
for line1 in file1:
line1=line1.strip()
map_id=line1.split("\t") [0].strip('"').split(":")[1]
descript=line1.split("\t") [1].strip('"')
out1.write(map_id+"\t"+descript+"\n")
输出文件为KEGG_KO_allspacies_description.txt、KEGG_pathway__allspacies_description.txt
image.png
image.png
这样我们GO、KEGG富集前准备文件就完成啦!
网友评论