Этногеномика беларусов — часть V

Обсуждение результатов и выводы

 

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

После внимательного изучения результатов нашего исследования,можно сказать, что оба из приведенных выше заключений представляют собой крайне упрощеные варианты сложного процесса формирования аутосомного генофонда беларусов. Хотя мы и не можем предоставить окончательных аргументов в пользу или опровержение каждой из этих версий, мы может предоставить более полное и подробное обозрение структуры аутосомного генофонда. В отличие от трех основых компонентов, упомянутых выше, в нашем исследовании мы выделили шесть основных компонентов, типичных для европейцев в целом. Основу генофонда составляет компонент, который мы обозначили как северо-восточно-европейский компонент. Именно этот компонент выделяет беларусов среди других восточных славян, приближая их к современным балтийским популяциям (у литовцев процент компонента составляет 81,9, у латышей — 79,5%, у беларусов -76,4%, у эстонцев — 75,2%). Примечательно, по мере удаления от территории Беларуси на север в с торону Латвии и Эстонии, увеличивается процент северо-европейского генетического компонента (как мы полагаем, этот компонент доминировал в генофонде доисторических жителей Скандинавии в эпоху до распространения финно-угоров и индо-европейцев). С другой стороны, беларусов и других восточных славян отдаляет от балтов и сближает друг к другу более высокий процент так называемого западно-азиатского или кавказского компонента (любопытно, что в этом случае эта закономерность может свидетельствовать в пользу западно-азиатской теории происхождения индо-европейцев).

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

Исходяизвышенаписанного,представляетсялогичнымсделатьвыводотом,чтоосновнойкритическийэтапстановленияаутосомногогенофондапришелсянапериодсмешиванияносителейсеверо-восточно-европейскогогенетическогокомпонентасносителямизападно-азиатского(кавказского)генетическогокомпонента,послечегопредковыйаутосомныйгенофондбеларусовприобрелотносительнуюстабильность.Разумеется,даннаямодельнеисключаетпозднейшиеэпизодысмешиванияпопуляций,ноониоставилименьшийследвструктуреаутосомногогенофондабеларусов.Вэтойсвязивозникаеточевидныйвопрос–вкакойименноисторическийпериодпроизошлосмешениеносителейсеверо-восточно-европейскогогенетическогокомпонентасносителямизападно-азиатского(кавказского)генетическогокомпонента,иктобылиихносителями?
В начале сентября 2012 года известная американская лаборатория популяционной генетики доктора Райха опубликовала альфа-версию программного продуктаADMIXTOOOLS1.0. Альфа-версия была разработана для внутреннего использования, поэтому modusoperandiэтого продукта вряд ли является кристально понятным для стороннего пользователя. Положительным аспектом на мой взгляд является то, что ADMIXTOOLSпакет обеспечивает полную совместимость с форматом другой очень популярной программыEIGENSOFT, которая была разработана в той же лаборатории. Это немаловажное обстоятельство намного упрощает процесс обучения в ADMIXTOOLS.

Вышеупомянутый пакет включает в себя 6 приложений, среди которых я считаю наиболее полезнойqp3Popи утилиты для вычисления частотной характеристики аллелей. Впрочем, я не собираюсь обсуждатьqp3popво всех деталях и в контексте данной заметки достаточно отметить, что эта программа реализует тест three_pop(F_3), подробно описанный в известной статье Рейха и соавт. 2009.

Однако другой имплементированный в пакете метод, – метод rolloff– нуждается в более пристальном внимании. Этот метод позволяет производить математическую оценку как времени, так как и уровня адмикса. Оценка производится на основании анализа неравновесия по сцеплению между SNP-ами. Тут необходимо вспомнить стандартное определение неравновесия по сцеплению.Неравновесием по сцеплению (часто используется английская аббревиатураLD) называется неслучайная связь между двумя аллелями, в силу которой определенные комбинации аллелей встречаются наиболее часть. В теории, чем дальше друг от друга находятся SNP-ы ,тем меньше будет уровень LD. Темп угасания снижения LDв адмиксе напрямую связана с числом поколений, прошедших с момента адмикса, так как cвозрастанием числа поколений увлечивается число рекомбинаций произошедших между двумя отдельными SNP-ами. Проще говоря: Rolloffсоответствует экспоненциальной кривой угасания уровня LDот расстояния, и эта скорость экспоненциального снижения как раз и используется для оценки числа поколений, так и уровня адмикса в анализируемой популяии. Учитывая, что одно поколение примерно равно 29 лет, можно преобразовать число поколений в года.

