构建ncbi——nr本地库

网上方案都是用ncbi-blast+程序包的updata_blastdb.pl自动下载，然而我下了两天，一直断···被逼无奈下，我只能用wget命令通过FTP地址去下载了

先附上几个有用的地址：
ncbi blast database 下载方式介绍
 ncbi blast database FTP地址
 wget高级用法
现在我们正式开始····

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 可持续断点下载，中途失败可续接；
-i nrftplist对应子文件；
-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.文件完整性下载

因为文件太大，下载时间太久，担心中间出现什么差错，用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结果一样，说明文件很完整，可以继续后续操作。

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子库构建

参考文章：
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的对应关系

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.判断文件完整性：如果两个结果一样，则文件完整
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

文件转入$HOME/bin能起作用的前提是$HOME/bin在你的环境变量里

6.2. 子库构建

以Lactobacillus为例

6.2.1 Taxonomy ID

NCBI物种数据库里搜Lactobacillus，有对应的Taxonomy ID
程序查找

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

Lactobacillus对应的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.应用

blast系列有以下几种

摘自徐洲更hoptop简书 o(*^＠^*)o

以npsA基因为例，寻找其在乳酸杆菌中的进化史。npsA protein fasta 序列下载网址

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的参数设置本地参数

7.2. 取出相似序列

7.2.1. 阈值选择

根据 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的相似蛋白质，我们需要将6.2.2中得到的taxid对应到taxName。

因为1)转换过程很快，2)该过程不需要重复操作，所以我认为没必要专门写转换程序

转换网址：https://www.ncbi.nlm.nih.gov/Taxonomy/TaxIdentifier/tax_identifier.cgi

示例

附录2. 蛋白质accession对应的具体strain

很多蛋白质对应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"