使用Diamond比对NR数据库获取物种注释_diamond nr

nohup wget -t 0 -c -b https://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nr.gz &
wget -c https://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nr.gz.md5
#检查下载的完整性（一样的话就是完整了）
md5sum nr.gz  ### 899ac219cac213c60fede9c3d9ef8f7b  nr.gz
cat nr.gz.md5 ### 899ac219cac213c60fede9c3d9ef8f7b  nr.gz
nohup gunzip  nr.gz  &
mv nr  nr.fa
 
###下载NR物种相关信息和taxid信息
wget -t 0 -c -b https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/accession2taxid/prot.accession2taxid.gz
wget -t 0 -c -b https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/accession2taxid/prot.accession2taxid.gz.md5
wget -t 0 -c -b https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz
wget -t 0 -c -b https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz.md5
#检查完整性
md5sum prot.accession2taxid.gz  #32b8ca7ea712e161c72af69135fc938e
cat prot.accession2taxid.gz.md5 #32b8ca7ea712e161c72af69135fc938e
md5sum taxdump.tar.gz           #26558b800bc1b3795c25d1f0ead65412
cat taxdump.tar.gz.md5          #26558b800bc1b3795c25d1f0ead65412  
#解压
tar -zxvf taxdump.tar.gz
gunzip prot.accession2taxid.gz
#移动需要的文件
cd ~
mkdir .taxonkit
cp *.dmp ~/.taxonkit

conda create -n NR_database_search
conda activate NR_database_search
conda install -c bioconda taxonkit=0.15.0
conda install -c bioconda diamond=2.1.8
conda install -c bioconda csvtk=0.28.0

nohup diamond makedb --threads 180 --in nr.fa --db NR_2023_07_23 &

diamond blastx --db NR_2023_07_23.dmnd --query nucleic_reads.fna\
-o nucleic_matches_fmt6.txt --threads 180 --evalue 0.00001 \
--max-target-seqs 5 --outfmt 6
 
diamond blastp --db NR_2023_07_23.dmnd --query protein_reads.fna\
-o protein_matches_fmt6.txt --threads 180 --evalue 0.00001 \
--max-target-seqs 5 --outfmt 6
## --outfmt 6 最好别改变

## 一些主要的物种编号
# 2     bacteria
# 2157  archaea
# 4751  fungi
# 10239 virus
# 2759 Eukaryota
 
#看taxnokit安好了么
taxonkit -h
 
#创建目录
cd /home/zhongpei/database/NR_2023_07_23
mkdir /home/zhongpei/database/NR_2023_07_23/NCBI_Main_tax
cd NCBI_Main_tax
 
#开始
taxonkit list -j 4 --ids 2,2157,4751,10239,2759 --indent "" > NCBI_Main.taxid.txt
# -j 是线程，软件说4个够了；--ids 是需要的物种编号，用逗号分隔
wc -l NCBI_Main.taxid.txt # 2708739 NCBI_Main.taxid.txt
head -n 5 NCBI_Main.taxid.txt #查看内容
 
# 提取taxid和taxonomy（界门纲目科属种）的对应信息到NCBI_Main.taxid.txt
less NCBI_Main.taxid.txt | taxonkit reformat -I 1 -r Unassigned -f "{k}\t{p}\t{c}\t{o}\t{f}\t{g}\t{s}\t{t}"| sed '1i\Taxid\tKingdom\tPhylum\tClass\tOrder\tFamily\tGenus\tSpecies\tStrain' > NCBI_Main.taxid_new.txt
# -I 1 一个制表符分隔；-r 没有找到的用什么字符去填充，这里用的“Unassigned” 
# -f 输出的格式；1i表示在第行之前插入文本(sed用法，不太会)

cat prot.accession2taxid | csvtk -t grep -f taxid -P NCBI_Main.taxid.txt | csvtk -t cut -f accession.version,taxid > NCBI_seqid_taxid.txt
# "cat prot.accession2taxid |" 把 prot.accession2taxid 的内容到下面的 csvtk
# -t 输入内容是制表符分隔；grep 这是csvtk的1个子命令，用于在文件中搜索匹配的行
# -P 搜索那些"taxid"字段的值出现在"NCBI_Main.taxid.txt"文件中的行
# cut 这是csvtk的1个子命令，用于从输入中选择特定的字段

#!/home/zhongpei/miniconda3/envs/py39/bin/python3.9
# ##########################################################
# match diamond blast NR result.txt and NCBI_seqid_taxid.txt
# written by PeiZhong in IFR of CAAS
 
import os
import argparse
 
parser = argparse.ArgumentParser(description='match diamond blast NR result.txt and NCBI_seqid_taxid.txt. '
'!!! all the file should end will .txt !!!')
parser.add_argument('diamond_result_folder_path',help='full Path of the folder that contain your diamond result txt')
parser.add_argument('result_files_mark',help='mark=The name of the mark specific to '
'your two-column diamond results file in this folder, e.g., clean, '
'for example,the mark of result file SY10_NR_diamond_clean_fmt6_taxid.txt is clean')
parser.add_argument('NCBI_acc_taxid',help='full Path that contain NCBI_seqid_taxid.txt')
args = parser.parse_args()
 
