|
0
|
1 #!/usr/bin/env python
|
|
|
2
|
|
|
3 import re
|
|
|
4 import sys
|
|
|
5 import argparse
|
|
|
6 from Bio import SeqIO
|
|
|
7 from Bio.Seq import Seq
|
|
|
8 from Bio.SeqUtils import nt_search
|
|
|
9 from Bio.Alphabet import generic_dna
|
|
|
10
|
|
|
11 choices = ['embl', 'fasta', 'fastq-sanger', 'fastq', 'fastq-solexa', 'fastq-illumina', 'genbank', 'gb']
|
|
|
12
|
|
|
13 def find_pattern(seqs, pattern, outfile_path):
|
|
|
14 """
|
|
|
15 Finds all occurrences of a pattern in the a given sequence.
|
|
|
16 Outputs sequence ID, start and end postion of the pattern.
|
|
|
17 """
|
|
|
18 pattern = pattern.upper()
|
|
|
19 rev_compl = Seq(pattern, generic_dna).complement()
|
|
|
20 search_func = simple_pattern_search
|
|
|
21 if set(pattern).difference(set('ATCG')):
|
|
|
22 search_func = complex_pattern_search
|
|
|
23
|
|
|
24 with open(outfile_path, 'w+') as outfile:
|
|
|
25 for seq in seqs:
|
|
|
26 search_func(seq, pattern, outfile)
|
|
|
27 search_func(seq, rev_compl, outfile, '-')
|
|
|
28
|
|
|
29
|
|
|
30 def simple_pattern_search(sequence, pattern, outfile, strand='+'):
|
|
|
31 """
|
|
|
32 Simple regular expression search. This is way faster than the complex search.
|
|
|
33 """
|
|
|
34 bed_template = '%s\t%s\t%s\t%s\t%s\t%s\n'
|
|
|
35 for match in re.finditer( str(pattern), str(sequence.seq) ):
|
|
|
36 outfile.write(bed_template % (sequence.id, match.start(), match.end(), sequence.name, '', strand))
|
|
|
37
|
|
|
38
|
|
|
39 def complex_pattern_search(sequence, pattern, outfile, strand='+'):
|
|
|
40 """
|
|
|
41 Searching for pattern with biopyhon's nt_search().
|
|
|
42 This allows for ambiguous values, like N = A or T or C or G, R = A or G ...
|
|
|
43 """
|
|
|
44 l = len(pattern)
|
|
|
45 matches = nt_search(str(sequence.seq), pattern)
|
|
|
46 bed_template = '%s\t%s\t%s\t%s\t%s\t%s\n'
|
|
|
47 for match in matches[1:]:
|
|
|
48 outfile.write(bed_template % (sequence.id, match, match+l, sequence.name, '', strand) )
|
|
|
49
|
|
|
50
|
|
|
51 if __name__ == "__main__":
|
|
|
52 parser = argparse.ArgumentParser()
|
|
|
53 parser.add_argument('-i', '--input' , required=True)
|
|
|
54 parser.add_argument('-o', '--output' , required=True)
|
|
|
55 parser.add_argument('-p', '--pattern' , required=True)
|
|
|
56 parser.add_argument('-f', '--format', default="fasta", choices=choices)
|
|
|
57 args = parser.parse_args()
|
|
|
58
|
|
|
59 with open(args.input) as handle:
|
|
|
60 find_pattern( SeqIO.parse(handle, args.format), args.pattern, args.output )
|
|
|
61
|
