Answers.Body: <p>Late to the party here, but I like to try to avoid writing scripts when some command line magic will do. It's good practice to index your <code>FASTA</code> so use it.</p>
  
  <pre><code>samtools faidx &lt;myfasta.fa&gt;
  </code></pre>
  
  <p>A little bit of <code>awk</code> and <code>sort</code> can determine the largest isoform of each gene (and they don't even have to be sorted by name in the source <code>FASTA</code>).</p>
  
  <pre><code>awk -F'[\t.]' '{print $1,$2,$3,$4}' &lt;myfasta.fa&gt;.fai | sort -k4nr,4 | sort -uk1,2 | cut -f1-3 -d' '| tr ' ' '.' &gt; selection.ls
  </code></pre>
  
  <ul>
  <li>The <code>-F[\t.]</code> specifies both the <code>fai</code> tab delimiters, and the full stops in your contig names (<code>Doug_NoIndex_L005_R1_001_contig_2.g7.t1</code>) as delimiters for <code>awk</code>.</li>
  <li>For each line we print four columns. <code>$1</code> to <code>$3</code> are the three components of your contig name (after splitting on the <code>.</code>), and column $4 is the column in the index for sequence length.</li>
  <li>We use sort to organise the lines by <code>$4</code> (the contig length). Reverse sort with <code>-r</code> so the longest genes appear first.</li>
  <li>We then subsort by columns <code>$1</code> and <code>$2</code> (the contig and gene name (<code>g1</code>...)). <code>-u</code> then drops any repeat of column <code>$1</code> and <code>$2</code> pairs (that is, all lines after the first pair -- the longest length).</li>
  <li>Finally, we reassemble <code>$1</code> to <code>$3</code> with <code>cut</code> and <code>tr</code> to give you a list of sequence names.</li>
  </ul>
  
  <p>A loop with our old friend <code>samtools faidx</code> can put your chosen sequences out to a new <code>FASTA</code>.</p>
  
  <pre><code>while read contig;
      do samtools faidx &lt;myfasta.fa&gt; $contig &gt;&gt; selection.fa;
  done &lt; selection.ls
  </code></pre>
  
  <p>Don't forget to index it, they are useful!</p>
  
  <pre><code>samtools faidx selection.fa
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-13T00:51:52.510Z
Answers.Text: Oh and if you wanted to remove the .tN part, just switch out that `cut -f1-3 -d' '` for `cut -f1-2 -d' '`
Answers.UserDisplayName: Sam Nicholls
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-13T00:28:05.247Z
Answers.LastActivityDate: 2017-06-13T00:51:03.047Z
Answers.LastEditDate: 2017-06-13T00:51:03.047Z
Answers.UserDisplayName: Sam Nicholls
Answers.Body: <p>Here is a solution in R. Could get really slow with big files. Works for the example you posted.</p>
  
  <pre><code>library(Biostrings)
  
  ## read your fasta in as Biostrings object
  fasta.s &lt;- readDNAStringSet("sample.fa")
  
  ## get the read names (in your case it has the isoform info)
  names.fasta &lt;- names(fasta.s)
  
  ## extract only the relevant gene and isoform id (split name by the period symbol)
  gene.iso &lt;- sapply(names.fasta,function(j) cbind(unlist(strsplit(j,'\\.'))[2:3]))
  
  ## convert to good data.frame = transpose result from previous step and add relevant column names
  gene.iso.df &lt;- data.frame(t(gene.iso))
  colnames(gene.iso.df) &lt;- c('gene','isoform')
  
  ## and length of isoforms
  gene.iso.df$width &lt;- width(fasta.s)
  
  ## split data.frame into list with entry for each gene
  gene.iso.df.split &lt;- split(gene.iso.df,gene.iso.df$gene)
  
  ## optional to keep all the information but really you just need indices
  ##gene.iso.df.split.best &lt;- lapply(gene.iso.df.split,function(x) x[order(x$width)[1],])
  
  ## pull out the longest isoform ID for each gene (in case of a tie just take the first one)
  best.id &lt;- sapply(gene.iso.df.split,function(x) row.names(x)[order(x$width)[1]])
  
  ## subset your original reads with the subset
  fasta.s.best &lt;- fasta.s[best.id]
  
  ## export new fastafile containing longest isoform per gene
  writeXStringSet(fasta.s, filepath='sample_best_isoform.fasta')
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-09T10:51:42.197Z
Answers.LastActivityDate: 2017-06-09T10:51:42.197Z
Answers.UserDisplayName: story
Answers.Body: <p>This solution is will work for your example and probably all Augustus-derived Fastas, but mileage will vary beyond that.</p>
  
  <pre class="lang-python prettyprint-override"><code>#!/usr/bin/env python
  from __future__ import print_function
  import sys
  
  def parse_fasta(data):
      """Stolen shamelessly from http://stackoverflow.com/a/7655072/459780."""
      name, seq = None, []
      for line in data:
          line = line.rstrip()
          if line.startswith('&gt;'):
              if name:
                  yield (name, ''.join(seq))
              name, seq = line, []
          else:
              seq.append(line)
      if name:
          yield (name, ''.join(seq))
  
  isoforms = dict()
  for defline, sequence in parse_fasta(sys.stdin):
      geneid = '.'.join(defline[1:].split('.')[:-1])
      if geneid in isoforms:
          otherdefline, othersequence = isoforms[geneid]
          if len(sequence) &gt; len(othersequence):
              isoforms[geneid] = (defline, sequence)
      else:
          isoforms[geneid] = (defline, sequence)
  
  for defline, sequence in isoforms.values():
      print(defline, sequence, sep='\n')
  </code></pre>
  
  <p>Save as <code>longest.py</code>, and invoke on the command line like so.</p>
  
  <pre><code>python longest.py &lt; intput.fasta &gt; output.fasta
  </code></pre>
  
  <hr>
  
  <p>The line <code>geneid = '.'.join(defline[1:].split('.')[:-1])</code> is key. Let me break it down.</p>
  
  <ul>
  <li><code>defline[1:]</code>: ignore the first character of the defline (the <code>&gt;</code> symbol)</li>
  <li><code>defline[1:].split('.')</code>: split the string on <code>.</code> symbol</li>
  <li><code>defline[1:].split('.')[:-1]</code>: ignore the last value after the split (the isoform name)</li>
  <li><code>'.'.join(defline[1:].split('.')[:-1])</code>: join the split values again</li>
  </ul>
  
  <p>This value will be the same for all isoforms for the same gene, so we use it in the script to keep track of the longest sequence associated with each gene ID.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-08T22:10:46.897Z
Answers.Text: This could save a lot of memory and a bit of time if you write out and reset `isoforms` whenever `geneid not in isoforms`. OP says isoforms are always next to each other.
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-08T22:14:29.683Z
Answers.Text: Hi Daniel, I tried running this script and received the following error message:  
  
      File "longest.py", line 29
           print(defline, sequence, sep='\n') 
                                  ^  
      SyntaxError: invalid syntax
Answers.UserDisplayName: ZincFingers
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-08T22:19:30.163Z
Answers.Text: @ZincFingers: Change it to `print(defline)` and `print(sequence)` on separate lines.
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-08T22:24:54.510Z
Answers.Text: Python 2 / 3 incompatibility. Will update.
Answers.UserDisplayName: Daniel Standage
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-08T22:06:35.127Z
Answers.LastActivityDate: 2017-06-09T08:07:00.293Z
Answers.LastEditDate: 2017-06-09T08:07:00.293Z
Answers.UserDisplayName: Daniel Standage
Answers.Body: <p>An ex-coworker (Josep Avril) has written a couple of very useful little scripts that convert fasta to tbl (<code>seqID&lt;TAB&gt;Sequence</code>) and back again. These are extremely handy for this sort of thing (and are included at the end of this answer). Using them, you can convert your fasta to a one sequence per line format, keep the longest sequence with a simple awk script and convert back to fasta.</p>
  
  <p>I am assuming that the part after the last <code>.</code> in the ID line should be removed (so that <code>Doug_NoIndex_L005_R1_001_contig_2.g7.t1</code> and <code>Doug_NoIndex_L005_R1_001_contig_2.g7.t2</code> both map to <code>Doug_NoIndex_L005_R1_001_contig_2.g7</code>). If that is a correct assumption, this should work for you:</p>
  
  <pre><code>$ FastaToTbl file.fa | sed 's/\.[^.]*\t/\t/' | 
      awk -F"\t" '{
                      if(length($2)&gt;length(a[$1])){
                          a[$1]=$0
                      }
                  }
                  END{
                      for(i in a){
                          print a[i]
                      }
                  }' | TblToFasta
  &gt;Doug_NoIndex_L005_R1_001_contig_2.g7 
  atggggcataacatagagactggtgaacgtgctgaaattctacttcaaagtctacctgat
  tcgtatgatcaactcatcattaatataaccaaaaacctagaaattctagccttcgatgat
  gttgcagctgcggttcttgaagaagaaagtcggcgcaagaacaaagaagatagaccg
  </code></pre>
  
  <p>One possible issue with that approach is that it needs to keep one sequence per ID in memory. If you have a huge fasta file and not much memory, that might be a problem. If that's the case, you can first sort the file to ensure that similar IDs always appear together and then print each line as soon as you find the next ID:</p>
  
  <pre><code>FastaToTbl file.fa | sed 's/\.[^.]*\t/\t/' | LC_ALL=C sort -t$'\t' -k1,1 |  
      awk -F"\t" '{
                      if(NR&gt;1 &amp;&amp; prev!=$1){
                          print a[$1]; 
                          prev=$1
                      } 
                      if(length($2)&gt;length(a[$1])){
                          a[$1]=$0
                      }
                  }
                  END{
                      print a[$1]
                  }' | TblToFasta
  </code></pre>
  
  <hr>
  
  <ul>
  <li><p>FastaToTbl</p>
  
  <pre class="lang-c prettyprint-override"><code>#!/usr/bin/awk -f
  {
          if (substr($1,1,1)=="&gt;")
          if (NR&gt;1)
                      printf "\n%s\t", substr($0,2,length($0)-1)
          else 
              printf "%s\t", substr($0,2,length($0)-1)
          else 
                  printf "%s", $0
  }END{printf "\n"}
  </code></pre></li>
  <li><p>TblToFasta</p>
  
  <pre class="lang-c prettyprint-override"><code>#! /usr/bin/awk -f
  {
    sequence=$NF
  
    ls = length(sequence)
    is = 1
    fld  = 1
    while (fld &lt; NF)
    {
       if (fld == 1){printf "&gt;"}
       printf "%s " , $fld
  
       if (fld == NF-1)
        {
          printf "\n"
        }
        fld = fld+1
    }
    while (is &lt;= ls)
    {
      printf "%s\n", substr(sequence,is,60)
      is=is+60
    }
  }
  </code></pre></li>
  </ul>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-08T22:03:44.497Z
