Getting sequences of one transcript per gene out of annotation and genome

Question

I got a genome (data/genome/genome.fasta) and braker-based genome annotation (data/genome/annotation.gff3), now I would like to get sequences of one transcript per gene (let's say the longest one).

Cufflinks has a gffread program to get the transcriptome out of the genome and the annotation. Something like

gffread data/genome/annotation.gff3 -g data/genome/genome.fasta -w data/genome/transcripts.fasta -y data/genome/proteins.faa

Will do the job. However, it is pulling multiple transcripts per gene. I read the help page twice, found parameters that could potentially do what I would like to (collapsing alternative transcripts into single one), but even when I specify my candidate parameters -M -K I get still .t1 and .t2 sequences for some genes.

The task is surprisingly simple, and I got already a custom script solution, but it just feels wrong. There must be a better way!

What is a reasonable/elegant way to get nucleotide sequences of genes out of a genome and annotation? Am I missing something with cufflinks?

Related, not exactly identical:

• How to get the longest fasta sequence including all possible switching isoforms of a gene out of isoforms (dealing more with terminology of Trinity)
• When I used Kallisto/Salmon for differential expression analysis, I collapsed transcripts into genes while importing to R via tximport.

Answer 1

Here is the custom python script solution:

#!/usr/bin/env python3

from collections import defaultdict
from Bio import SeqIO

# read file into dictionary of transcrpt->gene
# jg25690   jg25690.t1
# jg25690   jg25690.t2
# jg25691   jg25691.t1
transcripts2genes = dict()
gene2transcript_mapfilename = 'data/genome/transcripts2genes.map'
with open(gene2transcript_mapfilename, 'r') as g2t_file:
    for line in g2t_file:
        gene, trainscript = line.rstrip('\n').split('\t')
        transcripts2genes[trainscript] = gene

# construct a dictionary gene->list of transcripts 
# this fasta file is output of default gffread -w
genes2transcript_seuqneces = defaultdict(list)
transcripts_fasta_filename = 'data/genome/transcripts.fasta'
transcripts_fasta_file = SeqIO.parse(transcripts_fasta_filename, 'fasta')
for seq_record in transcripts_fasta_file:
    gene = transcripts2genes[seq_record.name]
    genes2transcript_seuqneces[gene].append(seq_record)

genes_fasta_filename = 'data/genome/genes.fasta'
with open(genes_fasta_filename, 'w') as out_file:
    # iterate through genes
    for gene in genes2transcript_seuqneces.keys():
        # get all the associated transcripts
        transcripts = genes2transcript_seuqneces[gene]
        # find the longest, by default the one with index 0
        lengest = 0
        for index, tr in enumerate(transcripts):
            if len(transcripts[lengest].seq) > len(tr.seq):
                lengest = index
        # this would be nicer with Bio interface (with wrapping of sequence lines)
        # but I was lazy and blast does not care, only thing I do is converting sequence to uppercase (not even sure if necessary)
        fasta_string = '>' + gene + '\n' + str(transcripts[lengest].upper().seq) + '\n'
        out_file.write(fasta_string)