最近做有關(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)的分類(lèi)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.txtcut -f 1
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
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ù)的頭文件，例如
所以使用的話，還需要寫(xiě)個(gè)perl把這些序列拆開(kāi)，最終形成nt.anmail.fa.gz
8 如果直接想構(gòu)建子庫(kù)，那么沒(méi)必要搞序列，直接運(yùn)行
blastdb_aliastool -gilist
NT -out nt_animal -title nt_animal
轉(zhuǎn)自：http://www.itdecent.cn/p/ec7773da27b9
