最近做有關(guān)小鼠腸道微生物宏基因組，遇到兩個(gè)問(wèn)題何暮，一是數(shù)據(jù)量太大奄喂，二是宿主污染嚴(yán)重。
估計(jì)宿主污染至少80%左右海洼，因而就想通過(guò)一些方法跨新，例如kraken、bowtie等把宿主污染去除坏逢。
那么就有一個(gè)問(wèn)題域帐，如何選擇去除污染的數(shù)據(jù)庫(kù)呢？
思來(lái)想去是整，還是從NT庫(kù)入手肖揣，打算把NT庫(kù)所有動(dòng)物的序列或者所有小鼠的序列提取出來(lái)，做成一個(gè)子庫(kù)浮入，用來(lái)去除宿主污染龙优。
百度了一下提取子庫(kù)的方法，大多都是人云亦云事秀，干脆還是自己整理整理彤断。下面是一些步驟
1 首先下載NCBI的taxonomy數(shù)據(jù)庫(kù)
下載完解壓縮，其中names.dmp和nodes.dmp兩個(gè)文件很重要易迹，是后續(xù)提取子庫(kù)的基礎(chǔ)
2 下載NCBI的TaxonKit軟件宰衙，http://bioinf.shenwei.me/taxonkit/download/，linux系統(tǒng)直接解壓
而后把names.dmp和nodes.dmp兩個(gè)文件直接cp到~/.taxonkit下睹欲，其余的.dmp也可一并cp到~/.taxonkit下
cp taxdump/* ~/.taxonkit
3?下載NCBI的csvtk軟件供炼，http://bioinf.shenwei.me/csvtk/download/，linux系統(tǒng)也是直接解壓窘疮，即可使用
4 （選擇性步驟）NCBI taxonomy數(shù)據(jù)庫(kù)下還有accession2taxid庫(kù)袋哼，這個(gè)庫(kù)里面也有蛋白以及核酸的accession以及對(duì)應(yīng)的分類id，但是經(jīng)過(guò)嘗試考余，采取這種方法提取的子庫(kù)序列往往出乎意料的少先嬉，很可能是該庫(kù)的accession與NT/NR庫(kù)的accession不一致轧苫，前者可能冗余更多楚堤，因此該方法可忽略疫蔓，見(jiàn)仁見(jiàn)智吧，下面給個(gè)例子身冬，例如：
#從taxonomy數(shù)據(jù)庫(kù)中的nucl_wgs.accession2taxid提取accession號(hào)
pigz -dc prot.accession2taxid.gz \
? ? | csvtk grep -t -f taxid -P $id.taxid.txt \
? ? | csvtk cut -t -f accession.version,taxid \
? ? | sed 1d \
? ? > $id.acc2taxid.txt
cut -f 1 $id.acc2taxid.txt > $id.acc.txt
5 查看衅胀，并提取動(dòng)物的taxid
grep -P "\|\s+[aA]nimal\w*\s*\|" ~/.taxonkit/names.dmp
33208
taxonkit list --ids 33208 --indent "" > $id.taxid.txt
6 從下載好的NT庫(kù)提取對(duì)應(yīng)的accession
blastdbcmd -db $NT -entry all -outfmt "%a %T" | csvtk grep -d ' ' -D ' ' -f 2 -P $id.taxid.txt \
? ? | cut -d ' '? -f 1 \
? ? > $id.acc.txt
7 從NT庫(kù)提取完整的nt序列，并提取子庫(kù)序列
blastdbcmd -db $NT -dbtype nucl -entry all -outfmt "%f" -out - | pigz -c > nt.fa.gz
time cat <(echo) <(pigz -dc nt.fa.gz) \
? ? | perl -e 'BEGIN{ $/ = "\n>"; <>; } while(<>){s/>$//;? $i = index $_, "\n"; $h = substr $_, 0, $i; $s = substr $_, $i+1; if ($h !~ />/) { print ">$_"; next; }; $h = ">$h"; while($h =~ />([^ ]+ .+?) ?(?=>|$)/g){ $h1 = $1; $h1 =~ s/^\W+//; print ">$h1\n$s";} } ' \
? ? | seqkit grep -f $id.acc.txt -o nt.$id.fa.gz
需要注意的是酥筝，這里又使用了seqkit軟件滚躯。這種從NT庫(kù)中還原的nt.fa序列里面有很多重復(fù)的頭文件，例如
所以使用的話嘿歌，還需要寫個(gè)perl把這些序列拆開(kāi)掸掏，最終形成nt.anmail.fa.gz
8 如果直接想構(gòu)建子庫(kù)，那么沒(méi)必要搞序列宙帝，直接運(yùn)行
blastdb_aliastool -gilist $id.acc.txt -db $NT -out nt_animal -title nt_animal
