Моделирование эволюции гена

1. Изображение дерева, описанного заданной формулой

Мне была выдана скобочная структура дерева, моделирующая судьбу гена, кодирующего белок CUER_ECOLI:
(((((А:23,В:23):20,С:43):7,D:40):20,Е:70):30,F:100);

Для данной скобочной формулы составляем 2 филогенетических дерева. Заметим по формуле, что дерево не ультраметрическое.

Укорененное:

и неукорененное:

2. Описание ветвей дерева как разбиения множества листьев

	A	B	C	D	E	F
alpha	.	.	*	*	*	*
beta	.	.	.	*	*	*
gamma	.	.	.	.	*	*

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

3. Получение искуственных мутантных последовательностей

Для получения "мутантов" воспользуемся программой msbar пакета EMBOSS. Этой программе нужно указать число мутаций, необходимое для "получения" данного потомка. А в дереве число мутаций рассчитано на 100 нуклеотидов. То есть надо пересчитать число мутаций для данного гена. Длина гена, кодирующего белок CUER_ECOLI, составляет 408 нуклеотидов. Соответственно, пересчитываем число мутаций для каждого потомка, в том числе, конечно, и для последовательностей, находящихся в узлах дерева (то есть общих предков наших "конечных" последовательностей. Например, общий предок последовательностей A, B и С назван ABC). x в расчетах - необходимое число мутаций:

F : 100-100
x - 408 x=408

ABCDE : 30-100
x - 408 x=122,4 (122)

E: 70-100
x - 408 x=285,6 (285)

ABCD: 20-100
x - 408 x=81,6 (81)

D : 40-100
x - 408 x=163,2 (163)

ABC: 7-100
x - 408 x=28,56 (28)

C: 43-100
x - 408 x=175,44 (175)

AB: 20-100
x - 408 x=81,6 (81)

A (B): 23-100
x - 408 x=93,84 (93)

Далее приведен текст скрипта для программы msbar, составленный на основе предыдущих расчетов:

msbar CUER_gene.fasta f.fasta -point 4 -count 408 -auto
msbar CUER_gene.fasta abcde.fasta -point 4 -count 122 -auto
msbar abcde.fasta e.fasta -point 4 -count 285 -auto
msbar abcde.fasta abcd.fasta -point 4 -count 81 -auto
msbar abcd.fasta d.fasta -point 4 -count 163 -auto
msbar abcd.fasta abc.fasta -point 4 -count 28 -auto
msbar abc.fasta c.fasta -point 4 -count 175 -auto
msbar abc.fasta ab.fasta -point 4 -count 81 -auto
msbar ab.fasta a.fasta -point 4 -count 93 -auto
msbar ab.fasta b.fasta -point 4 -count 93 -auto

4. Реконструкция дерева алгоритмами UPGMA, Neighbor-joining и максимального правдоподобия

Чтобы реконструировать дерево алгоритмом максимального правдоподобия, используем программу fdnaml :
fdnaml leaves.fasta -ttratio 1 -auto

В файл leaves.fasta помещены последовательности, соответсвующие листьям дерева. Программа построила следующее неукорененное дерево:

 
  +-----b         
  |  
  |     +----------d         
  |     |  
  1-----2     +-------------------------------f         
  |     |  +--4  
  |     +--3  +---------------------e         
  |        |  
  |        +-------c         
  |  
  +----a         

Чтобы реконструировать дерево алгоритмами UPGMA или Neighbor-joining, сначала надо посчитать попарные расстояния между последовательностями программой fdnadist:
fdnadist leaves.fasta -ttratio 1 -auto

Для реконструкции дерева алгоритмом Neighbor-joining файл leaves.fdnadist (выданный программой fdnadist), надо подать на вход программе fneighbor :
fneighbor leaves.fdnadist -auto

Получаем следующее дерево:

  +----b         
  ! 
  !     +----------d         
  1-----3 
  !     ! +--------c         
  !     +-4 
  !       !  +---------------------e         
  !       +--2 
  !          +---------------------------------f         
  ! 
  +-----a         

Из перечисленных программ только алгоритм UPGMA дает укорененное дерево:

                                         +----------a         
                              +----------1 
                              !          +----------b         
               +--------------3 
               !              ! +--------------------c         
  +------------4              +-2 
  !            !                +--------------------d         
