Во многих приложениях биоинформатики используются деревья, то есть связные графы без циклов. В частности, деревья — удобный инструмент для описания взаимоотношений родства между биологическими видами (так называемые филогенетические или эволюционные деревья). Для определения отношений используются морфологические признаки (размер индивидов, цвет волосяного покрова и тому подобное) и / или данные биохимических исследований (например, различающиеся нуклеотиды в генах или, соответственно, аминокислоты в аналогичных белках). Построением филогенетических деревьев занимается отдельная область биоинформатики — вычислительная филогенетика; методы создания деревьев обычно основаны на машинном обучении (например, иерархической кластеризации) или математической статистике (тестирование статистических гипотез, выбор наиболее вероятной вероятностной модели).
В системе Rosalind.info филогенетические деревья используются в нескольких задачах:
- построение дерева на основе признаков (CHBP);
- восстановление генетической информации предков на основе этих данных у потомков (ALPH);
- сравнение двух филогенетических деревьев с помощью расстояния на основе разделений (англ. split) и квартетов (англ. quartet) — SPTD и QRTD, соответственно.
Имеет смысл рассмотреть общие задачи, полезные для их решения:
- построение дерева на основе текстового представления;
- вычисление рекуррентных функций на деревьях;
- сравнение деревьев между собой.
Хранение деревьев
Первая задача, которая возникает при работе с деревьями — их хранение в памяти. В случае Python имеет смысл использовать объектно-ориентированный поход, представив все вершины дерева в виде объектов класса Tree, обладающего следующими полями:
- data — данные, связанные с вершиной, например название биологического вида, связанного с вершиной;
- parent — родительская вершина, None для корня дерева;
- children — дочерние вершины.
Помимо этого, для минимальной функциональности класс Tree содержит методы append() и prepend(), добавляющие поддерево в конец и начало списка дочерних вершин, соответственно.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
class Tree(object): def __init__(self, data): self.parent = None self.data = data self.children = [] def prepend(self, child): if child.parent is not None: child.parent.children.remove(child) self.children[0:0] = [ child ] child.parent = self def append(self, child): if child.parent is not None: child.parent.children.remove(child) self.children.append(child) child.parent = self |
Для текстового представления деревьев в филогенетике часто используется формат Newick. Согласно формату, дерево рекурсивно записывается как ограниченное скобками перечисление дочерних вершин через запятую, за которым могут следовать текст, связанный с корневой вершиной дерева. Запись завершается точкой с запятой. Формально с использованием БНФ-образной записи синтаксиса:
TREE ::= NODE ';'
NODE ::= '(' NODE ( ',' NODE )* ')' data
data ::= [^,();]*
((,A,)B,((C,D),E),(F,G))H;Задачу создания дерева на основе записи можно решить с помощью произвольного генератора LR- или LL-анализаторов (см. список инструментов для синтаксического анализа для Python). С другой стороны, поскольку грамматика Newick-записи настолько проста, анализатор легко пишется вручную; его можно выделить в отдельный метод класса Tree. Разбор удобнее проводить справа налево, с конца текстового представления в начало. Во время разбора хранятся два объекта класса Tree:
- node — последняя созданная вершина дерева;
- parent — вершина, в которую в данный момент добавляются дочерние вершины.
Первоначально node = None, а parent равняется вершине, для которой вызван метод.
На каждом этапе синтаксического разбора создается новая вершина node, после чего читаются связанные с ней данные. Эти данные соответствуют регулярному выражению [^,();]*, то есть корректно обрабатываются вершины без связанных данных. Новая вершина не создается, если последним обработанным символом была открывающаяся скобка, поскольку в этом и только в этом месте Newick-представления вершин быть не может. Последующие действия зависят от последнего необработанного символа текстового представления.
| символ | Действие |
|---|---|
| ) | Начинается список дочерних вершин для последней созданной вершины. parent := node. |
| ( | Список дочерних вершин для текущей родительской вершины закончился; новые вершины будут добавляться в родительскую для нее вершину. parent := parent.parent. |
| , | Игнорировать. |
| ; | Игнорировать. |
Итерации алгоритма продолжаются, пока не достигнуто начало Newick-представления.
Рекуррентные отношения на деревьях
Подобно рекуррентным отношениям, заданным на множестве натуральных чисел, в задачах на деревьях часто используются отношения, заданные на вершинах дерева. Обычно значение некоторой функции F в вершине дерева зависит от значения этой функции в дочерних вершинах v.children (хотя бывают и обратные случаи, когда F(v) определяется величиной F(v.parent)). Одно из таких отношений — преобразование дерева к Newick-представлению, то есть функция, обратная по своей сути синтаксическому разбору.
|
1 2 3 4 5 6 7 8 9 10 |
class Tree(object): ... # другие методы def __str__(self): return 'Tree(\'' + self.newick() + ';\')' def newick(self): s = ','.join([child.newick() for child in self.children]) if len(s) > 0: s = '(' + s + ')' return s + str(self.data) |
Простейший способ вычисления рекуррентных отношений — через рекурсивные вызовы — чреват теми же проблемами, что и, например, рекурсивное вычисление чисел Фибоначчи. Для более эффективного определения рекуррентных функций нужно поступить так же, как в динамическом программировании, — линейно упорядочить последовательность вызовов. Для этого годится обход дерева в ширину, начиная с корня; очередь, создаваемая при таком обходе, содержит все вершины дерева, причем родитель вершины всегда расположен до самой вершины. Для хранения вычисленных значений функции можно использовать словарь. При этом надо использовать функцию Python id(), так как деревья являются изменяемыми и не могут сами использоваться в качестве ключей:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
def newick(tree): queue = [ tree ]; ptr = 0 while ptr < len(queue): queue.extend(queue[ptr].children) ptr += 1 text = dict() for node in reversed(queue): s = ','.join([ text[id(child)] for child in node.children ]) if len(s) > 0: s = '(' + s + ')' s += str(node.data) text[id(node)] = s return text[id(tree)] + ';' |
Код создания очереди стоит выделить в отдельный метод класса Tree (этот метод достаточно полезен сам по себе). Логика прохода по очереди выделяется в декоратор, аналогичный декоратору динамического программирования:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 |
class Tree(object): ... # другие методы def queue(self): queue = [ self ]; ptr = 0 while ptr < len(queue): queue.extend(queue[ptr].children) ptr += 1 return queue def leaves_to_root(fn): cache = dict() def ltr_fn(tree): if id(tree) in cache: # значение функции уже вычислено result = cache[id(tree)] else: queue = tree.queue() # пройти по очереди от листьев к корню for node in reversed(queue): cache[id(node)] = fn(node) result = cache[id(tree)] cache.clear() return result return ltr_fn |
В отличие от декоратора @dp, кэш очищается после каждого вычисления функции, поскольку содержащиеся в нем значения могут устареть (например, если в дерево добавить новые вершины).
С помощью декоратора легко реализовать, например, подсчет количества листьев в дереве:
|
1 2 3 4 5 6 |
@leaves_to_root def nleaves(tree): if len(tree.children) == 0: # дерево состоит из одной вершины return 1 return sum([nleaves(ch) for ch in tree.children]) |
Равенство деревьев
Другой пример задачи на деревьях, которую можно решить при помощи рекуррентных функций — определение его разбиений (англ. split). Согласно определению, разбиение дерева T — это два непересекающихся множества его листьев, которые соответствуют двум поддеревьям, образуемым в результате удаления из T одного из ребер.
Каждое из ребер дерева можно представить в виде пары вершин (v, v.parent). Легко заметить, что одно из двух множеств в разбиении, которое порождается этим ребром, состоит из всех листьев, являющихся потомками v (включая саму вершину v). В дальнейшем будем подразумевать под произвольным разбиением именно это одно множество. Пусть S(v) — множество разбиений в поддереве с корнем в v, включая «полное» разбиение — множество всех листьев в этом поддереве. Для функции S существует рекуррентное отношение. В самом деле, множество разбиений для некоторой вершины равно объединению разбиений (включая полные) для дочерних вершин, плюс новое полное разбиение:
(Через max обозначен элемент S(u) с максимальным числом элементов, то есть «полное» разбиение для этой вершины.)
Значит, все разбиения можно вычислить за один проход по вершинам дерева:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 |
class Tree(object): ... # другие методы @leaves_to_root def all_splits(self): # вместо множеств используются списки, так как # деревья не хэшируются if len(self.children) == 0: return [ [ self ] ] splits = [] all_leaves = [] for child in self.children: child_splits = child.all_splits() splits += child_splits # последний элемент child_splits — полное разбиение all_leaves += child_splits[-1] splits.append(all_leaves) return splits |
Разбиения полезны для определения равенства между деревьями. Одно и то же дерево можно записать в формате Newick несколькими способами, в зависимости от выбранной в качестве корня вершины и порядка перечисления дочерних вершин. Например, приведенное выше на рисунке дерево можно записать как
((,A,)B,((C,D),E),(F,G))H; (,A,,(((C,D),E),(F,G))H)B; (C,D,(E,((F,G),(,A,)B)H)); ((G,((,A,)B,(E,(C,D)))H))F;
и еще кучей способов. Известно, что два дерева равны тогда и только тогда, когда равны их разбиения; значит, операцию равенства для класса Tree можно определить так:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
class Tree(object): ... # другие методы def __eq__(self, other): if not isinstance(other, Tree): return False to_data = lambda nodes: frozenset([ v.data for v in nodes ]) splits, osplits = map(to_data, self.all_splits()), \ map(to_data, other.all_splits()) leaves = splits[-1] # листья этого дерева for split in splits: if (split not in osplits) and \ (leaves - split not in osplits): return False return True |
(В приведенном ниже сценарии на Github используется немного более быстрая проверка равенства с использованием разбиений, содержащих более одного элемента.)