Answers.Body: <p>Here's one way using <a href="https://biopython.org/" rel="nofollow noreferrer">Biopython</a> and the <a href="https://biopython.org/wiki/SeqIO" rel="nofollow noreferrer">SeqIO</a> interface to read and write <a href="https://biopython.org/wiki/SeqRecord" rel="nofollow noreferrer">SeqRecord</a> objects. Using a generator, we can produce (yield) trimmed sequence records that can then be written to the output handle. This avoids storing an entire list of trimmed sequences in memory and is usually the way to go, assuming you're dealing with large files containing thousands of records:</p>
  <pre><code>import argparse
  import re
  
  from Bio import SeqIO
  </code></pre>
  <pre><code>def trim_reads(fileobj):
  
      pattern = (
          r'type=(?P&lt;type&gt;\w+),'
          r'start=(?P&lt;start&gt;\d+),'
          r'end=(?P&lt;end&gt;\d+),'
          r'length=(?P&lt;length&gt;\d+),'
          r'identity=(?P&lt;identity&gt;\d+(?:\.\d+)?)%'
      )
  
      for seq_record in SeqIO.parse(fileobj, &quot;fasta&quot;):
          match = re.search(pattern, seq_record.description)
          result = match.groupdict()
  
          start = int(result['start']) - 1
          end = int(result['end'])
  
          head = seq_record.seq[:start]
          tail = seq_record.seq[end:]
  
          seq_record.seq = head + tail
  
          yield seq_record
  </code></pre>
  <pre><code>def get_argument_parser():
  
      parser = argparse.ArgumentParser()
  
      parser.add_argument(
          &quot;--fasta&quot;,
           type=str,
           required=True,
           help=&quot;input fasta file&quot;,
      )
      parser.add_argument(
          &quot;--output&quot;,
          type=str,
          required=True,
          help=&quot;output tail.trimmed.fasta file&quot;,
      )
  
      return parser
  </code></pre>
  <pre><code>def main():
  
      parser = get_argument_parser()
      args = parser.parse_args()
  
      with open(args.fasta) as fasta, open(args.output, &quot;w&quot;) as out:
          SeqIO.write(trim_reads(fasta), out, &quot;fasta-2line&quot;)
  </code></pre>
  <pre><code>if __name__ == &quot;__main__&quot;:
      main()
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-12-24T15:21:40.630Z
Answers.Text: Its worth point out that `re` is notoriously slow and over large files parallelisation is needed across the Biopython object.
Answers.UserDisplayName: M__
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-12-27T10:41:36.597Z
Answers.Text: Hello Steve. Thank you so much for you reply but I tried running your code and it gave the following error: AttributeError: 'NoneType' object has no attribute 'groupdict'
Answers.UserDisplayName: pubsurfted
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-12-27T11:18:24.713Z
Answers.Text: No worries at all, @pubsurfted. Yes - [re.search](https://docs.python.org/3/library/re.html#re.search) will return None if the search pattern could not be found. You'll need to add some code to handle these cases. You might need a different search pattern, or maybe you just want to ignore these records? If, in your question above, you can provide a few examples of where the search pattern fails, I might be able to suggest a solution.
Answers.UserDisplayName: Steve
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-12-28T12:52:02.053Z
Answers.Text: Hello @Steve, I have edited the question. If you could please take a look, thank you.
Answers.UserDisplayName: pubsurfted
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-12-28T14:25:58.123Z
Answers.Text: @pubsurfted I've updated the regex above to handle the fact that the identity can also be a plain integer, i.e. the decimal places are optional. Also, FASTA files are often line wrapped at 60 characters, but this entirely optional. These are sometimes called multiline FASTA files. There's probably a way to change the line wrapping using Biopython.
Answers.UserDisplayName: Steve
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-12-28T14:32:19.817Z
Answers.Text: @pubsurfted There is indeed, using the ["fasta-2line"](https://biopython.org/docs/1.76/api/Bio.SeqIO.FastaIO.html#Bio.SeqIO.FastaIO.as_fasta_2line) format.
Answers.UserDisplayName: Steve
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-12-29T14:01:46.130Z
Answers.Text: Thank you so so much Steve! Also, the "fasta-2line" format is raising the error: ValueError: Unknown format 'fasta-2line'
Answers.UserDisplayName: pubsurfted
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-12-29T14:16:15.077Z
Answers.Text: It ran. I had to update the biopython package from 1.67 to 1.80.
Answers.UserDisplayName: pubsurfted
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-12-24T14:12:46.383Z
Answers.LastActivityDate: 2022-12-28T14:33:02.990Z
Answers.LastEditDate: 2022-12-28T14:33:02.990Z
Answers.UserDisplayName: Steve
Body: <p>I have a FASTA file that has reads like this:</p>
  <blockquote>
  <p>&gt;SRR5655563.745 745 length=126 (type=T,start=2,end=12,length=11,identity=81.8182%) CttttatttgttGTGTGTGAGGTTTGATTTGATGGGAAAATATCTTGAATCTGCGGCGAGGTTGGAAGAGCTATCGCGGATTGTGTCATCTGCTGCGAAGCCCAATAGGTCAAAGGGAACGCTACC</p>
  </blockquote>
  <p>As you can see, there is information about the start &quot;<code>start=2</code>&quot; and end &quot;<code>end=12</code>&quot; of a sequence within the header. I would like to slice the sequence like [start:end] and keep the rest of it. e.g</p>
  <p>The part I would like to trim:</p>
  <blockquote>
  <p>Ctttggtttcctttt</p>
  </blockquote>
  <p>And after trimming, I would like to keep the rest of the read:</p>
  <blockquote>
  <p>GTGTGTGAGGTTTGATTTGATGGGAAAATATCTTGAATCTGCGGCGAGGTTGGAAGAGCTATCGCGGATTGTGTCATCTGCTGCGAAGCCCAATAGGTCAAAGGGAACGCTACC</p>
  </blockquote>
  <p>So the new FASTA file will have reads like this:</p>
  <blockquote>
  <p>&gt;SRR5655563.745 745 length=126 (type=T,start=2,end=12,length=11,identity=81.8182%)
  GTGTGTGAGGTTTGATTTGATGGGAAAATATCTTGAATCTGCGGCGAGGTTGGAAGAGCTATCGCGGATTGTGTCATCTGCTGCGAAGCCCAATAGGTCAAAGGGAACGCTACC</p>
  </blockquote>
  <p>I have written a python script that takes a FASTA file as an input and outputs a trimmed FASTA file, like this:</p>
  <pre><code>tailtrim.py --fasta input.fasta --output output.tail.trimmed.fasta
  </code></pre>
  <p>The code does trim the reads where I want it to.</p>
  <pre><code>import numpy as np
  import argparse
  import re
  
  def main():
      parser = argparse.ArgumentParser()
  
      parser.add_argument(&quot;--fasta&quot;, type=str, required=True, help=&quot;input fasta file&quot;)
      parser.add_argument(&quot;--output&quot;, type=str, required=True, help=&quot;output tail.trimmed.fasta file&quot;)
  
      args = parser.parse_args()
  
      out_fasta = open(args.tail.trimmed.fasta, 'w')
  
      for line in open(args.fasta):
          header = line.rstrip()
          read = next(args.fasta).rstrip()
  
      #Extract header into a list
          header_lst = []
          header_lst= header.split(&quot; &quot;)
      
          tmp_start_lst = []
          tmp_end_lst = []
          tmp_lst = []
  
          for i in header_lst:
              if i.startswith(&quot;(&quot;) and i.endswith(&quot;)&quot;):
                  tmp_lst.append(i.split(&quot;,&quot;))
          #print(tmp_lst)
  
      # Extract the start and end numbers in a list
          for i in tmp_lst:
              for str in i: 
                  if str.startswith(&quot;start=&quot;):
                      tmp_start_lst.append(str)
                  if str.startswith(&quot;end=&quot;):
                      tmp_end_lst.append(str)
  
          pattern = re.compile(r'\d+$')
  
          start_lst = []
          for str in tmp_start_lst:
              match = re.search(pattern, str)
              start_lst.append(int(match.group()))
  
          end_lst = []
          for str in tmp_end_lst:
              match = re.search(pattern, str)
              end_lst.append(int(match.group()))
  
          #Extract reads into a list to manipulate i.e. trim
          read_lst = []
          read_lst.append(read)
  
          new_lst = []
          
          for str in read_lst:
              for num_start in start_lst:
                  for num_end in end_lst:
                      s = str.replace(str[(num_start-1):num_end],&quot;&quot;)
                      new_lst.append(s)
          
  </code></pre>
  <p>I do get the trimmed reads, but what I'm having trouble figuring out is how to write these trimmed reads with their corresponding headers to a new FASTA file.</p>
  <p>Also if there is any way that I could improve my code or employ better strategy to approach the problem, I would love to know.</p>
  <p>Best wishes.</p>
  <p>Edit: I tried running the solution by the super helpful user @Steve.</p>
  <p>The command I ran:</p>
  <pre><code>python3 tailTrimming.py --fasta tailtrim.fasta --output testtrim.trim.fasta
  </code></pre>
  <p>The command line error:</p>
  <pre><code>Traceback (most recent call last):
    File &quot;/home/x/Desktop/fastq/tailTrimming.py&quot;, line 61, in &lt;module&gt;
      main()
    File &quot;/home/x/Desktop/fastq/tailTrimming.py&quot;, line 58, in main
      SeqIO.write(trim_reads(fasta), out, &quot;fasta&quot;)
    File &quot;/home/x/.local/lib/python3.10/site-packages/Bio/SeqIO/__init__.py&quot;, line 478, in write
      count = writer_class(fp).write_file(sequences)
    File &quot;/home/x/.local/lib/python3.10/site-packages/Bio/SeqIO/Interfaces.py&quot;, line 209, in write_file
      count = self.write_records(records)
    File &quot;/home/x/.local/lib/python3.10/site-packages/Bio/SeqIO/Interfaces.py&quot;, line 193, in write_records
      for record in records:
    File &quot;/home/x/Desktop/fastq/tailTrimming.py&quot;, line 20, in trim_reads
      result = match.groupdict()
  AttributeError: 'NoneType' object has no attribute 'groupdict'
  </code></pre>
  <p>This is the whole input file I supplied:</p>
  <pre><code>&gt;SRR5655563.9676 9676 length=126 (type=T,start=1,end=11,length=11,identity=81.8182%)
  tttctttctttAGTTCCAATAGTCTTGACCACCTTTTATCAACAGAAGCTTGGCCTTGTCCAGTTATATATTGTATTGTACCCTCAGATTTCTCGCACAAT
  &gt;SRR5655563.10087 10087 length=126 (type=T,start=2,end=16,length=15,identity=80%)
  GtttgtttcattttttCCAAGGAAGAAATTCGACCTTCCAGCGATTTGTGGAGGTCATTCATCTGTTGAGTAAAGGATTGGTGCGCGACCTTGTAATCCTTGTGGAGGGTGTAG
  &gt;SRR5655563.13902 13902 length=126 (type=T,start=2,end=16,length=15,identity=80%)
  CtttttcttcttgtttGATGTAAAAAATCTTCGGTATCTAAAACGTATGAAACCTTAGATGAACTGCAAGCTATCTTACACTCTAGACCATTGATAATTT
  &gt;SRR5655563.18489 18489 length=126 (type=T,start=2,end=12,length=11,identity=81.8182%)
  CtttttcttcttCTTGGAAGTTTTGGATTCTGCAGGGACATTTTCTTTCTTTTCTCCATCTCCATCAGTCTCAACACCTGTCTTATCCTGAGCACCTTGAGACTCATCTTGACCATCTTTACTTT
  &gt;SRR5655563.20782 20782 length=126 (type=T,start=2,end=15,length=14,identity=71.4286%)
  CttcattttagttttCAACTCCTCTACCAACTGGGATGGTGCATTCAACTGTAGAACTGTTTCACGAATCTCCTTGAGAACACTGAAGTCTCCCTCAGTTAACTTCTTTTTCAAATTATTGA
  &gt;SRR5655563.23972 23972 length=126 (type=T,start=2,end=29,length=28,identity=75%)
  CttttttgttgttctcttctttcttctttGCATGGCTTGATGATGACTGATCCGATTGATTATTTCTTGATAGAATATCAATGTCATCACCAATTGTTGAGT
  &gt;SRR5655563.31739 31739 length=126 (type=T,start=1,end=10,length=10,identity=80%)
  ttttcttcttCTTGCATGAAAGTGGGTTTGGACCCTTAGCTCTTTTTCTCTTGAACTGAGGTTT
  &gt;SRR5655563.34244 34244 length=126 (type=T,start=2,end=14,length=13,identity=76.9231%)
  CtttcctttattttCTTGGAGGTTAGCCAGGATATTGGGATCAAAGCATGTAGCAAAACACACCTTGCTATTTCTTGCTCGTACAGATACATGGGATAGGA
  &gt;SRR5655563.38334 38334 length=126 (type=T,start=2,end=12,length=11,identity=81.8182%)
  CtttattgttttCTCATCACCCATATTATATTTCTTAAGGGCAGATTCCAACTCTTCCATTGTGATGTACCCGCTCTTATCATTGTCAAAATATTCAAAGGCTTTATATAGGCGATCCTCTCTTTC
  </code></pre>
  <p>And this is the whole output I have gotten:</p>
  <pre><code>&gt;SRR5655563.9676 9676 length=126 (type=T,start=1,end=11,length=11,identity=81.8182%)
  AGTTCCAATAGTCTTGACCACCTTTTATCAACAGAAGCTTGGCCTTGTCCAGTTATATAT
  TGTATTGTACCCTCAGATTTCTCGCACAAT
  </code></pre>
  <p>Just one read, that's all. Whereas I would like to keep all the sequences because all of them have this <code>start</code>, <code>end</code> attribute. The sequence being on two lines in the output is not something I did. It is like this in the output file. So it's no longer FASTA file, right?</p>
  <p>Thanks in advance for any help!</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-12-24T13:17:09.760Z
Comments.Text: I think there's a typo in your example. The part to trim does not exist in the input sequence. Also, if you want to keep the remaining sequence (i.e. "the rest of it"), shouldn't the output include the sequence before the starting position? Or would you like to extract the sequence after the end position?
Comments.UserDisplayName: Steve
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-12-24T16:45:08.610Z
Comments.Text: Just to add in terms of performance `Python` is slow (loops and stuff). `re` is slow even for Python. Although `Perl` is not popular it will outperform the solutions here on a single CPU. If there are performance issues it's worth posting back - ultimately it depends on the size of the file and sensitivity to runtime. Huge files can simply fall over (i.e. never complete). However, I agree with Steve's solution: what Python is very cool at is splicing a string and that is super important to the solution. This is the bit that is really missing in the above code IMO
Comments.UserDisplayName: M__
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-12-24T20:36:36.823Z
Comments.Text: You made mistake describing the problem (sequence does not match 2-12 position intervals). If you update it, maybe you get more help. Also, it's not clear what kind of header you want to print in the record, exact same FASTA header?
Comments.UserDisplayName: Supertech
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-02-10T21:19:35.993Z
Comments.Text: i have the same problem but looks like i found some interesting solution @Steve. However, i have want to ask... what does this mean? pattern = ( r'type=(?P<type>\w+),' r'start=(?P<start>\d+),' r'end=(?P<end>\d+),' r'length=(?P<length>\d+),' r'identity=(?P<identity>\d+(?:\.\d+)?)%' )
  r'type= ??? what should i put here? r'start = ??? refers to the forward primer?
  r'end = reverse primer?
  r'length = the length of the final sequence after the trimming? what happens if there are sequences with different length afterwards? r'identity = the similarity between the primer and its sequence?
Comments.UserDisplayName: Damianos
Tags: <python><fasta>
Title: How to manually trim FASTA file sequences with the information provided in the header and store it into a new FASTA file? - Python
ThreadText: # How to manually trim FASTA file sequences with the information provided in the header and store it into a new FASTA file? - Python
  
  (Asked by: username_0 on 2022-12-24)
  
  I have a FASTA file that has reads like this:
  
  > >SRR5655563.745 745 length=126 (type=T,start=2,end=12,length=11,identity=81.8182%) CttttatttgttGTGTGTGAGGTTTGATTTGATGGGAAAATATCTTGAATCTGCGGCGAGGTTGGAAGAGCTATCGCGGATTGTGTCATCTGCTGCGAAGCCCAATAGGTCAAAGGGAACGCTACC
  
  As you can see, there is information about the start "`start=2`" and end "`end=12`" of a sequence within the header. I would like to slice the sequence like [start:end] and keep the rest of it. e.g
  
  The part I would like to trim:
  
  > Ctttggtttcctttt
  
  And after trimming, I would like to keep the rest of the read:
  
  > GTGTGTGAGGTTTGATTTGATGGGAAAATATCTTGAATCTGCGGCGAGGTTGGAAGAGCTATCGCGGATTGTGTCATCTGCTGCGAAGCCCAATAGGTCAAAGGGAACGCTACC
  
  So the new FASTA file will have reads like this:
  
  > >SRR5655563.745 745 length=126 (type=T,start=2,end=12,length=11,identity=81.8182%)
  > GTGTGTGAGGTTTGATTTGATGGGAAAATATCTTGAATCTGCGGCGAGGTTGGAAGAGCTATCGCGGATTGTGTCATCTGCTGCGAAGCCCAATAGGTCAAAGGGAACGCTACC
  
  I have written a python script that takes a FASTA file as an input and outputs a trimmed FASTA file, like this:
  
  ```
  tailtrim.py --fasta input.fasta --output output.tail.trimmed.fasta
  
  ```
  
  The code does trim the reads where I want it to.
  
  ```
  import numpy as np
  import argparse
  import re
  
  def main():
      parser = argparse.ArgumentParser()
  
      parser.add_argument("--fasta", type=str, required=True, help="input fasta file")
      parser.add_argument("--output", type=str, required=True, help="output tail.trimmed.fasta file")
  
      args = parser.parse_args()
  
      out_fasta = open(args.tail.trimmed.fasta, 'w')
  
      for line in open(args.fasta):
          header = line.rstrip()
          read = next(args.fasta).rstrip()
  
      #Extract header into a list
          header_lst = []
          header_lst= header.split(" ")
      
          tmp_start_lst = []
          tmp_end_lst = []
          tmp_lst = []
  
          for i in header_lst:
              if i.startswith("(") and i.endswith(")"):
                  tmp_lst.append(i.split(","))
          #print(tmp_lst)
  
      # Extract the start and end numbers in a list
          for i in tmp_lst:
              for str in i: 
                  if str.startswith("start="):
                      tmp_start_lst.append(str)
                  if str.startswith("end="):
                      tmp_end_lst.append(str)
  
          pattern = re.compile(r'\d+$')
  
          start_lst = []
          for str in tmp_start_lst:
              match = re.search(pattern, str)
              start_lst.append(int(match.group()))
  
          end_lst = []
          for str in tmp_end_lst:
              match = re.search(pattern, str)
              end_lst.append(int(match.group()))
  
          #Extract reads into a list to manipulate i.e. trim
          read_lst = []
          read_lst.append(read)
  
          new_lst = []
          
          for str in read_lst:
              for num_start in start_lst:
                  for num_end in end_lst:
                      s = str.replace(str[(num_start-1):num_end],"")
                      new_lst.append(s)
          
  
  ```
  
  I do get the trimmed reads, but what I'm having trouble figuring out is how to write these trimmed reads with their corresponding headers to a new FASTA file.
  
  Also if there is any way that I could improve my code or employ better strategy to approach the problem, I would love to know.
  
  Best wishes.
  
  Edit: I tried running the solution by the super helpful user @username_1.
  
  The command I ran:
  
  ```
  python3 tailTrimming.py --fasta tailtrim.fasta --output testtrim.trim.fasta
  
  ```
  
  The command line error:
  
  ```
  Traceback (most recent call last):
    File "/home/x/Desktop/fastq/tailTrimming.py", line 61, in <module>
      main()
    File "/home/x/Desktop/fastq/tailTrimming.py", line 58, in main
      SeqIO.write(trim_reads(fasta), out, "fasta")
    File "/home/x/.local/lib/python3.10/site-packages/Bio/SeqIO/__init__.py", line 478, in write
      count = writer_class(fp).write_file(sequences)
    File "/home/x/.local/lib/python3.10/site-packages/Bio/SeqIO/Interfaces.py", line 209, in write_file
      count = self.write_records(records)
    File "/home/x/.local/lib/python3.10/site-packages/Bio/SeqIO/Interfaces.py", line 193, in write_records
      for record in records:
    File "/home/x/Desktop/fastq/tailTrimming.py", line 20, in trim_reads
      result = match.groupdict()
  AttributeError: 'NoneType' object has no attribute 'groupdict'
  
  ```
  
  This is the whole input file I supplied:
  
  ```
  >SRR5655563.9676 9676 length=126 (type=T,start=1,end=11,length=11,identity=81.8182%)
  tttctttctttAGTTCCAATAGTCTTGACCACCTTTTATCAACAGAAGCTTGGCCTTGTCCAGTTATATATTGTATTGTACCCTCAGATTTCTCGCACAAT
  >SRR5655563.10087 10087 length=126 (type=T,start=2,end=16,length=15,identity=80%)
  GtttgtttcattttttCCAAGGAAGAAATTCGACCTTCCAGCGATTTGTGGAGGTCATTCATCTGTTGAGTAAAGGATTGGTGCGCGACCTTGTAATCCTTGTGGAGGGTGTAG
  >SRR5655563.13902 13902 length=126 (type=T,start=2,end=16,length=15,identity=80%)
  CtttttcttcttgtttGATGTAAAAAATCTTCGGTATCTAAAACGTATGAAACCTTAGATGAACTGCAAGCTATCTTACACTCTAGACCATTGATAATTT
  >SRR5655563.18489 18489 length=126 (type=T,start=2,end=12,length=11,identity=81.8182%)
  CtttttcttcttCTTGGAAGTTTTGGATTCTGCAGGGACATTTTCTTTCTTTTCTCCATCTCCATCAGTCTCAACACCTGTCTTATCCTGAGCACCTTGAGACTCATCTTGACCATCTTTACTTT
  >SRR5655563.20782 20782 length=126 (type=T,start=2,end=15,length=14,identity=71.4286%)
  CttcattttagttttCAACTCCTCTACCAACTGGGATGGTGCATTCAACTGTAGAACTGTTTCACGAATCTCCTTGAGAACACTGAAGTCTCCCTCAGTTAACTTCTTTTTCAAATTATTGA
  >SRR5655563.23972 23972 length=126 (type=T,start=2,end=29,length=28,identity=75%)
  CttttttgttgttctcttctttcttctttGCATGGCTTGATGATGACTGATCCGATTGATTATTTCTTGATAGAATATCAATGTCATCACCAATTGTTGAGT
  >SRR5655563.31739 31739 length=126 (type=T,start=1,end=10,length=10,identity=80%)
  ttttcttcttCTTGCATGAAAGTGGGTTTGGACCCTTAGCTCTTTTTCTCTTGAACTGAGGTTT
  >SRR5655563.34244 34244 length=126 (type=T,start=2,end=14,length=13,identity=76.9231%)
  CtttcctttattttCTTGGAGGTTAGCCAGGATATTGGGATCAAAGCATGTAGCAAAACACACCTTGCTATTTCTTGCTCGTACAGATACATGGGATAGGA
  >SRR5655563.38334 38334 length=126 (type=T,start=2,end=12,length=11,identity=81.8182%)
  CtttattgttttCTCATCACCCATATTATATTTCTTAAGGGCAGATTCCAACTCTTCCATTGTGATGTACCCGCTCTTATCATTGTCAAAATATTCAAAGGCTTTATATAGGCGATCCTCTCTTTC
  
  ```
  
  And this is the whole output I have gotten:
  
  ```
  >SRR5655563.9676 9676 length=126 (type=T,start=1,end=11,length=11,identity=81.8182%)
  AGTTCCAATAGTCTTGACCACCTTTTATCAACAGAAGCTTGGCCTTGTCCAGTTATATAT
  TGTATTGTACCCTCAGATTTCTCGCACAAT
  
  ```
  
  Just one read, that's all. Whereas I would like to keep all the sequences because all of them have this `start`, `end` attribute. The sequence being on two lines in the output is not something I did. It is like this in the output file. So it's no longer FASTA file, right?
  
  Thanks in advance for any help!
  
  **Tags:** python, fasta
  
  > **Comments:**
  > **username_1**: I think there's a typo in your example. The part to trim does not exist in the input sequence. Also, if you want to keep the remaining sequence (i.e. "the rest of it"), shouldn't the output include the sequence before the starting position? Or would you like to extract the sequence after the end position?  
  > **username_2**: Just to add in terms of performance `Python` is slow (loops and stuff). `re` is slow even for Python. Although `Perl` is not popular it will outperform the solutions here on a single CPU. If there are performance issues it's worth posting back - ultimately it depends on the size of the file and sensitivity to runtime. Huge files can simply fall over (i.e. never complete). However, I agree with username_1's solution: what Python is very cool at is splicing a string and that is super important to the solution. This is the bit that is really missing in the above code IMO  
  > **username_3**: You made mistake describing the problem (sequence does not match 2-12 position intervals). If you update it, maybe you get more help. Also, it's not clear what kind of header you want to print in the record, exact same FASTA header?  
  > **username_4**: i have the same problem but looks like i found some interesting solution @username_1. However, i have want to ask... what does this mean? pattern = ( r'type=(?P<type>\w+),' r'start=(?P<start>\d+),' r'end=(?P<end>\d+),' r'length=(?P<length>\d+),' r'identity=(?P<identity>\d+(?:\.\d+)?)%' ) r'type= ??? what should i put here? r'start = ??? refers to the forward primer? r'end = reverse primer? r'length = the length of the final sequence after the trimming? what happens if there are sequences with different length afterwards? r'identity = the similarity between the primer and its sequence?  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2022-12-24. Score: 3 | Accepted answer
  
  Here's one way using [Biopython](https://biopython.org/) and the [SeqIO](https://biopython.org/wiki/SeqIO) interface to read and write [SeqRecord](https://biopython.org/wiki/SeqRecord) objects. Using a generator, we can produce (yield) trimmed sequence records that can then be written to the output handle. This avoids storing an entire list of trimmed sequences in memory and is usually the way to go, assuming you're dealing with large files containing thousands of records:
  
  ```
  import argparse
  import re
  
  from Bio import SeqIO
  
  ```
  
  ```
  def trim_reads(fileobj):
  
      pattern = (
          r'type=(?P<type>\w+),'
          r'start=(?P<start>\d+),'
          r'end=(?P<end>\d+),'
          r'length=(?P<length>\d+),'
          r'identity=(?P<identity>\d+(?:\.\d+)?)%'
      )
  
      for seq_record in SeqIO.parse(fileobj, "fasta"):
          match = re.search(pattern, seq_record.description)
          result = match.groupdict()
  
          start = int(result['start']) - 1
          end = int(result['end'])
  
          head = seq_record.seq[:start]
          tail = seq_record.seq[end:]
  
          seq_record.seq = head + tail
  
          yield seq_record
  
  ```
  
  ```
  def get_argument_parser():
  
      parser = argparse.ArgumentParser()
  
      parser.add_argument(
          "--fasta",
           type=str,
           required=True,
           help="input fasta file",
      )
      parser.add_argument(
          "--output",
          type=str,
          required=True,
          help="output tail.trimmed.fasta file",
      )
  
      return parser
  
  ```
  
  ```
  def main():
  
      parser = get_argument_parser()
      args = parser.parse_args()
  
      with open(args.fasta) as fasta, open(args.output, "w") as out:
          SeqIO.write(trim_reads(fasta), out, "fasta-2line")
  
  ```
  
  ```
  if __name__ == "__main__":
      main()
  
  ```
  
  > **Comments:**
  > **username_2**: Its worth point out that `re` is notoriously slow and over large files parallelisation is needed across the Biopython object.  
  > **username_0**: Hello username_1. Thank you so much for you reply but I tried running your code and it gave the following error: AttributeError: 'NoneType' object has no attribute 'groupdict'  
  > **username_1**: No worries at all, @username_0. Yes - [re.search](https://docs.python.org/3/library/re.html#re.search) will return None if the search pattern could not be found. You'll need to add some code to handle these cases. You might need a different search pattern, or maybe you just want to ignore these records? If, in your question above, you can provide a few examples of where the search pattern fails, I might be able to suggest a solution.  
  > **username_0**: Hello @username_1, I have edited the question. If you could please take a look, thank you.  
  > **username_1**: @username_0 I've updated the regex above to handle the fact that the identity can also be a plain integer, i.e. the decimal places are optional. Also, FASTA files are often line wrapped at 60 characters, but this entirely optional. These are sometimes called multiline FASTA files. There's probably a way to change the line wrapping using Biopython.  
  > **username_1**: @username_0 There is indeed, using the ["fasta-2line"](https://biopython.org/docs/1.76/api/Bio.SeqIO.FastaIO.html#Bio.SeqIO.FastaIO.as_fasta_2line) format.  
  > **username_0**: Thank you so so much username_1! Also, the "fasta-2line" format is raising the error: ValueError: Unknown format 'fasta-2line'  
  > **username_0**: It ran. I had to update the biopython package from 1.67 to 1.80.  
  
  ---
  
  