annotate find_subsequences.py @ 0:7f39014f9404 draft

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