|
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
|
|
1
|
13 def find_pattern(seqs, pattern, outfile_path, strand):
|
|
0
|
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:
|
|
1
|
26 if strand in ['both', 'forward']:
|
|
|
27 search_func(seq, pattern, outfile)
|
|
|
28 if strand in ['both', 'reverse']:
|
|
|
29 search_func(seq, rev_compl, outfile, '-')
|
|
0
|
30
|
|
|
31
|
|
|
32 def simple_pattern_search(sequence, pattern, outfile, strand='+'):
|
|
|
33 """
|
|
|
34 Simple regular expression search. This is way faster than the complex search.
|
|
|
35 """
|
|
|
36 bed_template = '%s\t%s\t%s\t%s\t%s\t%s\n'
|
|
|
37 for match in re.finditer( str(pattern), str(sequence.seq) ):
|
|
|
38 outfile.write(bed_template % (sequence.id, match.start(), match.end(), sequence.description, '', strand))
|
|
|
39
|
|
|
40
|
|
|
41 def complex_pattern_search(sequence, pattern, outfile, strand='+'):
|
|
|
42 """
|
|
|
43 Searching for pattern with biopyhon's nt_search().
|
|
|
44 This allows for ambiguous values, like N = A or T or C or G, R = A or G ...
|
|
|
45 """
|
|
|
46 l = len(pattern)
|
|
|
47 matches = nt_search(str(sequence.seq), pattern)
|
|
|
48 bed_template = '%s\t%s\t%s\t%s\t%s\t%s\n'
|
|
|
49 for match in matches[1:]:
|
|
|
50 outfile.write(bed_template % (sequence.id, match, match+l, sequence.description, '', strand) )
|
|
|
51
|
|
|
52
|
|
|
53 if __name__ == "__main__":
|
|
|
54 parser = argparse.ArgumentParser()
|
|
1
|
55 parser.add_argument('-i', '--input', required=True)
|
|
|
56 parser.add_argument('-o', '--output', required=True)
|
|
|
57 parser.add_argument('-p', '--pattern', required=True)
|
|
|
58 parser.add_argument('--strand', choices=['both', 'forward', 'reverse'], default='both')
|
|
0
|
59 parser.add_argument('-f', '--format', default="fasta", choices=choices)
|
|
|
60 args = parser.parse_args()
|
|
|
61
|
|
|
62 with open(args.input) as handle:
|
|
1
|
63 find_pattern( SeqIO.parse(handle, args.format), args.pattern, args.output, args.strand )
|
|
0
|
64
|
