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

[wertner]

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

Именно на этой идее построена моя программа, которую использую для построения деревьев. Но и в этом случае (как при использовании дерева с меньшим суммарным весом) необходимо использовать не максимально правдоподобное дерево, а семейство правдоподобных деревьев. Это делают и Мурка, и TNT, но не делает моя программа.

[Nimissin]

Уважаемый Овод, должен признать, что не могу понять формулу:

SUMi(LiKi/N) = T.

Из определения кратности Ki следует, что

SUMi(Ki) = N.

Подставляя вместо N сумму, получаем:

SUMi(LiKi)/SUMi(Ki) = T.

С другой стороны,

SUMi(LiKi)/SUMi(Ki) = <L>,

поскольку это общеизвестная формула для расчета среднего взвешенного. В нашем случае это среднее взвешенное значений Li с весами Ki.

Получается, что <L> = T.

Тогда выходит, что часть значений Li может быть меньше T, а другая часть - больше T.
Но из постановки задачи известно, что при любом i значение Li не может быть больше возраста предка T.
Следовательно, при любом i выполняется соотношение Li = T. Значит, все ветви генеалогического дерева растут от общего предка.

Что в моих рассуждениях не так?

[Nimissin]

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%.

Давайте посчитаем по "несмещенной" формуле. Рассматриваем случай L1=L2=T, K1=K2=N/2.

T - ASDн = L1*K1*(K1-1)/N/(N-1) + L2*K2*(K2-1)/N/(N-1) = (T/N/(N-1))*(N*(N-2)/4 + N*(N-2)/4) =
= T*(N-2)/(2*(N-1)).

Тогда T = ASDн только при N = 2.

Все-таки давайте уточним, что такое N. Я понимаю так, что это общая численность финальной популяции дерева.
Или это число ребер?

[Каржавин]

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

Что Вы предполагаете предпринять для оценки качества предлагаемой Вами формулы? А именно, насколько хорошо и при каких условиях компенсируется сдвиг, каков разброс оценки для разных скоростей мутаций, глубин генеалогических деревьев в поколениях и т.д. Если исследования покажут, что действительно по указанным параметрам Ваша формула дает выигрыш (особенно в разбросе (rms) оценки TMRCA), то это будет важный результат.
Как Вы знаете, в многие авторы вводят всякие поправки в различных процедурах статистического оценивания, и главная трудность - определение их качества.

[Овод]

Уважаемый Овод, должен признать, что не могу понять формулу:

SUMi(LiKi/N) = T.

Из определения кратности Ki следует, что

SUMi(Ki) = N.

Не следует. Возьмите простой пример: У деда два сына, каждый из которых родил ему два внука. Итого на дереве четыре финальных потомка и шесть ребер, два из которых (от деда к сыновьям) имеют кратность 2 (поскольку ведут к 2 внукам каждое), а четыре (от сыновей к внукам) - кратность по 1 (поскольку они ведут к одному внуку каждое). Итого по шести рёбрам имеем суммарную кратность 8, а не 4 (число финальной популяции внуков). Сумма кратностей вообще не имеет никакого физического смысла и, следовательно не может быть равна ни одному из остальных параметров топологии дерева.
 
Кстати убедитесь, что по первой формуле Вы и здесь получите возраст дерева 2 поколения (что соответствует истине). Но рисуйте любое другое мыслимое дерево, проставляйте на каждом его ребре длины и кратности и убедитесь в том, что первая формула (моя) верна всегда, а вторая (Ваша) - лишь в частном случае "звёздного" кластера.
 
И ещё раз внимательно прочитайте определение кратности ребра: это не число рёбер одинаковой длины, а число финальных потомков, к которым ведёт это ребро. И суммирование идёт по рёбрам, а не по потомкам.
 
Теперь понятно? Если - нет, я с удовольствием дам дальнейшие пояснения.
 

Все-таки давайте уточним, что такое N. Я понимаю так, что это общая численность финальной популяции дерева.

Это Вы правильно понимаете. Всё остальное - неверно.

[Nimissin]

Спасибо за разъяснение, уважаемый Овод.

В приведенном Вами примере L1=1, K1=2 - для ребра между дедом и его первым сыном,
L2=1, K2=1 - для ребра между дедом и его вторым сыном,
L3=L4=L5=L6=1, K3=K4=K5=K6=1 - для ребер между отцами (т.е. сыновьями деда) и их сыновьями.

T=(1*2+1*2+1*1+1*1+1*1+1*1)/4 = 2.
Все верно.

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

А как доказывается это замечательное соотношение в общем случае?

[Овод]

Что Вы предполагаете предпринять для оценки качества предлагаемой Вами формулы?

Поскольку правая часть формулы содержит только заранее заданные (определяющие) характеристики родового дерева, то статистическим колебаниям она не подвержена. Что же касается ASD (как оценки дисперсии), то разумеется, строго говоря, в формулах фигурирует её матожидание (иначе нельзя было бы ставить знак равенства). А стат. свойства диперсии  ASD уже были исследованы ранее (теми же Нордтведтом и Адамовым то же, если не ошибаюсь). А вообще-то стат. свойства АSD для гаплотипов были исчерпывающе исследованы ещё в прошлом веке Гольдштейном, Слаткиным и другими независимо.
 
