<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>

CC BY-SA 4.0
2018-10-09T17:36:13.970Z
This is exactly what I was looking for, I appreciate your assistance and thank you for taking your time to write this out.
hunter92
CC BY-SA 4.0
2018-10-05T03:02:39.043Z
2018-10-08T08:07:51.670Z
2018-10-08T08:07:51.670Z
user3479780
<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>

CC BY-SA 4.0
2018-10-04T00:58:43.780Z
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?
Bioathlete
<ngs><fasta><awk>
Extracting sequences from FASTA beginning with common 5' end
# 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 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.  

---

