Влияние топологии дерева на сдвиг оценки ТМRCA по внутрикладовой АSD.

[Овод]

Недавно, пытаясь (с подачи Вадима Веренича) "скрестить" ASD и ро-статистику в расчётах ВБОП по полученным в результате филогении деревьям, я "по-ходу" получил оценку для популяционного влияния дерева с заданной топологией на расчёт его возраста по ASD. Хотя полученная формула имеет лишь косвенное отношение к первоначальной затее, она может оказаться даже интересней в свете недавних разговоров о "популяционном смещении" в ASD-методах.
 
Итак, пусть задано произвольное дерево со следующими параметрами:
 
N - число современных гаплотипов на нём
Т - возраст дерева (число поколений от корня до вершины)
L_i и K_i - cоответственно длина (в поколениях) и кратность i-того ребра дерева.
 
Здесь под "ребром" понимается участок дерева без ветвлений на нём, а под "кратностью" - число финальных потомков, к которым ведёт этот участок.
 
Для начала рассмотрим формулу для "смещенной" ASDc = SUM_ij(a_i-a_j)²/N², переводя которую к гистограммному виду (по частоте аллелей в выборке) и рассматривая случай единичной мутации на любом i-том ребре с кратностью К_i получаем приращение дисперсии за счёт этой мутации равное 2*K_i(1-K_i/N)/N. Усредняя этот результат по всем рёбрам (с весом их длины), получаем, для средней  ASDс в результате ожидаемого числа мутаций на дереве со средней скоростью m мутаций/поколение ASDc = 2mSUM_i(L_iK_i(1-K_i/N)/N), и, учитывая, что для любого дерева SUM_i(L_iK_i/N) = T, имеем итоговую оценку смещения относительно истинного возраста Т:
 
Т - ASD_С/2m = SUM_i(L_iK_i²/N²)                            (1)
 
Однако, для "смещённой" дисперсии мы не можем чётко выделить именно "популяционное смещение" , так как чисто статистически эта оценка также смещена. Чтобы его всё же найти, рассмотрим формулу "несмещённой"  ASD:
 
ASD_Н = SUM_ij(a_i-a_j)²/(N(N-1)) = ASD_С*N/(N-1)
 
и, домножая обе части уравнения (1) на N/(N-1), после нехитрых преобразований, получаем:
 
Т - ASD_Н/2m  = SUM_i(L_iK_i(K_i-1))/(N(N-1))

что и является искомой формулой чисто "популяционного" смещения для внутрикладовой ASD. Интересны два факта, вытекающих из этой формулы:
 
1. Терминальные ветви дерева, кратность которых всегда равна 1 не влияют на смещение дисперсии, давая нулевой результат, а всё "зло" кроется исключительно в общих ветвях, имеющих двух или более современных потомков из выборки. Собственно, их наличие и является истинной причиной популяционного сдвига.
 
2. Смещение квадратично возрастает с увеличением кратности ветвей и линейно - с их длиной. То есть, наличие большего числа выживших потомков предка имеет гораздо больший эффект, чем его близость по возрасту к современникам, которая тоже увеличивает смещение - но в меньшей степени.
 
Надеюсь, что эта формула может помочь, если не в конкретных расчётах и прикидках, то в понимании сути "популяционного эффекта" - наверняка.
 
P.S. Eщё раз подчеркиваю, что величины ASD_C и ASD_H в приведённых здесь формулах следует понимать как матожидание (среднее) ASD, а не любое (выборочное) её значение.

[Nimissin]

Уважаемый Овод, полученная Вами формула мне очень нравится. Но пока чисто интуитивно. Возьму некоторое время для изучения.

[shekhol]

Спасибо, Овод!

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

[Овод]

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

Наоборот, омоложенное.

[Nimissin]

Появились следующие вопросы:

1. Вызвает сомнение верность формулы
SUMi(LiKi/N) = T.
Похоже, она верна только для случая, когда все ветви растут непосредственно от общего предка. В иных случаях после суммирования должен получаться средневзвешенный по кратности ветвей Ki возраст ветвей <L>, меньший чем T.

2. Рассмотрим самый простой случай двух ветвей, равных по кратности и происходящих от самого общего предка. Тогда L1 = L2 = T, K1 = K2 = N/2. Предположим также, что N велико и вместо (N-1) можно использовать N. Подставим значения в формулу для смещенного ASDc:

Т - ASDс/2m = T (1/4 + 1/4) = T/2.

Тогда ASDс/2m = T/2.

Результат неверный, так как из начальных условий должно быть

ASDс/2m = T.

Для того, чтобы "звездообразные" деревья (т.е. все ветви начинаются от общего предка) не имели какого-либо "популяционного смещения", по-видимому, в формуле должен содержаться множитель типа (T - Li).

[dima75]

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

Наоборот, омоложенное.

Звучит логично даже без формул: понятно, что без внесения поправок ВБОП будет омоложен в сторону самой большой ветви.

[Овод]

1. Вызвает сомнение верность формулы
SUMi(LiKi/N) = T.
Похоже, она верна только для случая, когда все ветви растут непосредственно от общего предка. В иных случаях после суммирования должен получаться средневзвешенный по кратности ветвей Ki возраст ветвей <L>, меньший чем T.

Формула верна абсолютно для любого мыслимого дерева. Ещё раз подумайте, и если придумаете хотя бы одно дерево, для которого это соотношение нарушено, то я Вам выдам персональную нобелевскую премию в виде ящика коньяка. Или хотите равноправное пари?

