Калибровка скоростей мутаций для целей ДНК-генеалогии

[kaa76]

Если я ничего не путаю, то с возрастом человека мутации по Y накапливаются. Интересно, появится ли со временем дополнительное поле в Y-гаплотипе - ВозрастОтца?

[Каржавин]

По теме использования известных генеалогий для калибровки частот ("скоростей") мутаций хотел бы добавить вот еще что.
1. Два года назад мы с Дмитрием провели сравнительное исследование классического ASD метода и линейного метода. Пожалуй, главным варьируемым параметром был объем выборки гаплотипов, который использовался для вычисления ВБОП. Результаты описаны в статье Д.Адамов, С.Каржавин. О влиянии численности популяции на точность оценки TMRCA, получаемой стандартными методами по совокупности аллелей STR локуса. RJGG, том2,№2, 2010.
Оказалось, что вычисляемый возраст начального предка СИСТЕМАТИЧЕСКИ ЗАНИЖАЕТСЯ, причем, данное занижение увеличивается с уменьшением количества гаплотипов. Напоминаю, что Рожанский-Клесов как раз и использовали малые совокупности гаплотипов для своих калибровок "скоростей" мутаций. Причину этого занижения мы полностью осознали несколько позднее, а заключается она в очень простом и очевидном факте. Чем меньше по численности выборка гаплотипов, тем вероятнее, что общий начальный предок этой совокупности гаплотипов моложе истинного родоначальника исследуемого субклада или предполагаемого начального предка. И чем меньше количество гаплотипов, тем В СРЕДНЕМ моложе общий начальный предок исследуемой совокупности гаплотипов. Например, имеем мы всего ДВА каких-то гаплотипа Рюриковичей. Ну и какова вероятность того, что каждый из них является потомком разных сыновей самого Рюрика? Конечно, эта вероятность очень мала. А теперь представим, что у нас миллион потомков Рюрика (ограниченное количество официальных потомков и подавляющее большинство от внебрачных связей Рюриковичей во всех поколениях  ). Очевидно, что из такой совокупности гаплотипов ПРАКТИЧЕСКИ ВСЕГДА найдется хотя бы пара гаплотипов, каждый из которых произошел от непосредственных сыновей (официальных или внебрачных) Рюрика. ТО есть, чем больше объем исследуемой выборки гаплотипов, тем "ближе" мы можем подобраться к возрасту истинного начального предка.
Дело усугубляется еще и тем, что чем старше используемая в калибровках знаменитая генеалогия, тем больше вероятность, что к настоящему времени большинство ветвей, "вырастающих" от непосредственно родоначальника, пресеклись, и все современные гаплотипы начали "ветвиться" от какого-то промежуточного потомка.
2. Что сделали в своих калибровках Рожанский-Клесов? Они ПОСТУЛИРОВАЛИ, что исследуемые ими совокупности гаплотипов знаменитых родов с вероятностью ЕДИНИЦА в качестве общего начального предка имеют истинного родоначальника династии, а не промежуточного потомка. Таким образом, они реально более короткую во времени генеалогию (от промежуточного предка) растянули на весь период времени от истинного родоначальника, и, тем самым, количество мутаций, которое произошло на более коротком промежутке времени. растянули на весь период времени. Очевидно, что при этом "скорость" мутирования получается НИЖЕ ИСТИННОЙ.
Очевидно, единственно правильным является использование частот ("скоростей") мутаций, полученных по парам гаплотипов "отец-сын". НУ НЕТ ПОКА ДРУГОГО ПРАВИЛЬНОГО способа!
3. Теперь замечание по частому вопросу. Факт периодического ЗАНИЖЕНИЯ результатов по линейному методу над классическим методом объясняется довольно просто: выборки гаплотипов для знаменитых генеалогий маленькие, а в этом случае разброс оценки ВБОП очень большой (см. рис.8,10 указанной выше статьи). Во-вторых, линейный метод дает лучшее приближение оценки к истинному возрасту начального предка, но расплата за это - больший разброс оценки (см. рис.8,10). Вот и получается, что периодически оценка по линейному методу ниже, чем у классического метода.

[Каржавин]

