-
Notifications
You must be signed in to change notification settings - Fork 3
/
BA1J.py
111 lines (99 loc) · 3.41 KB
/
BA1J.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
import rna
def hammingDistance(p, q):
ham = 0
for x, y in zip(p, q):
if x != y:
ham += 1
return ham
def immediateNeighbors(pattern):
neighborhood = [pattern]
for i in range(len(pattern)):
symbol = pattern[i]
if symbol != 'A':
neighborhood.append(pattern[0:i] + 'A' + pattern[i+1:])
if symbol != 'C':
neighborhood.append(pattern[0:i] + 'C' + pattern[i+1:])
if symbol != 'G':
neighborhood.append(pattern[0:i] + 'G' + pattern[i+1:])
if symbol != 'T':
neighborhood.append(pattern[0:i] + 'T' + pattern[i+1:])
return neighborhood
def neighbors(pattern, d):
if d == 0:
return [pattern]
if len(pattern) == 1:
return ['A', 'C', 'G', 'T']
neighborhood = []
sufneigh = neighbors(pattern[1:],d)
for x in sufneigh:
if hammingDistance(pattern[1:],x) < d:
for y in ['A', 'C', 'G', 'T']:
neighborhood.append(y + x)
else:
neighborhood.append(pattern[0] + x)
return neighborhood
def patternToNumber(pattern):
if len(pattern) == 0:
return 0
return 4 * patternToNumber(pattern[0:-1]) + symbolToNumber(pattern[-1:])
def symbolToNumber(symbol):
if symbol == "A":
return 0
if symbol == "C":
return 1
if symbol == "G":
return 2
if symbol == "T":
return 3
def numberToPattern(x, k):
if k == 1:
return numberToSymbol(x)
return numberToPattern(x // 4, k-1) + numberToSymbol(x % 4)
def numberToSymbol(x):
if x == 0:
return "A"
if x == 1:
return "C"
if x == 2:
return "G"
if x == 3:
return "T"
def approximatePatternCount(text, pattern, d):
count = 0
for i in range(len(text) - len(pattern) +1):
test = text[i:i+len(pattern)]
if hammingDistance(test, pattern) <= d:
count = count+1
return count
def frequencyArrayWithMismatches(text, k, d):
frequencyArray = [0] * 4**k
close = [0] * 4**k
for i in range(len(text)-k):
neighborhood = neighbors(text[i:i+k],d)
for x in neighborhood:
close[patternToNumber(x)] = 1
for i in range(4**k):
if close[i] == 1:
pattern = numberToPattern(i,k)
frequencyArray[i] = approximatePatternCount(text, pattern, d)
return frequencyArray
def frequentWordsWithMismatches(text, k, d):
frequentPatterns = []
frequent = frequencyArrayWithMismatches(text, k, d)
reverseFrequent = frequencyArrayWithMismatches(rna.reverseComplement(text), k, d)
frequencyArray = [0] * 4**k
for i in range(4**k):
frequencyArray[i] = frequent[i] + reverseFrequent[i]
maxCount = max(frequencyArray)
for i in range(4**k):
if frequencyArray[i] == maxCount:
frequentPatterns.append(numberToPattern(i,k))
return frequentPatterns
text = "ACGCCGAAGCGTGCAAGCGTGCAAGCGTGCAGCTGCAATTCTTTGCCCGGCGTTGCTGCAATCGGCGTTACGCCGATCTTTGCCTCTTTGCCGCTGCAATAGCGTGCACGGCGTTGCTGCAATGCTGCAATCGGCGTTTCTTTGCCTCTTTGCCACGCCGAAGCGTGCATCTTTGCCTCTTTGCCTCTTTGCCCGGCGTTAGCGTGCAGCTGCAATTCTTTGCCAGCGTGCAGCTGCAATACGCCGACGGCGTTCGGCGTTGCTGCAATACGCCGAAGCGTGCAGCTGCAATCGGCGTTACGCCGATCTTTGCCGCTGCAATCGGCGTTTCTTTGCCAGCGTGCAACGCCGAACGCCGAAGCGTGCATCTTTGCCCGGCGTTGCTGCAATACGCCGAAGCGTGCACGGCGTTGCTGCAATGCTGCAATGCTGCAATGCTGCAATACGCCGAGCTGCAATAGCGTGCAACGCCGAACGCCGATCTTTGCCGCTGCAATACGCCGATCTTTGCCTCTTTGCCGCTGCAATACGCCGAAGCGTGCATCTTTGCCGCTGCAATTCTTTGCCACGCCGATCTTTGCCAGCGTGCAGCTGCAATGCTGCAATCGGCGTTACGCCGACGGCGTTGCTGCAATACGCCGATCTTTGCCGCTGCAATCGGCGTTACGCCGAACGCCGAACGCCGAACGCCGAGCTGCAATAGCGTGCATCTTTGCCGCTGCAATACGCCGACGGCGTTCGGCGTTACGCCGACGGCGTTCGGCGTTAGCGTGCACGGCGTTAGCGTGCACGGCGTTAGCGTGCAAGCGTGCACGGCGTTTCTTTGCCTCTTTGCCGCTGCAAT"
k = 6
d = 2
a = frequentWordsWithMismatches(text, k, d)
p = ""
for x in a:
p += x + " "
print(p)