Answers.Body: <p>Here's one way using <a href="/questions/tagged/biopython" class="post-tag" title="show questions tagged &#39;biopython&#39;" aria-label="show questions tagged &#39;biopython&#39;" rel="tag" aria-labelledby="biopython-container">biopython</a>. Note that the <a href="https://biopython.org/wiki/Seq" rel="nofollow noreferrer"><code>Seq</code></a> object has a number of methods that act just like those of a Python string. One of these is the <a href="https://biopython.org/docs/1.80/api/Bio.Seq.html#Bio.Seq.SequenceDataAbstractBaseClass.find" rel="nofollow noreferrer"><code>find()</code></a> method, which returns the index of the first occurrence of the sequence if found. Replace with <a href="https://biopython.org/docs/1.80/api/Bio.Seq.html#Bio.Seq.SequenceDataAbstractBaseClass.rfind" rel="nofollow noreferrer"><code>rfind()</code></a> if you want the last occurrence if that makes more sense. This solution uses a generator to avoid storing the entire list of trimmed sequences in memory:</p>
  <pre><code>import argparse
  import contextlib
  import gzip
  import sys
  
  from Bio import SeqIO
  </code></pre>
  <pre><code>def trim_fasta(fileobj, sequence, num_bases):
      &quot;&quot;&quot;
      Trim num_bases upstream of the sequence provided
      &quot;&quot;&quot;
      for record in SeqIO.parse(fileobj, 'fasta'):
          index = record.seq.find(sequence)
          if index != -1 and index &gt; num_bases:
              record = record[index-num_bases:]
          yield record
  </code></pre>
  <pre><code>def open_file(filename, mode='r'):
      if filename.endswith('.gz'):
          return gzip.open(filename, mode)
      else:
          return open(filename, mode)
  </code></pre>
  <pre><code>def get_argument_parser():
  
      parser = argparse.ArgumentParser(add_help=False)
  
      group = parser.add_argument_group('input options')
      group.add_argument(
          &quot;fasta&quot;,
          type=str,
          metavar=&quot;FILE&quot;,
          help=&quot;The input FASTA file&quot;,
      )
      group.add_argument(
          &quot;sequence&quot;,
          type=str,
          metavar=&quot;STR&quot;,
          help=&quot;The sequence to find and trim&quot;,
      )
  
      group = parser.add_argument_group('trimming options')
      group.add_argument(
          &quot;-n&quot;,
          &quot;--num_bases&quot;,
          type=int,
          metavar=&quot;INT&quot;,
          default=1000,
          help=&quot;Trim INT bases upstream of input sequence&quot;,
      )
  
      group = parser.add_argument_group('output options')
      group.add_argument(
          &quot;-o&quot;,
          &quot;--output&quot;,
          type=str,
          metavar=&quot;FILE&quot;,
          default='-',
          help=&quot;Write the filtered output to FILE (default: stdout)&quot;,
      )
      group.add_argument(
          &quot;-f&quot;,
          &quot;--force&quot;,
          action='store_true',
          help=&quot;Overwrite the output file if it exists&quot;,
      )
  
      group = parser.add_argument_group('additional options')
      group.add_argument(
          &quot;-h&quot;,
          &quot;--help&quot;,
          action=&quot;help&quot;,
          help=&quot;Show this help message and exit&quot;,
      )
  
      return parser
  </code></pre>
  <pre><code>def main():
  
      parser = get_argument_parser()
      args = parser.parse_args()
  
      with contextlib.ExitStack() as stack:
          if args.output != '-':
              out = stack.enter_context(open_file(args.output, 'wt' if args.force else 'xt'))
          else:
              out = sys.stdout
  
          input_fasta = stack.enter_context(open_file(args.fasta))
          trimmed_reads = trim_fasta(input_fasta, args.sequence, args.num_bases)
  
          SeqIO.write(trimmed_reads, out, &quot;fasta-2line&quot;)
  
  
  if __name__ == '__main__':
      main()
  </code></pre>
  <hr />
  <hr />
  <pre><code>usage: trim_reads.py [-n INT] [-o FILE] [-f] [-h] FILE STR
  
  input options:
    FILE                  The input FASTA file
    STR                   The sequence to find and trim
  
  trimming options:
    -n INT, --num_bases INT
                          Trim INT bases upstream of input sequence
  
  output options:
    -o FILE, --output FILE
                          Write the filtered output to FILE (default: stdout)
    -f, --force           Overwrite the output file if it exists
  
  additional options:
    -h, --help            Show this help message and exit
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-26T12:31:19.623Z
Answers.LastActivityDate: 2023-01-26T12:44:12.530Z
Answers.LastEditDate: 2023-01-26T12:44:12.530Z
Answers.UserDisplayName: Steve
Body: <p>I need a quick way to trim multiple reads in a FASTA file. I need to trim everything that is 1kbp upstream of this sequence <code>AAGAGATGTTCAATCGTTTAAACAAATTCCAAGCTGCTTTAGCTTTGGCCCTTTACTCTCA</code>.</p>
  <p>I figure a quick python script might be the way to go but I'm not sure if there's a tool that already does that (i.e. Trimmomatic or something like it) that would be able to do that for me.</p>
  <p>Thank you!</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-01-27T19:42:46.470Z
