Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-10-09T17:36:13.970Z
Answers.Text: This is exactly what I was looking for, I appreciate your assistance and thank you for taking your time to write this out.
Answers.UserDisplayName: hunter92
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-10-05T03:02:39.043Z
Answers.LastActivityDate: 2018-10-08T08:07:51.670Z
Answers.LastEditDate: 2018-10-08T08:07:51.670Z
Answers.UserDisplayName: user3479780
Body: <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>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-10-04T00:58:43.780Z
Comments.Text: What is the common region that you are expecting to see?  Also do you want to allow for insertions and deletions or just base edits?
Comments.UserDisplayName: Bioathlete
Tags: <ngs><fasta><awk>
Title: Extracting sequences from FASTA beginning with common 5' end
ThreadText: # Extracting sequences from FASTA beginning with common 5' end
  
  (Asked by: username_0 on 2018-10-03)
  
  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.
  
  **Tags:** ngs, fasta, awk
  
  > **Comments:**
  > **username_1**: What is the common region that you are expecting to see?  Also do you want to allow for insertions and deletions or just base edits?  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_2 on 2018-10-05. Score: 2 | Accepted answer
  
  Here is a 'quick' `python3` solution, requiring installation of [BLASTN](https://www.ncbi.nlm.nih.gov/books/NBK52640/) 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
  
  ```
  
  > **Comments:**
  > **username_0**: This is exactly what I was looking for, I appreciate your assistance and thank you for taking your time to write this out.  
  
  ---
  
  