Прошедшие две недели я посветил отработке новой методики увеличения аккуратности определения вариантов снипов в геномах древних образцов. Я решил отказаться от предыдущих способов, когда с помощью samtools и GATK сначала генерировались файлы пайлапа, а потом из этой кучи возможных вариантов случайным образом выбирался аллелель и дублировался (т.е образец получал гомозиготные варианты). Проблема этого подхода выяснилась во время импутирования геномов, искусственная псевдогаплоидность древних геномов приводила к громадному искажению в сторону референсных геномов. Я решил упростить сложности и теперь вместо приведения генотипов к псевдогаплоидности, я определяю в GATK UnifiedGenotyper 38 миллионов известных снипов с таким условием, что алгоритм сам отбирает только те аллели, которые заданы в dbsnp как референсный и альтернативный аллель снипа.
В принципе, после долгих головоломок, удалось получить приемлимый алгоритм действий.
Я апробировал его на 55 опубликованных палеогеномах из балтийского региона (Литва, Латвия и Эстония) времен мезолита, раннего, среднего и позднего неолита, а также бронзового времени.
Для большой точности я ограничился только теми образцами, для которых удалось определить генотипы как минимум половины из 38 миллионов снипов dbsnp.

Sample Culture
Donkalnis6 Baltic_EMN
Gyvakarai1 Baltic_LN
Kivutkalns19 Baltic_BA
Kivutkalns207 Baltic_BA
Kivutkalns209 Baltic_BA
Kivutkalns215 Baltic_BA
Kivutkalns222 Baltic_BA
Kivutkalns25 Baltic_BA
Kivutkalns42 Baltic_BA
Kretuonas2 Baltic_EMN
Kretuonas4 Baltic_EMN
MA969 Baltic_BN
MA973 Baltic_LN
Plinkaigalis242 Baltic_LN
Spiginas1 Baltic_EMN
Spiginas2 Baltic_LN
Spiginas4 Baltic_Mesolithic
ZVEJ25 Baltic_Mesolithic
ZVEJ27 Baltic_Mesolithic
ZVEJ31 Baltic_EMN
ZVEJ32 Baltic_Mesolithic

Перед тем как использовать полученный набор в downstream aнализе, я решил посмотреть насколько точно определилось структурное разделение генофонда этих древних геномов.
Я использовал программы peddy, ATK, а также разбиение на фракции компонентов происхождения в программах iAdmix и fastNGSadmix.
На графиках видно, что в принципе основная масса этих геномов проецируется на то место в пространстве генетического разнообразия современных популяций людей, где оно и должно находится c точки зрения здрового смысла.

#family_id sample_id paternal_id maternal_id sex phenotype het_call_rate het_ratio het_mean_depth het_idr_baf ancestry-prediction PC1 PC2 PC3
Donkalnis6 Donkalnis6_Donkalnis6 0 0 0 -9 0.996 0.3029 -2 0 EUR -0.4471 -1.312 0.4822
Gyvakarai1 Gyvakarai1_Gyvakarai1 0 0 0 -9 0.9214 0.2377 -2 0 AMR -0.09174 -1.431 0.4644
Kivutkalns19 Kivutkalns19_Kivutkalns19 0 0 0 -9 0.9923 0.3483 -2 0 EUR -0.5558 -1.044 0.803
Kivutkalns207 Kivutkalns207_Kivutkalns207 0 0 0 -9 0.997 0.3443 -2 0 EUR -0.4681 -1.071 0.5988
Kivutkalns209 Kivutkalns209_Kivutkalns209 0 0 0 -9 0.9596 0.2518 -2 0 EUR -0.4277 -1.495 0.4507
Kivutkalns215 Kivutkalns215_Kivutkalns215 0 0 0 -9 0.973 0.2798 -2 0 EUR -0.2305 -1.201 0.901
Kivutkalns222 Kivutkalns222_Kivutkalns222 0 0 0 -9 0.8608 0.1615 -2 0 AMR -0.4777 -1.456 0.3636
Kivutkalns25 Kivutkalns25_Kivutkalns25 0 0 0 -9 0.8956 0.1933 -2 0 AMR -0.5087 -1.067 0.5996
Kivutkalns42 Kivutkalns42_Kivutkalns42 0 0 0 -9 0.8412 0.1575 -2 0 AMR -0.1253 -1.393 0.4066
Kreutonas2 Kreutonas2_Kreutonas2 0 0 0 -9 0.8462 0.1364 -2 0 EUR -0.4288 -1.337 0.6583
Kreutonas4 Kreutonas4_Kreutonas4 0 0 0 -9 0.9985 0.3136 -2 0 EUR -0.3243 -1.217 0.7842
MA969 MA969_MA969 0 0 0 -9 0.8092 0.1161 -2 0 AMR -0.2649 -1.263 -0.2799
MA973 MA973_MA973 0 0 0 -9 0.9482 0.2736 -2 0 EUR -0.3808 -1.319 -0.2429
Plinkgailis242 Plinkgailis242_Plinkgailis242 0 0 0 -9 0.9777 0.2811 -2 0 EUR -0.5622 -1.108 0.341
Spiginas1 Spiginas1_Spiginas1 0 0 0 -9 0.9943 0.3158 -2 0 EUR -0.4762 -1.402 0.7969
Spiginas2 Spiginas2_Spiginas2 0 0 0 -9 0.974 0.2945 -2 0 EUR -0.5128 -1.521 0.3943
Spiginas4 Spiginas4_Spiginas4 0 0 0 -9 0.8427 0.1399 -2 0 AMR -0.3 -1.208 0.6467
ZVEJ25 ZVEJ25_ZVEJ25 0 0 0 -9 0.969 0.2344 -2 0 EUR -0.2371 -1.254 1.072
ZVEJ27 ZVEJ27_ZVEJ27 0 0 0 -9 0.5763 0.0387 -2 0 UNKNOWN -0.2384 -1.622 -0.7302
ZVEJ31 ZVEJ31_ZVEJ31 0 0 0 -9 0.6926 0.06053 -2 0 UNKNOWN 0.04159 -1.332 -0.1725
ZVEJ32 ZVEJ32_ZVEJ32 0 0 0 -9 0.7095 0.06971 -2 0 UNKNOWN -0.06001 -1.699 -0.3068

Подробное разложени образцов по компонентамм можно посмотреть в этой таблице

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 не так сглаживает минорные компоненты происхождения (что позволяет легче выделить в пространстве главных компонент кластеры):

test (1)

Вторая фаза нового проекта

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

Как я уже рапортовал ранее, самой сложной из запланированных на втором этапе задач являлась импутирование (импутация) отсутствующих генотипов.  Читатели моего блога помнят, что две предыдущие экспериментальные попытки импутирования больших выборок     — в 2013  и в 2015  — закончились неудачно (или, если говорить точнее, качество импутированных генотипов не оправдало моих завышенных ожиданий). В предыдущих опытах я задействовал мощную комбинацию программ ShapeIT и IMPUTE и  метод импутирования снипов за счет использования большой референсной панели аутосомных гаплотипов (из 1000 genomes),  гарантирующей более аккуратное определение генотипов.

На этот раз, я решил не повторять ошибок, и обратился к использованию других программ — в частности , к  Minimac3, хорошо зарекомендовавшую себя в работе с геномами 1000G.  К моему счастью, я набрел на недавно появившиеся публичные сервера, работающие с «облачным» сервисом импутирования Cloudgene. геномов.
Серверы импутирования геномов позволяют использовать полную референсную панель гаплотипов для точного определения недостающих генотипов в анализируемых данных. Пользователи подобных серверов могут загружать (предварительно фазированные или несфазированные) данные генотипов на сервер. Процедура импутирования  будет осуществляться на удаленном сервере, и по окончанию этого процесса рассчитанные данные доступны пользователю для скачивания. Наряду с импутированием, подобные сервисы позволяют провести процедуру контроля качества (QC) и фазировки данных в качестве предварительного этапа процесса импутирования генотипов.

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

Основная рабочая лошадка сервиса — это комбинация двух или трех программ — две програмы для фазирования диплоидных генотипов в гаплоидную фазу  ShapeIT и Hapi-UR , а в качестве основного ПО для самого процесса импутирования (определения) недостающих генотипов — вышеупомянутую программу Minimac3.

Описание эксперимента с импутированием генотипов на удаленном сервере

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

  1. европейцы (европейцы + кавказцы) — 1715 образцов -87169 снипа
  2. азиаты (+американские аборигены и аборигены островов Тихого Океана) — 2356 образцов — 87044 снипа
  3. африканцы — 1054 образца — 86754 снипов
  4. палеогеномы древних жителей Евразии, Африки и Америки -340 — 594500 снипов
  5. смешанные группы — преимущественно мозабиты, пуэрто-риканцы и др.
QC-Report
На рисунке показана корреляция между частотами аллелей в изучаемоей выборке (здесь: европейская когорта) и частотами аллелей в референсной панели

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

Несмотря на значительную скорость обработки генотипов на удаленном сервере, к настоящему времени эксперимент еще не доведен до конца.  Пока я планирую ограничиться импутированием генотипов в 3 первых когортах (т.к. импутирование палеогеномов с помощью современных референсных панелей гаплотипов вероятнее всего приведет к искажению истинного разнообразия палеогеномов за счет проекции на современные группы населения, хотя авторы статьи Gamba et al. 2014 в сопроводительном материале к своей статье утверждают обратное).

После окончания фазирования и последующей обработки генотипов европейской когорты в программе Plink (были отсеяны все варианты с вероятностью ниже 0.9) ,  я получил выборку из 1715 европейцев с 25 215 169 снипами против изначальных 87169, т.е число снипов в выборке увеличилось в 290 раз!
В азиатской когорте соотношение импутированных генотипов к исходным составило чуть меньшую величину 19 048 308 / 87044 = 219.

Проверка результатов