2. Рассмотрим самый простой случай двух ветвей, равных по кратности и происходящих от самого общего предка. Тогда L1 = L2 = T, K1 = K2 = N/2. Предположим также, что N велико и вместо (N-1) можно использовать N. Подставим значения в формулу для смещенного ASDc:

Т - ASDс/2m = T (1/4 + 1/4) = T/2.

Тогда ASDс/2m = T/2.

Результат неверный, так как из начальных условий должно быть

ASDс/2m = T.

Для того, чтобы "звездообразные" деревья (т.е. все ветви начинаются от общего предка) не имели какого-либо "популяционного смещения", по-видимому, в формуле должен содержаться множитель типа (T - Li).

Результат верный. Для двух рёбер и гаплотипов может быть только два. И уж никак нельзя считать в этом случае, "что N велико". Вы употребили "смещённую" формулу и получили ожидаемое 50%-ное статистическое смещение, равное для двух гаплотипов T/2.
 
В таких случаях можно считать только по "несмещённой" формуле. Подставьте свои цифири в неё и убедитесь, что она даёт верный ответ, то есть ASD_H /2m = Т.
 
Я вообще рекомендую народу считать по "несмещённой" формуле всегда, хоть она и не слушком удобна. И уж, во всяком случае, для N<100, чтобы ошибка не превысила 1%.

[Yurgan]

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

Наоборот, омоложенное.

Согласен.

[Каржавин]

Итак, пусть задано произвольное дерево со следующими параметрами:
N - число современных гаплотипов на нём
Т - возраст дерева (число поколений от корня до вершины)
L_i и K_i - cоответственно длина (в поколениях) и кратность i-того ребра дерева.
Здесь под "ребром" понимается участок дерева без ветвлений на нём, а под "кратностью" - число финальных потомков, к которым ведёт этот участок.
Т - ASD_Н/2m  = SUM_i(L_iK_i(K_i-1))/(N(N-1))

Правильно ли я понимаю, что если все ветки начинаются с первопредка и не ветвятся, то кратность К становится равной N, суммирование вырождается в один член L=T, так как Ki(Ki-1)/(N(N-1) =1. При этом получается, что
T - ASD/2m = T

[Овод]

Итак, пусть задано произвольное дерево со следующими параметрами:
N - число современных гаплотипов на нём
Т - возраст дерева (число поколений от корня до вершины)
L_i и K_i - cоответственно длина (в поколениях) и кратность i-того ребра дерева.
Здесь под "ребром" понимается участок дерева без ветвлений на нём, а под "кратностью" - число финальных потомков, к которым ведёт этот участок.
Т - ASD_Н/2m  = SUM_i(L_iK_i(K_i-1))/(N(N-1))

Правильно ли я понимаю, что если все ветки начинаются с первопредка и не ветвятся, то кратность К становится равной N, суммирование вырождается в один член L=T, так как Ki(Ki-1)/(N(N-1) =1. При этом получается, что
T - ASD/2m = T

Нет, не правильно. Внимательно прочтите определение кратности ребра, приведённое в цитирумом Вами участке (выделено мной жирным шрифтом), и Вы поймёте, что "кратность" - это не количество рёбер равной длины, а количество современников, к которым ведёт данное ребро.
 
В Вашем примере каждое ребро ведёт к одному потомку  и кратность каждого равна 1, а длина - количеству поколений до предка. И, следовательно, правильный ответ:
 
Т-ASD_Н/2m = 0

[Каржавин]

Как я понимаю, Li берутся из филогенетического дерева, в котором эти Li на самом деле не соответствуют истинным длинам (короче).
Я нечто похожее вводил в оценку, причем, данные брал из дерева, которое мне известно было полностью (коли моделировал, значит, полностью известна структура дерева). Действительно, смещение оценки ЗНАЧИТЕЛЬНО уменьшалось, но увеличился разброс оценки. А ведь это самое главное - разброс. Смещение кое-как можно и калибровочными апроксимационными кривыми компенсировать.
Поскольку предлагаемая Вами формула (как и мои экзерсисы на данную тему) эмпирическая и не выведена из соответствующей статистической модели, хоть и тенденции она отслеживает в верном направлении, то увеличение разброса оценки я Вам гарантирую.
Я штуки 4-5 аналогичных формул учета филогении пробовал. Вывод однозначен: нужно, чтобы формула вытекала хоть из какой-то правдоподобной статистической модели и из соответствующих вероятностных распределений.

[VVR]

Недавно, пытаясь (с подачи Вадима Веренича) "скрестить" ASD и ро-статистику в расчётах ВБОП по полученным в результате филогении деревьям,

Не помню о каком сообщении Вадима речь, но если вдруг Вы перепутали В.В. с VVR, а ASD с ММП (http://forum.molgen.org/index.php/topic,1682.msg67422.html#msg67422) , то попробую объяснить, что я имел ввиду, тем более Каржавин тоже коснулся вопроса топологии.
Я имел ввиду при построении дерева применять ММП для проверки(уточнения) филогении, а затем на основе максимально правдоподобного дерева рассчитывать ро-статистику.

[Овод]

Нет, я имел в виду именно Вадима. Он просил меня здесь http://forum.molgen.org/index.php/topic,744.msg47029.html#msg47029 дополнить мой скрипт по ро-статистике расчётом ASD ещё весной. Но я что-то закопался и так его и не закончил (зато получил формулу сдвига).
 
Вы же говорили о ММП.

[VVR]

Понятно

[Овод]

Я имел ввиду при построении дерева применять ММП для проверки(уточнения) филогении, а затем на основе максимально правдоподобного дерева рассчитывать ро-статистику.

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