SNPweights: использование модели калькулятора K16 для анализа главных компонентов происхождения

Ранее я уже отрапортовал о создании двух новых моделей для стандартного этно-популяционного калькулятора, в разработке которых использовались геномы людей, cамостоятельно указавшими свое происхождение (self-reported ancestry).
К сожалению, очень часто субъективная оценка собственного происхождения (указываемого респондентами в опросниках) недостаточно надежна для статистических методов анализа происхождения, поскольку некоторые люди либо сообщают ложные сведения о своей родословной или же просто не знают о своем истинном происхождении. Что еще хуже, — во многих публичных популяционных выборках мы не находим никаких  сведений о точном этническом составе людей в выборке . Как многие из вас знают,  существует множество способов достаточно точной оценки происхождения индивида на основе данных SNP генотипирования.

Самый простой способ сводится к следующему: сначала исследователь объединяет генотипы из своего исследования с генотипами образцов в референсной панели (например: HapMap или 1000 геномов),  затем находит пересечение SNP-ов в каждом наборе данных, а затем запускает программу кластеризации, чтобы увидеть, каким образом образцы исследования группируются с популяциями референсных панелей.  В принципе,  сам процесс несложный, но требует немало времени

К счастью, в 2014 году лабораторией Alkes была предложена программа которая, по сути, значительно облегчает процесс, выполняя большую часть работу за вас. Программа называется SNPWEIGHTS и можно скачать здесь.  Говоря простым языком, программа принимает  в качестве входных данных генотипы SNP-ов, самостоятельно находит пересечение генотипов SNP с генотипами в эталонной выборке , рассчитывает веса SNP-ов на основе предварительно настроенных параметров, чтобы построить первую пару главных компонентов (иначе говоря,  cобственных векторов), а затем вычисляет процентное значение происхождения индивидуума из каждой предковой популяции (кластера).

Для того, чтобы запустить программу, необходимо убедится в том, что в вашей системе установлен Python, и что ваши данные генотипирования приведены в формате EIGENSTRAT. Краткую инструкции по преобразованию в формат EIGENSTRAT с помощью инструмента convertf можно почитать здесь.  Данные аутосомного генотипирования FTDNA или 23andme можно напрямую преобразовать в формат EIGENSTRAT с помощью утилиты aconv от Феликса Чандракумара (либо любого самописного софта).

Затем необходимо загрузить сам пакет SNPWEIGHTS и референтную панель с весами снипов.

• Панель весов SNP для популяций Европы и Западной Африки можно скачать здесь.
• SNP веса для населения Европы, Западной Африки и  Восточной Азии можно скачать здесь.
• SNP веса для населения Европы, Западной Африки, Восточной Азии и популяций американских индейцев можно скачать здесь.
• SNP веса для популяций северо-западной, юго-восточной части Европы, ашкеназских евреев и можно скачать здесь.

Затем необходимо создать файл параметров par.SNPWEIGHTS с названиями входных файлов EIGENSTRAT, референтной панели, и файл c результатами. Например:

input_geno: data.geno
input_snp: data.snp
input_ind: data.ind
input_pop: CO
output: ancestry.txt

И, наконец, нужно запустиь программу с помощью команды inferancestry.py —par par.SNPWEIGHTS. Для того чтобы программа работала, убедитесь, что inferancestry.info и  файл референтной панели  находятся в том же каталоге, что и файл inferancestry.py.

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

Первый график — обычный график PCA c двумя первыми компонентами (собственными векторами) и наложенный на график процентный расклад компонентов происхождения:

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

Вот простой код генерирования этих графиков в R. В программе R нет базовых пакетов для построения триангулярных графиков, поэтому  нужно будет сначала установить пакет plotrix. Ancestry.txt  — это файл полученный на выходе из SNPWEIGHTS:

# EV Plot with Percent Ancestry Overlay
data=read.table("ancestry.txt", as.is=T, header=F)
names(data)
plot(data$EV1, data$EV2, pch=20, col="gray", xlab="EV1", ylab="EV2")
text(data$EV1, data$EV2,labels=round(data$EUR,2)100, cex=0.4, offset=0.1, pos=3)
text(data$EV1, data$EV2,labels=round(data$AFR,2)100, cex=0.4, offset=0.1, pos=2)
text(data$EV1, data$EV2,labels=round(data$ASN,2)*100, cex=0.4, offset=0.1, pos=1)
#Triangle Plot
data$total=data$EUR+data$AFR+data$ASN # Need to account
data$European=data$EUR/data$total # for slight rounding
data$African=data$AFR/data$total # in the ancestry
data$Asian=data$ASN/data$total # estimation file for
data_p=data[c("European","Asian","African")] # triax.plot to work
library(plotrix)
triax.plot(data_p, pch=20, cc.axes=T, show.grid=T)

 

Разумеется, размещенные на сайте разработчика референтные панели носят ограниченный характер. Поэтому я решил заполнить пробелы, преобразовав аллельные частоты SNP-ов в 16 предковых компонентах в 16 синтетических «чистых» предковых популяций, каждая из которых состояла из 200 синтетических индивидов («симулянтов») состоящих на 100 процентов из одного компонента происхождения в модели K16). Файл с генотипами 3200 «симулянтов» я использовал для вычисления весов снипов в каждом компоненте. Продвинутые пользователи, желающие протестировать модель K16 до ее публичного релизма, могут скачать полученный файл с весами снипов  здесь, а затем, cледуя приведенным выше инструкциям, использовать его в качестве референтной панели (а затем сравнить свои результаты с усредненными результатами разных этнических популяций).

Я протестировал веса снипов в модели K16 (выражаю признательность автору программу Чену за помощь), и обнаружил, что между данными калькулятора и данными SNPWEIGHTS расхождения носят незначительный характер, хотя похоже, что SNPWEIGHTS не так сглаживает минорные компоненты происхождения (что позволяет легче выделить в пространстве главных компонент кластеры):

Древние геномы человека в перспективе генетического разнообразия современных популяций

Примерно месяц тому назад, один из замечательных представителей «гражданской науки» в области генетики, известный геномный блоггер Polako (Давид Веселовски) разместил в своем блоге заметку, в которой были приведены результаты самостоятельного изучения вариативности снип-мутаций в пяти наиболее известных  из отсеквенированных геномов древних людей.  Хотя, как мне представляется, основное внимание Давид уделил все же прояснению ответа на вопрос о расположении  древнего генома сибирского мальчика со стоянки Malta (13 тысяч снипов-вариантов в аутосомах) в пространстве главных компонентов генетического разнообразия (PCA) cовременных человеческих популяций. К слову, этот же образец (Malta-1) был на днях включен в новую таблицу откалиброванных процентных соотношений 13 конвенциональных генетических компонентов в популярном среди пользователей Gedmatch этно-популяционногенетическом калькуляторе Eurogenes K=13 .  Наряду с вышеназванным образцом, в отреферированном анализе использовались геномные снип-варианты древнего ДНК австралийского аборигена (46 тыс.снипов), Anzick-1 генома древнего индейца культуры Кловис (106 тыс.снипов), генома древнего экскимоса Saqqaq (68 тыс.снипов), геном обитателя мезолитической Испании La-Brana 1 (23 тыс.снипов).

Можно предположить, что при проведении статистических анализов PCA, Давид использовал в качества сравнительного эталона-референса известный график из статьи Lazaridis et al. 2013.

PCA из статьи-препринта Lazaridis et. al .2013.

К сожалению ,  Давид из Eurogenes по определенным причинам не включил в свой анализ варианты снипов остальных известных евразийских древних геномов задействованных в PCA-анализе статьи-препринта Lazaridis et al. 2013, в частности древние геномы неолитического периода — женщин  культур воронковидных кубков (Swedish_farmer) и культуры линейно-ленточной керамики Южной Германии (Stuttgart), а также неолитического жителя Тирольских Альп — Этци (Iceman). Нет в  анализе Давида и образцов мезолитического и эпинеолитического генофонда Европы — мезолитических охотников-собирателей Motala  и Losсhbour и неолитических охотников с острова Готланд (Skoglund_merge). C другой стороны, в широко обсуждаемой предварительной версии статьи Лазаридиса к анализу привлечены только актуальные в евразийской перспективе образцы, и поэтому на графике PCA отсутствуют геномы древнего аборигена Австралии и двух древних геномов из Северной Америки.

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

Ниже в таблице приведены сводные данные о древних геномах и размерности числа снипов  этих образцов, которые использовались в моем анализе

Аncient (Afontova Gora) 10965

Australian Aborigen 236880

Otzi_Tyrolean 171195

Swedish_merged_farmer 1600

Swedish_merged_HG 4053

La Brana  57050

Malta-1 44459

LBK_Stuttgart 54220

Motala12 54677

Loschbour 54591

Motala_merged 35010

Визуализация двух первых главных компонентов разнообразия в популяциях выборки

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

Ниже я разместил серию визуализаций графика PCA. Первая иллюстрация — визуализация двух первых главных компонентов разнообразия, ставшая уже классической форма V-образного клина.

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