Разумеется, все полученные результаты нуждались в дополнительной проверке качества генотипирования.
Cначала я объединил импутированную европейскую когорту с когортой палеогеномов (которая не была импутирована) и рассчитал в программе PLINK 1.9 матрицу IBS (т.е. сходства образцов в выборке между собой, эта метрика отдаленно напоминает Global Similarity в клиентских отчетах 23andme), а затем усреднил данные по популяциям и произвел по усредненным значениям иерархическую кластеризацию по признакам сходства (IBS, identity by state). Результат превзошел все мои пессимистические ожидания

 

 

Как становится очевидно из приведенной выше кластерограммы,  в целом взаимное расположение популяций в кластерах соответствует (в общих чертах) взаимному географическому положению. Присутствуют, правда, и некоторые огрехи. Так, например, венгры очутились в одном кластере с русскими из Курска,  норвежцы — с русскими из Смоленска, а усредненные «русские» — с американцами европейского происхождения из штата Юта и французами. Трудно сказать, в чем здесь причина, тем более что матрица была составлена по значениям IBS (идентичности по состоянию), а не IBD (идентичности по происхождению).  Более подробные данные о попарных значениях IBS между популяциями выборки можно посмотреть в этой таблице

Импутированная азиатская когорта (несмотря на расширение географии за счет включения образцов коренного населения Америки и аборигенов бассейна Тихого океана)  тоже  оказалась на удивление надежной. Я пока не буду останавливаться на подробностях изучения этой когорты, вместо этого я размещаю здесь результаты MDS- мультдименсионального шкалирования образцов выборки, образованной в ходе слияния 2 импутированных когорт (европейской и азиатской) с 1 неимпутированной (палеогеномы). Цветовое обозначение точек соответствует определенным кластерам, выявленных в выборке с помощью алгоритма MCLUST (cледуя рекомендациям Диенека Понтикоса). Всего этих кластеров 15 и они обозначены последовательностью чисел от 1 до 15, и каждый из этих кластеров имеет свою четкую географическую привязку:

  • 1 — кластер популяций ближнего Востока и  Анатолии
  • 2 — кластер популяций северного Кавказа
  • 3 — «индоевропейский» кластер древних популяций Синташта, шнуровой культуры, Ямной культуры и т.д.
  •  4 — кластер аборигенных жителей Америки (эскимосов и индейцев)
  • 5 — суперкластер популяций средиземноморского и восточноевропейского региона
  • 6 — сибирский кластер алтайских и самодийских популяций
  • 7 — кластер популяций западной и северной Европы
  • 8 — кластер палеосибирских популяций (таких как чукчи, ительмены и коряки)
  • 9 — кластер аборигенных (австронезийских и тай-кадайских) популяций юго-восточной Азии (даи, атаяла и ами)
  • 10 — кластер неолитических популяций
  • 11 — еще один ближневосточно-средиземноморский кластер (ашкеназим, сардинцы и так далее)
  • 12 — кластер североиндийских популяций
  • 13 — кластер центральноазиатских популяций
  • 14 — поволжские популяции
  • 15 — разные групп индусов

 

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

Несколькими постами ранее ув. Сергей Козлов подготовил замечательный по своей глубине русскоязычный обозор новой статьи Allentoft et al. 2015 (еще раз выражаю свою благодарность). В этом обзоре были затронуты преимущественно технические вопросы, в то время как в аналогичном разборе на сайте генофонд.ру было пересказано общее содержание статьи:  » Cтатья большого международного коллектива, опубликованная 11 июня в журнале Nature, посвящена исследованию геномов популяций Евразии в бронзовом веке (изучен период от 3000 до 1000 лет до н.э.). Первый автор Мортен Aллентофт (Morten E. Allentoft) и ведущий автор Эске Виллерслев (Eske Willerslev) представляют Центр географической генетики Музея естественной истории Университета Копенгагена, Дания. Эта статья вызвала огромный интерес у специалистов по истории популяций человека — ведь в ней представлен анализ самого большого массива древних геномов из разных археологических культур эпохи бронзы. На основании анализа древних геномов авторы пробуют реконструировать древние миграции и распространение археологических культур во времени и пространстве. В бронзовом веке, начало которого датируют временем 3500-3300 лет до н.э., в производстве орудий и оружия камень все больше уступает место металлу. Это сопровождается   радикальными культурными и социальными изменениями в жизни людей. Они касаются не только хозяйственного уклада – возникает новое понимание имущественных отношений, семьи и личности. Основной вопрос, на который попытались ответить авторы статьи — были ли эти изменения результатом передачи культурных навыков или результатом миграций населения. Иными словами, «была ли это циркуляция людей или идей». Важнейший вопрос — связаны ли эти события с распространением индоевропейских языков, на которых сейчас говорит большая часть человечества.»

Лавина публикаций древних геномов (кроме вышеупомянутых статей Allentoft et al. 2015, Haak et al. 2015, летом опубликовались статьи Pinhasi et al. Optimal Ancient DNA Yields from the Inner Ear Part of the Human Petrous Bone,  и Fu et al.  An early modern human from Romania with a recent Neanderthal ancestor, однако к сожалению, количество снипов в большинстве образцов палеогеномов недостаточно для проведения развернутых анализов вместе с палеогеномами из других статей) заставила меня ускорить подготовку своей сводной выборки референсных образцов популяций (об этом я писал в предыдущих записях).

Я решил отказаться от полного импутирования древних геномов (очевидно бессмысленного занятия, так как у нас нет надежной референсной панели для импутирования выборочно секвенированных палеогеномв), и вместо этого ограничился импутированным (с помощью панели 1000 Genomes) набором снипов в контрольном наборе популяций лаборатории Райха (Affymetrix Human Origins Fully Public Dataset), этот набор использовался в статье Lazaridis et al. 2014.  

Разумеется, ни одна процедура «импутирования генома» (imputation of genome -сложнопереводимый на русский язык термин) не обходится без ошибок. Поэтому перед тем как приступить к самому анализу, я провел проверку качества выборки. На этот раз, я использовал  инструментарий Python — PyGenClean. Этот инструментий существенно облегчает стандартизацию генетических данных и  контроль качества выходных данных платформы генотипирования. Он минимизирует ошибки манипулирования данными, и ускоряет процесс очистки данных от потенциальных ошибок генотипирования,  а также позволяет составлять информативные графики и автоматически оценивать предварительные параметры последующего статистического анализа.

После отсеивания снипов c низким качеством и индивидов с низкой степенью генотипирования, а также снипов с существенным отклонением от равновесия Харди-Вайнберга, я посмотрел оставшиеся образцы на предмет наличия в выборке «оutliers» (так называемых «статистических выбросов»). До процедуры нахождения выбросов график главных компонент выглядел следующим образом:

И после нахождения выбросов (т.е образцов со стандартным отклонением больше 5 сигм)

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

 

Поэтому такую очистку данных лучше проводить в два захода, один — для современных образцов, а другой — для палеогеномов. А затем полученные «качественные» выборки соединять в общую контрольную выборку.
В нашем случае, я так и поступил, получив выборку из 2250 этнопопуляционных образцов и 155 000 снипов.

Вот так выглядит взаимное расположение образцов геномов на PCA графике.

 

Caucasian, North-African, Afro-American, AG2, South-European, Alberstedt-LN, Native-American, Siberian, African, East-Asian, Near-Eastern, Atayal-Coriell, Native-Australian, Australian, Australian-ECCAC, East-European, Baalberge-MN, South-Asian, Volga-Ural, West-European, Bell-Beaker-LN, North-Indian, BenzigerodeHeimburg-LN, South-Indian, Ancient-African, American, Oceanian, South-East-Asian, Arctic, Corded-Ware, Near-East, Denisovan, Denmark-Carlstrup, Denmark-Falshoy, Denmark-Marbjerg, Denmark-Sebberskole, Esperstedt-MN, EuropeanIronAge, North-European, Halberstadt-LBA, Central-Asian, Hixton, Href, HungaryGamba-BA, HungaryGamba-CA, HungaryGamba-EN, HungaryGamba-HG, HungaryGamba-IA, Iceman, Karelia-HG, Karsdorf-LN, Kostenki14, LaBrana1, LateDorset, LBK-EN, WHG, MA1, Mezmaiskaya, MiddleDorset, North-Greek, South-Italian, Piramalai-Kallars, Poland-Polwice, Poland-Szczepankowice, Poland-Unetice, Poland-Chociwiel, Samara-HG, Saqqaq, East-Aasian, Spain-EN, Spain-EN-relative-of-I0410, Spain-MN, Starcevo-EN, Stuttgart, Sweden-Abekas, Sweden-Angamollan, Sweden-Visby, SwedenSkoglund-MHG, SwedenSkoglund-MN, SwedenSkoglund-NHG, Thule, Unetice-EBA, Ust-Ishim, Vindija, Yamnaya

Именно эту выборку я положил в основу своего нового тестового калькулятора K13 (о его создании я расскажу позже) — предназначенный для анализа «глубокого» происхождения популяций. Как всегда, модель нового калькулятора основана на базовой модели известного DIYDodecad калькулятора. Впервые я остался более или менее удовлетворен полученными результатами. Думаю, что от этой модели можно плясать дальше. И хотя модельная кластеризация с помощью алгоритма Mclust дает основание полагать, что используемая мной выборка из 2230 геномов наилучшим образом (т.е без неизбежного при больших значениях K вырождения компонентов) описывается моделью из 8 кластеров, я остановился на K=13 т.е 13 кластерах:

  1. Amerindian — модальный компонент американских индейцев

 

 

  • ANE — модальный компонент северных евразийцев, изолирован из общего с WHG кластера — наивысшие значения в древнесибирских образцах MA1, AG2, а также у андроновцев, синаштинцев, представителей ямной культуры, шнуровиков и т.д. Из ныне живущих популяций самый высокий процент у калашей. Практически совпадает с ANE в статье Lazaridis et al. 2014

 

 

  • Arctic — модальный компонент с пиком в популяциях коряков, чукчей, ительменов и эскимосов

 

 

  • ASI — модальный компонент южноиндийских популяций, у современных популяций наивысший процент у онге, идентичен ASI в работе Reich et al. 2009.

 

 

  • Caucas(us)-Gedrosia — идентичен кластеру, открытому в 2011 году Диенеком Понтикосом

 

 

  • EastAsian — модальный компонент жителей восточной Азии

 

 

  • ENF — компонент древних европейских земледельцев неолита, пик в образцах палеогеномов культуры линейно-ленточной керамики. Тождественен аналогичному компоненту в работах популяционных генетиков (Lazaridis et al. 2014, Haak et al. 2015). В современных этнопулах — наивысшие значения у сардинцев, корсиканцев и басков.

 

 

  • NearEast — модальный компонент жителей ближнего Востока

 

 

  • Oceanian — модальный компонент аборигенных жителей Океании, Австронезии, Меланезии и Микронезии — пик у современных папуасов и австралийских аборигенов

 

 

  • Paleo-African — модальный компонент африканских пигмеев и бушменов

 

 

  • Siberian — модальный компонент народностей юго-восточной Сибири

 

 

  • Subsaharian — второй африканских компонент — пик в популяциях мандинка, йоруба и ишан

 

 

  • WHG-UHG — компонент древних европейских мезолитических охотников-собирателей, пик в образцах палеогеномов мезолитических популяций европейских охотников-собирателей. Тождественен аналогичному компоненту в работах популяционных генетиков (Lazaridis et al. 2014, Haak et al. 2015). Из современных популяций — наивысший процент в популяциях эстонцев, литовцев, финнов и др.

 

 

