使用eggnog-mapper注释基因组

caigui

2023-04-10

基因组

462

参考资料：
github https://github.com/eggnogdb/eggnog-mapper
原始文献 https://doi.org/10.1093/molbev/msx148
一个中文教程 https://www.jianshu.com/p/9d38cb9f1d02

1. 实验室服务器下操作

1.1 默认参数运行

conda create -n eggnog
conda activate eggnog
conda install eggnog-mapper
mkdir /home/caigui/miniconda3/envs/eggnog/lib/python3.11/site-packages/data 
download_eggnog_data.py # 一路y下去，解压后有48G数据库
# 运行，这里使用每个基因里最长的那个蛋白作为代表，并去除转录本序号以简化基因名称
nohup emapper.py -i xiaocuiyun_protein_represent_rename.fa -o xiaocuiyun_20230404_eggnog --tax_scope Viridiplantae --cpu 40 --dbmem &
# 好像不到十分钟就注释完了？？？
# 结果文件
less -S xiaocuiyun_20230404_eggnog.emapper.annotations

1.2 与interproscan结果进行比较

选一个基因SK03G20400，对比一下eggnog-mapper和interproscan的注释结果

################################################################################################################################################
# eggnog
SK03G20400	88036.EFJ33044	1.36e-67	222.0	28MZX@1|root,2QVBF@2759|Eukaryota,37RA9@33090|Viridiplantae,3GF59@35493|Streptophyta	33090|Viridiplantae	K	wuschel-related homeobox	-	GO:0000003,GO:0001067,GO:0003006,GO:0003674,GO:0003676,GO:0003677,GO:0003700,GO:0005488,GO:0006355,GO:0007275,GO:0008150,GO:0009653,GO:0009790,GO:0009791,GO:0009793,GO:0009888,GO:0009889,GO:0009933,GO:0009987,GO:0010014,GO:0010016,GO:0010072,GO:0010087,GO:0010154,GO:0010468,GO:0010556,GO:0019219,GO:0019222,GO:0022414,GO:0031323,GO:0031326,GO:0032501,GO:0032502,GO:0043565,GO:0044212,GO:0048316,GO:0048367,GO:0048507,GO:0048508,GO:0048532,GO:0048608,GO:0048731,GO:0048856,GO:0050789,GO:0050794,GO:0051171,GO:0051252,GO:0051301,GO:0060255,GO:0061458,GO:0065007,GO:0080090,GO:0090421,GO:0097159,GO:0140110,GO:1901363,GO:1903506,GO:2000112,GO:2001141	-	ko:K18490	-	-	-	-	ko00000,ko03000	-	-	-	Homeobox
# eggnog的go term
GO:0000003,GO:0001067,GO:0003006,GO:0003674,GO:0003676,GO:0003677,GO:0003700,GO:0005488,GO:0006355,GO:0007275,GO:0008150,GO:0009653,GO:0009790,GO:0009791,GO:0009793,GO:0009888,GO:0009889,GO:0009933,GO:0009987,GO:0010014,GO:0010016,GO:0010072,GO:0010087,GO:0010154,GO:0010468,GO:0010556,GO:0019219,GO:0019222,GO:0022414,GO:0031323,GO:0031326,GO:0032501,GO:0032502,GO:0043565,GO:0044212,GO:0048316,GO:0048367,GO:0048507,GO:0048508,GO:0048532,GO:0048608,GO:0048731,GO:0048856,GO:0050789,GO:0050794,GO:0051171,GO:0051252,GO:0051301,GO:0060255,GO:0061458,GO:0065007,GO:0080090,GO:0090421,GO:0097159,GO:0140110,GO:1901363,GO:1903506,GO:2000112,GO:2001141
################################################################################################################################################

