構(gòu)建ncbi——nr本地庫

網(wǎng)上方案都是用ncbi-blast+程序包的updata_blastdb.pl自動下載天揖，然而我下了兩天夺欲，一直斷···被逼無奈下，我只能用wget命令通過FTP地址去下載了

先附上幾個有用的地址：
ncbi blast database 下載方式介紹
 ncbi blast database FTP地址
 wget高級用法
現(xiàn)在我們正式開始····

1.庫文件列表< nrftplist >準備

import os
from ftplib import FTP 
import re
workdir='/Users/gushanshan/server/ncbi_db'
os.chdir(workdir)


ftpsite="ftp.ncbi.nlm.nih.gov/blast/db/"
linepo=ftpsite.index("/")
ftp=FTP(ftpsite[0:linepo])
ftp.login()
ftp.cwd(ftpsite[linepo+1:])
dirs=ftp.nlst()
for dir_i in dirs:
    inter_i=re.match("^nr",dir_i)
    if inter_i is not None:
        print(dir_i)
        filepo='ftp://'+ftpsite+dir_i
        os.system('echo \"'+filepo+';type=i\"'+' >> nrftplist')

ftp.close()

type=i指定binary mode：Pre-formatted databases must be downloaded using the update_blastdb.pl script or via FTP in binary mode

2 庫文件下載

如果用wget -i nrftplist下載今膊，程序會按次序一個一個下些阅，太費時間。所以我按照文件大小把nrftplist拆成了6個子文件斑唬，每個子文件總下載量差不多

wget -b -c -I nrftplist_1
wget -b -c -I nrftplist_2
wget -b -c -I nrftplist_3
wget -b -c -I nrftplist_4
wget -b -c -I nrftplist_5
wget -b -c -I nrftplist_6

-c 可持續(xù)斷點下載市埋，中途失敗可續(xù)接；
-i nrftplist對應(yīng)子文件恕刘；
-b 后臺下載缤谎。

