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

Я поработал тут над протоколом определения мужских Y-гаплогрупп в палеоДНК. В конце концов — через пару дней — я остановился над следующим варианте.
Протокол содержит две части — первая для геномов с высоким покрытием, вторая для геномов с низким качеством и малым покрытием.


1) Для геномов с высоким покрытием варианты Y определяются в программе GATK и выводятся в формат VCF
Файл VCF вводится в программу yHaplo (написанную Позником на основании алгоритма определения Y-гаплогруппы в 23andme)
2) Для геномов с низким покрытием используется программа samtools mpileup c параметрами -B -q30 -Q30 -C50. Файл пайлапа преобразуется в формат 23andme и вводится в ту же программу (yHaplo)


Я проверил работоспособность протокола на нескольких примерах, похоже все работает (варианты гаплогрупп в таблице совпадают с теми что были опубликованы в статьях)
Сначала геномы с высоким покрытием — 2 генома древних гладиаторов из Йорка

3DT26 J-CTS8938 J-M304 J
6DT3 R-L52 R-P311 R1b1a2a1a

Теперь геномы с низким покрытием — 2 древнеегиптских образца

ERR1654486 J-P58 J-P58 J1a2b
ERR1654487 E-V22 E-L677 E1b1b1a1b2

Теперь еще более экстремальный случай (качество и покрытие плохое) — геномы римского периода с территории Польши (предположительно готы из Вельбарской культуры)

kow45 I-L35 I-M436 I2a2
kow55 I-L80 I-M253 I1

Геном англосакса из Йорка

NO3423 I-DF29 I-DF29 I1a

Геном неолитического периода с территории Польши

pl-7 R-S24902 R-S24902 R1a1a1b1a2c

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

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

Данная заметка представляет собой критический анализ методологических предпосылок создания неандертальского калькулятора, имплементированного в 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%

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

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

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

В качестве примера в нашем туториале мы будем использовать данные, любезно предоставленные авторами работы 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, о котором мы поговорим в следущей части туториала.

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

Еще раз о экзомном тестировании — пример отчета Exome 80x 23andme

Пару дней тому назад в своем блоге я изложил

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

  554597_4120657144542_494648084_n 303382_4120657544552_2124700585_n 252231_4120657784558_1755564213_n 479720_4120658184568_761081981_n 285695_4120658504576_1081385431_n 527651_4120658824584_192594454_n 551743_4120659024589_475852363_n 599435_4120659384598_2023625463_n 527612_4120659544602_1322921037_n 250942_4120659744607_1560066672_n 542556_4120660024614_715032428_n 250944_4120660464625_514191875_n 224935_4120660584628_446354628_n

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

В апреле прошлого года начал работать над добавлением в своей BGA проект SNP-ов извлеченных из «ископаемых» геномов трех неолитических охотников-собирателей и одного неолитического земледельца из недавней, но уже нашумевшей статьи Понтуса Скоглунда и др. (2012):

«Новые исследования ДНК древних останков помогли исследователям доказать, что сельское хозяйство пришло в Европу вместе с древними переселенцами. Это ставит под сомнение альтернативную теорию, которая предполагает, что сельское хозяйство было развито уже существующими в Европе племенами, которые до этого занималось охотой и собирательством, а затем стали делиться новым опытом друг с другом. Команда исследователей сравнила ДНК скелета древнего фермера со скелетами трех охотников. Они обнаружили, что «фермер» отличался на генетическом уровне от охотников. Понтус Скоглунд (Pontus Skoglund) из Университета Упсалы, Швеция, и его коллеги выделили генетический материал, который содержался в ядре клеток, из скелета человека, который жил 5 тысяч лет назад в южной части современной Швеции. ДНК, полученное из этих останков было на самом деле очень древним и не имело примесей современной ДНК, что позволило исследователям использовать молекулярные и статические техники. Они сравнили генетический профиль фермера каменного века и охотника-собирателя с современным населением. Женщина-фермер, как оказалось, родилась в этих местах, так как ее геном очень похож на геном современного человека из юго-восточной Европы.»

В настоящее время Диенек Понтикос поделился со мной ссылкой на BAM-элайнменты геномов древних шведских фермера и древних шведских охотников (Ajv52, Ajv70, Ire8, Gok4 и Ste7 ). Планируется извлечь известные (dBSNP-cовместимые) снипы с помощью алгоритма UnifiedGenotyper, имплементированного в коде G.A.T.K.

Полученные на выходе VCF файлы были конвертированы в PED/MAP формат с помощью vcftools.Как это уже стало привычным, в случае определения генетических вариантов в древних геномах нужно быть предельно осторожным и применять дополнительные QC-фильтры.  Понтус Скоглунд дал следующие рекомендации по работе с древней ДНК эпохи раннего шведского неолита.

По этой ссылке  можно скачать BAM-файлы с данными  генома неолитических жителей Швеции. Согласно установленному протоколу, эти данные приведены в соответствие с геномом современного человека. Что касается извлечения SNP-ов, то тут нет универсального простого рецепта.  Самым простым способом, вероятно, является циклическое (‘loop’) извлечение SNP-ов из исходных методом одиночного извлечения (случайных) аллель  из ‘древней’ ДНК с использованием следующей команды:

samtools mpileup <bam_file>-Q 30-Q 15-е hg18.fasta-р <chr>: <start> — <stop>

, где hg18.fasta — это референс человеческого генома версии UCSC в формтае FASTA.

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

Третье важное замечание. В связи с низким охватом (low coverage) генома по многим геномным позициям были получены только гаплоидные данные. Поэтому перед использованием данных с целью сравнения генетического разнообразия крайне важно убедиться, что данные древних геномов сравниваются с соответствующими данными современных популяций предварительно приведенным в гаплоидную фазу. В противном случае, могут возникнуть странные эффекты в большинстве анализов структуры генофонда населения. Кроме того, в полученных данных был отмечен  избыток C>T и G>A. Ошибки в определении  древних данных из-за посмертных (postmortem) повреждений, так что при сохранении этих полиморфизмов было бы неплохо, во избежании неожиданности, исключить их. Наконец, нужно быть осторожным при слияния SNP-ов древнего ДНК лица (приведенных в соответствие с  hg18) с данных где генотипами выравненными по референсному геному hg19, поскольку если оба аллеля присутствуют в обоих наборах данных, то flip — формальная команда Plink для замены снипов одной цепи (strand) ДНК на снипы другой цепи  в этом случае работать не будет.