Обратите внимание, что это часть заметки, взятой из этого курса Coursera по алгоритму ДНК.
Проблема: при сборке генома нам дается набор прочтений, которые являются подстроками исходной последовательности, и наша цель — воссоздать исходную последовательность. Например. задано ATCG, TCGC, GCAA => ATCGCGCAA
Один из способов сделать это — предположить, что у нас есть все до единого k-меры (чтения длины k), взятые из чтения. Например. взятие всех 4-меров ATCGCGCAA дает ATCG, TCGC, CGCG, GCGC, CGCA, GCAA. Конечно, в реальном мире нет гарантии, что предположение верно, но это хорошая отправная точка.
График Де Брёйна – удобный способ реконструкции исходной последовательности. Чтобы построить такой граф, мы должны выполнить следующие шаги:
Шаг 1: взять все (k-1)-меры из набора k-меров, например ATCG, TCGC -> ATC, TCG, TCG, CGC. У нас должен быть удвоенный размер чтения k-mer.
Шаг 2: построить мультиграф с узлами, являющимися k-1-мерами; провести ребро между двумя мерами k-1 только в том случае, если два меров k-1 взяты из одного и того же чтения. Например. ATCG, TCGC дает
Шаг 3: граф, построенный таким образом, гарантированно имеет эйлеров след, мы идем по следу и соединяем узлы, чтобы сформировать нашу исходную последовательность.
Подробнее: почему алгоритм работает? Второй k-1-мер каждого k-мера является началом k-1-мера в следующем k-мере, продолжаясь таким образом до конца последовательности.
Когда алгоритм не сработает? когда у нас есть повторы в последовательности. Например, ниже приведен график, полученный путем взятия всех 3-меров ATGATATG.
В этом случае мы могли бы построить TGATGATG из графа, выполнив эйлеров цикл (точнее, контур), который отличается от исходной последовательности. Чтобы преодолеть двусмысленность из-за коротких повторов, мы можем увеличить длину чтения, чтобы у нас было больше различных k-1-меров.
На пути к реальному алгоритму сборки
Мы видели алгоритм для сборки чтений k-mer. Мы можем ослабить условие, чтобы принимать чтения выше заданной длины, и разбить каждый k-мер на (k-n)-мер, чтобы учитывать чтения, которые имеют k-n перекрытий вместо k-1 перекрытия.
Кроме того, мы можем попытаться реконструировать части, которые легче собрать (контиги), и исключить неоднозначные части. Например, ATGGCCGAAGCAACGGGTTTGTGACGCGCGCGCGCGACCAGG имеет длинные повторы CG в середине, но два конца довольно легко собрать.
Фрагмент кода здесь:
# coding: utf-8
# In[2]:
def debruijnize(reads):
nodes = set()
not_starts = set()
edges = []
for r in reads:
r1 = r[:-1]
r2 = r[1:]
nodes.add(r1)
nodes.add(r2)
edges.append((r1,r2))
not_starts.add(r2)
return (nodes,edges,list(nodes-not_starts))
# In[3]:
def build_k_mer(str,k):
return [str[i:k+i] for i in range(0,len(str)-k+1)]
# In[4]:
reads = build_k_mer("ATCGTTGCGCGACCG",4)
print(reads)
# In[5]:
G = debruijnize(reads)
print(G)
# In[6]:
def make_node_edge_map(edges):
node_edge_map = {}
for e in edges:
n = e[0]
if n in node_edge_map:
node_edge_map[n].append(e[1])
else:
node_edge_map[n] = [e[1]]
return node_edge_map
m = make_node_edge_map(G[1])
print(m)
# In[7]:
def eulerian_trail(m,v):
nemap = m
result_trail = []
start = v
result_trail.append(start)
while(True):
trail = []
previous = start
while(True):
if(previous not in nemap):
break
next = nemap[previous].pop()
if(len(nemap[previous]) == 0):
nemap.pop(previous,None)
trail.append(next)
if(next == start):
break;
previous = next
# completed one trail
print(trail)
index = result_trail.index(start)
result_trail = result_trail[0:index+1] + trail + result_trail[index+1:len(result_trail)]
# choose new start
if(len(nemap)==0):
break
found_new_start = False
for n in result_trail:
if n in nemap:
start = n
found_new_start = True
break # from for loop
if not found_new_start:
print("error")
print("result_trail",result_trail)
print(nemap)
break
return result_trail
start = G[2][0] if (len(G[2]) > 0) else G[0][0]
print m
t = eulerian_trail(m,start)
print(t)
# In[8]:
def visualize_debruijn(G):
nodes = G[0]
edges = G[1]
dot_str= 'digraph "DeBruijn graph" {\n '
for node in nodes:
dot_str += ' %s [label="%s"] ;\n' %(node,node)
for src,dst in edges:
dot_str += ' %s->%s;\n' %(src,dst)
return dot_str + '}\n'
# In[9]:
get_ipython().magic(u'load_ext gvmagic')
# In[10]:
get_ipython().magic(u'dotstr visualize_debruijn(G)')
# In[11]:
def assemble_trail(trail):
if len(trail) == 0:
return ""
result = trail[0][:-1]
for node in trail:
result += node[-1]
return result
assemble_trail(t)
# In[12]:
def test_assembly_debruijn(t,k):
reads = build_k_mer(t,k)
G = debruijnize(reads)
v = visualize_debruijn(G)
nemap = make_node_edge_map(G[1])
print(G)
print(v)
start = next(iter(G[2])) if (len(G[2]) > 0) else next(iter(G[0]))
trail = eulerian_trail(nemap,start)
return assemble_trail(trail)
# In[17]:
test_assembly_debruijn("ATGGCCGAAGCAACGGGTTTGTGACGCGCGCGCGCGACCAGG",4)