速度似乎并沒有上去······，從中午下到了午夜 /("▔□▔)/

可以用下方程序檢查下載到哪里了

tail wget-log
tail wget-log.1
tail wget-log.2
tail wget-log.3
tail wget-log.4
tail wget-log.5

其實我覺得連寫6個wget的方式很笨雪营，但還不知道怎么能更簡潔些···還在學習中

3.文件完整性下載

因為文件太大弓千，下載時間太久，擔心中間出現(xiàn)什么差錯献起，用md5檢查下：

import os
import re
workdir='/Users/gushanshan/server/ncbi_db/nr'
os.chdir(workdir)
files=os.listdir()
for file in files:
    if re.search("gz$",file) is not None:
        print(file)
        os.system('md5sum '+file+' >> md5_file')
    elif re.search("md5$",file) is not None:
        print(file)
        os.system('cat '+file+' >> md5_md5')

太慢了洋访，我好焦灼···

如果md5_file和md5_md5結(jié)果一樣，說明文件很完整谴餐，可以繼續(xù)后續(xù)操作姻政。

4. 庫文件解壓

nohup gzip *.tar.gz -d &
nohup tar -xvf nr....tar &

5.本地庫配置——.ncbirc文件

NCBI官方文檔

以分號開頭的為注釋句；
可以放在當前工作路徑或用戶的home目錄岂嗓，我偏向于后者汁展；

; Start the section for BLAST configuration
[BLAST]
; Specifies the path where BLAST databases are installed
BLASTDB=/picb/evolgen/users/gushanshan/ncbi_db/nr
; Specifies the data sources to use for automatic resolution 
; for sequence identifiers 
DATA_LOADERS=blastdb 
; Specifies the BLAST database to use resolve protein sequences 
BLASTDB_PROT_DATA_LOADER=/picb/evolgen/users/gushanshan/ncbi_db/nr
; Specifies the BLAST database to use resolve nucleotide sequences 
BLASTDB_NUCL_DATA_LOADER=/picb/evolgen/users/gushanshan/ncbi_db/nt

; Windowmasker settings
[WINDOW_MASKER]
WINDOW_MASKER_PATH=/picb/evolgen/users/gushanshan/ncbi_db/windowmasker
; end of file

6.nr子庫構(gòu)建

參考文章：
https://www.bioinfo-scrounger.com/archives/207/
http://www.chenlianfu.com/?p=2691

6.1.軟件安裝及文件下載

6.1.1. TaxonKit

下載地址：https://bioinf.shenwei.me/taxonkit/download/

wget https://github.com/shenwei356/taxonkit/releases/download/v0.6.0/taxonkit_linux_amd64.tar.gz
tar -zxvf taxonkit_linux_amd64.tar.gz 
mkdir -p $HOME/bin/
cp taxonkit $HOME/bin/

6.1.2. csvtk

下載地址：https://bioinf.shenwei.me/csvtk/download/

wget https://github.com/shenwei356/csvtk/releases/download/v0.20.0/csvtk_linux_amd64.tar.gz
tar -zxvf csvtk_linux_amd64.tar.gz
cp csvtk $HOME/bin/

6.1.3. prot.accession2taxid

該文件記錄NCBI的accession與taxid的對應(yīng)關(guān)系

wget ftp://ftp.ncbi.nih.gov/pub/taxonomy/accession2taxid/prot.accession2taxid.gz
gzip -dc prot.accession2taxid.gz > prot.accession2taxid

6.1.4. taxdump

文件包含 taxonomic lineage of taxa， type strains和material的信息，以及host information
下載地址

#1.下載
wget https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/new_taxdump/new_taxdump.tar.gz
wget https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/new_taxdump/new_taxdump.tar.gz.md5

#2.判斷文件完整性：如果兩個結(jié)果一樣食绿，則文件完整
md5sum new_taxdump.tar.gz#7458b27039e743d8bb4fcbf71ccaf4ac  new_taxdump.tar.gz
cat new_taxdump.tar.gz.md5#7458b27039e743d8bb4fcbf71ccaf4ac  new_taxdump.tar.gz

#3.解壓
mkdir ~/.taxonkit
tar zxf new_taxdump.tar.gz -C ~/.taxonkit
#里面包括citations.dmp侈咕，delnodes.dmp，division.dmp器紧，gencode.dmp耀销，typeoftype.dmp，merged.dmp铲汪，names.dmp熊尉，nodes.dmp，host.dmp掌腰，typematerial.dmp狰住，rankedlineage.dmp，fullnamelineage.dmp齿梁，taxidlineage.dmp催植，readme.txt
#詳細文檔：https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/new_taxdump/taxdump_readme.txt

文件轉(zhuǎn)入$HOME/bin能起作用的前提是$HOME/bin在你的環(huán)境變量里

6.2. 子庫構(gòu)建

以Lactobacillus為例

6.2.1 Taxonomy ID

NCBI物種數(shù)據(jù)庫里搜Lactobacillus，有對應(yīng)的Taxonomy ID
程序查找

grep -P "\|\s+[lL]actobacillus\w*\s*\|" ~/.taxonkit/names.dmp

Lactobacillus對應(yīng)的Taxonomy ID為1578

6.2.2.提取特定taxons下的所有taxid

taxonkit list --ids 1578 --indent "" > Lactobacillus.taxid.txt
wc -l Lactobacillus.taxid.txt #共2733個taxid

$nrdir是我存放prot.accession2taxid的地方

--ids是指定taxid士飒，--indent ""是將所列出的taxid左邊的空格去除查邢，以左對齊排列

6.2.3. 提取 Lactobacillus.taxid所有的accession

csvtk在prot.accession2taxid文件中提取Lactobacillus.taxid所有的accession，共6694759條entry

cat /picb/evolgen/users/gushanshan/ncbi_db/prot.accession2taxid |csvtk -t grep -f taxid -P Lactobacillus.taxid.txt  |csvtk -t cut -f taxid,accession.version > Lactobacillus.taxid.acc.new.txt
wc -l Lactobacillus.taxid.acc.new.txt

cat /picb/evolgen/users/gushanshan/ncbi_db/prot.accession2taxid |csvtk -t grep -f taxid -P nr_Lactobacillus/Lactobacillus.taxid.txt  |csvtk -t cut -f taxid,accession.version > all_accession_taxid

6.2.4. 建立Lactobacillus子庫

blastdb_aliastool -seqidlist Lactobacillus.taxid.acc.new.txt -db nr -out nr_Lactobacillus -title nr_Lactobacillus

blastdb_aliastool必須是ncbi toolkit里的酵幕，比如ncbi-blast-2.10.1+。conda安裝的不行缓苛。

7.應(yīng)用

blast系列有以下幾種

摘自徐洲更hoptop簡書 o(*^＠^*)o

以npsA基因為例芳撒，尋找其在乳酸桿菌中的進化史。npsA protein fasta 序列下載網(wǎng)址

7.1. blast

#!/bin/bash

db_Lactobacillus=/picb/evolgen/users/gushanshan/probiotics/13tree/npsABC/genetree/nr_Lactobacillus

blastp -num_threads 2 -word_size 6 -gapopen 11 -gapextend 1 -threshold 21 -window_size 40 -db $db_Lactobacillus/nr_Lactobacillus -query npsa_wcfs1.fasta -out npsA_blast_result.txt -outfmt 6

得到2,459個hit

按照ncbi blastp的參數(shù)設(shè)置本地參數(shù)

7.2. 取出相似序列

7.2.1. 閾值選擇

根據(jù) An Introduction to Sequence Similarity (“Homology”) Searching的介紹未桥，將閾值選為bit score> 50笔刹。

得到2,458個hit。

500條同屬序列+3條外群序列+npsA=504條

7.2.2. 抽取序列

epost -db protein -input high_similarity_acc.txt | efetch -format fasta > high_similarity.fasta

epost -db protein -input allacc.txt | efetch -format ipg > allan_tax

7.3. 多序列比對

mafft --auto --thread -1 high_similarity.fasta > mafft_npsA.fasta

附錄1. taxId2taxName

為了了解哪些物種中存在npsA的相似蛋白質(zhì)冬耿，我們需要將6.2.2中得到的taxid對應(yīng)到taxName舌菜。

因為1)轉(zhuǎn)換過程很快，2)該過程不需要重復(fù)操作亦镶，所以我認為沒必要專門寫轉(zhuǎn)換程序

轉(zhuǎn)換網(wǎng)址：https://www.ncbi.nlm.nih.gov/Taxonomy/TaxIdentifier/tax_identifier.cgi

示例

附錄2. 蛋白質(zhì)accession對應(yīng)的具體strain

很多蛋白質(zhì)對應(yīng)species日月，但實際只屬于某些特定strain＄凸牵可以通過identical protein

#單個id
esearch -db protein -query 'WP_011101060.1' | efetch -format ipg| grep "RefSeq"
#多個id
epost -db protein -input test | efetch -format ipg| grep "RefSeq"