生物信息学原理作业第四弹：DNA序列组装（非循环子图）
原理：生物信息学（孙啸）
大致思想：
1. 这个算法理解细节理解比较困难，建议看孙啸的生物信息学相关章节。
2. 算法要求所有序列覆盖整个目标DNA，并保证相邻片段有足够的覆盖连接（引自孙啸 生物信息学）。
3. 最后推导出符合条件的序列构成的有向图没有回路，并有哈密顿路径。
4. 利用拓扑排序，得到顶点的有序排列。
5. 组装。
贴上Python代码，发现问题我会及时更正。
转载请保留出处！
# -*- coding: utf-8 -*-
"""
Created on Sat Dec 2 16:09:14 2017
@author: zxzhu
python3.6
"""
from functools import reduce def get_weight(s1,s2): #通过两条序列的overlap计算出权值
l = min(len(s1),len(s2))
while l>0:
if s2[:l] == s1[-l:]:
return l
else:
l-=1
return 0 def print_result(s1,s2): #将两条序列去除首尾overlap后合并
weight = get_weight(s1,s2)
s = s1 + s2[weight:]
#print(s)
return s def dir_graph(l,t=3): #得到满足条件的有向图(权值大于等于t)
graph = {}
for i in l:
VW = []
for j in l:
if i!=j:
weight = get_weight(i,j)
if weight >= t:
VW.append(j)
graph[i] = VW
#print(graph)
for i in graph.keys(): #不能有孤立顶点
if not graph[i]:
count = get_in_V(graph,i)
if count ==0:
graph.clear()
print('The sequence:\n"{0}"\n can\'t align with others!'.format(i))
break
return graph def get_in_V(graph,v): #得到某顶点入度
count = 0
all_in = reduce(lambda x,y:x+y,graph.values())
for i in all_in:
if i == v:
count+=1
return count def aligner(graph,topo=[]): #得出顶点顺序
while graph:
V = graph.keys()
for i in V:
flag = 1
in_num = get_in_V(graph,i)
if in_num ==0:
topo.append(i)
graph.pop(i)
flag = 0
break
if flag: #存在环
#print('The t score is too small!')
return None
else:
aligner(graph,topo)
return topo x = 'CCTTTTGG'
y = 'TTGGCAATCACT'
w = 'AGTATTGGCAATC'
u = 'ATGCAAACCT'
z = 'AATCGATG'
v = 'TCACTCCTTTT'
graph = dir_graph([x,y,z,w,u],t=3)
topo = aligner(graph)
if topo:
result = reduce(print_result,topo)
else:
result = topo
print(result)