宏基因组测序分析(十)功能注释
eggnog注释
使用在线版 eggnog-mapper 进行 eggnog 数据库注释。
eggnong数据库全称为evolutionary genealogy of genes: Non-supervised Orthologous Groups,是对NCBI的COG数据库的拓展,提供不同分类水平蛋白的直系同源组(Orthologous Groups,OG)。它扩展了COG数据库的分类方法,并且提供进行eggnog注释的工具:eggNOG-mapper
在线版 eggnog-mapper 单次输入序列条数不能超 10 万条,可以对氨基酸序列进行切分。
seqkit split --by-size 100000 --out-dir split unigene_pep.fasta
数据准备好以后,登录在线网站上传蛋白序列, 在线进行蛋白注释,完成后下载注释结果 *emapper.annotations,并将其合并。
uniprot 注释
是一个提供蛋白质序列和注释数据的综合性数据库。对于序列数据库,我们会使用blastp/diamond这类软件进行比对。
由于 eggnog 数据库对细菌的 GO 注释信息包含很少,这里我们使用uniref90 数据库对氨基酸序列进行比对注释,并基于 idmapping 信息提取 GO 注释。
# 下载uniref90数据库
wget https://ftp.uniprot.org/pub/databases/uniprot/uniref/uniref90/uniref90.fasta.gz
# 下载id mapping文件 11G
wget https://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/idmapping/idmapping_selected.tab.gz
# 构建diamond database
diamond makedb --in ./uniref90.fasta.gz --db uniref90
# 进行diamond比对
diamond blastp \
--db ./uniref90 \ # db名称
--query unigene_pep.fasta \ # 输入,氨基酸fasta文件
--out unigene_pep.uniref90.m6 \ # 输出,表格格式
--threads 8 \ # 线程
--outfmt 6 \ # 输出blast表格格式结果
--max-target-seqs 5 \
--evalue 1e-5
# 基于idmapping 提取GO注释信息
perl uniref90_idmapping.pl \
unigene_pep.uniref90.m6 \ # 比对结果
./idmapping_selected.tab.gz \ # idmapping文件
> unigene_pep.uniref90.GOanno # 输出GO注释信息
# 构建orgdb,并统计GO 注释
Rscript GOmapperx.R unigene_pep.uniref90.GOanno
抗生素耐药基因注释
目前常用的耐药数据库:
-
ARDB数据库:
最早的耐药数据库,其核心架构包含完整的耐药基因序列以及,除此之外,还包括部分抗生素靶位点序列及其相关信息。目前ARDB已不再维护,相关信息已经被整合进CARD数据库,因此,目前常用CARD进行耐药性研究。
-
CARD数据库
CARD是一个基于志愿者贡献数据的耐药性研究共享平台,包括ARDB数据库的所有耐药基因信息,且更新很及时,从而保证了数据的时效性。CARD数据库以ARO为分类单位,在一个term中关联抗生素、作用靶位点、作用机制、序列变异等信息。除此之外,还包括一些毒力基因的信息、相关的可移动元件信息等。分析模式包括BLAST和RGI模式,除了可以对已知耐药基因进行注释外,还可以通过RGI预测潜在的耐药基因。RGI提供了3种预测标准,即Perfect、Strict和Loose;通过选择同源比对的判定标准,可以得到不同可信度和数量的潜在耐药基因,有助于发现新的耐药基因。
RGI 可以网页版运行,也可以安装在本地服务器,在线地址https://card.mcmaster.ca/analyze/rgi
# 下载CARD数据库参考数据
wget https://card.mcmaster.ca/latest/data
tar -xvf data ./card.json
# 准备输入文件
unigene_pep.fasta
# 加载数据库
rgi load \
-i ./card.json \ # 输入
--local # 数据库存放在当前目录
# 运行主程序
rgi main \
--input_sequence unigene_pep.fasta \ # 输入,氨基酸序列
--output_file rgi_out \ # 输出文件前缀
--input_type protein \ # 输入数据类型contig|protein
--alignment_tool DIAMOND \ # 比对软件
--num_threads 8 \ # 线程数
--local \ # 数据库在当前目录
--clean \ # 删除临时文件
--include_loose # 保留宽松比对结果
碳水化合物酶注释
碳水化合物酶注释的常用数据库为CAZy 数据库。
CAZy预测工具dbCAN2
这里使用 dbCAN2 提供的 run_dbCAN 本地版进行碳水化合物酶注释,dbCNA 也可以使用在线版,网址如下:dbCAN3 server (unl.edu)
# 使用diamond比对方法
run_dbcan \
--out_dir cazy_diamond \ # 输出目录
--db_dir 你的路径 \ # 数据库路径
--tools diamond \ # 使用diamond比对进行注释
--dia_cpu 8 \ # 线程数
unigene_pep.fasta \ # 输入,蛋白fasta序列
protein # 输入数据格式
# 使用hmmer结构域搜索方法
run_dbcan \
--out_dir cazy_hmmer \
--db_dir /app/db \ # 数据库路径
--tools hmmer \ # 使用结构域搜索进行注释
--dia_cpu 8 \ # 线程数
--hmm_cpu 8 \
--eCAMI_jobs 5 \
unigene_pep.fasta \ # 输入,蛋白fasta序列
protein
主要输出结果: