text: user3451: I am trying to figure out the best way to extract sequences from a FASTA file which begin with a common 5' region of 43 nucleotides. Preferably, I would like to to allow for "fuzziness" in this region to allow for mutations or read overlaps. The idea is that any sequences that begin with this region, regardless of size, could be extracted into a new file.
  Example of the FASTA file:
  >1-69050-454.08
  GTACGGGGAAGGACGTCAATAGTC
  
  >2-65989-433.95
  AATCTGTAGTACGACTCACTATAGC
  
  >3-62181-408.91
  GATCTGTAATACGACTCACTATAGG
  
  >4-49959-328.53
  GGGGAAGGACGTCAATAGTCACAC
  
  And what I would like to get from the code is such:
  >2-65989-433.95
  AATCTGTAGTACGACTCACTATAGC
  
  >3-62181-408.91
  GATCTGTAATACGACTCACTATAGG
  
  The FASTA file is quite large and I have tried using a couple grep and awk methods retrieve the sequences. Any help you could provide is much appreciated.
  
  user3381: Here is a 'quick' python3 solution, requiring installation of BLASTN and a pip3 install biopython.  The results are output in FASTA format to stdout.  
  Save the following python3 code in hit_finder.py
  #!/usr/bin/env python3
  
  import sys
  from Bio import SeqIO
  from Bio.Blast.Applications import NcbiblastnCommandline
  
  cline = NcbiblastnCommandline(query=sys.argv[1],
                                subject=sys.argv[2],
                                outfmt="'6 qseqid length qstart'",
                                max_hsps=1,
                                word_size=4)
  
  result = [j[0] for j in [i.split('\t') for i in
                           cline()[0].rstrip().split('\n')] 
            if int(j[1]) >= int(sys.argv[3]) and int(j[2]) <= int(sys.argv[4])]
  
  for record in SeqIO.parse(sys.argv[1], 'fasta'):
      if record.id in result:
          print(record.format('fasta'))
  
  Make it executable: chmod +x hit_finder.py.
  This script takes four arguments (which in future versions could be enhanced using argparse instead of sys.argv):
  arg1 = query_fasta_file (file to search in)
  arg2 = subject_fasta_file (database file to contain the sequence to match at 5')
  arg3 = min_alignment_length (this means the length of the hit, allows mismatching in string with blastn and a word_size of 4)
  arg4 = start_of_hit_at_5prime (where should the alignment hit start?)
  
  Let's say you wanted to search for AATCTGTAGTACGACTCACTATAGC in your query seqeuences.  Save this sequence in subject.fa as:
  >subject
  AATCTGTAGTACGACTCACTATAGC
  
  Your query sequences are stored in queries.fa as:
  >1-69050-454.08
  GTACGGGGAAGGACGTCAATAGTC
  
  >2-65989-433.95
  AATCTGTAGTACGACTCACTATAGC
  
  >3-62181-408.91
  GATCTGTAATACGACTCACTATAGG
  
  >4-49959-328.53
  GGGGAAGGACGTCAATAGTCACAC
  
  Now, to search for hits of at least 23 bp starting no further than 1 base from the 5' end of your sequence, you would run this search using:
  ./hit_finder.py queries.fa subject.fa 23 1
  
  To search for hits of at least 30 bp no more than 10 bp from the 5' end of your queries, do:
  ./hit_finder.py queries.fa subject.fa 30 10
  
  To save your results, output stdout to file, e.g.,:
  ./hit_finder.py queries.fa subject.fa 30 10 > results.fa
  
  
original_text: 3451: <p>I am trying to figure out the best way to extract sequences from a FASTA file which begin with a common 5' region of 43 nucleotides. Preferably, I would like to to allow for "fuzziness" in this region to allow for mutations or read overlaps. The idea is that any sequences that begin with this region, regardless of size, could be extracted into a new file.</p>
  
  <p>Example of the FASTA file:</p>
  
  <pre><code>&gt;1-69050-454.08
  GTACGGGGAAGGACGTCAATAGTC
  
  &gt;2-65989-433.95
  AATCTGTAGTACGACTCACTATAGC
  
  &gt;3-62181-408.91
  GATCTGTAATACGACTCACTATAGG
  
  &gt;4-49959-328.53
  GGGGAAGGACGTCAATAGTCACAC
  </code></pre>
  
  <p>And what I would like to get from the code is such:</p>
  
  <pre><code>&gt;2-65989-433.95
  AATCTGTAGTACGACTCACTATAGC
  
  &gt;3-62181-408.91
  GATCTGTAATACGACTCACTATAGG
  </code></pre>
  
  <p>The FASTA file is quite large and I have tried using a couple grep and awk methods retrieve the sequences. Any help you could provide is much appreciated.</p>
   3381: <p>Here is a 'quick' <code>python3</code> solution, requiring installation of <a href="https://www.ncbi.nlm.nih.gov/books/NBK52640/" rel="nofollow noreferrer">BLASTN</a> and a <code>pip3 install biopython</code>.  The results are output in FASTA format to <code>stdout</code>.  </p>
  
  <p>Save the following <code>python3</code> code in <code>hit_finder.py</code></p>
  
  <pre><code>#!/usr/bin/env python3
  
  import sys
  from Bio import SeqIO
  from Bio.Blast.Applications import NcbiblastnCommandline
  
  cline = NcbiblastnCommandline(query=sys.argv[1],
                                subject=sys.argv[2],
                                outfmt="'6 qseqid length qstart'",
                                max_hsps=1,
                                word_size=4)
  
  result = [j[0] for j in [i.split('\t') for i in
                           cline()[0].rstrip().split('\n')] 
            if int(j[1]) &gt;= int(sys.argv[3]) and int(j[2]) &lt;= int(sys.argv[4])]
  
  for record in SeqIO.parse(sys.argv[1], 'fasta'):
      if record.id in result:
          print(record.format('fasta'))
  </code></pre>
  
  <p>Make it executable: <code>chmod +x hit_finder.py</code>.</p>
  
  <p>This script takes four arguments (which in future versions could be enhanced using <code>argparse</code> instead of <code>sys.argv</code>):</p>
  
  <pre><code>arg1 = query_fasta_file (file to search in)
  arg2 = subject_fasta_file (database file to contain the sequence to match at 5')
  arg3 = min_alignment_length (this means the length of the hit, allows mismatching in string with blastn and a word_size of 4)
  arg4 = start_of_hit_at_5prime (where should the alignment hit start?)
  </code></pre>
  
  <p>Let's say you wanted to search for <code>AATCTGTAGTACGACTCACTATAGC</code> in your query seqeuences.  Save this sequence in <code>subject.fa</code> as:</p>
  
  <pre><code>&gt;subject
  AATCTGTAGTACGACTCACTATAGC
  </code></pre>
  
  <p>Your query sequences are stored in <code>queries.fa</code> as:</p>
  
  <pre><code>&gt;1-69050-454.08
  GTACGGGGAAGGACGTCAATAGTC
  
  &gt;2-65989-433.95
  AATCTGTAGTACGACTCACTATAGC
  
  &gt;3-62181-408.91
  GATCTGTAATACGACTCACTATAGG
  
  &gt;4-49959-328.53
  GGGGAAGGACGTCAATAGTCACAC
  </code></pre>
  
  <p>Now, to search for hits of at least 23 bp starting no further than 1 base from the 5' end of your sequence, you would run this search using:</p>
  
  <pre><code>./hit_finder.py queries.fa subject.fa 23 1
  </code></pre>
  
  <p>To search for hits of at least 30 bp no more than 10 bp from the 5' end of your queries, do:</p>
  
  <pre><code>./hit_finder.py queries.fa subject.fa 30 10
  </code></pre>
  
  <p>To save your results, output <code>stdout</code> to file, e.g.,:</p>
  
  <pre><code>./hit_finder.py queries.fa subject.fa 30 10 &gt; results.fa
  </code></pre>
  