MDS plot - K13 ancestral population

Как я и предполагал, модель калькулятора оказалась особенно хороша в применении к анализу древних геномов. И на самом деле, на нижеприведенном графике PCA (пространстве 2 главных компонент результатов анализа древних геномов в моем бета-калькуляторе K13) видны замечательные вещи. Расположение геномов хорошо вписывается в треугольник, один из углов которого образуют геномы древних «ямников» (из работы Haak et al. 2015), причем геномы «русских»ямники из работы Allentoft at al.2015 чуть-чуть сдвинуты в сторону древних мезолитических геномов древних европейских охотников-собирателей. За ними (в направлении «неолитического» угла) следуют представители шнуровой культуры, еще дальше — геномы представителей геномов унетицкой культуры и т.д. Второй угол треугольника образован неолитическим геномами, причем если более поздние неолитические геномы сдвигаются ближе к представителям линейно-ленточной культуры (англ. Linear Pottery culture, фр. Culture rubanée, нем. Linearbandkeramische Kultur, LBK — наиболее распространенная неолитическая культура Центральной Европы 5500—4500 гг. до н. э.), то более ранние геномы — геном представительницы более ранней фазы этой культуры (Stuttgart-LBK), а также геномы представителей балканских неолитических культур — Старчево и Винча — очень близки к палеогеному из Barcin (культура Чатал-Хююк, cамые ранние найденные культурные слои относятся к 7400 г. до н. э.). Таким образом генетика подтверждает утверждения археологов о близости неолитических культур Балкан и Анатолии. Более того — данные генетики свидетельствуют о том, что во времена т.н «неолитической революции» происходила не только и не столько миграция технологий (как считали некоторые археологи), но и миграция населения (из Анатолии на Балканы). Причем, судя по моему графику PCA, миграция происходила в несколько, хронологически удаленных, этапов, и — скорее всего — из разных мест. Крайную точку в этому угле треугольника я обозначил как «анатолийские земледельцы» (ближайший к этой точки геном — геном «земледельца» из культуры Старчево — взят из работ Haak et al. 2015).

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

Образно говоря, вектор градиента начинается в геномах ямников (больше половины генома которых состояла из этого компонента) и затем идет к геномам представителей синташтинской, афанасьевской, андроновской, окуневской и карасукской культур.
Пару слов об этих культурах (положение геномов представителей которых можно посмотреть на графике).
1) Синташтинская культура формировалась из древнеямных и катакомбных племён и местного населения. Синташтинцев связывают с индоиранскими племенами.
2) Андроновская культура также развивается на базе ямной. На западе она доходила до района Урала и Волги, где контактировала со срубной культурой. На востоке андроновская культура распространилась до Минусинской котловины, частично включив в себя территорию ранней афанасьевской культуры. Андроновцев (также как и синаштинцев) относят к индоиранской сообщности.
3) Афанасьевская культура была создана мигрантами из Восточной Европы, в частности, носителями древнеямной культуры, ассимилировавшими местное население. Сменилась карасукской и окуневской культурами.Наследниками афанасьевцев были племена тагарской культуры, дожившей до III в. до н. э., по другой версии, тагарцы были скифами, а потомки афанасьевцев — тохарами, которых именно скифы-тагарцы вытеснили в Синьцзян.Большинство исследователей ассоциируют афанасьевскую культуру с (прото-)тохарами.
4) Окуневская культура — фнтропологический тип населения этой эпохи был смешанного европеоидно-монголоидного происхождения, с преобладанием монголоидного. Как отмечает А. В. Громов, бросается в глаза их морфологическая разнородность — встречаются как чисто монголоидные черепа, так и типично европеоидные, не обнаруживающими никаких следов монголоидной примеси. Проведя обстоятельный анализ антропологических особенностей населения неолита и ранней бронзы, А. А. Громов пришел к выводу, что физический тип окуневцев сложился в результате смешения местного неолитического населения с выходцами из территории Средней Азии и Казахстана (афанасьевцами)
5) Карасукская культура — развилась на основе окуневской культуры под влиянием андроновской культуры.

Интересно, что геном мальчика с южносибирской палеолитической стоянки MA-1 как раз проецируется между центроидами геномов представителей синташтинской, афанасьевской, андроновской, окуневской и карасукской культур. Эти геномы (вернее их центроиды) занимают на графике значительное место. Самый дальний из них — геном алтайца из эпохи железного века (примерно 50 год до нашей эры). Сразу за ним идут все из имеющихся у меня палеогеномов жителей Америков (палеоэскимосы — в том числе и Saqqaq; и «палеоиндейцы» — Clovis, древние жители Перу и палеогеномы Botocudo). Любопытно что последние — геномы Botocudo — хотя и являются самыми современными (1600 год нашей эры), однако в них хорошо заметен «океанский компонент», именно поэтому они смыкаются на графике с палеогеномом австралийского аборигена. В этой связи я вспоминаю оригинальную теорию Тура Хейердала о наличии доисторических контактов между жителями островов Тихого океана и жителями Южной Америки.

Особое место на графики занимают «живые реликты» — онге, один из коренных андаманских народов (адиваси), геномы так называемого «усть-ишимца» (возраст 45000 лет), костенковца (Kostenki-14, возраст 38 700 -36 200 лет), и недавно опубликованный палеогеном Oase из Румынии (возрастом 37000-42000 лет). Они образуют отдельную группу (особенно близки друг к другу румынский палеогеном Oase и усть-ишимец), однако я терясь в догадках о том, что именно означает столь заметная близость этих геномов.

 

 

Eurasian and American paleogenomes

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

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

Ранее эта тема поднималась в описании программы SupportMix, а также в сжатом изложении методологии оценки происхождения хромосомных сегментов (инструмент PCAdmix).  Данная заметка будет посвящена третьему инструменту — LAMP (Local Ancestry in adMixed Populations) (Sankararaman et al.2008).

Очевидно, что алгоритмы определения локального происхождения отдельных сегментов человеческих хромосом могут дать неплохие результаты при комбинированном использовании программ PLINK /ADMIXTUIRE/LAMP: например, комбинация этих программ позволяет довольно точно определить не только стратификацию отдельных этно-популяционных групп,  но также и уровень «адмикса» у отдельных людей. Поскольку одна из задач нашего проекта MDLP состояла в определении практических и теоретических преимуществ и/или ограничений конкретных методологий биоинформатического анализа полных генома, я провел эксперимент, позволяющий прояснить ряд ограничений, которые значительно уменьшают уровень достоверности результатов  субструктуры аутосомного генофонда населения Европы.

В качестве инструмента контроля качества комбинированного набора данных (аутосомных SNP-ов 22 хромосом) я использовал Plink, с помощью которого я выбрал для последующего анализа только качественные снипы (99% генотиприрования),  частоты минорных аллелей которых превышают 1%.

Поскольку этно-популяционный фон неравновесного линикиджа марекеров (LD) может существенным образом влиять на основные компоненты субструктуры популяции, я исключил из выборки маркеры, характеризующиеся статистически значимым уровнем LD (с коэффициентом попарной корреляции r2 Пирсона > 0,4) в «скользящем окне» из 100 снипов  с пошаговым сдвигом на 10 снипов. Кроме этого, я также использовал  другие методы Plink для получения однородной выборки  — например, кластеризации на основе IBS для обнаружения пары индивидов (outliers) с  уровнем «родства», значительно более высоким, чем у пары выбранных случайным образом индивидов в однородной популяции.  Под более высоким родством здесь понимается  резко отклоняющиеся значения (более 3 стандартных отклонений) парных значнений IBS по отношению к остальной части выбаки, а также случаи с высоким значения PIHAT (более 0,05) и  высокой степень инбридинга (гомозиготности*). Индивиды с подобными аномальными значениями («выбросы») были удалены из  «обучающего» подмножества нашей выборки .


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

 

homozyg
Стратификация образцев в соответствии с уровнями гомозиготности. Вдоль оси Х отображена общая сумма гомозиготных сегментов в килобазах; вдоль Y-оси — средний размер гомозиготных сегментов в килобазах

 

 

homozyg2
Уровни индивидуальной гомозиготности в выбороке: вдоль ости X отложено количество сегментов NSEG. Общая длина гомозиготных сегментов отображается осью Y

 

По окончанию описанных выше процедур фильтрации снипов и удаления «выбросов», окончательный набор данных представлял собой набор данных из 90 455 снипов и 317 человек (289 мужчин, 82 женщин). Эти данные были использованы в последующем анализе.

Прежде всего, мы использовали программу ADMIXTURE (Alexandre, Novembre, Lange 2009), в которой реализована модель оценки максимального правдоподобия (ML), т.е алгоритм кластеризации и оценки структуры популяции в наборе генетических данных (снипов).

