Monthly Archives: January 2013

Реконструкция филогенетического дерева в BEAST/BEAUTi с использованием разных геномных участков и вторичной структуры генов

Рассмотрим пример построения филогенетического дерева, когда в качестве исходных данных используется выровненная нуклеотидная или аминокислотная матрица, состоящая из более чем одного гена.

Исходным файлом для такой задачи является файл с данными в формате nexus, где помимо блока данных и их описания (BEGIN DATA;) присутствует блок разбивки матрицы на блоки (BEGIN MRBAYES;). Разбивка нуклеотидной матрицы на блоки может осуществляться по целом ряду признаков – например, по типам ДНК (митохондриальная, ядерная, рибосомная), по функциям ДНК, по их вторичной структуре (интроны и экзоны), по генами.  Структура такого файла будет выглядеть следующим образом:

#NEXUS

BEGIN DATA;

DIMENSIONS NTAX=3 NCHAR=14;

FORMAT DATATYPE=DNA GAP=- MISSING=? ;

MATRIX

‘Ctenotus angusticeps’  TAAAACATTGGGTT [14 bp]

‘Ctenotus astarte’            ACCCTTCACTTTTAA[14 bp]

‘Ctenotus brooksi’           AATGATTGAATCTG[14 bp]

;

END;

BEGIN MRBAYES;

CHARSET 12S = 1-7;

CHARSET 16S = 8-14;

partition favored = 2: 12S, 16S;

set partition = favored;

END;

Обратите внимание, как именно в блоке BEGIN MRBAYES задается разбивка данных.

Во-первых, мы прописываем начальную и конечную позиции для каждого гена в нуклеотидной матрице. Например, в примере выше, у  нас два кодирующих гена (12S и 16S). Ген 12S начинается с 1го нуклеотида в матрице и заканчивается 7ой позицей, ген 16S начинается с 8ой позиции и заканчивается на 14ой. Каждый ген имеет свое имя и задается командой CHARSET. Во-вторых, мы в обязательно порядке должны указать разбивку с помощью команды partition. В этой команде мы указывем имя разбивки (favored), сколько блоков будет в данной разбивке, и их имена. В-третьих, с помощью команды set partition мы прописываем ту разбивку, которая будет активной в данном nexus файле. Таким образом, очевидно, что в блоке BEGIN MRBAYES может быть прописано несколько разных разбивок (команда partition), но быть активной должна быть только одна (команда set partition).

Про то, как именно можно получить такой файл читайте в заметке про SequenceMatrix. В сегодняшнем упражнении, мы воспользуемся набором данных, полученным в предыдущих упражнениях (скачать файл для обработки;логин – public@phylogenetics.ru; пароль: cornus).

Запустим программу BEAUTi и загрузим в нее файл ctenotus_4genes.nexus.

pr23_1

Так как мы использовали файл с указанной разбивкой – а в нашем случае этой разбивкой являеются 4 гена (12S, 16S, ATP, NADH) – BEAUTi автоматически разбиала набор данных на четыре строки, по одной для каждого блока данных. Однако, до тех пор, пока параметры этих четырех блоков «залинкованы» (linked) информация о структуре матрицы не будет использована при реконструкции филогенетического дерева.

У нас есть три опции разлинковки блоков:

  • Разлинковка эволюционных моделей нуклеотидных замен (для каждого гена – свои параметры и своя эволюционная модель)
  • Разлинковка модели эволюционных часов (при датировке дерева)
  • Разлинковка топологий деревьев

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

pr23_2

Чаще всего, используется разлинковка эволюционных моделей. Для того, чтобы разлинковать наши блоки, выделите правой кнопкой мыши все четыре строки с данными. Кнопки разлинковки станут активными. Нажмите кнопку Unlink Subst. Models. Посмотрите как поменялась колонка Subst.Models.

pr23_3

Теперь в колонке Subst.Models каждая строка с данными имеет свое уникальное название. Таким образом, параметры будут оцениваться отедельно для  гена 12S, гена 16S, гена ATP, и гена NADH.

pr23_5

Оставим остальные две колонки нетронутыми. При этом, неважно по имени какого из генов (блоков) будут называться параметры в колонках Clock Model и Partition Tree – важно то, что везде стоит одинаковое название, а значит параметры будут общими для всех генов (блоков матрицы).