Посмотрел на "Родстве" последние посты Рожанского и Клесова по проблеме калибровки частот ("скоростей") мутаций: http://www.rodstvo.ru/forum/index.php?showtopic=127&pid=92008&st=460&#entry92008
Меня заинтересовало следующее утверждение А.Клесова
Скорости мутаций откалиброваны не ДЛЯ линейного метода, но С ИСПОЛЬЗОВАНИЕМ линейного метода, но если взять для калибровки квадратичный, то будет то же самое. Для любой не "кривой" выборки линейный и квадратичный метод обязаны совпадать. В их основе - те же самые принципы, то же распределение Пуассона, те же принципы химической кинетики.
Вся эта возня вокруг "разницы" отражает низкую квалификацию тех, кто возится, их нежелание посмотреть на основы методов шире, и их нежелание читать литературу. Всё это было проанализировано Д.Адамовым и мной еще в далеком 2008 году, и опубликовано в серии статей в Вестнике.

Во-первых, очередное заклинание уважаемого А.Клесова про распределение Пуассона и принципы химической кинетики становится уже даже не смешным.
Во-вторых, видимо, в категорию "низкоквалифицированных", получается, попал не только я, но и ... Д.Адамов. Еще ДВА года назад мы именно Д.Адамовым совместно написали статью (название ее в моем предыдущем посте), в которой проанализированы эти самые классический квадратичный метод и линейный (автор - Д.Адамов). В этой работе выводы получены ОДНОЗНАЧНЫЕ: и квадратичный метод и классический имеют систематическое СМЕЩЕНИЕ оценки ВБОП, которое зависит от объема родовой популяции. Игнорирование уважаемыми И.Рожанским и А.Клесовым полученных ОДНОЗНАЧНЫХ и ОЧЕВИДНО ИНТЕРПРЕТИРУЕМЫХ результатов - это просто глупо.
Приходится еще раз повторить, что проблема калибровки по известным генеалогическим родам (Рюриковичи, Гедеминовичи и пр.) в том, что с учетом даже их неизвестных внебрачных детей их современные  родовые популяции относительно немногочисленны. Конечно, определить их численность не представляется возможным, но вряд-ли какая-либо из них превышает 100 тыс. мужчин. Исключение можно сделать, пожалуй, для Аднанов, существующих довольно давно. А для генеалогий, начавшихся позднее XIII в., их численность просто не успела вырасти. Что из этого следует? А следует то, что в соответствии с рис.7 и рис.9 из указанной статьи необходимо для объема родовой популяции менее 100 тыс. вводить соответствующие поправки по апроксимационным формулам (6) и (9), полученным Д.Адамовым (на которого так любит ссылаться А.Клесов). Только для популяций объемом не менее 100 тыс. и классический квадратичный и линейный методы дают близкий результат с практическим отсутствием смещения оценки ВБОП. Но для молодых генеалогий исторических персонажей (коих довольно много в списке у Рожанского) вряд-ли численность превышает тысячу. Очевидно, что в этом случае и линейный метод дает довольно приличное смещение (но значительно меньше, чем квадратичный, поэтому он и был Д.Адамовым придуман).
Вообще говоря, Рожанский с Клесовым попали в ловушку, которую поставили сами себе, опубликовав и распропагандировав весьма несовершенный метод калибровки частот ("скоростей") мутаций. Причем, не просто проигнорировав, а одновременно с этим весьма отрицательно публично охарактеризовав метод калибровки по парам галпотипов "отец-сын". Как они в дальнейшем будут выпутываться из этой дурацкой ситуации, я лично буду с удовольствием наблюдать.

Далее из комментария Рожанского:
И еще раз о калибровках линейного и квадратичного методов, раз уж эта тема приобрела какой-то нездоровый характер...
А вот здесь я не могу согласиться с уважаемым И.Рожанским. Эта тема как раз стала носить здоровый характер после того, как началась разумная критика используемых методов калибровки 
Вместе с тем, хочется отметить, что идея И.Рожанского калибровать по генеалогиям известных исторических персонажей весьма оригинальна и интересна. Но она была очень хороша несколько лет назад, когда еще статистики по парам гаплотипов "отец-сын" было еще мало. Но в настоящее время эта методика значительно хуже по качеству. Сейчас, несомненно, уже надо переходить на "прямой" метод оценивания частот ("скоростей") мутаций по парам гаплотипов "отец-сын".

