【自用推文】NCBI blast物种注释2

最近公司2个转录组样本mapping 比例在10%，公司的NT数据库比对和结果物种注释脚本稀烂，并不能很好的看出物种比例。于是就有了之前的推文《NCBI blast物种注释》。今天终于有空，将其封为了脚本。

本文中各脚本均不公开，仅对脚本过程和特点进行描述，可以结合之前的推文写一下自己的过程。

示例

zcat SampleID.xml.gz |grep 'Iteration_query-def' -A 7 |grep 'Hit_accession' |sed -e 's/.*>\(\S\+\)<.*/\1/g' > SampleID_accession.txt

Rscript run_accession2taxid.r -i SampleID_accession.txt -o SampleID_accession2taxid.xls

cat SampleID_accession2taxid.xls |cut -f 3 |sed '1d' |sort |uniq > SampleID_uniq_taxid.xls

taxonkit lineage SampleID_uniq_taxid.xls  -j  20 |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' > SampleID_species.xls

Rscript run_accession2species.r --accession SampleID_accession.txt --spices SampleID_species.xls --accession2taxid SampleID_accession2taxid.xls -o SampleID_species_anno.xls

步骤

第一步：从NT数据库的Blast XML结果中提取每个reads的第一个比对结果条目的accession。

SampleID_accession.txt

通过定位搜索序列ID的位置来定位第一个比对结果，提取对应的accession。

$ zcat SampleID.xml.gz |grep 'Iteration_query-def' -A 7 |grep 'Hit_accession' |sed -e 's/.*>\(\S\+\)<.*/\1/g' > SampleID_accession.txt

$ head SampleID_accession.txt
XM_001166213
XM_008958708
NG_042093
XM_011725543
5A2Q_2
XM_012746077
NM_001271972
XR_930872
XM_011728545
NG_027821

第二步：从NT数据库的nucl_gb.accession2taxid.gz 中匹配，获取每条accession 对应的分类IDtaxid。

这一步比较吃内存资源，我分析时耗费内存大概在90G左右，而且是这4步注释中的限速步骤。

$ Rscript run_accession2taxid.r -i SampleID_accession.txt -o SampleID_accession2taxid.xls

$ head SampleID_accession2taxid.xls
accession	accession.version	taxid	gi
AB041222	AB041222.1	9606	446512455
AB046838	AB046838.1	9606	10047312
AB064199	AB064199.1	9606	21669604
AB064207	AB064207.1	9606	21669620
AB074982	AB074982.1	9606	19912273
AB208902	AB208902.1	9606	62087383
AB209021	AB209021.1	9606	62087621
AB209477	AB209477.4	9606	349500966
AB209585	AB209585.1	9606	62088749

第三步：利用taxonkit ，将taxid 注释到物种。

$ cat SampleID_accession2taxid.xls |cut -f 3 |sed '1d' |sort |uniq > SampleID_uniq_taxid.xls
$ taxonkit lineage SampleID_uniq_taxid.xls  -j  20 |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' > SampleID_species.xls

$ head SampleID_species.xls
Taxid	Kingdom	Phylum	Class	Order	Family	Genu	Species
10020	Eukaryota	Chordata	Mammalia	Rodentia	Heteromyidae	Dipodomys	Dipodomys ordii
10036	Eukaryota	Chordata	Mammalia	Rodentia	Cricetidae	Mesocricetus	Mesocricetus auratus
10090	Eukaryota	Chordata	Mammalia	Rodentia	Muridae	Mus	Mus musculus
10160	Eukaryota	Chordata	Mammalia	Rodentia	Octodontidae	Octodon	Octodon degus
10181	Eukaryota	Chordata	Mammalia	Rodentia	Bathyergidae	Heterocephalus	Heterocephalus glaber
112262	Eukaryota	Chordata	Mammalia	Artiodactyla	Bovidae	Ovis	Ovis canadensis
1155699	Eukaryota	Mollusca	Gastropoda	unclassified Gastropoda order	Potamididae	Cerithideopsis	Cerithideopsis scalariformis
1298539	unclassified other entries superkingdom	unclassified other entries phylum	unclassified other entries class	unclassified other entries order	unclassified other entries family	unclassified other entries genus	Cloning vector pDAL649
143302	Eukaryota	Chordata	Mammalia	Eulipotyphla	Talpidae	Condylura	Condylura cristata

第四步：合并上述各步骤结果，并作简单统计。

这一步主要就是数据的merge，没什么特殊的。

• SampleID_accession.txt 主要就是通过accession提供参与分析的reads数量，便于后续统计
• SampleID_accession2taxid.xls 是accession2taxid的结果。
• SampleID_species.xls 是taxid2species 的结果。

上述3个文件就能统计出比对到每个物种的reads数量和比例。

$ Rscript run_accession2species.r --accession SampleID_accession.txt --spices SampleID_species.xls --accession2taxid SampleID_accession2taxid.xls -o SampleID_species_anno.xls

$ head SampleID_species_anno.xls
Taxid	Number	Percent%	Kingdom	Phylum	Class	Order	Family	Genu	Species
9606	3635	61.205590166694726	Eukaryota	Chordata	Mammalia	Primates	Hominidae	Homo	Homo sapiens
61853	520	8.755682774877926	Eukaryota	Chordata	Mammalia	Primates	Hylobatidae	Nomascus	Nomascus leucogenys
9545	471	7.930628051860583	Eukaryota	Chordata	Mammalia	Primates	Cercopithecidae	Macaca	Macaca nemestrina
32630	414	6.970870516922041	unclassified other entries superkingdom	unclassified other entries phylum	unclassified other entries class	unclassified other entries order	unclassified other entries family	unclassified other entries genus	synthetic construct
9597	209	3.5191109614413203	Eukaryota	Chordata	Mammalia	Primates	Hominidae	Pan	Pan paniscus
9598	167	2.811921198855026	Eukaryota	Chordata	Mammalia	Primates	Hominidae	Pan	Pan troglodytes
9531	118	1.986866475837683	Eukaryota	Chordata	Mammalia	Primates	Cercopithecidae	Cercocebus	Cercocebus atys
336983	80	1.3470281192119886	Eukaryota	Chordata	Mammalia	Primates	Cercopithecidae	Colobus	Colobus angolensis
37293	66	1.1112981983498906	Eukaryota	Chordata	Mammalia	Primates	Aotidae	Aotus	Aotus nancymaae

• Number 列是比对到该物种的reads数。结果中针对该列进行了降序排列。

• Percent% 列是比对到该物种的reads数占注释到结果的物种的比例。并不是每个reads都能注释到结果，但是注释无结果的reads仅占极少比例，因此这个数值也是近似于比对到该物种的reads数占参与比对的reads的比例。

可能存在的问题

• 不用在笔记本上尝试，内存不够。我测试时的内存峰值大概在90G左右。
• blast步骤可以优化下，仅保留最佳匹配，但是这对新手可能有一点点难度。diamond能否轻松实现，我没测试过。
• 因为NT数据库与nucl_gb.accession2taxid.gz版本原因，结果中可能存在NA，也是正常的。具体处理可以参考taxonkit 的说明文档。