如何利用NR庫(kù)快速進(jìn)行物種鑒定

<meta charset="utf-8">

引言

在做基因或者蛋白質(zhì)功能注釋的時(shí)候，NR數(shù)據(jù)庫(kù)包含非常重要的信息。但是有沒(méi)有人跟我一樣遇到相同的苦惱，就是在運(yùn)行的時(shí)候特別耗時(shí)，經(jīng)過(guò)一番搜索后，發(fā)現(xiàn)可以根據(jù)實(shí)際數(shù)據(jù)建立自己的NR子庫(kù)，這樣運(yùn)行效率將大幅提高。那么怎么創(chuàng)建NR子庫(kù)呢，跟著我一起來(lái)做吧。

01 NR 數(shù) 據(jù) 庫(kù) 介 紹

目前我們可以從很多數(shù)據(jù)庫(kù)中獲取蛋白質(zhì)的序列信息，比如GenPept, Swissprot, PIR, PDF, PDB 和 NCBI RefSeq等。但是在這些數(shù)據(jù)庫(kù)之間，蛋白質(zhì)序列存在冗余性，為了解決這個(gè)問(wèn)題，NCBI構(gòu)建了一個(gè)非冗余的蛋白質(zhì)序列數(shù)據(jù)庫(kù)，即NR（non-redundant proteins）。

我們可以看到最新版的NR庫(kù)是剛剛更新的，壓縮后的文件存儲(chǔ)已經(jīng)達(dá)到了73G，這是相當(dāng)恐怖的。跟著來(lái)建立自己的NR子庫(kù)吧~

image

02 準(zhǔn) 備 工 作

首先我們下載NR數(shù)據(jù)庫(kù)及需要用到的軟件。

2.1 下載NR數(shù)據(jù)庫(kù)

https://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nr.gz
下載后需要解壓，fasta格式文件，存放路徑 mypath/nr 。
既然是非冗余數(shù)據(jù)庫(kù)，那么NR數(shù)據(jù)庫(kù)中序列合并必須同時(shí)滿足以下兩個(gè)原則：

  • 兩條序列有相同的長(zhǎng)度；
  • 兩條序列上的每個(gè)殘基都是相同的。

2.2 下載分類(lèi)數(shù)據(jù)庫(kù)

1 mkdir  mypath/taxdump
2 tar  -zxvf  mypath/taxdump.tar.gz   -C mypath/taxdump
  • taxdump 目錄中有兩個(gè)重要文件：
    names.dmp：記錄物種名及其分類(lèi)編號(hào)
    nodes.dmp：記錄分類(lèi)編號(hào)的分類(lèi)節(jié)點(diǎn)信息

2.3 下載accession與taxid的對(duì)應(yīng)關(guān)系

下載鏈接為：
https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/accession2taxid/prot.accession2taxid.gz
1gzip -dc prot.accession2taxid.gz >prot.accession2taxid

2.4 其他工具

03 構(gòu) 建 命 令

我們以構(gòu)建病毒子庫(kù)為例，首先要知道病毒的tax id為10239，這個(gè)需要要先驗(yàn)知識(shí)，或者在taxdump/names.dmp文件中查找。

1.提取病毒包含的所有物種編號(hào)

1 mypath/taxonkit/taxonkit list -j 2 --ids 10239 --indent "" --data-dir mypath/taxdump/ >mypath/virus.list
這行命令運(yùn)行很快，大概十秒鐘，看下結(jié)果，總共提取出來(lái)205975個(gè)taxid。

-j ，線程，給2足以；
--ids，需要提取的分類(lèi)的taxid；
--data-dir，該目錄下必須包含文件names.dmp和nodes.dmp；
--indent，提取的物種編號(hào)縮進(jìn)位置，這個(gè)參數(shù)很重要，記得一定要設(shè)置為空 ""

2.根據(jù)taxid獲取對(duì)應(yīng)accession.version

