De novo 基因组GO 注释小窗口

关于GO 注释的心得体会

目前对于GO功能注释的思路有 以下常见的三种:

1、BLAST+InterProScan => BLAST2GO

可以先通过使用 blastp 进行nr注释,那么问题就来了,nr 数据库截至2018年12月已经有120G左右的数据量了,根据我的个人使用经验,对25000条蛋白进行32线程的并行化比对,差不多要半个月左右(心中一万匹cnm飘过)。很显然,在计算资源有限的条件下这是不现实的,那么大概就会出现两种思路:

(1)缩小数据库:一种是使用Swiss-Prot等小的数据库进行注释;另一种就是根据基于taxid构建 nr子数据库了,具体的构建可见 基于taxid构建Blast database_bioinfomatics2medicine_新浪博客 (2)能不能使用一个更高效的比对工具呢,答案是肯定的 那就是 GitHub – bbuchfink/diamond: Accelerated BLAST compatible local sequence aligner.

其主页介绍如下:

DIAMOND is a sequence aligner for protein and translated DNA searches, designed for high performance analysis of big sequence data. The key features are:

Pairwise alignment of proteins and translated DNA at 500x-20,000x speed of BLAST.

Frameshift alignments for long read analysis.

Low resource requirements and suitable for running on standard desktops or laptops.

Various output formats, including BLAST pairwise, tabular and XML, as well as taxonomic classification.

对的,它就是那么快,相对于 blastp 的半个月,在 ‘–more-sensitive’ 模式下也仅需要4-5个小时。

在nr注释完之后,就到使用 Blast2GO mapping GO id 了。Blast2GO 有 PRO (学术版一年的费用为850€,就看你的大老板可不可以支持了;当然你也可以申请个7天试用体验一下) 和 Basic (免费申请,期限好像为一年有效期;其对于 PRO 来说自然少了很多功能,但最尴尬的是 mapping 是真的有点慢,我的25000条蛋白 mapping 了接近两天的时间,当然也会有人为了加快速度split蛋白文件,多个PC进行跑 )两个版本,但不管怎么说,对于资金缺乏的人来说,使用 Basic 进行mapping 还是有点慢,那么不如就来个本地化吧,对于本地化网上也有很多的教程了:

blast2go本地化2017教程:blast2go本地化2017教程 – 生信技能树

非root权限的blast2go的安装和使用(二)· blast2go的数据和软件准备及使用:https://mp.weixin.qq.com/s?__biz=MzI1MjU5MjMzNA==&mid=2247486358&idx=1&sn=4b095e9a200a419079d947d930c38abd&chksm=e9e02237de97ab21f4227ee76eca0aa94433fe6cc374e030093eb3acc2094a546110763bb411&mpshare=1&scene=23&srcid=0101kMbzD2el0RPASBiC9pFi#rd

陈连福 centos6.9 代码:

“# installing Blast2go Databases

mkdir /opt/biosoft/blast2go

cd /opt/biosoft/blast2go

wget http://archive.geneontology.org/latest-full/go_monthly-assocdb-data.gz

wget ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/gene_info.gz ### ascp -k 1 -T -l 200M -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh –host=ftp-private.ncbi.nlm.nih.gov –user=anonftp –mode=recv /gene/DATA/gene_info.gz 。/

wget ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/gene2accession.gz  #### ascp -k 1 -T -l 200M -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh –host=ftp-private.ncbi.nlm.nih.gov –user=anonftp –mode=recv /gene/DATA/gene2accession.gz ./

wget ftp://ftp.pir.georgetown.edu/databases/idmapping/idmapping.tb.gz

## 对于不能使用aspera加速下载的文件可以试试 lftp -e ‘pget -n NUM -c url; exit’实行多线程下载

gzip -dv go_monthly-assocdb-data.gz

gzip -dv gene_info.gz

gzip -dv gene2accession.gz

gzip -dv idmapping.tb.gz

tar zxf ~/software/local_b2g_db.tar.gz