В целях сохранения совместимости с MDLP калькулятором, я остановился  на модели, в которой выборка представлена в виде комбинации 7 предковых компонентов (K=7).  Индивидуальные значения процентной составляющей каждого компонента в индивидуальном геноме (матрица Q), была визуализированы в R (ниже приведен график с результатами участников проекта MDLP, полный список  доступен в этой таблице).

Результаты K=7

Полученные предковые компоненты (K=7) я обозначил следующими названиями (с сопутствующей цветовой легендой)**:

  • Транс-кавказский — красный
  • Балканском / средиземноморском -желтый
  • Северо-кавказский -зеленый
  • Западно-европейский
  • Алтайский — светло-голубой
  • Балто-славянский — темно-синий
  • Прибалтийско-финский / Северо-европейский -фиолетовый

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

 

 

 

 

MDLP v4 components

 

 

 

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

Я предположил, что все образцы в моей выборке (представленной образцами VID)  проекта MDLP возникли в в результате смешивание 7 отдельных предковых групп населения. Данное предположение означает, что «чистые» референсные группы населения тесно связаны с истинными предковыми популяциями. Исходя из этого предположения мы снова задействовали программное обеспечение Admixture,  на этот раз с целью определения предковых компонентов в фазированном наборе данных из отдельных неполовых (аутосомных) хромосом.

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

Я  расчитал в программном обеспечении Plink частоты аллелей (в стратифицированных по этническим признакам кластерах), и добавил в файл фиксированные частоты рекомбинации (определяются отдельно для каждой из 22 хромосом). Для моделирования динамического процесса смешивания предковых компонентов, я использовал различное количество поколений G ( 5, 10,25 поколений),  предполагая 3 хронологически разных варианта, в которых при  K = 7  предковые популяции A1, …, Ak,  перемешивались в течение G = 5,10,25 поколений.

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

1) результаты Admixture для фазированных генотипов хромосомы (Chr * -phased)
2) результаты Admixture для нефазированных генотипов хромосомы (Chr * -unphased)
3) результаты LAMP для G = 5 (Chr * -lamp-GEN5)
4) результаты LAMP для G = 10 (Chr * -lamp-GEN5)
5) результаты LAMP для G = 25 (Chr * -lamp-GEN5)

Образец этих выходных данных можно посмотреть в файле Excel с результатами анализа хромосомы 1 (Chr1).

PCAdmix: инструмент и методология для оценки происхождения хромосомных сегментов

В марте прошлого года  Сергей Козлов — один из соавторов данного блога, — опубликовал важную с точки зрения методологии генетико-генеалогического анализа заметку о принципах оценки вероятности определения времени жизни последнего общего предка при попарном сравнении аутосомных данных двух или более сравниваемых индивидов.  Действительно, в последние годы среди людей, интересующихся генеалогией, приобрели заметную популярность сервисы, производящие поиск генетических родственников по всем линиям, а не только по прямой мужской и прямой женской. В качестве примера можно привести Family Finder от FTDNA и DNA relatives от 23andMe. Участник получает достаточно длинный список так называемых «совпаденцев» — людей, имеющих с ним один или более участок половинного совпадения (УПС) на аутосомах (неполовых хромосомах). Если участок достаточно длинный (а его длина измеряется в сантиморганидах, обозначающих вероятность разрыва участка при каждой передаче в следующее поколение), то это говорит о наличии общего предка (от которого участок и получен).
Для значительной части клиентов сервисов персональной коммерческой геномики, интересующихся исключительно вопросами своего происхождения, вопрос о достоверном определении времени жизни общих предков имеет первостепенное значение. И вместе с тем, именно проблема с получением четкого ответа на этот краеугольный вопрос служит одной из главных причин недовольства и раздражения клиентов компаний вроде FTDNA или 23andme.

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

Однако наряду с  географически близкими (в пределях 50-100 км)  «совпаденцами», нередко в списках «совпаденцев», предоставляемых в 23andme или FTDNA появляются совершенно экзотические «совпаденцы». Например, у финна может появится совпаденец из Италии, а у корейца — из  Великобритании. Совершенно очевидно, что подобные случаи очень сложно объяснить не только простым сопоставлением сведений о географическом происхождении предков, но даже и безотказной в простых случаях  моделью наложения «этнопопуляционного аутосомного фона в виде коротких реликтовых  IBD сегментов».

В этой связи возникает практический вопрос — как интерпретировать подобные случаи, при условии что подобные сегменты представляют собой не «ложно-позитивные», а вполне достоверные совпадения, указываюшие на существование в неопределенный момент прошлого некоего общего предка. И подобные случаи характерны не только для коммерческих «выборок», но и для вполне серьезных научных баз данных, например 1000 Genomes. В частности,  в этой базе данных при сравнении редких снипов у 89 британцев и 97 китайцев были обнаружены три англо-китайские пары с отдаленным генеалогическим родством ( в геноме этих пар были обнаружены идентичные по происхождению фрагменты (IBD сегменты) ДНК,  которые составляют 0,001%, 0,004% и 0,01%  их геномов).

Самое простое решение этой проблемы некоторые из любителей генетической генеалогии пытались найти в обращении к сервисам главного инструмента аутосомной генетической генеалогии  Gedmatch. В частности, как известно, данный сервер содержит онлайн-версии практически всех популярных среди любителей модификаций DIYDodecad калькуляторов. Например, выбрав разработанный мною калькулятор MDLP K23b в режиме Chromosome painting: Paint differences between 2 kits, 1 chromosome   и сравнив характер распределения предковых компонентов на гомологичных хромосомах у двух сравниваемых людей, можно получить примерное представление о географическом ареале, в котором мог жить общий предок этих людей (вероятно, на этот ареал будет указывать доминирующий на совпадающем сегменте компонент). Логика простая. Предположим, например, что мы сравниваем  сегменты хромосомы X в данных индивида A этнического происхождения D c данными индивида В этнического происхождения С. Здесь возможны три варианта

  • С-происхождение предка или предков индивида A
  • D-происхождение предка или предков индивида B
  • Y-происхождение подмножества предков обоих индивидов

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

 


Пример I.

В этом примере мы будем использовать свои данные и данные женщины, с которой у нас был обнаружен подтвержденный генеалогией общий предок, живший в середине 19 века.  При сравнении наших данных, алгоритм поиска достоверных генеалого-генетических совпадений обнаружил три сегмента с генетической дистанцией > 7 cантиморганов, cостоящих в блочной записи из более чем 700 последовательно совпадающих снипов

Start Location End Location Centimorgans (cM) SNPs
4 32232224 42421625 13.2 1115
7 8295405 13845989 9.8 885
11 36784445 45084878 8.0 881

Самый большой сегмент = 13.2 cM
Общий размер сегментов с сантиморганах > 7 cM = 30.9 cM
Приблизительное число поколений до общего предка  = 4.4

Задетектированные  сегменты хромосом идеографически отображаются при попарном сравнении в цветовой гамме — черный цвет означает несовпадающие сегменты, другие цвета — компонентную привязку к одному из компонентов моего калькулятора MDLP K23b.  Ниже приведены фрагменты идеографического отображения 2 из 3 вышеуказанных совпадающих сегментов на кариограмму 4 и 7 хромосомы.:

M051225_F298455_4_D64088
Сегмент на 4 хромосоме
M051225_F298455_7_BC1A38
Сегмент на 7 хромосоме

Самый значительный сегмент (13.2 сM) на 4 хромосоме имеют хорошо заметную привязку к северо-восточно-европейскому компоненту [зеленый цвет], в исторической перспективе связанному с наследием мезолитического населения этого региона. А вот сегмент на 7 хромосоме имеет более сложную структуру, в которой характерно преобладание кавказского компонента [голубой цвет]. Таким образом можно уверено утверждать, что общий предок (или предки) могли жить в регионе восточной Европы.

К сожалению, данный инструмент сегментного сравнения на  Gedmatch хотя и прост в обращении (в силу интуитивной понятности), однако  далек от совершенства. В первую очередь, на аккуратность определения «генографического»происхождения сегмента влияет отсутствие на сервере  гаплоидных фаз похромосомных данных. В результате, сравнение ведется не по конкретной фазе (т.е по конкретной хромосоме доставшейся ребенку от каждого из родителей), а по диплоидному составному блоку, т.е вместо настоящих IBD мы можем оперировать half-IBD (HBD), которые на слэнге русскоязычных любителей именуются УПС-ами. Во вторых, аккуратность генографического определения  зависит от аккуратности определения предковых компонентов в используемом варианте калькулятора, но это отдельная тема для разговора.


