Методы построения деревьев

Основной задачей филогенетического анализа является восстановление истинного филогенетического дерева. Эта задача включает в себя подбор топологии и оценивание длин ветвей дерева.

Статистические методы для построения филогенетического дерева по данным молекулярной генетики можно разделить на три группы:

1) методы расстояний;

2) методы парсимонии;

3) методы максимального правдоподобия.

В этой группе существует несколько методов для построения филогении. В них для каждой пары изучаемых таксонов или других форм вычисляются эволюционные расстояния, и по ним строится филогенетическое дерево.

Обобщенного статистического критерия для выбора подходящей меры расстояний не существует. Но в работе М. Нея, С. Кумара (2004) сформулированы эмпирические правила – рекомендации, которым можно следовать при построении филогении:

1. Если d – оценка числа нуклеотидных замен на сайт, определенная по модели Джукса – Кантора меньше или равна 5%, то, вне зависимости от наличия предпочтения транзиций/трансверсий и (или) разной скорости замен для разных сайтов, следует использовать расстояние р или модель Джукса – Кантора. Расстояние р особенно подходит для коротких последовательностей.

2. Если 5% < d < 30% и отношение числа транзиций/трансверсий (R) не слишком велико (R < 5), то для достаточно длинных последовательностей можно использовать расстояние, определенное по модели Джукса – Кантора. Если доли транзиций и трансверсий существенно различны, то для достаточно длинных последовательностей следует использовать расстояние, определенное по модели Кимуры.

3. Если 30% < d < 100% и нуклеотидные частоты (А, Т, С, G) значительно различны, но нет сильно выраженных предпочтений транзиций/трансверсий, следует использовать расстояние, определенное по модели Таджимы – Нея. Если есть выраженные предпочтения транзиций/трансверсий и содержания GC, следует использовать расстояние, определенное по модели Тамуры – Нея.

4. Если d > 100% для большого числа попарных сравнений, восстановленное дерево филогении обычно не достоверно по целому ряду причин (например, большая дисперсия и ошибки выравнивания последовательностей). Такие данные лучше не использовать. Хотя можно попробовать не рассматривать фрагменты, эволюционирующие с повышенной скоростью, и проанализировать только медленно эволюционирующую часть данных; так часто поступают с вариабельным районом генов иммуноглобулина. Также можно попытаться найти другие гены, эволюционирующие медленнее.

5. Многие меры расстояний для оценки числа нуклеотидных замен на сайт (d) нельзя применять, когда последовательности находятся на больших расстояниях. Это происходит потому, что формула для вычисления расстояния, как правило, содержит логарифм, и аргумент логарифма становится отрицательным. Следовательно, лучше избегать строить филогению по сильно дивергировавшим последовательностям, а если это все же необходимо, то в качестве меры расстояния лучше использовать расстояние р с меньшей дисперсией и без логарифма.

6. Для эволюционно удаленных видов лучше использовать расстояние между аминокислотами.
Невзвешенный парно–групповой метод расстояний.

Самый простой метод в этой группе – это невзвешенный метод попарной группировки с использованием среднего арифметического (НПГМ или UPGMA).

В этом методе сначала для каждой пары таксонов (i, j), согласно выбранной модели замен, определяются эволюционные расстояния (d_ij). Затем используют кластерный анализ по d_ij, алгоритм которого состоит в следующем.

Два таксона, имеющие минимальное расстояние объединяют, создавая новый узел (кластер). Расстояние между любыми двумя кластерами (А и В) вычисляется по формуле

где r и s – число таксонов в кластерах А и В, a d_ij –расстояние между таксоном i кластера А и таксоном j кластера В. Точка ветвления (уровень, время расхождения) двух кластеров вычисляется как d_AB / 2. После этого матрица расстояний пересчитывается, и процесс кластеризации повторяется, пока не произойдет объединение всех таксонов. В результате получаем дерево с корнем.

Рассмотрим алгоритм этого метода на модельном примере для пяти таксонов.

На первом этапе в методе UPGMA строится матрица попарных эволюционных расстояний (d) в процентах между каждой парой i, j (табл 6.2.).

Таксон	1	2	3	4
2	d12=5
3	d13=7	d23=8
4	d14=10	d24=12	d34=11
5	d15=13	d25=15	d35=15	d45=13

