Question: Hi guys,

I have a couple thousands sequences that I translated using getorf from EMBOSS. I now need to select the the longest sequence out of many sequences. How can we do this? For example I have translated these two sets of contigs:

Contig: 01003765.1 Contig562 mRNA sequence

and

Contig: 5553765.1_1 Contig562 mRNA sequence

```
>01003765.1_1 [227 - 259] Contig562 mRNA sequence
ATGTCCTCCAGAAATGTTCGACGTTATTTTGGA
>01003765.1_2 [259 - 291] Contig562 mRNA sequence
ATGAGAATCAACTGGACGATGCCTGTGAACATA
>01003765.1_3 [427 - 459] Contig562 mRNA sequence
ATGTCTACGGATACAGTCACGATTCTAGAGGTA
>01003765.1_4 [201 - 476] Contig562 mRNA sequence
ATGGTAGCTGCCGACAAATTGGCTCAATGTCCTCCAGAAATGTTCGACGTTATTTTGGAT
GAGAATCAACTGGACGATGCCTGTGAACATATAGCGGAATATCTGGAGGCATATTGGAGA
GCTACTCATCCAGAAATTGTAACGAGTACAACACGACAGATCGGTAGTCCACCTCAAGCA
TCGCCTAGTGGAGACATGGGAGAAACGACGCTTCCAGCTCAGCACGATGTCTACGGATAC
AGTCACGATTCTAGAGGTATAAATTCAGGGTTCAGC
>5553765.1_1 [227 - 259] Contig562 mRNA sequence
ATGTCCTCCAGAAATGTTCGACGTTATTTTGGACAACTGGACGATGCCTGTGAACATATAGCGGAATATCTGGAGGCATATTGGAGAGCTACTCATCCAGAAATTGTAACGAGTACAACACGACAGATCGGTAGTCCACCTCAAGCATCGCCTAGTGGAGACATGGGAGAAACGACGCTTCCAGCTCAGCACGATGTCTACGGATACAGTCACGATTCTAGAGGTATAAATTCAGGGTTCAGC
>5553765.1_2 [259 - 291] Contig562 mRNA sequence
ATGAGAATCAACTGGACGATGCCTGTGAACATA
>5553765.1_3 [427 - 459] Contig562 mRNA sequence
ATGTCTACGGATACAGTCACGATTCTAGAGGTA
>5553765.1_4 [201 - 476] Contig562 mRNA sequence
ATGGTAGCTGCCGACAAATTGGCTCAATGTCCTCCAGAAATGTTCGACGTTATTTTGGAT
GAGAAT
```

As a result, I only want the longest sequence for each contig translated in the same file:

```
>01003765.1_4 [201 - 476] Contig562 mRNA sequence
ATGGTAGCTGCCGACAAATTGGCTCAATGTCCTCCAGAAATGTTCGACGTTATTTTGGAT
GAGAATCAACTGGACGATGCCTGTGAACATATAGCGGAATATCTGGAGGCATATTGGAGA
GCTACTCATCCAGAAATTGTAACGAGTACAACACGACAGATCGGTAGTCCACCTCAAGCA
TCGCCTAGTGGAGACATGGGAGAAACGACGCTTCCAGCTCAGCACGATGTCTACGGATAC
AGTCACGATTCTAGAGGTATAAATTCAGGGTTCAGC

>5553765.1_1 [227 - 259] Contig562 mRNA sequence
ATGTCCTCCAGAAATGTTCGACGTTATTTTGGACAACTGGACGATGCCTGTGAACATATAGCGGAATATCTGGAGGCATATTGGAGAGCTACTCATCCAGAAATTGTAACGAGTACAACACGACAGATCGGTAGTCCACCTCAAGCATCGCCTAGTGGAGACATGGGAGAAACGACGCTTCCAGCTCAGCACGATGTCTACGGATACAGTCACGATTCTAGAGGTATAAATTCAGGGTTCAGC
```

I would really appreciate your help if you could provide me the solution to this problem. Thank you.

Answer: Here is a dirty python script you can try to adapt. I haven't rigorously tested it, so spot check the results and make sure it is correct. This script assumes the output sequences are in blocks, meaning all possible ORFs of a transcript are grouped together in the fasta file. You need biopython.

```
import sys
from Bio import SeqIO

currentCid = ''
buffer = []

for record in SeqIO.parse(open(sys.argv[1]),"fasta"):
cid = str(record.description).split('.')[0][1:]

if currentCid == '':
currentCid = cid
else:
if cid != currentCid:
buffer.sort(key = lambda x : len(x[1]))
print '>' + buffer[-1][0]
print buffer[-1][1]
currentCid = cid
buffer = [(str(record.description),str(record.seq))]
else:
buffer.append((str(record.description),str(record.seq)))

buffer.sort(key = lambda x : len(x[1]))
print '>' + buffer[-1][0]
print buffer[-1][1]
```

Save it and use it by:

python longestORF.py yourORF.fasta > output.fasta