Этот метод открывает интересные перспективы. Для целей этого анализа, я создал специальный набор SNP-данных, который включает в себя около 750 000 cнипов, частично или полностью в 250 различных популяциях человека. Далее, я разбил популяции 3 * 62 000 трио в следующем виде (X, Y, Z), где X и Y – пара рефренсных групп, а Z – белорусы из коллекцииBehar et al.2010. После этого я провел q3Pop анализ этих трио.

Результаты изложены в нижеприведенной таблице

Indian Polish Belarusian -0.000736 0.000251 -2.935
Polish Indian Belarusian -0.000736 0.000251 -2.935
Karitiana Sardinian Belarusian -0.001278 0.000517 -2.471
Sardinian Karitiana Belarusian -0.001278 0.000517 -2.471
Otzi North_Amerind Belarusian -0.002556 0.001126 -2.271
Cirkassian Polish Belarusian -0.000488 0.000231 -2.113
Polish Cirkassian Belarusian -0.000488 0.000231 -2.113
Pima Otzi Belarusian -0.002727 0.00137 -1.99
Pima Sardinian Belarusian -0.000794 0.000431 -1.843
Sardinian Pima Belarusian -0.000794 0.000431 -1.843
Otzi Surui Belarusian -0.002938 0.001931 -1.522
Surui Otzi Belarusian -0.002938 0.001931 -1.522

 

На первый взгляд, результаты нашего эксперимента с 3qPop, кажется, неплохо согласуются с выводами, содержащимися в работеПаттерсон и др. 2012: “Самый поразительный вывод состоит в обнаружени четкого сигнала адмикса в северной Европе, один из элементов которого связан с предками населения наиболее близкого по своей генетике к баскам и жителям Сардинии, а другой – с предками современного населения северо-восточной Азии и Америки. Этот явный сигнал, вероятно, отражает историю смешивания неолитических мигрантов с коренным населением Европы, что подтверждается недавним генетическим анализом древних костей Швеция и секвенированием полного генома Отци Тирольца”. Что касается собственно белорусов, то источники сигнала смешивания с посторонними популяцими менее ясны и расплывчаты. Как было показано ранее, с точки зрения формального анализа примесей (f3 статистики), белорусы могут быть представлены в виде популяционного микса поляков и индусов / черкессов. Первый компонент смеси может быть связан с носителями культуры шнуровой керамики/боевых топоров и культуры колоковидных кубков; второй, в соответствии с результатами, должен быть общим для индусов и черкесов.

 

Белорусы = ((неолитические культуры Европы) + “носители культуры колоковидных кубков”) + (мезолитическое население Европы) + компонент носителей культуры шнуровой керамики)) + скифо-сарматский тип

 

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

 

rolloff

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

154,158 + -87,024 поколений назад (или, 4470 + -2523 года до настоящего времени / 2510 – +2523 лет до н.э.).

 

Исходя из этого, мы решили модифицировать Rolloff-анализ генофонда белорусов, используя на этот раз в качестве референсов литовцев и пуштунов. Следуя этому совету, я решил предпринять вторую попытку формального анализа адмикса в двух имеющихся у нас выборках беларусов ( выборка беларусов из статьи Behar et al. 2011), и выборка беларусов, собранная в нашем проекте.Ниже приведены результаты эксперимента с двумя этими группам (в отличие результатов нашей предыдущей попытки, результаты данного эксперимента менее “зашумленные”):

rolloff2

 

Интервал числа поколений, прошедших со времен анализируемого адмикса (105.086+-52.59) или 3069 +- 1525 лет до настоящего времени, что соответствует временном интервалу 2 тыс. до нашей эры – 6 век нашей эры. Принимая во внимание эти выводы, мы можем предположить, что основной аутосомный эпизод смешивания предковых популяций беларусов произошел в течении довольно таки продолжительного времени, охватывающего несколько тысяч лет. В этой связи, вопрос о том, кто именно был носителями северо-восточно-европейскогогенетическогокомпонентасносителямизападно-азиатского(кавказского)генетическогокомпонента, остается открытым.

Реклама

Этногеномика беларусов — часть III

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

 

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

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

Pygmy

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

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

simul

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

fst dist

Рисунок 1. FST-дистанции между компонентами

 

Приведенная выше матрица Fstдистанций была использована для определения наиболее вероятной топологии NJ-дерева всех 22 предковых компонентов ( примечание: в качестве outgroup-таксона использовался South-Africancomponent).