Центроиды популяций

 Как видно из второго графика, мировый популяции равномерно распределились по углам триангуляции. Африканские популяции длинным шлейфом-вектором  от пигмеев до фулани, cахарцев и эфиопских этносов распредились в левой части V-клина. Между ними и европейцами находится большая группа смешанных рассовых групп — пуэрто-риканцы, доминиканцы, афроамериканцы Карибского региона и Северной Америки, морокканцы, мозабиты и жители Туниса. В вершине угла V клина находятся все классические европейские этнические группы и народности. Они образуют внутренний европейский градиент генетической вариативности, уменьшающийся по мере удаления на север.  Северные популяции европейцев (особенно в Скандинавии и Прибалтике) смыкаются с находящимися на самой веришине угла древними геномами европейцев времен мезолита (Motala, Loschbour, La Brana,и перехода к неолита. Эта картина соответствует тому, что мы наблюдаем на графике Lazaridis et al. 2013.  Наблюдаемая на моем графике более значительная дистанция шведских охотников-собирателей шведской культуры ямочной керамики от современных популяций северной Европы объясняется только тем, что в работе Lazaridis et al. 2013 использовалась большее количество тех снипов древних геномов, которые встречаются и в современных популяциях (т.е находятся в пределах современной вариативности генов жителей современной северной Европы). Поэтому дистанция в узказанной работе между древними и современными популяцими ниже (тот же феномен наблюдается и в неолитическом векторе). Неолитический «вектор» представлен шведским неолитическим фермером, Этци Тирольцем, женщиной из неолитического поселения возле современного Штуттгарта. Из современных популяций к этому вектору находятся близко сардинцы и баски.
 

Однако наиболее интересная картина наблюдается в правой части графика, где мы наблюдаем наложение сразу нескольких клинов-градиентов разнообразия. Наиболее сложная структура наблюдается в том месте правого «крыла» графика, куда проецируются геномы двух палеолитических жителей Сибири (Malta-1 и AG). В этом месте график начинает ветвиться на три тесно переплетенные вектора-градиенты. Один уходит через Средную Азию-Непал-Северную Индию на юг, где встречается в двигающимся ему навстречу вектору-градиенту представленному австралийскими аборигенами, онге, папуасами, меланизийцами, андаманцами и дравидами.  Второй вектор ведет через Алтай-Монголию и Китай в Индокитай и юго-восточную Азию.

Третий вектор разделяется сразу на две части — одна ведет к палеосибирским народами и далее к алеутам и экскимосам. Этот вектор заканчивается древним геномом Saqqaq, который видимо является самым чистым «образчиком» генома древних людей, связанных с этими группами. Второй уходит через группу североамериканских индейских народов на юг, в Мезоамерику и далее к индейцам южной Америки. Вектор заканчивается на Anzick-1, и — по аналогии c Saqqaq, — можно сделать вывод о том, что этот геном является квинтэссенцией «чистого америндского компонента» без позднейших вкраплений в ходе контактов с европейцами.

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

Трехмерная перспектива PCA

О так называемом «эффекте калькулятора»

Уважаемый Сергей Козлов написал короткую заметку посвященную так называемому эффекту калькулятора.

«Эффект калькулятора» и его влияние на результаты оракулов аутосомных этно-калькуляторов.

«Эффект калькулятора» и его влияние на результаты оракулов аутосомных этно-калькуляторов.

На Молгене и других форумах периодически всплывает тема несоответствия результатов реальных людей и предсказаний встроенных в этно-калькуляторы оракулов. Недоброжелатели обычно приводят в качестве объяснений плохие исходные выборки либо неправильность самого подхода. Я же считаю, что многие аутосомные калькуляторы работают отлично. Проблема кроется в оракулах — а именно, в «эффекте калькулятора».

Понятие «Calculator effect», которое принято переводить на русский язык, как «эффект калькулятора», было введено известным геномным блоггером Polako. На почве этого у него возникла длительная перепалка с не менее известным блоггером Dienekes Pontikos .
Суть эффекта в том, что результаты людей, включенных в расчет для выделения предковых компонентов (то есть участников его проекта), отличались от результатов людей, которые в проекте не участвовали.
На этом основании он отказывался делать оракулы на базе своих калькуляторов, так как их результат должен был неизбежно искажаться. Однако позже он все-таки выложил новый парный калькулятор Eurogenes JTest/EUtest (14 и 13 компонентов соответственно) вместе с оракулом, который давал хорошие предсказания для большинства европейцев. Чтобы добиться этого, Поляко использовал для расчета предковых компонентов исключительно научные выборки, а таблицы эталонов для оракула строил на основании данных участников своего проекта. Таким образом, ключ к решению проблемы — разделение выборок,
использованных для выделения компонентов, и выборок, использованных для построения таблицы эталонов.

В чем же причина «эффекта калькулятора»? Изучая таблицы эталонов разных калькуляторов и сравнивая их с данными реальных людей, я пришел к выводу, что различия между исходными выборками («эталонами») показываются более резкими, чем они есть между представляемыми популяциями на самом деле.
Содержание главного компонента в эталонах завышается, а содержание второстепенных, особенно близких к главному — занижается. Это естественно, если вспомнить, что сами различия выводятся на основании этих выборок. Приведу в пример один из лучших калькуляторов: MDLP World-22. Содержание основных европейских компонентов North-East European и Atlantic-Mediterranean Neolithic в эталонах и у реальных русских (в среднем):

Русские центра эталон 70,4 16
Русские центра реальные 63 19
Северные русские эталон 63,5 11,8
Северные русские реальные 60 15
Украинцы эталон 64,2 19,7

Реальные русские не могут достичь показателей своего эталона и показываются более удаленными украинцами. В смешанном режиме оракул может показать их, к примеру, на 2/3 русскими и на 1/3 немцами — у немцев North-East European ниже, а Atlantic-Mediteranean выше.
Аналогично, например, у финнов специфически финский компонент North-European mesolithic по таблице эталонов должен достигать 24-37 (в зависимости от местности), обычно же у финнов он в районе 16-18.

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

Мной был проведен эксперимент с бета-версией нового 27-компонентного калькулятора Вадима Веренича, который идеологически заменяет его предыдущий этно-калькулятор World-22. Были собраны данные по восточноевропейцам, не участвовавшим в исходном расчете — в общей сложности использованы результаты 68 человек. В основном это были русские из разных регионов, также собраны результаты эрзя, чувашей, балтов, беларусов, ашкенази, бойко, поляка, финнов, норвежцев, турок.
Закономерности полностью подтвердились. Протестированные со схожим происхождением оказались близки друг к другу по содержанию основных компонентов. Популяции уверенно выделяются и каждый попадает в нужный кластер. Некоторое исключение составляют сложно-смешанные случаи — там результат искажается либо неравномерной передачей генов, либо играют роль какие-то другие эффекты. Несколько хуже других группируются также жители Юга России — очевидно, из-за того, что регион был заселен сравнительно недавно и переселенцы происходили из разных местностей.

Таким образом, «эффект калькулятора» хорошо поддается устранению и в будущем, как я надеюсь, исчезнет из оракулов.

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

В целях более полного освещения проблемы, необходимо изложить предысторию вопроса.  Известный геномный блоггер Давид Polako Веселовский в свой записи на блоге Eurogenes (май 2012 года) указывает на то обстоятельство, что  при использовании геномных калькуляторов и оракулов на основе инструментов Диенека Понтикоса, многие люди получают искаженные результаты, несмотря на то, что все делается по инструкции. Например, пользователи из Великобритании часто оказываются в этих калькуляторах гораздо ближе к выходцам из континентальной Европы, чем  этого следовало бы ожидать. Некоторые из них на самом деле полагали, что эта странная картина обусловлена  тем, что они генетически гораздо более похожи на «норманнов» или англосаксов, нежели средний британец. Polako полагает, что истинная причина кроется в том, что он называл «эффектом калькулятора». Эффект калькулятора  состоит в том, когда результаты  людей, которые входят в состав выборки использованной в анализе Admixture при определении частот аллелей калькулятора значительно отличаются от аналогичных результатов людей, которые не были включены в первоначальную выборку, несмотря на то, что обе группы пользователей имеют одинаковое происхождение, и следовательно в обеих случаях можно ожидать идентичные результаты.

В своей ответной статье-апологии «On the so called «Calculator Effect»» (август 2012 года), автор программы — Диенек Понтикос — приводит сокрушительную аргументацию против негативных и поверхностных замечаний Polako. Во-первых, он убедительно показывает, что хотя эффект отклонения результатов сторонних людей от эталонных результатов и имеет место быть, однако он не имеет никакого отношения к алгоритмам калькулятора и оракула. Поэтому, в корне неверно полагать, что расхождение результатов между описанными выше группами является артефактом алгоритма программа. Из этого вытекает очевидное заключение о неверном описании проблемы как следствия работы программы.  Следовательно, сам термин эффект калькулятора неправилен.

Далее Диенек, c присущей ему убедительной и взвешенной аргументацией, показывает что эта проблема или аномалия была описана им задолго до Polako, в заметке «Further caution on admixture estimates: at the edges of variation», датированной октябрем 2011 года (!).

В ходе анализа причин сильного расхождения результатов выборки армян из статьи Юнусбаева и др. (2011) и результатов армян-участников-проекта и армян из выборки в статье Дорона Бехара (2011),  Понтикос использует наглядную геометрическую визуализацию вариации в виде клинов.

Диенек приходит к очевидным заключениям.

При определении относительного  (в отношению к другим индивидам, в том числе и индивидам из референсных выборок) положения индивида на клиньях (clines)  необходимо помнить, что положение  индивида наиболее точно определяется в тех случаях когда края клина надежно «закреплены». На практике это означает,  что при анализе населения определенного региона необходимо использовать как можно больше групп популяций и отдельных лиц по всему периметру изучаемого региона. Далее, позиция  индивида наиболее точно определяется в тех случаях, когда «клин» длинный, т.е небольшие отклонения из-за дрейфа или неполноты выборки по краям  пренебрежимо малы в сравнении с длиной «клина». Наконец,  компоненты маркирующие значительные меж- и внутриконтинентальные расстояния на (например, европейские компоненты по отношению компонентов из  восточной Азии) оцениваются более точно, чем те,  которые маркируют короткие расстояния (например, южноевропейский и западно-азиатские компоненты). От себя могу добавить к выводам Диенека одно важное наблюдение:  точность результатов Admixture (и основанных на этих результатах входных данных соответствующих калькуляторов и оракулов) сильно зависит от двух факторов — степени (Fst)  дивергенции компонентов и числа снипов, использованных при анализе.  Авторы Admixture отмечают, что для компонентов со значительной степенью относительной дивергенции достаточно 10 000 снипов для выявления генетической структуры выборки, для компонентов с относительно низкой дивергенцией (то есть с короткой дистанцией или, в терминах Диенека, коротким «клином») необходимое число снипов должно быть значительно выше (минимум 100 000 снипов). Значительное влияние на величину оказывает и присутствие в изучаемых выборках так называемых редких вариантов, то есть одиночных нуклеотидных вариантов, с частотой распространения менее 1 процента (смотри статью Estimating and interpreting FST: The impact of rare variants).

 

Решение проблемы.

Как справедливо отмечает Понтикоса, у этой проблемы нет простого решения. Для коротких клинов (например, между популяциями в Европе), которые не надежно закреплены ( так как каждая популяция представлена несколькими людьми из отдаленных групп), можно ожидать относительно большие систематические ошибки. Понтикос считает, что самое доступное в ближайшее время решение — это увеличение числа и качества выборок. Другого обходного пути не существует. Хотя общие детали генетического ландшафта (например,  основные континентальные группы-компоненты генетического разнообразия) уже сейчас легко вывести, но детали всегда можно будет улучшить.

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

 

 

Новая работа по геному неандертальцев и денисовцев.

В июле этого года в журнале Gene появилась замечательная работа Neanderthal and Denisova genetic affinities with contemporary humans:
Introgression versus common ancestral polymorphisms/Robert K. Lowery, Gabriel Uribe , Eric B. Jimenez , Mark A. Weiss, Kristian J. Herrera, 
Maria Regueiro, Rene J. Herrera. Gene . Особого внимания в этой работе заслуживает постановка вопроса в исследовании вопроса о схожести геномов денисовцев и неандертальцев с современными человеческими популяциями.  В этой связи я позволю себе удовольствие процитировать краткий реферат этой статьи в изложении профессионального русскоязычного генетика Людмилы Р.:

Авторы решали вопрос — являются ли те 1-4% генетического сходства между
архаичными гоминидами и современными людьми результатом имевшего место смешивания или общего наследственного полиморфизма, который сохранился в человеческой популяции?
Авторы сравнили 5 млн.SNPs (финальный набор 37,758 SNPs) ныне живущих людей (n=827 из 11 популяций) и архаичных гоминид. Они разделили снипы на 4 группы, которые, логично предположить, происходили в разные отрезки времени –
NdDa –у неандертальцев –derived (мутировавшие) и ancestral (предковые) – у денисовцев,
NaDd — у неандертальцев – ancestral и derived у денисовцев,
NdDd — derived у неандертальцев и денисовцев,
NaDa – ancestral у неандертальцев и денисовцев.
Ancestral и derived снипы определялись по сравнению с шимпанзе.
Ранее предполагалось, что не-Африканские популяции содержат 1-4% генома, доставшегося им от неандертальцев, в отличие от популяций Sub-Saharan-Africans, за счет того, что было смешивание с неандертальцами после выхода человека из Африки. По этому сценарию, все потомки древней человеческой популяции должны содержать равное количество неандертальской ДНК. При этом отличия Sub-Saharan-Africans и non- Sub-Saharan-Africans приписываются gene flow от неандертальцев. Но то, что какой-то SNP у человека, найден у неандертальцев, но не найден у шимпанзе, не означает, что он появился у неандертальца. Такая мутация могла произойти от времени разделения линий шимпанзе и гоминид ( 4-7 млн.лет назад) до времени разделения ветвей человека и неандертальца (400-800 тыс.лет назад). Т.е. выявленные общие SNPs у человека и неандертальцев могут означать их общий предковый полиморфизм.
Авторы не отказываются от “выхода человека из Африки”, но предполагают, что региональные различия в Африке внутри общей предковой популяции были уже на ранних стадиях, и люди, которые мигрировали из Африки, могли представлять собой субпопуляции с более высоким сродством с неандертальцами или денисовцами.
11 популяций, которые участвовали в сравнении:
Abbreviation n Region Populations included
a 123 Sub-Saharan Africa Yoruba, Mandenka, San, Bantu,
Biaka Pygmy, Mbuti Pygmy
b 41 Northern Africa Ethiopians, Egypt, Morocco
c 68 Caucasus Georgia, Armenians, Lezgins, Adygei
e 124 Europe Lithuanian, Belorussian, Romanian,
Cypriot, Hungarian, Basque, Russian,
Spanish, Chuvash
m 33 Melanesian Papuan, Bouganville
n 31 Amerindian Pima, Piapoco, Curripaco, Mayan
s 67 South Central Asia Paniya, Kannadi, Sakilli, Kalash, Uygur,
Barusho, Balochi
r 35 SouthWest Asia Iranian, Uzbekistan
d 30 South East Asia Yizu, Cambodian, Lahu, Malayan
t 34 North East Asia Yakut, Mongolian, Daur
z 241 Near-East Jordan, Samaritan, Syrian, Druze,
Bedouin, Mozabite, Palestinian,
Turkey, Lebanon, Saudi, Yemen
В работе использовали методы популяционной генетики — Principal component (PC) и Structure analyses, D-statistics. Авторы делают выводы, что присутствие 3,6 % неандертальских генов в европейских геномах более похоже на полиморфизм нашего общего предка, чем на результат спаривания видов. % общих генов уменьшается с продвижением на восток в Евразию. Предполагаемая примесь у меланезийцев денисовских генов может также свидетельствовать об их общем предке.

Примечательно, что задолго до публикации этой интересной статьи, к аналогичным выводам пришли любители — антрополог Джон Хоукз (анализ интрогрессии геномов в выборке 1000genomes) и уже ставший живой легендой геномный блогер Диенек Понтикос ( пост о вопросе наличия неандертальский/денисовский адмикса) . Эти выводы противоречат широко растиражированному в масс-медиа выводу о том что «неандертальцы занимались сексом с предками современных людей, за исключением африканцев из региона Суб-Сахары». Этот фривольный медиа-мем возник на основании вольной интерпретации серьезного исследования коллектива под руководством Сванте Паабо. Позже появилось еще одно исследование «The Shaping of Modern Human Immune Systems by Multiregional Admixture with Archaic Human», в котором было показано, что  вклад денисовца в евразийские гены оказался более скромным, однако его доля, как выяснилось, достигает 6% у современных меланезийцев и населения Новой Гвинеи. Соответственно, в средства массовой информация прошла очередная ‘сексуальная’ новость — оказывается, «cпособность успешно противостоять евразийским микробам мы обрели благодаря бракам с неандертальцами и денисовцами». Причем никто из журналистов, похоже не вникал в технические особенности этих работ, в которых ascertainment (установление) снипов производилось по субсахарской популяции бушменов.

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

Существуют два возможных сценария генезис снипов, которые обнаружены у неандертальцев, денисовских людей и современных людей: 1) они возникли у общих предков всех трех групп или 2) снипы могут являться следствием  обмена генами между эти тремя группами.

