Python. 用BioPython将genbank文件中三个最长的基因核苷酸序列排序并写入fasta文件

Question

我对Python还比较陌生，所以请原谅我这个问题的愚蠢之处。我有一个genbank文件，写了一段代码，可以找出最长的三个基因，并把它们放到一个新生成的fasta文件里。

from Bio import SeqIO
file="sequence.gb"
output=open("Top3.faa", "w")
record=SeqIO.parse(file, "genbank")
rec=next(record)
print('The genes with the top 3 longest lengths have beens saved in Top3.faa')
for f in rec.features:
        end=f.location.end.position
        start=f.location.start.position
        length=end-start
        bug=(rec.seq)
        if f.type=='CDS':
            if 'gene' in f.qualifiers:
                        if length>7000:
                                geneName=f.qualifiers['gene']
                                name=str(geneName)
                                lenth=str(length)
                                seq=str(bug[start:end])
                                output.write('>')
                                output.write(lenth)
                                output.write('\n')
                                output.write(seq)
                                output.write('\n')
output.close()

我想做的是，不再手动检查基因长度是否超过7kb，而是让代码自己完成这个检查，自动找到最长的三个基因。如果能给我一些建议，告诉我该往哪个方向去做，我会非常感激。谢谢！

Answer 1

你可以保存一个包含N个最大元素及其大小的列表。

可以像这样做（虽然我不能测试，所以可能会出错，但这个思路是对的）：

from Bio import SeqIO
file="sequence.gb"
output=open("Top3.faa", "w")
record=SeqIO.parse(file, "genbank")
rec=next(record)
print('The genes with the top 3 longest lengths have beens saved in Top3.faa')

# Largest genes and their size, sorted from the shortest to the longest.
# size first, gene name next, then seq.
largest_genes = [ (0, None, None) ] * 3;  # initialize with the number of genes you need.
for f in rec.features:
  end = f.location.end.position
  start = f.location.start.position
  length = end-start
  bug = (rec.seq)
  if f.type=='CDS' and 'gene' in f.qualifiers:
    if length > largest_genes[0][0]:  # [0] gives the first, [0] gives the length.
      # This one is larger than the smallest one we have.
      largest_genes = largest_genes[1:]  # drop the smallest one.
      # add this one
      largest_genes.append((length, f.qualifiers['gene'], str(bug[start:end])))  
      largest_genes.sort()  # re-sort.

for length, name, seq in largest_genes:   
  # name is not used but available.
  output.write('>')
  output.write(str(lenth))
  output.write('\n')
  output.write(seq)
  output.write('\n')
output.close()