text: user16067: 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!
  
  user10775: Here's one way using biopython. Note that the Seq object has a number of methods that act just like those of a Python string. One of these is the find() method, which returns the index of the first occurrence of the sequence if found. Replace with 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
  
  
original_text: 16067: <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>
   10775: <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>
  