################################################################################################################################################
# interproscan
SK03G20400	05c99823c90e86309c61782e8ea36ed9	323	CDD	cd00086	homeodomain	102	164	1.04591E-11	T	07-01-2022	IPR001356	Homeobox domain	GO:0003677
SK03G20400	05c99823c90e86309c61782e8ea36ed9	323	Gene3D	G3DSA:1.10.10.60	-	101	168	3.1E-17	T	07-01-2022	-	-
SK03G20400	05c99823c90e86309c61782e8ea36ed9	323	MobiDBLite	mobidb-lite	consensus disorder prediction	163	212	-	T	07-01-2022	-	-
SK03G20400	05c99823c90e86309c61782e8ea36ed9	323	ProSiteProfiles	PS50071	'Homeobox' domain profile.	98	163	15.256396	T	07-01-2022	IPR001356	Homeobox domain	GO:0003677
SK03G20400	05c99823c90e86309c61782e8ea36ed9	323	PANTHER	PTHR46777	WUSCHEL-RELATED HOMEOBOX 13	36	217	9.1E-61	T	07-01-2022	IPR044559	WUSCHEL-related homeobox 13-like	GO:0003700
SK03G20400	05c99823c90e86309c61782e8ea36ed9	323	Pfam	PF00046	Homeodomain	102	162	3.7E-17	T	07-01-2022	IPR001356	Homeobox domain	GO:0003677
SK03G20400	05c99823c90e86309c61782e8ea36ed9	323	SUPERFAMILY	SSF46689	Homeodomain-like	94	166	4.28E-17	T	07-01-2022	IPR009057	Homeobox-like domain superfamily	-
SK03G20400	05c99823c90e86309c61782e8ea36ed9	323	SMART	SM00389	HOX_1	100	167	5.6E-13	T	07-01-2022	IPR001356	Homeobox domain	GO:0003677
SK03G20400	05c99823c90e86309c61782e8ea36ed9	323	MobiDBLite	mobidb-lite	consensus disorder prediction	168	197	-	T	07-01-2022	-	-
SK03G20400	05c99823c90e86309c61782e8ea36ed9	323	Coils	Coil	Coil	60	80	-	T	07-01-2022	-	-
# interproscan的go term
GO:0003677,GO:0003700
################################################################################################################################################

怎么说呢，interproscan的go，eggnog那边都有。而且eggnog那边的注释非常简洁，就只有一句“wuschel-related homeobox”

1.3 与拟南芥同源基因注释比较

看一下这个基因的拟南芥同源基因在TAIR上的注释

# AtWOX13在org.At.tair.db的GO注释
"GO:0000976" "GO:0003700" "GO:0005634" "GO:0006355" "GO:0006631" "GO:0006996" "GO:0009888" "GO:0010025" "GO:0010087" "GO:0010629" "GO:0018193" "GO:0019748" "GO:0032787" "GO:0042546" "GO:0044255" "GO:0044550"

AtWOX13有以上16个go term，和eggnog的59个有4个交集。为什么eggnog可以自信地给出这么多GO注释？这就需要看eggnog的原始文献了。原始文献的摘要就说了，平均每个蛋白多41个term，确实。摘要还说它各个方面都比我之前用的interProScan好。

原始文献 https://doi.org/10.1093/molbev/msx148

1.4 hmmer代替默认的dimand

用hmmer代替默认的dimand进行搜索试试。hmmer速度很慢，dimand速度很快，看看注释结果有什么区别。

# 用hmmer代替默认的dimand试一试
nohup emapper.py -i xiaocuiyun_protein_represent_rename.fa -o xiaocuiyun_20230404_eggnog --tax_scope Viridiplantae --cpu 40 --dbmem -m hmmer &
# 报错
# emapper.py: error: HMMER mode requires a target database (-d, --database).
nohup emapper.py -i xiaocuiyun_protein_represent_rename.fa -o xiaocuiyun_20230404_eggnog --tax_scope Viridiplantae --cpu 40 --dbmem -m hmmer -d euk &
# 还是报错

# 解决办法：需要下载数据库
# http://eggnog5.embl.de/#/app/downloads 查看taxID，下面的33090是Viridiplantae的ID
# 下载数据库
download_eggnog_data.py -H -d 33090
# 运行注释程序
nohup emapper.py -i xiaocuiyun_protein_represent_rename.fa -o xiaocuiyun_20230404_eggnog --tax_scope Viridiplantae --cpu 40 --dbmem -m hmmer -d 33090 &
################################################################################################################################################
SK03G20400	88036.EFJ33044	2.2e-85	280.0	28MZX@1|root,2QVBF@2759|Eukaryota,37RA9@33090|Viridiplantae,3GF59@35493|Streptophyta	33090|Viridiplantae	K	wuschel-related homeobox	-	GO:0000003,GO:0001067,GO:0003006,GO:0003674,GO:0003676,GO:0003677,GO:0003700,GO:0005488,GO:0006355,GO:0007275,GO:0008150,GO:0009653,GO:0009790,GO:0009791,GO:0009793,GO:0009888,GO:0009889,GO:0009933,GO:0009987,GO:0010014,GO:0010016,GO:0010072,GO:0010087,GO:0010154,GO:0010468,GO:0010556,GO:0019219,GO:0019222,GO:0022414,GO:0031323,GO:0031326,GO:0032501,GO:0032502,GO:0043565,GO:0044212,GO:0048316,GO:0048367,GO:0048507,GO:0048508,GO:0048532,GO:0048608,GO:0048731,GO:0048856,GO:0050789,GO:0050794,GO:0051171,GO:0051252,GO:0051301,GO:0060255,GO:0061458,GO:0065007,GO:0080090,GO:0090421,GO:0097159,GO:0140110,GO:1901363,GO:1903506,GO:2000112,GO:2001141	-	ko:K18490	-	-	-	-	ko00000,ko03000	-	-	-	Homeobox
################################################################################################################################################

