-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path全局比对算法.py
295 lines (251 loc) · 11.2 KB
/
全局比对算法.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
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
# -*- coding: utf-8 -*-
"""
全局最优算法的实现尝试
有一个A,C,G,T(U),-的五阶矩阵
有两个长度不超过三的碱基序列
PPT上的写法有点反人类，行列在这个作业里还是正常点好了
使得程序可以满足如
输入AGG AGC
可以得到AGG- / -AGC
和AAG- / A-GC
两组结果
暂时不考虑输入格式异常的问题
chris will
xiaozhi
2019/4/5
"""
import numpy as np
import copy
dic={"A":1,"C":2,"G":3,"T":4,"U":4}
dicc={0:1,1:4,2:2}
def compute(i,j,s,SS,SF):
#i为行号，j为列号，s为距离矩阵,本函数返回SS矩阵位于i行j列的值
#考虑是构建两个矩阵还是改变当前矩阵中每个元素的值
if(i==0 and j==0):
SF[i][j]=0
return 0
elif(j==0):
SF[i][j]=1
return i*int(s[i][-1])
elif(i==0):
SF[i][j]=2
return j*int(s[-1][j])
else:
#矩阵生成完毕
t1=SS[i-1][j]+int(s[dic[A[i-1]]][-1])#来自上方：1
t2=SS[i-1][j-1]+int(s[dic[A[i-1]]][dic[B[j-1]]])#来自斜方：4
t3=SS[i][j-1]+int(s[-1][dic[B[j-1]]])#来自左方：2
#现在需要找到哪个是最大的，是否有一样大的
#三个一样大/两个一样大/一个最大
mn=max(t1,t2,t3)
flag=[i for i,j in enumerate([t1,t2,t3]) if j==mn]
#print(flag,[t1,t2,t3])
if(len(flag)==1):
#只有一条路径
SF[i][j]=dicc[flag[0]]
else:
#有多条路径
temp=0
for k in range(0,len(flag)):
temp+=dicc[flag[k]]
SF[i][j]=temp
return mn
#生成目标矩阵，定义结构体推到方向分为三个，斜向下，向右，向下，分别附权重(4,2,1)
def job(A,B):
#该函数结束之后可以得到SS和SF即全局比对的结果和替代方向矩阵
if(len(A)==len(B)):
#满足最基本的条件，其实好像不用
SS=np.eye(len(A)+1)#SS矩阵
SF=np.eye(len(A)+1)#SF矩阵，用来记录演变方向
for i in range(0,len(A)+1):
for j in range(0,len(B)+1):
SS[i][j]=compute(i,j,s,SS,SF)
#print(SS," SS")
#print(SF," SF" )
return SS,SF
def multipath(SF,ans,stack,pathstack):
#每个位置的权重，正在处理的是哪一条路线，堆栈的情况
#print(ans,"ans!!")
#print(pathstack," pathstack!!")
#print(stack," stack!")
while(len(stack)>0 and len(pathstack)>0):
pathnum=pathstack.pop()
if(pathnum in ans):
if(len(ans[pathnum])==len(SF)):
#这条路已经走到尽头，堆栈中的是留给下一条路的
multipath(SF,ans,stack,pathstack)
break
#开始操作
#每次循环后匹配链长度加一
temp1=stack.pop()
#print(temp1,": temp1")
#获取当前这个元素，判断前进方向
#temp1是含有三个元素的列表，分别为行号，列号，和属性值
if(pathnum not in ans):
#初始化每一条匹配链
ans[pathnum]=[]
if(temp1[2]==1):
#向上,且没有第二条分支
ans[pathnum].append((A[temp1[0]-1],'-'))
stack.append([temp1[0]-1,temp1[1],SF[temp1[0]-1][temp1[1]]])
pathstack.append(pathnum)
multipath(SF,ans,stack,pathstack)
elif(temp1[2]==2):
#向左，且没有第二条分支
ans[pathnum].append(('-',B[temp1[1]-1]))
stack.append([temp1[0],temp1[1]-1,SF[temp1[0]][temp1[1]-1]])
pathstack.append(pathnum)
multipath(SF,ans,stack,pathstack)
elif(temp1[2]==4):
#向斜，且没有第二条分支
ans[pathnum].append((A[temp1[0]-1],B[temp1[1]-1]))
stack.append([temp1[0]-1,temp1[1]-1,SF[temp1[0]-1][temp1[1]-1]])
pathstack.append(pathnum)#俗话说的一条路走到黑
multipath(SF,ans,stack,pathstack)
else:
#多路分支,每个分支先走权重小的,且只要分路，路线条数必增加
if(temp1[2]==3):
#上/左
ans[pathnum+1]=copy.deepcopy(ans[pathnum])
#先拷贝,确保分路前之前的路径内容都一样
ans[pathnum].append(('-',B[temp1[0]-1]))#大的,应为左
ans[pathnum+1].append((A[temp1[0]-1],'-'))
stack2=copy.deepcopy(stack)#分两条路
stack.append([temp1[0]-1,temp1[1],SF[temp1[0]-1][temp1[1]]])
stack.append([temp1[0],temp1[1]-1,SF[temp1[0]][temp1[1]-1]])
pathstack.append(pathnum+1)
pathstack.append(pathnum)
multipath(SF,ans,stack,pathstack)
#multipath(SF,ans,stack2,pathstack)
elif(temp1[2]==5):
#上/斜
ans[pathnum+1]=copy.deepcopy(ans[pathnum])
#先拷贝
#print(ans,"haha")
ans[pathnum].append((A[temp1[0]-1],B[temp1[1]-1]))
#print(ans,"haha1")
ans[pathnum+1].append((A[temp1[0]-1],'-'))
#print(ans,"haha2")
#stack2=copy.deepcopy(stack)#分两条路,但实际上发现拷贝出错
pathstack.append(pathnum+1)
pathstack.append(pathnum)#大的先走，先走就走老路
stack.append([temp1[0]-1,temp1[1],SF[temp1[0]-1][temp1[1]]])
stack.append([temp1[0]-1,temp1[1]-1,SF[temp1[0]-1][temp1[1]-1]])
multipath(SF,ans,stack,pathstack)
#multipath(SF,ans,stack2,pathstack)
elif(temp1[2]==6):
#左/斜
ans[pathnum+1]=copy.deepcopy(ans[pathnum])
#先拷贝
#print(ans,"haha")
ans[pathnum].append((A[temp1[0]-1],B[temp1[1]-1]))
#print(ans,"haha1")
ans[pathnum+1].append(('-',B[temp1[1]-1]))
#print(ans,"haha2")
#stack2=copy.deepcopy(stack)#分两条路,但实际上发现拷贝出错
pathstack.append(pathnum+1)
pathstack.append(pathnum)#大的先走，先走就走老路
stack.append([temp1[0],temp1[1]-1,SF[temp1[0]][temp1[1]-1]])
stack.append([temp1[0]-1,temp1[1]-1,SF[temp1[0]-1][temp1[1]-1]])
multipath(SF,ans,stack,pathstack)
#multipath(SF,ans,stack2,pathstack)
else:
#三路分支
if(temp1[2]==0):
pathstack.append(pathnum)
multipath(SF,ans,stack,pathstack)
else:
#三条路开始
ans[pathnum+1]=copy.deepcopy(ans[pathnum])
ans[pathnum+2]=copy.deepcopy(ans[pathnum])
#先拷贝
#print(ans,"haha")
ans[pathnum].append((A[temp1[0]-1],B[temp1[1]-1]))#斜着
#print(ans,"haha1")
ans[pathnum+1].append(('-',B[temp1[1]-1]))#向左
ans[pathnum+2].append((A[temp1[0]-1],B[temp1[1]]))#向上
#print(ans,"haha2")
#stack2=copy.deepcopy(stack)#分两条路,但实际上发现拷贝出错
pathstack.append(pathnum+2)
pathstack.append(pathnum+1)
pathstack.append(pathnum)#大的先走，先走就走老路
stack.append([temp1[0]-1,temp1[1],SF[temp1[0]-1][temp1[1]]])
stack.append([temp1[0],temp1[1]-1,SF[temp1[0]][temp1[1]-1]])
stack.append([temp1[0]-1,temp1[1]-1,SF[temp1[0]-1][temp1[1]-1]])
multipath(SF,ans,stack,pathstack)
#multipath(SF,ans,stack2,pathstack)
#print(ans,"ans")
def findpath(SS,SF):
#开始处理SF，从SF矩阵最后一个元素开始判断
#且每个位置如果是1,2,4则说明只有一条匹配路线，该点匹配之后则跳过
#如果是3,5,6则是两个点共同的结果要分两条路，如果是7则是三条路替代后的结果
#每个结果输出两个字符串，即匹配结果匹配的长度应为碱基序列长度加一
row=len(SF)#当前主要判断元素的行号+1
col=len(SF[0])
stack=[]#每条路向下走的时候压栈
#先判断第一个点位，并将数据导入堆栈。然后进行判断
pathstack=[]
if(SF[row-1][col-1]==1 or SF[row-1][col-1]==2 or SF[row-1][col-1]==4):
pathstack.append(1)
stack.append([row-1,col-1,SF[row-1][col-1]])
elif(SF[row-1][col-1]==3):
pathstack.append(2)
pathstack.append(1)
stack.append([row-1,col-1,1])
stack.append([row-1,col-1,2])
elif(SF[row-1][col-1]==5):
pathstack.append(2)
pathstack.append(1)
stack.append([row-1,col-1,1])
stack.append([row-1,col-1,4])
elif(SF[row-1][col-1]==6):
pathstack.append(2)
pathstack.append(1)
stack.append([row-1,col-1,2])
stack.append([row-1,col-1,4])
else:
pathstack.append(3)
pathstack.append(2)
pathstack.append(1)
stack.append([row-1,col-1,1])
stack.append([row-1,col-1,2])
stack.append([row-1,col-1,4])
#print(stack)
#最初的哪那一个碱基位置已经被处理
#假定最初进行替换只有一条路径
#计划采用字典来存储不同的路径
#对应节点采用元组来存储
ans={}
multipath(SF,ans,stack,pathstack)
#print(ans,": ans")
#对最后的ans进行处理
final={}
for p,i in enumerate(ans):
#print(ans[i])
temp1=[]
temp2=[]
for k,j in enumerate(ans[i]):
#print(k,j)
temp1.append(j[0])
temp2.append(j[1])
temp1.reverse()
temp2.reverse()#注意反转的方式，否则返回的东西很奇葩
final[p]=["".join(temp1),"".join(temp2)]
print(' {} answers has been found for you!.'.format(len(ans)))
print(final)#finall 即为所求
#先生成矩阵,让矩阵的行列标可以从1开始
s=np.array([
['$','A','C','G','T','-'],
['A','2','-7','-5','-7','-5'],
['C','-7','2','-7','-5','-5'],
['G','-5','-7','2','-7','-5'],
['T','-7','-5','-7','2','-5'],
['-','-5','-5','-5','-5','0']
])
A=input("Please input your first serial which divided by ','.\n")
B=input("Now, input your second serial which divided by ','.\n")
A=A.split(',')
B=B.split(',')
AB=job(A,B)
#已经得到SS,SF
findpath(AB[0],AB[1])