result_file_folder_path = args.diamond_result_folder_path
NCBI_file = args.NCBI_acc_taxid
file_mark = args.result_files_mark
 
os.chdir(result_file_folder_path)
files = os.listdir(result_file_folder_path)
print(files)
 
db={}
 
print("start db read")
f_table = open("%s" % (NCBI_file), 'r')
print("start db build")
for line in f_table.readlines():
    line=line.split('\t')
    acc_num = line[0].strip()
    tax_num = line[1].strip()
    db[acc_num] = tax_num
print("finish db build")
 
file_ls = []
for result_file in files:
    if file_mark in result_file:
        file_ls.append(result_file)
file_ls.sort()
print(file_ls)
 
header = "qseqid\taccession.version\tpident\tlength\tmismatch\tgapopen\tqstart\tqend\tsstart\tsend\tevalue\tbitscore"
 
for result_file in file_ls:
    print(result_file)
    f_result = open("%s" % (result_file), 'r+')
    content = f_result.read()
    f_result.seek(0, 0)
    f_result.write(header.rstrip('\r\n') + '\n' + content)
    f_result.close()
    out_name = str(result_file).strip('txt')
    out_name = out_name+'_taxid.txt'
    f_result = open("%s" % (result_file), 'r')
    f_out = open("%s" % (out_name), 'a')
    for line in f_result.readlines():
        line = line.split('\t')
        query_num = line[0].strip()
        acc_q_num = line[1].strip()
        if acc_q_num in db:
            print(query_num,end="\t",file=f_out)
            print(acc_q_num, end="\t", file=f_out)
            print(db[acc_q_num], file=f_out)
    f_out.close()
    f_result.close()
f_table.close()

chmod +x diamond_NR_tax.py(完整地址)
diamond_NR_tax.py(完整地址) -h
diamond_NR_tax.py(完整地址) diamond_result_folder_path result_files_mark NCBI_acc_taxid
diamond_result_folder_path：你存放上面处理完的diamond比对文件，txt结尾的文件，的目录（完整地址）
result_files_mark：这个地址中，你的这些文件独有的标识字符串
#for example,the mark of result file 'SY10_clean.txt' and 'SY11_clean.txt' is 'clean'
NCBI_acc_taxid：你的NCBI_seqid_taxid.txt文件的完整地址

#!/home/zhongpei/miniconda3/envs/py39/bin/python3.9
# ##########################################################
# match diamond blast NR taxid result.txt and NCBI_Main.taxid_new.txt
# written by PeiZhong in IFR of CAAS
 
import os
import argparse
 
parser = argparse.ArgumentParser(description='match diamond blast NR taxid result.txt and NCBI_Main.taxid_new.txt. '
'!!! all the file should end will .txt !!!')
parser.add_argument('diamond_result_folder_path',help='full Path of the folder that contain your result txt')
parser.add_argument('result_files_mark',help='mark=The name of the mark specific to '
'your two-column diamond results file in this folder, e.g., clean, '
'for example,the mark of result file SY10_NR_diamond_clean_fmt6_taxid.txt is clean')
parser.add_argument('NCBI_taxid_tax',help='full Path that contain NCBI_Main.taxid_new.txt')
args = parser.parse_args()
 
result_file_folder_path = args.diamond_result_folder_path
NCBI_file = args.NCBI_taxid_tax
file_mark = args.result_files_mark
 
os.chdir(result_file_folder_path)
files = os.listdir(result_file_folder_path)
print(files)
 
db={}
 
print("start db read")
f_table = open("%s" % (NCBI_file), 'r')
print("start db build")
for line in f_table.readlines():
    line=line.split('\t')
    taxid = line[0].strip()
    tax_anno = line[1:8]
    db[taxid] = tax_anno
print("finish db build")
 
file_ls = []
for result_file in files:
    if file_mark in result_file:
        file_ls.append(result_file)
file_ls.sort()
print(file_ls)
 
for result_file in file_ls:
    print(result_file)
    out_name = str(result_file).strip('txt')
    out_name = out_name+'_tax.txt'
    f_result = open("%s" % (result_file), 'r')
    f_out = open("%s" % (out_name), 'a')
    for line in f_result.readlines():
        line = line.split('\t')
        query_num = line[0].strip()
        acc_q_num = line[1].strip()
        taxid_1 = line[2].strip()
        if taxid_1 in db:
            print(query_num,end="\t",file=f_out)
            print(acc_q_num, end="\t", file=f_out)
            print(taxid_1, end="\t", file=f_out)
            tax_in_db = db[taxid_1]
            str_ls = map(str, tax_in_db)
            tax = '\t'.join(str_ls)
            print(tax, file=f_out)
    f_out.close()
    f_result.close()
f_table.close()