Кластеризация таксонов начинается с пары, имеющей наименьшее расстояние. Для рассмотренной матрицы наименьшее расстояние d₁₂=5. Тогда объединяются таксоны 1 и 2, а точка ветвления (узел) для них находится на уровне b = d₁₂ / 2=5 / 2=2,5. Предполагается, что таксоны 1 и 2 одинаково удалены от точки ветвления. Таким образом, таксоны 1 и 2 объединяются в один составной таксон, или кластер. Обозначим его u [и = (1 – 2)].

Матрица расстояний пересчитывается по формуле (6.1). Расстояние между этим кластером и и каждым из остальных таксонов с номером k (k ≠ 1, 2) вычисляется как d_uk= (d_1k + d_2k) / 2. В частности: d_u3= (d₁₃ + d₂₃) / 2=(7+8) / 2=7,5; d_u4= (d₁₄ + d₂₄) / 2=(10+12) / 2=11; d_u5= (d₁₅ + d₂₅) / 2=(13+15) / 2=14. Получаем новую матрицу попарных эволюционных расстояний (табл. 6.3.):
Таблица 6.3. Матрица расстояний для 4–х кластеров.

Таксон	u=(1–2)	3	4
3	du3=7,5
4	du4=11	d34=11
5	du5=14	d35=15	d45=13

Теперь наименьшее расстояние в матрице d_u3=7,5. Тогда в новый кластер [v = (1 – 2 – 3)] объединяются таксоны и (1, 2) и 3, а точка ветвления вычисляется как b = d_u3/2 = [(d₁₃ + d₂₃)/2] / 2=7,5 / 2=3,75. Расстояние между кластером v и оставшимися таксонами (k) пересчитывается: d_v4 = (d₁₄ + d₂₄ + d₃₄) / 3= (10+12+11) / 3=11; d_v5 = (d₁₅ + d₂₅ + d₃₅) / 3= (13+15+15)=14,33.

Получаем новую матрицу (табл. 6.4.)
Таблица 6.4. Матрица расстояний для 3–х кластеров.

Таксон	u=(1–2–3)	4
4	dv4=11
5	dv5=14,33	d45=13

Наименьшее расстояние d_v4=11. Поэтому объединяем v = (1–2–3) и таксон 4 с точкой ветвления b = d_v4/2 = [(d₁₄ + d₂₄ + d₃₄) /3] / 2)=11/2=5,5. Последний таксон 5 – объединяется со всеми остальными, и точка ветвления вычисляется как b = [(d_l5 + d₂₅ + d₃₅ + d₄₅) / 4] / 2)=(13+15+15+13) / 8=7.

Расстояние между кластером и (1 – 2) и таксоном 3 равно 7,5. Взяв половину этого значения 7,5 / 2=3,75, получаем время до их общего предка. Расстояние до кластера u (1–2), равно 3,75 – 2,5=1,25.

Расстояние от таксона 4 до кластера v (1 – 2 – 3) равно 5,5 – 3,75=1,75.

Расстояние от таксона 5 до кластера, объединяющего 4 таксона (1–2–3–4) равно 7 – 5,5=1,5.

Полученное дерево представлено на рис 6.3.

2,5

1

Рис 6.3. Дерево, построенное по методу UPGMA с указанием длин ветвей.

Для деревьев, восстановленных методом UPGMA, подразумевается аддитивность длин ветвей. Длины ветвей данного дерева называются аддитивными, если расстояние между любыми двумя таксонами равно сумме длин ветвей на пути, соединяющем эти таксоны.

Так как построение дерева основано на ограниченном наборе данных о последовательностях нуклеотидов, то важно оценить достоверность построенного дерева. Для этого применяется тест бутстреп, о котором будет сказано ниже. Метод UPGMA может использоваться для анализа молекулярных данных, если скорость нуклеотидных замен для гена более–менее постоянна. Если скорость замен варьирует, или когда длина и число сравниваемых последовательностей невелико, то может быть получена неправильная топология дерева.

Для дерева, построенного по методу UPGMA длину ветвей можно рассматривать как время расхождения таксонов, отмеренное «молекулярными часами», идущими с постоянной скоростью. То есть сумма времен по любому «пути» от любого узла к терминальным таксонам слева – направо на рис.6.3 постоянна и не зависит от выбора пути.