Результаты D-статистических анализов демонстрируют более высокое присутствие NdDd (derived у неандертальцев и денисовцев) аллелей в Африке к югу от Сахары относительно всех евразийцев и населения Северной Африки. Конечно, сочетание этих двух сценариев может объяснить происхождение подмножества снипов в наборе NdDd подмножество. В модели европейской примеси (адмикса),  у африканцев Субсахары должно быть меньше NdDd аллелей, чем у евразийцев и населения Северной Африки. Исходя из этого, высокий процент NdDd аллелей в субсахарских популяциях является решающим  аргументом в пользу происхождения этих аллелей от древних гоминид, а не в пользу версии смешивания с архаичными людими. Кроме того, к югу от Сахары доля 3 из 5 NdDd компонентов адмикса , включает в себя примерно 30% от общего числа снипов в NdDd  (1 компонент, 6 и 10 на рис. 10), что предполагает общее происхождение предков, а не трехстороннее смешивание для снипов, включенных в панель NdDd . Дополнительным аргументов в пользу сценария общего предкового полиморфизма является расположение снипов NdDd неандертальцев и Денисова NdDd в пространстве первого главного компонента разнообразия PC1 (рис. 5) рядом с субафриканскими популяциями. Этот результат является неожиданным, учитывая, что генотипы NdDd состоят только из деривативных (derived) аллелей. Интересно, что снипы NdDd демонстрируют сопоставимые характеристики D-статистики  в популяциях меланезийцев и африканцев Субсахары. Обе группы — меланезийцы и субсахарские африканцы —  генерируют более высокие показатели D-статистики на основании снипов NdDd основана D-статистики, — примерно на 0,7% выше, чем у выходцев из Северной Африки и на 1,5% выше, чем у жителей Северо-Восточной Азии (табл. 2;. Рис 11). Кроме того, если мы предположим, что смешивание между тремя группами (предками современных людей, неандертальцами и денисовцами) является важным источником снипов NdDd , то можно  было бы ожидать более высокую степень генетического обмена между европейцами и меланезийцами, однако, во всех наших анализах Structure, ни один такой обмен не наблюдается.

Эксперимент

С целью проверки выводов этого исследования я провел дополнительное исследование этого вопроса с использованием большого количества современных популяций (более трехста популяций), большего числа снипов (примерно 300 000 снипов) и альтернативных методов — программы Admixture и MDS (мультидименсионального скалирования).

Результаты нашего анализа замечательно согласуются с результатами в обсуждаемой статье. Например, на графике PCA (в пространстве двух первых главных компонентов генетического разнообразия), неандертальцы и денисовцы оказываются в окружении трех субсахарских популяций — гумуз, ари из Восточной Африки и лемба из южной Африки. 

Результаты Admixture также мало чем отличаются от результатов аналогичного анализа Structure в работе профессиональных попгенетиков. Снипы неандертальцев и денисовского человека (взятые из кураторского набора данных лаборатории Райха (SNP ascertainment panel)) образуют особый компонент вместе со снипами субсахарских популяций бушменов, коса и сандаве.

Denisova	Denisova	99,98%
Neander	Vindija	99,98%
San	HGDP00991	99,98%
San	HGDP01032	99,98%
San	HGDP01036	99,98%
San	SA36	99,98%
San	SA34	99,98%
San	SA52	99,98%
San	SA19	99,98%
San	HGDP00988	99,54%
San	HGDP01029	99,19%
San	HGDP00992	98,47%
San	SA53	97,53%
San	SA47	93,98%
San	SA41	93,28%
San	SA22	92,13%
San	SA32	91,10%
Neander	Clint	90,75%
San	SA48	89,58%
San	SA30	89,40%
San	SA55	88,93%
San	SA35	88,18%
San	SA61	85,45%
San	SA50	83,92%
San	SA29	81,14%
San	SA04	78,20%
San	SA37	74,40%
San	SA56	74,34%
San	SA38	74,17%
San	SA21	70,00%
San	SA06	69,85%
San	SA28	61,13%
San	SA03	57,39%
San	SA40	56,62%
San	SA49	54,89%
San	SA45	47,39%
San	SA58	43,01%
San	SA39	41,33%
San	SA59	34,80%
Bantu	HGDP01030	33,37%
Xhosa	XH4	26,85%
Xhosa	XH20	25,99%
Xhosa	XH14	24,78%
Bantu	HGDP00993	23,99%
Bantu	HGDP00994	23,02%
Bantu	HGDP01034	21,48%
San	SA25	21,28%
Bantu	HGDP01033	15,40%
Sandawe	HG43	14,20%
Sandawe	HG60	14,04%
Sandawe	HG40	13,77%
Sandawe	HG35	13,56%
Sandawe	HG44	13,51%
Sandawe	HG56	13,37%
Sandawe	HG46	13,26%
Sandawe	HG41	13,25%
Sandawe	HG66	13,18%
Sandawe	HG47	13,09%
Sandawe	HG49	12,93%
Sandawe	HG67	12,75%
Sandawe	HG55	12,63%
Sandawe	HG45	12,43%
Sandawe	HG63	12,14%
Aricultivator	Aricultivator11	12,13%
Ariblacksmith	Ariblacksmith2	12,13%
Sandawe	HG42	12,10%
Ariblacksmith	Ariblacksmith3	11,92%
Sandawe	HG38	11,85%
Ariblacksmith	Ariblacksmith7	11,83%
Sandawe	HG53	11,76%
Ariblacksmith	Ariblacksmith6	11,70%
Aricultivator	Aricultivator2	11,67%
AricultivatorIbd	Aricultivator23Ibd	11,54%
Ariblacksmith	Ariblacksmith10	11,49%
Ariblacksmith	Ariblacksmith8	11,48%
Aricultivator	Aricultivator17	11,46%
Aricultivator	Aricultivator4	11,42%
AricultivatorIbd	Aricultivator24Ibd	11,28%
Sandawe	HG48	11,22%
Aricultivator	Aricultivator15	11,18%

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

Neander	Clint	100,00%
Denisova	Denisova	100,00%
Neander	Vindija	100,00%
Papuan	HGDP00542	100,0000%
Papuan	HGDP00554	100,0000%
NAN_Melanesian	HGDP00662	100,0000%
NAN_Melanesian	HGDP01027	100,0000%
Papuan	HGDP00543	100,0000%
Papuan	HGDP00555	100,0000%
NAN_Melanesian	HGDP00663	100,0000%
Papuan	HGDP00544	100,0000%
Papuan	HGDP00556	100,0000%
NAN_Melanesian	HGDP00664	100,0000%
Papuan	HGDP00545	100,0000%
NAN_Melanesian	HGDP00490	100,0000%
NAN_Melanesian	HGDP00787	100,0000%
Papuan	HGDP00546	100,0000%
NAN_Melanesian	HGDP00491	100,0000%
NAN_Melanesian	HGDP00788	100,0000%
Papuan	HGDP00547	100,0000%
NAN_Melanesian	HGDP00655	100,0000%
NAN_Melanesian	HGDP00789	100,0000%
Papuan	HGDP00548	100,0000%
NAN_Melanesian	HGDP00656	100,0000%
NAN_Melanesian	HGDP00823	100,0000%
Pima	HGDP01048	100,0000%
Papuan	HGDP00541	100,0000%
Papuan	HGDP00553	100,0000%
NAN_Melanesian	HGDP00661	100,0000%
NAN_Melanesian	HGDP00979	100,0000%
Karitiana	HGDP00998	100,0000%
Karitiana	HGDP01011	100,0000%
Surui	HGDP00833	100,0000%
Surui	HGDP00846	100,0000%
Karitiana	HGDP01010	100,0000%
Surui	HGDP00832	100,0000%
Surui	HGDP00845	100,0000%
Papuan	HGDP00550	100,0000%
NAN_Melanesian	HGDP00658	100,0000%
NAN_Melanesian	HGDP00825	100,0000%
Karitiana	HGDP00999	100,0000%
Karitiana	HGDP01012	100,0000%
Surui	HGDP00834	100,0000%
Surui	HGDP00847	100,0000%
Papuan	HGDP00540	100,0000%
Papuan	HGDP00552	100,0000%
NAN_Melanesian	HGDP00978	100,0000%
Karitiana	HGDP01000	100,0000%
Karitiana	HGDP01013	100,0000%
Surui	HGDP00835	100,0000%
Surui	HGDP00848	100,0000%
Karitiana	HGDP01001	100,0000%
Karitiana	HGDP01014	100,0000%
Surui	HGDP00837	100,0000%
Surui	HGDP00849	100,0000%
Karitiana	HGDP01003	100,0000%
Karitiana	HGDP01015	100,0000%
Surui	HGDP00838	100,0000%
Surui	HGDP00850	100,0000%
Karitiana	HGDP01004	100,0000%
Karitiana	HGDP01016	100,0000%
Surui	HGDP00839	100,0000%
Surui	HGDP00851	100,0000%
Surui	HGDP00843	100,0000%
Pima	HGDP01050	100,0000%
Karitiana	HGDP01005	100,0000%
Karitiana	HGDP01017	100,0000%
Surui	HGDP00840	100,0000%
Surui	HGDP00852	100,0000%
Karitiana	HGDP01006	100,0000%
Karitiana	HGDP01018	100,0000%
Surui	HGDP00841	100,0000%
Karitiana	HGDP00996	100,0000%
Karitiana	HGDP01008	100,0000%
Papuan	HGDP00551	100,0000%
Karitiana	HGDP00995	100,0000%
Karitiana	HGDP01007	100,0000%
Karitiana	HGDP01019	100,0000%
Surui	HGDP00842	100,0000%
Karitiana	HGDP00997	100,0000%
Surui	HGDP00830	100,0000%
Surui	HGDP00844	100,0000%
NAN_Melanesian	HGDP00657	100,0000%
NAN_Melanesian	HGDP00824	100,0000%
Papuan	HGDP00549	100,0000%
Hadza	BAR01	100,0000%
Hadza	BAR04	100,0000%
Hadza	BAR07	100,0000%
Hadza	BAR08	100,0000%
Hadza	BAR10	100,0000%
Hadza	BAR11	100,0000%
Hadza	BAR13	100,0000%
Hadza	END08	100,0000%
Hadza	END09	100,0000%

Перспективы изучения линкаджа в плане определения генеалогической наследственности в изолированных популяциях (заметки доктора К.Булаевой)

По просьбе уважаемой К.Булаевой, я произвел анализ линкаджа в отдельном регионе 6 хромосомы в одной популяционной выборке (какая именно это была выборка, я точно не могу сказать).

Kazima Bulayeva:

Привет Вадим, LD? а admixture ? Мы же как договорились -результаты вместе смотрим ваши и мои -решаем их совместную публикацию. О моих линкадже я говорила. Что LD показало? По идее более узкий регион? Но этот метод-ассоциативный, а у меня нет выборки здоровых….не соображу что может дать этот метод нам. Расскажите плиз что получилось и далее детально обсудим, идет?

Vadim Viarenič-Stachowski: Просто Вы ничего не говорили про admixture

Vadim Viarenič-Stachowski: А сейчас я обработал Ваши данные в программе Haploview.

Kazima Bulayeva: Вадим, я же почти не знаю этот метод. То, что я знаю -это когда изучают в популяциях -можно определить степень геномной гетерогенности в популяции и даже у каждого члена.

Vadim Viarenič-Stachowski: Ее отличие — она позволяет показать блоки с высоким сцеплением наглядно.
http://www.broadinstitute.org/scientific-community/science/programs/medical-and-population-genetics/haploview/screenshots-0

Vadim Viarenič-Stachowski: То есть выявить блоки LD или гаплоблоки

Kazima Bulayeva: Я думаю как раз сейчас -что может дать LD по снипам в хр 6 в дополнении к линкадже? Прежде всего, линкадже я делала на основе STR сканированных по 10 сМ по всему геному каждого…но как понимаете -это too spread. LD может уловить тоньше локус такого сцепления…единственно —как нам сравнить с нормой? Может быть HapMap для контроля?

Vadim Viarenič-Stachowski: Ok, но для вычисления геномной гетерогенности нужны GWAS-данные. Одной хромосомы маловато будет.

Kazima Bulayeva: по популяциям? Да, согласна

Vadim Viarenič-Stachowski: Так Вам нужны результаты анализа в Haploview?

Kazima Bulayeva: Блок 1 -какие снипы включает?

Vadim Viarenič-Stachowski: Я напримре видел такие вот треугольные плоты в презентациях Степанова

Vadim Viarenič-Stachowski: В графике все подписано

Vadim Viarenič-Stachowski: с обозначением снипа в rs-формате.

Kazima Bulayeva: снипы какие-то другие названия…rs….по идее должны быть ?

Vadim Viarenič-Stachowski: Так ведь это не мой график )), а в качестве примера с сайта программы на Broad Instutute )

Vadim Viarenič-Stachowski: Я справшивал другое — Вам нужны графики такого формата?

Kazima Bulayeva: Нет. Давайте сформулируем задачу: у нас есть данные из 4-этнически разных изолятов с высоким сцеплением с SCZ в 6p21. В сцепленном регионе локализовано около 25 генов…много генов-большой отрезок генома -около 10 сМ т.к STRs/ Что позволят определить снипы? Не все эже 25 генов связаны -а какие-то 1-2 гена из общего числа. Поэтому снипы и LD могут помочь выявить из числа 25 те именно гены которые сцеплены с заболеванием. Согласен с задачей?

Kazima Bulayeva: мне кажется логично поставленный вопрос. и LD вполне адекватный инструмент даже без контроля, т.к. мы его используем как 2-й этап углубления в мезанизм установленного в родословных сцепления

Kazima Bulayeva: permutation p -недостоверен нигде?

Kazima Bulayeva: Вадим, далее: если ы ЛД мы установили внутри сцепленного с STRs региона блоки снипов у больных-мы можем проверить функциональную роль снипов-типа интрон или экзон и в каких генах…т.е. выявляем конкретные гены и геномнын механизмы

