На главную страницу четвертого семестра.

Занятие 3. Дополнительные задания.


Использование статистического способа оценки jackknife, вместо bootstgrap.


В строке Unix вводилась следующая команда, для построения дерева методом jackknife.:

fseqboot argbz.fasta -test j


Метод jackknife - несколько отличный способ размножения выборки, чем метод bootstrap. Метод "складного ножа" (jackknife), предложен М.Кенуем еще в 1949 г. "Размножение выборок" при этом осуществляется путем исключения одного наблюдения. При этом для выборки объема n получаем n "похожих" на нее выборок объема (n - 1) каждая. Если же исключать по 2 наблюдения, то число "похожих" выборок возрастает до n * (n - 1) / 2 объема (n - 2) каждая. Применительно к биоинформатическим способам анализа филогении последовательностей ДНК, метод модифицируется следующим образом: подаваемое на вход программе fseqboot выравнивание последовательностей ДНК, разбивается случайным образом на несколько частей (если не ошибаюсь, в варианте данной программы пакета PHYLIP, общая выборка (выравнивание) разбивается на три части), которые затем случайным образом совмещаются вместе, давая новое выравнивание, меньшего размера. Всего создается 100 таких выравниваний. Затем, для таких кусочных выравниваний с помощью программы fdnaml строятся соответствующие филогенетические деревья, с указанием для каждого из них:
Данная команда вводилась в строке Unix со следующим синтаксисом:

fdnaml jackknife.fseqboot -ttratio 1 -auto


Результаты программы представленны в двух файлах: jackknife.fdnaml, где приведены вышеперечисленные расчеты, и jackknife.treefile, где приведены скобочные формулы построенных деревьев. Теперь требуется из всего набора возможных направлений эволюций, выбрать одно, описываемое консенсусным филогенетическим деревом. Для его построения использовалась программа fconsense, подробно описанная на странице прошлого занятия. Программа запускалась с теми же установками, что и в обязательном задании, чтобы можно было сравнить результаты двух методов: jackknife и bootstrap.

fconsense jackknife.treefile -auto


В результате было получено консенсусное дерево, приводимое ниже. Затем для полученного выравнивания строилось филогенетическое дерево, для которого также анализировалась его топология. В целом, создавалось 100 таких выборок, соответственно 100 филогенетических деревьев, помещенных в файл с расширением .fseqboot.
Рис.1. изображение консенсусного дерева, созданного методом jackknife 
 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