Таким образом, если все же гипотеза молекулярных часов справедлива, то, при достаточном объеме данных для анализа, метод UPGMA правильно восстановит дерево.

Если рассматривать небольшой промежуток времени эволюции, то ее скорость конечно не будет постоянна. Действительно, новые мутации возникают непрерывно, но случайно. Одни из них могут элиминироваться отбором, другие фиксироваться в популяции случайно (из-за дрейфа генов) или под давлением позитивного отбора.

Если на дереве можно найти такую точку, что расстояния от нее до всех терминальных таксонов одинаковы, то такое дерево называется ультраметрическим. Ультраметрическое дерево можно «укоренить» в эту точку. Таким образом, если алгоритм предполагает точность «молекулярных часов», то реконструированное им дерево укорененное ультраметрическое. Таковым является дерево рис.6.3. Если же алгоритм не опирается на «молекулярные часы», то полученное дерево всегда не ультраметрическое и, как правило, не укорененное (например, метод NJ, описанный ниже, не предполагает точности «молекулярных часов» и строит дерево без корня).
Метод наименьших квадратов (НК).

Если для разных таксонов скорость замен варьирует, то метод UPGMA использовать нельзя, так как он может привести к неправильной топологии. Для этого случая подходит метод наименьших квадратов, поскольку он допускает разную скорость эволюции для ветвей филогенетического дерева.

Итак, оценивается родство между группой таксонов на основе последовательностей, полученных из этих таксонов. В стандартном методе наименьших квадратов, для каждой из возможных топологий филогенетического дерева, вычисляется остаточная сумма квадратов:

где d_ij – наблюдаемое расстояние между последовательностями i, j, а е_ij – ожидаемое расстояние между таксонам i, j. е_ij являются суммой оценок всех длин ветвей дерева, которые соединяют два таксона.

Алгоритм вычисляет значения R_s для всех вероятных топологий. Наиболее вероятной считается та топология, для которой значение R_s минимально.

Используют также взвешенный метод наименьших квадратов, где

Как стандартный, так и взвешенный методы наименьших квадратов часто приводят к явно неправильной топологии, поскольку при наименьшем R_s некоторые из подобранных чисто математически длин ветвей могут оказаться отрицательными. Поэтому один из возможных путей улучшения этих методов – при вычислении ввести запрет на отрицательную длину ветвей.

Итак, чтобы для каждой топологии вычислить R_s необходимо подобрать длины ветвей и вычислить е_ij – их сумму для каждой пары таксонов (i, j). Простой метод оценки был предложен Фитчем и Марголишем.

Р
x

