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:

  • $\begingroup$ By gene sequences, you mean the entire sequence including the introns? I have been in similar situations before and my approach was to create a BED file from the GFF3 using some set of awk commands followed by bedtools getfasta for sequence. I am curious what your custom script solution was. $\endgroup$
    – vkkodali
    Apr 8, 2020 at 13:22
  • $\begingroup$ I made an edit, I mean to have one transcript per gene. So I suppose the longest transcript would be the most appropriate term. I just want to avoid multiple isoforms from the same gene. $\endgroup$ Apr 8, 2020 at 14:49
  • $\begingroup$ I posted the script here too... $\endgroup$ Apr 8, 2020 at 14:57

1 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]

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'
  • $\begingroup$ Do you know which one of the 3 steps takes the longest in your case? UCSC Kent Utils have some tools that you can use as substitution. One change I would suggest is that when you are constructing the genes2transcript_seuqneces dictionary, just keep the longest transcript and nothing else. That way, you are not storing the entire transcripts.fasta in memory. $\endgroup$
    – vkkodali
    Apr 9, 2020 at 11:57
  • $\begingroup$ Not sure which of the steps takes the longest, but all of this takes only a few seconds, so I don't think it's really worth any optimisation effort. I will take a look at UCSC Kent Utils. $\endgroup$ Apr 9, 2020 at 15:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.