Comments.Text: What is the source of these reads (i.e. sequencing platform)? Is it necessary to account for some error in the recognition sequence?
Comments.UserDisplayName: gringer
Tags: <python><trimming><multi-fasta>
Title: Trim reads 1kb upstream of sequence
ThreadText: # Trim reads 1kb upstream of sequence
  
  (Asked by: username_0 on 2023-01-25)
  
  I need a quick way to trim multiple reads in a FASTA file. I need to trim everything that is 1kbp upstream of this sequence `AAGAGATGTTCAATCGTTTAAACAAATTCCAAGCTGCTTTAGCTTTGGCCCTTTACTCTCA`.
  
  I figure a quick python script might be the way to go but I'm not sure if there's a tool that already does that (i.e. Trimmomatic or something like it) that would be able to do that for me.
  
  Thank you!
  
  **Tags:** python, trimming, multi-fasta
  
  > **Comments:**
  > **username_1**: What is the source of these reads (i.e. sequencing platform)? Is it necessary to account for some error in the recognition sequence?  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_2 on 2023-01-26. Score: 2 | Accepted answer
  
  Here's one way using [biopython](/questions/tagged/biopython "show questions tagged 'biopython'"). Note that the [`Seq`](https://biopython.org/wiki/Seq) object has a number of methods that act just like those of a Python string. One of these is the [`find()`](https://biopython.org/docs/1.80/api/Bio.Seq.html#Bio.Seq.SequenceDataAbstractBaseClass.find) method, which returns the index of the first occurrence of the sequence if found. Replace with [`rfind()`](https://biopython.org/docs/1.80/api/Bio.Seq.html#Bio.Seq.SequenceDataAbstractBaseClass.rfind) if you want the last occurrence if that makes more sense. This solution uses a generator to avoid storing the entire list of trimmed sequences in memory:
  
  ```
  import argparse
  import contextlib
  import gzip
  import sys
  
  from Bio import SeqIO
  
  ```
  
  ```
  def trim_fasta(fileobj, sequence, num_bases):
      """
      Trim num_bases upstream of the sequence provided
      """
      for record in SeqIO.parse(fileobj, 'fasta'):
          index = record.seq.find(sequence)
          if index != -1 and index > num_bases:
              record = record[index-num_bases:]
          yield record
  
  ```
  
  ```
  def open_file(filename, mode='r'):
      if filename.endswith('.gz'):
          return gzip.open(filename, mode)
      else:
          return open(filename, mode)
  
  ```
  
  ```
  def get_argument_parser():
  
      parser = argparse.ArgumentParser(add_help=False)
  
      group = parser.add_argument_group('input options')
      group.add_argument(
          "fasta",
          type=str,
          metavar="FILE",
          help="The input FASTA file",
      )
      group.add_argument(
          "sequence",
          type=str,
          metavar="STR",
          help="The sequence to find and trim",
      )
  
      group = parser.add_argument_group('trimming options')
      group.add_argument(
          "-n",
          "--num_bases",
          type=int,
          metavar="INT",
          default=1000,
          help="Trim INT bases upstream of input sequence",
      )
  
      group = parser.add_argument_group('output options')
      group.add_argument(
          "-o",
          "--output",
          type=str,
          metavar="FILE",
          default='-',
          help="Write the filtered output to FILE (default: stdout)",
      )
      group.add_argument(
          "-f",
          "--force",
          action='store_true',
          help="Overwrite the output file if it exists",
      )
  
      group = parser.add_argument_group('additional options')
      group.add_argument(
          "-h",
          "--help",
          action="help",
          help="Show this help message and exit",
      )
  
      return parser
  
  ```
  
  ```
  def main():
  
      parser = get_argument_parser()
      args = parser.parse_args()
  
      with contextlib.ExitStack() as stack:
          if args.output != '-':
              out = stack.enter_context(open_file(args.output, 'wt' if args.force else 'xt'))
          else:
              out = sys.stdout
  
          input_fasta = stack.enter_context(open_file(args.fasta))
          trimmed_reads = trim_fasta(input_fasta, args.sequence, args.num_bases)
  
          SeqIO.write(trimmed_reads, out, "fasta-2line")
  
  if __name__ == '__main__':
      main()
  
  ```
  
  ---
  
  ---
  
  ```
  usage: trim_reads.py [-n INT] [-o FILE] [-f] [-h] FILE STR
  
  input options:
    FILE                  The input FASTA file
    STR                   The sequence to find and trim
  
  trimming options:
    -n INT, --num_bases INT
                          Trim INT bases upstream of input sequence
  
  output options:
    -o FILE, --output FILE
                          Write the filtered output to FILE (default: stdout)
    -f, --force           Overwrite the output file if it exists
  
  additional options:
    -h, --help            Show this help message and exit
  
  ```
  
  
  ---
  
  