[Agni]

Согласен, а вот какие всё-таки последние результаты по парам "отец-сын"?

[Alexander]

На самом же деле, как Вы видите на последнем дереве N1 для реальных генеалогий пришлось не понизить, а ещё больше увеличить скорости - практически вдвое выше клёсовских.

Попробовал прикинуть возраста МВП по общему дереву, но не скажу, что удачно получилось.

1. Гедиминовичей всего четверо. И по дереву ТНТ они разделяются, как 1+3 гаплотипа. Мурка делит их по другому - как 4 отдельные ветви (от узла - Гедимина?). На ТНТ-дереве Ро-возраст не проставлен (интересно совпадает ли с Ро-Мурки?). На мурочном дереве в узле стоит 496+/-138 (лет?). Если в годах, то маловато - явно промежуточный потомок. Оценка по калькулятору-МВП для ТНТ-шного расклада =0,0522, пересчёт по Клёсову даёт TMRCA=(0,0522/0,00179)*25=725 лет.а вот по мурочному - не уверен, что получилось - оценку-то дал (=0,070), но в результате оказалось прописано всего 3 ветви (а не 4, как я задавал) - это сбой калькулятора (?). Если принять эту оценку за верную, то получается чересчур много TMRCA=975 лет. Явно - не то.
Конечно, в этом случае очень уж мало гаплотипов.
Если взять mrate=0,00319 (как на дереве), не написано - сколько лет на поколение, то если 25 лет, то TMRCA=(0,0522/0,00319)*25=400 лет, если 30 лет, то TMRCA=480 лет (эта последняя величина близка к Ро-Гедимина на дереве).

2. Рюриковичей - оставил, как у меня было 8 гаплотипов. Структуры в ТНТ и Мурке - совпадают. В одной ветке - 7 шт., в другой - 1 шт. (47176 Gagarin). В паре у него стоит 167330 Volkonsky (Ольгович?), но у меня его не было, а во-вторых в "таблице князей" его гаплотипа почему-то не нашёл (все другие - в наличии. Есть пара Волконских R1a). Эти 2 линии сходятся на схеме условно на Ярославе Мудром, у которого Ро-возраста нет, но ниже его имеется узел с возрастом 1400+/-196 лет. Тем не менее результат получился - оценка=0,1023. В пересчёте по Клёсову - явный перебор TMRCA=(0,1023/0,00179)*25=1425 лет - хотя к узлу 1400 лет совсем близко. Перехлёст, или действительно первопредок не Ярослав, а более старый.
Если взять mrate=0,00319, то при 25 лет TMRCA=800, при 30 годах TMRCA=960 лет. Несколько меньше желаемого. Возможно, МВП было бы точней с учётом Волконского 167330 - подскажете, как отыскать его гаплотип?

Общее резюме - нужно бы доработать калькулятор, имеют место сбои в расчётах.

[Alexander]

Согласен, а вот какие всё-таки последние результаты по парам "отец-сын"?

По парам не хватает данных по частотам нескольких маркеров ftdna (порядка 7 шт. из 67), но для 60 маркеров вполне можно сделать калибровку, вряд ли она будет сильно отличаться от 67. Впрочем, это отличие тоже придётся как-то оценить.

[Mukovnikov]

... Возможно, МВП было бы точней с учётом Волконского 167330 - подскажете, как отыскать его гаплотип?

http://www.familytreedna.com/public/RussianNobilityDNA/default.aspx?section=ycolorized
номер 167330 ; в графе "Name" стоит только литера "V"

[Agni]

Вспоминаю свои старые знания по теории вероятностей, поэтому не взыщите ежели что не так.

