通过blast比对NCBI的nt库来鉴定物种类别

2022-03-02 本文已影响0人莫讠

美吉生物信息云流程ktps://cloud.majorbio.com/page/project/task.html?project_id=0nccl4374u5lopn15761jdrrhl

有一个CVTree的工具可以参考一下
http://cvtree.online/v3/php/projectpage.php?project=20211228_1525_59913&ptype=cvtree&db=20160125

知乎上面建议是用nt库，然后再用megan
参考链接：https://zhuanlan.zhihu.com/p/91633444

强哥给出的建议是blast以及nt库进行比对：

以下是Linux中Blast+的使用：
https://huans.github.io/2017/12/16/blast/

https://huans.github.io/2017/12/16/blast/

对于组装的文献参考：
https://genomebiology.biomedcentral.com/articles/10.1186/s13059-020-1938-2#Sec16

对于Blast的安装

对于NCBI blast的下载地址：
https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/

wget -c https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ncbi-blast-2.12.0+-x64-linux.tar.gz

安装完成之后在bin目录下还有一个updata_blastdb.pl的一个工具，但是执行这个工具的话，还需要再安装一个执行perl脚本的软件

sudo apt-get install perl-doc

之后依照 https://zhuanlan.zhihu.com/p/91633444 中的方法使用update_blastdb.pl脚本显示NCBI里面所有的库并安装nt库

显示所有可安装的库

(base) qianwj@ubuntu-NF5280M5:~/software/blast/ncbi-blast-2.12.0+/bin $ ./update_blastdb.pl --showall

Connected to NCBI
16S_ribosomal_RNA
18S_fungal_sequences
28S_fungal_sequences
Betacoronavirus
ITS_RefSeq_Fungi
ITS_eukaryote_sequences
LSU_eukaryote_rRNA
LSU_prokaryote_rRNA
SSU_eukaryote_rRNA
env_nt
env_nr
human_genome
landmark
mito
mouse_genome
nr
nt
pataa
patnt
pdbaa
pdbnt
ref_euk_rep_genomes
ref_prok_rep_genomes
ref_viroids_rep_genomes
ref_viruses_rep_genomes
refseq_select_rna
refseq_select_prot
refseq_protein
refseq_rna
swissprot
tsa_nr
tsa_nt
taxdb

有个+比较烦，所以对文件名进行修改

mv ncbi-blast-2.12.0+ ncbi-blast-2.12.0

如果想直接用update_blastdb.pl进行下载的话，会出现下载中断的问题，比如下面这种的

(base) qianwj@ubuntu-NF5280M5:/data/nt $ ~/software/blast/ncbi-blast-2.12.0/bin/update_blastdb.pl nt
Connected to NCBI
Downloading nt (53 volumes) ...
Downloading nt.00.tar.gz...corrupt download, trying again.
Downloading nt.00.tar.gz...Timeout at /usr/share/perl/5.26/Net/FTP.pm line 583.

尝试了两次还是不行，于是决定放弃这种方法

而且对于这次下载的nt库来说，总计大小是150个G，解压之后会更大，使用wget进行下载的时候速度贼慢，所以还是不建议用这种方法进行传输。

wget -c https://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nt.gz

使用Aspera进行数据传输文件

下面我们介绍一种非常高效的的数据传输工具Aspera

1、安装以及解压

wget -c https://download.asperasoft.com/download/sw/connect/3.8.1/ibm-aspera-connect-3.8.1.161274-linux-g2.12-64.tar.gz
tar -zxvf  ibm-aspera-connect-3.8.1.161274-linux-g2.12-64.tar.gz
bash ibm-aspera-connect-3.8.1.161274-linux-g2.12-64.sh

Note: 有时候sh命令对于有些sh脚本不一定适效，当sh不管用的时候就将其替换bash再试一下，简而言之，bash是sh的超集，具有更优雅的语法和更多的功能。几乎在所有情况下使用bash

sh 脚本运行完了之后，整个安装过程就完成了，其会在~目录下生成.aspera目录，其所有文件均在这个文件夹下。我们这次的是在/home/qianwj/.aspera/目录下

使用的方式也比较简单，而且非常方便的一点是这种传输工具支持断点传输

/home/qianwj/.aspera/connect/bin/ascp -k 1 -i /home/qianwj/.aspera/connect/etc/asperaweb_id_dsa.openssh -QT -l 500m anonftp@ftp.ncbi.nlm.nih.gov:blast/db/FASTA/nt.gz ./

