NCBI的NT库比对——blastn
NCBI的NT库比对——blastn步骤一：NT库下载步骤二：随机提取10000条或5000条序列情况1：二代数据，取双端不配对数据。各自取5000条序列情况2：二代数据，取双端配对数据。取5000条序列
步骤三：对提取10000条序列进行全NT库blastn比对补充：staxid的scinames的关系信息步骤四：对blastn结果进行统计总结
NCBI的NT库比对——blastn
NT库比对，包括对测序数据和组装后的基因组序列进行NT库比对，查看是否存在菌污染以及是否是自己的数据。这里我提供这一部分的具体操作过程。
步骤一：NT库下载
前面有一篇文章NCBI下载nt/nr/swissprot库中介绍了下载NT库的全过程，查看完成即可。
步骤二：随机提取10000条或5000条序列
具体数目，自行决定 二代双端数据文件：FDSW220005290-2r_1.fq.gz和FDSW220005290-2r_2.fq.gz
情况1：二代数据，取双端不配对数据。各自取5000条序列
cd /home/zhaohuiyao/Genome_survey/02NT /home/zhaohuiyao/Biosoft/seqkit fx2tab -n …/01cleandata/FDSW220005290-2r_1.fq.gz | shuf -n 5000 > sample_1.name /home/zhaohuiyao/Biosoft/seqkit fx2tab -n …/01cleandata/FDSW220005290-2r_2.fq.gz | shuf -n 5000 > sample_2.name /home/zhaohuiyao/Biosoft/seqkit grep -n -f sample_1.name …/01cleandata/FDSW220005290-2r_1.fq.gz > sample_1.fq /home/zhaohuiyao/Biosoft/seqkit grep -n -f sample_2.name …/01cleandata/FDSW220005290-2r_2.fq.gz > sample_2.fq #将两个fastq文件整合并输出到fasta格式 cat sample_1.fq sample_2.fq > sample.fq /home/zhaohuiyao/Biosoft/general/seqkit fq2fa sample.fq > sample.fa
情况2：二代数据，取双端配对数据。取5000条序列
cd /home/zhaohuiyao/Genome_survey/02NT /home/zhaohuiyao/Biosoft/seqkit fx2tab -n …/01cleandata/FDSW220005290-2r_1.fq.gz | shuf -n 5000 > sample.name /home/zhaohuiyao/Biosoft/seqkit grep -n -f sample.name …/01cleandata/FDSW220005290-2r_1.fq.gz > sample_1.fq /home/zhaohuiyao/Biosoft/seqkit grep -n -f sample.name …/01cleandata/FDSW220005290-2r_2.fq.gz > sample_2.fq #将两个fastq文件整合并输出到fasta格式 cat sample_1.fq sample_2.fq > sample.fq /home/zhaohuiyao/Biosoft/general/seqkit fq2fa sample.fq > sample.fa
步骤三：对提取10000条序列进行全NT库blastn比对
保证软件blast已安装成功且测试通过 保证步骤一和步骤二完成
#全NT库已经是构建好的，直接进行blast
cd /home/zhaohuiyao/Genome_survey/02NT/
nohup blastn -num_threads 32 -max_target_seqs 10 -evalue 1e-05 -db /home/zhaohuiyao/Database/nt/nt/nt -outfmt "7 qseqid sseqid evalue pident ppos length mismatch gapopen qstart qend sstart send bitscore staxid sscinames stitle" -query ./XXX_10000.fa -out ./XXX_10000.blastn.out &
#具体的参数含义，自行查阅即可
#参数-max_target_seqs：最多比对序列数目，但是一条query和一条reference之间有多个比对情况，因此会出现多条比对情况
#可能会遇到这样的报错：BLAST Database error: Error: Not a valid version 4 database。原因可能是blast版本低，使用高版本的blast
#报警告，不影响blast结果，只是因为没有物种信息，无法用staxid获取sscinames信息
Warning: [blastn] Taxonomy name lookup from taxid requires installation of taxdb database with ftp://ftp.ncbi.nlm .nih.gov/blast/db/taxdb.tar.gz
#更新命令
cd /home/zhaohuiyao/Database/nt/
wget http://ftp.ncbi.nlm.nih.gov/blast/db/taxdb.tar.gz
tar -zxvf ./taxdb.tar.gz -C ./nt/
#更新后没有用，因此这里我整理一下staxid、scinames的信息，最后补充到blast结果中。如果你没有出现这个问题，则忽略
#我的blast结果（缺失scinames的信息，显示为N/A，倒数第二列）(可以依据倒数第三列Taxid，得到第二列信息)
# BLASTN 2.10.1+
# Query: A00808:1021:H5HHFDSX3:3:2478:32217:24627 1:N:0:GACCAAGCTT+AGTGGAGTCA
# Database: /home/zhaohuiyao/Database/NT/nt/nt
# Fields: query id, subject id, evalue, % identity, % positives, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, bit score, subject tax id, subject sci name, subject title
# 46 hits found
A00808:1021:H5HHFDSX3:3:2478:32217:24627 gi|2268683873|emb|OX155639.1| 2.56e-54 93.421 93.42 152 8 2 1 150 15472611 15472762 224 218720 N/A Carterocephalus palaemon genome assembly, chromosome: 2
补充：staxid的scinames的关系信息
这一步只有你的blastn的结果中scinames这一栏是N/A的情况下进行，若是没有，则跳过
#下载最新的taxdmp.zip（2022-09-28）
mkdir -p /home/zhaohuiyao/Database/taxdmp/ && cd /home/zhaohuiyao/Database/taxdmp/
wget https://ftp.ncbi.nih.gov/pub/taxonomy/new_taxdump/new_taxdump.zip
wget https://ftp.ncbi.nih.gov/pub/taxonomy/new_taxdump/new_taxdump.zip.md5
md5sum -c ./new_taxdump.zip.md5
unzip ./new_taxdump.zip
#解压后的文件中names.dmp，是我们需要的staxid-scinames关系文件
#执行脚本staxid-scinames.py，拿到两者关系文件staxid-scinames.txt，第一列：staxid；第二列：scinames
python3 ./staxid-scinames.py -i ./names.dmp -o ./
#staxid的scinames关系文件位置：/home/zhaohuiyao/Database/taxdmp/staxid-scinames.txt
步骤四：对blastn结果进行统计
这一步执行的前提的是你已经拿到了blastn的结果文件。因为我的脚本里面需要补充scinames信息，所以要提供staxid-scinames.txt文件，若你不需要，则需要将脚本略微修改。（两种我都会提供）
#对blast结果补充scinames信息，并进行总结分类
python3 blastn_stat.py -i ./XXX_10000.blastn.out -f /home/zhaohuiyao/Database/taxdmp/staxid-scinames.txt -o ./XXX_10000.blastn.out.stat -num 10000 -level all
#参数-num：指定进行Blastn的序列数目
#参数-level：两种情况。all：表示对所有比对结果进行统计。one：表示仅对一条比对结果进行统计。
#补充：blastn比对时，一个query序列，会对应多个比对结果的。在one情况下，我只看第一个比对结果，这个结果是什么物种，我就定义query序列是什么物种。在all情况下，会对所有比对结果进行统计，做一个字典。这个query序列可以有多个物种相对应。
#结果文件如下，cat XXX_10000.blastn.out.stat
Species_name mapping_reads total_reads ratio
unknown 9265 10000 92.65%
Wolbachia pipientis 207 10000 2.07%
Wolbachia endosymbiont of Ostrinia scapulalis 99 10000 0.99%
Wolbachia endosymbiont of Chrysomya megacephala 92 10000 0.92%
Opisthograptis luteolata 91 10000 0.91%
Wolbachia endosymbiont of Ostrinia furnacalis 90 10000 0.90%
Wolbachia endosymbiont of Corcyra cephalonica 76 10000 0.76%
Wolbachia pipientis wAlbB 76 10000 0.76%
Wolbachia endosymbiont of Drosophila mauritiana 75 10000 0.75%
Wolbachia endosymbiont of Diaphorina citri 73 10000 0.73%
#取blast前10的物种输出。第一列：物种名；第二列：比对到该物种的reads数；第三列：总参与比对reads数；第四列：占比
总结
因为考虑到全库比对会比较费时间，所以计划构建子库。但是尝试后发现，全库比对的时间可以接受，所以放弃制作子库，但是可以作为拓展练习。（正在努力，完成后会和大家分享学习~~ 这里提供别人的一篇文章，供大家参考，这是一篇nr子库构建：https://cloud.tencent.com/developer/article/1943973）
还有上面提到的四个脚本，我放在了gitee上，大家可以下载学习，也可以自己写。https://gitee.com/zhao-huiyao/python_script/tree/master/NT%E5%BA%93%E6%AF%94%E5%AF%B9