text: user3876: 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?
  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  
  
  
  user1310: 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
  
  
original_text: 3876: <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>
   1310: <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>
  