(числа, расположенные на ветвях, обозначают количество 
деревьев, в которых данная ветвь встретилась, в расчете
 на 100 филогенетических деревьев (размер выборки).

                       +------mD.
         +-------100.0-|
         |             +------mC.
  +------|
  |      |             +------mB.
  |      |      +-98.0-|
  |      +100.0-|      +------mA.
  |             |
  |             +-------------mF.
  |
  +---------------------------mE.


Полученное дерево легко можно сравнить с результатами метода bootstrap. Действительно, для запуска программы fseqboot использовались одинаковые параметры, как при использовании метода статистической обработки выборки (в нашем случае общего выравнивания) bootstrap, так и jackknife.
Рис.2. Сравнение консенсусных деревьев, полученных методами (слева направо): jackknife, bootstrap. 
                       +------mD.
         +-------100.0-|
         |             +------mC.
  +------|
  |      |             +------mB.
  |      |      +-98.0-|
  |      +100.0-|      +------mA.
  |             |
  |             +-------------mF.
  |
  +---------------------------mE.


                              +------mD.
                       +-96.0-|
                +-95.0-|      +------mC.
                |      |
         +-66.0-|      +-------------mE.
         |      |
  +------|      +--------------------mF.
  |      |
  |      +---------------------------mA.
  |
  +----------------------------------mB.


Также примечательно сравнить количество включенных и не включенных в консенсусное дерево ветвей:
Рис.3. Сравнение количества включенных и не включенных ветвей, рассчитанных методами (слева направо): jackknife, bootstrap. 
Порядок следования листьев в формулах
разбиения множества ветвей.
 1. mE.
 2. mD.
 3. mC.
 4. mB.
 5. mA.
 6. mF.

Ветви, включенные в дерево:

Тип ветви     частота встречаемости ветви в 
                    выборке деревьев
...***                     100.00
.**...                     100.00
...**.                     98.00


Ветви, не включенные в дерево:

Тип ветви     частота встречаемости ветви в 
                    выборке деревьев
...*.*                      2.00

Порядок следования листьев в формулах
разбиения множества ветвей. (Изначально 
он был иной, но для удобства сравнения
изменен)
1. mE.
2. mD.
3. mC.
4. mB.
5. mA.
6. mF.
 
Ветви, включенные в дерево:
Тип ветви     частота встречаемости ветви в 
                    выборке деревьев
...***                     96.00
**....                     95.00
***..*                     66.00


Ветви, не включенные в дерево:

Тип ветви     частота встречаемости ветви в 
                    выборке деревьев
***.*.                     32.00
*.*..*                      5.00
..*..*                      3.00
.*...*                      2.00
*....*                      1.00


Как видно из приведенных выше рисунков, деревья получились полностью одинаковыми по топологии, и совпадающими с деревом реальной эволюции (см. выполнение обязательных заданий). Но показатели встречаемости расчетных ветвей в построенных деревьях, и также количество невключенных ветвей, сильно различаются. Поэтому, чтобы действительно как-то оценить эффективность данных методов при создании консенсусных деревьев, предлагаю поступить следующим образом: если метод правильно идентифицировал какую-либо ветвь, то число, соответствующее частоте встречаемости данной ветви среди набора деревьев, берем со знаком "+", если же метод идентифицировал ветвь, не наблюдаемой в реальном филогенетическом дереве, и приписал ей какое-либо значение частоты встречаемости, то возьмем это число в общую сумму со знаком "-". Полученное значение суммы частот встречаемости каждой ветви консенсусного дерева отнормируем на n*100%, где n - количество ветвей, присутствующих в реальном дереве. Так можно поступить, так как если сложить все проценты встречаемости по каждой ветви (и включенной и не включенной в консенсусное дерево), то получим как раз значение n*100. Итак, общая формула расчетов выглядит так:

eff = Σ xi / [n*100] * 100%


Итак, расчеты показали следующее:

метод оценки jackknife bootstrap
эффективность 98.7% 42.7%


То есть, фактически метод jackknife эффективнее в 2,3 раза!! Конечно, это очень приближенный и грубый способ оценки эффективности, так как в нем не учитывается то, что в обеих программах установлен порог для значений частот встречаемости ветвей, равный 50%: если значение ветви выше 50%, то оно включается в консенсусное дерево, если наоборот - то не включается. Но тогда эффективность была бы одинакова для обоих методов. Поэтому решено было оценивать методы по "общей дисперсии" всех возможных топологий деревьев и потенциальных ветвей. Но также есть и более объективное преимущество метода jackknife перед bootstrap. Предположим, что по выборке делаются какие-либо статистические выводы. Также необходимо знать, насколько эти выводы устойчивы. Если есть другие выборки, описывающие то же явление, то можно применить к ним ту же статистическую процедуру и сравнить результаты. То есть, фактически применить метод bootstrap.
(В терминологии биоинформатики, это означало бы общее выравнивание последовательностей случайным образом "нарезать" на лоскуты, рассчитать все необходимые данные по этим "лоскутам", а затем сравнить полученные данные.)
А если таких выборок нет? Тогда можно их построить искусственно. Берется исходную выборку и исключается один элемент. Получается похожая выборка. Затем возвращается этот элемент и исключается другой. Получается вторая похожая выборка. Поступив так со всеми элементами исходной выборки, получается столько выборок, похожих на исходную, каков ее объем. Остается обработать их тем же способом, что и исходную, и изучить устойчивость получаемых выводов - разброс оценок параметров, частоты принятия или отклонения гипотез и т.д. То есть использовать метод jackknife.

Укоренение дерева в среднюю точку с помощью программы fretree.


На вход программе fretree подавалось неукорененное дерево, построенное методом Neighbor-joining. Так как в дереве эволюции последовательности гена argb всего шесть листьев, то параметр при запуске равнялся шести. Итак, команда, задаваемая в строке Unix, выглядит так:

fretree 6 argbz_fneighbor.treefile


После исполнения этой записи, программа вывела на stdout изображение неукорененного дерева из файла argbz_fneighbor.treefile и потребовала ввести определенные однобуквенные команды. Набрав "?", было выяснено, что означает каждая однобуквенная команда:
Как видно из приведенного списка описаний односимвольных программ, для создания укорененного дерева в середину из неукорененного, необходимо было провести следующие операции: после получения запроса на ввод однобуквенных команд, ввести:
  1. M - укоренить в среднюю точку
  2. X - выход из программы
  3. Y - на вопрос о сохранении файла ответить кратко "да" (Y)
  4. R - дерево сохранить укорененным

Итак, был получен файл argbz_fretree.treefile, содержащий скобочную формулу укорененного в середину дерева.

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


Полученный файл со скобочной формулой укорененного в середину филогенетического дерева можно визуализировать с помощью соответствующей программы из пакета PHYLIP, которой является программа fdrawgram. Так, для визуализации неукорененных деревьев следует использовать программу fdrawtree, тогда как для создания изображения укорененного дерева следует запускать программу fdrawgram. Последняя, кроме необходимых условий запуска: имена входного и выходного файлов, также имеет множество других опциональных определителей. Среди них есть определитель plotter, параметрами которого задается тип графического файла, в котором будет сохранено изображение дерева; -previewer, с параметрами чего становится возможным быстрый просмотр результатов; также есть некоторое количество статистических определителей, связанных с особенностями работы программы и также возможности расположить итоговое укорененное дерево в определенном месте и положении листа. Но на мой взгляд, наиболее полезными оакзались следующие определеители:
Итак, для получения изображения дерева в виде фенограммы, необходимо использовать следующую команду строки Unix:

fdrawgram fretree.ps fretree1.ps -style p


На вход программа подавалась скобочная формула укорененного в середину дерева из файла fretree.ps, и полученное изображение сохранялось в файле fretree1.ps. Итак, на рисунке ниже изображена данная фенограмма:

Рис.4. Изображение укорененного в середину дерева в виде фенограммы.

Восстановление методом максимального правдоподобия предковой последовательности сравнение её с реальной предковой последовательностью.


Так как метод наибольшего правдоподобия является действительно одним из лучших способов оценки филогении последовательностей, то безусловно затем требуется попытаться восстановить хотя бы приблизительно предковую последовательность. Для выполнения поставленной задачи использовалась программа dnamlk пакета PHYLIP. В качестве исходных данных, программа использует выравнивание последовательностей мутантных генов (в нашем случае), или последовательности генов, для которых исследуется филогения и требуется получить предковую последовательность, или просто восстановить филогенетическое дерево. В отличие от другой программы пакета PHYLIP, fdnaml, в опциях fdnamlk также можно задать другое филогенетическое дерево, полученное иными методами, а не реконструировать его заново. Среди других установок (опциональных, то есть тех, которые можно менять по усмотрению) есть как относящиеся к статистической обработке данных, так и к виду результатов, выдаваемых программой в выходном файле. В числе первых установок, на мой взгляд, интересно изменять следующие определители:
Итак, для начала использовалась следующая команда, вводимая в строке Unix:

fdnamlk argbz.fasta fretree.ps argbz.fdnamlk -ttratio 1 -hypstate Y


После её исполнения в файле argbz.fdnamlk содержалось введенное дерево из файла fretree.ps (расширение здесь представлено .ps, но на самом деле в нем содержалось скобочная формула в формате PHYLIP); для ветвей были пересчитанны значения длин согласно алгоритму ML; отношение частот транзиций и трансверсий принималась равной единице, как и в предыдущем занятии по реконструкции деревьев; и также были реконструированы последовательности предшественников, находящихся в узлах дерева, и также "корневая" последовательность. Из этого файла затем извлекалась предковая (корневая) последовательность и записывалась в формате .FASTA вместе с последовательностью гена argB. Собственно, после всех проделанных операций, как раз эта последовательность и должна была "восстановиться", согласно тому, что именно из последовательности argB были произведены все шесть мутантных последовательностей. Затем это, фактически, выравнивание последовательностей в формате .fasta импортировалось в файл программы GENEDOC для визуализации выравнивания. В результате получено следующее изображение (приведена лишь часть выравнивания):

Рис.5. Изображение участка выравнивания предсказанной предковой последовательности и реальной последовательности гена argB


На рисунке завездочкой "*" обозначена третья позиция кодона (на каждый десятый раз приводится значение длины последовательности в цифрах), и сравниваемые последовательности объединены в группу, выделенную красным цветом; в названии предковой последовательности также указано значение ttratio (ttr1), и файл использовавшегося дерева (fretree). Как видно из этого участка, обнаружены некоторые участки сходства, но они довольно разрозненны и если перевести в белковую последовательность, то эти участки "уменьшатся" ещё в три раза. Также (на рисунке не показано) программа fdnamlk при реконструкции предковой последовательности ДНК, часто отмечает позиции, в которых алгоритм не уверен за определенный нуклеотид, следующими символами:
отчего информативность выравнивания снижается. Так же снижается содержательность самой последовательности: фактически, мы получаем не точный порядок нуклеотидов в предковой последовательности, а приблизительную оценку того, как она должна выглядеть. Но стоит отметить, что в целом эта оценка неплохая: даже сайты с обозначениями из списка, часто точно совпадают с природой нуклеотида в реальной последовательности, то есть если в реальной последовательности на n-ой позиции стоит нуклеотид G, то в реконструированной последовательности на месте этого сайта стоит R или др. символ. Таким образом, неопределенность в последовательности снижается от четырех возможных нуклеотидов на данный сайт, до двух нуклеотидов. В целом, такие участки довольно многочисленны, и их можно видеть на рисунке ниже (отмечены голубым):

Рис.6. Изображение участка выравнивания предсказанной предковой последовательности и реальной последовательности гена argB с выделенными сайтами неустановленных нуклеотидов.


Уменьшение информативной содержательности сопровождает рисунок №7, где последовательности генов транслированны в последовательности белков, и приведен именно исследуемый участок. Но я полагаю, этот непривлекательный факт не должен повергать в уныние, так как всегда можно провести поиск реального предка с помощью реконструированной последовательности программой BlastN, или другими вариантами, использующими на входе нуклеотидную последовательность, и работающих с базами данных как нуклеотидов, так и белков. И в принципе, таких разрозненных участков сходства должно хватить для обнаружения возможных гомологов.

Рис.7. Уменьшение информативности выравнивания при транслировании в белковое выравнивание.


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

fdnamlk argbz.fasta fretree.ps argbzm.fdnamlk -ttratio n -hypstate Y


Где n - [1 и 0,5], m - [2 и 3]. То есть было получено два файла: argbz2.fdnamlk, argbz3.fdnamlk, из которых извлекались реконструированные последовательности. Затем с ними проводились те же операции, описанные выше. В итоге, получено выравнивание этих последовательностей с реальной предковой последовательностью, которое затем также импортировалось в msf-формат. Поэтому, имеем следующее изображение того же участка выравнивания, который иллюстрировал результаты программы выше:

Рис.8. Сравнение результатов программы fdnamlk при разных значениях ttratio.


На данной картинке изображены три выравнивания, каждое состоит из реконструированной последовательности и реального предка. Каждое выравнивание собрано в группы, классифицированные по значению ttratio, с помощью которого были построены предположительные предки. Красным выделено выравнивание, где у построенной последовательности ttratio - 1, желтым - ttratio - 0.5, оранжевым - ttratio - 2. Как видно из выравниваний, наилучшим оказалось та группа, где ttratio = 0.5, при реконструкции предковой последовательности. Чтобы показать этот факт более наглядно, дополнительно верно идентифицированные нуклеотидные остатки покрашены красным. Но здесь приведен лишь участок выравнивания, на самом деле дополнительно определенных остатков ещё больше. Стоит отметить, что это довольно сильно сказывается на информативности восстановленной последовательности, так как "заполняются" бреши в числе недостаяющих остатков в кодонах, кодирующих аминокислотную последовательность белка. Возможно, это связано с тем, что ttratio устанавливает именно отношение частот транзиций и трансверсий, а не отношение вероятностей. Чем чаще происходят замены между практически одинаковыми по физико-химическим свойствам нуклеотидными основаниями, тем "точнее" программа идентифицирует остатки в предковой последовательности, тем больше в ней будет символов A, G, T, C, R, Y. Что впрочем и наблюдается на практике.
Для подтверждения увеличения информативности при использовании ttratio = 0.5, можно оттранслировать последовательности в выравнивании в белковую последовательность и сравнить количество идентичных остатков в последовательности реального предка и реконструированного:

Рис.9. Сравнение результатов программы fdnamlk при разных значениях ttratio.


Также можно видеть, что при ttratio = 2 вообще ни одного остатка не было идентифицировано точным и единственным образом.
Вообще, был создан файл root.msf, в котором содержатся данные про всё вышесказанное, но только для целой последовательности.
И, наконец, было интересно проверить предположение об идеальных частотах нуклеотидных оснований в уникальных участках ДНК - кодирующих участках. В строке Unix вводилась следующая команда:

fdnamlk argbz.fasta argbz_fdnaml.treefile argbz8.fdnamlk -freqsfrom N -basefreq 0.25_0.25_0.25_0.25 -ttratio 0.5 -hypstate Y


С помощью этой команды ttratio задается равным 0.5, а частоты каждого нуклеотида в выборке установлены на идеальные значения: 0.25 для каждого нуклеотида. Итак, проделывая с полученным файлом argbz8.fdnamlk ту же последовательность действий, что делалось для каждого такого файла выше, удалось "вытащить реконструированную последовательность предка и импортировать её в .Fasta формате в файл root.msf. Но в данном случае интереснее сравнить выравнивание этой последовательности (обозначена на рисунке и в файле как 1ttr0.5_freq0.25) и гена argB с выравниванием, которое признано лучшим в предыдущем рисунке (группа последовательностей, выделенная желтым цветом). Так в итоге получается, что если использовать идеальные частоты нуклеотидов в последовательности, тем лучше получится результат - намного больше было идентифицировано одинаковых нуклеотидов с реальным предком (см. Рис.10).

Рис.10. Сравнение результатов программы fdnamlk при разных значениях частот нуклеотидов: в желтой группе, программа, при реконструкции предка, использовала свои, расчетные частоты нуклеотидов, в фиолетовой группе программе задавались идеальные частоты нуклеотидов.


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




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