Mercurial > repos > bgruening > find_subsequences
comparison find_subsequences.py @ 0:7f39014f9404 draft
Imported from capsule None
| author | bgruening |
|---|---|
| date | Fri, 20 Mar 2015 06:23:17 -0400 |
| parents | |
| children | d882a0a75759 |
comparison
equal
deleted
inserted
replaced
| -1:000000000000 | 0:7f39014f9404 |
|---|---|
| 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 | |
| 13 def find_pattern(seqs, pattern, outfile_path): | |
| 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: | |
| 26 search_func(seq, pattern, outfile) | |
| 27 search_func(seq, rev_compl, outfile, '-') | |
| 28 | |
| 29 | |
| 30 def simple_pattern_search(sequence, pattern, outfile, strand='+'): | |
| 31 """ | |
| 32 Simple regular expression search. This is way faster than the complex search. | |
| 33 """ | |
| 34 bed_template = '%s\t%s\t%s\t%s\t%s\t%s\n' | |
| 35 for match in re.finditer( str(pattern), str(sequence.seq) ): | |
| 36 outfile.write(bed_template % (sequence.id, match.start(), match.end(), sequence.description, '', strand)) | |
| 37 | |
| 38 | |
| 39 def complex_pattern_search(sequence, pattern, outfile, strand='+'): | |
| 40 """ | |
| 41 Searching for pattern with biopyhon's nt_search(). | |
| 42 This allows for ambiguous values, like N = A or T or C or G, R = A or G ... | |
| 43 """ | |
| 44 l = len(pattern) | |
| 45 matches = nt_search(str(sequence.seq), pattern) | |
| 46 bed_template = '%s\t%s\t%s\t%s\t%s\t%s\n' | |
| 47 for match in matches[1:]: | |
| 48 outfile.write(bed_template % (sequence.id, match, match+l, sequence.description, '', strand) ) | |
| 49 | |
| 50 | |
| 51 if __name__ == "__main__": | |
| 52 parser = argparse.ArgumentParser() | |
| 53 parser.add_argument('-i', '--input' , required=True) | |
| 54 parser.add_argument('-o', '--output' , required=True) | |
| 55 parser.add_argument('-p', '--pattern' , required=True) | |
| 56 parser.add_argument('-f', '--format', default="fasta", choices=choices) | |
| 57 args = parser.parse_args() | |
| 58 | |
| 59 with open(args.input) as handle: | |
| 60 find_pattern( SeqIO.parse(handle, args.format), args.pattern, args.output ) | |
| 61 |