Занятие 2. Моделирование эволюции гена.

(((А:60,В:65):10,((С:25,D:25):15,Е:40):30):20,F:120);

Задание №1. Создание изображения дерева, описанного заданной скобочной формулой.

Рис.1. изображение дерева, восстанавливаемое по скобочной формуле.

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

     A B C D E F     
     * * . . * *
     * * . . . *
     * * . . . . 

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

• Создана поддиректория "Mutant", куда занесены последовательности мутантов, в рабочей директории Practice2.
• Длина моего гена: 777 п.н. Поэтому необходимо пересчитать число нуклеотидных замен на всю длину гена. Но сперва надо локализовать корень на бескорневом дереве (см. Рис.1.). Предложено расположить его в ветви F <->BAECD, так как здесь наиболее большие эволюционные расстояния между листом F и предшественником листьев BAECD. Также Итак, укорененное дерево выглядит так:
  
  

  Рис.3. изображение укорененного дерева

• Соответственно, пересчет эволюционных расстояний как число нуклеотдных замен:
  Для этого надо знать длину гена argB: 777 п.н. Тогда число замен можно рассчитать по общей формуле:

  N = 777/100*n = 7,77*n

  
  где n - число замен на 100 нуклеотидов.
• В результате были рассчитаны все эволюционные расстояния, выраженные в количествах нуклеотидных замен в последовательности гена argB:
  
  

  Рис.4. изображение дерева с числом замен в последовательности гена argB. Последовательности промежуточных мутантов в узлах не подписаны.

• Был написан скрипт, генерирующий последовательности всех мутантов, представленных на рисунках с помощью программы msbar. Последовательности промежуточных мутантов, соответствующих узлам на филогенетическом дереве, обозначены mXYZ, где XYZ - названия мутантов, соответствующих листьям и которые происходят из данного мутанта. Он выглядит следующим образом:
  
  

  Рис.5. Скрипт-файл m2.script (сохранен в Unix-формате).

  msbar argb.fasta mABCDEF.fasta -point 4 -count 155 -auto msbar argb.fasta mF.fasta -point 4 -count 932 -auto msbar mABCDEF.fasta mAB.fasta -point 4 -count 78 -auto msbar mABCDEF.fasta mCDE.fasta -point 4 -count 233 -auto msbar mAB.fasta mA.fasta -point 4 -count 466 -auto msbar mAB.fasta mB.fasta -point 4 -count 505 -auto msbar mCDE.fasta mCD.fasta -point 4 -count 117 -auto msbar mCDE.fasta mE.fasta -point 4 -count 311 -auto msbar mCD.fasta mC.fasta -point 4 -count 194 -auto msbar mCD.fasta mD.fasta -point 4 -count 194 -auto cat mF.fasta >> argbz.fasta cat mA.fasta >> argbz.fasta cat mB.fasta >> argbz.fasta cat mC.fasta >> argbz.fasta cat mD.fasta >> argbz.fasta cat mE.fasta >> argbz.fasta

  
  
• В результате исполнения данного скрипта были получены 10 мутантных последовательностей:mABCDEF.fasta, mF.fasta, mAB.fasta, mCDE.fasta, mA.fasta, mB.fasta, mCD.fasta, mC.fasta, mD.fasta, mE.fasta и один файл argbz.fasta, содержащий последовательности мутантов гена argB, соответствующие листьям на филогенетическом дереве.

Реконструкция дерева эволюции гена argB алгоритмами UPGMA, Neighbor-joining и Maximum Likelihood (максимального правдоподобия).