Vadim Viarenič-Stachowski: Разумно.

Vadim Viarenič-Stachowski: Хорошо, я перешлю Вам выходные данные из своего анализа, а потом подумаем каких применить и что ценного можно извлечь.

Kazima Bulayeva: статическая достоверность есть у блоков ? Всего 3 блока выявляются? и наверняка мы сможем определить их цитобанды и гены в блоках? Еще-я сделал CNV и LOH в этой же хромосоме. Мы сможем посмотреть эти блоки в LD на предмет указанный аберраций

Kazima Bulayeva: по-моему должно что-то быть выявлено интересное с добавлением LD по снипам—

Vadim Viarenič-Stachowski: Вот и ладненько. Сегодня или завтра перешлю.

Kazima Bulayeva: Вадим, а вы в Stanley Center работаете?

Kazima Bulayeva: там по писихиатрической генетике работы давно проводятся….Не смогли бы узнать-есть ли у них возможность для типирования снипов? у меня есть ДНК из родословных с психопатологией и с STRs

За кулисами: как создавался этно-популяционный калькулятор World-22

Летом 2011 года я создал целый рядсобственных модификаций получившего широкую известность калькулятора DIY Dodecad гениального грека Диенека Понтикоса. К моему приятному удивлению, за прошедшее время калькулятором успело воспользоваться несколько тысяч людей, некоторые из которых даже выложили свои результаты в Интернете.  Разумеется, многие также разместили и свои собственные интерпретации полученных результатов. Некоторые из приведенных в комментариях интерпретации выделялись (в хорошем смысле этого слова) высоким академическим уровнем, но мне попадались и такие комментарии, при чтении которых становилось понятно, что авторы не только не понимают принципов и сути парадигмы анализа, предложенного Понтикосом, но и — что гораздо хуже — выдавали свои фантазии за действительности. Особенно часто мне попадались подобные фантастические рассуждения в русскоязычном секторе Интернета.Пример такого невежества можно найти в рассуждениях само-провозглашенного академика ДНК-генеалогии Анатолия  Клесова:

Но и в этом случае различия все равно будут между русскими и монголами. Качественно и как-то полуколичественно его можно рассматривать, но не в виде профанации, как это делает Понтикос. Более того, это рассмотрение – если правильно – надо проводить не на выбранных маленьких фрагментах, а действительно по всему геному. На маленьких фрагментах будут вылезать отдельные особенности – то присущие в основном, например, гаплогруппам Y-I2 и мтДНК-Н, то кому-то еще. И это еще будет зависеть от разрешения, которые и обозначают индексами К=4, К=8 и другими. То есть берут маленький фрагмент генома, да еще с малым (или бóльшим) разрешением, стягивают в точку, и все равно получают в целом ерунду. Но для коммерции годится. Годятся для коммерции и вот такие, в частности, «открытия» того же Понтикоса: Перевод: Интересно то, что европейская популяция показывает присутствие американских индейцев, что показывает и f-статистика, и она же показывает присутствие компонента с Сардинией. Как видим, Понтикос уже забыл, что названия им придуманы как попало, и уже придает им абсолютные значения. Про Сардинию Понтикос уже вошел в состояние экзальтации. Он придает Сардинии некую пра-европейскую значимость, на основании, конечно, этой ерунды с «геномом», который анализирует как хочет. Пример – он трубил по всему свету, что Отци, «ледовый человек», имел геном «Сардинии». Однако только что опубликована статья о том, что Отци – никакая не Сардиния, а типичная Центральная Европа. Ну, и что делать будем? Понтикос, с его страстным желанием сенсаций, каждый раз наступает на одни и те же грабли. Впрочем, фарс продолжается. Теперь тем же занялся некто российский Веренич, а именно тоже насчитывает «польскую компоненту», пользуясь подходом своего гуру-Понтикоса.

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

Предварительные замечания

Как вы возможно знаете, MDLP блог не обновлялся с февраля 2012 года.  Полгода тому назад я пообещал себе, что я не буду писать новые сообщения на MDLP блоге до те пор пока я не напишу краткую научный отчет о проделенной работе. Так как приоритеты завершения научной работы были важнее рутиного обновления блога,  то  в связи с нехваткой времени, я был не в состоянии продолжать обновление блога на регулярной основе, в связи с нехваткой времени, я должен был внести изменения в свой исследовательский график. Поэтому я решил воздерживался от размещения новых данных на блоге в течение нескольких месяцев, фокусируясь на более важных вопросах. Несмотря на все ограничения, я продолжал втайне работать  на проектом MDLP, сбором необходимых данных и выполением различных ‘геномных’ экспериментов в целях достижения своей конечной цели. Однако с течением времени, некоторые результаты секретных экспериментов с новыми полногеномными популяционными выборками и инструментами в конечном итоге просочились в Интернет,  порождая огромный интерес к моему проекту. После выпуска новой версии моей собственной модификации DIYDodecad калькулятор на сайте Gedmatch.com, я был буквально завален письмами пользователями сервиса Gedmatch.com.

Тогда я осознал свою основную стратегическую ошибку, которая заключалась в  отсутствии подробной документации к выпущенными мной данными и результатам анализа, и почувствовал себя обязанным разместить более подробные разъяснения. Очевидно, я начну новую серию публикацию в своем блоге,  которая будет тесным образом связанна с теми аспектами моей работы, которая наиболее интересует общественность, то есть с калькулятором MDLP World22.

Основы отбора референсных популяций калькулятора MDLP World22.

Референсный набор  популяций в этом калькуляторе был собран в программе PLINK   методом «intersection&thinning» ( дословно «пересечением и истончением») образцов из различных источников данных: HapMap 3 (отфильтрованный набор данных КЕС, YRI, JPT, CHB), 1000genomes,   Rasmussen et al. (2010),   HGDP (кураторская база данных Стэнфордского университета), Metspalu et al. (2011),  Yunusbayev et al (2011), Chaubey et al. (2010) и т.д.

Кроме того, я отобрал произвольным образом по 10 сэмплов (или максимальное количество доступных сэмплов в тех случаях, когда общее число сэмплов в популяции было меньше 10) от каждой европейской страны, представленной в панеле базе данных POPRES. Наконец, для того чтобы оценить степень корреляции между современным и древним генетическим разнообразием населения Европы, я также включил в выборку образцы древней ДНК Эци (Keller et al. (2012)) ,  образцы житлей шведского неолита Gök4, Ajv52, Ajv70, Ire8, STE7 ( Skoglund et al. (2012)) и 2 образца La Braña  — останков мезолитических жителей Пиренейского полуострова (Sánchez-Quinto et al.(2012)).

Затем я добавил 90 образцов — анонимизированных данных — участников моего проекта. После слияния вышеупомянутых наборов данных и истончения набора SNP с  помощью особой команды PLINK, я исключил SNP-ы с  более чем 1% минорных аллелей. После чего я отфильтровал дубликаты, лиц с высоким уровнем общих по происхождению идентичных сегментов (IBD). В качестве критерия фильтрации  были использованы расчеты IBD в Plink, где IBD представлена как средняя доля аллелей общих между двумя людьми по всем анализируемым локусам.  Затем я удалил из выборки лиц с  высоким коэффициентом предпологаемого родства (коэффициенты родства были вычислены в программном обеспечении King). Для получения более стабильных результатов, я также отфильтровал сэмплы с более чем 3 стандартными отклонениями от средних данных  по популяции. Поскольку коэффициент родства может быть надежно определен с помощью оценки HWE (ожидания, вытекающего из закона Харди-Вайнберга) между SNP-ами с той же базовой частотой аллелей, то SNP-ы с существенным отклонением (p < 5.5 x10⁻⁸) от  ожидания Харди-Вайнберга были удалены из объединенного набора данных. После этого я выделил те SNP-ы, которые присутствовали в чипах Illumina / Affymetrix, и затем произвел фильтрацию снипов на основе расчетов степени неравновесного сцепления  (в этой я использовал хромосомное ‘окно’ размером в 50 базовых пар, с шагом 5 базовых пар и пороговым значением уровня сцепления R ^ 2, равным 0,3).

По окончанию этой сложной последовательности операций, я получил окончательноый набора данных, который включал в себя 80 751 снипов,  2516 человек и 225  референсных популяций.

Анализ этно-популяционного адмикс

 В ходе следующенго этапа, окончательный набор данных по референсным популяциям (которые я храню в linkage-формате PLINK) был обработан в программе Admixture.  Во время выбора подходящей модели проведения теста на этно-популяционный адмикс, я столкнулся с крайне трудной задачей: как было показано в профильных научных исследованиях (Patterson et al.2006) количество маркеров, необходимых для надежной стратификации популяций в анализе обратно пропорциональна генетическому расстоянию (фСТ) между популяциями. Согласно рекомендациям пользователей программы Admixture, считается что примерно 10 000 генетических SNP-маркеров достаточно для выполнения интер-континентальной GWAS -коррекции обособленных популяций (например, уровень дивергенции между африканскими, азиатскими и европейскими популяциями  FST > 0.05), в то время как для аналогичной коррекции между внутриконтинентальными популяциями требуется более чем 100000 маркеров (в Европе, например, ФСТ < 0.01). Для повышения точности результатов Admixture я решил использовать метод, предложенный Dienekes. Этот метод позволяет преобразовать частот аллелей в «синтетические» индивиды (см. также пример Зака Аджмала из проекта HarappaDNA). Идея метода довольно проста: сначала необходимо запустить unsupervised анализ Admixture с целью вычисления частот аллелей в так называемых предковых компонентов, а затем на основании аллельных частот сгенерировать «фиктивные популяции».  Именно эти фиктивные популяции и индивиды будут использоваться в ходе чистых референсов в ходе последующего анализа этно-популяционного анализа.
Впрочем, как и любой другой исследователь, работающий над четким решением проблемы этно-популяционного адмикса, я вынужден считаться с ограничениями этого подхода. Хотя я и отдаю себе отчет в существовании явных методологических подвохов в использовании смоделированных искусственных индивидов для определения адмикса в реальной популяции, я все же скорее склонен согласиться с Понтикосом, которые считаeт полученных в ходе  аллельно-частотного моделирования «фиктивных индивидов» лучшей аппроксимацией древних генетических компонентов мирового народонаселения.Как бы то не было, моделирующий подход, предложенный Диенеком и Заком, сослужил свою хорошую службу, поскольку были мной были получены  значимые результаты в ходе создания нового калькулятора. Сначала я произвел unsupervised Admixture (при значении К = 22, т.е 22 кластера частот аллель или предковых компонентов). По выполнению анализа нами были получены оценки коэффициентов адмикса в каждой из этих 22 аллельных кластеров, а также частоты аллелей для всех SNP-ов в каждой из 22 родовых популяций.

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

Pygmy

West-Asian
North-European-Mesolithic
Tibetan
Mesomerican
Arctic-Amerind
South-America_Amerind
Indian
North-Siberean
Atlantic_Mediterranean_Neolithic
Samoedic
Proto-Indo-Iranian
East-Siberean
North-East-European
South-African
North-Amerind
Sub-Saharian
East-South-Asian
Near_East
Melanesian
Paleo-Siberean
Austronesian

Вышеупомянутые частоты аллель, вычисленные в ходе unsupervised (безнадзорного) анализа (Admixture K = 22) объединенного набора данных, были затем использованы для симуляции синтетических индивидов, по 10 индивидов на каждую из 22 предковых компонент.  Это симуляционное моделирование проводилось с помощью PLINK команды -simulate Когда моделирование было закончено, я сделал визуализацию расстояния между симулированными индивидами с использованием многомерного масштабирования.

На следущем этапе, я включил группу смоделированных индивидов (220 индивидов) в новую эталонную популяцию. После чего я запустил новый анализ А, на этот раз в полном «поднадзорном» режиме для K = 22, причем полученные в ходе симуляционного моделирования фиктивные популяции фиктивных индивидов использовались в качестве новых референсных эталонных групп.  На конвергенцию 22 априорно заданых предковых компонентов было затрачено  31 итераций (3 7773,1 сек) с окончательным loglikelihood: -188032005,430318 (ниже приведена таблица значений Fst  между расчетными ‘предковыми’ популяциями):

Приведенная выше матрица  Fst дистанций  была использована для определения наиболее вероятной топологии NJ-дерева всех 22 предковых компонентов ( примечание: в качестве outgroup-таксона использовался South-African component). Индивидуальные результаты ‘поднадзорного’ анализа этно-популяционных миксов (в формате Excel) для участников проекта были загружены на GoogleDrive.

MDLP World22 DIYcalculator

Выходные файлы «поднадзорного» анализа  Admixture K=22 (средние значения коэффициентов адмикса в референсных популяциях и значения Fst) были использованы для разработки новой версии DIYcalculator MDLP, который более известен под кодовым названием «World22» (онлайн версия доступна разделе Admixture-утилит на сервисе Gedmatch в рамках проекта MDLP). Как я уже упоминал выше, MDLP DIYcalculator работает на коде Dodecad DIY calculator (c) Dienekes Pontikos.

В свою очередь,  реализованная на сервисе  Gedmatch модификация DIYcalculator ‘World22’ комбинирована с  Oracle ‘World22’ MDLP, который также работает на коде Диенека и Зака Аджмала ​​(Хараппа/DodecadOracle). Программа «Oracle» работает в двух режимах. В режиме single population программа определяет ближайщие (к анализируемому геному) референсные популяции калькулятора Word22. В смешанном режиме, Oracle рассматривает все пары населения, и для каждой из пар вычисляет минимальное Fst-взвешенное расстояние между парой и анализируемым геномом, а также  коэффициенты сходства.

