之前无论是做blastn或者是blastp都是仅仅在blast结果中显示了种属信息,没有添加详细的信息,添加种属信息很简单,只要结果里生成XML格式,之后转化成xls格式就OK:
blastn -query 16S_seq/996-SLM.new.fa -db $dbNT -out 996-SLM.new.fa.nt.xml -outfmt 5 -evalue 1e-05 -max_target_seqs 1 -num_threads 28 #以blastn 为例
但是如果想要添加详细的物种信息,例如界门纲目科属种,就需要利用到NCBI taxonomy数据库。
还是以blastn,比对NT库为例。
首先,我们需要找到accession号对应的ncbi里面的taxonomy id号,也就是taxid。
打开NCBI taxonomy地址,ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/,在accession2taxid下找到nucl_gb.accession2taxid.gz,这个文件记录了所有核酸序列(非WGS或TSA)asscession对应的taxid。这里要说明一下,在18年甚至更久以前,NCBI是通过gi号与taxid对应起来,而如今,gi号已经慢慢淡出视野。
有了blastn结果,注意输出-outfmt选择6,我们可以使用grep命令抓取对应的taxid号:
grep '\bNR_043136.1\b' nucl_gb.accession2taxid | awk '{print $3}' #这里为了保证结果唯一,可添加个head -1输出一条结果
有了taxid号后,借助一个很好用的TaxonKit工具,来获取详细信息。该工具之前记录过,可在“如何从NT/NR数据库中提取子库”中找到。
与上文抓取到的taxid号,整合taxonkit命令生成详细的物种信息:
grep '\bNR_043136.1\b' nucl_gb.accession2taxid | awk '{print $3}' | taxonkit lineage | taxonkit reformat -F | cut -f 3
最终得到:
Bacteria;Proteobacteria;Alphaproteobacteria;Sphingomonadales;Erythrobacteraceae;Qipengyuania;Qipengyuania vulgaris
该结果可以自行写个perl或者python添加到blastn结果中。
taxonkit工具是用GO语言写的,很好用,还会自动获取更新,作者真大神。
网友评论