--5            ! 
  !            +------------------------------------e         
  ! 
  +-------------------------------------------------f         


  

5. Сравнение все полученных деревьев между собой и с исходным деревом

	a	b	c	d	e	f	исходное дерево	реконструкция алгоритмом максимального правдоподобия	реконструкция алгоритмом Neighbor-joining	реконструкция алгоритмом UPGMA
1	.	.	*	*	*	*	+	+	+	+
2	.	.	.	*	*	*	+	-	-	-
3	.	.	.	.	*	*	+	+	+	+
4	.	.	*	.	*	*	-	+	+	-
5	*	*	.	.	*	*	-	-	-	+

номерами отмечены ветви (учтены все, существующие во всех 4 деревьях). Следует отметить, что ветвь ?1 присутствует во всех реконструкциях, а также в реальном дереве. Это, вероятно, связано с тем, что длина ветви, отделяющей общего предка последовательностей А и В, относительно велика. Как следствие этого, все программы явно "видят" эту ветвь.
Интересно также то, что реконструкция алгоритмом максимального правдоподобия совпала с реконструкцией алгоритмом Neighbor-joining. При этом они отличаются от реального дерева присутсвием ветви 4. Это также можно объяснить большой длиной ветви, одним из концов которой является общий предок A, B, C и D. Странно то, что при этом последовательность С оказывается отдельно от четырех вышеперечисленных. Возможно, это происходит вследствие большой длины ветви, отделяющей С от всех остальных. С точки зрения теории эволюции это означает, что данная последовательность наиболее "далека" от остальных мутантов (наибольшее количество мутаций должно было произойти для ее появления). При этом в реконструкциях алгоритмом максимального правдоподобия и Neighbor-joining отсутствует ветвь ?2. Ее нормированная длина в скобочной формуле - всего 7. Возможно, она очень мала, чтобы быть "замеченной" и учтенной данными алгоритмами. Зато множество [ABCD] выделяется хорошо. Но внутри него все как раз не так определенно. Это может быть следствием небольшой длины ветви, отделяющей D от [ABC].
Алгоритм UPGMA основан на "молекулярных часах" - количество мутаций прямо пропорционально времени эволюции. Этот алгоритм выделил ветвь, отделяющую мутантов C и D от остальных. Могу предположить, что это произошло, потому что расстояние между предками этих последовательностей небольшое, но при этом количество мутаций, необходимое для "возникновения" данных последовательностей велико, что "относит" мутантов C и D далеко от остальных. Также в качетсве объяснения надо учесть, что с помощью псевдографики невозможно нарисовать длин меньше одной палочки.
Все восстановленнные с помощью алгоритмов деревья отмечают, что мутанты E и F (особенно F) наиболее далеки от остальных.

Bootstrap и drawtree

6. Бутстреп-анализ

Полученные 100 выравниваний подаем на вход программе fdnaml . В выходном файле leaves.treefile получаем 100 скобочных формул, соответствующих реконструкциям, сделанным по каждому из выравниваний.
Запустив программу fconsense, получим результаты бутстреп-анализа. Ниже приведено полученное консенсусное дерево:

  +---------------------------c
  |
  |                    +------a
  |             +-89.0-|
  |             |      +------b
  |      +-34.0-|
  |      |      |      +------f
  +------|      +-55.0-|
         |             +------e
         |
         +--------------------d

A	B	C	D	E	F	реальное дерево	консенсусное дерево	бутстреп-значения
.	.	*	*	*	*	+	+	89
.	.	.	*	*	*	+	-	-
.	.	.	.	*	*	+	+	55
.	.	*	*	.	.	-	+	34

Из таблицы видно, что 2 из 3 ветвей дерева найдены верно. Программа fconsense "отдаляет" мутант С от остальных. В то время как с точки зрения эволюции (реального дерева) мутант F отделился от исходной предковой последовательности раньше, чем С. Но, следует отметить, что бутстреп-значение ветви, выделяющей множество [ABEF] мало (всего 34). Значит, и ветвь эта статистически ненадежна.

7. Изображение исходного дерева