Answers.LastActivityDate: 2017-06-10T07:35:23.737Z
Answers.LastEditDate: 2017-06-10T07:35:23.737Z
Answers.UserDisplayName: terdon
Answers.Body: <p>While the solution from <a href="https://bioinformatics.stackexchange.com/users/96/daniel-standage">https://bioinformatics.stackexchange.com/users/96/daniel-standage</a> should work (after adjusting for possible python3 incompatibility), the following is a shorter and less memory demanding method that uses biopython:</p>
  
  <pre class="lang-python prettyprint-override"><code>#!/usr/bin/env python
  from Bio import SeqIO
  import sys
  
  lastGene = None
  longest = (None, None)
  for rec in SeqIO.parse(sys.argv[1], "fasta"):
      gene = ".".join(rec.id.split(".")[:-1])
      l = len(rec)
      if lastGene is not None:
          if gene == lastGene:
              if longest[0] &lt; l:
                  longest = (l, rec)
          else:
              lastGene = gene
              SeqIO.write(longest[1], sys.stdout, "fasta")
              longest = (l, rec)
      else:
          lastGene = gene
          longest = (l, rec)
  SeqIO.write(longest[1], sys.stdout, "fasta")
  </code></pre>
  
  <p>If you saved this as <code>filter.py</code>, then <code>filter.py original.fa &gt; subset.fa</code> would be the command to use.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-09T07:17:47.383Z