# 将--tax_scope Viridiplantae改为--tax_scope 33090，试试有什么不同
nohup emapper.py -i xiaocuiyun_protein_represent_rename.fa -o xiaocuiyun_20230404_eggnog --tax_scope 33090 --cpu 40 --dbmem -m hmmer -d 33090 &
# 运行正常，估计和上面没什么不同
# 经过比对，结果确实一模一样

# 试试把-d也换成Viridiplantae
nohup emapper.py -i xiaocuiyun_protein_represent_rename.fa -o xiaocuiyun_20230404_eggnog --tax_scope Viridiplantae --cpu 40 --dbmem -m hmmer -d Viridiplantae &
# 报错，看来-d只能选33090，但是tax_scope可以选33090或者Viridiplantae

1.5 比较一下web版本的注释结果

# 试试官方的web版本的注释情况
# http://eggnog-mapper.embl.de/
# 显示注释完成了，但是结果文件下载全部失败，显示为空
# 过了几个小时再去看，发现文件又都有了，看来只是反应有点慢？
# 同一个基因的结果文件
################################################################################################################################################
SK03G20400	88036.EFJ33044	1.21e-68	224.0	28MZX@1|root,2QVBF@2759|Eukaryota,37RA9@33090|Viridiplantae,3GF59@35493|Streptophyta	35493|Streptophyta	K	wuschel-related homeobox	-	GO:0000003,GO:0001067,GO:0003006,GO:0003674,GO:0003676,GO:0003677,GO:0003700,GO:0005488,GO:0006355,GO:0007275,GO:0008150,GO:0009653,GO:0009790,GO:0009791,GO:0009793,GO:0009888,GO:0009889,GO:0009933,GO:0009987,GO:0010014,GO:0010016,GO:0010072,GO:0010087,GO:0010154,GO:0010468,GO:0010556,GO:0019219,GO:0019222,GO:0022414,GO:0031323,GO:0031326,GO:0032501,GO:0032502,GO:0043565,GO:0044212,GO:0048316,GO:0048367,GO:0048507,GO:0048508,GO:0048532,GO:0048608,GO:0048731,GO:0048856,GO:0050789,GO:0050794,GO:0051171,GO:0051252,GO:0051301,GO:0060255,GO:0061458,GO:0065007,GO:0080090,GO:0090421,GO:0097159,GO:0140110,GO:1901363,GO:1903506,GO:2000112,GO:2001141	-	ko:K18490     	-	-	-	-	ko00000,ko03000	-	-	-	Homeobox
################################################################################################################################################
可以看到这里自动识别为了Streptophyta，而不是前面用的Viridiplantae。
看了看，上面那个查看taxID的链接内，也有Streptophyta的选项

2. 最后画个维恩图看一下本地的dimand，本地的hmmer，和web版本默认参数下的注释区别

结果就是一模一样

3. 写一个小脚本将GO term提取出来

# 输出所有基因到一个文件
grep -o "SK[0-9Z]\{1,2\}G[0-9]\{5\}" xiaocuiyun_20230404_eggnog.emapper.annotations > sk.eggnog.id

# 跑个循环，找出每个基因的goterm，并按照"gene，go_term"的格式输出
for i in `cat sk.eggnog.id`
do
go_list=($(grep $i xiaocuiyun_20230404_eggnog.emapper.annotations | grep -o "GO:[0-9]\{7\}"))

for go in "${go_list[@]}"
do
echo "$i,$go" >> sk.eggnog.emapper.annot.GO.trem.20230616.csv
done

done

4. 总结

总的来说，用默认参数就可以了。
运行示例nohup emapper.py -i xiaocuiyun_protein_represent_rename.fa -o xiaocuiyun_20230404_eggnog --tax_scope Viridiplantae --cpu 40 --dbmem &
其中--tax_scope这个参数可以在这里查询
ps：--tax_scope这个参数好像不用也行哈