0

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

我想要做的不是手动输入检查是否超过 7kb 来找到代码本身的方法并自动找到 3 个热门点击。任何有关我可以去哪里的帮助将不胜感激。谢谢

4

1 回答 1

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()
于 2014-03-02T21:00:04.830 回答