注意 -i 表示引入私钥，一定要从绝对路径引入私钥asperaweb_id_openssh,-T最好加上，取消加密。否则很容易下载失败。-k表示断点续传，通常设置为1即可。-l 500m表示设置最大下载速度为每秒500Mb。在此设置下载下，平均速可达400Mb/s。

当然如果要下载的文件比较大的话，会经常出现中断的情况，比如

(base) qianwj@ubuntu-NF5280M5:/data/nt $ /home/qianwj/.aspera/connect/bin/ascp -k1  -i /home/qianwj/.aspera/connect/etc/asperaweb_id_dsa.openssh --overwrite=diff -QTr -l200m anonftp@ftp.ncbi.nlm.nih.gov:blast/db/FASTA/nt.gz ./
nt.gz                                                                                                                                                            24%   36GB  194Mb/s  1:26:30 ETA
Session Stop  (Error: Session shutdown failed, Unable to send channel data)

这个时候因为选了断点上传的功能，于是只需要重新执行一下之前的代码就好，非常方便好用

那么接下来我们就等待下载结束了

nt.gz全部下载结束之后我们需要对其进行解压缩

gunzip -c nt.gz > nt

解压完成之后就开始建立索引，目录是在/data/nt下

makeblastdb -in nt -dbtype nucl -title make_nt -parse_seqids -out /data/nt/make_nt/mk_nt -logfile make_nt.log

等这边建立索引结束之后就开始就行比对操作了（执行代码）

(roary) qianwj@ubuntu-NF5280M5:/data/nt/maknnt $ blastn -query ~/project/ONT/multi_sample_basecalling_gpu_sup/3_Error_correction/barcode03/barcode03_racon_round1.fasta -db mk_nt  -out /data/nt/barcode03_nt_blast/blastn_result.xls -evalue 1e-6  -outfmt 0

最后等待结果

另外补充一下对于blastn的使用参数

*** Formatting options
 -outfmt <String>
   alignment view options:
     0 = pairwise,
     1 = query-anchored showing identities,
     2 = query-anchored no identities,
     3 = flat query-anchored, show identities,
     4 = flat query-anchored, no identities,
     5 = XML Blast output,
     6 = tabular,
     7 = tabular with comment lines,
     8 = Text ASN.1,
     9 = Binary ASN.1,
    10 = Comma-separated values,
    11 = BLAST archive format (ASN.1) 

   Options 6, 7, and 10 can be additionally configured to produce
   a custom format specified by space delimited format specifiers.

balst的参考操作
https://www.cnblogs.com/leezx/p/6425620.html

这条路一直显示走不通，一直报错，显示不对，于是准备再次试试使用自带的`update_blastdb.pl`去下载`nt`库，参考链接：https://genehub.wordpress.com/page/8/

以及NCBI blast官方文档： https://www.ncbi.nlm.nih.gov/books/NBK279690/

执行代码如下：（需要conda activate roary）

nohup time update_blastdb.pl nt > nt_download.log &

然后就放到后台让其自动运行好了，到时候再来试一下

还是失败了，于是继续尝试其他的方法：https://www.jianshu.com/p/c4b6723fd786 这篇文档里面用的是比较传统的wget方法
下载地址： https://ftp.ncbi.nih.gov/blast/db/

成功下载！！开心，不知道这次的行不行~

还是失败了！！同样的原因，我怀疑可能是内存带不动导致的

于是试一下直接下载大肠杆菌的全部的数据再进行比对试一下

参考工具 https://github.com/kblin/ncbi-genome-download/blob/master/README.md 的用法

使用NCBI-genome-download 工具

直接安装

 pip install ncbi-genome-download

之后开始下载流程

ncbi-genome-download --assembly-levels complete bacteria

最后还是失败了

鉴于多次的失败，怀疑可能是比对的数据库太大导致内存不够进而运行失败
于是换个思路，blast还是继续照常，但是数据库换成自己之前下载构建的Ecoil全基因组序列，看看这次行不行得通

运行代码如下：

blastn -query ~/project/ONT/multi_sample_basecalling_gpu_sup/3_Error_correction/barcode03/barcode03_racon_round1.fasta -db ./Blast_ref/all_Escherichia_coli_complete_genome -out ./barcode03_blast/blast.out -evalue 1e-6  -outfmt 0

这次终于成功运行，对最后显示的blast.out结果进行查看

(blast) qianwj@ubuntu-NF5280M5:/data/barcode03_blast $ head blast.out  -n 100
BLASTN 2.12.0+