Предковые популяции (т.е. полученные в ходе симуляционное моделирования популяции — см. выше) обозначены в результатах Oracle суффиксом anc, в то время реальные современные и древние популяции обозначены суффиксом der.

Если у Вас возникли проблемы с пониманием/интерпретацией результатов Oracle и DIYcalculcator,  то я настоятельно рекомендую обратится к соответствующим темам в блогах  Dodecad и НаrappaWorld . Я полагаю, что не имеет особого практического смысла заново изобретать велосипед и слово в слово повторять то, что уже было написано более компетентными в этом вопросе людьми.

Что представляют собой компоненты MDLP World-22?

Один из наиболее частых вопросов, которые задают мне пользователи калькулятора, напрямую касается практической интерпретации референсных популяций и предковых компонентов в моих калькуляторах K = 12 и World-22 анализов в виду. Чуть выше по тексту я уже привел часть ответа на этот вопрос , но — как гласит старинная китайская пословица — одна картинка стоит десять тысяч слов. Вот почему я решил визуализировать компоненты на поверхности земного шара путем отображения коэффициентов адмикса. Избегая излишних премудростей, я воспользовался готовым рецептом Франсуа Оливье, который предложал  использовать графическую библиотеку статистического программного обеспечения R для отображения пространственной интерполяции  коэффициентов адмикса (Q матрица) в двух измерениях (где пространственные координаты записываются как географические долгота и широта).  Благодаря этому решению, мне удалось создать по 2 контурные карты на каждый из предковых компонентов.Pygmy (модальный компонент в  популяциях африканских пигмеев Biaka и Mbuti)

West-Asian (бимодальный компонет с пиком на Кавказе и юго-восточной части Ирана, приблизительно идентичен компонентам Caucasian/Gedrosia Диенека Понтикоса)

North-European-Mesolithic (локальный архаичный компонент с пиком в популяции древних европейских жителей Иберийского полуострова La_Brana и современной популяции саамов).

Tibetan (Indo-Burmese) component (Гималаи-Тибет)

Mesomerican (главный генетический компонент  у мезоамериканских америндов)

 

North-Amerind (нативный компонент северо-американских америндов)

South-Amerind (нативный компонент южно-американских индейцев)

  Atlantic-Mediterranean-Neolithic (доминируюший компонент  в западной и юго-западной Европе)

Контурные карты прочих компонентов можно скачать здесь.

Практические рекомендации по работе с данными древней ДНК

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

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

В качестве примера в нашем туториале мы будем использовать данные, любезно предоставленные авторами работы P Skoglund, H Malmström, M Raghavan, J Storå, P Hall, E Willerslev, MTP Gilbert, A Götherström* & M Jakobsson* (2012) Origins and genetic legacy of Neolithic farmers and hunter-gatherers in Europe, и данные работы  Federico Sánchez-Quinto, Hannes Schroeder, Oscar Ramirez, María C. Ávila-Arcos, Marc Pybus, Iñigo Olalde, Amhed M.V. Velazquez, María Encina Prada Marcos, Julio Manuel Vidal Encinas, Jaume Bertranpetit, Ludovic Orlando, M. Thomas P. Gilbert, Carles Lalueza-Fox Genomic Affinities of Two 7,000-Year-Old Iberian Hunter-Gatherers.

Для успешного прохождения туториала нам потребуется:
1) наличие мотивации и желание изучить основы практической геномики

2) посколько большинство инструментов задействованных в данном туториале написны под Unix, то необходимо наличие опыта работы с Unix shell: желательно также иметь доступ к значительным вычислительным мощностями (некоторые из операций описанных ниже я производил в вычислительном кластере Тартуского университета).

3) пакет samtools последней версии

4) пакет snpEFF/snpSift

5) пакет vcftools и программа Plink

6) FASTA-файл с человеческий референсным геномом в  версии билда hg18: я рекомендую использовать модифицированную версию файла, в котором старый референсный митосиквенс заменен на новый референсный митосиквенс (rCRS:NC_012920 gi:251831106).

7) каталог генетических полиморфизмов dbSNP в версии билда hg18.

I этап — bamtools.

Используемые в нашем туториале исходные файлы представлены в формате bam — бинарном варианте стандартного файла SAM используемого для хранения элайнментов сиквенсов.
В нашем случае исходные файлы представляют собой конечный продукт,  в котором уже удалены дупликаты и артифакты клонирования в ходе PCR. Поэтому мы можем сразу же приступить к следущему этапу — объединения файлов bam в один общий файл:

samtools merge AjvIre.bam SNPs_Ajv52_r1_hits_rmdup.bam SNPs_Ajv52_r2_hits_rmdup.bam SNPs_Ajv70_r1_hits_rmdup.bam SNPs_Ajv70_r2_hits_rmdup.bam SNPs_Ire8_r1_hits_rmdup.bam SNPs_Ire8_r2_hits_rmdup.bam  
samtools merge GokSte.bam SNPs_Ste7_r1_hits_rmdup.bam SNPs_Ste7_r2_hits_rmdup.bam SNPs_Ste7_r3_hits_rmdup.bam SNPs_Ste7_r4_hits_rmdup.bam\

Далее, мы провидем сортировку файлов по контигам (контиг — это набор упорядоченных перекрывающихся клонов ДНК, охватывающих всю хромосому ил и какой-либо ее участок):

 samtools sort AjvIre.bam AjvIre.sorted.bam
 samtools sort GokSte.bam GokSte.sorted.bam
 samtools sort BRA.bam BRA.sorted.bam

Скачиваем референсный файл билда hg18

wget http://hgdownload.cse.ucsc.edu/goldenPath/hg18/bigZips/hg18.2bit
wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/twoBitToFa
./twoBitToFa hg18.2bit hg18.fa

Производим индексацию референсного файла человеческого генома (билд hg18) и сравниваем систему обозначения хромосомных контигов с аналогичной системой в наших образцах древних геномов:

samtools faidx hg18.fa
chr10    135374737    7    50    51
chr10_random    113275    138082253    50    51
chr11    134452384    138197801    50    51
chr11_random    215294    275339247    50    51
chr12    132349534    275558854    50    51
chr13    114142980    410555386    50    51
chr13_random    186858    526981240    50    51
chr14    106368585    527171843    50    51
chr15    100338915    635667807    50    51
chr15_random    784346    738013515    50    51
chr16    88827254    738813555    50    51
chr16_random    105485    829417369    50    51
chr17    78774742    829524971    50    51
chr17_random    2617613    909875222    50    51
chr18    76117153    912545195    50    51
chr18_random    4262    990184706    50    51
chr19    63811651    990189061    50    51
chr19_random    301858    1055276960    50    51
chr1    247249719    1055584862    50    51
chr1_random    1663265    1307779589    50    51
chr20    62435964    1309476127    50    51
chr21    46944323    1373160818    50    51
chr21_random    1679693    1421044042    50    51
chr22    49691432    1422757336    50    51
chr22_random    257318    1473442611    50    51
chr22_h2_hap1    63661    1473705091    50    51
chr2    242951149    1473770032    50    51
chr2_random    185571    1721580217    50    51
chr3    199501827    1721769506    50    51
chr3_random    749256    1925261383    50    51
chr4    191273063    1926025631    50    51
chr4_random    842648    2121124169    50    51
chr5    180857866    2121983676    50    51
chr5_random    143687    2306458713    50    51
chr5_h2_hap1    1794870    2306605288    50    51
chr6    170899992    2308436062    50    51
chr6_random    1875562    2482754067    50    51
chr6_cox_hap1    4731698    2484667156    50    51
chr6_qbl_hap2    4565931    2489493503    50    51
chr7    158821424    2494150759    50    51
chr7_random    549659    2656148625    50    51
chr8    146274826    2656709284    50    51
chr8_random    943810    2805909620    50    51
chr9    140273252    2806872313    50    51
chr9_random    1146434    2949951044    50    51
chrM    16571    2951120413    50    51
chrX    154913754    2951137322    50    51
chrX_random    1719168    3109149365    50    51
chrY    57772954    3110902923    50    51

samtools view -H AjvIre.sorted.bam

@HD    VN:1.0    SO:unsorted@PG    ID:dvtgmlqtca    PN:stampy    VN:1.0.10_(r854)    CL:-g hg18 -h hg18 --solexa --sensitive -f sam -o output/stampy_Ajv52_r1_aln1.sam -M /bubo/proj/b2010050/private/seqdata/neolitikum/Neolitisar/pruned/Ajv52_r1_trimmed.txt@CO    TM:Tue, 30 Nov 2010 12:35:43 CET    WD:/bubo/proj/b2010050/private/program/stampy-1.0.10    HN:q207.uppmax.uu.se    UN:pontuss@SQ    SN:NC_000001.9    LN:247249719    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000002.10    LN:242951149    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000003.10    LN:199501827    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000004.10    LN:191273063    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000005.8    LN:180857866    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000006.10    LN:170899992    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000007.12    LN:158821424    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000008.9    LN:146274826    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000009.10    LN:140273252    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000010.9    LN:135374737    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000011.8    LN:134452384    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000012.10    LN:132349534    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000013.9    LN:114142980    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000014.7    LN:106368585    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000015.8    LN:100338915    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000016.8    LN:88827254    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000017.9    LN:78774742    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000018.8    LN:76117153    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000019.8    LN:63811651    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000020.9    LN:62435964    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000021.7    LN:46944323    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000022.9    LN:49691432    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000023.9    LN:154913754    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000024.8    LN:57772954    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_012920.1    LN:16560    AS:hg18_ncbi36_rCRS    SP:human

Итак,  при сравнении вышеупомянутых двух файлов, мы видим что обозначение в референсном генома  отличается от обозначения начала хромосом в файле AjvIre (вместо традиционного обозначения chr1…chrM в этом файле используется номер сиквенса хромосомы в Генбанке, например сиквенс первой хромосомы — NC_000001.9, и т.д.).
Эта проблема решается сравнительно легко с помощью редактирования заголовка bam файла (заменой SO:unsorted@PG  на SO:sorted@PG и номеров Генбанка на порядковый номер хромосом) и следущих комбинаций директив samtools:

samtools view -H AjvIre.sorted.bam > originalheader
gedit originalheader
samtools reheader newheader AjvIre.reh.sorted.bam

Аналогичные операции производим и с файлом GokSte.sorted.bam. Файл Bra.sorted.bam редактировать нет надобности, поскольку обозначения хромосом соответствуют обозначению хромосом в референсном файле.

Таким образом, после выполнения означенных выше операций, мы подошли к самой важной процедуре — snp and indel calling, то есть «вызову» (определению) снипов и инделов в наших отсортированных и модифицированных bam файлах.

Нужно сразу отметить, что процедура нахождения генетических вариантов в древней ДНК существенно отличается от аналогичной процедуры в случае с современной ДНК.  Поэтому приходится применять фильтры samtools, которые в большинстве рутинных анализов просто не используются. Я не буду объяснять, что означает каждый из используемых фильтров. Достаточно будет сказать, что я следую рекомендациям профессора Понтуса Скоглунда.  Принимая во внимание ресурсоемкость операции нахождения генетических вариантов, я задействовал возможности тартуского вычислительного центра (ниже приведен пример с BRA.srt.bam):

qsub runSamstools.sh

#!/bin/bash
# This file is runSamtools
#
#PBS -N Samtools
#PBS -m be
#PBS -k oe
#PBS -l walltime=01:30:00
#PBS -l nodes=4:ppn=8
#PBS -l vmem=4gb
#PBS -d /storage/hpchome/vadim78

cd /storage/hpchome/vadim78/conversion/ancient
module load storage_software
samtools mpileup BRA.srt.bam -q 30 -Q 15 -uf hg18.fa |
/storage/hpchome/vadim78/samtools/bcftools/bcftools view -vcg - > BRA.vcf

II. Аннотация VCF файлов — snpSift.

Итак,  мы получили три файла VCF, которые содержат в себе информацию о найденных генетических вариантах — инделах и снипах.  При визуальном осмотре файлов сразу же бросается в глаза отсутствие идентификаторов снипов/инделов. Вместо привычных rs-id, варианты индексированы с помощью точек . Поскольку нам необходима для дальнейшего анализа традиционная система обозначения, мы должны произвести аннотирование файлов. Путем метода проб и ошибок я выбрал самую удобный для начинающих геномиков пакет snpEff.

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

Аннотирование  индексов снипов  в VCF-файлах выполняется с помощью  несложной командой (ниже приведен пример командной строки для файла GokSte.vcf).

java -Xmx2g -jar SnpSift.jar annotateMEM -id  ../dbsnp.vcf ../GokSte.vcf > GokSte.annotate.vcf

Очевидно, что древняя ДНК содержит значительное число новых истинных и ложных снипов, которых нет в индексах dbSNP.  В нашем туториале мы ограничимся лишь известными снипами, и поэтому отфильтруем «новельные» снипы.

java -Xmx2g -jar SnpSift.jar filter "(exists ID) & ( ID =~ 'rs' )"  GokSte.annotate.vcf > GokSte.snp.vcf

III.  Фильтрация снипов в vcftools

