Answers.Body: <p>Here a python solution:</p>
  
  <pre><code>#!/usr/bin/env python
  from __future__ import print_function
  from Bio import SeqIO
  from argparse import ArgumentParser
  
  
  def main():
      args = get_args()
      last_gene = None
      transcripts = []
  
      for record in SeqIO.parse(args.input, "fasta"):
          gene = record.id.split("_")[0]
  
          if last_gene is None:
              last_gene = gene
  
          if last_gene == gene:
              transcripts.append(record)
          else:
              print_longest_transcripts(transcripts, args.number)
              last_gene = gene
              transcripts = [record]
  
      print_longest_transcripts(transcripts, args.number)
  
  
  def print_longest_transcripts(transcripts, number):
      longest = sorted(transcripts, key=lambda x: len(x.seq), reverse=True)
  
      for record in longest[:number]:
          print("&gt;"+record.id, record.seq, sep="\n")
  
  
  def get_args():
      parser = ArgumentParser(description="")
      parser.add_argument("-i", "--input", help="input fasta", required=True)
      parser.add_argument("-n", "--number", help="number of longest transcript per gene", type=int, required=True)
      return parser.parse_args()
  
  
  if __name__ == '__main__':
      main()
  </code></pre>
  
  <p>Run it like this:</p>
  
  <pre><code>python longest.py -i input.fasta -n 3
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-12-21T14:01:51.660Z
Answers.Text: wow, thanks, finswimmer, for this script; although I liked the samtools solution, this python code is superior because it runs on one of my input files in 0m 1.857s while the samtools commands took almost 2h for writing the list "longestID.txt" for the same input file. A small suggestion: Because we want fasta formatted output, we could add to line no. 33 a fasta hook, like this: print(">"+record.id, record.seq, sep="\n")
Answers.UserDisplayName: Pedro
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-12-21T14:06:34.103Z
Answers.Text: You're right with your suggestion. I missed it and edited my answer. I'm surprised that the python solution is faster here. Would be interesting to find the bottleneck in the `samtools` version. You could try increasing the number of parallel jobs by applying the `-j` parameter for `parallels` e.g. `-j 8`.
Answers.UserDisplayName: finswimmer
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-12-21T12:25:01.593Z
Answers.LastActivityDate: 2018-12-21T14:04:11.763Z
Answers.LastEditDate: 2018-12-21T14:04:11.763Z
Answers.UserDisplayName: finswimmer
Answers.Body: <p>I know you are asking for a python solution. But maybe you are interested in other ways.</p>
  
  <p>You can use <code>samtools faidx</code> to index your fasta file. The resulting file contains the information about the length of every sequence. We can sort by length and extract the first n transcripts per gene.</p>
  
  <p><strong>1. Index fasta</strong></p>
  
  <pre><code>samtools faidx input.fa
  </code></pre>
  
  <p><strong>2. Extract IDs of the n longest transcripts per gene</strong></p>
  
  <pre><code>sed 's/_.\+//g' input.fa.fai | sort -u | 
      parallel -k 'grep {} -m 3 &lt;(sort -k2,2 -r --sort g input.fa.fai) |
          cut -f1' &gt; longestID.txt
  </code></pre>
  
  <p>The <code>sed</code> command replaces everything from <code>_</code> onwards in each line, leaving behind only the gene name. <code>sort -u</code> is used to get a list of distinct gene names. For each of theses gene names we <code>grep</code> for the first 3 matches in the descending sorted by length index file. As we just want the sequence ids, we select only the first column with <code>cut</code>. Adopt the <code>-m</code> parameter if you want more or fewer transcripts.</p>
  
  <p><strong>3. Get sequences</strong></p>
  
  <pre><code>samtools faidx input.fa -r longestID.txt &gt; output.fa
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-12-21T10:35:23.037Z
Answers.Text: Thanks for your suggestion, fin swimmer; this solution is elegant and completely fine with me (and I can understand it!). I just thought that modifying the python code would be easier. I didn't know that the samtools index can be used like this, and also learned from your answer about the parallel -k and the grep -m parameters. Thank you for your help and solving my problem.
Answers.UserDisplayName: Pedro
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-12-21T08:24:54.040Z
Answers.LastActivityDate: 2018-12-21T11:02:40.580Z
Answers.LastEditDate: 2018-12-21T11:02:40.580Z
Answers.UserDisplayName: finswimmer
Body: <p>The question of how one can extract the longest isoform per gene from a multi-fasta file has been answered expertly in this thread: 
  <a href="https://bioinformatics.stackexchange.com/questions/595/how-can-longest-isoforms-per-gene-be-extracted-from-a-fasta-file">How can longest isoforms (per gene) be extracted from a FASTA file?</a></p>
  
  <p>I wonder: how can the python script "filter.py" by Devon Ryan, that is marked as best answer, be modified so that it extracts the N longest isoform sequences per gene instead of only the single longest isoform? </p>
  
  <p>Any help on this to a python novice is very much appreciated. </p>
  
  <p>Here is an example of how the data looks (two genes with 5 isoformes each, and I need the 3 longest isoforms for each gene):</p>
  
  <pre><code>&gt;Animal1|gene000001_1  
  TTLQLLNELVLLLHYTTLLLHTKKTKESTINNYFKQKEQK  
  &gt;Animal1|gene000001_2  
  CYYYTTQHYYYIQKKQKNLQSTITSNKKNKNKHAK  
  &gt;Animal1|gene000001_3  
  ARDMLDYTTTTKRVSATTTLHNTTTTYKKNKRIYNQQLLQTKRTKISMQNKNLNQINYKK  
  LQETTEKHAK  
  &gt;Animal1|gene000001_4  
  HHAVNFQPMSEATGIYLNSTEIRNSLLQYSSWSWFTFHTCKLWNVTF  
  &gt;Animal1|gene000001_5  
  EILCYSIPPGVGLLFTPVNFGMLRSKLIKVRILDIINFNCSLLICLPPYLFQIDVGEDKD  
  DPVVEIVLYLACSI  
  &gt;Animal1|gene000002_1  
  TTLQLLNELVLLLHYTTLLLHTKKTKESTINNYFKQKEQK  
  &gt;Animal1|gene000002_2  
  CYYYTTQHYYYIQKKQKNLQSTITSNKKNKNKHAK  
  &gt;Animal1|gene000002_3  
  ARDMLDYTTTTKRVSATTTLHNTTTTYKKNKRIYNQQLLQTKRTKISMQNKNLNQINYKK  
  LQETTEKHAK  
  &gt;Animal1|gene000002_4  
  LCSAVARQKTKSRKYDATVDTCFSNIKININVLSFSLIVNPKYSKVSYQYGTNDAFRKFF  
  ITSRSLR  
  &gt;Animal1|gene000002_5  
  ETLEYFGLTINEKESTLILILMLLKQVSTVASYFLDLVFCLATALHS  
  </code></pre>
  
  <p>As output, I would like to have a fasta file like the following: </p>
  
  <pre><code>&gt;Animal1|gene000001_5  
  EILCYSIPPGVGLLFTPVNFGMLRSKLIKVRILDIINFNCSLLICLPPYLFQIDVGEDKD  
  DPVVEIVLYLACSI  
  &gt;Animal1|gene000001_3  
  ARDMLDYTTTTKRVSATTTLHNTTTTYKKNKRIYNQQLLQTKRTKISMQNKNLNQINYKK  
  LQETTEKHAK  
  &gt;Animal1|gene000001_4  
  HHAVNFQPMSEATGIYLNSTEIRNSLLQYSSWSWFTFHTCKLWNVTF  
  &gt;Animal1|gene000002_3  
  ARDMLDYTTTTKRVSATTTLHNTTTTYKKNKRIYNQQLLQTKRTKISMQNKNLNQINYKK  
  LQETTEKHAK  
  &gt;Animal1|gene000002_4  
  LCSAVARQKTKSRKYDATVDTCFSNIKININVLSFSLIVNPKYSKVSYQYGTNDAFRKFF  
  ITSRSLR  
  &gt;Animal1|gene000002_5  
  ETLEYFGLTINEKESTLILILMLLKQVSTVASYFLDLVFCLATALHS  
  </code></pre>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-12-20T15:21:40.250Z
