--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/find_subsequences.py	Fri Mar 20 06:23:17 2015 -0400
@@ -0,0 +1,61 @@
+#!/usr/bin/env python
+
+import re
+import sys
+import argparse
+from Bio import SeqIO
+from Bio.Seq import Seq
+from Bio.SeqUtils import nt_search
+from Bio.Alphabet import generic_dna
+
+choices = ['embl', 'fasta', 'fastq-sanger', 'fastq', 'fastq-solexa', 'fastq-illumina', 'genbank', 'gb']
+
+def find_pattern(seqs, pattern, outfile_path):
+    """
+    Finds all occurrences of a pattern in the a given sequence.
+    Outputs sequence ID, start and end postion of the pattern.
+    """
+    pattern = pattern.upper()
+    rev_compl = Seq(pattern, generic_dna).complement()
+    search_func = simple_pattern_search
+    if set(pattern).difference(set('ATCG')):
+        search_func = complex_pattern_search
+
+    with open(outfile_path, 'w+') as outfile:
+        for seq in seqs:
+            search_func(seq, pattern, outfile)
+            search_func(seq, rev_compl, outfile, '-')
+
+
+def simple_pattern_search(sequence, pattern, outfile, strand='+'):
+    """
+    Simple regular expression search. This is way faster than the complex search.
+    """
+    bed_template = '%s\t%s\t%s\t%s\t%s\t%s\n'
+    for match in re.finditer( str(pattern), str(sequence.seq) ):
+        outfile.write(bed_template % (sequence.id,  match.start(), match.end(), sequence.description, '', strand))
+
+
+def complex_pattern_search(sequence, pattern, outfile, strand='+'):
+    """
+    Searching for pattern with biopyhon's nt_search().
+    This allows for ambiguous values, like N = A or T or C or G, R = A or G ...
+    """
+    l = len(pattern)
+    matches = nt_search(str(sequence.seq), pattern)
+    bed_template = '%s\t%s\t%s\t%s\t%s\t%s\n'
+    for match in matches[1:]:
+        outfile.write(bed_template % (sequence.id, match, match+l, sequence.description, '', strand) )
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser()
+    parser.add_argument('-i', '--input' , required=True)
+    parser.add_argument('-o', '--output' , required=True)
+    parser.add_argument('-p', '--pattern' , required=True)
+    parser.add_argument('-f', '--format', default="fasta", choices=choices)
+    args = parser.parse_args()
+
+    with open(args.input) as handle:
+        find_pattern( SeqIO.parse(handle, args.format), args.pattern, args.output )
+