# HG changeset patch
# User bgruening
# Date 1426846997 14400
# Node ID 7f39014f9404ff57eadd8d518d7eaf0eb6e97f96
Imported from capsule None
diff -r 000000000000 -r 7f39014f9404 find_subsequences.py
--- /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 )
+
diff -r 000000000000 -r 7f39014f9404 find_subsequences.xml
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/find_subsequences.xml Fri Mar 20 06:23:17 2015 -0400
@@ -0,0 +1,291 @@
+
+ providing regions in BED format
+
+ biopython
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ ^[a-zA-Z]+$
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+`_.
+
+This tool is searching on both strands.
+
+]]>
+
+
+ 10.1093/bioinformatics/btp163
+
+
diff -r 000000000000 -r 7f39014f9404 test-data/find_subsequences_advanced_result1.bed
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/find_subsequences_advanced_result1.bed Fri Mar 20 06:23:17 2015 -0400
@@ -0,0 +1,2 @@
+forward_advanced 9 15 forward_advanced +
+reverse_advanced 9 15 reverse_advanced -
diff -r 000000000000 -r 7f39014f9404 test-data/find_subsequences_input1.fasta
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/find_subsequences_input1.fasta Fri Mar 20 06:23:17 2015 -0400
@@ -0,0 +1,11 @@
+>forward_simple
+ACCTGTAGCATCGGGGCTCGTA
+
+>reverse_simple
+AAAATTTAGCGATTGTAGTGTAGATAGTA
+
+>forward_advanced
+AAAAAAAAAATGGCGTTTTTTTTTT
+
+>reverse_advanced
+AAAAAAAAATAGGGCTTTTTTTTTTT
diff -r 000000000000 -r 7f39014f9404 test-data/find_subsequences_simple_result1.bed
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/find_subsequences_simple_result1.bed Fri Mar 20 06:23:17 2015 -0400
@@ -0,0 +1,3 @@
+forward_simple 9 13 forward_simple +
+forward_simple 5 9 forward_simple -
+reverse_simple 6 10 reverse_simple -
diff -r 000000000000 -r 7f39014f9404 test-data/find_subsequences_simple_result2.bed
diff -r 000000000000 -r 7f39014f9404 tool_dependencies.xml
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tool_dependencies.xml Fri Mar 20 06:23:17 2015 -0400
@@ -0,0 +1,6 @@
+
+
+
+
+
+