Постановка задачи: По мне мутации STR  это процесс Бернулли.
Тогда имеем, что за дискретное время t процесс может оказаться в целочисленной точке х, причем с вероятностью p он прирастает  на 1, с q уменьшает на 1, с r остается на месте. p+q+r=1. Следовательно, надо найти такой набор чисел n,m,s соответствующих шагов чтоб было n-m=x. Тогда вероятность будет P(x,t) = [n m s]pⁿq^mr^s. Следовательно, максимум этой вероятности для известного x по t=n+m+s и будет наиболее вероятным временем (в поколениях) до исходной точки. Это для одного маркера. Для многих маркеров нужно найти максимум P₀(x₀,t)P₁(x₁,t)...P_i(x_i,t), где вероятности собственные для каждого маркера вычисленного независимо по парам "отец-сын", то есть просто количество людей у которых произошло данное изменение деленное на общее количество людей.
Чтоб найти максимумы как известно надо взять производную P по t и приравнять ее нулю, но можно и просто численно, все равно нужно искать максимум в целых точках. Если возникнут вероятности увеличения шага до 2 и более, то просто надо добавить их вероятности и увеличить мультиномиал.

Этот способ независит от количества маркеров в смысле его не надо калибровать для разного количества, на графике можно увидеть максимумы, мне кажется что он правильней биологически и математически. Извините, если открываю секрет полишенеля.

[Alexander]

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

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

[Alexander]

номер 167330 ; в графе "Name" стоит только литера "V"

Вот, спасибо, эту таблицу шерстил вдоль и поперёк, и поиском, и так - прям, затмение какое нашло.

Однако добавление 167330 V к расчёту МВП очень сильно увеличивает возраст - оценка МВП 0,1429 (это порядка 2000 лет). А только с одним V (без Гагарина) и вообще 0,1834. Явно что-то не то. На дереве этот V сидит очень хорошо, одинаково и в ТНТ и в Мурке. А расчёт выбивает...

Заодно по калькулятору:

прошу пока не использовать этот калькулятор в части линейного метода. Он еще весьма сырой, не доведенный до рабочего состояния.В линейном методе, как правило, не допускаются дробные значения базового гаплотипа. Расчет опирается на целочисленный аллель, который был принят за предковый.
В приведенном Вами модальном гаплотипе есть дробные значения. Семаргл обещал мне внести поправку.

В самом алгоритме всё верно - используются целочисленные значения (проверил тест в экселе). Просто нужно будет исправить написание базового в листке результатов.

И ещё вопрос к mouglley:
На дереве 12 Рюриковичей-ромбов, однако двоих из них нет в таблице гаплотипов (229414 Rus и 200703 Poland). Зато в таблице в разделе Рюриковичей есть 158597    Myshetsky  и 167322 Tbd. Но они на других страницах дерева, в совсем других ветвях. В чём тут дело?

[Alexander]

Нашёл эту пару Рюриковичей.
Почему-то в Russian Nobility DNA Project их всего 12 (включая 158597 Myshetsky и 167322 Tbd).
А в Rurikid Dynasty DNA Project их всего 11 (включая 229414    Lobanov-Rostovsky и 200703    Korybut Woroniecki), но без вышенаписанной пары и без V 167330.
Подскажите всё-таки, кто есть кто? Что за разночтение в проектах?
Или вернее всего верить дереву, исключив первую пару?

[mouglley]

Однако добавление 167330 V к расчёту МВП очень сильно увеличивает возраст - оценка МВП 0,1429 (это порядка 2000 лет). А только с одним V (без Гагарина) и вообще 0,1834. Явно что-то не то. На дереве этот V сидит очень хорошо, одинаково и в ТНТ и в Мурке. А расчёт выбивает...

Вот-вот.
И какая скорость теперь получается, близкая к той, что у меня?

Причём ветка довольно-таки отделённая от других.

[mouglley]

На дереве 12 Рюриковичей-ромбов, однако двоих из них нет в таблице гаплотипов (229414 Rus и 200703 Poland). Зато в таблице в разделе Рюриковичей есть 158597    Myshetsky  и 167322 Tbd. Но они на других страницах дерева, в совсем других ветвях. В чём тут дело?

229414 и 200703 в проекте уважаемого Andr
http://www.familytreedna.com/public/rurikid/

158597 Myshetsky не относится к основной линии Рюриковичей и его генеалогия подвергается сомнениям.
Скорее всего, общий предок с Рюриком у него ещё дальше, чем общий предок меня и Рюрика.

То же самое относится и к 167322 Tbd.
Хотя мой с 167322 предок жил относительно недавно.

[Каржавин]