mv local_b2g_db/* ./  && rm -rf local_b2g_db/

perl -p -i -e ‘s/go_201512-assocdb-data/go_2017-assocdb-data/’ install_blast2goDB.sh

./install_blast2goDB.sh

#一个自己的例子 #

##下载NR 数据库到自己的Blast+ db 中

wget ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/2.2.28/ncbi-blast-2.2.28+-x64-linux.tar.gz

tar -zxf ncbi-blast-2.2.28+-x64-linux.tar.gz -C /opt/biosoft/

/opt/biosoft/ncbi-blast-2.2.28+/bin/blastdbcmd -db nr -entry all -out nr.faa

### 值得注意的是这里我用的blastdbcmd是ncbi-blast-2.2.28+的,输出的 nr.faa 文件中序列名是 “>gi|66816243|ref|XP_642131.1| hypothetical protein DDB_G0277827” ,但如果用的是 ncbi-blast-2.6.0+的 blastdbcmd,输出没有GI号,原因可能是:As of September 2016, the integer sequence identifiers known as “GIs” will no longer be included in the GenBank, GenPept, and FASTA formats supported by NCBI for the display of sequence records.In addition, the FASTA format will no longer include the database source abbreviation. Please refer to the NCBI News Announcement posting for more detail.

wget https://github.com/bbuchfink/diamond/releases/download/v0.9.24/diamond-linux64.tar.gz

### 同样注意版本,因为diamond输出的XML与 b2g4pipe_v2.5 不兼容而出现 “Annotation of 0 seqs with 0 annots finished. Now searching for orfan IPRs…” 这个issues的讨论,你可以在到下面的链接去看。

tar -zxf diamond-linux64.tar.gz

mv diamond ~/bin/

cp ../interpro/proteins.fasta ./

ascp -k 1 -T -l 200M -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh –host=ftp-private.ncbi.nlm.nih.gov –user=anonftp –mode=recv /pub/taxonomy/accession2taxid/prot.accession2taxid.gz ./prot.accession2taxid.gz

ascp -k 1 -T -l 200M -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh –host=ftp-private.ncbi.nlm.nih.gov –user=anonftp –mode=recv /pub/taxonomy/taxdmp.zip ./

unzip taxdmp.zip

/opt/biosoft/ncbi-blast-2.2.28+/bin/blastdbcmd -db nr -entry all -out nr.faa

diamond makedb –in nr.faa -d nr –taxonmap prot.accession2taxid.gz –taxonnodes nodes.dmp #### DIAMOND v0.9.24

diamond blastp –query proteins.fasta –more-sensitive –db nr –evalue 1e-5 –salltitles –threads 64 –outfmt 5 –out nr.xml

perl -p -e ‘s/diamond 0.9.24/BLASTP 2.2.26/’ nr.xml > nr_new.xml

perl -p -i -e ‘s/^Dbacces.dbname=.*/Dbacces.dbname=b2gdb/’ b2gPipe.properties

perl -p -i -e ‘s/^Dbacces.dbhost=.*/Dbacces.dbhost=127.0.0.1/’ b2gPipe.properties

