昨天下午搗鼓了一下宏基因組物種注釋過(guò)程(基于nr庫(kù))，現(xiàn)在將整個(gè)流程記錄一下。
軟件需求：blast，diamond赛不，taxonkit(安裝自行百度)
構(gòu)建細(xì)菌子庫(kù)
blast方法可能會(huì)準(zhǔn)確點(diǎn)，但是它的速度簡(jiǎn)直讓我懷疑人生罢洲，倆種軟件的方法我都說(shuō)下吧踢故，因?yàn)槲冶葘?duì)的主要是細(xì)菌文黎，我首先想到是干脆按照網(wǎng)上的方法構(gòu)建一個(gè)細(xì)菌的子庫(kù)可能速度會(huì)更快點(diǎn)~
說(shuō)干就干
參考連接:https://www.bioinfo-scrounger.com/archives/207/;
https://bioinf.shenwei.me/taxonkit/tutorial/
我最初想法是構(gòu)建一張表，第一列是taxid殿较，后面7列跟著門綱目科屬種的名稱耸峭，這樣的話我們根據(jù)比對(duì)結(jié)果中的taxid，就直接可以得到物種各層級(jí)的信息淋纲，相信這也是大多數(shù)人的想法劳闹，就像下圖這樣：
這里我們就可以使用taxonkit這個(gè)好輪子去方便的達(dá)成這一目的，代碼如下：
1.#輸出細(xì)菌的taxid
taxonkit list --ids 2 --indent "" > bacteria.taxid.txt
#33090為植物 2為細(xì)菌 2157為Archaea洽瞬；10239 為Viruses本涕；Eukaryota為2759；Fungi 4751
#有時(shí)間的話可以將這幾大類的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,現(xiàn)在已經(jīng)得到了這張表菩颖，接下來(lái)我們就要去比對(duì)，得到每個(gè)基因的taxid
先來(lái)構(gòu)建細(xì)菌子庫(kù)吧
需要文件：
ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/accession2taxid/prot.accession2taxid.gz 蛋白acc號(hào)和taxid對(duì)應(yīng)文件
方法參考自爪哥文檔：https://bioinf.shenwei.me/taxonkit/tutorial/
#還是像上面那樣先得到細(xì)菌的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 :
#直接構(gòu)建(速度慢)
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版本 低版本沒(méi)有seqidlist) 這步結(jié)果就只有一個(gè)pal文件
blastdb_aliastool -seqidlist bacteria.taxid.acc.txt -db /home/software/nr-2019-12-18/nr -out nr_bac -title nr_bac
Option 3:
#采用了分割并行的策略加上爪哥自己寫的腳本为障，速度更加快晦闰，建議大家去看爪哥的源文檔，這里不再贅述.
序列比對(duì)
blastp
我們拿上面Option2得到的文件去進(jìn)行比對(duì)鳍怨，命令行如下：
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
#記得結(jié)果加上 staxid 得到物種taxid信息
diamond
由于索引庫(kù)不兼容呻右，我們將blastcmd抽提出來(lái)的nr庫(kù)，用diamond先構(gòu)建索引庫(kù)
要想得到taxid和種名信息鞋喇，需要構(gòu)建的時(shí)候額外增加倆個(gè)參數(shù)--taxonmap和--taxonnodes
1是我們上述說(shuō)的 蛋白acc號(hào)和taxid的對(duì)應(yīng)文件prot.accession2taxid.gz
2是存儲(chǔ)有taxonomy數(shù)據(jù)庫(kù)的層級(jí)文件taxdmp.zip
注意wget的時(shí)候會(huì)出現(xiàn)連接超時(shí)情況声滥，多等幾次即可.
diamond makedb --in nr.fa -d Dimond_nr --taxonmap prot.accession2taxid.gz --taxonnodes nodes.dmp
(#nodes.dmp文件是taxdmp.zip壓縮包內(nèi)的文件,這個(gè)地方直接跟這個(gè)文件即可，否則建庫(kù)不完整)
#比對(duì)
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
補(bǔ)充一點(diǎn)：
我們可以先將有taixd也有種名信息的提取出來(lái)侦香，過(guò)濾掉沒(méi)有taxid和種名信息的基因醒串，因?yàn)槲募Y(jié)果中可能會(huì)出現(xiàn)一種情況是該基因比對(duì)結(jié)果沒(méi)有taxid，但是后面匹配到了種名信息鄙皇，這時(shí)候我們?cè)賮?lái)用taxonkit這個(gè)軟件根據(jù)種名去得到這些物種的taixid，這樣我們的結(jié)果里會(huì)多一點(diǎn)仰挣，這時(shí)候還沒(méi)有匹配到taxid的再將其過(guò)濾掉伴逸。
這里我將腳本貼出來(lái)
#!usr/bin env python
import re
import sys
#根據(jù)nr結(jié)果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()
好了膘壶，上面兩種軟件得到結(jié)果都會(huì)有taxid信息和種名了错蝴，由于是比對(duì)的子庫(kù)文件，速度也會(huì)有所提升颓芭，我們可以拿著taxid去匹配最后那張表的信息就能匹配種水平上面幾個(gè)層級(jí)的taxonomy的信息顷锰，或者好好利用taxonkit這個(gè)好輪子也是可以得到的~~~
下面是剛創(chuàng)建個(gè)人公眾號(hào)，會(huì)定時(shí)更新R亡问、linux官紫、python肛宋，組學(xué)方面的學(xué)習(xí)內(nèi)容，請(qǐng)多多支持呦~~
