<p>Here is the custom python script solution:</p>

<pre><code>#!/usr/bin/env python3

from collections import defaultdict
from Bio import SeqIO

# read file into dictionary of transcrpt-&gt;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-&gt;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) &gt; 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 = '&gt;' + gene + '\n' + str(transcripts[lengest].upper().seq) + '\n'
        out_file.write(fasta_string)
</code></pre>

CC BY-SA 4.0
2020-04-09T11:57:21.410Z
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
CC BY-SA 4.0
2020-04-09T15:12:56.170Z
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
CC BY-SA 4.0
2020-04-08T14:57:21.883Z
2020-04-08T14:57:21.883Z
Kamil S Jaron
<p>I got a genome (<code>data/genome/genome.fasta</code>) and braker-based genome annotation (<code>data/genome/annotation.gff3</code>), now I would like to get sequences of one transcript per gene (let's say the longest one).</p>

<p>Cufflinks has a <code>gffread</code> program to get the transcriptome out of the genome and the annotation. Something like</p>

<pre><code>gffread data/genome/annotation.gff3 -g data/genome/genome.fasta -w data/genome/transcripts.fasta -y data/genome/proteins.faa
</code></pre>

<p>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 <code>-M -K</code> I get still <code>.t1</code> and <code>.t2</code> sequences for some genes.</p>

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

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

<p>Related, not exactly identical: </p>

<ul>
<li><a href="https://bioinformatics.stackexchange.com/questions/5585/how-to-get-the-longest-fasta-sequence-including-all-possible-switching-isoforms">How to get the longest fasta sequence including all possible switching isoforms of a gene out of isoforms</a> (dealing more with terminology of Trinity)</li>
<li>When I used Kallisto/Salmon for differential expression analysis, I collapsed transcripts into genes while importing to R via <code>tximport</code>.</li>
</ul>

CC BY-SA 4.0
2020-04-08T13:22:38.363Z
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
CC BY-SA 4.0
2020-04-08T14:49:11.530Z
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
CC BY-SA 4.0
2020-04-08T14:57:39.153Z
I posted the script here too...
Kamil S Jaron
<genome><sequence-annotation><isoform>
Getting sequences of one transcript per gene out of annotation and genome
# Getting sequences of one transcript per gene out of annotation and genome

(Asked by: username_0 on 2020-04-08)

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.

**Tags:** genome, sequence-annotation, isoform
> **Comments:**
> **username_1**: 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.  
> **username_0**: 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.  
> **username_0**: I posted the script here too...  

---

## 1 answers

### Answer by username_0 on 2020-04-08. Score: 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)

> **Comments:**
> **username_1**: 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.  
> **username_0**: 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.  

---