Answers.LastActivityDate: 2017-06-09T08:06:41.350Z
Answers.LastEditDate: 2017-06-09T08:06:41.350Z
Answers.UserDisplayName: Devon Ryan
Body: <p>Is there a convenient way to extract the longest isoforms from a transcriptome fasta file? I had found some <a href="https://www.biostars.org/p/107759/" rel="nofollow noreferrer">scripts on biostars</a> but none are functional and I'm having difficulty getting them to work.</p>
  
  <p>I'm aware that the longest isoforms aren't necessarily 'the best' but it will suit my purposes.</p>
  
  <p>The fasta was generated via Augustus. Here is what the fasta file looks like currently (sequence shortened to save space)</p>
  
  <pre><code>&gt;Doug_NoIndex_L005_R1_001_contig_2.g7.t1
  atggggcataacatagagactggtgaacgtgctgaaattctacttcaaagtctacctgattcgtatgatcaactcatca
  ttaatataaccaaaaacctagaaattctagccttcgatgatgttgcagctgcggttcttgaagaagaaagtcggcgcaagaacaaagaagatagaccg
  &gt;Doug_NoIndex_L005_R1_001_contig_2.g7.t2
  atggggcataacatagagactggtgaacgtgctgaaattctacttcaaagtctacctgattcgtatgatcaactcatca
  </code></pre>
  
  <p>The format is as such:</p>
  
  <pre><code>Gene 1 isoform 1  
  Gene 1 isoform 2  
  Gene 2 isoform 1  
  Gene 2 isoform 2   
  </code></pre>
  
  <p>and so forth. There are several genes that have more than one pair of isoforms (up to 3 or 4). There are roughly 80,000 total transcripts, probably 25,000 genes. I would like to extract the single longest isoform for each gene.</p>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-08T21:38:20.837Z
