0

I am relatively new to python so please forgive the idiocies that come with this question. I have a genbank file and have written a piece of code that will take the top 3 longest genes and place them into a newly generated fasta file.

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

What i'm trying to do is instead of manually imputing a check if it's over 7kb to find a way of the code do that itself and find the 3 top hits automatically. Any sort of help with direction of where i could go with this would be much appreciated. Thanks

Peepo
  • 3
  • 2

1 Answers1

1

You can keep the list of the N biggest (with their sizes).

Something like this (this might crash as I can't test it but the idea is there:

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()
Nicolas Defranoux
  • 2,646
  • 1
  • 10
  • 13
  • This is spot on what i was needing. It works great apart from matching it up with the name but hopefully i can sort that! Thanks for the help – Peepo Mar 02 '14 at 21:20
  • I did not include the name in my example because you were not outputing the name to the file... I can add it if you need it. – Nicolas Defranoux Mar 02 '14 at 22:00
  • Sorry about that...new to this whole site. It would be much appreciated if you could help me out with that one thank you – Peepo Mar 02 '14 at 22:03
  • No problem. I just added the name in my response, let me know if it does not work, I might have missed something ;) – Nicolas Defranoux Mar 02 '14 at 22:03