Grab Only Genes from Trinity Assembly FASTA Output

[eagb]

I'm doing some transcriptome assembly work with Trinity, which outputs a FASTA file full with ID numbers that look something like this: c0000_g1. Only the c0000 part is the ID number, the g bit refers to other transcripts for that gene. I want to go through and grab only the longest sequences per ID number and save that file as a new FASTA. I've been trying to use dictionaries with Python/Biopython but having trouble (new to using python for scientific computing). Here's a made up data set to explain what I'm talking about: https://gist.github.com/anonymous/69b46ddd94a91d30555d. So here, my new FASTA file would contain c00001_g1, c00002_g2, c00003_g1, and c00004_g2 and their respective sequences. Any ideas/code/scripts that would be of use here?

[bio_boris]

Hi Eagb,

I hacked something together in perl, it should work for smaller datasets, but it loads the entire fasta into memory, instead of just sequence lengths.

https://gist.github.com/bio-boris/91d7e627af46b451d8be

If two sequences are the same length, the second one would be thrown away.

[Apexy]

Hi eagb,

Is this exercise an attempt to reduce redundancy or complexity of the transcriptome assembly? Out of curiosity, what are the odds that the largest sequence per 'gene' is a chimera or that the 'c0000 part' is a collection of paralogs? The last time I checked, the largest sequence per 'c0000 part' is likely the one that is highly expressed. This stems from the behaviour of 'Inchworm'. Perhaps you are interested in achieving something other than reduction of redundancy but I thought I let you know of the possibility of loosing genuine transfrags.

[eagb]

Reduce redundancy. For the record, this is the solution I hacked together https://gist.github.com/ethanabaker1...38a614ada09e85 .

Probably not going to actually end up using these results - adopting a different analysis pipeline...