К счастью, парадокс «экзотических» совпаденцев имеет более точное решение с помощью одной из программ, позволяющих определять геногеографическую структуру или «локальное происхождение» совпадающих сегментов.  Можно использовать разные программы, HAPMIX, LAMP , HAPAA, ANCESTRYMAP — так как несмотря на ряд принципиальных отличий, все они используют алгоритмы моделнй скрытых марковских цепей (HMM) и поэтому выдают в целом схожие результаты. К этому же классу программ относится и более новая програма PCAdmix, которую я буду использовать в своем втором примере, в котором я задействую фазированные в BEAGLE генотипы.  В целях разжевывания принципов работы программы, следует вкратце описать рабочий процесс PCAdmix.
PCAdmix являет cобой метод, который оценивает локальное происхождение хромосомных сегментов с помощью анализа главных компонентов (PCA)  фазированных гаплотипов. В самом начале выполняется анализ главных компонентов в 2-3 референсных панелех, необходимых доя построения пространства главных компонентов, например, для хромосомы 22 . Поскольку метод использует фазированные данные, каждая копия хромосомы 22 в референсных панелях рассматривается как отдельная точка в пространстве главных компонентов. Первые две главные компоненты, как правило, представляют собой оси «предкового» расхождения популяций референсных панелей, что хорошо заметно на графиках. Если подобного рассхождения не наблюдается,  то скорее всего в популяциях референсных панелей «маскируется» присутствие неявной популяционной субструктуры. В построенное таким способом пространство главных компонентов в дальнейшем проецируется группа лиц «смешанного» происхождения, и затем определяется значение нагрузки главных компонентов для каждого снипа.  После этого метод переходит к анализу коротких «окон» снипов — для каждого из этих окон вычисляются  вероятности того, что данное окно в гаплотипе человека «смешанного» происхождения происходит от одной из референсных популяций. Вычисленные таким образоом вероятности различных вариантов происхождения каждого окна снипов, используются на заключительном этапе метода в  скрытой моделе Маркова (HММ) для сглаживания шума в определении происхождения «окон» снипов. Таким образом, данная скрытая модель Маркова НММ зависит от значений главных компонентов, доли каждого «компонента происхождения» на заданной хромосоме, а также матрицы перехода, которая, в свою очередь, зависит от числа поколений прошедших с момента смешивания популяций и генетического расстояния (сM) между двумя окнами снипов. В текущей версии метода, рекомбинаторные расстояния и число поколений определяются параметрами.
Конечным результатом рабочего процесса PCAdmix является матрица состяний скрытой модели Маркова, содержащая апостериорную вероятность каждого из возможных вариантов происхождения для данного «окна снипов», и эта вероятность обусловлена остальной частью данных для хромосомы. Важно отметить, что происхождение каждого окна снипов определяется только в том случае если апостериорная вероятность для одного из возможных происхождений > = 0,8. Любое окно, для которого максимальная апостериорная вероятность любого варианта происхождения <0,8, считается «неопределенным».


Пример 2

Данный пример основан на реальном случае, когда ко мне обратился человек, чьи предки происходят из центральных регионов Азии. Смущенный наличием в списке своих совпаденцев в сервисе Relative Finder 23andme  человека с корейскими и японскими корнями, а также  семейными легендами о «восточноазиатской»прабабушке, он попросил меня определить вероятность присутствия японцев в числе своих ближайших (в пределах 5 поколений) предков, опираясь исключительно на аутосомные данные.

В этом эксперименте, я решил скурпулезно следовать инструкциям разработчиков PCAdmix, и для начала произвел фазирование (биоинформатическую реконструкцию гаплотипных фаз аутосомных хромосом) в программе BEAGLE. Данные тестанта (ок 400 тыс. снипов) были фазированы в присутствии 3 контрольных референсных групп популяций — британцев GBR, китайцев CHB и японцев JPT — поскольку эти группы были позднее задействованы мной в качестве 3 референсных панелей. В целях уменьшения количества ошибок, которые неизбежно появляются в результате импутации пропущенных «генотипов» снипов, я использовал только те общие снипы, которые были определены как в аутосомных данных клиента 23andme, так и в трех референсных группах.

Затем фазированные данные тестанта были похромосомно обработаны в рабочих циклах программы PCAdmix. Программа отфильтровала cнипы с низким значением MAF и высоким значением LD, в результате чего число снипов уменьшилось почти вдвое. Оставшиеся снипы были разбиты на «окна снипов», каждое из которых состяло из 20 снипов.  При расчете по всем 22 хромосомах, общее количество полученных таким разбиением «окон» составило 11 997. В конце рабочего цикла (метод главных компонентов + HMM) программа выдала для каждой парной аутосомной хромосомы A и B  файл в формате bed, удобном для отображения дополнительной информации в аннотации генома (номер хромосомы, начало и конец сегмента, наиболее вероятный регион происхождения сегмента, cM, максимальная вероятность и апостериорная вероятность одного из трех вариантов происхождения — JPT, GBR, CHB, непоказана в таблице). В конечном отчете GBR используется как индикатор сегментов не-восточноазиатского происхождения (nEA), JPT — японского происхождения (JPA), CHB — неспецифичных сегментов восточноазиатского происхождения (EA) :

10 111955 468599 GBR 0.004885 0.134147 GBR* 0.636943
10 521723 811876 GBR 0.142147 0.582463 GBR* 0.646868
10 815149 1151723 GBR 0.585829 0.898724 GBR* 0.676252
10 1156487 1335849 GBR 0.901503 1.23673 GBR 0.925059
10 1337709 1449849 GBR 1.24246 1.60705 GBR 0.99999
10 1454864 1510208 GBR 1.61249 1.76798 GBR 0.999506
10 1512546 1623734 GBR 1.77039 2.12653 GBR 0.999647
10 1624900 1669347 GBR 2.13038 2.25357 GBR 0.999778


Выбор формата BED в качестве формата выходных в моем случае также был далеко неслучайным. C помощью одной из библиотеки платформы Bioconductor формат BED легко отображается в кариограмме 22 пар аутосомных хромосом человека (я использовал координаты геномного билда b37). Чтобы было понятно, что именно изображают эти «кариоплоты» (идеографические изображения хромосом), необходимо пояснить, что  «японское происхождение» (JPA) приписывалась 20-сниповому сегменту только в том случае, если апостериорная вероятность японского происхождения данного «окна из 20 снипов» составляла > = 0,8. Любое окно, для которого максимальная апостериорная вероятность любого варианта составляля <0,8, засчитывалось как окно  с «неопределенным» происхождением (UND).Chromosomes A

Chromosomes A

 

Chromosomes B
Chromosomes B

Эксперимент показал, что среди 11997 «окн» число  «окон» не-восточноазиатского (nEA) происхождения (7650) почти в два раза больше чем число «восточноазиатских» сегментов. Происхождение 2750 геномных «окон» снипов невозможно определеить, и только 965 «окна» могут быть определены как «японские по происхождению». Вместе с 617 окнами «китайского» (EA),  восточно-азиатские сегменты составляют меньше, чем 10% генома.
Не менее важно и то обстоятельства, что значительная доля этих сегментов-окон пришлась на низких «консервативные, низкорекомбинантные» области хромосом,  — такие, как  например, теломеры, центромеры и регионы с низкой плотностью снипов: сегменты в таких регионах могут переходить от одного поколения к другому фактически в неизменном виде. Наконец, те же закономерности распределения родословной были отмечены в обеих фазированных наборах аутосомных хромосом, что опровергает версию о недавной «восточноазиатской» примеси со стороны одного из родитедей и скорее  свидетельствует о древнем эпизоде смешивание определенных центрально- и юго-западноазиатских групп с группами восточноазиатского происхождения (например, в ходе монгольских или тюркских нашествий).

Разумеется, как и во многих других моделях анализа, основанных на вероятностях, наше заключение нельзя считать окончательным вердикторм. Вместо этого, лучше сказать, что шансы в пользу существования «недавнего японского предка» против шансов отсутствия такого, составляют 10 к 90. Другими словами, вариант с недавней японской «примесью» нельзя полностью исключить, поскольку вероятность такого сценария  составляет 11%.

 

Расширенные карты для палеогеномов

Обновлено 30.11.2014

Этот пост также продолжает один из предыдущих, а именно визуализацию суммы IBD-сегментов (а возможно, это и не IBD — вопрос остается открытым) двух палеоевропейцев и мальчика со стоянки Мальта с современными  выборками. С тех пор в открытом доступе появилось еще несколько обработанных палеогеномов — «усть-ишимец«, «Костенки-14» («человек с Маркиной горы») и два палеогенома хорошего качества из Венгрии.

Результаты собраны мной в онлайн-таблицу, а также отрисованы на расширенных картах. Поскольку усть-ишимец явно тяготел к восточноазиатам, пришлось добавить к сравнению выборки из Южной и Восточной Азии. Ну а после этого логика подсказывала, что неолитических земледельцев Европы неплохо бы сравнить с жителями Ближнего Востока. Таким образом, все карты перерисованы.

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

«Неолитическая фермерша» )) Stuttgart/LBK

«Неолитический земледелец» NE1:

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

«Охотник-собиратель» Loschbour:

Разница между «охотником-собирателем» и усреднением по двум земледельцам. Красное — больше сегментов с Loschbour, зеленое — c Stuttgart и NE1

Европеец позднего бронзового века BR2 из Венгрии:

«Усть-ишимец»:

Костенки-14 (ослабленные настройки):

Мальтинец (аналогично):

И наконец, в качестве примера результата нашего современника, моя собственная карта:

 

 

Методологические заметки к созданию неандертальского калькулятора

Данная заметка представляет собой критический анализ методологических предпосылок создания неандертальского калькулятора, имплементированного в cоответствующем сервисе 23andme (Neanderthal lab). В основу заметки положен перевод технического документа 23andme (white paper), описывающего процесс создания неандертальского калькулятора.

Существует несколько методологических подходов  к созданию неандертальского калькулятора (т.е инструмента для оценки того, сколько процентов ДНК  в геноме анализируемого индивида имеет  неандертальское происхождение).   Есть несколько способов прямой экспериментальной оценки величины процента «неандертальской » ДНК с помощью ресеквенирования ДНК клиента в тех регионах, в которых ученые обнаружили возможные варианты, имеющие предполагаемое неандертальское происхождение. Но в силу технической сложности реализации этих способов и необходимости каждый раз заново производить секвенирование в полном объеме  регионов неандертальского происхождения, нет особой нужды рассматривать их в этой записи. Вместо этого я предлагаю рассмотреть две оставшиеся методики определения вклада неандертальского ДНК.  Хотя оба метода не без своих изъянов, они позволяют существенно снизить влияние неопределенности (ascertainment bias) в оценке вклада неандертальского ДНК, и в принципе,  других приемлемых альтернатив этим методам не существует, так как в противном случае получаемый другими методами (например, Dstatisticsили ABBABABA) разброс оценки величины неандертальского вклада будет в несколько раз отличаться от тех величин, которые получаются на выходе соответствующих программ, используемых в  NationalGeographicGeno и 23andme (обе программы основаны на одном из двух нижеописанных методов).Именно по этой причине, каждая из нижеприведенных методик заслуживает отдельного рассмотрения. 

  1. Метод PCA