2 [b = 2] [b' = 3]

3 [c = (3–4–5)] [c' = (4–5)]

Рис. 6.4. Дерево для трех таксонов примера. Пояснения в тексте.

где x, y и z – неизвестные длины ветвей для таксонов 1, 2 и 3.

Решив эти уравнения, получаем оценки для длин ветвей:

x=(d₁₂+d₁₃ – d₂₃)/2

y=(d₁₂ – d₁₃+d₂₃)/2 (6.4)

z=(–d₁₂+d₁₃+d₂₃)/2

Итак, даны три таксона. Пусть эволюционные расстояния между каждой парой таксонов: d₁₂=5; d₁₃=7; d₂₃=8.

Использую формулу (6.4) получаем оценки длин ветвей x, y и z.

x=(d₁₂+d₁₃ – d₂₃)/2=(5+7–8)/2=2

y=(d₁₂ – d₁₃+d₂₃)/2=(5–7+8)/2=3

z=(–d₁₂+d₁₃+d₂₃)/2=(–5+7+8)/2=5

После того как оценены длины ветвей x, y, z можно определить (е_ij) – ожидаемые расстояния между таксонами i, j. Это расстояние определяется через сумму оценок всех длин ветвей, соединяющих два таксона. Затем можно вычислить значение R_s по формуле (6.2) или (6.3).

При m = 3 R_s = 0. Полученное решение – единственное. Но при m > 3 число расстояний d_ij становится больше числа ветвей дерева без корня (табл. 6.1). Число линейных уравнений, связывающих d_ij и длины ветвей, становится больше числа неизвестных – длин ветвей. Следовательно, при m > 3 необходимо применять специальные алгоритмы и методы оценки длины ветвей.

Рассмотрим такой алгоритм на примере для 5 таксонов. Матрицу попарных расстояний возьмем из примера построения дерева UPGMA (табл.6.2).

Наименьшее расстояние в этой матрице между таксоном 1 и 2 (d₁₂=5). Обозначим таксоны 1 и 2 через а и b соответственно, а 3 остальных – как с (рис. 6.4). Тогда d_ab=5; d_ас=(7+10+13)/3=10; d_bc=(8+12+15)/3=11,67. Теперь, используя уравнения (6.5), можно вычислить значения x, y и z (рис. 6.4):

x=(d_ab+d_ac – d_bc)/2=(5+10–11,67)/2=1,67

y=(d_ab – d_ac+d_bc)/2=(5–10+11,67)/2=3,34

z=(–d_ab+d_ac+d_bc)/2=(–5+10+11,67)/2=8,34

x и y – оценки числа нуклеотидных замен – длины ветвей до таксонов 1 и 2 (а и b); z – расстояние между составным таксоном с и точкой ветвления для кластера 1–2. Теперь объединим таксоны 1 и 2 и обозначим составной таксон (ab). Пересчитаем расстояния между таксоном (аb) и каждым из остальных (входивших в с) и снова выберем два таксона с наименьшим расстоянием. Обозначим их через а' и b' и объединим оставшиеся таксоны в составной таксон с'. Новые значения x', y' и z' вычисляются как раньше (x, y, z).

Мы уже вычисляли расстояние между кластером (1–2) и остальными таксонами при построении дерева UPGMA (7,5;11;14 – табл.6.3). Наименьшее расстояние в новой матрице между кластером (ab) и таксоном 3. Тогда (ab) и таксон 3 будут представлять новые таксоны а' и b', а в таксон с' объединятся таксоны 4 и 5 (рис.6.4.). Теперь d_a'b'=7,5; d_а'с'=(10+12+13+15)/4=12,5; d_b'c'=(11+15)/2=13. Определяем значения x', y' и z':

Теперь можно «развернуть» полное дерево для 5 таксонов. На рисунке 6.5 представлена полученная топология. В скобках указаны буквенные обозначения ветвей (a, b, c, d, e, f ,g).

1

(c)1

(a)1.67


2

(e)2,75

(b)3,34


3

(d)4


(f)4,75


4

(g)8,25


5

Рис. 6.5. Полученное дерево без корня для 5 таксонов модельного примера.

Уже известно, что (a+b)/2=2,5 и z'=9. Тогда d=13–9=4; с=7,5–2,5–4=1.

Объединим таксоны a'b' Наименьшее расстояние в новой матрице между (a'b') и таксоном 4. Тогда (a'b') и таксон 4 будут представлять новые таксоны а'' и b'', а в таксон с'' попадает таксон 5. Теперь d_a''b''=11; d_а''с''=(14+15)/2=14,5; d_b''c''=13. Определяем значения x'', y'' и z'':

x''=(d_a''b''+d_a''c'' – d_b''c'')/2=(11+14,5–13)/2=6,25

y''=(d_a''b'' – d_a''c''+d_b''c'')/2=(11–14,5+13)/2=4,75

z''=(–d_a''b''+d_a''c''+d_b''c'')/2=(–11+14,5+13)/2=8,25

Длины ветвей e и f определим, используя формулы:

Откуда f=13–8,25=4,75; е=14,5–8,25–1–2,5=2,75

Длина ветви g = z''

Теперь можно вычислить е_ij – ожидаемые расстояния между таксонами i, j, используя рис 6.5.

Далее можно определить R_s по формуле (6.2) и (6.3).

1,0936

Для метода взвешенных наименьших квадратов (формула 6.3):

Следует обратить внимание на несходство полученных оценок длин ветвей дерева, построенного методами UPGMA (рис. 6.3.) и НК (рис. 6.5.).

Чтобы найти дерево методом наименьших квадратов, нужно рассмотреть все возможные топологии. Однако на практике их число огромно и для вычисления R_s рассматривается лишь небольшой процент всех возможных топологий. В методе Фитча и Марголиаша первая топология выбирается как было описано выше (например, для 5 таксонов рис. 6.5), а дальнейшие получаются из нее с использованием различных алгоритмов обмена ветвей. Эти алгоритмы будут описаны при рассмотрении метода парсимонии.