Перейдем во вкладку Taxon Sets. Нажмите на плюс «+» расположенный в левой нижней части окна. В правом подокне появится запись Untitled1 (в колонке Taxon Sets), а в левом подокне в разделе Excluded Taxa появится список всех видов из файла nexus. Переимнуйте untitled1 в ctenotus_prasino и поставьте галочку в колонке Monophyletic?. Выделите все виды в колонке Excluded Taxa кроме Sphenomorphus_muelleri. После того, как вы выберете виды в колонке, небольшая стрелочка справа от колонки станет активной. Нажмите на эту стрелочку. Все выделенные виды перемещены в колонку Included Taxa. В колонке Excluded Taxa остался только один вид, который в нашем случае не относится к роду Ctenotus. В результате этих действий мы сможем задавать датировку дереву на узле между двумя родами Sphenomorphus и Ctenotus (~75 миллионов лет).

pr23_6

Пропустите закладку Tip Dates. В данном упражнении у нас нет информации для этой закладки, поэтому оставьте все параметры неизменными. Точно также ничего не меняйте в закладке Traits. Она может быть использована в тех случаях, когда есть дополнительная информация о морфологии вида, но ее использование должно быть строго оправдано, так как морфологические признаки (особенно конвергентные) могу нарушать процесс реконструкции филогенетического дерева.

Перейдите во вкладку Site Models. Обратите внимание на то, что в списке справа в подокне Substitution Models у нас есть 4 гена (четыре блока матрицы). Для какждого гена мы можем выбрать свою эволюционную модель. Из списка Substitution Model выберите модель GTR для первого гена, из списка base frequencies выберите Estimated, а из списка Site Heterogeneity Model выберите Gamma и Number of Gamma Categories оставьте равной 4.

pr23_8

Последовательно выделяя одну за другой строчки в правом подокне, повторите шаги выше для трех остальных генов.

pr23_9

Несмотря на то, что мы выбрали идентичную эволюционную модель для всех генов (GTR+GAMMA) параметры будут оцениваться независимо для каждого гена. Чтобы выяснить какая эволюционная модель оптимальная для каждого гена прочитайте статью Тестирование эволюционных моделей.

В закладке Clock Model в колонке Model выберите Relaxed Clock: uncorrelated lognormal и в колонке Estimate поставьте галочку (таким образом мы позволяем модели варьировать скорость нуклеотидных замен в разных частях дерева)

pr23_10

Перейдите во вкладку Trees. Оставьте Link tree prior for all trees включенным, а в разделе Tree Prior выберите из списка Speciation: Yule Process

pr23_11

Перейдите во вкладку Priors. В колонке Parameter выберите tmrc(ctenotus_prasino) – то есть априорное распределение для узла между группой Сtenotus и Sphenomorphus. Нажмите в колонке Prior на Using tree prior. В появившемся окне выберите из списка априорное распределение Normal, укажите mean 78 и stdev 2.8. Нажмите OK.

pr23_12

Таким образом мы задаем возраст для узла соединяющего две клады – Сtenotus и Sphenomorphus. Этот узел позволяет задать референсную точку всему филогенетического дереву и произвести его датировку. Одного узла для точной датировки дерева почти всегде не достаточно, особенно если этот узел расположен ближе к вершине дерева. Наиболее оптимальным является наличие трех-четырех дат расположенных равномерно по всему дереву (например, один узел ближе к корню дерева, один узел в центре, и два – ближе к вершине дерева). Источником информации для датирования деревьев являются палеонтологические данные. В данном случае мы воспользовались информации полученной в результате датировки более крупного филогенетического дерева для нескольких таксономических групп ящериц из научной статьи Rabosky 2007. Этот метод не идеален, но может быть использован для примерной датировки.

Теперь выберите строку ulcd.mean нажмите в колонке Priors на ?Not yet specified и в открывшемся окне и укажите Uniform, Initial Value = 0.00001, Lower =0, Upper = 0. Нажмите OK.

pr23_13

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

Перейдите в последнюю вкладку MCMC. Укажите следующие параметры: Length of chain 100 000 000, Echo state to screen every = 10000, Log parameters = 10000. Поставьте галочку напротив Create tree log file with branch length in substitutions.

pr23_14

Нажмите кнопку Generate BEAST file и укажите названия xml-файлу ctenotus_dated.xml. Файл сгенерирован.

Теперь запустите программу BEAST. Нажмите на кнопку Choose file… и выберите созданный нами xml-файл. Нажмите RUN. В зависимости от мощности компьютера расчет MCMC может занять до нескольких суток.

Признак-зависимая диверсификация: практикаCharacter-Associated Diversification: practice

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