Comments.Text: It depends on the exact format of your fasta. Where did you get it? What do the fasta header lines look like? Can you paste a small example?
Comments.UserDisplayName: user172818
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-08T21:38:28.420Z
Comments.Text: >Doug_NoIndex_L005_R1_001_contig_4.g13.t1
  
  (sequence data here)
  
  >Doug_NoIndex_L005_R1_001_contig_4.g13.t2
  
  (sequence data here)
Comments.UserDisplayName: ZincFingers
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-08T21:38:54.943Z
Comments.Text: 1. What exactly have you tried? 2. What does them not being functional mean? 3. Do you just want the single longest one or the top N longest?
Comments.UserDisplayName: Devon Ryan
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-08T21:41:54.887Z
Comments.Text: 1. I've tried all of the solutions here: https://www.biostars.org/p/107759/
  
  2. Non functional in that I couldn't get it work
  
  3. I would like the single longest read for each set of isoforms
Comments.UserDisplayName: ZincFingers
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-08T21:44:02.940Z
Comments.Text: @ZincFingers: Please put all of this in your original post and include whether isoforms are ever separated from each other by unrelated genes.
Comments.UserDisplayName: Devon Ryan
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-08T21:58:31.180Z
Comments.Text: You don't have the same sequence IDs there. Should we assume that the part after the last `.` in the ID line needs to be discarded?
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-08T22:02:11.683Z
Comments.Text: I'm not sure what you mean. The ".t1" or ".t2" indicates which isoform the read is from, I would like to discard the shorterst isoforms and retain the longest isoform.
Comments.UserDisplayName: ZincFingers
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-08T22:03:16.100Z
Comments.Text: OK, so they do need to be removed.
Comments.UserDisplayName: terdon
Tags: <fasta><isoform><filtering>
Title: How can longest isoforms (per gene) be extracted from a FASTA file?
ThreadText: # How can longest isoforms (per gene) be extracted from a FASTA file?
  
  (Asked by: username_0 on 2017-06-08)
  
  Is there a convenient way to extract the longest isoforms from a transcriptome fasta file? I had found some [scripts on biostars](https://www.biostars.org/p/107759/) but none are functional and I'm having difficulty getting them to work.
  
  I'm aware that the longest isoforms aren't necessarily 'the best' but it will suit my purposes.
  
  The fasta was generated via Augustus. Here is what the fasta file looks like currently (sequence shortened to save space)
  
  ```
  >Doug_NoIndex_L005_R1_001_contig_2.g7.t1
  atggggcataacatagagactggtgaacgtgctgaaattctacttcaaagtctacctgattcgtatgatcaactcatca
  ttaatataaccaaaaacctagaaattctagccttcgatgatgttgcagctgcggttcttgaagaagaaagtcggcgcaagaacaaagaagatagaccg
  >Doug_NoIndex_L005_R1_001_contig_2.g7.t2
  atggggcataacatagagactggtgaacgtgctgaaattctacttcaaagtctacctgattcgtatgatcaactcatca
  
  ```
  
  The format is as such:
  
  ```
  Gene 1 isoform 1  
  Gene 1 isoform 2  
  Gene 2 isoform 1  
  Gene 2 isoform 2   
  
  ```
  
  and so forth. There are several genes that have more than one pair of isoforms (up to 3 or 4). There are roughly 80,000 total transcripts, probably 25,000 genes. I would like to extract the single longest isoform for each gene.
  
  **Tags:** fasta, isoform, filtering
  
  > **Comments:**
  > **username_1**: It depends on the exact format of your fasta. Where did you get it? What do the fasta header lines look like? Can you paste a small example?  
  > **username_0**: >Doug_NoIndex_L005_R1_001_contig_4.g13.t1  (sequence data here)  >Doug_NoIndex_L005_R1_001_contig_4.g13.t2  (sequence data here)  
  > **username_2**: 1. What exactly have you tried? 2. What does them not being functional mean? 3. Do you just want the single longest one or the top N longest?  
  > **username_0**: 1. I've tried all of the solutions here: https://www.biostars.org/p/107759/  2. Non functional in that I couldn't get it work  3. I would like the single longest read for each set of isoforms  
  > **username_2**: @username_0: Please put all of this in your original post and include whether isoforms are ever separated from each other by unrelated genes.  
  > **username_3**: You don't have the same sequence IDs there. Should we assume that the part after the last `.` in the ID line needs to be discarded?  
  > **username_0**: I'm not sure what you mean. The ".t1" or ".t2" indicates which isoform the read is from, I would like to discard the shorterst isoforms and retain the longest isoform.  
  > **username_3**: OK, so they do need to be removed.  
  
  ---
  
  ## 5 answers
  
  ### Answer by username_2 on 2017-06-09. Score: 6 | Accepted answer
  
  While the solution from <https://bioinformatics.stackexchange.com/users/96/username_6-username_6> should work (after adjusting for possible python3 incompatibility), the following is a shorter and less memory demanding method that uses biopython:
  
  ```python
  #!/usr/bin/env python
  from Bio import SeqIO
  import sys
  
  lastGene = None
  longest = (None, None)
  for rec in SeqIO.parse(sys.argv[1], "fasta"):
      gene = ".".join(rec.id.split(".")[:-1])
      l = len(rec)
      if lastGene is not None:
          if gene == lastGene:
              if longest[0] < l:
                  longest = (l, rec)
          else:
              lastGene = gene
              SeqIO.write(longest[1], sys.stdout, "fasta")
              longest = (l, rec)
      else:
          lastGene = gene
          longest = (l, rec)
  SeqIO.write(longest[1], sys.stdout, "fasta")
  
  ```
  
  If you saved this as `filter.py`, then `filter.py original.fa > subset.fa` would be the command to use.
  
  
  ---
  
  ### Answer by username_4 on 2017-06-13. Score: 3
  
  Late to the party here, but I like to try to avoid writing scripts when some command line magic will do. It's good practice to index your `FASTA` so use it.
  
  ```
  samtools faidx <myfasta.fa>
  
  ```
  
  A little bit of `awk` and `sort` can determine the largest isoform of each gene (and they don't even have to be sorted by name in the source `FASTA`).
  
  ```
  awk -F'[\t.]' '{print $1,$2,$3,$4}' <myfasta.fa>.fai | sort -k4nr,4 | sort -uk1,2 | cut -f1-3 -d' '| tr ' ' '.' > selection.ls
  
  ```
  
  * The `-F[\t.]` specifies both the `fai` tab delimiters, and the full stops in your contig names (`Doug_NoIndex_L005_R1_001_contig_2.g7.t1`) as delimiters for `awk`.
  * For each line we print four columns. `$1` to `$3` are the three components of your contig name (after splitting on the `.`), and column $4 is the column in the index for sequence length.
  * We use sort to organise the lines by `$4` (the contig length). Reverse sort with `-r` so the longest genes appear first.
  * We then subsort by columns `$1` and `$2` (the contig and gene name (`g1`...)). `-u` then drops any repeat of column `$1` and `$2` pairs (that is, all lines after the first pair -- the longest length).
  * Finally, we reassemble `$1` to `$3` with `cut` and `tr` to give you a list of sequence names.
  
  A loop with our old friend `samtools faidx` can put your chosen sequences out to a new `FASTA`.
  
  ```
  while read contig;
      do samtools faidx <myfasta.fa> $contig >> selection.fa;
  done < selection.ls
  
  ```
  
  Don't forget to index it, they are useful!
  
  ```
  samtools faidx selection.fa
  
  ```
  
  > **Comments:**
  > **username_4**: Oh and if you wanted to remove the .tN part, just switch out that `cut -f1-3 -d' '` for `cut -f1-2 -d' '`  
  
  ---
  
  ### Answer by username_5 on 2017-06-09. Score: 3
  
  Here is a solution in R. Could get really slow with big files. Works for the example you posted.
  
  ```
  library(Biostrings)
  
  ## read your fasta in as Biostrings object
  fasta.s <- readDNAStringSet("sample.fa")
  
  ## get the read names (in your case it has the isoform info)
  names.fasta <- names(fasta.s)
  
  ## extract only the relevant gene and isoform id (split name by the period symbol)
  gene.iso <- sapply(names.fasta,function(j) cbind(unlist(strsplit(j,'\\.'))[2:3]))
  
  ## convert to good data.frame = transpose result from previous step and add relevant column names
  gene.iso.df <- data.frame(t(gene.iso))
  colnames(gene.iso.df) <- c('gene','isoform')
  
  ## and length of isoforms
  gene.iso.df$width <- width(fasta.s)
  
  ## split data.frame into list with entry for each gene
  gene.iso.df.split <- split(gene.iso.df,gene.iso.df$gene)
  
  ## optional to keep all the information but really you just need indices
  ##gene.iso.df.split.best <- lapply(gene.iso.df.split,function(x) x[order(x$width)[1],])
  
  ## pull out the longest isoform ID for each gene (in case of a tie just take the first one)
  best.id <- sapply(gene.iso.df.split,function(x) row.names(x)[order(x$width)[1]])
  
  ## subset your original reads with the subset
  fasta.s.best <- fasta.s[best.id]
  
  ## export new fastafile containing longest isoform per gene
  writeXStringSet(fasta.s, filepath='sample_best_isoform.fasta')
  
  ```
  
  
  ---
  
  ### Answer by username_6 on 2017-06-08. Score: 1
  
  This solution is will work for your example and probably all Augustus-derived Fastas, but mileage will vary beyond that.
  
  ```python
  #!/usr/bin/env python
  from __future__ import print_function
  import sys
  
  def parse_fasta(data):
      """Stolen shamelessly from http://stackoverflow.com/a/7655072/459780."""
      name, seq = None, []
      for line in data:
          line = line.rstrip()
          if line.startswith('>'):
              if name:
                  yield (name, ''.join(seq))
              name, seq = line, []
          else:
              seq.append(line)
      if name:
          yield (name, ''.join(seq))
  
  isoforms = dict()
  for defline, sequence in parse_fasta(sys.stdin):
      geneid = '.'.join(defline[1:].split('.')[:-1])
      if geneid in isoforms:
          otherdefline, othersequence = isoforms[geneid]
          if len(sequence) > len(othersequence):
              isoforms[geneid] = (defline, sequence)
      else:
          isoforms[geneid] = (defline, sequence)
  
  for defline, sequence in isoforms.values():
      print(defline, sequence, sep='\n')
  
  ```
  
  Save as `longest.py`, and invoke on the command line like so.
  
  ```
  python longest.py < intput.fasta > output.fasta
  
  ```
  
  ---
  
  The line `geneid = '.'.join(defline[1:].split('.')[:-1])` is key. Let me break it down.
  
  * `defline[1:]`: ignore the first character of the defline (the `>` symbol)
  * `defline[1:].split('.')`: split the string on `.` symbol
  * `defline[1:].split('.')[:-1]`: ignore the last value after the split (the isoform name)
  * `'.'.join(defline[1:].split('.')[:-1])`: join the split values again
  
  This value will be the same for all isoforms for the same gene, so we use it in the script to keep track of the longest sequence associated with each gene ID.
  
  > **Comments:**
  > **username_2**: This could save a lot of memory and a bit of time if you write out and reset `isoforms` whenever `geneid not in isoforms`. OP says isoforms are always next to each other.  
  > **username_0**: Hi username_6, I tried running this script and received the following error message:        File "longest.py", line 29          print(defline, sequence, sep='\n')                                  ^       SyntaxError: invalid syntax  
  > **username_2**: @username_0: Change it to `print(defline)` and `print(sequence)` on separate lines.  
  > **username_6**: Python 2 / 3 incompatibility. Will update.  
  
  ---
  
  ### Answer by username_3 on 2017-06-08. Score: 1
  
  An ex-coworker (Josep Avril) has written a couple of very useful little scripts that convert fasta to tbl (`seqID<TAB>Sequence`) and back again. These are extremely handy for this sort of thing (and are included at the end of this answer). Using them, you can convert your fasta to a one sequence per line format, keep the longest sequence with a simple awk script and convert back to fasta.
  
  I am assuming that the part after the last `.` in the ID line should be removed (so that `Doug_NoIndex_L005_R1_001_contig_2.g7.t1` and `Doug_NoIndex_L005_R1_001_contig_2.g7.t2` both map to `Doug_NoIndex_L005_R1_001_contig_2.g7`). If that is a correct assumption, this should work for you:
  
  ```
  $ FastaToTbl file.fa | sed 's/\.[^.]*\t/\t/' | 
      awk -F"\t" '{
                      if(length($2)>length(a[$1])){
                          a[$1]=$0
                      }
                  }
                  END{
                      for(i in a){
                          print a[i]
                      }
                  }' | TblToFasta
  >Doug_NoIndex_L005_R1_001_contig_2.g7 
  atggggcataacatagagactggtgaacgtgctgaaattctacttcaaagtctacctgat
  tcgtatgatcaactcatcattaatataaccaaaaacctagaaattctagccttcgatgat
  gttgcagctgcggttcttgaagaagaaagtcggcgcaagaacaaagaagatagaccg
  
  ```
  
  One possible issue with that approach is that it needs to keep one sequence per ID in memory. If you have a huge fasta file and not much memory, that might be a problem. If that's the case, you can first sort the file to ensure that similar IDs always appear together and then print each line as soon as you find the next ID:
  
  ```
  FastaToTbl file.fa | sed 's/\.[^.]*\t/\t/' | LC_ALL=C sort -t$'\t' -k1,1 |  
      awk -F"\t" '{
                      if(NR>1 && prev!=$1){
                          print a[$1]; 
                          prev=$1
                      } 
                      if(length($2)>length(a[$1])){
                          a[$1]=$0
                      }
                  }
                  END{
                      print a[$1]
                  }' | TblToFasta
  
  ```
  
  ---
  
  * FastaToTbl
  
    ```c
    #!/usr/bin/awk -f
    {
            if (substr($1,1,1)==">")
            if (NR>1)
                        printf "\n%s\t", substr($0,2,length($0)-1)
            else 
                printf "%s\t", substr($0,2,length($0)-1)
            else 
                    printf "%s", $0
    }END{printf "\n"}
  
    ```
  * TblToFasta
  
    ```c
    #! /usr/bin/awk -f
    {
      sequence=$NF
  
      ls = length(sequence)
      is = 1
      fld  = 1
      while (fld < NF)
      {
         if (fld == 1){printf ">"}
         printf "%s " , $fld
  
         if (fld == NF-1)
          {
            printf "\n"
          }
          fld = fld+1
      }
      while (is <= ls)
      {
        printf "%s\n", substr(sequence,is,60)
        is=is+60
      }
    }
  
    ```
  
  
  ---
  
  