1. Сперва был воссоздано дерево эволюции методом максимального правдоподобия. Это достигалось посредством использования команды

   fdnaml argbz.fasta -ttratio 1 -auto

   
   Значение ttratio устанавливает отношение частот транзиций и трансверсий. С математической точки зрения, частота транзиций в два раза меньше частоты трансверсий - так как возможных комбинаций замен в два раза меньше (пиримидин <-> пиримидин, пурин <-> пурин, в отличие от пиримидин <-> пурин). Правда это верно в том случае, если все замены равновероятны. Таким образом, если задать значение ttratio = 2, то можно учесть различную частоту трансверсий и транзиций - фактически восстанавливая алгоритм двупараметрической модели Кимуры. И программа fdnaml, в качастве установок по умолчанию, предлагает работать именно с значением ttratio = 2. Но в задании предложено использовать значение ttratio = 1, то есть не разграничивать типы замен, а приписать всем им одинаковую частоту. Но тогда мы имеем дело с алгоритмом однопараметрической модели Джукса-Кантора.
   В результате были получены два файла: argbz.fdnaml и argbz.treefile. argbz.treefile содержит скобочную формулу реконструированного дерева,а в файле argbz.fdnaml представлено довольно много интересных данных: эмпирические частоты встречаемости нуклеотидных оснований в мутантных последовательностях, визуализация дерева, данные о длинах ветвей, и также подсчеты границ доверительных интервалов для значений длин ветвей (подробнее см. ниже). Итак, дерево выглядит следующим образом:
   
   

   Рис.6. изображение неукорененного дерева, созданного методом максимального правдоподобия

   +--------mE. +------4 | | +-----mD. | +--3 | +-----mC. | | +--------------mA. 2---1 | +---------------mB. | +------------------------------------mF.



2. Для реконструкции дерева алгоритмами UPGMA или Neighbor-joining, сначала были подсчитаны попарные расстояния между последовательностями программой fdnadist:

   fdnadist argbz.fasta -ttratio 1 -auto

   
   Программа fdnadish использует множественное выравнивание последовательностей для подсчета попарных расстояний. Но так как работа проводится в предположении, что эволюция происходит путем точечных замен (без учета возможных инверсий, транспозиций, кроссинговера и др.), то выравнивание которое, нужно подавать на вход программе, совпадает с содержимым файла argbz.fasta. Подсчеты программы сохранены в файле argbz.fdnadist. В этой матрице расстояний эволюционные расстояния приведены в количестве замен, приходящемся на 1 н.п (а не на 100 н.п.). В более удобном виде матрица попарных расстояний представлена ниже (расстояния переведены в количество замен на 100 нуклеотидов):
   
   
   
   Вычисленные расстояния в целом соответствуют истинным эволюционным расстояниям, но все же есть некоторые различия. Так наибольшие различия (до 200 замен на 100 н.п.) свойственны последовательностям mutF и mutA, хотя в истинном дереве наиболее удаленными являются мутанты mutB и mutF (теоретически, расстояние между ними равно 65 + 10 + 20 + 120 = 215 нуклеотидных замен), хотя для них рассчитано 177 нуклеотидных замен. Видимо это является следствием инструментальной ошибки алгоритма джукса-Кантора, проявляетмой при достаточно высоких эволюционных расстояний.
   
   Для получения дерева с помощью алгоритма Neighbor-joining, этот файл подавался на вход программе fneighbor, используя следующую команду :

   fneighbor argbz.fdnadist -auto

   
   В результате были получены два файла: файл argbz.fneighbor и файл argbz.treefile. Первый из них содержит информацию о длинах ветвей реконструированного дерева и визуализацию дерева в текстовом формате. Файл argbz.treefile содержит скобочную структуру реконструированного дерева. Итак, дерево выглядит так:
   
   

   Рис.7. изображение неукорененного дерева, созданного методом Neighbor-joining

   +----------------mA. +--3 ! +-------------mB. ! ! +-----mC. ! +---1 4-------2 +-----mD. ! ! ! +--------mE. ! +--------------------------------------mF.

   
   
