ncbi_nr数据库blast
1. 数据库下载
# 下载nr数据库
# 原始fasta文件
https://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/
# 构建好的db数据库（网络不好的话推荐下载这个，因为已经分成了很多小份，方便重新下载）
# 推荐下载这个，因为还省去了重新makeblastdb的步骤
https://ftp.ncbi.nlm.nih.gov/blast/db/
# 我这里下载了原始fasta文件
wget https://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nt.gz # 311 GB
wget https://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nt.gz.md5
remind.sh nt.download.end
cd /mnt/caigui/77_QHT_clear_assemble/7_blast_ntnr/taxid
# 下载分类名称id
# https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/accession2taxid/
wget -c https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/accession2taxid/nucl_gb.accession2taxid.gz # 2.2 GB
# accession accession.version taxid gi
# A00001 A00001.1 10641 58418
# A00002 A00002.1 9913 2
# A00003 A00003.1 9913 3
# A00004 A00004.1 32630 57971
wget -c https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz # 61 MB
# 需要其中的names.dmp文件
# 1 | all | | synonym |
# 1 | root | | scientific name |
# 2 | Bacteria | Bacteria <bacteria> | scientific name |
# 2 | bacteria | | blast name |
# 2 | eubacteria | | genbank common name |
# 2 | Monera | Monera <bacteria> | in-part |
# 2 | Procaryotae | Procaryotae <bacteria> | in-part |
# 2 | Prokaryotae | Prokaryotae <bacteria> | in-part |
# 比如这一行，分别是ACCESSION, VERSION, taxon, gi(好像是在网址上体现)
A80507 A80507.1 39946 6731332
# 对应https://www.ncbi.nlm.nih.gov/nuccore/6731332
# 39946就可以对应到Oryza sativa Indica Group
grep ^39946$'\t' names.dmp | grep "scientific name"
# 39946 | Indian rice | | common name |
# 39946 | Indica rice | | common name |
# 39946 | long-grained rice | | genbank common name |
# 39946 | Oryza sativa (indica cultivar-group) | | synonym |
# 39946 | Oryza sativa Indica Group | | scientific name |
# 39946 | Oryza sativa (indica group) | | synonym |
# 39946 | Oryza sativa subsp. indica Kato | | authority |
# 39946 | Oryza sativa subsp. indica | | synonym |
# 39946 | Oryza sp. Poi-6 | | includes |
grep ^39946$'\t' names.dmp | grep "scientific name"
2. 数据库问题
在下载了两天的nt.gz后，明明硬盘还剩1.1T的空间，它解压的时候还是报告"gzip: nt: No space left on device"
因此我只能切换下载建立好索引的nt数据库。
cd /mnt/caigui/77_QHT_clear_assemble/7_blast_ntnr/wget
vim alldb.txt #将https://ftp.ncbi.nlm.nih.gov/blast/db/这一页的内容全部复制上来
cat alldb.txt |grep "^nt\.[0-9]" | awk '{print $1}' > download.txt
cat alldb.txt |grep "^nt\.[0-9]" | awk '{print $4}' | grep G | sed 's/G//g' | awk '{sum+=$1} END {print sum}' # 总共259G，估计得下载个一天半吧 # 实际下载时间，一晚上，23小时
for i in `cat download.txt`
do
wget -c "https://ftp.ncbi.nlm.nih.gov/blast/db/${i}" & #引入变量时要用引号括起来
done
wait
remind.sh nt.db.download.end
md5sum -c *.md5
for i in `cat download.txt | grep "gz$"`
do
tar -zxvf $i -C ./ntdb/
echo "$i done"
done
wait
remind.sh nt.db.tar.end
3. 比对测试
cd /mnt/caigui/77_QHT_clear_assemble/7_blast_ntnr/blast
nohup blastn -query test100bp.modif.fa -db ../wget/ntdb/nt -evalue 1e-5 -outfmt 6 -num_threads 40 -num_alignments 1 > test100bp.result.only1.txt &
cat test100bp.result.only1.txt
# ptg000005l:35000-35500 NC_077626.1 90.683 161 3 6 1 161 237448 237596 2.11e-47 204
# ptg000005l:72500-73000 LC706752.1 100.000 32 0 0 89 120 341465 341496 4.84e-04 60.2
# ptg000005l:63000-63500 LC706752.1 97.967 492 1 1 18 500 175070 174579 0.0 845
# ptg000005l:44500-45000 LC706752.1 92.621 515 3 4 1 500 221762 222256 0.0 708
# ptg000005l:52000-52500 LC706752.1 92.969 512 1 10 1 500 60084 59596 0.0 713
# ptg000005l:52000-52500 LC706752.1 92.969 512 1 10 1 500 376914 377402 0.0 713
# ptg000005l:51500-52000 LC706752.1 94.106 509 7 8 1 500 60579 60085 0.0 752
# ptg000005l:51500-52000 LC706752.1 94.106 509 7 8 1 500 376419 376913 0.0 752
4. blast结果注释
- 对blast结果的第一列去重复，取第二列。（虽然设置了
-num_alignments 1，但如果是比对到一条序列的不同位置的话，还是会出现多行） - 在nucl_gb.accession2taxid中根据第二列，取第三列的第一个。（有时候第三列会出现两行，可以取第一行来避免后续代码出错）
- 根据上一行取到的结果，在names.dmp中查找后，在查找"scientific name"行，取第三列作为最终注释的名称。