1 zcat mypath/prot.accession2taxid.gz | mypath/csvtk/csvtk -t grep -f taxid -P mypath/virus.list | mypath/csvtk/csvtk -t cut -f accession.version > mypath/virus.taxid.acc.txt
這一步運(yùn)行也比較快，我這不到10分鐘就運(yùn)行結(jié)束了，這一步提取出來(lái)的accession就能對(duì)應(yīng)到NR數(shù)據(jù)庫(kù)中的標(biāo)識(shí)符，總共提取出來(lái)了6065494行數(shù)據(jù)。

很顯然這個(gè)數(shù)目比上一步taxid的數(shù)目多了很多，這是因?yàn)橐粋€(gè)accession可能會(huì)對(duì)應(yīng)多個(gè)taxid，這個(gè)我們?cè)趍ypath/prot.accession2taxid.gz 中就可以看到。

3.NR數(shù)據(jù)庫(kù)建索引

1 mypath/ncbi-blast/bin/makeblastdb -in nr -dbtype prot -parse_seqids
建索引的步驟非常耗時(shí)，可能會(huì)花費(fèi)好幾個(gè)小時(shí)，需要耐心等待。

4.提取NR子庫(kù)

1 mypath/ncbi-blast/bin/blastdb_aliastool -seqidlist mypath/virus.taxid.acc.txt -db nr -out nr_virus -title nr_virus
最后這一步也會(huì)比較耗時(shí)，大概運(yùn)行了2個(gè)多小時(shí)。運(yùn)行結(jié)束后只生成了文件 mypath/nr_virus.pal，不要懷疑自己，確實(shí)只生成一個(gè)文件。

04 功 能 注 釋

對(duì)于NR數(shù)據(jù)庫(kù)，一般都是blast直接比對(duì)，獲取基因或者蛋白質(zhì)的功能注釋信息，更近一步，如果有Taxonomy 數(shù)據(jù)庫(kù)的數(shù)據(jù)，我們還可以做物種注釋。

接下來(lái)我們?cè)囋囉锰崛〉牟《咀訋?kù)做功能注釋吧，測(cè)試數(shù)據(jù)我選取了100條蛋白質(zhì)序列，命名為test.fa。

運(yùn)行命令：
1 mypath/ncbi-blast/bin/blastp -query mypath/test.fa -db mypath/nr_virus -evalue 1e-5 -num_alignments 5 -outfmt 6 -out mypath/test.blastp.out

運(yùn)行結(jié)果是m8格式的文件如下圖：

image

該文件包含12列，每列信息簡(jiǎn)單介紹如下：

  1. qseqid：查詢序列ID標(biāo)識(shí)
  2. sseqid：比對(duì)上的目標(biāo)序列ID標(biāo)識(shí)
  3. pident：序列比對(duì)的一致性百分比
  4. length：符合比對(duì)的比對(duì)區(qū)域的長(zhǎng)度
  5. mismatch：比對(duì)區(qū)域的錯(cuò)配數(shù)
  6. gapopen：比對(duì)區(qū)域的gap數(shù)目
  7. qstart：比對(duì)區(qū)域在查詢序列(qseqid)上的起始位點(diǎn)
  8. qend：比對(duì)區(qū)域在查詢序列(qseqid)上的終止位點(diǎn)
  9. sstart：比對(duì)區(qū)域在目標(biāo)序列(sseqid)上的起始位點(diǎn)
    10.send：比對(duì)區(qū)域在目標(biāo)序列(sseqid)上的終止位點(diǎn)
    11.evalue：比對(duì)結(jié)果的期望值，evalue越小，越有可能是真實(shí)的相似序列
    12.bitscore：比對(duì)結(jié)果的bit score值

注 意 ! ! !

  • 因?yàn)槲覀冇玫膎cbi-blast+中的blastp，-outfmt 要設(shè)置為6，如果你用的是blast中的blastp工具，要設(shè)置 -m 8；
  • 一般我們比較關(guān)注第3、11、12列信息；
  • 一條序列比對(duì)可能有多條結(jié)果，我們可以按照bitscore篩選得分最高的一條即可。

