NCBI的NT库比对 |
您所在的位置：网站首页 › blast序列比对结果 › NCBI的NT库比对 |
NCBI的NT库比对
|
NCBI的NT库比对——blastn
NCBI的NT库比对——blastn步骤一：NT库下载步骤二：随机提取10000条或5000条序列情况1：二代数据，取双端不配对数据。各自取5000条序列情况2：二代数据，取双端配对数据。取5000条序列
步骤三：对提取10000条序列进行全NT库blastn比对补充：staxid的scinames的关系信息步骤四：对blastn结果进行统计总结
NCBI的NT库比对——blastn
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 |
【本文地址】
今日新闻 |
推荐新闻 |