Вспоминаю свои старые знания по теории вероятностей, поэтому не взыщите ежели что не так.

Постановка задачи: По мне мутации STR  это процесс Бернулли.
Тогда имеем, что за дискретное время t процесс может оказаться в целочисленной точке х, причем с вероятностью p он прирастает  на 1, с q уменьшает на 1, с r остается на месте. p+q+r=1. Следовательно, надо найти такой набор чисел n,m,s соответствующих шагов чтоб было n-m=x. Тогда вероятность будет P(x,t) = [n m s]pⁿq^mr^s. Следовательно, максимум этой вероятности для известного x по t=n+m+s и будет наиболее вероятным временем (в поколениях) до исходной точки. Это для одного маркера. Для многих маркеров нужно найти максимум P₀(x₀,t)P₁(x₁,t)...P_i(x_i,t), где вероятности собственные для каждого маркера вычисленного независимо по парам "отец-сын", то есть просто количество людей у которых произошло данное изменение деленное на общее количество людей.

В нашем форумном журнале RJGG в ряде статей (Д.Адамов, С.Каржавин) подробно рассмотрены вероятностные математические модели, на которых основаны (из которых ВЫВОДЯТСЯ) соответствующие статистические решающие правила по оцениванию ВБОП, Так что, все у нас стоит на твердой теоретической базе. Очень хорошо, что Вы разбираетесь в этих вопросах, будете в нашей команде. У нас много еще нерешенных вопросов в ДНК-генеалогии, которые ждут своих рабочих математических рук. Давайте, осваивайтесь потихоньку и включайтесь.

[Agni]

Вспоминаю свои старые знания по теории вероятностей, поэтому не взыщите ежели что не так.

Постановка задачи: По мне мутации STR  это процесс Бернулли.
Тогда имеем, что за дискретное время t процесс может оказаться в целочисленной точке х, причем с вероятностью p он прирастает  на 1, с q уменьшает на 1, с r остается на месте. p+q+r=1. Следовательно, надо найти такой набор чисел n,m,s соответствующих шагов чтоб было n-m=x. Тогда вероятность будет P(x,t) = [n m s]pⁿq^mr^s. Следовательно, максимум этой вероятности для известного x по t=n+m+s и будет наиболее вероятным временем (в поколениях) до исходной точки. Это для одного маркера. Для многих маркеров нужно найти максимум P₀(x₀,t)P₁(x₁,t)...P_i(x_i,t), где вероятности собственные для каждого маркера вычисленного независимо по парам "отец-сын", то есть просто количество людей у которых произошло данное изменение деленное на общее количество людей.

В нашем форумном журнале RJGG в ряде статей (Д.Адамов, С.Каржавин) подробно рассмотрены вероятностные математические модели, на которых основаны (из которых ВЫВОДЯТСЯ) соответствующие статистические решающие правила по оцениванию ВБОП, Так что, все у нас стоит на твердой теоретической базе. Очень хорошо, что Вы разбираетесь в этих вопросах, будете в нашей команде. У нас много еще нерешенных вопросов в ДНК-генеалогии, которые ждут своих рабочих математических рук. Давайте, осваивайтесь потихоньку и включайтесь.

Посмотрел, конечно поверхностно, Вашу статью 1 от 2009 года. Извините, если что-то недопонял
1. Конечно, схема Бернулли что Вы используете все таки не совсем корректна - там нет отсутствия мутации в поколении. У вас мутация есть всегда в дискретном времени, нет вероятности отсутствия мутации (но только для оценки вероятности мутаций, без поколений, она верна).
Но
2. Вы это компенсируете переходом к Пуассоновскому потоку, который как известно работает с непрерывным временем, и при выводе самой формулы используется именно вероятность отсутствия события.
Ординарностью в непрерывном времени Вы компенсируете появления нескольких мутаций в поколение. Поэтому, в принципе, использование формулы Пуассона оправданно.

Единственно, извините если у Вас это есть, а я не увидел, скорость мутации которая вроде у Вас обозначена тетой или лямбдой, это есть -ln p, где p - вероятность не наступления события в единицу времени т.е. отсутствия мутации в поколении (у меня обозначена как r).

Об остальном как нибудь буду думать.