java -cp *:ext/*: es.blast2go.prog.B2GAnnotPipe -in ../nr_new.xml -out go -annot -dat -annex

# 至于整合InterProScan的结果可以使用下面的脚本,当然也可以在java -cp *:ext/*: es.blast2go.prog.B2GAnnotPipe -ips 来实现,或者将b2g4pipe的输出文件加载到 Blast2GO Basic 中进行后续可视化分析。

merge_interpro_to_go.pl b2g4pipe/go.annot ../interpro/interpro.tsv > go.annot

# Anyone use diamond output xml file  for b2gpipe? · Issue #79 · bbuchfink/diamond · GitHub

# Problem in uploading the Diamond blastx results into Blast2GO · Issue #159 · bbuchfink/diamond · GitHub

# diamond to blast2go erro  · Issue #165 · bbuchfink/diamond · GitHub

2、InterProScan

对于 InterProScan 可以本地化,也可以用 interProScan5.pl 将序列发送到官方服务器进行注释,资源允许,还是推荐本地化,网络版感觉还是有点慢~一个本地化教程:

一个用interproscan做基因注释的简易教程

3、eggNOG-Mapper (包括网页版和本地化)

(1)关于这个注释方法网页版的推荐自己去摸索吧(就是和NCBI提交类似)=> http://eggnogdb.embl.de/#/app/emapper

(2)本地化推荐两个教程吧:

序列功能注释神器:eggNOG-mapper,KEGG/COG/KOG/GO/BiGG 一网打尽 : 序列功能注释神器:eggNOG-mapper,KEGG/COG/KOG/GO/BiGG 一网打尽 « Biostack.org

应该是最好的eggnog-mapper功能注释教程:应该是最好的eggnog-mapper功能注释教程 – 简书(生信媛)

关于eggNOG-Mapper 和 BLAST2GO的小心得:

这些方法我都用过,eggNOG-Mapper 是很多人推荐的,我当然也喜欢网页版的,不耗资源,数度快。相比其他的都快。eggNOG-Mapper自己的文章评测认为其准确性高于BLAST2GO,原理上我也这么认为。但有可能是我物种的特异性还是eggNOG-Mapper目前数据库还不够大,自己注释出的蛋白只有6000左右(25000),而BLAST2GO注释有12000多,所以两者目前我觉得都有优点及不足,关键还是看自己的选择。

前方高能:

对于准确性的个人思考,基于比对,也就是基因相似性的注释方法(blastp,diamond)可能会取决于数据库的大小,因为在blast2go中有涉及到相似度cut这一步,如果数据库不够大,那么相似度有的可能只为50(在cut值之上),所以只使用Swiss-Prot等小数据库注释时可能需要提高cut值;当然也会听到直系同源和旁系同源的说法,所以在新版Blast2go中增加了filter GO with Taxid (可能是看到eggnog的文章慌了)这一新模块,减少旁系同源的注释结果。那么如果要在老版的Blast2go中运用这个模块,或许可以换个思路,在比对时,通过Taxid来过滤比对结果。

eggnog的注释,存在两种模式diamond和hummer,前者官方推荐为存在接近物种,也就是在直系同源注释中有更好的表现,hummer则反之. 同样的,只要是比对,就会有cut相似度这一说,数据库的大小及深度(涉及的物种),就会影响结果。其一个特点在于直系同源注释。这点Blast2go中出了filter with Taxid 来补充。

居于以上的讨论,我个人觉得目前还是blast+interpro的准确性会更高一些,当然我觉得eggnog如果持续开发,以及整合diamond和hummer的结果将会成为新的宠儿。

小提示:

首先,第一次写东西,难免会有粗心的地方,望同行能批评指正;对于行文格式,以后会努力。其次,我个人不太会写程序,所以喜欢收集别人的程序,可以在 Zhanmengtao (zhanmengtao) · GitHub 找到一些可能对你来说感兴趣的程序,包括上面提到的一些。

参考文献

http://blog.sina.com.cn/s/blog_16152d7d70102xnc3.html

https://github.com/bbuchfink/diamond

https://mp.weixin.qq.com/s?__biz=MzI1MjU5MjMzNA==&mid=2247486358&idx=1&sn=4b095e9a200a419079d947d930c38abd&chksm=e9e02237de97ab21f4227ee76eca0aa94433fe6cc374e030093eb3acc2094a546110763bb411&mpshare=1&scene=23&srcid=0101kMbzD2el0RPASBiC9pFi#rd

http://www.jinciwei.cn/a342744.html

https://github.com/bbuchfink/diamond/issues/79

https://github.com/bbuchfink/diamond/issues/159

https://github.com/bbuchfink/diamond/issues/165

序列功能注释神器:eggNOG-mapper,KEGG/COG/KOG/GO/BiGG 一网打尽

https://www.jianshu.com/p/e646c0fa6443

NGS 生物信息学分析 V6.0 陈连福 郑越

https://www.jianshu.com/p/6296385adf21

    原文作者:詹梦涛
    原文地址: https://www.jianshu.com/p/2f8ec50020ec
    本文转自网络文章,转载此文章仅为分享知识,如有侵权,请联系博主进行删除。
点赞