На мой личный взгляд, наилучшим  (как в плане аккуратности, так и в плане легкости реализации) методом оценки величины неандертальца в ДНК клиентов является метод главных компонент PCA, так как он представляет собой очень мощный инструмент для представления корреляции данных высокой размерности (порядка миллионов снипов и даже больше) в виде гораздо меньшего, некоррелирующего набора переменных, которые носят название «главные компоненты». Итак, метод главных компонент — это один из способов понижения размерности, состоящий в переходе к новому ортогональному базису, оси которого ориентированы по направлениям максимальной дисперсии набора входных данных (в нашем случае это набор генотипов снипов). Вдоль первой оси нового базиса дисперсия максимальна, вторая ось максимизирует дисперсию при условии ортогональности первой оси, и т.д., последняя ось имеет минимальную дисперсию из всех возможных. Такое преобразование позволяет понижать информацию путем отбрасывания координат, соответствующих направлениям с минимальной дисперсией. Можно отметить, что в основе метода главных компонент лежат следующие допущения: (a) допущение о том, что размерность данных может быть эффективно понижена путем линейного преобразования, и  (b)  допущение о том, что больше всего информации несут те направления, в которых дисперсия входных данных максимальна.

 

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

Сначала скачиваем полные геномы неандертальца, денисовского человека, и шимпанзе Clint. Затем с помощью программы samtools генерируем для каждого из трех геномов файлы с геномными вариантами (vcf), отфильтровываем из полученных файлы инделы, таким образом чтобы на выходе остались только снипы и проводим аннотацию  снипов с использованием базы данных dbSNP; при аннотации находятся те варианты, которые присутствуют в базе данных и им назначается соответствующий индекс, например rs4213456 (это условный пример). Затем необходимо выбрать из это файла только те cнипы, которые присутствуют в контрольной выборке с референсными популяциями современного человека. Описание примерного порядока выполнения этой задачи можно найти в двух записях в моем блоге (здесь и здесь).

В конечном итоге, по окончанию первого этапа,  мы получаем три файла VCF c аннотированным снипами, которые необходимо соединить в один файл либо в vcftools, либо в Plink. Затем провести анализ PCA с двумя заданными главными компонентами (K2) в самом Plink, либо конвертировать данные в формат Eigenstrat и провести в программе Eigensoft анализ PCA (также с двумя заданными главными компонентами). Последний вариант предпочтителен, так алгоритм Eigensoftдает более точные данные за счет kernel-преобразований данных. В конечном результате проведенного анализа двух основных компонентов должны получится нормированный лист cобственных векторов — эйгенвекторов так называемый лист факторной загрузки –factor loading) для каждого из индивидуальных образцов, входящих в анализируемый набор. Первый главный компонент, PC1 , чьи значения отображаются вдоль первой оси ортогонального  базиса, характеризуется максимальной дисперсией набора снипов входящих данных, эта ось отображает общее генетическое сходство архаичных людей (неандертальца и денисовского человека). Ось второго компонента , PC2 , оптимизирует дисперсию при условии ортогональности первой оси (т.е, PC1), и  отображает генетическое расхождение между неандертальцами и денисовским человеком. 

pca

 

На следующем этапе генотипы клиентыпроецируются на плоскость, образованную двум яосями PC1 и PC2.  Я полагаю, что на этом этапе в самом PCA анализе нет необходимости, вместо этого можно имплементировать метод с использованием высчитанного в первом анализе PCA листа загрузки компонентов (loadings). Подобный подход реализован, например, в программе shellfish. 

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

Как видно из иллюстрации к работе (Reich et al.2011), европейцы и жители Восточной Азии существенно сдвинуты в сторону неандертальцев по сравнению с афро-американцами (как видно из приведенного ниже графика,  расстояние между неандертальским «углом» и положением афроамериканцеввесьма значительно, это следствие неопределенности определения предковых аллелей неандертальца по африканским популяциям, поэтому для коррекции этой дистанции в 23andme высчитали центроид генетического положения африканцев с использование данных проекта 1000G, и расчет дистанции вели от него).

reich

 

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

  1. Методтеговых (маркерных) снипов— NAIM (Neanderthal Ancestry Informative Markers)

Существует более прямой и простой способ  вычисления неандертальского вклада в геном клиентов. Простота метода обусловлена отсутствием надобности в сравнительно сложных алгоритмах вычисления главных компонентов. Согласно известной публикации драфтовой версии генома неандертальца (Green et al., 2010), в геномах современных людей были обнаружены 13 геномных регионов, которые, как предполагают авторы, имели неандертальское происхождение.  Эти регионы генома  современных людей  были маркированы с помощью маркерных (теговых) снипов – то есть таких снипов, в которых неандертальский вариант часто встречается в современных неафриканских популяциях людей, но отсутствует в коренных африканских популяциях.

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

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

  1. Во-первых, не существует никаких формальных гарантий того, что эти варианты действительно  имеют неандертальское происхождение.
  2. Во-вторых, даже в том идеальном случае, когда все эти 180 вариантов действительно имеют неандертальское происхождение, они охватывают только 13 геномных регионов, самый длинный из которых представляет собой сегмент длиной всего лишь в 160 000 базовых пар. Эта длина на два порядка величин ниже, чем среднестатистические 2,5% неандертальского вклада в среднестатистическом геноме современного человека неафриканского происхождения . Поэтому простой подсчет числа неандертальских вариантов в маркерных снипах, где встречается будет в 2-3  раза занижать реальный процент неандертальского вклада в клиентском геноме.
  3. В-третьих, существует еще несколько трудных моментов, связанных с практической реализацией этого метода.

3.1.     Списка вышеупомянутых 180 снипов нет в открытом доступе, и так как в оригинальной статье было упомянуто другое количество снипов (166), похоже на то, что это число снипов варьируется в зависимости от использованного чипсета (поэтому и число снипов разное).

3.2.     Технически  эту проблему можно решить следующим образом. Самый простой способ состоит в определении того, какие снипы из используемого компанией чипсета попадают в эти сегменты. Например, берется первый сегмент на хромосоме 1 (начало 168 110 000 – конец 168 220 000, длина в базовых парах – 110 000) и выбираются снипы попадающие в этот регион, и так далее по всем регионам. При этом сначала надо узнать какой билд используется в контрольной выборке популяций современных людей. Если используется build 37, тогда необходимо конвертировать координаты сегментов в более ранний build 36. После того, как будут определены все снипы попадающие в эти 13 сегментов, нужно найти неандертальские варианты этих файлов (это можно сделать в базе данных неандертальских снипов) и составить список, который затем использовать в качестве затравки при сравнении с значениями снипов у современных людей.

3.3.     Другой вариант более сложный, но очевидно более точный. Список снипов найденных в ходе сравнения геномов шимпанзе, 5 референсных популяций современных людей и неандертальца  выгружен на сайте геномного браузера UCSC. Это большой файл (в распакованном виде 363 Mb), общее количество снипов 5 615 438. Формат файла следующий:

971    chr1       50600811             50600812             AA_AAD:0D,1A  0             +             50600811             50600812             0

971    chr1       50603655             50603656             AAD_AA:0D,2A    0             +             50603655             50603656             0

971    chr1       50604033             50604034             AADAA_:0D,1A    0             +             50604033             50604034             0

971    chr1       50605949             50605950             AAA_DA:0D,1A    0             +             50605949             50605950             0
Первая колонка представляет собой номер сегмента чтения, вторая – название хромосомы, вторая и третья – физическое положение снипа, далее идет длинная колонка с указанием характера варианта в  шимпанзе, 4 популяций людей и неандертальца. «A» обозначает предковое значение аллеля, «D» — derived, т.е мутировавшее значение. После двоеточия идет специфическая неандертальская колонка (например, :0D,1A)с указанием того сколько предковых и сколько мутировавших значений снипа обнаружено в исследованных геномах неандертальцев. В данном случае, в первом снипе обнаружено 0D (0 мутировавших) и 1A (1 предковое значение). Трудность задачи состоит в определении только тех снипов, в которых  у неандертальцев нет предковых значений, а встречаются только мутировавшие значения. Эти снипы — кандидаты на неандертальский вклад в человеческий геном. Затем сравнить отфильтрованный список со списком снипов в  контрольной выборке (опять-таки, надо знать какой билд используется, координаты этого списока  приведен по билду 36) и выбрать только те, что имеются в чипсете компании. Далее алгоритм тот же, что и выше – определяется значение снипа у неандертальца и  сравнивается с соответствующим значением у современных людей. Совпадающие у неандертальца и современных людей варианты подсчитываются и определяется конечный процент неандертальского вклада.

 

Эксперимент.

 

Я решил проверить эфективность первого метода (метода PCA) на своей контрольной выборке (2778 образцов современных людей, шимпанзе, денисовского человека и неандертальского человека и 142429 снипа). В качестве рабочей программы я использовал новую версию Plink, которая позволяет использовать в анализе PCA заданные контрольные кластеры, в которые проецируются исследуемые индивиды. В качестве трех контрольных групп я выбрал, следуя рекомендациям авторов обсуждаемого исследования,  геномы шимпанзе, неандертальца из Vindja и денисовского человека. Однако число априорных главных компонентов я намерено изменил,  с 2 на 3 (K3), таким образом на выходе я получил эйгенвекторы трех главных компонентов.  По этой причине, полученный мной график PCA несколько отличается от вышеприведенного графика 23andme (вместо PC1 и PC2 я использовал PC2 и PC3, то есть второй и третьи главные компоненты, более точно описывающие в данном случае сходство/различие геномов архаичных и современных людей).

R Graphics Output
Как видно из наших результатов, все популяции современных людей разместились внутри условного треугольника образованного дисперсией геномов денисовского человека, неандертальца и шимпанзе.
Впрочем, на графике нельзя разглядеть, какие именно популяции сдвигаются в сторону неандертальца, а какие — в сторону денисовского человека (такой сдвиг свидетельствовал бы о наличии адмикса).  Чтобы устранить этот досадный артефакт графика, придется убрать с графика геномы денисовца, неандертальца и шимпанзе (из-за значительной генетической дистанции популяции современных людей сдвигаются в одну кучу).

 

R Graphics Output
R Graphics Output

 