Запустим программу R и установим рабочую директорию, в которую нужно сохранить все рабочие файлы (загрузить файлы; логин – public@phylogenetics.ru; пароль: cornus) и куда по умолчанию, будут сохраняться результаты нашей работы

setwd(“D:\\site\\practicle22\\”)

Загрузим пакет R diversitree – этот пакет содержит набор функций, позволяющих оценить скорость диверсификации и протестировать разные модели. Если этот пакет неустановлен, воспользуйтесь командой install.packages(“diversitree”)

require(diversitree)

Загрузим филогенетическое дерево. Дерево должно быть доступно в формате nexus.

phy<-read.nexus(“tree.nex”)

Визуализируем дерево

plot(phy)

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

traits<-read.csv(“traits.csv”,stringsAsFactors=F)

Трансформируем загруженый файл из формата data.frame в формат vector

data<-traits[,2]

Присвоим имена вектору. Это обязательный шаг. Таким образом функции пакета diversitree будут знать соответствие между значениями в векторе data и филогенетическим деревом

names(data)<-traits[,1]

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

lik <- make.bisse(phy, data)

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

p <- starting.point.bisse(phy)

В пакете diversitree параметры модели (скорость видообразования, вымирания и переход из состояние 1 в 2) кодируются в следующей последовательности: скорость видообразования для видов, с признаком 1 (lambda0), скорость видообразования для видов с признаком 2 (lambda1), скорость вымирания для видов, с признаком 1 (mu0), скорость вымирания для видов с признаком 2 (mu1), скорость перехода из состония 1 в 2 (q01) и скорость перехода из состояния 2 в 1 (q10). Эти начальные параметры выбираются случайным образом, обычно имеют порядок меньше 1. В процессе оптимизации модели эти параметры будут постепенно меняться, и в результате процесса оптимизации методом максимального правдоподобия будут найдены новые, оптимальный параметры для конкретного филогенетического дерева и морфологического признака.

Теперь оптимизируем модель с помощью метода максимального правдоподобия

fit <- find.mle(lik, p)

Посмотрим результат

fit

Итак, перед нами результат оптимизации модели. Нас интересует несколько параметров. Во-первых, logLik модели равен -230.1526, мы должны записать этот параметр, чтобы потом иметь возможность сравнить его с logLik других моделей.

Эта модель является наиболее сложной, так как одноврменно оценивает все шесть параметров. Это можно увидеть в блоке fit$par – где мы увидем различные значения для скорости видообразования, вымирания и переходов для обоих состояний (с и без шпорца).

Как мы можем видеть, для видов без шпорца скорость видообразования равна lambda0 = 0.12028, в то время как для видов со шпорцем скорость видообразования lambda1 в 7 раз выше и равна 0.73893. Кроме того, виды без шпорца вымирают (mu0) со скоростью в 66 раз (!) выше (0.29267), чем виды несущие шпорец (mu1 = 0.00374). Таким образом мы можем расчитать скорость диверсификации. Она равна разнице между скоростью видообразования и вымирания:

lambda – mu = ro

Для видов со шпорцем, скорость диверсификации равна lambda1 – mu1 = 0.73893 – 0.00374 = 0.73519, а для видов без шпорца диверсификация равна lambda0 – mu0 = 0.12028 – 0.29267 = -0.17239. Диверсификация для видов без шпорца негативная, и мы можем сделать вывод о том, что данные виды находятся под угрозой постепенного вымирания, в то время как виды со шпорцем диверсифицируют более чем в 7 раз быстрее.  Кроме того, скорость перехода из состояния без шпорца в состояние со шпорцем превышает скорость обратных переходов в десятки раз (q01 и q10)

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

lik.l <- constrain(lik, lambda1 ~ lambda0, mu1 ~ mu0, q01~q10)

Таким образом, в модели будут оценены только три параметра, вместо шести и будет соблюдено условие равенства между парами параметров для видов с и без шпорца

Теперь оценим параметры модели с помощью метода максимального правдоподобия, как ранее мы оценили для более сложной модели

fit.l <- find.mle(lik.l, p[argnames(lik.l)])

Посмотрим результат анализа.

Наиболее интересным для нас является logLik модели, который равен -250.2319. Сравним этот logLik с тем, который мы получили для более сложной модели (-230.1526). В целом, модель у которой logLik больше, лучше описывает данные. В данном случае, logLik больше для более сложной модели. Однако, мы помним, что в одной модели оценивались шесть параметров, а в другой – только три. Таким образом, сравнение logLik напрямую не является корректным и нам нужно пересчитать logLik в AIC  для того, чтобы учесть число параметров (сложность) модели.

