列注释_宏基因组之物种注释(基于nr库)

点击上方「蓝字」关注我们

宏基因组物种注释的方法有很多,可以基于质控后的cleandata直接用metaphlan2，kraken2进行序列比对，还可以拿组装去冗余后的基因集翻译得到蛋白序列后拿去和Nr蛋白库比对,比对的工具有主流的blast和diamond,前几天刚好捣鼓了一下物种注释过程(基于nr库)，现在将整个流程和注意的细节记录一下，以供参考。

软件需求：blast，diamond，taxonkit(安装自行百度)

构建细菌子库

blast方法可能会准确点，但是它的速度简直让我怀疑人生，俩种软件的方法我都说下吧，因为我比对的主要是细菌，我首先想到是干脆按照网上的方法构建一个细菌的子库可能速度会更快点~ 说干就干

我最初想法是构建一张表，第一列是taxid，后面7列跟着门纲目科属种的名称，这样的话我们根据比对结果中的taxid，就直接可以得到物种各层级的信息，相信这也是大多数人的想法，就像下图这样：

这里我们就可以使用taxonkit这个好轮子去方便的达成这一目的，代码如下：

1.#输出细菌的taxid
 taxonkit list --ids 2 --indent "" > bacteria.taxid.txt 
#33090为植物 2为细菌 2157为Archaea；10239 为Viruses；Eukaryota为2759；Fungi 4751
#有时间的话可以将这几大类的id全都整合在一起，形成一张表
2.less bacteria.taxid.txt|taxonkit lineage |taxonkit reformat -f "{k}\t{p}\t{c}\t{o}\t{f}\t{g}\t{s}" -F |cut -f 1,3- | sed '1i\Taxid\tKingdom\tPhylum\tClass\tOrder\tFamily\tGenu\tSpecies' > Bacteria_taxid_Ano.txt

OK,现在已经得到了这张表，接下来我们就要去比对，得到每个基因的taxid

先来构建细菌子库吧

需要文件：

ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/accession2taxid/prot.accession2taxid.gz 蛋白acc号和taxid对应文件

方法参考自爪哥文档：https://bioinf.shenwei.me/taxonkit/tutorial/

#还是像上面那样先得到细菌的taxid
taxonkit list --ids 2 --indent "" > bac.taxid.txt
#得到bac.taxid.acc.txt文件  
zcat prot.accession2taxid.gz |csvtk -t grep -f taxid -P bac.taxid.txt |csvtk -t cut -f accession.version >bac.taxid.acc.txt 

Option 1 :
#直接构建(速度慢)
blastdbcmd -db /home/software/nr-2019-12-18/nr -entry_batch bac.taxid.acc.txt -out - | pigz -c > nr.$id.fa.gz

Option 2: 
#用blast配套工具(这里要注意blastdb_aliastool版本 低版本没有seqidlist) 这步结果就只有一个pal文件
blastdb_aliastool -seqidlist bacteria.taxid.acc.txt -db /home/software/nr-2019-12-18/nr -out nr_bac -title nr_bac  

Option 3:
#采用了分割并行的策略加上爪哥自己写的脚本，速度更加快，建议大家去看爪哥的源文档，这里不再赘述.

序列比对

blastp

我们拿上面Option2得到的文件去进行比对，命令行如下：

blastp -query NR100pro.fasta -db /home/pub_guest/db/nr/nr_bac -out D1_nr.out -outfmt "6 qseqid qgi qacc qaccver qlen sseqid qseq sseq evalue  score length pident  staxids sscinames salltitles " -num_threads 16 -evalue 1e-5 -num_alignments 5
 #记得结果加上 staxid 得到物种taxid信息

diamond

由于索引库不兼容，我们将blastcmd抽提出来的nr库，用diamond先构建索引库 要想得到taxid和种名信息，需要构建的时候额外增加俩个参数--taxonmap和--taxonnodes

1.第一个参数需要的是我们上述说的 蛋白acc号和taxid的对应文件prot.accession2taxid.gz

2.第二参数需要的是存储有taxonomy数据库的层级文件taxdmp.zip

注意wget的时候会出现连接超时情况，多等几次即可.

diamond makedb --in nr.fa -d Dimond_nr --taxonmap prot.accession2taxid.gz --taxonnodes nodes.dmp  
(#nodes.dmp文件是taxdmp.zip压缩包内的文件,这个地方直接跟这个文件即可，否则建库不完整)
#比对
diamond blastp -p 8 --db Dimond_nr  -q NR100pro.fasta --outfmt 6 qseqid sseqid pident length qlen mismatch gapopen qstart qend evalue bitscore staxids stitle salltitles --max-target-seqs 5 -e 1e-5 -o NR.out

补充一点：

我们可以先将有taixd也有种名信息的提取出来，过滤掉没有taxid和种名信息的基因，因为文件结果中可能会出现一种情况是该基因比对结果没有taxid，但是后面匹配到了种名信息，这时候我们再来用taxonkit这个软件根据种名去得到这些物种的taixid，这样我们的结果里会多一点，这时候还没有匹配到taxid的再将其过滤掉。

这里写个小python脚本贴出来供参考：输入文件为比对结果的三列信息:基因名，taxid，种名信息，结果就是筛选完后的信息了~~

#!usr/bin env python
import re
import sys
#根据nr结果1,12,13列信息提取taxid，种名
input_file = sys.argv[1]
fr = open(input_file,'r')

pattern = "\[[^\]]+"
regexp = re.compile(pattern)

for line in fr:
        line = line.strip()
        gene = line.split('\t')[0]
        taxid = line.split('\t')[1]
        if '[' not in line.split('\t')[2]:
                sth = ''
        else:
                sth = regexp.findall(line.split('\t')[2])[0]
        print(f"{gene}\t{taxid}\t{sth}")
fr.close()

好了，根据blast或者diamond得到结果中都会有taxid信息和种名，由于是比对的子库文件，速度也会有所提升，我们可以拿着taxid去匹配最后那张表的信息就能匹配种水平上面几个层级的taxonomy的信息，或者好好利用taxonkit这个好轮子也是可以得到的~~~

参考连接 

1.https://www.bioinfo-scrounger.com/archives/207/; 

2.https://bioinf.shenwei.me/taxonkit/tutorial/