text: user779: Is there a convenient way to extract the longest isoforms from a transcriptome fasta file? I had found some scripts on biostars 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.
  
  user298: 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
  #!/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
  #! /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
    }
  }
  
  
  
  user96: This solution is will work for your example and probably all Augustus-derived Fastas, but mileage will vary beyond that.
  #!/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.
  
  user383: 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')
  
  
  user215: 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
  
  
  user77: While the solution from https://bioinformatics.stackexchange.com/users/96/daniel-standage should work (after adjusting for possible python3 incompatibility), the following is a shorter and less memory demanding method that uses biopython:
  #!/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.
  
original_text: 779: <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>
   298: <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>
   96: <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>
   383: <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>
   215: <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>
   77: <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>
  