Mercurial > repos > bgruening > find_subsequences
annotate find_subsequences.py @ 8:e614ec7a43b8 draft default tip
planemo upload commit caa6f11d80692d7e03f2224f69a535ddd577746e
| author | bgruening |
|---|---|
| date | Wed, 12 Jul 2017 11:51:17 -0400 |
| parents | 99f356eeba15 |
| children |
| rev | line source |
|---|---|
| 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 | |
| 6 | 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: | |
| 6 | 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' | |
|
7
99f356eeba15
planemo upload commit 7ff6d582464c6b7183e1acd421272167f8a2433f
bgruening
parents:
6
diff
changeset
|
37 for match in re.finditer( str(pattern), str(sequence.seq), re.IGNORECASE ): |
| 4 | 38 outfile.write(bed_template % (sequence.id, match.start(), match.end(), sequence.description, '', strand)) |
| 0 | 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:]: | |
| 4 | 50 outfile.write(bed_template % (sequence.id, match, match+l, sequence.description, '', strand) ) |
| 0 | 51 |
| 52 | |
| 53 if __name__ == "__main__": | |
| 54 parser = argparse.ArgumentParser() | |
| 6 | 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: | |
| 6 | 63 find_pattern( SeqIO.parse(handle, args.format), args.pattern, args.output, args.strand ) |
| 0 | 64 |