DNA序列组装（贪婪算法）

生物信息学原理作业第四弹：DNA序列组装（贪婪算法）

原理：生物信息学（孙啸）

大致思想：

　　　　　　1. 找到权值最大的边；

　　　　　　2. 除去以最大权值边的起始顶点为起始顶点的边；

　　　　　　3. 除去以最大权值边为终点为终点的边；

　　　　　　4. 重复上述步骤，得到所有符合条件的边；

　　　　　　5. 拼接得到的边；

　　　　　　6. 加入孤立点（如果有）。

附上Python代码，如果有问题我会及时更正（确实不太熟算法）

DNA序列组装（贪婪算法）

转载请保留出处！

  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)

 

 

posted @ 2017-12-04 21:33  orange1002  阅读(2445)  评论(4)  编辑  收藏  举报