BLAST及NT库本地化

Log

下载NCBI的NT数据库到服务器本地
Blast比对结果中输出物种信息
用TaxonKit、blastdbcmd、makeblastdb建立子库

NT库本地化

下载文件

  1. 可选文件格式

    1. blast index文件
      https://ftp.ncbi.nlm.nih.gov/blast/db中的nt.*.tar.gz 。
      下载完成后，MD5验证，解压文件，就可以直接用于blast。index中已经包含物种Taxonomy信息。但是压缩包中没有Fasta文件，需要的话可以用blastdbcmd提取。
    2. Fasta文件
      https://ftp.ncbi.nlm.nih.gov/blast/db/FASTA中的nt.gz 。
      下载完成并解压后，需要自建blast index；如果结果中还需要物种Taxonomy信息，建index前还需要准备好TaxIDMapFile文件。

  2. 下载blast index文件
    wget下载有点问题，不知道是服务器还是网络问题，wget下载的文件总是MD5校验失败；有时候重新下又能校验成功。
    aspera下载不行，链接NCBI失败，可能是必须开放某个指定端口。
    最后是用的rsync下载。

    1. 文件列表下载 
      1
      2
      3
      rsync --no-motd --files-from=<文件列表> <文件源> <下载存放路径>
      # 示例
      rsync --no-motd --files-from=***/File.list rsync://ftp.ncbi.nlm.nih.gov/blast/db/ ***/NT_Database
      File.list（每个文件1列）
      File.list
    • 实际下载文件路径是<文件源>加上File.list中的内容，完整路径如：rsync://ftp.ncbi.nlm.nih.gov/blast/db/nt.22.tar.gz。
    • 当File.list中的内容包含子目录，<下载存放路径>中会自动生成对应层级目录，例如File.list中的内容是“blast/db/nt.22.tar.gz”，则下载的nt.22.tar.gz会保存在<下载存放路径>/blast/db中。
    1. 并行下载
      推荐这种方式，文件列表逐个下载太慢，这个每行加“&”后台并行下载文件比较快。 
      1
      2
      3
      4
      rsync --no-motd <文件源> <下载存放路径>
      # 示例
      rsync --no-motd rsync://ftp.ncbi.nlm.nih.gov/blast/db/nt.22.tar.gz ***/NT_Database &
      rsync --no-motd rsync://ftp.ncbi.nlm.nih.gov/blast/db/nt.23.tar.gz ***/NT_Database &

MD5验证

https://ftp.ncbi.nlm.nih.gov/blast/db中还有nt.*.tar.gz.md5文件，用于MD5验证。
全部nt..tar.gz和nt..tar.gz.md5文件都下载完后，md5文件合并一起验证。
nt.all.md5.check中，文件名后是OK即文件完整；如果是Fail，则不完整，需要重新下载该文件。 

1
2
3
4
# 合并所有md5文件
cat nt.*.tar.gz.md5 > nt.all.md5
# md5验证
md5sum -c nt.all.md5 > nt.all.md5.check

解压文件

tar -zxvf nt.*.tar.gz
解压后，nt.*.tar.gz文件删除或者备份都可以。
一般是删除？除了需要迁移到别的服务器，好像也没有能再用到的地方了。

Blast比对

blast结果中Subject Seq-id（NT库的序列ID）格式为“gi|384474605|emb|HE793683.1|”，含有GI numbers和GenBank accession numbers；不包含Taxid信息。
如果blast index中有Taxonomy信息，可以在输出格式增加staxids sscinames，使结果额外输出物种的Taxid和物种名称。 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# Shell

sample="Input"
query_fa="Input.fa"

database="<NT库目录>/nt"
identity="98"
cov_hsp="90"
evalue="1e-5"

threads="8"

# 如果不export BLASTDB的话，要在bashrc声明BLASTDB，否则报错
# 在-outfmt添加staxids sscinames，使结果额外输出物种Taxid、物种名称
export BLASTDB=<NT库目录> && blastn -query ${query_fa} -db ${database} -out ${sample}.noheader -outfmt "6  qaccver saccver pident length mismatch gapopen qstart qend sstart send evalue bitscore staxids sscinames stitle" -num_threads ${threads} -evalue ${evalue} -perc_identity ${identity} -qcov_hsp_perc ${cov_hsp}

# 表头
echo -e "Query Seq-id\tSubject Seq-id\tIdentity\tAlignment length\tMismatch\tGap Open\tQuery Start\tQuery End\tSubject Start\tSubject End\tEvalue\tBitscore\tTaxonomy ID\tScientific Name\tSeq Title" > ${sample}.tsv

