我基于tblastn点击提取核苷酸序列的python代码也是s

Question

我使用独立的blast从python脚本在本地数据库上执行了多个查询tblastn搜索，结果是一个大的XML文件。我需要解析只接受前4次点击的输出，因为有些查询返回的结果超过4次。但是，由于输出是作为氨基酸序列比对记录给出的，因此对于文件中的每个tblastn记录，解析原始核苷酸序列中hsp开始和结束的命中主题的高分对是很重要的。在

所以我有了这段代码，但它的速度非常慢，而且考虑到数据量，它甚至需要一个多月才能完成它的工作。有人能帮我改进一下吗？在

> from Bio import SeqIO 
> from Bio.Blast import NCBIXML
> 
> infile_path = '/home/edson/ungulate/ungulate.fa'    # this is a file
> which contain unaligned nucleotide sequences  outfile_path =
> '/home/edson/ungulate/tblastn_result.fa'
> 
> for seq_record in SeqIO.parse(infile_path, 'fasta'):
>     flag = seq_record.description            # a flag is sequence identifier in a fasta file format

      with open(outfile_path, 'a') as outfile:
          with open('/home/edson/ungulate/tblastn_result.xml') as tblastn_file:
             tblastn_records = NCBIXML.parse(tblastn_file)
             for tblastn_record in tblastn_records:
                 for alignment in tblastn_record.alignments[:4]:
                     for hsp in alignment.hsps:
                         if flag in alignment.title:         
                          # this cross check if sequence identifier is present in an XML file
>                            sub_record = seq_record.seq[hsp.sbjct_start:hsp.sbjct_end] 
                             # this takes sequences in an infile path and slice them based on tblastn output
>                            outfile.write('>' + seq_record.description + '\n')
>                            outfile.write(str(sub_record + '\n'))

谢谢。在

Answer 1

至少有两个明显的瓶颈-对于外部循环的每个迭代，您

1. 重新打开outfile
2. 重新打开并重新解析tblastn_file

只需将这些操作移到外循环之外，就可以显著提高性能（当然，如果有多个外循环迭代）。在

另一个可能的改进：在每次迭代̀alignment.hsps上测试{}。对于相同的hsps，这个测试对于所有hsps都是常量，所以最好放在前面，即：

for alignment in tblastn_record.alignments[:4]:
    if flag in alignment.title:  
        for hsp in alignment.hsps:
           # etc...