子庫(kù)運(yùn)行時(shí)間不到50分鐘，而選擇整個(gè)NR數(shù)據(jù)庫(kù)比對(duì)時(shí)，將近4個(gè)小時(shí)才結(jié)束，很顯然，選擇子庫(kù)比對(duì)大大提升了比對(duì)效率。因此，如果我們對(duì)研究的物種比較確定，提取子庫(kù)再進(jìn)行比對(duì)，將會(huì)提升運(yùn)行效率降低集群壓力。

05 物 種 注 釋

5.1 合并物種分類(lèi)文件

names.dmp，nodes.dmp 這兩個(gè)文件的內(nèi)容我們上文已經(jīng)介紹過(guò)，現(xiàn)在我們先來(lái)處理一下。將這兩個(gè)文件的信息整合在一起，便于后續(xù)應(yīng)用。

這里介紹網(wǎng)上一個(gè)腳本，首先運(yùn)行
1 git clone https://gitee.com/wangshun1121/TaxonomyPickUp.git
命令后會(huì)在當(dāng)前路徑下生成目錄mypath/TaxonomyPickUp/，然后將names.dmp，nodes.dmp 兩個(gè)文件拷貝到目錄mypath/TaxonomyPickUp/。
最后運(yùn)行
1 perl NodeExtract.pl >TaxInfo.txt
大概一兩分鐘就得到結(jié)果，它存儲(chǔ)了NCBI的Taxonomy物種分類(lèi)的詳細(xì)層級(jí)信息，格式如下：

image
  • 第一列為T(mén)axID，
  • 第二列為該TaxID的拉丁學(xué)名。
  • 第三列表示該TaxID所處的層級(jí)：即經(jīng)過(guò)追溯多少次Parents到達(dá)分類(lèi)級(jí)別的最高層級(jí)。
  • 第四列則是TaxID向上追溯的詳細(xì)信息。不同分類(lèi)層級(jí)之間以“-”間隔，最左邊為該TaxID代表的層級(jí)，例如7 “Azorhizobium caulinodans”表示species，其上一級(jí)的TaxID為6，層級(jí)為genus，再向上為335928的family，依次類(lèi)推。

5.2 進(jìn)行物種注釋

比對(duì)之后我們可以繼續(xù)做物種分類(lèi)注釋?zhuān)枰玫缴弦徊絙lastp的結(jié)果，prot.accession2taxid.gz，names.dmp，nodes.dmp 這幾個(gè)文件了。接下來(lái)，我們只需要寫(xiě)個(gè)簡(jiǎn)單的腳本，就能從blastp比對(duì)結(jié)果mypath/test.blastp.out獲取每條reads對(duì)應(yīng)的物種分類(lèi)信息了。這里簡(jiǎn)單說(shuō)下邏輯：

    1. 從mypath/test.blastp.out文件前兩列獲取reads id 和accession對(duì)應(yīng)關(guān)系；
    1. prot.accession2taxid.gz文件包含accession和taxid的對(duì)應(yīng)關(guān)系，結(jié)合1可以得到每條reads id和taxid的對(duì)應(yīng)關(guān)系；
    1. 上述TaxInfo.txt包含了taxid和物種的對(duì)應(yīng)關(guān)系，結(jié)合2就能獲取reads id和物種的關(guān)系了，還可以方便篩選所有層級(jí)的分類(lèi)單位。

通過(guò)這個(gè)詳細(xì)的描述，相信大家都已經(jīng)掌握了整個(gè)過(guò)程，快動(dòng)手試試吧。

作者：生信阿拉丁
鏈接：http://www.itdecent.cn/p/45fdb5cf930a
來(lái)源：簡(jiǎn)書(shū)
著作權(quán)歸作者所有。商業(yè)轉(zhuǎn)載請(qǐng)聯(lián)系作者獲得授權(quán)，非商業(yè)轉(zhuǎn)載請(qǐng)注明出處。

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書(shū)系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容