Is there a convenient way to extract the longest isoforms from a transcriptome fasta file? I had found some scripts on biostars but none are functional and I'm having difficulty getting them to work.I'm aware that the longest isoforms aren't necessarily 'the best' but it will suit my purposes.The fasta was generated via Augustus. Here is what the fasta file looks like currently (sequence shortened to save space)>Doug_NoIndex_L005_R1_001_contig_2.g7.t1atggggcataacatagagactggtgaacgtgctgaaattctacttcaaagtctacctgattcgtatgatcaactcatcattaatataaccaaaaacctagaaattctagccttcgatgatgttgcagctgcggttcttgaagaagaaagtcggcgcaagaacaaagaagatagaccg>Doug_NoIndex_L005_R1_001_contig_2.g7.t2atggggcataacatagagactggtgaacgtgctgaaattctacttcaaagtctacctgattcgtatgatcaactcatcaThe format is as such:Gene 1 isoform 1  Gene 1 isoform 2  Gene 2 isoform 1  Gene 2 isoform 2   and so forth. There are several genes that have more than one pair of isoforms (up to 3 or 4). There are roughly 80,000 total transcripts, probably 25,000 genes. I would like to extract the single longest isoform for each gene.
How can longest isoforms (per gene) be extracted from a FASTA file?
fasta;isoform;filtering
While the solution from https://bioinformatics.stackexchange.com/users/96/daniel-standage should work (after adjusting for possible python3 incompatibility), the following is a shorter and less memory demanding method that uses biopython:#!/usr/bin/env pythonfrom Bio import SeqIOimport syslastGene = Nonelongest = (None, None)for rec in SeqIO.parse(sys.argv[1], fasta):    gene = ..join(rec.id.split(.)[:-1])    l = len(rec)    if lastGene is not None:        if gene == lastGene:            if longest[0] < l:                longest = (l, rec)        else:            lastGene = gene            SeqIO.write(longest[1], sys.stdout, fasta)            longest = (l, rec)    else:        lastGene = gene        longest = (l, rec)SeqIO.write(longest[1], sys.stdout, fasta)If you saved this as filter.py, then filter.py original.fa > subset.fa would be the command to use.