Так что исследуемый мною сдвиг будет иметь те же самые свойства.

Другое дело, если дерево для сравнения не задано, а вычислено (скажем, в результате филогении). В этом случае стат.свойства сдвига будут зависеть от конкретного метода построения дерева.

[Овод]

А как доказывается это замечательное соотношение в общем случае?

Оно не может быть доказано, поскольку именно оно и определяет смысл входящих в него величин в алгебраической форме. И отражает соотношение между этими четырьмя независимыми параметрами топологии дерева по определению. Его просто нужно понять.
 
Для понимания оного приведу следующее объяснение. Представьте звёздный кластер - от одного первопредка к каждому потомку ведёт одна независимая линия-нить, длиной T, всего таких нитей - N. Общая длина всех нитей будет N*T. Теперь мысленно соберите все линии в "пучки", так чтобы получить дерево произвольной структуры, чтобы эти пучки образовали рёбра этого дерева. При этом очевидно, что ни длина каждой нити, ни общая их длина (в поколениях) не изменятся.
 
Поскольку, количество "нитей" в пучке, это и есть - кратность. то суммарная длина всех нитей и будет равна сумме длин нитей во всех пучках-рёбрах, помноженном на количество нитей в каждом пучке, то есть SUMi(LiKi)=N*T. Поэтому поделив эту величину на N, мы неизбежно получим T.

[Nimissin]

Да, это действительно так. Где бы мы не пережимали пучки нитей, общая длина будет N*T.

[Овод]

Ой не всегда. Есть тут кое-какие расчёты. Как-нибудь приведу. Там есть нюансы, но в подавляющем большинстве случаев - занижение.

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

[Овод]

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

[VVR]

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

[VVR]

Итак, модель. Популяция от общего предка развивается 35 поколений без популяционного эффекта, достигая 10 тыс.чел. Далее бутылочное горлышко, выживает 2 случайных потомка. Первый, через 65 поколений имеет 1млн. потомков, второй - 100 тыс. Поп. эффект за эти 65 поколений также отсутствует. Общий возраст всей конечной популяции численностью 1 100 000 чел. от общего предка 100 поколений. Таким образом, я пытался рассмотреть и бут.горлышко, и неравномерность развития дерева.
Рассматривал на одном маркере с интенсивностью(скоростью) мутаций 0,002 за поколение.
На первых 35 покол. распределение мутаций от предкового значения получилось таким.
+2 - 6 чел.
+1 - 327
0   - 9 334
-1 - 327
-2 - 6
Итого 10 тыс.чел.
Для проверки просчитал возраст ASD-методом-35,1 покол., линейным -35,09.
После бут. горлышка из них остаётся случайных 2 чел. Рассмотрев все возможные варианты и просчитав их вероятность, учитывая. что для дальнейших расчётов важно только расстояние между ними в мутациях, я свёл их к пяти вариантам:
4 мутации - вероятность 0,0000007
3- 0,0000785
2- 0,0043792
1- 0,1221794
0- 0,8733622
Замечу, что два последних варианта дают в данном примере более 99,5%.
Распределение на последующих 65 поколения такие:
Для большей части 1 000 000 чел.
+4-1
+3-38
+2-1 836
+1-57 299
0-881 652
-1-57 299
-2-1 836
-3-38
-4-1
Проверка возраста расчётом ASD-65,00 покол. ,лин-65,06
Меньшая часть - 100 000 чел.
+3-4
+2-184
+1-5 730
0-88 164
-1-5 730
-2-184
-3-4
Проверка расчётом возраста.ASD-65,02;лин-65,07.

[VVR]

Считаем возраст всей популяции 1 100 000чел. методом внутрикладовой ASD.Напомню, что реальный возраст до общего предка, который и должен получиться 100 покол.
4мут. между выжившими -726,17 покол.
3-436,91
2-230,29
1-106,33
0-65,00
Линейный метод -
4-309,45
3-236,79
2-171,37
1-112,51
0-65,01
Средний возраст(средневзвешенный) получается
ASD-70,8 покол.
лин.-71,34
Средневзвешенный по двум наиболее частым(более 95,5%) вариантам(0 и 1 мутация)
ASD- 70,07
лин.-70,88
В большинстве случаев, в данном примере 87,3%, оба метода дают возраст бутылочного горлышка - 65 покол.
Надо заметить, что в реальных расчётах разброс результатов будет компенсироваться количеством маркеров. При большом их количестве большая часть маркеров будет давать занижение, меньшая завышение. По всему гаплотипу, опять-таки в большинстве случаев, результат будет занижен из-за популяционного эффекта.

[VVR]

Теперь представим, что два выживших после бутылочного горлышка потомка имеют разные снипы, что позволяет их потомков надёжно разделить на субклады и посчитать интеркладовыми методами.
Кросс-ASD
4-4065 покол.
3-2315
2-1065
1-315
0-65
Разброс результатов огромный, но в большинстве случаев те же 65 покол. - возраст бутылочного горлышка.
А вот средневзвешенный получился 100,1 покол. Это совершенно не значит, что метод не подвержен популяционному влиянию, просто в данном приведенном примере заданный популяционный эффект не оказал популяционного влияния на оценку TMRCA этим методом.
Средневзвешенный по двум наиболее частым разностям(0 и 1 мутация) дал 95,7 покол.