Comments.Text: Please [edit] your question and include an example of your fasta file so we don't need to guess from the other question.
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-12-20T23:16:06.423Z
Comments.Text: I can do this but it would be entirely different code, you are grouping objects by name and reporting 3 of each name that are the longest. Would you like me to write you a code that handled this a different way? Also, can you use the pandas library?
Comments.UserDisplayName: d_kennetz
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-12-21T10:30:10.137Z
Comments.Text: Thanks for your answer, d_kennetz; I did not expect that entirely different code would be necessary, but simply thought that adding one or two statements to the existing python script would be the easiest way to do the job. I can understand the suggestion of fin swimmer (see below) and it does exactly what I need, so it is not necessary that you write extra code for me. Thanks for your help.
Comments.UserDisplayName: Pedro
Tags: <rna-seq><filtering><isoform>
Title: How can I extract the longest N isoforms per gene from a fasta file?
ThreadText: # How can I extract the longest N isoforms per gene from a fasta file?
  
  (Asked by: username_0 on 2018-12-20)
  
  The question of how one can extract the longest isoform per gene from a multi-fasta file has been answered expertly in this thread:
  [How can longest isoforms (per gene) be extracted from a FASTA file?](https://bioinformatics.stackexchange.com/questions/595/how-can-longest-isoforms-per-gene-be-extracted-from-a-fasta-file)
  
  I wonder: how can the python script "filter.py" by Devon Ryan, that is marked as best answer, be modified so that it extracts the N longest isoform sequences per gene instead of only the single longest isoform?
  
  Any help on this to a python novice is very much appreciated.
  
  Here is an example of how the data looks (two genes with 5 isoformes each, and I need the 3 longest isoforms for each gene):
  
  ```
  >Animal1|gene000001_1  
  TTLQLLNELVLLLHYTTLLLHTKKTKESTINNYFKQKEQK  
  >Animal1|gene000001_2  
  CYYYTTQHYYYIQKKQKNLQSTITSNKKNKNKHAK  
  >Animal1|gene000001_3  
  ARDMLDYTTTTKRVSATTTLHNTTTTYKKNKRIYNQQLLQTKRTKISMQNKNLNQINYKK  
  LQETTEKHAK  
  >Animal1|gene000001_4  
  HHAVNFQPMSEATGIYLNSTEIRNSLLQYSSWSWFTFHTCKLWNVTF  
  >Animal1|gene000001_5  
  EILCYSIPPGVGLLFTPVNFGMLRSKLIKVRILDIINFNCSLLICLPPYLFQIDVGEDKD  
  DPVVEIVLYLACSI  
  >Animal1|gene000002_1  
  TTLQLLNELVLLLHYTTLLLHTKKTKESTINNYFKQKEQK  
  >Animal1|gene000002_2  
  CYYYTTQHYYYIQKKQKNLQSTITSNKKNKNKHAK  
  >Animal1|gene000002_3  
  ARDMLDYTTTTKRVSATTTLHNTTTTYKKNKRIYNQQLLQTKRTKISMQNKNLNQINYKK  
  LQETTEKHAK  
  >Animal1|gene000002_4  
  LCSAVARQKTKSRKYDATVDTCFSNIKININVLSFSLIVNPKYSKVSYQYGTNDAFRKFF  
  ITSRSLR  
  >Animal1|gene000002_5  
  ETLEYFGLTINEKESTLILILMLLKQVSTVASYFLDLVFCLATALHS  
  
  ```
  
  As output, I would like to have a fasta file like the following:
  
  ```
  >Animal1|gene000001_5  
  EILCYSIPPGVGLLFTPVNFGMLRSKLIKVRILDIINFNCSLLICLPPYLFQIDVGEDKD  
  DPVVEIVLYLACSI  
  >Animal1|gene000001_3  
  ARDMLDYTTTTKRVSATTTLHNTTTTYKKNKRIYNQQLLQTKRTKISMQNKNLNQINYKK  
  LQETTEKHAK  
  >Animal1|gene000001_4  
  HHAVNFQPMSEATGIYLNSTEIRNSLLQYSSWSWFTFHTCKLWNVTF  
  >Animal1|gene000002_3  
  ARDMLDYTTTTKRVSATTTLHNTTTTYKKNKRIYNQQLLQTKRTKISMQNKNLNQINYKK  
  LQETTEKHAK  
  >Animal1|gene000002_4  
  LCSAVARQKTKSRKYDATVDTCFSNIKININVLSFSLIVNPKYSKVSYQYGTNDAFRKFF  
  ITSRSLR  
  >Animal1|gene000002_5  
  ETLEYFGLTINEKESTLILILMLLKQVSTVASYFLDLVFCLATALHS  
  
  ```
  
  **Tags:** rna-seq, filtering, isoform
  
  > **Comments:**
  > **username_1**: Please [edit] your question and include an example of your fasta file so we don't need to guess from the other question.  
  > **username_2**: I can do this but it would be entirely different code, you are grouping objects by name and reporting 3 of each name that are the longest. Would you like me to write you a code that handled this a different way? Also, can you use the pandas library?  
  > **username_0**: Thanks for your answer, username_2; I did not expect that entirely different code would be necessary, but simply thought that adding one or two statements to the existing python script would be the easiest way to do the job. I can understand the suggestion of fin swimmer (see below) and it does exactly what I need, so it is not necessary that you write extra code for me. Thanks for your help.  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_3 on 2018-12-21. Score: 0 | Accepted answer
  
  Here a python solution:
  
  ```
  #!/usr/bin/env python
  from __future__ import print_function
  from Bio import SeqIO
  from argparse import ArgumentParser
  
  def main():
      args = get_args()
      last_gene = None
      transcripts = []
  
      for record in SeqIO.parse(args.input, "fasta"):
          gene = record.id.split("_")[0]
  
          if last_gene is None:
              last_gene = gene
  
          if last_gene == gene:
              transcripts.append(record)
          else:
              print_longest_transcripts(transcripts, args.number)
              last_gene = gene
              transcripts = [record]
  
      print_longest_transcripts(transcripts, args.number)
  
  def print_longest_transcripts(transcripts, number):
      longest = sorted(transcripts, key=lambda x: len(x.seq), reverse=True)
  
      for record in longest[:number]:
          print(">"+record.id, record.seq, sep="\n")
  
  def get_args():
      parser = ArgumentParser(description="")
      parser.add_argument("-i", "--input", help="input fasta", required=True)
      parser.add_argument("-n", "--number", help="number of longest transcript per gene", type=int, required=True)
      return parser.parse_args()
  
  if __name__ == '__main__':
      main()
  
  ```
  
  Run it like this:
  
  ```
  python longest.py -i input.fasta -n 3
  
  ```
  
  > **Comments:**
  > **username_0**: wow, thanks, username_3, for this script; although I liked the samtools solution, this python code is superior because it runs on one of my input files in 0m 1.857s while the samtools commands took almost 2h for writing the list "longestID.txt" for the same input file. A small suggestion: Because we want fasta formatted output, we could add to line no. 33 a fasta hook, like this: print(">"+record.id, record.seq, sep="\n")  
  > **username_3**: You're right with your suggestion. I missed it and edited my answer. I'm surprised that the python solution is faster here. Would be interesting to find the bottleneck in the `samtools` version. You could try increasing the number of parallel jobs by applying the `-j` parameter for `parallels` e.g. `-j 8`.  
  
  ---
  
  ### Answer by username_3 on 2018-12-21. Score: 1
  
  I know you are asking for a python solution. But maybe you are interested in other ways.
  
  You can use `samtools faidx` to index your fasta file. The resulting file contains the information about the length of every sequence. We can sort by length and extract the first n transcripts per gene.
  
  **1. Index fasta**
  
  ```
  samtools faidx input.fa
  
  ```
  
  **2. Extract IDs of the n longest transcripts per gene**
  
  ```
  sed 's/_.\+//g' input.fa.fai | sort -u | 
      parallel -k 'grep {} -m 3 <(sort -k2,2 -r --sort g input.fa.fai) |
          cut -f1' > longestID.txt
  
  ```
  
  The `sed` command replaces everything from `_` onwards in each line, leaving behind only the gene name. `sort -u` is used to get a list of distinct gene names. For each of theses gene names we `grep` for the first 3 matches in the descending sorted by length index file. As we just want the sequence ids, we select only the first column with `cut`. Adopt the `-m` parameter if you want more or fewer transcripts.
  
  **3. Get sequences**
  
  ```
  samtools faidx input.fa -r longestID.txt > output.fa
  
  ```
  
  > **Comments:**
  > **username_0**: Thanks for your suggestion, fin swimmer; this solution is elegant and completely fine with me (and I can understand it!). I just thought that modifying the python code would be easier. I didn't know that the samtools index can be used like this, and also learned from your answer about the parallel -k and the grep -m parameters. Thank you for your help and solving my problem.  
  
  ---
  
  