Пересчет делается по формуле:

-2*(число параметров) – 2*logLik

Таким образом, получаем:

AIC1 =  -2*6 – 2*-230.1526 = 494.4638

AIC2 =  -2*3 – 2*-250.2319 = 448.3052

Наименьшее (в отличие от наибольшего для logLik!) значение AIC обозначает наилучшую модель. Статистически значимой считается разница между моделями как минимум в четыре единицы. В нашем случае, наименьшее AIC имеет более сложная модель (шесть параметров, полная модель) и она отличается более чем 45 пунктов от следующей модели (нулевой). Таким образом, мы отвергаем нулевую гипотезу о равенстве скорости диверсификации для видов с и без шпорцами, и говорим о том, что виды со шпорцами диверсифицируют примерно в 7 раз быстрее, чем виды без шпорца. 

Признак-зависимая диверсификация: методология анализаCharacter-Associated Diversification: theory

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

Признак-зависимая диверсификация  (Maddison, Midford & Otto, 2007) выражается в наличии различных скоростей видообразования и вымирания видов в зависимости от наличия или отсутствия у них некоторого морфологического или экологического признака. Признаки могут быть как дискретными, так и непрерывными.

Для простоты рассмотрим пример бинарного дискретного признака.

Представим филогенетическое дерево, в котором для каждого ныне живущего вида известен некоторый дискретный бинарный признак. В качестве примера такого признака используем гипотезу S. Hodges в 1997 г.  Исследователь предположил, что шпорец у рода Аквилегий (семейство Лютиковые) является ключевой инновацией, предопределившей ускоренное видообразование и диверсификацию в группе. Нашим бинарным признаком становится “шпорец” и “отсутствие шпорца”.

Аквилегия - пример цветка со шпорцами

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

Пример филогенетического дерева с ярко выраженной разницей в скоростях видообразования между кладами

При оценке скорости диверсификации клады оцениваются как минимум два параметра: скорость видообразования и скорость вымирания (обозначаются буквами λ и μ соответственно). Разница между скоростью видообразования и вымирания назвается диверсификацией (r) и их дробь – видовым обновлением (t).  При сравнительном анализе клад с различными дискретными бинарными признаками оценивается статистическая вероятность того, что скорости видообразования и вымирания неравны между кладами, а также третий параметр – скорость перехода из одного бинарного состояния в другое (обозачается как σ).

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

Максимально сложная модель диверсификации при наличии дискретного бинарного признака выглядит следующим образом:

λ1 ≠ λ2 AND μ1 ≠ μ2 AND  σ12 ≠  σ21

где 1 и 2 обозначают состояние бинарного признака – в случае со шпорцем его наличие или отсутствие;

λ1 – скорость видообразования группы, где виды несут шпорец;

λ2 – скорость видообразования группы, где виды не несут шпореца;

μ 1 – скорость вымирания группы, где виды несут шпорец;

μ 1 – скорость вымирания группы, где виды не несут шпорец;

σ 12 – скорость потери шпорца (переход из состояния 1 в 2)

σ 21 – скорость приобретения шпорца (переход из состояния 2 в 1)

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

λ1 = λ2 AND μ1 = μ2 AND  σ12 = σ21

Такая модель сводится к оценке всего трех параметров вместо шести, как, например в модели выше.

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

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

λ1 = λ2 AND μ1 = μ2 AND  σ12 = σ21

λ1 ≠ λ2 AND μ1 = μ2 AND  σ12 = σ21

λ1 = λ2 AND μ1 ≠ μ2 AND  σ12 = σ21

λ1 = λ2 AND μ1 = μ2 AND  σ12 ≠ σ21

λ1 = λ2 AND μ1 ≠ μ2 AND  σ12 ≠ σ21

λ1 ≠ λ2 AND μ1 ≠ μ2 AND  σ12 = σ21

λ1 ≠ λ2 AND μ1 = μ2 AND  σ12 ≠ σ21

λ1 ≠ λ2 AND μ1 ≠ μ2 AND  σ12 ≠ σ21

Сравнение моделей должно проводиться с использованием параметра AIC, а не Log-Likelihood, так как число параметров отличается между моделями и необходима коррекция.

В следующей заметке мы рассмотрим, каким образом можно протестировать данные модели с использование R и пакета DiversiTree.