Python. Попытка сортировки файла для 3 самых длинных нуклеотидных последовательностей гена из файла genbank в файл fasta с использованием BioPython

Question

Я относительно новичок в Python, поэтому, пожалуйста, прости идиотизм, который идет с этим вопросом. У меня есть файл genbank и я написал фрагмент кода, который возьмет 3 самых длинных гена и поместит их во вновь сгенерированный файл 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()

То, что я пытаюсь сделать, это вместо того, чтобы вручную вводить проверку, если она превышает 7 КБ, чтобы найти способ, которым код делает это сам, и автоматически находит 3 главных хита. Любая помощь с указанием того, куда я мог бы пойти с этим, была бы очень признательна. Спасибо

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()