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

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:

• 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. – vkkodali Apr 8 at 13:22
• 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. – Kamil S Jaron Apr 8 at 14:49
• I posted the script here too... – Kamil S Jaron Apr 8 at 14:57

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)

• 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. – vkkodali Apr 9 at 11:57
• 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. – Kamil S Jaron Apr 9 at 15:12