Географическое распространение компонентов нового калькулятора проекта MDL K27

Не успела бета-версия моего нового этно-популяционного калькулятора и сопутствующего ему геномного оракула (Dodecad oracle) пойти в массы, как один талантливый россиянин Сергей Козлов из Новосибирска (о котором я уже много раз упоминал в блоге) написал программу, позволяющую довольно точно проецировать/предсказывать ареал происхождения анализируемого человека по мере степени увеличения или убывания  процентов предковых компонентов (или аллельных частот)  в отношении к априори заданым точкам на контурным картам (эти точки на контурной карте соответствуют контрольным группам референсных популяций).

Отклоняясь в сторону от темы, хочу отметить что два года назад, когда я начал работу над проектом MDL, я не надеялся найти активных последователей среди русскоговорящего населения, хотя задекларированный в анонсе ареал проекта частично охватывал часть современной западной европейской части России.  Причина моего пессимизма была очевидна — современные русские (впрочем как и 90% прочего человеческого населения) ленивы, глупы и любят бесплатно паразитировать на результатах труда других людей.  К началу 2011 года можно было по пальцам пересчитать тех русскоязычных людей, которые занимались  практическим изучением аутосомного родства и изучения происхождения, или создавали соответствующее программное обеспечение. По прошествии 2 лет,  я должно признать, что в своих мрачных прогнозах немного ошибался.   К счастью, не перевелись еще в известных российских IT-селениях вроде Новосибирска энтузиасты-кулибины.  А это означает, что надежда на пробуждение массового  интереса к ДНК-генеалогии в РФ, так же как это произошло уже в США, где уже сейчас можно наблюдать геномную революцию  во всей ее динамике.

Но вернемся к теме.

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

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

В связи  с этим, нужно отметить два важных нюанса.
Во-первых,  число реперных точек (т.е «реперных» популяций) по европейской части РФ у Сергея  гораздо выше, чем в оригинальной бета-версии моего калькулятора K27 . Число точек в модификации Сергея было увеличено за счет включения фиксирующих дополнительных групп народонаселения РФ.
Во-вторых, cама идея визуализации  геномных данных на географической карте далеко не нова.  Весной этого года, в своей большой обзорной статье о принципах созданиях этно-популяционных калькуляторов на примере MDL World K22,  я указал на возможность визуализации коэффициентов адмикса в географическом пространстве:

… я решил визуализировать компоненты на поверхности земного шара путем отображения коэффициентов адмикса. Избегая излишних премудростей, я воспользовался готовым рецептом Франсуа Оливье, который предложал  использовать графическую библиотеку статистического программного обеспечения R для отображения пространственной интерполяции  коэффициентов адмикса (Q матрица) в двух измерениях (где пространственные координаты записываются как географические долгота и широта).

При вдумчивом прочтении подобных методов, встает неизбежный вопрос — почему градиенты аллельных частот в геноме людей являются крайне информативными при определении места их происхождения?   Частичный ответ на этот вопрос можно найти в другой моей заметке «О новых перспективах геномной геногеографии: SPA анализ участников проекта MDL«. В этой заметке я обсуждал перспективу расширения традиционных геногеографических методов, так как эксплицитное пространственное моделирования частот аллелей позволяет достаточно точно локализовать положение отдельно взятого человеческого индивида  на географической карте только на основании генетической информации. Если географическое происхождение лиц известно априори, то можно использовать эту информацию для определения функции частот аллелей в каждом SNP. Однако, если таковая информация отсутствует, то наша модель позволяет определить географическое происхождение физических лиц, используя только их генетические данные аналогично более известному методу  многомерного скалирования , основанному на определении пространственных координат статистических параметров.
Это заключение подтверждается в независимом исследовании компании 23andme, согласно которому анализ главных компонентов генетического разнообразия в геноме человека позволяет точно определить его место происхождения в Европе.

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

Легенда

Результаты «типичного русского»

Результаты карпатского русина

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

Ancestral-South-Indian Ancestral-Yayoi Arabic Australo-Melanesian Austronesian Baltic-Finnic Bantu Bushmen Caucasian-Near-Eastern Central-African-Hunter-Gatherers Central-African-Pygmean Congo-Pygmean Cushitic East-Siberean Gedrosia-Caucasian Kalash Nilo-Saharian Nilotic-Omotic North-African North-Amerindian North-Circumpolar North-European-Baltic Papuan-Australian South-Meso-Amerindian South-West-European Tibeto-Burman Uralic

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

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

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

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

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