Reference: Zheng Zhang, Scott Schwartz, Lukas Wagner, and Webb
Miller (2000), "A greedy algorithm for aligning DNA sequences", J
Comput Biol 2000; 7(1-2):203-14.



Database: all_Escherichia_coli_complete_genome
           6,225 sequences; 9,814,017,588 total letters



Query= contig_1 LN:i:4473778 RC:i:47929 XC:f:1.000000

Length=4473778
                                                                      Score        E
Sequences producing significant alignments:                          (Bits)     Value

NZ_CP053607.1 Escherichia coli strain NEB5-alpha_F'Iq chromosome,...  2.398e+06  0.0  
NZ_CP017100.1 Escherichia coli strain K-12 NEB 5-alpha chromosome...  2.398e+06  0.0  
NZ_CP047127.1 Escherichia coli K-12 chromosome, complete genome       2.081e+06  0.0  
NZ_CP025520.1 Escherichia coli strain DH5alpha chromosome, comple...  1.957e+06  0.0  
NZ_CP045741.1 Escherichia coli strain DH5alpha chromosome, comple...  1.920e+06  0.0  
NC_017625.1 Escherichia coli DH1, complete sequence                   1.690e+06  0.0  
NZ_CP012125.1 Escherichia coli strain DH1Ec095 chromosome, comple...  1.584e+06  0.0  
NZ_CP016018.1 Escherichia coli strain ER1821R chromosome, complet...  1.584e+06  0.0  
NC_017638.1 Escherichia coli DH1, complete sequence                   1.584e+06  0.0  
NZ_CP042867.1 Escherichia coli strain ATCC BAA-196 chromosome, co...  1.584e+06  0.0  
NZ_CP069132.1 Escherichia coli BW25113 substr. CHXR1G20 chromosom...  1.584e+06  0.0  
NZ_CP064682.1 Escherichia coli K-12 strain BZKR2G40 chromosome, c...  1.584e+06  0.0  
NZ_CP064681.1 Escherichia coli K-12 strain BZKR3G40 chromosome, c...  1.584e+06  0.0

(blast) qianwj@ubuntu-NF5280M5:/data/barcode03_blast $ tail blast.out -n 100

Query  5354     CCAGGCATCAAATAAAACGAAAGGCTCAGTCGAAAGACTGGGCCTTTCGTTTTATCTGTT  5413
                ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  3548651  CCAGGCATCAAATAAAACGAAAGGCTCAGTCGAAAGACTGGGCCTTTCGTTTTATCTGTT  3548710

Query  5414     GTTTGTCGGTGAACGCTCTC  5433
                ||||||||||||||||||||
Sbjct  3548711  GTTTGTCGGTGAACGCTCTC  3548730


 Score = 148 bits (80),  Expect = 2e-31
 Identities = 80/80 (100%), Gaps = 0/80 (0%)
 Strand=Plus/Plus

Query  10912    CCAGGCATCAAATAAAACGAAAGGCTCAGTCGAAAGACTGGGCCTTTCGTTTTATCTGTT  10971
                ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  3548651  CCAGGCATCAAATAAAACGAAAGGCTCAGTCGAAAGACTGGGCCTTTCGTTTTATCTGTT  3548710

Query  10972    GTTTGTCGGTGAACGCTCTC  10991
                ||||||||||||||||||||
Sbjct  3548711  GTTTGTCGGTGAACGCTCTC  3548730


 Score = 148 bits (80),  Expect = 2e-31
 Identities = 80/80 (100%), Gaps = 0/80 (0%)
 Strand=Plus/Plus

Query  13695    CCAGGCATCAAATAAAACGAAAGGCTCAGTCGAAAGACTGGGCCTTTCGTTTTATCTGTT  13754
                ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  3548651  CCAGGCATCAAATAAAACGAAAGGCTCAGTCGAAAGACTGGGCCTTTCGTTTTATCTGTT  3548710

Query  13755    GTTTGTCGGTGAACGCTCTC  13774
                ||||||||||||||||||||
Sbjct  3548711  GTTTGTCGGTGAACGCTCTC  3548730

同样的操作方式，对于barcode01的数据也进行这样的比对

(blast) qianwj@ubuntu-NF5280M5:/data/barcode01_blast $ blastn -query ~/project/ONT/multi_sample_basecalling_gpu_sup/3_Error_correction/barcode01/racon/racon_round1.fasta -db ../Blast_ref/all_Escherichia_coli_complete_genome -out ./barcode01_Ecoil_blast.out  -evalue 1e-6  -outfmt 0