Положение удаленных денисовца, неандертальца и шимпанзе размечено на новом графике буквенными обозначениями — D, N, Chimp. Из человеческих популяций я разметил группы африканских популяций (Africans), и коренных американцев (Native Americans). Европейские и азиатские популяций смещены в одну общую группу, с сильным креном в сторону неандертальца. Судя по всему, мои результаты, в общих чертах, практически не отличаются от результатов исследований Грина и Райха. Как отмечает  Дробышевский: » «денисовские гены», несмотря на свою экзотичность, обнаружились у современных людей. Первоначально они были найдены у папуасов Новой Гвинеи и меланезийцев острова Бугенвиль (Reich et al., 2010), затем – у австралийских аборигенов (Gibbons, 2011), а полнейшее исследование вопроса констатировало наличие их у огромного числа популяций (Reich et al., 2011). Они были выявлены в тридцати трёх популяциях Океании и Юго-Восточной Азии, в том числе у папуасов Новой Гвинеи, австралийских аборигенов (даже больше, чем у папуасов), полинезицев, фиджийцев, восточных индонезийцев с разных островов, филиппинцев и у филиппинских аэта-маманва.»

Что касается неандертальца, то уже с 2010 года известно, что в целом неандертальская ДНК составляет 1-4% генома нынешних людей, живущих за пределами Африки. Авторы двух исследований, опубликованных в среду журналах Science и Nature, выяснили, что чаще всего неандертальская наследственность присутствует в нескольких генах, связанных с выработкой кератина, присутствующего в коже, волосах и ногтях. В этой части генома неандертальские аллели обнаружены у 70% европейцев и 66% азиатов.

Гораздо интереснее те мои результаты, которые отличаются от общепринятых. Так например, довольно неожиданным результатом является наблюдаемое на графике значительное смещение южноамериканских индейцев в сторону денисовского человека, причем это смещение гораздо значительнее смещения папуасов и меланезийцев, у которых были найдены «денисовские гены» в наибольшем количестве. Что это означает, трудно сказать — наличие реального сигнала смешивания в данном случае равновероятен обнаружению статистического артефакта.  Впрочем, если верить работам Скоглунда этот результат может быть правдоподобным — моделирование миграций генов показало, что «денисовские» гены должны встречаться не только в Юго-Восточной Азии, но даже в некоторых группах Южной Америки (Skoglund et Jakobsson, 2011)

Оставим в стороне этот вопрос, который нуждается в более детальном изучении, и передем к расчетам процентной величины вклада неандертальских генов в популяции современных людей. Очевидно, что средняя величина этого вклада по каждой из популяций может дать только приблизительное представление о характере архаичной интрогресси неандертальских генов. Индивидуальный уровень вклада в каждой популяции может иметь большую частотную амплитуду в интервале между 1 и 6% процентами. Тем не менее, представляется возможным апроксимировать эти значения путем умножения собственного вектора (eigenvector) главных компонентов каждого индивида каждой популяции на собственное число линейного преобразования (eigenvalue), и последующим усреднением по популяции.

Ниже приведены эти усредненные значения в процентах (неандертальских генов), в порядке уменьшения. Вызывают сомнения ультра-высокие значения в первых десяти популяциях — скорее всего это результат комплексного воздействия статистических эфектов недостаточной представленности выборки, а также высокой степени гомозиготности, характерной для изолированных популяций (исландцев, албанцев и басков). Довольно высок уровень неандертальского вклада в образцах древних европейцев, хотя это и логично с точки зрения исторической модели адмикса. С другой стороны, средние значения (2-2.7%) неандертальского адмикса в популяциях Восточной Европы выглядят реалистичными. Так, например, по расчетам 23andme у меня уровень «неандертальских генов» составляет 2.67% :

Icelandic 10.50%
Norwegian 9.00%
1_Motala12 8.00%
Spain_BASC 8.00%
Albanian 7.00%
Korean 7.00%
Tiwari 5.11%
1_LBK380 5.00%
1_Loschbour 5.00%
French_South 4.00%
Kashmiri 4.00%
Tubalar 4.00%
Atayal_Coriell 3.60%
Ami_Coriell 3.10%
1_Motala_merge 3.00%
Bolivian 3.00%
Croatian 3.00%
Totonac 2.80%
Qatari 2.71%
Mixed_East_Slav 2.57%
Gujarati 2.43%
Ulchi 2.39%
North-Russian 2.36%
Center-Russian 2.36%
Aonaga 2.33%
British 2.33%
Chenchu 2.33%
East-Belarusian 2.33%
Ukrainian 2.33%
Finn 2.29%
Latvian 2.29%
Mixed_European 2.28%
South-Russian 2.27%
Pole 2.26%
Lithuanian 2.25%
West-Belarusian 2.25%
Belarusian 2.23%
Vepsa 2.23%
Bosnian 2.22%
Cree 2.20%
Georgian_Imereti 2.20%
Polish 2.20%
Orcadian 2.15%
Russian 2.15%
Karelian 2.13%
Welsh 2.12%
Swede 2.11%
Ukranians 2.11%
Greek 2.10%
Lithuanians 2.10%
Gagauz 2.09%
Croat 2.08%
Slovak 2.08%
Estonians 2.08%
Adygei 2.07%
Serb_Serbia 2.07%
Toscani 2.07%
French 2.06%
Komi 2.06%
1_LaBrana 2.00%
Algonquin 2.00%
Avar 2.00%
Azeri_Dagestan 2.00%
Azov_Greek 2.00%
Bashkir 2.00%
Belgian 2.00%
Bulgarians 2.00%
Central-Greek 2.00%
CEU 2.00%
Cirkassian 2.00%
Cochin_Jew 2.00%
Corsican 2.00%
Cretan 2.00%
Croat_BH 2.00%
Don_cossack 2.00%
Eskimo 2.00%
Haida 2.00%
Hungarian 2.00%
Hungarians 2.00%
Inkeri 2.00%
Inkeri-Finn 2.00%
Italian_Abruzzo 2.00%
Kets 2.00%
Kosovar 2.00%
Kryashen 2.00%
Kuban_cossack 2.00%
Lezgin 2.00%
Macedonian 2.00%
Meghawal 2.00%
Mishar 2.00%
Mixed_CEU 2.00%
Mixed_East_European 2.00%
Mixed_German 2.00%
Mixed_Slav 2.00%
Montenegrian 2.00%
Mordovian 2.00%
Mordovians 2.00%
North_Italian 2.00%
Occitan 2.00%
Roma_Bulgarian 2.00%
Roma_Macedonian 2.00%
Romanian_Jew_2 2.00%
Russian_South 2.00%
Saami 2.00%
Selkup 2.00%
Serb_BH 2.00%
Slovenian 2.00%
South_Greek 2.00%
Swedish 2.00%
Tabassaran 2.00%
Tatar_Lithuanian 2.00%
Velama 2.00%
West_Greenland 2.00%
French_Basque 1.95%
Chechens 1.94%
Iberian 1.94%
Chuvash 1.94%
Tatar 1.93%
Balkars 1.92%
German 1.92%
North-Ossetian 1.92%
Hant 1.89%
North_Greek 1.89%
Georgians 1.88%
Lak 1.88%
Abhkasians 1.85%
Sardinian 1.84%
Udmurd 1.84%
Maris 1.82%
Romanians 1.82%
Georgian_Laz 1.80%
Kumyks 1.80%
Lodi 1.80%
Mansi 1.77%
Chukchis 1.75%
Crimean_Tatar 1.75%
Italian_Piedmont 1.75%
Ket 1.75%
Moldavian 1.75%
Vaish 1.75%
Hallaki 1.67%
Lezgins 1.67%
Ossetian 1.67%
Tlingit 1.67%
Greek-Islands 1.63%
Turks 1.63%
Armenians 1.60%
Nogais 1.60%
Selkups 1.60%
Hakas 1.57%
Ashkenazy_Jews 1.56%
Apache 1.50%
Jew_Tat 1.50%
Kabardin 1.50%
Karitiana 1.50%
Kurds 1.50%
Nenets 1.50%
Samaritians 1.50%
Santhal 1.50%
Srivastava 1.50%
Syrian_Jew 1.50%
Tuva 1.50%
Uygur 1.50%
Mexican 1.45%
Italian_Jew 1.40%
Portugese 1.40%
Tajiks 1.40%
Kyrgyzians 1.38%
Roma_Slovenian 1.38%
Altaians 1.36%
Koryaks 1.33%
Pashtun 1.33%
Satnami 1.33%
Sicilian 1.33%
Yakut 1.31%
Cypriots 1.30%
Spaniards 1.30%
Turkmen 1.30%
French_Jew 1.29%
Iraqi_Jews 1.29%
Sephardic_Jews 1.29%
Turkmens 1.29%
Parsi 1.28%
Buryats 1.27%
Pathan 1.27%
Tadjik 1.27%
Athabask 1.25%
Iran_Jew 1.25%
Kurd_Jew 1.25%
Nganassans 1.25%
Nysha 1.25%
Azeri 1.22%
Mixtec 1.22%
Tharu 1.20%
Tunisian_Jew 1.20%
Uzbek 1.20%
Evenkis 1.18%
Kazakhs 1.18%
Roma 1.17%
Tuvinians 1.17%
Druze 1.16%
Karakalpak 1.14%
Mongolians 1.14%
Uzbeks 1.13%
Ojibwa 1.10%
Buryat 1.00%
Cochimi 1.00%
Cucupa 1.00%
Dolgan 1.00%
Dolgans 1.00%
Even 1.00%
Evenk 1.00%
Hazara 1.00%
Huichol 1.00%
Kalash 1.00%
Kalmyk 1.00%
Kamsali 1.00%
Koryak 1.00%
Kumiai 1.00%
Lambadi 1.00%
Luiseno 1.00%
Maya 1.00%
Mongol_Halha 1.00%
Nganassan 1.00%
Oroqen 1.00%
Pima 1.00%
Roma_BH 1.00%
Romanian_Jew_1 1.00%
Romanian_Jew_3 1.00%
Shor 1.00%
Surui 1.00%
Tharus 1.00%
Tsimsian 1.00%
Uyghur 1.00%
Uzbekistan_Jew 1.00%
Uzbekistani_Jews 1.00%
Vysya 1.00%
Yukaghirs 1.00%
Sindhi 0.91%
Hezhen 0.86%
Xibo 0.80%
Navajo 0.78%
Bhil 0.75%
Brahmins_UP 0.75%
Burusho 0.75%
Mongola 0.75%
Naga 0.75%
Iranians 0.71%
Daur 0.67%
Kshatriya 0.67%
Mala 0.67%
Moroccan_Jews 0.67%
Japanese 0.58%
Chinese_Dai 0.53%
Evens 0.50%
Kol 0.50%
Morocco_Jew 0.50%
Mumbai_Jews 0.50%
Scheduled_Caste_UP 0.50%
South_Han 0.50%
Tu 0.50%
North_Han 0.45%
Brahui 0.45%
She 0.44%
Tujia 0.44%
Iraki 0.43%
Naxi 0.43%
Dharkars 0.40%
Han 0.40%
Kanjars 0.40%
Miaozu 0.40%
Velamas 0.38%
Balochi 0.33%
Chenchus 0.33%
Dusadh 0.33%
Hakkipikki 0.33%
Lahu 0.33%
Piramalai_Kallars 0.33%
Yizu 0.33%
Colombian 0.25%
Chamar 0.22%
Syrians 0.22%
Dai 0.20%
Libyan_Jew 0.17%
Makrani 0.08%