Как я указывал в предыдущем разделе, файлы VCF содержат в себе информацию о всех найденных генетических вариантах — инделах и снипах.  Несмотря на всю важность инделов в определение вариативности генофонда популяций,  во многих популяционно-генетических исследованиях явное предпочтение отдается снипам. Принимая это во внимание, я решил отсеять инделы в сторону и трансформировать файл VCF в более традиционный формат Plink PED:

./vcftools/bin/vcftools --vcf  GokSte.snp.vcf --remove-indels  --plink --out GokSte

На выходе, мы получил файл Plink PED, о котором мы поговорим в следущей части туториала.

Продолжение следует.

Открытая полемика с историком Носевичем по вопросу о происхождении динарской субклады I2a -часть 3. Интермеццо.

На форуме AFB разместили важную новость которая касается пресловутого «динарского» кластера I2a1b1-L621 (снип включает в себя также и кластер Disles).

Речь идет о результатах теста Geno 2.0 одного поляка из южной Польши.

29.03.2013 известный активист Берни Куллен отрапортовал в Rootsweb-группе Y-DNA-HAPLOGROUP-I:

Я получил результаты Geno 2.0 одного поляка (родом из южной Польши). Он сегодня присоединился к
Проекту I2a на FTDNA, но поскольку он еще не тестировался в FTDNA, его результаты еще не появились на странице нашего проекта.

Но он прислал мне CSV файл со своими результатами, и я сравнил его результаты с результатами 5 других представителей «динарского» кластера I-L621/L147 и одним представителем кластера I-L621 L147-«Disles».

Его результат идентичен результатов других «динариков», с двумя важными исключениями:

у него предковое значение аллели (А) в снипе CTS10228
у него предковое значение аллели (T) в снипе CTS5966

По трем другим снипам CTS10936, CTS11768, CTS4002, его результаты идентичны результатам других динариков:

Примечательно, что у представителя кластера Disles во всех 5 снипах предковые аллели.

Более того, у этого поляка установлено:
1) предковое значение аллеля в снипе, обнаруженном в кластере Disles
2) предковое значение аллеля в снипе, обнаруженном у одного «греческого» динарика,
3) предковое значение аллеля снипа, найденного у одного представителя еврейского субкластера динарского кластера.

Таким образом, поляк принадлежит к ветви динариков, которая наиболее близка к терминальной точки соединяющей генетические линии кластеров Disles и Dinaric.

Что из всего следует?

На мой взгляд эти результаты являются еще одним надежным аргументом в пользу моей версии происхождения динарского кластера (которую я впервые озвучил еще в 2009 году). В свете новых фактов она представляется более предпочтительной, чем классическая балканская версия традиционной попгенетики и «карпийская» гипотеза Носевича.

Новые перспективы коммерческого экзомного тестирования/секвенирования

Поскольку слово «экзом» является совсем свежым заимствованием из английского языка. Наиболее простое определеие: экзом состоит из совокупности экзонов, а экзон — это участок гена (ДНК) эукариот, несущий генетическую информацию, кодирующую синтез продукта гена (белка).  Соответствующие экзонам участки ДНК, в отличие от интронов, полностью представлены в молекуле информационной РНК, кодирующей первичную структуру белка. По мнению некоторых исследователей  соответствуют доменам (структурно автономным областям) в белке и являются первичными генетическими единицами, рекомбинация которых приводит к возникновению в ходе эволюции новых генов и соответственно новых белков. Э. чередуются в структуре гена с другими фрагментами — интронами.
Иными словами, экзом — это совокупность всех участков ДНК, несущих информацию, определяющую экспрессию белка.

Здесь уместно вспомнить недавную видео-лекцию Павла Певзнера «Персональная медицина и ассемблирование геномов: паззл с миллиардом частей» (где-то на мордокніге я давал ссылку). Певзнер, в числе прочего, мимоходом упоминул про недавную работу одного из ведущих сотрудников института Сангера (ведущего центра персональной геномики, одним из исследовательских направлений как раз и является -The 500 Exome Project with collaborators from WTSI, GSK and Lausanne University ).
Речь идет о нашумевшей работе, в котором ученый описывает как на протяжении полугода он в рутинном режиме ежемесячно «проверял экзомы» на предмет анализа экспрессии белков. В ходе работы был не только выявлен целый ряд ранее неизвестных вариантов генов, ответственных за предрасположенность к диабету второго типа, но и произведен анализ динамики белковых изменений.Этот тщательный анализ позволил «излечить» пациентов от диабета.

После того как медиа взбудоражила общественное сознание этой новостью, целая группа коммерческих компаний обратилась к этому, ранее коммерчески неосвоенному типу генотипирования (хотя некоммерческіе исследования ведутся уже не менее десятка лет). Чутко следящая за коньюктурой геномного рынка компания Illumina сразу опустила планку цен на «экзомы» до 200 долларов. Чем не приминула воспользоваться компания 23andme, предлагающая (в качестве посредника, т.к. само типирование проводится в лабах Иллюмины) конечному потребителю продукт по цене 999 долларов.Легкая доступность экзомного тестирования будет иметь свои преимущества, поскольку позволит не только проводить анализ генетических маркеров, определяющие риски, но и анализировать экспрессию белков под воздействием определенных эпигенетических факторов (тип приминяемых медикаментов, питания и т.д.).

Сейчас врачи посылают пациентов на анализ крови, cлюны, мочи и прочих биологических субстанций.Лет этак через 10 врач будет писать в истории болезнии: «Пациенту назначено прохождение годового курса экзомного генотипирования и анализа»

Перед тем как перейти к конкретному примеру, немного сухой теории.

Принципы и платформы экзомного сивенирования

Авторы ряда недавних исследований пытались сравнить коммерческие технологии экзомного секвенирования. Предложенные в течение последних нескольких лет коммерческие платформы секвенирования следующего поколения секвенирования, разработаы с целью секвенирования кодирующей части генома — так называемого экзома. В недавнем исследовании лаборатории Майкла Снайдера (Стэнфорд), авторы сравнили три основных платформы экзомного секвенирования  созданные компаниями Agilent’s SureSelect Human All Exon (50 Mbp), Roche/Nimblegen’s SeqCap EZ v2.0, and Illumina’s TruSeq Exome Enrichment. Платформы сравнивались как между собой, так и с платформой полногеномного сиквенирования  (35x) на примере сиквенирования генома одного и того же человека.

Рис.1. Основные экзомные платформы (Clark и др., Nat Biotech, 2011 год.).

Различия платформ

Для  начала, следует отметить различия между платформами.

Фото: Кларк и др., Nat. Biotech. 2011

Большое количество базовых пар генома (29,45 Мбит), которое (предположительно) составляет «ядро» покрывается всеми  сравниваемыми платформами. В индивидуальном плане, каждая из платформ имеет от 4 до 28 Мбит уникального  «таргетного пространства покрытия». Платформа Agilent лучше подходит для транскриптов Ensembl, в то время как платформа NimbleGen имеет более широкий охват микроРНК.  На приведенной выше диаграмме Венна видно что пересчение множеств покрываемых платморфмами базовых пар генома («таргетное пространство») гораздо больше, чем пересечение каждой из платформ  с множеством таргентного пространства платформы Illumina. Это в первую очередь  объясняется тем, что Illumina в первую очередь нацелена  на секвенирование нетранслируемых областей генома (НТО). Трудно сказать,  является ли это преимуществом платформаы Illumina, или нет. С одной стороны, этот факт, конечно же, представляет определенный интерес исследователям, которые ищут  геномные варианты в регионах НТО. С другой стороны, эта характерная особенность экзомной платформы  Illumina приводит к неизбежному ухудшению качества покрытия генома.  И, действительно, как отмечают авторы отмечают, при секвенировании генома в режиме 50 миллионов 2 × 100 б.п. покрытие генома на платформе Illumina составляет 30x,  для сравнения 60x на платформе Agilent -60x и на платформе NimbleGen — 68x.

Определение таргетной эффективности платформы  и содержание GC

Авторы произвели определение и секвенирование экзома на одном образце — взятом у добровольца европейского происхождения — с использованием всех трех экзомных решений. Каждая экзомная библиотека  получила одну полосу  прочтения (2 × 100 б.п.)  на  Illumina HiSeq 2000 (11-18 Гигобаз в каждой библиотеке).  С помощью программы BWA 99%  полученных геномных «ридов» было отображено на референсную последовательность человеческого генома, примерно 10-15% геномных ридов оказались ПЦР-дубликатами. Затем была вычислена общая таргетная эффективность прочитки  ( при расчете исходили из того, что на каждый экзом приходится примерно 80 миллионов ридов,  и вычисляли процентные доли нуклеотидов покрытых при прочитке 10x, 20x, 30x ). Как пишут авторы исследования: «при всех аналогичных параметров покрытия и числа ридов, платформа NimbleGеn дала более высокий процент опеределения  своих «»таргетных» (целевых) баз, в сравнении с  другими платформами.» Авторы объясняют эффективность платформы использованием высокой плотности перекрывающихся биотинилированных олигов, используемых в платформе NimbleGen для «захвата» экзонов.

Неудивительно, что все платформы продемонстрировали заметное снижение охвата при  повышении и понижении GC целей. Однако, при низкой GC (40% до 20%),  платформа Agilent показал лишь небольшое cнижение глубины прочитки, что, возможно, связано с меньшим числом циклов ПЦР,  и большим числом биолитинированных пробов для захвата РНК, а также с использованием уникальных РНК-пробов.

Обнаружение единичных нуклеотидных вариантов (SNV-ов) и небольших инделов

Как известно, обнаружение небольших вариантов последовательности, особенно SNV-ов, является одной из основных целей секвенирования экзома. Используя нормированные  наборы геномных «ридов» размером в ~ 80 мегабаз, авторы  обнаружили в каждом наборе  экзома SNV-ы  (авторы использовали программное обеспечения GATK). Все три платформы показали высокую степень корреляции между «вызовами»  SNV-ов и высокой плотностью массива генотипов SNP-ов. Как показал эксперимент, референсный аллель имел небольшое преимущество (0.53-0.55) в позициях SNP-ов, и это дает повод утверждать о небольшом отклонении результатов в случае  с ридами, которые содержат геномные варианты. Вместе с тем, не было найденно никаких отклонений в отношении определенных типов субституций.   Как и следовало ожидать, на всех платформах  количество обнаруженных SNV-ов увеличивалось при увеличении охвата генома. Этот рост, однако,  не носил линейный характе; при прочитке 30 миллионов баз,  не было обнаружено более 95% SNV-ов. В общих регионах секвенированных экзомов, платформа NimbleGen  обнаружила большинство SNVs при наименьшем количестве прочтений.

Кроме того, платформа NimbleGen (за счет более эффективного захвата и, следовательно, более глубокого охвата генома) обнаружила большинство инделов  — как в общих регионах экзомов, так и в регионах рефересной последовательности генома. При более низком уровне прочитке, платформа Agilent обнаружила больше инделов в общих регионах экзомов, но при  50 миллионах прочтений, платформа Illumina превзошла Agilent (и, что неудивительно, обнаружила еще много инделов в нетранслируемых регионах). Размер большинства инделов был равне 1 базовой паре, хотя авторы отметили небольшое увеличение размера инделов  до 4-8 базовых пар (что подтверждается путем сравнения генома человека с геномом приматов), а также увеличение за счет отбора против  сдвига рамки мутаций.

Сравнение  результатом экзомного секвенирования с результатами полногеномного секвенирования

Основным достоинством данного исследования можно считать то, что авторы  произвели и  полногеномное секвенирование образца при средней величине покрытия 35х.  Корреляция по гетерозиготным позициям SNP-ов между результатами экзомного секвенирования и полногенмного секвенирования составила 98,5%. Для имитации мультиплексного секвенирования 3 или 6  экзомных библиотек на одну полосу (GAIIx или HiSeq, соответственно), авторы нормировали экзомные продукту  до 50 миллионов  ридов на каждой из платформ. В каждом попарном сравнении продуктов  полногеномного и экзомного секвкенирования, набор данных полногеномного секвенирования был ограничен теми же таргетными регионами, что и сравниваемый экзомный продукт.  Хотя этот шаг  и представляется необходимым для сравнения по методу «подобное к подобному», нужно отметить, что он сводит к минимуму мощь полногеномного секвенирование, которое обеспечивает относительно объективное освещение всех кодирующих областей. Другими словами, это ограничение дает определенное преимущество экзомным продуктам,  поскольку сравниваются только те таргетные регионы, под которые заточена платформа экзомного продукта.

Перекрытие SNV-ов в экзомных продуктах и полно-геномном продукте (Clark и др., Nat Biotech 2011 года.)

Подавляющее большинство SNV-ов в таргетных  регионах было обнаружено как в экзомных, так и в полногеномных регионах, с небольшими различиями. Примечательно, что при попарном сравнении специфические SNV в  экзомных и полногеномных продуктах, ка правило, имеют (1) низкий  доверительный порог, (2) более высокую долю новельных (по отношению к dbSNP) вариантов, и (3) лучшее покрытие в детектирующей платформе. Специфические для полногеномного продукта SNV-ы часто  имеют нулевые риды в экзомных продуктах (вероятно, из-за проблемы с гибридизацией). Напротив, большинство SNV-ов, специфических для экзомных продуктов, было покрыто в полногеномном продукте, хотя их общее число в полногеномном продукте все равно остается ниже, чем в экзомном продукте.

