Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-04-09T11:57:21.410Z
Answers.Text: 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.
Answers.UserDisplayName: vkkodali
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-04-09T15:12:56.170Z
Answers.Text: 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.
Answers.UserDisplayName: Kamil S Jaron
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-04-08T14:57:21.883Z
Answers.LastActivityDate: 2020-04-08T14:57:21.883Z
Answers.UserDisplayName: Kamil S Jaron
Body: <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>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-04-08T13:22:38.363Z
Comments.Text: 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.
Comments.UserDisplayName: vkkodali
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-04-08T14:49:11.530Z
Comments.Text: 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.
Comments.UserDisplayName: Kamil S Jaron
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-04-08T14:57:39.153Z
Comments.Text: I posted the script here too...
Comments.UserDisplayName: Kamil S Jaron
Tags: <genome><sequence-annotation><isoform>
Title: Getting sequences of one transcript per gene out of annotation and genome
ThreadText: # 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](https://bioinformatics.stackexchange.com/questions/5585/how-to-get-the-longest-fasta-sequence-including-all-possible-switching-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.  
  
  ---
  
  