# 合并blast结果和表头
cat ${sample}.noheader >> ${sample}.tsv

# 删除没有表头的blast结果
rm -f ${sample}.noheader

构建子库

  • 用blastdbcmd可以根据Taxid从NT库抽取相关序列，从而构建子库。
  • 用TaxonKit可以获取某个Taxid节点以下所有Taxid，也就是说不止可以抽取单个物种的序列，还可以提供Taxid列表来抽取整个目、科、属…的序列。
  • TaxonKit需要taxdump.tar.gz。NT库和taxdump.tar.gz最好日期一致，因为Taxonomy会更新，可能导致NT库的Taxid和TaxonKit提取的有区别。具体见TaxonKit的TaxID changelog。
  • 如果要用Taxid统计blast结果中的物种信息，需要注意有些亚种/株的Taxid与物种（Species）等级Taxid不一样。如Streptococcus parasanguinis ATCC 903（Taxid 888048）的Rank是strain；它是属于Streptococcus parasanguinis（Taxid 1318）这个物种的。用taxonkit lineage可以获取某个Taxid在物种等级的名称和Taxid。如echo 888048 | taxonkit lineage | taxonkit reformat -t -f \"{s}\" | cut -f 3,4。
  • 用子库作为参考序列blast时，将-db参数换成{Index_Dir}/${Name}.fa即可。

以链球菌属举例，在NCBI的Taxonomy数据库可以找到Streptococcus的Taxid是1301。 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
# Shell

Name="Streptococcus"
Taxid="1301"
Index_Dir="Streptococcus_Index"

# 获取该Taxid节点以下所有子节点的Taxid
taxonkit list -I "" --ids $Taxid > ${Index_Dir}/${Name}.taxid

# 提取相关序列
blastdbcmd -db <NT库目录>/nt -taxidlist ${Index_Dir}/${Name}.taxid > ${Index_Dir}/${Name}.fa

# 提取GenBank accession numbers（GI号也可以）和Taxid的对应关系文件
blastdbcmd -db <NT库目录>/nt -taxidlist ${Index_Dir}/${Name}.taxid -outfmt "%a %T" > ${Index_Dir}/taxid_map.txt

# 建blast index，通过-parse_seqids -taxid_map在index中附加物种Taxonomy信息
makeblastdb -dbtype nucl -in ${Index_Dir}/${Name}.fa -out ${Index_Dir}/${Name}.fa -parse_seqids -taxid_map ${Index_Dir}/taxid_map.txt

# 提取序列的GenBank accession numbers、Taxid、物种名称（比对用不上，就是看看参考里有什么）
blastdbcmd -db <NT库目录>/nt -taxidlist ${Index_Dir}/${Name}.taxid -outfmt "%a  %T  %S" > ${Index_Dir}/${Name}.Name.tsv

taxid_map.txt（第1列GenBank accession numbers，第2列Taxid）
taxid_map.txt
Name.tsv（第1列GenBank accession numbers，第2列Taxid，第3列物种名称）
Name.tsv

假如用nt.gz建库

  • nt.gz解压后是NT库的fasta文件。如果不需要物种信息的话，直接makeblastdb建库就行。如果balst结果要有物种信息，则需要提供GenBank accession numbers和Taxid对应关系，类似上面的子库构建的-parse_seqids -taxid_map方式。
  • GenBank accession numbers和Taxid对应关系文件可以用nucl_gb.accession2taxid.gz。和taxdump.tar.gz一样原因，accession2taxid最好与NT库日期一致。
  • nucl_gb.accession2taxid.gz格式是tsv，第一列Accession，第二列Accession.version，第三列TaxId，第四列GI；整理成taxid_map格式后，用makeblastdb加-parse_seqids -taxid_map建库就行。

一些blastn参数

-task：默认是megablast，可选blastn、blastn-short、dc-megablast、megablast、rmblastn
-query：需要比对的fasta序列
-db：数据库路径，写到fasta前缀为止，如***/NT_Database/nt
-out：输出文件
-outfmt：输出格式，6是tab分割且无表头的格式，其他格式见blastn帮助文档
-num_threads：线程数，默认1；与-mt_mode相关，具体见这里
-evalue：e值过滤阈值，默认10
-perc_identity：identity百分比过滤阈值，取值0 ~ 100
-qcov_hsp_perc：hsp覆盖度百分比过滤阈值，取值0 ~ 100
-max_target_seqs：输出结果保留多少比对上的subject序列，默认500，建议大于等于5
-max_hsps：每个query序列比对上的每个subject序列保留多少个HSP，需要填写大于等于1，无默认值