Как становится ясно из вышеприведенного рисунка, число SNV-ов, обнаруженых в экзомных продуктах и полногеномном продукте отлично коррелирует с «»потолком» каждой экзомной платформы. Illumina, которая имеет наибольшее таргетное пространство (особенно в нетранслируемых регионах), имеет наибольшее количество общих SNV-оd.  Число общих SNV-ов в Agilent больше, чем NimbleGen, однако чувствительность NimbleGen в определении истинно-положительных результатов в целевых регионах  гораздо выше, чем на двух других платформах.

Как выбрать экзомную платформу

Авторы приходят к выводу, что все три платформы экзомного секвенирования по-своему  хороши. Выбор платформы, вероятно, зависит от целей, приоритетов и бюджета исследователя. Для малобюджетных проектов, NimbleGen предлагает наиболее эффективное обогащение экзонов (а также микро-РНК). Agilent подходит для охоты за вариантами генома, поскольку обеспечивает более широкий охват, но требует больше  секвенированных данных. Illumina наиболее требовательна в плане секвенированных даннных, но зато обследует нетранслируемые области, могут заинтересовать некоторых исследователей.

Мои собственные практические выводы на основе анализа данных

Летом я освоил навыки анализоа результатов экзомного генотипирования, любезно представленные одним из немногих россиян, участвовавших в пилотном проекте экзомного генотипирования в компании 23andme.

Поскольку слово «экзом» является совсем свежым заимствованием из английского языка. Наиболее простое определеие: экзом состоит из совокупности экзонов, а экзон — это участок гена (ДНК) эукариот, несущий генетическую информацию, кодирующую синтез продукта гена (белка).  Соответствующие экзонам участки ДНК, в отличие от интронов, полностью представлены в молекуле информационной РНК, кодирующей первичную структуру белка. По мнению некоторых исследователей Э. соответствуют доменам (структурно автономным областям) в белке и являются первичными генетическими единицами, рекомбинация которых приводит к возникновению в ходе эволюции новых генов и соответственно новых белков. Э. чередуются в структуре гена с другими фрагментами — интронами.
Иными словами, экзом — это совокупность всех участков ДНК, несущих информацию, определяющую экспрессию белка.

Здесь уместно вспомнить недавную видео-лекцию Павла Певзнера «Персональная медицина и ассемблирование геномов: паззл с миллиардом частей» (где-то на мордокніге я давал ссылку). Певзнер, в числе прочего, мимоходом упоминул про недавную работу одного из ведущих сотрудников института Сангера (ведущего центра персональной геномики, одним из исследовательских направлений как раз и является -The 500 Exome Project with collaborators from WTSI, GSK and Lausanne University ).
Речь идет о нашумевшей работе, в котором ученый описывает как на протяжении полугода он в рутинном режиме ежемесячно «проверял экзомы» на предмет анализа экспрессии белков. В ходе работы был не только выявлен целый ряд ранее неизвестных вариантов генов, ответственных за предрасположенность к диабету второго типа, но и произведен анализ динамики белковых изменений.
Этот тщательный анализ позволил «излечить» пациентов от диабета.
После того как медиа взбудоражила общественное сознание этой новостью, целая группа коммерческих компаний обратилась к этому, ранее коммерчески неосвоенному типу генотипирования (хотя некоммерческіе исследования ведутся уже не менее десятка лет). Чутко следящая за коньюктурой геномного рынка компания Illumina сразу опустила планку цен на «экзомы» до 200 долларов. Чем не приминула воспользоваться компания 23andme, предлагающая (в качестве посредника, т.к. само типирование проводится в лабах Иллюмины) конечному потребителю продукт по цене 999 долларов.

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

Сейчас врачи посылают пациентов на анализ крови, cлюны, мочи и прочих биологических субстанций.
Лет этак через 10 врач будет писать в истории болезнии: «Пациенту назначено прохождение годового курса экзомного генотипирования и анализа»

Как выглядит конечный продукт экзомного тестирования предлагаемый 23andme за 999 зеленых американских рублей?

Это набор из четырех файлов:
1) x.bam
2) x.bai
3) x.pdf
4) x.vcf.

X -это кодовый номер участника. BAM файл являющийся бинарной версией формата SAM (формата множественного выравнивания ДНК по референсному сиквенсу), BAI — индекс контигов в BAM файле. Наконец, VCF — это файл содержащий все «задетектированные» в BAM файле варианты (прежде всего SNPs и INDELs)

Но вернемся к экзомному тестированию.
Cуществует определенная группа лиц, которых больше интересуют вопросы происхождения и генеалогии. Медицинские аспекты, как правило, им неинтересны.

Что интересного могут излечь из экзомных данных ДНК-генеалогии? Не трудно ответить. Большинство ДНК-генеалогов знает принципы наследования ДНК, в первую очередь Y-хромосомы и митохондриального генома (которые наследуются соответственно строго по мужской и женской линии).

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

Первая и основная проблема — это характер экзомного типировния. При экзомном типировании определяются только те снипы и инделы, которые находятся в экзонах. С митохондрионом здесь проблем особых нет — в [человеческом] митогеноме практически все вариативные позиции, являются экзомными, т.е несут генетическую информацию («код» синтеза белка). Поэтому фактически данные экзомного тестирования уже содержат полный сиквенс генома (аналог FGS от FTDNA). Остается только их извлечь. И вот тут появляется другая проблема. Для определения генетических вариантов (т.е различий нуклеотидов в локусах) необходмо провести «выравнивание» анализируемого по референсному сиквенсу. Как известно, в митохондрионе для этих целей используется «классический» сиквенс rCRS (Cambridge Reference Sequence, GenBank:NC_012920.1). Однако в геномных билдах-ассамблеях hg18 и hg19_Chr37, этот референс заменен другим. Поэтому результаты выравнивания митогенома по дефольтным вариантам вышеуказанных билдов дают результаты, сильно отличающиеся от привычного формата.

После замены дефолтного сиквенса на rCRS все получилось. Вот фрагмент из VCF файла, содержащий интересующие нас отличия от rCRS:
#CHROM POS REF ALT
MT 73 A G
MT 195 T C
MT 263 A G
MT 709 G A
MT 750 A G
MT 1438 A G
MT 1888 G A
MT 2141 T C
MT 2706 A G
MT 3106 CN C
MT 4216 T C
MT 4917 A G
MT 5894 A G
MT 7028 C T
MT 8697 G A
MT 8860 A G
MT 9117 T C
MT 10463 T C
MT 11191 C T
MT 11251 A G
MT 11719 G A
MT 11812 A G
MT 12741 C T
MT 13260 T C
MT 13368 G A
MT 13965 T C
MT 13966 A G
MT 14233 A G
MT 14687 A G
MT 14766 C T
MT 14905 G A
MT 15326 A G
MT 15452 C A
MT 15607 A G
MT 15928 G A
MT 16126 T C
MT 16294 C T
MT 16296 C T
MT 16324 T C
MT 16519 T C

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

Убедился и я в этом на примере реальных данных (это представитель Y хромосомной гаплогруппы R1a1).
samtools view -h x.bam Y > Y.sam
samtools view -h -b -S Y.sam > Y.bam
samtools/samtools mpileup -C 50 -ugf chrY.fa Y.bam | /samtools/bcftools/bcftools view -vcg — > Y.raw.vcf

Данный подход позволил обнаружить у тестанта около сотни генетических полиморфизмов (координаты данные по билду hg19):
Y 4058546 0 A C
Y 4058566 0 ta t
Y 4457069 0 tctctcct tct
Y 6028350 0 A T
Y 8149348 0 G A
Y 8566853 0 GCCC GCCCC
Y 8783761 0 C T
Y 8881927 0 GGTGT GGTGTGT
Y 9198243 0 T A
Y 9304866 0 G A
Y 9368340 0 tg tGNg
Y 9384631 0 A C
Y 9385720 0 CGG CG
Y 9909058 0 T A
Y 9930114 0 C A
Y 9931330 0 T A
Y 9938790 0 C A
Y 9938851 0 A T
Y 9938982 0 T C
Y 9939117 0 T A
Y 9952497 0 A G
Y 9982892 0 G A
Y 9982917 0 C A
Y 10007709 0 C A
Y 10007727 0 G A
Y 10007741 0 G A
Y 10011344 0 A G
Y 10011487 0 A G
Y 10011498 0 G C
Y 10011502 0 A G
Y 10011545 0 T G
Y 10011604 0 C CTT
Y 10011648 0 T G
Y 10011673 0 G A
Y 10011677 0 G A
Y 10011698 0 A G
Y 10011878 0 G A
Y 10011935 0 C CT
Y 10011960 0 T C
Y 10011966 0 ATT AT
Y 10012012 0 T A
Y 10013318 0 A G
Y 10028123 0 C T
Y 10028180 0 A G
Y 10029163 0 A G
Y 10029228 0 G A
Y 10029308 0 A T
Y 10029322 0 T C
Y 10029340 0 T C
Y 10029485 0 G C
Y 10029487 0 T A
Y 10029513 0 A G
Y 10029610 0 G A
Y 10029616 0 G T
Y 10029623 0 C T
Y 10029629 0 A G
Y 10029649 0 C G
Y 10029711 0 A C
Y 10043269 0 C T
Y 13241432 0 G T
Y 13241656 0 G A
Y 13243050 0 C G
Y 13243352 0 G A
Y 13244666 0 C T
Y 13244690 0 A G
Y 13254228 0 C T
Y 13262943 0 ACCC ACC
Y 13263091 0 G A
Y 13263304 0 C T
Y 13263364 0 A G
Y 13263374 0 C G
Y 13266266 0 G A
Y 13266286 0 C T
Y 13266301 0 A G
Y 13266368 0 T G
Y 13266377 0 G C
Y 13266499 0 A G
Y 13266520 0 G T
Y 13266556 0 T G
Y 13266560 0 C T
Y 13266587 0 C G
Y 13268187 0 T C
Y 13268361 0 T C
Y 13268377 0 A G
Y 13268521 0 C T
Y 13307425 0 G T
Y 13307562 0 G A
Y 13309174 0 A T
Y 13309226 0 A C
Y 13309239 0 G C
Y 13309262 0 T C
Y 13309348 0 C T
Y 13311223 0 T A
Y 13311491 0 C T
Y 13311501 0 G A
Y 13312579 0 G A
Y 13312666 0 G C
Y 13312729 0 C T
Y 13312756 0 A G
Y 13312789 0 A G
Y 13332277 0 C T
Y 13357224 0 C T
Y 13370991 0 C A
Y 13445929 0 G C
Y 13445957 0 C G
Y 13463779 0 A C
Y 13463831 0 T A
Y 13463837 0 G A
Y 13463860 0 C G
Y 13465055 0 A G
Y 13470805 0 G A
Y 13470834 0 T C
Y 13470855 0 T G
Y 13470880 0 G A
Y 13470897 0 G A
Y 13475849 0 C T
Y 13476553 0 T C
Y 13478387 0 A T
Y 13478445 0 G C,A
Y 13478569 0 T G
Y 13478583 0 T G
Y 13478613 0 A G
Y 13485671 0 T G
Y 13488312 0 C A
Y 13488330 0 A G
Y 13488337 0 C T
Y 13488370 0 G A
Y 13488395 0 A G
Y 13488410 0 A T
Y 13488429 0 A G
Y 13488601 0 A C
Y 13488621 0 A G
Y 13488946 0 A C
Y 13488952 0 T C
Y 13488972 0 C G,T,A
Y 13488988 0 A G
Y 13488992 0 T C
Y 13489043 0 G A
Y 13489069 0 A C,G
Y 13489077 0 T C
Y 13489206 0 C G
Y 13489220 0 T C
Y 13489234 0 T C
Y 13489255 0 A G
Y 13489292 0 A G
Y 13489300 0 A G
Y 13492264 0 C A
Y 13500410 0 T G
Y 13500424 0 T C
Y 13500443 0 T C
Y 13502048 0 C T
Y 13524378 0 T C
Y 13524752 0 G T
Y 13524761 0 C T
Y 13524873 0 T C
Y 13537129 0 G A
Y 13537569 0 A T
Y 13537581 0 C T
Y 13541022 0 C A
Y 13541053 0 CA CATA
Y 13541068 0 T C
Y 13541199 0 A G
Y 13541232 0 A T
Y 13541288 0 G A
Y 13541293 0 ATTT ATT
Y 13541420 0 A C
Y 13541454 0 T C
Y 13541478 0 G T
Y 13541520 0 C T
Y 13541556 0 A C
Y 13541561 0 T G
Y 13541584 0 C G
Y 13572922 0 A C
Y 13572932 0 T C
Y 13572999 0 A G
Y 13573033 0 A C
Y 13573108 0 G C
Y 13573152 0 C A
Y 13573216 0 G A
Y 13573240 0 C T
Y 13573271 0 G T
Y 13595280 0 T C
Y 13687807 0 T G
Y 13688825 0 C G
Y 13689634 0 T C
Y 13689668 0 C G
Y 13689755 0 G C
Y 13690562 0 C T
Y 13694899 0 G A
Y 13694929 0 G A
Y 13694956 0 C G
Y 13694983 0 T A
Y 13695051 0 T G
Y 13726074 0 T A
Y 13726129 0 C G
Y 13842718 0 G C
Y 14482235 0 C A
Y 14485120 0 G A
Y 14498990 0 C T
Y 14771478 0 A T
Y 14898094 0 A G
Y 14958218 0 C T
Y 15026424 0 A C
Y 15027529 0 T G
Y 15930958 0 ccttcttcctc cCTTCTTCCTCCTcttcttcctc
Y 16751825 0 A G
Y 16832517 0 T C
Y 17231616 0 A G
Y 21154004 0 A C
Y 21154323 0 G A
Y 21154426 0 G A
Y 21154466 0 T A
Y 21208056 0 A G
Y 21208066 0 C G
Y 22260237 0 C T
Y 22510104 0 G A
Y 22510163 0 T A
Y 23473201 0 T A
Y 23800360 0 T G
Y 23805478 0 C A
Y 24008079 0 T A
Y 28582510 0 G C
Y 28582566 0 C G
Y 28582605 0 T C
Y 28582622 0 G A
Y 28582676 0 G A
Y 28582685 0 C A
Y 28582863 0 A G
Y 28582865 0 A G
Y 28582921 0 A G
Y 28582932 0 G A
Y 28583310 0 C T
Y 28583314 0 A G
Y 28583382 0 G C
Y 28583394 0 T C
Y 28583410 0 C G
Y 28583415 0 T C
Y 28583431 0 A T
Y 28583432 0 A G
Y 28583590 0 A C
Y 28586782 0 G A
Y 28586959 0 T C
Y 28587232 0 T C
Y 28689055 0 G T
Y 28709343 0 A G
Y 28780767 0 A C
Y 28780823 0 T A
Y 28780883 0 G A
Y 28815270 0 C A
Y 28815656 0 T C
Y 28816806 0 T C
Y 28816831 0 C T
Y 28816870 0 T G
Y 28816948 0 C G
Y 28817276 0 T G
Y 28817286 0 T G
Y 28817559 0 T G
Y 28817636 0 G A
Y 58856145 0 G C
Y 58883603 0 A T,C
Y 58883784 0 T A
Y 58883834 0 A T
Y 58893627 0 A T
Y 58968939 0 G A
Y 58975896 0 T C
Y 58981639 0 cctccactcca cCTCCActccactcca
Y 58982160 0 G T
Y 58982559 0 A C
Y 58982671 0 tcttccttc tcttc
Y 58985524 0 T G
Y 58996230 0 G A
Y 58996257 0 G T
Y 58999765 0 C T
Y 58999773 0 G A
Y 59001429 0 G A
Y 59001608 0 C T
Y 59001620 0 A C
Y 59001647 0 G A
Y 59001685 0 G C
Y 59001722 0 G A
Y 59001753 0 T C
Y 59001773 0 A C
Y 59001782 0 C A
Y 59001792 0 T C
Y 59001960 0 T A
Y 59002047 0 C G
Y 59002139 0 G T,A
Y 59005179 0 C A
Y 59010280 0 A G
Y 59015256 0 T A
Y 59017005 0 A G
Y 59017181 0 T A
Y 59017206 0 A G
Y 59017378 0 T G
Y 59017384 0 ag aGg
Y 59018341 0 C G
Y 59020728 0 A G
Y 59022718 0 A G
Y 59022723 0 C T
Y 59022734 0 C T
Y 59022768 0 A G
Y 59027525 0 A G
Y 59027700 0 A C
Y 59027882 0 T G
Y 59029728 0 C T

