简单DNA序列组装（贪婪算法）
Posted orange1002
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了简单DNA序列组装（贪婪算法）相关的知识，希望对你有一定的参考价值。
生物信息学原理作业第四弹：DNA序列组装（贪婪算法）
原理：生物信息学（孙啸）
大致思想：
1. 找到权值最大的边；
2. 除去以最大权值边的起始顶点为起始顶点的边；
3. 除去以最大权值边为终点为终点的边；
4. 重复上述步骤，得到所有符合条件的边；
5. 拼接得到的边；
6. 加入孤立点（如果有）。
附上Python代码，如果有问题我会及时更正（确实不太熟算法）
转载请保留出处！
1 # -*- coding: utf-8 -*-
2 """
3 Created on Mon Dec 4 15:04:39 2017
4 @author: zxzhu
5 python3.6
6 """
7 from functools import reduce
8 def get_weight(s1,s2): #通过两条序列的overlap计算出权值
9 l = min(len(s1),len(s2))
10 while l>0:
11 if s2[:l] == s1[-l:]:
12 return l
13 else:
14 l-=1
15 return 0
16
17 def print_result(s1,s2): #将两条序列去除首尾overlap后合并
18 weight = get_weight(s1,s2)
19 s = s1 + s2[weight:]
20 #print(s)
21 return s
22
23 def dir_graph(l,t=3): #得到满足条件的有向图(权值大于等于t)
24 graph = {}
25 for i in l:
26 for j in l:
27 if i!=j:
28 weight = get_weight(i,j)
29 if weight >= t:
30 key = (i,j)
31 graph[key] = weight
32 return graph
33
34 def rm_path(graph,path): #贪婪算法加入一条边后应该去除与该边首尾顶点相同的边
35 key = graph.keys()
36 rm_key = []
37 for i in key:
38 if i[1] == path[1] or i[0] == path[0]:
39 rm_key.append(i)
40 for i in rm_key:
41 graph.pop(i)
42 return graph
43
44 def get_path(graph,path = []): #得到满足条件的所有边
45 while graph:
46 max_weight = 0
47 for i in graph.keys():
48 if graph[i] > max_weight:
49 max_weight = graph[i]
50 cur_path = i
51 path.append(cur_path)
52 graph = rm_path(graph,cur_path)
53 get_path(graph,path)
54 return path
55
56 def out_num(path,V): #计算某顶点的出度
57 count = 0
58 for i in path:
59 if i[0] == V:
60 count+=1
61 return count
62
63 def get_last_V(path,last_V = None): #得到最后一条边
64 index = 0
65 if last_V: #非随机寻找出度为0的顶点
66 for i in path:
67 if i[1] == last_V:
68 return i,index
69 else:
70 index+=1
71 return None #没有找到指向last_V的顶点(一条路径结束)
72 else: #随机寻找出度为0的顶点
73 for i in path:
74 if out_num(path,i[1]) == 0:
75 return i,index
76 else:
77 index+=1
78 return -1 #首尾相连
79
80
81 def assemble(cur_V,path,new_path = []): #给满足条件的边排序
82 while path:
83 path.pop(cur_V[1])
84 new_path.insert(0,cur_V[0])
85 cur_V = get_last_V(path,last_V = cur_V[0][0])
86 if cur_V:
87 assemble(cur_V,path,new_path)
88 else:
89 cur_V = get_last_V(path)
90 assemble(cur_V,path,new_path)
91 return new_path
92
93 def align_isolated(path,sequence): #加入孤立顶点
94 new_path = reduce(lambda x,y:x+y,path)
95 for i in sequence:
96 if i not in new_path:
97 new_path.append(i)
98 return new_path
99
100 x = \'CCTTTTGG\'
101 y = \'TTGGCAATCACT\'
102 w = \'AGTATTGGCAATC\'
103 u = \'ATGCAAACCT\'
104 z = \'AATCGATG\'
105 v = \'TCACTCCTTTT\'
106 a = w
107 b = y
108 c = \'TCACTAGTA\'
109 sequence = [x,y,w,u,z]
110 sequence1 = [a,b,c]
111 graph = dir_graph(sequence1,t=3)
112 print(graph)
113 path = get_path(graph)
114 path = [list(i) for i in path] #将path中的tuple元素换成list
115 #print(path)
116 start = get_last_V(path) #起始出度为0的顶点所在的边
117 if start == -1: #序列首尾相连
118 new_path = reduce(lambda x,y:x+y, path)
119 new_path = new_path[:-1]
120 result = reduce(print_result,new_path)
121 else:
122 new_path = assemble(start,path) #排序后的边
123 new_path = align_isolated(new_path,sequence1) #加入孤立顶点
124 #print(new_path)
125 result = reduce(print_result,new_path) #组装
126 #print(new_path)
127 print(result)
以上是关于简单DNA序列组装（贪婪算法）的主要内容，如果未能解决你的问题，请参考以下文章