Деанонимизация через генетическую информацию

Обычно на Хабрахабре не пишут интересных сообщений по геномике, но вот этот вот пост заслуживает определенного внимания.

В cети в открытом доступе некоторое время назад появились базы данных с генетической информацией людей (информацией разного уровня детализации — от полных последовательностей (сиквенса) всего генома до ограниченной информации по коротким тандемным повторам Y-хромосомы (Y-STRs). Например, энтузиасты делятся информацией о своих Y-STRs (гаплотип) на генеалогических сайтах для выяснения родственных связей и поиска дальних родственников, эти данные не анонимны. Так же в свободном доступе находится анонимная медицинская генетическая информация, например из научного проекта «1000 геномов человека» (проект по полной расшифровке геномов тысячи разных людей), где анонимность доноров ДНК поддерживается по этическим причинам.

Здесь начинается самое интересное. Генеалогические базы данных (даже весьма плохо заполненные, но тем не менее) позволяют деанонимизировать людей. Например, показано, что в случае искусственного оплодотворения спермой от анонимного донора, использование генеалогических баз данных позволяет узнать по крайней мере фамилию настоящего биологического отца ребенка (то есть через очень дальних родственников, засветившихся в базе, узнать, из какой семьи был донор), а при наличии дополнительной информации, такой, как место проживания и т.п., позволяет однозначно идентифицировать биологического отца. Недавно было показано, что находящиеся в свободном доступе анонимные генетические данные, плюс дополнительная информация о возрасте и т.п., позволяют точно установить личности примерно 50 анонимных доноров ДНК из проекта «1000 геномов человека». Это весьма настораживающее достижение, так как полная генетическая информация этих людей, находящаяся в открытом доступе, содержит данные об их предрасположенности к тем или иным заболеваниям и т.д., она может быть корыстно использована страховыми компаниями и подобными организациями.

Генетические компоненты в геномном калькуляторе MDLP-World22: взгляд со стороны. Анализ Евгения Марчева.

В качестве краткого вступления.

В первую очередь позвольте мне представиться терпеливому и мудрому читателю моего блога — Евгению Марчеву и выразить ему свое почтение и большое удовольствие, которое мне доставляет обращение к нему. Я вполне знаком с его характером, так как, судя по предмету и стилю моего блога, можно поручиться, что он один из миллионов. Подобный критически настроенный читатель готов столкнуться с утверждениями, от которых он будет склонен первоначально отступить; но в конце концов он ждет, чтобы его убедили в истинности каких-либо из них. Он также догадывается, что размышления над предметом и попытка изложить размышления в письменном виде заняло не скажу сколько времени, но уж определенно больше, чем четверть часа, и, следовательно, фундаментальные возражения, столь очевидные, что могли бы ненароком поразить каждого, пришли автору на ум неожиданно, хотя ответы на них не из тех, полное значение которых можно оценить сразу.

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

Итак, представляю читателю самостоятельный анализ Е.Марчева, предметом разбора которого являются рассчитанные  в моем проекте аллельные частоты отдельных предковых компонентов.

I. Теория

Решил немного проанализировать данные из этой таблицы применительно к европейцам с уклоном на север, центр и восток Европы.

У европейцев 7 наиболее частых компонентов:

West Asian, North-European-Mesolithic, Atlantic_Mediterranean_Neolithic, Samoedic, Indo-Iranian, North-East-European и Near_East. 

Разберем вкратце каждый из них.

Итак, первый West Asian. Наибольший процент на Западном Кавказе- абхазы (64, 9%), имеретинцы (63,7), лазы (56,6), аварцы (56,8), лезгины (55,4).
На интересуемой нас территории этот компонент чаще встречается у казанских татар (9,9 %), южных немцев (8,4), украинцев (от 6,6 до 7,7 %), южных русских (6,2%). На западе высок процент у итальянцев (21,5 % у центральных итальянцев), французов (6,7 %).

Второй компонент North-European-Mesolithic. Условно его можно назвать саамско-финским.
Саамы (76,4 %), Финны (от 30,1 до 37,3 %), вепсы (24,1), карелы (23,2), ижорцы (22, 7).
Част этот компонент у северных русских (10,5 %), норвежцев (9,8 %), шведов (7,8 %), эстонцев (7,1 %).

Третий- Atlantic_Mediterranean_Neolithic. Юго-западно-европейский или просто западно-европейский.
Больше всего у сардинцев (68,1 %), басков (59,2 %), иберийцев (48,8), корсиканцев (47,8), португальцев (46,6), северных итальянцев (44,3), французов (43,5 %).
Вообще, данный компонент достаточно выражен у всех западноевропейцев- более 30 %.
У восточноевропейцев ниже- чехи (27,8), поляки (18,4), украинцы ( от 17 до 21%), русские (от 11 у северных до 17,3 у южных), у коми (8,9 %), манси (8,8 %).

Следующий Samoedic. Уральский.
Значительно присутствует у селькупов (68,1%), хантов (64,6), ненцы (37,1), манси (30,9 %- ожидал больше), удмурты (29,6), марийцы (27, 8), шорцы (22,0 %), башкиры (21,7%), чуваши и хакассы по 17,6 %, коми- 16,4 %, казанских татар (11,9 %).
У западноевропейцев этот компонент практически не встречается. Есть у русских (от 1,0% у центральных до 4,7 % у северных), у карел (1,6%), словаков (1,4%), западных украинцев (1,7 %).

Пятый. Indo-Iranian.
Очень интересный. Чрезвычайно большой процент больше у калашей- 84,4 %. Вроде и у соседей должно быть прилично, ан нет: буриши (8,7 %), пуштуны (8,5 %), белуджи (6,9 %). Почему так- не знаю.
Что у европейцев? Встречается у норвежцев (3,8%), валлийцев (3,4), бритишей (2,8), французов (2,1). В Восточной Европе как правило до 1%.

Следующий North-East-European занимает значительное место у европейцев, почти у всех.
Больше всего у балтов и славян. Литовцы (81,9), латыши (79,5), белорусы (76,4), эстонцы (75,2), поляки (70,2), русские (67- 70,4), украинцы (62,1- 67,1), сорбы (65,9), карелы (60,2), вепсы (62,5).
Но много и у других: чехи (57,4), северные немцы (54,6), южные- 42,6, у британцев от 46 до 49, норвежцы- 48,1, шведы- 53,7! Низок у саамов- 15,5%.
У южан сравнительно меньше- французы (38,4), центр. итальянцы (18,9%).

И последний- седьмой. Near_East. Условно «семитский» (по имени Сима) или арабский, кому как угодно.
Евреи Йемена (60,9 %), Сауд. аравия (59,5), бедуины (56,7), евреи Эфиопии (52,5), египтяне (43,8).
В Европе:
Относительно много у итальянцев (центр- 17,4), португальцев (11,9).
У южных немцев (3,5), украинцы (от 2,3 у восточных до 3,8 % у западных), чехи (3,0), белорусы (3,4), словаки (3,2), у русских от 1,0 до 1,5%, у литовцев- 1,4%, у поляков- 1,3 %. У других почти нет.

И две диаграммы для наглядности:

II. Практический пример

Рассмотрим частный случай на примере ув. exorcio. Происхождение- северный русский/немец.

Основные компоненты (8):

West-Asian 4.12%- промежуточное значение
North-European-Mesolithic 6.87%- промежуточное значение
Arctic-Amerind 0.99%- градиент к сев. русским
Atlantic_Mediterranean_Neolithic 22.20%- градиент к немцам
Samoedic 2.77%- градиент к русским
East-Siberean 1.49%- градиент к русским
North-East-European 57.77%- промежут.
Near_East 1.15%- промежут.

Хорошо видно, что градиент в сторону северо-восточных популяций (русские, финны) за счет трех компонентов- Arctic-Amerind, Samoedic, East-Siberean, а в сторону немцев за счет Atlantic_Mediterranean_Neolithic.*

В конечном результате в пространстве PCA генетического разнообразия, excorcio находится рядом с чехами, то есть между немцами и восточными славянами.

———————

* Комментарий уважаемого А.Бурнашева: Я хочу заметить, что с математической точки зрения калькуляторы Веренича и Понтикоса не страдают специфическим дефектом калькулятора Веселовского (Поляко) — у них все компоненты независимы.В то же время, в EU-test Поляко появляется один зависимый компонент (т.е. принадлежность к нему куска генома вычисляется на основе других компонент), а в J-test таких зависимых уже 2 (очевидно, ashkenazi компонент тоже вычисляется как суперпозиция).При проведении PCA анализа это все видно.