В продолжение о Y-хромосоме. Совершенно ясно, что большинство из снипов, обнаруженных у протестированного ранее не были известны, и поэтому отсутсвуют в официальном списке ISOGG.
C помощью незамысловатой комманды grep -f snps ISOGGsnps я нашел лист известных ISOGG-снипов, которые также присутствуют и в данных тестанта>

L146 R1a M420 rs17250535 21882589 23473201 T->A
L265 R1b1a2 rs9786882 8209348 8149348 A->G
L269 G 13467612 14958218 T->C
M173 R1 P241; Page29 rs2032624 13535818 15026424 A->C
M201 G rs2032636 13536923 15027529 G->T
M379 I2a2a2 rs2032636 13536923..13536924 15027529..15027530 GT->del
M420 R1a L146 rs17250535 21882589 23473201 T->A
P241 R1 M173; Page29 rs2032624 13535818 15026424 A->C
Page7 R1a1a1 rs34297606 13008998 14498990 C->T
Page29 R1 M173; P241 rs2032624 13535818 15026424 A->C
Page83 P rs35361051 13407488 14898094 A->G

Исходя из вышеприведенной таблицы, Y-хромосомная сигнатура тестируемого в классическом ФТДНА-шном виде будет выглядеть следущим образом:

L269-:M201-:M379-:Page83+:L265-:Page29(M173+,P241+):M420+:Page7+.

Из чего еrgo (следует), что тестант принадлежит к группе R1a1a1

ALDer: анализ генеографии дагестанских народов в эволюционной перспективе

В свете наших споров с уважаемой Казимой Булаевой  (один из ведущих российских генетиков), я решил продемонстрировать робастность метода ALDer, предложенного в статье Loh et al.2012 в анализе демографически сложных популяций Дагестана.

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

Программа ALDer использовалась в двух режимах.

Первый режим — 1-reference weighted LD curve, второй режим — 3+ reference weighted LD curve. Термины нуждаются в грамотном переводе в русскоязычную терминологию, так что пока я оставил их в том виде, в каком они есть.

В режиме 3+ references для вычисления значения корреляции использовались следущие популяции:

ItalianCenter;Sicilian;Sardinian;German;Lithuanian_V;Lithuanian;Latvian;Belarusian;Swedish;Polish_V;Russian_V;Russian_Center;Latvian_V;Inkeri;Russian_South;Ukrainian_V;Slovakian;Czech;Sorb;Estonian;Ukrainian;Belarusian_V;UkrainianEast;UkrainianWest;Mordovian;CEU;CEU_V;British;French;Orcadian;GermanSouth;GermanNorth;German_V;Bulgarian;FinnishNorth;Cirkassian;Russian_cossack;Saami;Udmurd;Komi;Karelian;Vepsa;Mari;Bashkir;Nenets;Hant;Chuvash;Mansi;FinnishSouth;Polish;Gagauz;Moldavian;Romania;Bosnian;Adygei;Croatian;Serbian;Slovenian;Montenegrin;Macedonian;Kosovar;Austrian;Greek_Azov;Greek_Center;Greek_North;Greek_South;Tatar_Crim;Azeri;Tadjik;Kyrgyz;Kazakh;Georgian;Georgian_Imereti;Georgian_Laz;NorthOssetian;Armenian;Kumyk;Chechen;Turk;Turkmen;Uzbek;Mongol;Karakalpak;Lak;Balkarian;Lezgin;Abhkasian;Kalmyk;Syrian;Kurd;Tabassaran;Hakas;Altaic;Tatar_Kryashen;Tartar_Mishar;Parsi;Avar;Nogai;Italian-North;Hungarian

I.

Итак, начну с результатов ногайцев в тесте «3+ reference populations».
Результаты свидетельствуют о наличии синхронного адмикса у предков современных дагестанских ногайцев, имевшего место быть в интервалме между 17.20 +/- 3.32 и 12.49 +/- 2.55 поколениями до настоящего времени.

DATA: success 0.00052 Nogai Sorb Uzbek 5.17 2.15 2.63 15% 17.20 +/- 3.32 0.00006274 +/- 0.00000886 19.29 +/- 8.98 0.00002377 +/- 0.00000606 20.01 +/- 7.60 0.00001471 +/- 0.00000378
DATA: success 0.0022 Nogai Ukrainian-West Karakalpak 4.89 3.08 4.01 17% 12.49 +/- 2.55 0.00006670 +/- 0.00000793 14.31 +/- 3.68 0.00000975 +/- 0.00000316 14.85 +/- 3.70 0.00003094 +/- 0.00000499

Адмикс был двухкомпонентный — преобладающий центральноазиатский, и восточноевропейский. Внизу приведены данные по нижнему значению величины адмикса

Сорбы Mixture fraction % lower bound (assuming admixture): 44.9 +/- 7.2
Каракалпаки Mixture fraction % lower bound (assuming admixture): 55.8 +/- 4.4

Узбеки Mixture fraction % lower bound (assuming admixture): 72.1 +/- 16.1
Западные украинцы Mixture fraction % lower bound (assuming admixture): 52.0 +/- 11.7

II.

Cледущий пример — кумыки.Из всех возможных 2-референсных комбинаций кривых взвешенной LD статистически значимой оказалась только одна комбинация:

DATA: success 3.3e-06 Kumyk Italian-Center Hakas 6.07 2.71 5.82
23% 18.15 +/- 2.99 0.00004734 +/- 0.00000710 15.15 +/- 4.49 0.000004
61 +/- 0.00000170 19.14 +/- 3.29 0.00003496 +/- 0.00000519

Это весьма примечательный результат. Как видно из результатов, кривые угасания LD (cцепления по неравновесию) обеих популяций имеют положительную корреляцию между собой. Время двухстороннего адмикса — 18.15 +/- 2.99, то есть интервал между серединой 14 века и cерединой 16 века.

Примечательно, что величина нижнего порога «cредиземноморского» (Italian-Сenter) компонента адмикса выше чем аналогичная величина «тюркского» (Hakas) компонент адмикса у кумыков (см.ниже):

«итальянцы» -Mixture fraction % lower bound (assuming admixture): 47.9 +/- 8.3

хакасы — mixture fraction % lower bound (assuming admixture): 12.9 +/- 1.4

Можно поспекулировать по поводу исторических интерпретаций данных результатов. Если мы вслед за некоторыми генетиками будем рассматривать хакасов как наиболее близкую к древним тюркам популяцию, то можно предположить общее происхождение тюркского «компонента» кумыков с кыпчаками, либо (что менее вероятно) хазарами. Здесь много свободного места для спекуляций.

«Итальянская» часть адмиксf вызывает больше вопросов, чем ответов. Освежив свои неглубокие познания в истории дагестанского региона, смог вспомнить лишь смутные упоминания о присутствии итальянцев в Дагестане в 14-15 веках. «Согласно Фануччи*, генуэзцы выстроили и заселили поселение Кубачи в Дагестане …»; (Исторические записки. Том 3.,1938 г., Зевакин Е. С. и Пенчко Н. А. «Очерки по истории генуэзских колоний на Западном Кавказе в XIII—XV вв.» (72-129)).
Сами кубачинцы именуют свой аул грозным именем Угбуг, что означает «убийцы, губители». Но, впрочем, имеется в виду не свирепость кубачинцев, а непобедимое кубачинское оружие. В Кубачи живут мастера, которым нет равных. Они исполняют любые тонкие работы с металлом, но их призванием, прежде всего, всегда было оружие и доспехи; но это не просто ремесло, а сакральное искусство, философия, если угодно — магия. Слава кубачинского оружия — распространилась на весь мир.

С VI века в арабских источниках упоминается название Зирихгеран. Это название на фарси означает «страна тех, кто делает доспехи», по-русски получается нечто вроде «бронники» или «кольчужники». Около 1467 года впервые упоминается имя Кубачи (или Гюбечи), слово это тюркское, означает «бронники, изготовители доспехов».

Все три имени аула и живущего в нем народа означают одно — ремесло. И в этом странность: получается, что кубачинцы — народ без древнего, исконного наименования.

Но есть еще одно имя, четвертое. Соседи (лакцы, кайтаги и лезгины) кубачинцев называют странно — пранг-капур, то есть франки. Более того, сами угбуги-кубачи считают, что их род из Франции.

Первым о потомках европейцев на Кавказе сообщил полковник Иоганн Густав Гербер (умер в 1734 году), — он побывал в тех краях в 1727 году. Спустя полвека академик Иоганн Антон Гильденштедт (1745–1781) в описании своего путешествия по Кавказу сообщил, что в Кобачи живут потомки генуэзцев. Они бежали в горы от войск Чингисхана в 1220–1230-х годах, долго оставались христианами, потерпели гонения, скрывались и только после долгих преследований приняли ислам.

Не попал ли «итальянский» адмикс к кумыкам от кубачинцев?

III.

Наконец, лезгины. Точно также как и в случае с лакцами, поражает отсутствие статистически значимой двух-референсной корреляции кривой экспонентного угасания LD (неравновесного сцепления).

Согласно результатам предварительного теста, только две популяции — башкиры и казахи — имеют однорефренсную взвешенную LD-кривую с лезгинами

Pre-test: Does Lezgin have a 1-ref weighted LD curve with…… Bashkir: YES (z = 1.99) Kazakh: YES (z = 2.12)

Датировка казахского и башкирского адмиксов в популяции лезгинов также представляется мне логичной 8+-4 поколений [башкиры] и 9+-3 поколений тому назад [казахи]:

DATA: failure (warning: decay rates inconsistent) 2.6e+03 Lezgin Bashkir Kazakh 0.00 1.99 2.12 194% 500.00 +/- inf -341600.24428451 +/- inf 7.89 +/- 3.96 0.00000246 +/- 0.00000113 9.02 +/- 3.24 0.00000325 +/- 0.00000154

Поскольку decay rates (скорости угасания) казахского и башкирского адмикса не коррелируют между собой, то их источник был разный.Поскольку оба адмикса недавние — в пределах 100-300 лет назад — то как и ожидалось, %-доля этих адмиксов в генофонде лезгинов невелика.

[башкиры] Mixture fraction % lower bound (assuming admixture): 2.7 +/- 1.1
[казахи] Mixture fraction % lower bound (assuming admixture): 3.6 +/- 1.5

Вывод — смешение башкиров/казахов с лезгинами носило случайный и несистематический характер.