1 """
2 nussinov.py
3 NussinovAlgorithm code python
4
5 input.txt:
6 GGGAAAUCG
7
8 c:\workplace> python nussinov.py <input.txt >output.txt
9
10 output.txt (which will be created in the same directory):
11 (2,9)
12 (3,8)
13 (4,7)
14 """
15 """nussinov"""
16
17 class NussinovMatrix:
18 def __init__(self, rnaseq):
19 self.rnaseq = rnaseq
20 self.L = len(rnaseq)
21 self.matrix = []
22 for i in range(self.L):
23 self.matrix.append(range(self.L))
24 for i in range(self.L):
25 for j in range(self.L):
26 self.matrix[i][j] = 'n'
27 self.fillInitialMatrix()
28
29 def fillInitialMatrix(self):
30 for i in range(1, self.L):
31 self.matrix[i][i-1] = 0
32 for i in range(self.L):
33 self.matrix[i][i] = 0
34
35 def fillMatrix(self):
36 endOfSeq = self.L + 1
37 for oneDiagonal in range(1, endOfSeq):
38 i = 1
39 for oneCell in range(1, endOfSeq):
40 j = i + oneDiagonal
41 if j > self.L: break
42 else:
43 self.writeValue(i, j)
44 i += 1
45 endOfSeq -= 1
46
47 def writeValue(self, i, j):
48 unpairi = self.getValue(i+1, j)
49 unpairj = self.getValue(i, j-1)
50 pair = self.getValue(i+1, j-1) + self.score(i, j)
51
52 ks = range(i+1, j)
53 if ks == []: bifurcation = 0
54 else:
55 bifurcations = []
56 for k in ks:
57 case = self.getValue(i, k) + self.getValue(k+1, j)
58 bifurcations.append(case)
59 bifurcation = max(bifurcations)
60 result = max(pair, unpairi, unpairj, bifurcation)
61 self.setValue(i, j, result)
62
63 def score(self, ai, aj):
64 i = self.rnaseq[ai-1]
65 j = self.rnaseq[aj-1]
66 RNAComplementary = {'A':'U', 'U':'A', 'G':'C', 'C':'G'}
67 return RNAComplementary[i] == j
68
69 def getInitialMatrix(self):
70 return self.matrix
71 def getFinalMatrix(self):
72 self.fillMatrix()
73 return self.matrix
74 def getValue(self, i, j):
75 return self.matrix[i-1][j-1]
76 def setValue(self, i, j, aValue):
77 self.matrix[i-1][j-1] = aValue
78
79 class Nussinov(NussinovMatrix):
80 def __init__(self, aRnaseq):
81 NussinovMatrix.__init__(self, aRnaseq)
82 matrix = self.getFinalMatrix()
83 self.stack = []
84 self.stack.append((1, self.L))
85 self.basepair = []
86 self.iterateTraceback()
87
88 def iterateTraceback(self):
89 while 1:
90 try: popedItem = self.stack.pop()
91 except IndexError: break
92 self.doRecursion(popedItem)
93
94 def doRecursion(self, popedItem):
95 i, j = popedItem
96 if i >= j: pass
97 elif self.getValue(i+1, j) == self.getValue(i,j):
98 self.stack.append((i+1, j))
99 elif self.getValue(i, j-1) == self.getValue(i,j):
100 self.stack.append((i, j-1))
101 elif self.getValue(i+1, j-1) + self.score(i,j) == self.getValue(i,j):
102 self.basepair.append((i,j))
103 self.stack.append((i+1, j-1))
104 else:
105 ks = range(i+1, j)
106 for k in ks:
107 if self.getValue(i,k) + self.getValue(k+1,j) == self.getValue(i,j):
108 self.stack.append((k+1,j))
109 self.stack.append((i,k))
110 break
111
112 def getBasePair(self):
113 return self.basepair
114
115 def writeResult():
116 inputRnaTotal = ''
117 while 1:
118 try: inputRna=raw_input()
119 except EOFError: break
120 if not inputRna: break
121 inputRnaTotal += inputRna.strip()
122
123 a = Nussinov(inputRnaTotal)
124 for each in a.getBasePair():
125 print each
126
127 if __name__=='__main__':
128 writeResult()