3. Для получения дерева с помощью алгоритма UPGMA, файл argbz.fdnadist подавался на вход программе fneighbor, используя следующую команду :

   fneighbor argbz.fdnadist -treetype u -auto

   
   
   В результате были получены два файла: файл argbz_upgma.fneighbor и файл argbz_upgma.treefile. Первый из них содержит информацию о длинах ветвей реконструированного дерева и визуализацию дерева в текстовом формате. Причем построено было укорененное дерево. Файл argbz.treefile содержит скобочную структуру реконструированного дерева.
   Итак, укорененное дерево выглядит так:
   
   

   Рис.8. изображение укорененного дерева, созданного методом UPGMA

   +-------------------------------------------------------mF. ! --5 +------------------------------mA. ! +---3 ! ! +------------------------------mB. +--------------------4 ! +----------mC. ! +------1 +----------------2 +----------mD. ! +-----------------mE.

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

   тип ветви: A  B  C  D  E  F	True	Maximum likelihood	Neighbor-joining	UPGMA
   *  *  .  .  *  *	X	X	X	X
   *  *  .  .  .  *	X	X	X	X
   .  .  .  .  .  *	X	X	X	X
   .  .  *  *  *  *	X	X	X	X

   
   
   Как можно видеть из данных таблицы, все встреченные ветви содержатся как в реальном, так и реконструированных деревьях. Если воспринимать неукорененные деревья как множество всех возможных укоренений, то все построенные разными методами неукорененные деревья (Nj, Ml) также полностью совпадают с реальным. Таким образом, топология полученных деревьев одинакова и полностью совпадает с исходным "истинным" филогенетическим деревом. Но тогда для оценки эффективности каждого из алгоритмов придется сравнивать данные вычисленных расстояний. Итак, для возможности сравнивать все изображения имеющихся деревьев, они все приведены в нижеприводимой таблице. Все изображения были получены с помощью программ drawgram (строит укорененные деревья по скобочной формуле), и drawtree, (строит неукорененные деревья). Этими программами довелось пользоваться при выполнении заданий второго семестра, которые доступны на сайте Пастеровского института:
   
   

   Рис.9. Изображения всех деревьев, реконструированных разными методами и также реального, описываемого данной мне скобочной формулой. Сверху вниз, слева направо представлены результаты программ: ML, NJ, UPGMA, True

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

   S = [1/(n - 1)*Σ(X_i - X_p)²]

   
   Где Х_i - значение расчетной длины ветви одной из программ; X_p - теоретическое значение длины той же ветви одного из деревьев; суммирование по количеству ветвей в деревьях (соответственно в неукорененных будет девять ветвей, а в укорененном - десять); тогда n в случае неукорененного дерева равно n = 9; а в случае укорененного n = 10. Конечно, данная формула верна для последовательности случайных событий, с большим размером выборки и "хороших" распределений. Но в принципе, в эту формулу можно вложить несколько иной смысл, нежели воспринимать её как точечная оценка дисперсии случайной величины. Действительно, достаточно представить гипотетическую модель: расположим все рассчитанные ветви каждого дерева (а также истинного) одна за другой. Тогда сумма квадратичных отклонений длин каждой ветви по сравнению с соответствующей "истинной" веткой истинного дерева, нормированная на "(количество ветвей) - 1" и приведенная к одной размерности длины ветви (то есть извлечен корень из полученного числа), покажет на сколько полученное программой значение длин ветвей (фактически, их сумма) отличается от истинного значения (в плюс/минус сторону). Тогда по итогам нехитрых расчетов, получились следующие оценки:
   
   

   Рис.10. Оценки отклонения "размера" расчетного дерева от "истинного".

   программа	NJ	ML	UPGMA
   оценка	9.25	10.52	22.26

   
   
   Как видно из полученных оценок, наименьшим отклонением от длины истинного дерева отличается дерево, построенное по алгоритму Neighbor-joining (NJ). Фактически, примерно на 9-ть мутаций расчетное дерево "отклонилось" от длины истинного дерева. Но не далеки оценки и другого метода: Maximum likelihood (ML): всего в 11 мутаций отличается от истинного дерева. Фактически, это показывает, что предпринятая оцека не позволяет разграничить результаты этих методов, общее для которых является то, что воспроизводят они неукорененное дерево. Но зато эта оценка позволяет наглядно показать несостоятельность укорененного дерева, полученного программой UPGMA. На целых 22 мутации отличается длина истинного дерева от реконструированного. Видимо, это объясняется особенностями алгоритма программы и исходными данными, которые подавались на вход. Действительно, UPGMA восстанавливает укорененное, ультраметрическое дерево. Ультраметрия в основном верна только для случая выполнения гипотезы "молекулярных часов" для эволюционирующих последовательностей. Но этот экспериментальный факт верен не для всяких биологических текстов, а только для достаточно "близких" и довольно консервативных последовательностей. Но уже в скобочной формуле, описывающей реальное дерево, четко видно, что для данных мутантов эта теория не применима!!! - действительно, длины всех ветвей от корня до каждого листа не равны друг другу (120 ≠ 95 ≠ 90). Поэтому, опираясь на применимость теории "молекулярных часов", UPGMA создает неверное дерево. Видимо, это также связано с теми мутантами, подаваемыми на вход программе: количество мутаций в них было довольно большое (если учесть, что эволюционные рассотояния измеряются 10⁵ лет, и за это время действительно происходят видимые изменения в биологическом тексте (максимум 10 мутаций), то в моем случае для некоторых мутантов (mutF) длина в эволюционном времени огромна. Поэтому, приводимые для реконструкции филогении мутанты можно считать довольно далекими, потерявшими какое-либо сходство. Но отличные результаты показали алгоритмы, воспроизводящие неукорененные филогенетические деревья: Maximum Likelihood, Neighbor-joining. Очевидно, что идея поиска пару листьев, для которых сумма длин веток такого дерева минимальна для довольно отдаленных последовательностей, очень плодотворна и показывает высокие результаты.

Задание №2. Bootstrap и drawtree.

1. Программой fseqboot создано 100 бутстреп-реплик выравнивания из файла argbz.fasta:

   fseqboot argbz.fasta -auto

   
   Каждая реплика - это множественное выравнивание в формате PHYLIP (в верхней строке указано количество последовательностей в выравнивании и длина выравнивания, ниже расположено само выравнивание). Поэтому, подав на вход выравнивание в формате fasta, которое не может быть корректно обработано программами пакета PHYLIP, программа fseqboot выдает реплики в другом формате, удобном для анализа программами пакета PHYLIP. Это помогает избежать трудностей по преобразованию формата реплик из одного типа файлов в другой. Помимо формата fasta программа fseqboot способна обрабатывать исходные выравнивания и в других форматах, например, в формате ClustalW. Таким образом, возможность анализа выравниваний последовательностей в различных форматах является одним из преимуществ версий программ пакета PHYLIP.
2. Полученные 100 выравниваний поданы на вход программе fdnaml. В результате программа произвела 100 скобочных формул, представленных в файле argbz_100.treefile, соответствующих реконструкциям, сделанным по каждому из выравниваний. ttratio положено то, которое предлагается программой по умолчанию (то есть отношение частот транзиций и трансверсий равно двум).
3. Затем, была запущена программа fconsense. Синтаксис данной программы в командной строке Unix таков:

   fconsense argb_100.treefile -auto

   
   В качестве входящего файла подан результат предыдущей программы (argbz_100.treefile). Программа fconsense реконструирует консенсусное дерево с помощью одного из четырех возможных методов: (ml) minimum fraction (от 0.5 до 1.0), (mj) majority rule, (mre) extended majority rule и (s) strict consensus tree. Используемый метод задается параметром программы ,b>method.Первый метод, minimum fraction,включает в состав консенсусного дерева только те ветви, которые присутствуют более чем в 50% количестве реплик.То есть те ветви, которые обнаруживаются в более чем половины деревьев, реконструированных по репликам,включаются в консенсусное дерево, а ветви, встречаемыме меньше, чем в половины деревьев, игнорируются. Методы majority rule и strict являются частными случаями метода minimum fraction: для первого из них значение параметра mlfrac принято равным 0.5, для второго — 1.0. Метод extended majority rule, используемый программой fconsense по умолчанию, сначала включает в состав консенсусного дерева те ветви, которые присутствуют более чем в половине реконструированных по репликам деревьев (как и метод majority rule), а затем в порядке уменьшения процента встречаемости добавляет оставшиеся ветви, которые совместимы с уже включенными. Под совместимостью ветвей понимается то, сто не происходит пересечений новых ветвей с более достоверными, включенными ранее. При этом возможно получение нерахрешенного дерева, которое понимается так: "лишняя ветвь" (то есть та, что является третьей при любом ходе эволюции, при условии, что в один момент не могло произойти образования трех видов) может занять лябое из доступных ветвей, образующих этот узел. (По данным использования параметра -help к программе fconsense.
   В результате получено два файла: argbz.fconsense и argbz_consense.treefile (он был переименован, чтобы не перезаписывать имеющиеся данные в одноименном файле).Первый файл содержит список всех ветвей (в виде разбиения листьев), присутствующих в деревьях; для каждой из них указано количество реплик [в процентах от общего объема реплик], и соответствующих деревьев, в составе которых были выявлены данные ветви; список ветвей,вошедших в состав консенсусного дерева. Приведена визуализация консенсусного дерева - дерево неукорененное, так как реконструкция деревьев по репликам проводилась с помощью метода maximum likelihood. Второй файл (argbz_100.treefile) содержит скобочную формулу консенсусного дерева, в которой вместо длин ветвей указано количество деревьев в составе которых эти ветви были выявлены (длины ветвей, отделяющих один лист от всех остальных, составляют 100%, так как такие ветви присутствуют в любом дереве, реконструированном для данных последовательностей).
   
   

   Рис.9. изображение консенсусного неукорененного дерева, созданного на основе 100 бутстреп-реплик

   Extended majority rule consensus tree (метод построения (по умолчанию) CONSENSUS TREE: the numbers on the branches indicate the number of times the partition of the species into the two sets which are separated by that branch occurred among the trees, out of 100.00 trees +------mD. +-96.0-| +-95.0-| +------mC. | | +-66.0-| +-------------mE. | | +------| +--------------------mF. | | | +---------------------------mA. | +----------------------------------mB.

   
   
   По топологии реконструированное дерево полностью совпало с реальным деревом. Можно отметить интересную закономерность: если эволюционное расстояние между двумя мутантами (например, входящих в состав узлов, между внутренними ветвями) небольшое (количество мутаций, отделяющих их, находится в пределах 50-ти), то процент встречаемости данной ветви в полученных бутстреп-репликах повышается. Так длина ветви между мутантами mutCD - mutCDE равна 15 нуклеотидных замен, то она встречается в 96% деревьях, реконструированных по соответсвующим бутстреп-репликам, тогда как длина ветви mutCDE - mutCDEF равна 30-ти, то она встречается уже в 95% деревьев.
   
   

   Рис.10. Выдача файла argbz.fconsense

   Species in order (мутантные последвательности, соответствующие листьям): 1. mB. 2. mA. 3. mE. 4. mD. 5. mC. 6. mF. Sets included in the consensus tree (ветви, включенные в консенсусное дерево): Set (species in order) How many times out of 100.00 (частота встречаемости) ...**. 96.00 ..***. 95.00 ..**** 66.00 Sets NOT included in consensus tree (ветви, не включенные в консенсусное дерево): Set (species in order) How many times out of 100.00 .****. 32.00 ...*** 5.00 ....** 3.00 .*...* 2.00 ...*.* 1.00

   
   
   
   
   

   Рис.10. изображение дерева и скобочной формулы (повторное).

   
   
   Или можно сравнить множества ветвей деревьев:
   
   

   тип ветви: A  B  C  D  E  F	реальное дерево	реконструированное дерево
   *  *  .  .  *  *	X	X
   *  *  .  .  .  *	X	X
   .  .  .  .  .  *	X
   .  .  *  *  *  *	X	X

   
   
   Как видно, реконструированное дерево содержит все истинные ветви и полностью совпадает с истинным деревом.
   
   

   Задание №3. Изображение реального дерева программой fdrawtree.

   
   Для реализации программы, в строке Unix вводилась следующая команда:

   fdrawtree tree.txt tree.ps -auto

   
   С помощью программы fdrawtree получено следующий вариант визуализации истинного дерева:
   
   

   Рис.11. Изображение дерева программой fdrawtree.

   
   Программа строит неукорененное дерево, в формате postscript (по умолчанию); длины ветвей дерева пропорциональны эволюционным расстояниям.
   
   
   
   

   ---

   ©Володя Рудько