5
$\begingroup$

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.

$\endgroup$
  • 4
    $\begingroup$ 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? $\endgroup$ – Bioathlete Oct 4 '18 at 0:58
2
$\begingroup$

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
$\endgroup$
  • 1
    $\begingroup$ This is exactly what I was looking for, I appreciate your assistance and thank you for taking your time to write this out. $\endgroup$ – hunter92 Oct 9 '18 at 17:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.