【数据库】本地NR数据库如何按物种拆分？

2022-11-26 08:00:05

1.准备本地数据库文件
- 1.1 NR库下载
- 1.2 Taxonomy数据库下载
2.按物种拆分NR库

1.准备本地数据库文件

NR(Non-Redundant Protein Sequence Database)非冗余蛋白库，是所有GenBank+EMBL+DDBJ+PDB中的非冗余蛋白序列。Taxonomy物种分类数据库，包括大于7万余个物种的名字和系谱，这些物种都至少在遗传数据库中有一条核酸或蛋白序列。NR和Taxonomy数据库都是NCBI的子数据库，会提供比较全面的对应关系。在本地数据库按物种拆分的话，必须下载这两个数据库的文件。

1.1 NR库下载

ftp下载地址：ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA
NR数据库更新是相当频繁的，如果追求新，估计每个月甚至每周就重新下一次，但它又非常大，对于商业流程使用不可能更新得这么频繁，可以半年或一年更新一次。

1.2 Taxonomy数据库下载

ftp下载地址：ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/
同样，taxonomy更新也很快。

我们分库需要用到两个文件，一个是accession2taxid中的prot.accession2taxid文件：

该文件将accession和taxid关系对应起来（也有GI号，2016年以前大家用的是GI和taxid的对应文件，现在该文件已淘汰）。其格式为：

另一个是taxdump文件，里面包含了物种层级和物种名称等文件。解压后文件：

readme.txt文件中解释了每个文件的每一列信息（注意|是列间隔，而非列本身）：

*.dmp files are bcp-like dump from GenBank taxonomy database.

General information.
Field terminator is "\t|\t"
Row terminator is "\t|\n"

nodes.dmp file consists of taxonomy nodes. The description for each node includes the following
fields:
	tax_id					-- node id in GenBank taxonomy database （Taxonomy记录号）
 	parent tax_id				-- parent node id in GenBank taxonomy database （上一层分类级别的tax_id）
 	rank					-- rank of this node (superkingdom, kingdom, ...)  该tax_id所处的分类层级）
 	embl code				-- locus-name prefix; not unique
 	division id				-- see division.dmp file
 	inherited div flag  (1 or 0)		-- 1 if node inherits division from parent
 	genetic code id				-- see gencode.dmp file
 	inherited GC  flag  (1 or 0)		-- 1 if node inherits genetic code from parent
 	mitochondrial genetic code id		-- see gencode.dmp file
 	inherited MGC flag  (1 or 0)		-- 1 if node inherits mitochondrial gencode from parent
 	GenBank hidden flag (1 or 0)            -- 1 if name is suppressed in GenBank entry lineage
 	hidden subtree root flag (1 or 0)       -- 1 if this subtree has no sequence data yet
 	comments				-- free-text comments and citations

Taxonomy names file (names.dmp):
	tax_id					-- the id of node associated with this name （为taxonomy的记录号）
	name_txt				-- name itself  （即对应tax_id号的物种名称）
	unique name				-- the unique variant of this name if name not unique
	name class				-- (synonym, common name, ...)

Divisions file (division.dmp):
	division id				-- taxonomy database division id
	division cde				-- GenBank division code (three characters)
	division name				-- e.g. BCT, PLN, VRT, MAM, PRI...
	comments

Genetic codes file:
	genetic code id				-- GenBank genetic code id
	abbreviation				-- genetic code name abbreviation
	name					-- genetic code name
	cde					-- translation table for this genetic code
	starts					-- start codons for this genetic code

Deleted nodes file (delnodes.dmp):
	tax_id					-- deleted node id

Merged nodes file (merged.dmp):
	old_tax_id                              -- id of nodes which has been merged
	new_tax_id                              -- id of nodes which is result of merging

Citations file (citations.dmp):
	cit_id					-- the unique id of citation
	cit_key					-- citation key
	pubmed_id				-- unique id in PubMed database (0 if not in PubMed)
	medline_id				-- unique id in MedLine database (0 if not in MedLine)
	url					-- URL associated with citation
	text					-- any text (usually article name and authors).
						-- The following characters are escaped in this text by a backslash:
						-- newline (appear as "\n"),
						-- tab character ("\t"),
						-- double quotes ('\"'),
						-- backslash character ("\\").
	taxid_list				-- list of node ids separated by a single space

其中最关键的是names.dmp和nodes.dmp文件。names.dmp示例（共四列，重要的也就taxid和物种名的前两列信息）：

nodes.dmp示例（共13列，重要的也就taxid、上层级taxid、分类层级这前三列信息）：

为了让分类更简单，我们按taxonomy数据库本身分类的分法，即division.dmp文件，共12类物种。

2.按物种拆分NR库

2.1 第一步：获得Aceesson和分类物种的对应关系

根据以上的prot.accession2taxid.gz、nodes.dmp和division.dmp文件，可通过编写脚本来获得accession和以上12类物种的对应关系。脚本略，自己写。假设结果文件命名为acc2sp.xls，格式如下：

2.2 第二步：获得分类物种的序列

根据acc2sp.xls这个文件以及NR总库序列文件nr.gz，我们就可以获得各类物种的序列信息了。当然除了taxonomy数据库本身分的这12类，我们也可以将它们合并来自定义子库。比如这12类中没有动物，我们可以将Invertebrates.fa、 Mammals.fa、 Primates.fa、 Rodents.fa 和Vertebrates.fa合并为动物作为一类，也可以将"Bacteria"、"fungi"、"Viruses"、"Phages"和"Environmental.samples"等合并为微生物作为一类（这在宏组学注释中常用）。当然NR中也有这12类中没包含的序列，我们可将其归为unknown.fa（不同于Unassigned.fa，它是没有物种信息）。

脚本自己写，最后得到的是各个子数据库的fasta序列文件。

2.3 第三步：建库和比对

blast或diamond比对工具进行序列数据库建库，后面比对选择对应的字库就可。
blastall：

formatdb -p T  -i Plants.fa
blastall -i query.fa -d Plants.fa -o blastout.nr -p blastp -F F -m 7 -e 1e-5 -b 10 -v 10 -a 5

或diamond：

diamond makedb --in Plants.fa -d Plants.fa
diamond blastp --evalue 1e-5 --threads 4 --outfmt 5 -q query.fa  -d Plants.fa.dmnd -o blastout.nr --seg no --max-target-seqs 20 --more-sensitive -b 0.5 --salltitles

码农公寓