Способ определения этно-географической группы происхождения и территории происхождения индивидуума и панель однонуклеотидных полиморфизмов
METHOD FOR DETERMINING ETHNO-GEOGRAPHICAL GROUP OF ORIGIN AND TERRITORY OF ORIGIN OF AN INDIVIDUAL AND PANEL OF SINGLE-NUCLEOTIDE POLYMORPHISMS

Область техники, к которой относится изобретение
Настоящее изобретение относится к области прогностических тестов определения происхождения человека в судебно-медицинской экспертизе и используется для определения этно-географической группы и территории происхождения предков индивидуума.
Уровень техники
Возможность определения вероятного этно-географического (биогеографического) происхождения индивида по его ДНК перспективна для решения нескольких задач, важных как для отдельного человека, так и для общества. Возможность узнать место происхождения индивида и принадлежность к той или иной популяции по анализу ДНК очень ценна для криминалистов, так как позволяет значительно сузить поиск потенциального преступника или жертвы преступления. Этно-географическое происхождение индивидов важно знать для персонализированного подхода в медицине, так как в разных популяциях мира значительно различается частота генетических вариантов, связанных с теми или иными болезнями и переносимостью лекарств. Наконец, все больше людей хотят знать, где жили их предки, или к какой этнической группе они принадлежали, что сформировало область генетической генеалогии, в том числе.
Установлено, что генетическое разнообразие народонаселения мира находится в зависимости от географического расположения популяций. Иными словами, генетический ландшафт планеты проецируется на географическую карту. Эта идея постулировалась еще классиками геногеографии и лежит в основе возникновения геногеографии как науки. Она воплощается в корреляции между генетическими расстояниями между популяциями (которые характеризуют их генетическое сходство) и географическими расстояниями между ними. Эта закономерность проявляется на всех уровнях исследования, с использованием разных генетических систем. Этническое происхождение является более сложным предметом для исследований, поскольку этническая идентичность определяется самосознанием человека, и поэтому может не соответствовать его биологическому происхождению. Сопряжение этнического и биологического происходит только через популяции – стабильные группы населения, заключающие большинство браков в своих пределах; большинство этносов являются популяциями, поскольку справедливо правило этнической эндогамии (более половины браков заключается с представителем того же этноса). Поэтому на практике определяется популяция, из которой происходят предки обследованного индивида, далее определяется, к какому этносу относится данная популяция, и таким образом определяется этническое происхождение изучаемого индивида. Эти исследования проводятся в рамках популяционной генетики человека. Популяционные генетики изучают генофонды популяций, и основной объект их исследований  – популяция, сопряженная с этносом [Юсупов и др., 2017].
Географическое происхождение индивида отражается в его геноме, а поскольку его геном складывается из фрагментов, унаследованных от родителей, то характерные участки укажут на происхождение родителей и других предков, вглубь поколений. В геноме существуют изменчивые участки, которые служат генетическими маркерами происхождения, так как передаются из поколения в поколение. Среди них есть однородительские маркеры – которые передаются от отца к сыну (это маркеры мужской Y-хромосомы) или от матери к ее детям (это маркеры на митохондриальной ДНК – мтДНК). Отслеживание этих маркеров реконструирует происхождение человека прицельно только по отцовским или только по материнским линиям наследования. Но есть и маркеры, которые разбросаны по всему геному – они называются аутосомными (генотипирование сотен тысяч таких маркеров, разбросанных по всему геному, называется широкогеномным анализом) и передаются человеку от обоих родителей; исследование этих маркеров дает усредненную информацию о всех предках человека. Для предсказания места происхождения человека на основе его генетической информации ключевым является поиск ДНК-маркеров, информативных для определения происхождения индивида с перспективами применения в криминалистике.
Эпоха широкогеномных и полногеномных исследований привнесла новые возможности поиска таких маркеров. Ниже упомянуты несколько работ, в которых впервые на широкогеномных данных были выявлены географические закономерности генетического ландшафта. Они были выявлены по данным анализа главных компонент (РСА), на основе которого и возник первый алгоритм, который по широкогеномным SNP маркерам может вычислять принадлежность индивида к определенной популяции или группе популяций в масштабах континента или внутри континента.
В одной из исходных для этого направления работ [Novembreetal., 2008] авторы описали генетическую вариабельность на выборке из 3000 европейцев, которых они генотипировали по более чем 500 тыс. маркеров ДНК на панели Affymetrix (Фигура 3). На основе генетических расстояний между популяциями они применили анализ главных компонент (РСА), на графике которого популяции расположились в двухмерном пространстве. Индивиды из одного и того же географического региона группировались в кластеры. Эти кластеры по форме напоминали очертания частей Европы. Такой географически зависимый характер генетического разнообразия соответствует идеальной теоретической модели, в которой генетическое расстояние между популяциями коррелирует с географическим расстоянием.
В работе сделан вывод, что метод анализа главных компонент позволяет на основании генетической информации поместить конкретного индивида в географическую точку его происхождения. Важной характеристикой метода является его точность. Как утверждают авторы статьи, используя метод множественной регрессии, в пределах Европы они смогли поместить 50% индивидов в пределах 540 км от места рождения и 90% индивидов – в пределах 840 км от места рождения. Они считают, что метод анализа главных компонент адекватен для предсказания географических координат места происхождения индивида и может быть использован для генетического тестирования.
SPA (Spatial Ancestry Analysis)
Вслед за этим стали появляться другие методы для предсказания биогеографического происхождения индивида по его ДНК. Одна группа методов относится к локальным, которые предназначены для определения происхождения отдельных сегментов хромосом, например, методы Lanc-CSV, LAMP-LD, MULTIMIX, ADMIXTURE, STRUCTURE, другие – глобальными, нацеленными на предсказания происхождения собственно индивида, например, HAPMIX, LAMP, SPA. Среди них выделяют метод , в основе которого лежит модель пространственного распределения для каждого SNP, использующая частоту аллелей как непрерывную функцию в географическом пространстве. Все упомянутые выше методы решают задачу определения происхождения индивида с приемлемой точностью на уровне континентов и субконтинентальных регионов, но на уровне отдельных стран их точность оставляет желать лучшего. Это послужило стимулом для работы над новыми методами.
GPS (Geographic Population Structure)
Шагом вперед в увеличении точности стала разработка метода , представленный в статье [Elhaiketal., 2016]. Этот метод основан на анализе смешения предковых компонентов популяций (на анализе ADMIXTURE). Метод GPS состоит из нескольких стадий. На первой стадии собирается панель популяций мира, которые анализируются на предковые компоненты методом ADMIXTURE (Фигура 4) - формируется база референсных популяций (референсные популяции – это популяции, живущие в течение долгого времени в данном географическом регионе). На втором этапе определяется пропорция генетических вкладов референсных популяций в геноме индивида. Наконец, происходит конвертирование генетического расстояния от индивида до референсных популяций в географическое расстояние и вычисляются географические координаты происхождения индивида.
Точность GPS положительно выделялась среди предшествующих методов. С помощью GPS авторам удалось предсказать происхождение 83% индивидов с точностью до страны и для 65% - с точностью до конкретного региона страны (Фигура5). Место происхождения для 50% индивидов из популяций всего мира удалось предсказать в пределах 87 км от реального места их рождения и для 90% индивидов - в пределах 645 км. Примененный к жителям Сардинии, метод GPS продемонстрировал еще большую точность: он поместил 25% жителей в свою родную деревню и 50% - в радиус 50 км от нее. В случае смешанных популяций, по словам Татьяны Татариновой, одного из авторов метода, GPS определяет происхождение индивида с точностью 320 км от места его рождения.
Как подчеркивают авторы статьи, основные ограничения подходов SPA и PCA для определения биогеографического происхождения индивидов связаны с тем, что они адекватно работают на относительно гомогенных популяциях, таких как европейские. В то же время подход GPS продемонстрировал успешность использования на смешанных популяциях. Благодаря своим преимуществам метод GPS имеет массу потенциальных применений. Например, в генеалогических исследованиях он может определить место рождения индивида, в криминалистических исследованиях он может дать информацию по ДНК о происхождении индивида (например, был использован в работе ряда компаний, предлагающих платное генетическое тестирование - ProsapiaGenetics, GPS Origins и др.). В то же время, как показали некоторые другие работы авторов, метод не всегда приводит к адекватным результатам, если с его помощью пытаются определить не индивидуальное происхождение индивида, а прародину того или иного народа, и проводят лингвистико-исторические интерпретации.
reAdmix
Другой метод, созданный для решения задач в этом направлении, можно считать логическим продолжением и развитием метода GPS, поскольку его разработали те же авторы, среди которых проф. Татьяна Татаринова, проф. Университета Южной Калифорнии, и Эран Элхаик (Университет Шеффилда, Великобритания) [Kozlov et al., 2015]. Это метод , предназначенный для того, чтобы в геноме индивида вычислить пропорцию предковых вкладов, соответствующих современным популяциям. Иными словами, его предназначение – определение этно-географического происхождения индивида, происходящего из смешанной популяции.
reAdmix
reAdmix
reAdmix –
reAdmix
mSpectrum, HAPMIX
LAMP
Суть  состоит в том, что том, что происхождение индивида представляется как взвешенная сумма современных референсных популяций. требует построения базы данных по популяциям мира, с применением анализа ADMIXTURE [Alexander, Lange, 2011] с разным числом компонент. После выбора оптимального числа предковых популяций (К) из частот аллелей для каждой предковой популяции формируют базу для последующих шагов. Индивиды проецируются на эту референсную базу К предковых популяций. Иными словами, индивидуальный генотип «разбивается» на определенную совокупность предковых компонентов. Например, в случае К=9 предковые популяции для индивида-это Северо-Восточная Азия, Средиземноморье, Южная Африка, Юго-Западная Азия, Америка, Океания, Юго-Восточная Азия, Северная Европа, Африка южнее Сахары. Конечная задача представить происхождение индивида как взвешенную сумму современных референсных популяций (например, 25% французов, 25% немцев, 50% японцев). Точность метода характеризуется так: в случае несмешанного происхождения для 94% индивидов была правильно определена этническая группа, и среднее расстояние до точного места происхождения составило 87 ± 21 км. В случае двух компонентов в геноме среднее расстояние до точного места происхождения составило 521 ± 74 км. В случае трех компонентов в геноме среднее расстояние до точного места происхождения составило 556 ± 55 км. Метод показывает большее соответствие реальному месту происхождения индивида, чем методы и (Фигура 6).
Существуют и средства для вычисления «этнических компонентов» в геноме, которые были созданы и создаются представителями «гражданской науки» - генетической генеалогии и геномными блогерами. Это так называемые «этнокалькуляторы» (это название условное, поскольку, как указывалось выше, они вычисляют не этнос, а популяцию, сопряженную с этносом). Базу ссылок на разные этнокалькуляторы можно найти на научно-популярном сайте «МолГен» http://forum.molgen.org/index.php/topic,5736.0.html. Обзор этнокалькуляторов дает на своем сайте один из известных российских генетических генеалогов и автор нескольких десятков калькуляторов Вадим Веренич: https://verenich.wordpress.com/.
Этнокалькуляторы также работают на базе программы ADMIXTURE. Автор калькулятора задает количество предковых компонентов (K), закладывает в него выборки по популяциям и режим работы. На базе этих данных программой создается файл аллельных частот для каждого компонента. В дальнейшем любой желающий, используя этот файл и некоторые дополнительные утилиты, может «спроецировать» имеющиеся у него широкогеномные данные, свои собственные или чужие, на модель и получить «раскладку» своего генома по предковым компонентам.
Поскольку любой «этнокалькулятор» работает на основе собранной базы данных в референсных популяциях, от качества этой базы данных, в первую очередь, зависит качество его работы. Основная претензия ученых (популяционных генетиков) к тестированию, которое проводится в некоторых коммерческих компаниях, состоит именно в том, что они не раскрывают методику своего тестирования и непонятно, на каких базах данных делают свои оценки. Зарубежные «этнокалькуляторы» не содержат в своей основе репрезентативные выборки по популяциям коренных народов нашей страны, поэтому заведомо не могут давать правильные результаты, а данные по таким выборкам еще находятся в стадии накопления.
Работы, связанные с решением проблемы определения этно-географического происхождения индивида по его ДНК, встречаются и в криминалистических научных журналах. Одна группа исследований посвящена сравнению баз маркеров, использующихся для ДНК-идентификации в криминалистике, в разных популяциях мира. Так, например, в работе [Песик и др., 2014] впервые по панели маркеров, используемых при ДНК-идентификации, охарактеризован широкий круг сельских русских популяций Европейской части России. Изучено коренное население Архангельской, Белгородской, Воронежской, Курской, Ростовской, Рязанской и Орловской областей. Создана база данных по частотам маркеров, определено распределение частот аллелей и генотипов в семи сельских популяциях, проанализированы генетические взаимоотношения между популяциями. Сравнение изученных популяций как друг с другом (по всем 19 маркерам), так и с опубликованными данными по городским русским популяциям (по 12 маркерам), используемым в большинстве работ) не выявило достоверных различий в частотах аллелей, хотя на графиках многомерного шкалирования проявляются отличия генофонда северных русских от основной массы русских  популяций (Фигура 7). Как пишут авторы, полученный результат о почти полной идентичности генофондов различных региональных русских популяций по «криминалистической» панели маркеров создает возможность использования при ДНК-идентификации единой референсной базы данных для всех русских популяций по всей территории России. При этом авторы подчеркивают, что эта база применима только к «русским», а не «российским» популяциям, и ее применение для других народов России может привести к существенным искажениям результатов.
Существуют работы и в другом направлении - поиск ДНК-маркеров, информативных для определения происхождения индивида, с перспективами применения в криминалистике. Они обозначаются аббревиатурой AIMs (AncestryInformativeMarkers), причем среди них есть как SNP маркеры (AIM-SNPs), так и STR маркеры (AIM-STRs). Панели таких специализированных на определение этно-географического происхождения маркеров специалисты разрабатывают на базе исследований популяций из различных регионов мира, причем среди них есть панели как SNP, так и STR-маркеров.
Разные панели маркеров отличаются по трем основным параметрам:
1) по типу маркеров - SNP и STR панели;
2) по количеству маркеров;
3) по уровню, на котором они адекватно определяют происхождение индивидов (то есть разрешающей способности). Постепенно от уровня различения континентальных групп популяций исследователи переходят на уровень различения отдельных популяций в пределах регионов.
Например, исследование информативных маркеров происхождения (AIMs) для различения популяционной принадлежности индивидов на уровне континентов проводили авторы работы [Kosoy et al., 2009]. Они генотипировали 825 человек с разным происхождением, включая Европу, Восточную Азию, Америку (америнды), Африку, Южную Азию, Мексику и Пуэрто-Рико.
Исходя из признаков информативности и широкогеномного распределения, они проводили генотипирование на двух платформах (TaqMan и Illumina). Исследователи создали и протестировали полную панель из 128 информативных по происхождению (AIM – SNP) маркеров и панели из неполных наборов: 96, 64 и 24 маркера. На уровне континентов панели на 128, 96 и 64 маркера работали одинаково эффективно (Фигура 8), они успешно распределяли индивидов на основании их происхождения на четыре континентальные группы (Европа, Америка (америнды), Западная Африка и Восточная Азия). При использовании 24-маркерной панели точность определения континентальной группы снижалась.
Для определения разных по происхождению компонентов в геноме индивидов смешанных популяций эффективным оказалось использование 128-маркерной панели.  Она показала адекватность для определения западноафриканского компонента в AFA, европейского компонента у пуэрториканцев, компонента америндов у американских мексиканцев и мексиканцев. В работе делается вывод о том, что исследователи могут использовать TaqMan зонды для выбранных AIMs как простой и относительно недорогой инструмент для контроля различий в происхождении на уровне континентов при проведении исследований на этнически различных популяциях. Для различения индивидов из несмешанных популяций на уровне отдельных континентов практически одинаково эффективны панели на 128, 96 и 64 AIM – SNP маркеров. Для исследования индивидов из популяций смешанного происхождения предпочтительнее использование 128-маркерной панели.
Итак, в международной научной литературе накоплено немало данных об этнических генофондах и даже разработаны некоторые криминалистические панели, оптимизированные для разных стран или крупных регионов мира. Тем не менее для востребованной в отечественной криминалистике задачи определения происхождения индивида из популяций мультиэтничной России и сопредельных стран надежной и экономически приемлемой (для массового тестирования) системы до сих пор не известно. Отсюда естественным образом вытекает потребность в разработке нового, более точного способа прогнозирования этно-географического происхождения индивида по его ДНК, разработанного на основе изучения генофонда народонаселения России и оптимизированного для применения на российских популяциях.
Сущность изобретения
Настоящее изобретение основано на панели генетических маркеров, позволяющих определить этно-географическую группу и территорию, из которых происходят предки индивидуума, и способе определения этно-географической группы и территории.
1. Более конкретно, изобретение относится к способу определения этно-географической группы населения, из которой происходят предки индивидуума, на основе генотипа по однонуклеотидным полиморфизмам в образце, полученном от индивидуума, где полиморфизмы определены в таблице 1, и вероятность происхождения предков индивидуума из этно-географической группы населения вычисляется по формуле:
, ,
i,
I
j
j
где  и  определены в таблице 1 для каждого полиморфизма где представляет собой целое число от 1 до 5229,и для каждой из этно-географических групп  , где  представляет собой целое число от 1 до 29,
i
– генотип образца по полиморфизму  в формате 0, 1 или 2, где 0 указывается для гомозиготного генотипа по референсному аллелю, 1 указывается для гетерозиготного генотипа, а 2 указывается для гомозиготного генотипа по альтернативному аллелю,
j
при этом этно-географическую группу  населения, из которой происходят предки индивидуума, определяют на основании максимального полученного значения .
Таким образом, после расчета значений вероятности происхождения из каждой из 29 этно-географических групп производится определение наибольшего из полученных значений вероятности, и индивид признается происходящим из этно-географической группы, для которой получено это значение.
2. Изобретение также относится к способу определения территории происхождения предков индивидуума, основанному на картографировании вероятности происхождения, где вероятность происхождения предков индивидуума для каждой этно-географической группы вычисляется по формуле:
, ,
i
I
j
j
I
где  и  определены в таблице 1 для каждого полиморфизма , где представляет собой целое число от 1 до 5229, и каждой из этно-географических групп, где представляет собой целое число от 1 до 29,  – генотип данного образца по полиморфизму в формате 0, 1 или 2, где 0 указывается для гомозиготного генотипа по референсному аллелю, 1 указывается для гетерозиготного генотипа, а 2 указывается для гомозиготного генотипа по альтернативному аллелю,
и по полученным 29 значениям  вычисляется сеточная цифровая модель вероятности происхождения как функции географических координат по интерполяционной формуле:
,  ,
G
R
j
где  – значение вероятности происхождения в точке  с географическими координатами (), () – географические координаты точек привязки этно-географических групп,  – расстояние между точками  и  на сфере радиуса =6371 км, суммирование ведется по всем этно-географическим группам , точки привязки которых отстоят от точки G не более чем на 1000 км,
значения цифровой модели вычисляются в узлах квадратной сетки в проекции карты, координаты узлов сетки перед подстановкой в интерполяционную формулу пересчитываются в долготу и широту по обратным формулам проекции,
территория, занятая на карте узлами цифровой модели со значениями  выше 0,2, признается территорией происхождения индивида,
и цифровая модель визуализируется на географической карте способом изолиний с послойной окраской, при этом территории происхождения соответствует зона на карте, ограниченная изолинией со значением 0,2.
Таким образом, если способ, описанный в п. 1, позволяет определить наименование этно-географической группы, из которой происходит индивид, то для определения географической территории, с которой происходит индивид, производится картографирование по способу, описанному в п.2.
Для картографирования на первом этапе вычисляется цифровая модель. Для этого из значений вероятностей для 29 этно-географических групп вычисляются значения вероятности для каждого узла равномерной сетки, покрывающей территорию России и сопредельных государств. Как этно-географические популяции, используемые для расчета, так и узлы, для которых проводится расчет, имеют географическую привязку, то есть географические координаты, а именно широту и долготу, обозначаемые (), поэтому каждый узел цифровой модели соотносится с определенной точкой географического пространства по правилам картографической проекции.
На втором этапе территория, занятая узлами цифровой модели карты со значениями  выше 0,3 признается территорией происхождения индивида.
На третьем этапе цифровая модель визуализируется в виде географической карты. На карте территории, соответствующие различным интервалам значений вероятности, окрашиваются различным цветом и разделяются изолиниями. Зона на карте, ограниченная изолинией со значением 0,2, показывает территорию происхождения индивидуума.
3. Также изобретение относится к панели для определения происхождения из одной из этно-географических групп, включающей в себя 5229 однонуклеотидных полиморфизмов (Фигура 2), где указанные однонуклеотидные полиморфизмы определяют в образце ДНК, полученном от индивидуума.
Под термином «Россия и сопредельные страны» понимается совокупность стран Северной Евразии, с территории которых происходит большинство лиц, проживающих в России, а именно следующие государства: Азербайджан, Армения, Белоруссия, Грузия,  Молдавия, Монголия, Казахстан, Киргизия, Россия, Таджикистан, Узбекистан, Украина. Включение стран ближнего зарубежья (сопредельных стран) важно и должно учитываться для целей практического применения изобретения в практике криминалистических лабораторий России, поскольку в их работе производится анализ ДНК любых лиц, находящихся на территории России, включая как граждан России, так и граждан указанных стран – основных источников миграций.
Под термином «этно-географические группы» (ЭГГ) понимаются группы коренного населения, совокупность которых охватывает все население России и сопредельных стран таким образом, что каждая ЭГГ генетически гомогенна внутри себя, но при этом обладает генофондом, отличающимся от других ЭГГ. Указание и на этническую, и на географическую привязку в термине «этно-географическая» отражает тот факт, что различия между популяциями человека (т.е. существующими более одного поколения группами людей, заключающих большинство браков  в своих пределах) определяются в наибольшей степени двумя факторами: географическим фактором (коренное населения разных территорий всегда генетически различно) и этническим фактором (группы населения разной этнической идентичности, даже проживающие на одной территории, генетически различаются).
Задачей изобретения является разработка нового способа определения этно-географической группы, из которой происходят предки индивидуума, с выявлением основных этно-географических групп населения России и сопредельных стран – основных источников миграций.
Достигаемым техническим результатом является возможность определения этно-географической группы, из которой происходят предки индивидуума, за счет оценки комплекса генетических данных по геномным сайтам, высоко информативным для дифференцировки генетически различающихся групп населения России и сопредельных стран.
Техническая проблема, решаемая настоящим изобретением, заключается в расширении арсенала технических средств, которые используются для определения происхождения индивидуума.
Таким образом, изобретение состоит в обнаружении генетических полиморфизмов (маркеров), ассоциированных с этно-географической группой происхождения человека, формировании панели (набора) таких маркеров, совместное использование которых позволяет определить этно-географическую группу происхождения и территорию происхождения, и в разработке метода определения происхождения по генотипам, определенным по данной панели маркеров, причем данная панель и метод разработаны на основе изучения российских популяций и оптимизированы для населения России.
Краткое описание чертежей
0
ij
На фиг.1 представлена Таблица 1: Значения коэффициентов θ и θ.
На фиг.2 представлена панель генетических полиморфизмов (генетических маркеров). Каждый маркер однозначно идентифицируется номером хромосомы, номером позиции на хромосоме и нуклеотидом альтернативного аллеля, для каждого маркера эти параметры указаны через знак подчеркивания. Номера позиций указаны согласно референсному геному человека Build 37.
На фиг.3 представлена структура генофонда Европы по широкогеномным маркерам. Обозначения на фиг.3:
а. Статистическое суммирование генетических данных от 1387 европейцев на основе анализа главных компонент (компоненты PC1 и PC2). Маленькие цветные значки обозначают индивидов, большие цветные значки обозначают средние значения PC1 и PC2 для каждой страны. Цвета показаны на встроенной карте. Оси PC повернуты так, чтобы выявить сходство с географической картой Европы. AL, Албания; AT, Австрия; BA, Босния-Герцеговина; BE, Бельгия; BG, Болгария; CH, Швейцария; CY, Кипр; CZ, Чехия; DE, Германия; DK, Дания; ES, Испания; FI, Финляндия; FR, Франция; GB, Великобритания; GR, Греция; HR, Хорватия; HU, Венгрия; IE, Ирландия; IT, Италия; KS, Косово; LV, Латвия; MK, Македония; NO, Норвегия; NL, Нидерланды; PL, Польша; PT, Португалия; RO, Румыния; RS, Сербия и Черногория; RU, Россия, Sct, Шотландия; SE, Швеция; SI, Словения; SK, Словакия; TR, Турция; UA, Украина; YG, Югославия;
b. Увеличение области вокруг Швейцарии с выделением языковой дифференциации;
с. Генетическое сходство в сопоставлении с географическими расстояниями. Средняя генетическая корреляция между парами индивидов как функция географического расстояниями между популяциями.
На фиг.4 представлен анализ предковых компонентов по популяциям и субпопуляциям мира методом ADMIXTURE, Анализ выполнен для заданного числа предковых популяций K=9. По оси х представлены индивиды из популяций, рассортированные в соответствии с местом рождения предков (по опросам). Каждый индивид представлен вертикальной линией. Разными цветами обозначены генетические вклады предковых популяций.
На фиг.5 представлены диаграммы, на которых точность предсказания места происхождения в популяциях обозначена ярко-голубым цветом для стран и светло-голубым для регионов в пределах стран. Средняя величина точности предсказания обозначена красной линией.
На фиг.6 показано сравнение точности методов reAdmix, mSpectrum, HAPMIX и LAMP для индивидов с двумя (вверху) и четырьмя (внизу) компонентами в геноме. Первый столбик обозначает реальное соотношение предковых компонентов. Цвета обозначают: красный – Европу, зеленый – Африку, желтый – Америку, голубой – Восточную Азию, белый – неопределенный компонент.
На фиг.7 показано положение семи русских популяций в пространстве двух первых размерностей многомерного шкалирования матрицы генетических расстояний (по 19 STR-маркерам).
На фиг.8 показан анализ популяционной генетической структуры по панели из 128 AIMs. Каждая вертикальная линия представляет одного индивида. Показаны группа популяций: европейские американцы (EURA, 188 человек), Западная Африка (AFR, 98 человек), америнды (AMI, 88 человек), Восточная Азия (105 человек), Южная Азия (SAS, 64 человека), афроамериканцы (88 человек), пуэрториканцы (PRA, 28 человек), мексиканские американцы (MAM, 40 человек) и мексиканцы (MXN, 26 человек). А, В – анализ с использованием панели на 128 маркеров; C, D – анализ с использованием панели на 64 маркера.
На фиг.9 показан пример распределения по длинам волн средних интенсивностей флуоресценции зондов для одного полиморфизма и 48 образцов. Каждая точка на графике – единичный образец, каждый кластер сформирован из образцов, обладающих одинаковым генотипом (слева направо: гомозигота, гетерозигота, гомозигота).
На фиг.10 показаны TOP, BOTTOM, FORWARD схемы на примере двух различных полиморфизмов, находящихся внутри экзонов различных генов.
На фиг.11 показан график первой и второй главных компонент изменчивости генофонда населения России и сопредельных стран с разделением на группы методом K-means.
На фиг.12 показан график главных компонент изменчивости для региона Кавказа.
На фиг.13 показан график главных компонент изменчивости с итоговым разделением на этно-географические группы.
На фиг.14 показана итоговая карта выделенных этно-географических групп.
Цифрами в кружках на карте подписаны номера групп (совпадают с номерами в таблице 3).
На фиг.15 представлены генотипы образцов №1, №2, №3, использованные в описании примеров использования изобретения.
На фиг.16 показана территория происхождения образца №1, определенная по способу, предложенному в изобретении.
На фиг.17 показана территория происхождения образца №1, определенная по способу, предложенному в изобретении.
На фиг.18 показана территория происхождения образца №1, определенная по способу, предложенному в изобретении.
Осуществление изобретения
Для разработки изобретения были проведены следующие работы:
Этап 1: Формирование коллекций образцов ДНК, представляющих основных группы народонаселения России и стран – основных источников миграций;
Этап 2: Генотипирование панели из 4,5 млн. генетических маркеров;
Этап 3: Выявление оптимального набора этно-географических групп, для которых может производиться определение происхождения;
Этап 4: Идентификация высокоинформативных маркеров происхождения и построение алгоритма определения этно-географической группы и территории.
Ниже описаны каждый из этих этапов, после этого следуют примеры применения предложенного способа.
Этап 1: Формирование коллекций образцов ДНК, представляющих основных группы народонаселения России и стран – основных источников миграций
Для успешного применения системы прогнозирования в географическом масштабе России и сопредельных стран потребовалось охватить практически всё популяционное разнообразие данного региона, сочетая спектр этнических и субэтнических групп с минимально требуемым размером каждой выборки (а также - с условиями отбора образцов в эти выборки, о чём подробнее изложено ниже). Также, для итогового выбора только 5 тысяч наиболее информативных SNP-маркеров потребовалось генотипировать почти в 1000 раз больше таких полиморфизмов. Для столь масштабного молекулярно-генетического анализа потребовалась не менее масштабная подготовка коллекции образцов ДНК.
Коллекция включила более 2500 образцов ДНК из более 150 популяций России и сопредельных стран. В каждом случае учитывались этническая и географическая принадлежность индивида на глубину не менее трёх поколений, для образцов ДНК были проведены измерения концентрации и объема ДНК с последующим включением в коллекцию только аликвот, соответствующих требуемым параметрам; в процессе работы весь массив информация вносился в базу данных. Подробнее эти этапы описаны ниже.
В коллекцию образцов были включены популяции, представляющие все крупные географические регионы Северной Евразии: Восточная Европа и Урал, Кавказ и Центральная Азия, Сибирь и Дальний Восток. Внутри каждой крупной региональной группы учитывалось этническое и субэтническое разнообразие популяций, показанное в широком спектре предшествующих исследований. Например, располагая данными о структуре генофонда Кавказа (подразделение генофонда на четыре основных кластера: Западный, Центральный, Восточный и Закавказье, что чётко показано в [Balanovskyetal., 2011; Jeongetal., 2019]), в коллекцию обязательно включали популяции каждого из этих кластеров, причем из каждого кластера – по несколько этнических групп, а где известно субэтническое разнообразие – учитывали и его (чтобы охватить максимальное число предковых линий, представленных в геноме каждого индивида). В Западной Сибири, формируя коллекции образцов из популяций сибирских татар, принимали во внимание показанные в [Агджоян и др., 2016] значительные отличия генофондов пяти субэтнических групп тоболо-иртышских татар; в Центральной Азии и Южной Сибири – учитывали родоплеменные подразделения казахов, монголов, бурят, и т.д. Таким образом было охвачено генетическое разнообразие населения Северной Евразии на различных иерархических уровнях – от регионального до родоплеменного.
С другой стороны, для репрезентативности выборок (соответственно, и статистической надежности будущих расчетов частот маркеров) учитывались и количественные показатели. Например, популяция этнической группы была представлена не менее 10 индивидами, а в случае известной подразделённой этнического генофонда – не менее 5 индивидов для каждой субэтнической или родоплеменной группы.
При отборе образцов строго соблюдались три необходимых условия формирования коллекции: «биоэтика», «три поколения», «неродственность».
«Биоэтика»: для каждого образца существует письменное информированное согласие на предоставление биологического образца для популяционно-генетических исследований, на хранение ДНК и индивидуальной информации по форме, одобренной Этической комиссией ФГБУН МГНЦ. Информированные согласия были получены исследователями, сформировавшими первоначальные коллекции по итогам экспедиционных обследований популяций коренного населения Северной Евразии [Балановская и др., 2016].
«Три поколения». Индивиды, чьи образцы представляют данную популяцию, соответствуют международному правилу «трех поколений» (и оба его родителя, и все четверо его бабушек и дедов родились в данной популяции и все четверо относили себя к данной этнической или субэтнической группе). Данные для отбора по описанным выше параметрам были получены в ходе анкетирования обследуемых лиц исследователями во время экспедиционных обследований популяций [Балановская и др., 2016].
«Неродственность». Индивиды, чьи образцы были включены в коллекцию, не являются родственными друг другу до третьей степени родства, что позволяет избежать уменьшения генетически эффективного размера выборки и смещения частоты генетических маркеров.
Описанные подходы полностью согласуются с логикой поискового исследования, конечной целью которого представляется отбор сравнительно небольшого числа (5 тыс. из 4,5 млн.) наиболее информативных SNP-маркеров для определения происхождения из населения столь этнически пестрой территории, как Северная Евразия. Однако генотипирование 4,5 млн. геномных точек для одного образца требует достаточное количество и качество материала (ДНК). И эти требования привели к включению в план работ ещё одной стадии отбора образцов в коллекцию по молекулярным критериям качества и количества ДНК.
Проведение этой стадии отбора образцов в коллекцию включало оценку качества образцов ДНК, для чего были проведены измерения концентрации ДНК и объема растворов. По итогам измерений в состав коллекции были отобраны образцы с концентрацией не менее 40 нг/мкл и объемом аликвот не менее 10 мкл.
Образцы ДНК были помещены на хранение в SBS формат (штативы или планшеты, содержащие 96 образцов, 8×12 рядов) для обеспечения возможности их обработки с помощью роботизированных станций аликвотирования или пробоподготовки).
Подготовленная таким образом коллекция позволила провести генотипирование 4,5 млн. SNP-маркеров.
Этап 2: Генотипирование панели из 4,5 млн. генетических маркеров
Генотипирование маркеров проведено с использованием панели Infinium Omni5Exome-4 v1.3 BeadChip (Illumina, США). Данная панель является наиболее обширной из представленных на рынке и наиболее современной, и поэтому была выбрана для генотипирования.
Данный чип используется в полногеномных исследованиях и содержит в своем составе набор зондов для детекции 4 559 465 генетических маркеров человека. Среди них 544 817 экзонных и 1 886 632 интронных маркера, соответствующих базе данных RefSeq; 6 040 нонсенс-маркеров, 275 657 миссенс-маркеров, 76 120 маркеров синонимичных замен и 186 900 молчащих маркеров, соответствующих базе данных GenomeBrowser; 408 митохондриальных маркеров, 4 766 маркеров инсерций и делеций, а также маркеры половых хромосом: 118 204 для X-хромосомы, 2 574 для Y-хромосомы и 5 364 для их псевдоаутосомальных и гомологичных регионов, в соответствии с базой данных NCBI GenomeReferenceConsortium GRCh37.
Генотипирование образцов ДНК осуществляли в соответствии с предварительно отработанным протоколом. Пробоподготовка была проведена в 8 мкл ДНК с концентрацией 50 нг/мкл. Нанесение образцов ДНК на биочипы было проведено с использованием роботизированной системы TecanFreedomEvo (Tecan, Швейцария) и программного обеспечения Hyb Multi-BC2 AMP5. Гибридизация была проведена в камере HybChamber, промывка – с помощью Te-Flow LCG Flow-ThroughChambers (Illumina, США). После промывки и обработки (также с использованием роботизированной системы TecanFreedomEvo) было проведено сканирование результатов с использованием сканера биочипов iScan (Illumina, США).
Анализ полученных данных (файлы формата *.idat) проведен в программе GenomeStudio v2011.1 для проведения первичного анализа и оценки качества полученных данных. Для исследованных образцов показатель CallRate (доля успешно генотипированных маркеров) был определён на уровне не менее 99%, что говорит о пригодности полученных результатов для дальнейшей обработки.
При работе с чипами высокоплотного кварцевого нанесения (технология BeadChips от Illumina) статус полиморфизма определяется по интенсивности свечения флуорофора, ассоциированного с аллельным вариантом. При этом следует учитывать, что уровни интенсивности нормализуются как по всем бусинам чипа (порядка 30 повторов), так и по всем образцам чипа или запуска (не менее 12 образцов).
Обычно полученные интенсивности на графике группируются в 3 кластера: 2 гомозиготы и гетерозигота (Фигура 9). Для каждого полиморфизма вычисляют вероятность ошибки при отклонении нулевой гипотезы (P-значение, P-value) классификации в ожидаемый кластер (альтернативная гипотеза) против попадания в какой-либо другой (нулевая гипотеза). Если полученное значение P-value>0,05, то считается, что точка не попадает в указанный кластер и данный полиморфизм у анализируемого образца считается определившимся некорректно. Для каждого прошедшего первичный контроль полиморфизма вычисляется внутренний параметр GC-score. Эта величина является непосредственной характеристикой яркости флуоресценции зонда и измеряется как относительная интенсивность в сравнении со специальной технической бусиной, яркость которой всегда максимальная, так как все зонды этой бусины флуоресцируют. Поэтому, GC-score может принимать значения от 0 (свечение отсутствует) до 1 (все зонды всех бусин одного полиморфизма светятся максимально). Производителем установлен минимальный уровень флуоресценции в 0,15 единиц GC-score для того, чтобы учитывать наличие соответствующего сигнала.
Для дополнительного повышения точности генотипирования была изучена и снята проблема ошибочной записи генотипа полиморфизма по комплементарной нити ДНК. Эта проблема типична для ДНК-чипов, так как для разработчиков в зависимости от нуклеотидного контекста вокруг полиморфизма иногда удобнее подбирать зонд, комплементарный обратной цепи ДНК.
Научным сообществом принято, что та цепь, которая со стороны начала нумерации нуклеотидов в хромосомах имеет 5’ конец называется TOP (верхняя), а комплементарная ей – BOTTOM (нижняя) (Фигура 10). Но, к сожалению, не всегда удобно указывать значения нуклеотидов по схеме TOP/BOTTOM, поэтому были разработаны альтернативные схемы обозначения. Схема FORWARD подразумевает, что нуклеотид будет записан согласно расположению гена, последовательность которого была отсеквенирована. Таким образом, при схеме FORWARD требуется указать ген и интересующее нас положение, а эту информацию, зная на какой из цепей TOP или BOTTOM расположен ген, можно легко перевести в более простые форматы. К сожалению, до появления 38-й версии сборки человеческого генома (HumanGenome 38, hg38) не существовало единого стандарта на запись нуклеотидных состояний по всему миру, поэтому различные компании-производители для каждого из форматов (TOP, FORWARD) использовали геномную последовательность, которая имела координатную структуру генома hg19 (или hg18), но информация об "основной" цепи различалась. Поэтому очень частой является ситуация несоответствия результатов чип-гибридизации и секвенирования именно в части направления цепей.
Несмотря на то, что в аннотационном файле к ДНК-чипу указана правильная ориентация для всех зондов, зачастую ошибочное позиционирование не определяется в ходе предварительной подготовки данных и вносит погрешность на последующих этапах анализа. Для исключения этого фактора в качестве исходных данных для биоинформатической обработки использовали текстовые файлы, содержащие в каждой из строк наименование зонда, наименование образца, определившийся генотип, также показатели GC-score, оценивающие эффективность работы зонда. Известно, что в технологии чип-гибридизации Illumina зонды могут быть комплементарны любой из цепей ДНК. Для стандартизации получаемых результатов с данными секвенирования в качестве референсной последовательности была использована сборка генома человека HumanGenome 19 (hg19, GRCh37), размещенная в базе данных NCBI. Компания Illumina для каждого из чипов поставляет текстовый файл, который содержит последовательности зондов и их координаты в той версии генома, в которой они были разработаны. Для исключения описанных погрешность были использованы оригинальные  алгоритм и программный продукт, который на основании предоставленного от компании-производителя манифест-файла (содержащего списки зондов, находящихся на стекле чипа) формирует список позиций и определяет, возможно ли доверять производителю по соответствию зонд-координата и на каком из направлений находится прочитанный нуклеотид (+ или -). На первом шаге производится поиск последовательностей зондов, полученных от производителя, по их координатам в референсном геноме (hg19) как в прямом, так и в комплементарно-обратном вариантах. В качестве критерия соответствия было установлено, что допустимое максимальное количество несоответствующих нуклеотидов при выравнивании не может превышать 2 на зонд. После первого этапа формировался список зондов, которые были обнаружены находящимися в соответствии с манифест-файлом производителя, вместе с указанием цепи их нахождения, а также список тех зондов, координаты которых не соответствуют заявленным.
На втором этапе зонды, координаты которых не соответствовали сборке hg19, с помощью программного пакета blastn выравнивались на референсный геном. Считалось, что поиск был произведен успешно в том случае, если при таком выравнивании с учетом не более 2-х несоответствий находился только один уникальный фрагмент ДНК. Результаты такой проверки представлены в таблице 2.
Таблица 2. Результаты анализа структуры зондов чипов компании Illumina, использованных в данной работе
Чип
InfiniumOmni 5 Exome 4
Общее число полиморфизмов
4559465
Число корректно определяемых полиморфизмов
3803518
Число тестируемых геномных состояний, координаты для которых указаны верно, а цепь - ошибочно
710353
Число тестируемых геномных состояний, координаты которых определены путем выравнивания blastn
13678
Число геномных состояний, исключенных из анализа по различным причинам
31916
Сборки референсного генома, на основе которой производился дизайн зондов
hg19 / hg38
В итоге, если обработка результатов генотипирования с использованием чипа InfiniumOmni 5 Exome 4 была бы проведена по стандартным протоколам Illumina, и в качестве источника данных был бы выбран файл-отчет с указанием генотипов в виде TOP (как подразумевается компанией Illumina), то в 74243 геномных состояниях из 710353 был бы неточно определен генотип без возможности последующего исправления методами биоинформатики. Указанные полиморфизмы являются транзициями (A/T, G/C), которые, согласно документации производителя, находятся на TOP цепи их сборки генома hg19, тогда как они расположены на BOTTOM цепи основной версии генома hg19 NCBI. Кроме того, у 636110 из 710353 полиморфизмов, являющихся трансверсиями, возникла бы ситуация несоответствия нуклеотидов результатам секвенирования из-за того, что они также находятся на BOTTOM цепи, но совпадают с направлением гена (т.е. FORWARD). Нуклеотиды таких полиморфизмов находятся в согласии с базами данных (например, dbSNP), но являются комплементарными тем парам, которые могут быть обнаружены NGS-методами. Для таких состояний потребовалось бы написать программу-обработчик, которая бы определяла такие состояния и автоматически заменяла значения генотипов на комплементарные пары у всех образцов.
Кроме того, со времени разработки чипа InfiniumOmni 5 Exome 4 v1.3 (2016 год) некоторые полиморфизмы изменили свой индекс rs и на следующих чипах находятся с новой его версией. Именно такие зонды не находились по координатам, но легко обнаруживались сплошным поиском (blastn), и, в случае если бы интеграция данных производилась по индексам rs, такие точки были бы утеряны. В процессе анализа было обнаружено 1985 таких полиморфизма.
Таким образом, следует заметить, что при использовании штатного программного обеспечения (GenomeStudio версии 2.0), процент неправильно идентифицированных полиморфизмов может достигать 2%. При использовании высокоплотных ДНК-чипов, таких как InfiniumOmni 5 Exome 4, в эту категорию попадут более 50000 полиморфизмов. Поэтому, применение разработанного протокола обработки данных генотипирования на ДНК-чипах явилось важным этапом, дополнительно повысившим точность генотипирования.
Далее были подготовлены стандартные файлы формата среды Plink: .lgen, .map и .fam. Файл .lgen содержал 5 полей: идентификатор семьи, идентификатор индивида, идентификатор варианта, аллель 1, аллель 2. Файл .map содержал 4 поля: номер хромосомы, идентификатор варианта, расстояние сцепления (в сантиморганах), координата варианта в хромосоме. Для создания .map файла был использован манифест-файл соответствующего чипа Illumina. Файл .fam содержал 6 полей: идентификатор семьи, идентификатор индивида, идентификатор отца, идентификатор матери, пол, фенотип.
Далее набор файлов .lgen, .map, .fam был преобразован в набор файлов .bed (в бинарном формате), .bim, .fam в среде Plink, команда lfile. С помощью команд geno и mind в среде Plink были удалены варианты и образцы с уровнем пропусков генотипирования выше 10%. С помощью команд list-duplicate-vars и exclude в среде Plink были удалены дубликаты.
Таким образом, в соответствии с разработанным алгоритмом этап предварительной биоинформатической предобработки данных состоял из следующих шагов:
1) Формирование выгрузки из GenomeStudiov2.0 в формате txt с полями snp_id, user_id, allele_top_1, allele_top_2, gc_score;
2) Преобразование аллелей в сформированном файле с использованием собственного файла-схемы. На выходе формат файла сохраняется;
3) Замена аллелей на N в случае, если gc_score< 0.35;
4) Формирование собственным скриптом ped и map файлов на основании полученной модифицированной выгрузки из GenomeStudio;
5) Запуск Plink и преобразование текстовых ped и map в бинарные bed, bim, fam.
Результатомгенотипирования явился массив данных в формате PLINK, содержащий генотипы 1883 образцов по 4,5 млн. генетических маркеров.
Этап 3: Выявление оптимального набора этно-географических групп, для которых может производиться определение происхождения.
Население России характеризуется исключительно высокой гетерогенностью не только в этническом и лингвистическом плане, но и в плане разнообразия генофонда [Triska et al., 2019; Jeong et al., 2019]. В частности, велики генетические различия населения Центральной России не только от коренных народов Сибири, но и от народов Кавказа, а также от народов Волго-Уральского региона. Таким образом, базовой задачей является выделение в народонаселении России набора конкретных популяций (этно-географических групп), которые будут включены в разрабатываемый способ определения происхождения, стандартизованную панель скрининга. Этот набор популяций должен удовлетворять следующим требованиям:
- охватывать население по всей территории России;
- включать территории сопредельных государств, демографически тесно связанных с Россией, в том числе уроженцы которых часто встречаются на территории России;
- на каждой территории быть ориентированным на коренное население;
- каждая выделенная популяция должна быть генетически достаточно гомогенна, чтобы определенные для нее частоты маркеров были репрезентативны для каждого из народов или регионов, входящих в состав популяции;
- Число популяций должно быть минимизировано, чтобы типичный объем выборки составлял не менее 25 образцов, поскольку при меньших объемах выборок (при большом числе признаков, как в нашем случае) машинное обучение классифицирующих алгоритмов малоэффективно.
Для формирования набора популяций ключевыми являются требования (3) и (4). Их направленность противоположная: требование гомогенности увеличивает число выделяемых групп, а требование объемов выборки их уменьшает. Поэтому нужен математически корректный анализ для выявления оптимального баланса между ними. На практике, принятой в международных исследованиях (например, проекте “1000 Genomes”) стандартом минимального рекомендуемого объема выборки является 50 человек (100 хромосом), хотя нередко используют и выборки меньшего объема. Полученные к настоящему времени широкогеномные массивы данных по российским популяциям включают от 1 до 3 тысяч образцов [Triska et al., 2019; Jeong et al., 2019]. Проект “Российские геномы” также ставит своей целью секвенировать 2-3 тысячи полных геномов представителей коренных народов России. Такие объемы доступных общих выборок (2-3 тысячи образцов) и требования к среднему объему одной выборки (около 50 образцов) определяют, что максимальное число выделяемых популяций может быть 40-60. При этом для удобства использования и для дополнительного увеличения объемов выборок чисто популяций желательно минимизировать.
Основные проведенные к настоящему времени геномные исследования российских популяций включают порядка 100-200 популяций, изученная каждая по выборке всего в несколько образцов. Таким образом, возникает задача объединять несколько сотен популяций в 3-6 десятков более крупных групп, которые были бы гомогенны внутри себя, но тем не менее отражали бы основное разнообразие генофонда различных народов и регионов.
Поэтому целью данного этапа работы являлось разделение народонаселения России и сопредельных стран на группы популяций, удовлетворяющие вышеперечисленным условиям, а также имеющие репрезентативную выборку в существующих данных и биобанках.
Если переформулировать цель исследования на языке математических моделей, возникнет задача кластеризации популяций, представленных в имеющихся широкогеномных массивах данных, на основе генетического разнообразия популяций. Существует много различных подходов к кластеризации, в том числе специально для генетических данных [Alhusain et al., 2018]. Методы глобально делятся на два подхода - параметрический и непараметрический. Параметрические методы опираются на статистические модели, в основе которых лежат несколько генетических допущений. К таким методам относятся STRUCTURE, ADMIXTURE и другие. Однако данные методы являются ресурсоемкими и поэтому трудно применимы к разностороннему анализу больших датасетов с большим количеством снипов.
К непараметрическим методам относятся методы либо на основе уменьшения размерности и дальнейшей кластеризации [Liu et al., 2006], либо на основе расчета попарных расстояний между образцами и дальнейшей кластеризации. Непараметрические методы менее ресурсоемки и не требуют допущений.
Для данной работы был выбран непараметрический подход с уменьшением размерности методом главных компонент с применением программного пакета EIGENSTRAT\smartpca [Patterson et al., 2006]. В качестве метода кластеризации мы использовали K-means [Hartigan et al., 1979] на главных компонентах, который дает результаты не хуже, чем STRUCTURE [Lee et al., 2009].
Чтобы избежать терминологической путаницы при использовании понятий “население”, “народонаселение”, “народ”, “субэтническая группа”, “этнографическая группа”, “регион” и т.д., мы используем термин “этно-географические группы” (ЭГГ) для обозначения групп популяций коренного населения, совокупность которых охватывает все население изучаемого крупного региона таким образом, что каждая ЭГГ генетически сравнительно гомогенна внутри себя, но при этом обладает генофондом, отличающимся от других ЭГГ.
Для работы использовались данные генотипирования 4,5 млн. генетических полиморфизмов (маркеров) в различных популяциях Северной Евразии, полученные на предыдущем этапе исследования (описан выше). Анализ проведен на массиве данных, полученном при генотипировании панели производства Illumina (InfiniumOmniExomeBeadChipKit, Illumina; США). При этом выделяемые группы зависят от структуры генофонда самого населения, а не от используемой панели маркеров, поэтому очень близкие результаты могут быть получены и при использовании других широкогеномных панелей (например HumanOrigin производства Affimetrix, США) или полногеномногосеквенирования.
Общий объем используемого в анализе массива данных: 1883 образца, представляющие 247 популяций России и сопредельных стран. Минимальным размером ЭГГ был выбран порог в 25 образцов, оптимальным признавалась ЭГГ размером от 50 до 100 образцов. Таким образом, задачей было разделение данных 1883 образцов на гомогенные ЭГГ размером не менее 25 образцов каждая.
При биоинформатическом анализе для фильтрации использовалось программное обеспечение PLINK [Changetal., 2015]. Для первичной фильтрации применялась команда mind 0.1. Для фильтрации перед построением применением метода главных компонент последовательно применялись команды geno 0.05, maf 0.01, mind 0.1, indep-pairwise 1500 150 0.2. Для фильтрации перед расчетом значений FST последовательно применялись команды geno 0.1, maf 0.01, mind 0.05, indep-pairwise 1500 150 0.5.
Для поиска родственников ближе второй ступени применялось программное обеспечение KING 2.2.4 [Manichaikul et al., 2010] в режиме related с настройками по умолчанию.
Для расчета значений главных компонент всех образцов применялась программа smartpca, входящая в пакет eigensoft. Данные конвертировались из формата plink (bed-bim-fam) в формат eigensoft (eigenstratgeno-snp-ind) с помощью программы convertf из того же пакета с параметрами по умолчанию. Конвертированные данные подавались на вход программе smartpca с параметрами по умолчанию, кроме параметра, задающего количество итераций исключения outliers (numoutlieriter). Данный параметр имел значение 3 при построении промежуточных графиков, и значение 0 при построении итоговых. Результаты работы smartpca визуализировались с помощью языка программирования Python 3, в том числе с использованием библиотек pandas, matplotlib и seaborn.
Кластеризация проводилась на результатах работы smartpca на трех первых компонентах. Значения этих компонент для всех образцов считывались библиотекой pandas и подавались на вход методу sklearn.cluster.KMeans [Pedregosa et al., 2011] с фиксированным количеством кластеров и фиксированным генератором псевдослучайных чисел для воспроизводимости результатов. Для итогового деления было выбрано количество кластеров, равное 30.
ST
ST
Для расчета значений генетических расстояний между парами популяций F применялась программа smartpca, входящая в пакет eigensoft. Данные конвертировались из формата plink (bed-bim-fam) в формат eigensoft (eigenstratgeno-snp-ind) с помощью программы convertf из того же пакета с параметрами по умолчанию. Конвертированные данные подавались на вход программе smartpca с параметром “fstonly: YES” для расчета F. Результаты работы smartpca визуализировались с помощью языка программирования Python 3, в том числе с использованием библиотек pandas, matplotlib и seaborn.
Пространственное распределение выделенных этно-географических групп было картографировано в программе GeneGeo [Koshel, 2012] по данным о географических координатах 247 исходных популяций, каждая из которых была отнесена к одной из 29 результирующих ЭГГ. Области на карте, раскрашенные в индивидуальный цвет для каждой группы, получены как объединение ячеек диаграммы Вороного с одинаковым номером группы, построенной на множестве точек, соответствующих исходным популяциям (на карте отмечены звездочками). Расстояния между точками при построении диаграммы Вороного вычислялись на сфере с использованием их географических координат. Важно понимать, что полученные таким образом области условны и границы между ними не следует трактовать как ареалы расселения соответствующих популяций. Однако можно говорить о том, что ареалы содержатся внутри этих областей, что позволяет наглядно представить географические закономерности в их пространственном распределении.
В соответствии с описанным подходов, первым шагом работы было проведение фильтрации данных с помощью программного обеспечения PLINK. Был отфильтрован 51 плохо прочитанный образец (с покрытием менее 90%). Далее были исключены 19 образцов, для которых анализ с помощью программного обеспечения KING показал родство ближе второй степени с другим образцом из выборки. После данной фильтрации осталось 1813 образцов.
Для проведения всех анализов методом главных компонент применялся следующий алгоритм. Изначально выбирались образцы, для которых мы хотели применить метод (либо все образцы, прошедшие фильтрацию, либо только образцы из определенного региона). Затем исключались полиморфизмы, прочитанные у менее чем 95% образцов и полиморфизмы, частота встречаемости минорного аллеля которых была ниже 1%. После этого снова исключались образцы, у которых прочитано менее 90% позиций из уже отфильтрованного набора. После этого исключались сцепленные полиморфизмы с коэффициентом r2 более 0.2. Заключительным этапом построения было применение smartpca с тремя итерациями исключения outliers (кроме итоговых графиков, на которых pca проводился без исключения outliers).
Изначально были рассчитаны главные компоненты по всем отфильтрованным 1813 образцам. Затем были рассчитаны координаты центроидов популяций, как среднее значение по компонентам всех образцов данной популяции. Центроидам был присвоен вес, пропорциональный размеру выборки образцов из данной популяции. К имеющим вес центроидам по первым 3 компонентам был применен метод кластеризации k-means из программного пакета sklearn. Таким образом была выполнена кластеризация, имитирующая кластеризацию по отдельным образцам, однако обязательно относящая популяции к ЭГГ целиком.
Значение количества кластеров k подбиралось исходя из условия на размеры кластеров (не менее 25 образцов в каждом кластере и минимальное количество крупных кластеров). Оптимальным значением k было 30, то есть метод давал разделение популяций на 30 ЭГГ. При данном k из 30 групп всего 8 имели менее 25 образцов, и всего 3 группы имели более 100 образцов. Увеличение k увеличивало количество малых групп, не разделяя крупные, а уменьшение k увеличивало размеры крупных групп, не уменьшая количество малых. Полученное деление на графике первых двух главных компонент приведено на фигуре 11.
Неравномерность разделения далее была исправлена вручную с применением дополнительных построений графиков главных компонент на подмножествах образцов.
Малые группы включали в себя группы состоящие из следующих популяций:  “алеуты”, “болгары, гагаузы и молдаване”, “ительмены”, “ительмены-камчадалы”,  “нанайцы и орочи”, “нивхи и ульчи”, “чукчи”, “эвенки” и “эвены”. Учитывая внутреннее разнообразие данных популяций, объединение “эвенов” с “эвенками”, “нивхов и ульчей” с “нанайцами и орочи”, “ительменов” с “ительменами-камчадалами” и “чукчами” практически не нарушило гомогенность данных групп. Также к группе “ительмены и чукчи” были добавлены “коряки”, так как после объединения этой группы их внутреннее разнообразие было достаточно большим, чтобы включение “коряков” не нарушало гомогенность этой группы.
Малый размер и генетическое разнообразие “болгар, молдаван и гагаузов” и “алеутов” не позволил выделить их в отдельные группы или присоединить к какой-либо существующей группе. Тем не менее группу из болгар и гагаузов стоило бы выделить отдельно, если бы она имела достаточное количество образцов, в то время как молдаване имеют слишком большое внутреннее разнообразие, и поэтому могут относиться как к группе с болгарами и гагаузами, так и к группе с украинцами.
В свою очередь крупные группы включали в себя группы “популяции Кавказа”, “русские” и “буряты, дунганы, калмыки, монголы и хамнигане”. Все эти группы включают в себя достаточно разнообразные субпопуляции, поэтому они были разделены на более мелкие группы на основании графиков главных компонент для субпопуляций из данных групп. В качестве примера приведем график первых двух главных компонент с итоговым разделением на ЭГГ для региона Кавказ (Фигура 12).
В результате группа “буряты, дунганы, калмыки, монголы и хамнигане” была разделена на “калмыки и монголы” и “буряты, хамнигане и якуты” (якуты перемещены сюда, так как в результате этого деления они оказались ближе к этой группе, чем к изначальной “тофалары, тувинцы и якуты”), а популяция “дунгане” исключена, так как сильно отличалась от остальных популяций этих групп. При наличии достаточного количества образцов “дунгане” составили бы отдельную группу.
Группа “русские” была разделена на три группы: “южные русские”, “северные русские” и “русские севера Архангельской области”. После этого разделения популяция водь оказалась ближе к группе “северные русские”, чем к карелам и вепсам, к которым ее отнесла кластеризация k-means изначально, в связи с чем она была перенесена.
Группа “Кавказ” была разделена на четыре группы: “Дагестан”, “западный Кавказ”,  “ингуши и чеченцы” и “осетины”, а также некоторая часть популяций была перенесена в группу “Закавказье и Крым”. Группа “ингуши и чеченцы” включила в себя 23 образца, однако должна была быть больше, так как часть образцов была исключена из нее во время фильтраций. В связи с этим мы оставили ее, как отдельную ЭГГ, хотя формально она не подходит из-за размера выборки.
Также достаточно крупной являлась группа “алтайцы, хакасы и шорцы” (73 образца). При отдельном рассмотрении она отлично разделилась на “хакасы и южные алтайцы” и “шорцы и северные алтайцы”.
Также удмурты и бесермяне были перенесены из группы с мари и чувашами в группу с коми, а саамы исключены из этой группы, как вносящие излишнюю гетерогенность.
В связи с невозможностью определенно отнести к какой-то группе или выделить в отдельную ЭГГ по причине малых объемов выборок также были исключены греки, литовцы, караногайцы, фарси и кумыки.
ST
При распределении популяций по группам в случае спорных моментов об отнесении популяции к одной или другой группе использовались также значения F между популяциями.
Итоговое разделение включило в себя 29 групп. Их состав приведен в таблице 3.
Таблица 3. Итоговый набор этно-географических групп.
№ЭГГ
Этно-географическая группа (ЭГГ)
Народы, включенные в ЭГГ
Кол-во образцов (в использованном массиве данных)
1
Чукчи, коряки и ительмены
коряки, чукчи, ительмены, камчадалы
67
2
Эвенки и эвены
эвены, эвенки
45
3
Нанайцы, нивхи, орочи и ульчи
нанайцы, нивхи, ульчи, орочи
39
4
Буряты и хамнигане
буряты, хамнигане, якуты
57
5
Ханты, манси и ненцы
ненцы, ханты, манси
53
6
Монголы и калмыки
калмыки, монголы
126
7
Тувинцы и тофалары
тофалары, тувинцы, монголы-тувинцы
56
8
Хакасы и Южный Алтай
алтайцы, хакасы
42
9
Шорцы и Северный Алтай
алтайцы, шорцы
35
10
Сибирские татары
татары сибирские
68
11
Казахи, каракалпаки, уйгуры и ногайцы
уйгуры, казахи, каракалпаки, ногайцы астраханские, ногайцы ставропольские
33
12
Киргизы
киргизы
35
13
Таджики, памирцы и ягнобцы
таджики, памирские народы, ягнобцы
72
14
Узбеки и туркмены
узбеки, туркмены
45
15
Башкиры
башкиры
43
16
Татары
татары, татары мишари, татары астраханские, татары кряшены, татары казанские
52
17
Марийцы и чуваши
чуваши, марийцы
53
18
Коми и удмурты
коми-пермяки, удмурты, бесермяне, коми-зыряне
84
19
Карелы и вепсы
карелы, вепсы, карелы тверские
38
20
Северные русские Архангельской области
русские
35
21
Северные русские
русские, ижора, водь
76
22
Южные русские
русские
198
23
Мордва
мордва-мокша, мордва-эрзя, мордва-шокша
40
24
Украинцы
украинцы
79
25
Дагестан
табасаранцы, аварцы, кубачинцы, даргинцы, лакцы, лезгины, рутульцы
68
26
Чеченцы и ингуши
чеченцы, ингуши
28
27
Осетины
осетины
36
28
Западный Кавказ
карачаевцы, абхазы, адыгейцы, кабардинцы, шапсуги, балкарцы, черкесы, абазины
87
29
Закавказье и Крым
армяне, азербайджанцы, караимы, турки-месхетинцы, курды, татары крымские
83
Итоговый график первых двух главных компонент приведен на Фигуре 13. Можно видеть, что выделенные ЭГГ охватывают все генетическое разнообразие населения России сопредельных стран, при этом перекрывание между генетической изменчивостью групп минимальное.
Итак, в результате применения методов биоинформатического анализа к массиву широкогеномных данных о популяциях Северной Евразии, был сформирован набор из 29 этно-географических групп (ЭГГ). Эти группы характеризуются генетической гомогенностью внутри себя, генетически различаются друг от друга, охватывают все генетическое разнообразие населения России и сопредельных стран, то есть соответствуют требованиям, сформулированным в начале нашего исследования. На Фигуре 14 представлена карта расположения выделенных групп.
Построенная карта фактически представляет собой схему районирования территории России и сопредельных стран по признаку генетического разнообразия населения. Важно, что карта получена по результатам наиболее подробного массива данных, охватывающего большое число конкретных групп населения (247 исходных популяций) и основанного на типировании большого числа маркеров, покрывающих весь геном. Поэтому проведенное районирование является если не финальным и не единственно возможным, то приближается к этой цели, поскольку опирается на репрезентативную выборку популяций и маркеров и получено путем применения объективных математических методов. Отметим, что хотя критерий географического соседства не применялся при выделении ЭГГ, практически все выделенные группы имеют целостные, а не разорванные ареалы.  Полученная карта ЭГГ может использоваться при определении происхождения человека по его ДНК в генетической генеалогии и криминалистике.
Этап 4: Идентификация высокоинформативных маркеров происхождения и построение алгоритма определения этно-географической группы и территории
Для отбора наиболее информативных ДНК маркеров, позволяющих определить происхождение (популяционную принадлежность к этно-географическим группам Северной Евразии) были проанализированы индивиды, принадлежащие к различным географическим популяциям и этническим группам. Вторая группа индивидов (из общего числа 2 400 образцов была зарезервирована для построения и тестирования классификатора. Для каждой этно-географической группы соотношение между группами для отбора наиболее информативных ДНК маркеров и для построения и тестирования классификатора было сделано в пропорции от 30% до 40%.
Все указанные образцы были генотипированы на панели Infinium Omni5Exome-4 v1.3 BeadChip (Illumina, США), в результате чего имелись генотипы по 4 559 465 полиморфизмам.
Для наиболее эффективного применения методов машинного обучения для отбора наиболее информативных ДНК маркеров данные были отфильтрованы с помощью программного пакета PLINK 1.9, а именно:
1. Отобраны аутосомные полиморфизмы и полиморфизмы X хромосомы (4 549 591 полиморфизмов из 4 559 465).
2. Отобраны полиморфизмы, прочитанные у более 95 процентов индивидов (4 198 721 полиморфизмов из 4 549 591).
3. Отобраны полиморфизмы, у которых встречались хотя бы два вхождения альтернативного аллеля среди всех индивидов (частота альтернативного аллеля более 0.0005) (4 140 504 полиморфизмов из 4 198 721).
4. Отбор индивидов, у которых прочитано (есть данные генотипов) более 90 процентов полиморфизмов.
5. Исключение полиморфизмов, сцепленных с другими полиморфизмами панели (критерий отбора - коэффициент детерминации R2 не превышает 0.2 для полиморфизмов, находящихся на расстоянии ближе 1 500 000 пар оснований к исследуемому) (817120 полиморфизмов из 4140504)
После применения данных фильтров для анализа остались 817 120 полиморфизмов.
Для корректного применения методов машинного обучения необходимо анализировать наиболее дифференцированные между собой популяции, объем которых не менее чем 20 образцов. В связи с этим популяции были объединены в более крупные этно-географические группы (далее ЭГГ). Выбор ЭГГ и их состав описаны выше (описание этапа 3).
На первом шаге отбора наиболее информативных ДНК маркеров были отобраны по 2 000 полиморфизмов на каждую ЭГГ, показавшие наиболее высокое значение в статистическом тесте по критерию Хи-квадрат между данным полиморфизмом и принадлежностью индивидов к данной ЭГГ. Таким образом выбираются полиморфизмы, имеющие наибольшую связь с принадлежностью индивидов к ЭГГ. На данном этапе после исключения повторов были отобраны 54 522 полиморфизма (часть полиморфизмов входила в наилучшие для нескольких ЭГГ).
На втором этапе отбора была построена модель классификатора на основе случайного леса решающих деревьев Extra-treesclassifier. Данная модель строит множество решающих деревьев, используя случайные наборы признаков (в нашем случае полиморфизмов) и образцов (в нашем случае индивидов), и усредняет результаты предсказаний этих деревьев. Это позволяет избежать переобучения модели (сильное ухудшение предсказаний на новых данных по сравнению с данными, на которых обучалась модель). Важной для нас особенностью данной модели является возможность получить численные характеристики важности отдельных полиморфизмов для итогового предсказания ЭГГ. Задавая границу значения важности, мы можем отобрать наиболее информативные для предсказания ЭГГ полиморфизмы. Для компенсации неравного количества индивидов в различных ЭГГ при обучении каждой ЭГГ был присвоен вес, обратно пропорциональный размеру ЭГГ.
Также для предотвращения переобучения применялся прием “кросс-валидация”. Набор данных был разделен на пять равных частей, в которые в равной степени входят представители каждой ЭГГ (метод Stratified K-Folds). Модель обучается на каждой из комбинаций из четырех частей, в результате чего получаются пять наборов полиморфизмов. Из этих пяти наборов берется пересечение, как итоговые наиболее важные для предсказания ЭГГ полиморфизмы.
Для подбора наилучших гиперпараметров модели (характеристик модели, задающихся до обучения, и не изменяющихся в процессе) для каждого набора гиперпараметров было произведено обучение описанной модели, отбор лучших полиморфизмов, обучение модели логистической регрессии на каждой из комбинаций четырех из пяти частей кросс-валидации, после чего на каждой из оставшихся неиспользованных частей кросс-валидации получена F1 метрика качества предсказания. Оптимальными гиперпараметрами считались те, при которых среднее значение метрики F1 было максимальным.
Метрика F1 представляет собой среднее гармоническое между точностью (precision) и полнотой (recall) предсказания и позволяет оценить, насколько правильны предсказания классификатора для каждого класса (в нашем случае ЭГГ).
После отбора наиболее информативных ДНК маркеров с оптимальными гиперпараметрами были отобраны 4 851 полиморфизм.
Для улучшения качества предсказания ЭГГ в итоговый список также были добавлены по 100 полиморфизмов для пар или троек ЭГГ, генетически близких между собой, и различие происхождения из которых вызывало наибольшие затруднения у алгоритма (например, мордва и украинцы, казахи и киргизы). Для этого для каждой пары или тройки ЭГГ был построен график главных компонент, включающий только индивидов этих ЭГГ, была найдена компонента, разделяющая данные ЭГГ, и были выбраны по 100 полиморфизмов, вносящие наибольший вклад в данную компоненту (наибольшие по модулю).
Объединенный список состоял из 5229 аутосомных ДНК маркеров, наиболее информативных для определения происхождения индивида (Фигура 1).
Наконец, финальной задачей было построение алгоритма определения вероятности происхождения из различных этно-географических групп с использованием выбранной панели маркеров.
Для этого был взят полный массив генотипированных образцов. Образцы были распределены по ЭГГ (этно-географическим группам) согласно их популяционный принадлежности. Часть образцов, относящихся к популяциям, не включенным ни в одну из ЭГГ, исключили из анализа. Полученный список из 1773 образцов - обучающая выборка.
Для каждой ЭГГ с помощью библиотеки pandas для python 3 рассчитали средние координаты и записали в таблицу. Генотипы по 5229 маркерам были переведены сначала в формат .vcf, а затем с помощью pandas для python 3 перевели из формата 0/0, 0/1, 1/1 в формат записи генотипов 0, 1, 2 соответственно. Было рассчитано среднее значение генотипа для каждого маркера по каждой ЭГГ. Усреднением по всем ЭГГ (с равным весом каждой ЭГГ) для каждого маркера, получили по одному значению для каждого маркера – эти значения сохранили и в дальнейшем использовали для импутации отсутствующих генотипов. Импутировали непрочитанные генотипы в извлеченных данных для обучающей выборки.
Для выбора модели разбили обучающую выборку на две — пробную и валидирующую. В валидирующую выборку вошло 30% образцов. Разбиение было выполнено с помощью метода sklearn.model_selection.train_test_split с параметрами test_size=0.3 (доля образцов во второй выборке, в нашем случае валидирующей), shuffle=True (случайное перемешивание образцов внутри каждой ЭГГ перед разделением), stratify=y (сохранение в обеих полученных выборках пропорций ЭГГ как в выборке до разбиения; в переменной y записана информация о принадлежности образцов к ЭГГ) и фиксированным генератором случайных чисел для воспроизводимости результата.
Далее перебирались модели логистической регрессии (для этого использовался класс sklearn.linear_model.LogisticRegression()). Для уменьшения эффекта переобучения применялась L2-регуляризация, уменьшающая значения коэффициентов при маркерах, но не зануляющая их, как L1-регуляризация (гиперпараметр penalty=’l2’). Изменялись следующие гиперпараметры:
Подход к количеству классов более двух (гиперпараметр multi_class). Были испробованы подходы ‘ovr’ (один против остальных) и ‘multinomial’ (мультиномиальный), первый из которых означает сравнение каждого класса со всеми оставшимися, то есть по сути обучение 29 бинарных предикторов на основе обычной логистической регрессии, а второй - применение многомерных формул для регрессии внутри одного предиктора.
Вес классов при обучении (гиперпараметр class_weight). Были испробованы подход ‘None’, дающий равный вес всем ЭГГ, и подход ‘balanced’, рассчитывающий веса классов обратно пропорционально их размеру, что сглаживает неравномерность в размере ЭГГ.
Сила регуляризации (гиперпараметр C - численный коэффициент, обратный силе регуляризации). Позволяет контролировать величину коэффициентов, избегая переобучения или уменьшая его последствия. Перебирались значения от 0.1 до 100.
Каждая модель с различным набором гиперпараметров обучалась на пробной выборке с помощью метода sklearn.linear_model.LogisticRegression().fit(), после чего предсказывала ЭГГ для образцов из валидирующей выборки с помощью метода sklearn.linear_model.LogisticRegression().predict(). По полученным предсказаниям рассчитывались метрики качества предсказания - precision, recall и f1-score (среднее гармоническое precision и recall) по каждой ЭГГ, а также их усредненные и взвешенные значения по всем классам. Также рассчитывалось значение accuracy. Сравнивалось среднее значение f1-score по всем классам. Наилучший результат показывали модели на основе подхода multinomial, однако проверка предсказания данных моделей на образцах, не вошедших ни в одну ЭГГ, показала, что данные модели предсказывают ошибочные ЭГГ со слишком большой вероятностью (>0.95), что является минусом для предикторов с нашей целью. Поэтому был выбран подход ‘ovr’, и он показал наилучшие значения при обучении с class_weight=’balanced’ и C=1.
Далее модель с выбранными параметрами была обучена на всем обучающем датасете. Полученные коэффициенты являются итоговыми и используются для предсказания в программе. Значения коэффициентов представлены в Таблице 1.
Примеры использования
Пример применения 1:
определение этно-географической группы населения и территории происхождения для образца №1.
Для образца №1 был известен генотип по 5229 маркерам, входящим в разработанную панель (Фиг.15). По формуле
, ,
i,
I
j
j
где  и  определены в таблице 1 для каждого полиморфизма где представляет собой целое число от 1 до 5229,и для каждой из этно-географических групп , где  представляет собой целое число от 1 до 29,
i
– генотип образца по полиморфизму , определенный на Фиг.15 в формате 0, 1 или 2, где 0 указывается для гомозиготного генотипа по референсному аллелю, 1 указывается для гетерозиготного генотипа, а 2 указывается для гомозиготного генотипа по альтернативному аллелю,
были рассчитаны значения  (вероятности происхождения) для каждой из 29 этно-географических групп. Полученные значения вероятностей представлены в Таблице 4.
Таблица 4. Значения вероятности происхождения из различных этно-географических групп, вычисленные для образца №1.
Наименования этно-географической группы (ЭГГ)
<img file="00000006.jpg" he="6" wi="4" img-content="undefined" img-format="jpg" />, вычисленное для образца №1
Широта
Долгота
Нанайцы, нивхи, орочи и ульчи
0
50,97
138,72
Башкиры
0,0025
54,15
56,45
Буряты и хамнигане
0
55,74
114,28
Чеченцы и ингуши
0,0003
43,33
44,95
Чукчи, коряки и ительмены
0
59,38
162,26
Дагестан
0,0002
41,98
47,39
Эвенки и эвены
0
56,69
141,16
Карелы и вепсы
0,0005
60,84
33,77
Казахи, каракалпаки, уйгуры и ногайцы
0,0026
45,14
57,48
Хакасы и южные алтайцы
0,0001
52,29
89,04
Ханты, манси и ненцы
0,0001
63,68
63,74
Коми и удмурты
0,0023
59,22
54,34
Киргизы
0,001
40,54
75,46
Марийцы и чуваши
0,0002
56,05
47,47
Монголы и калмыки
0
47,35
93,44
Мордва
0,0007
54,59
42,77
Осетины
0,002
42,81
44,15
Северные русские
0
59,87
40,25
Южные русские
0
54,94
35,97
Северные русские Архангельской области
0
65,19
43,74
Шорцы и северные алтайцы
0,0001
52,61
87,71
Сибирские татары
0
57,25
70,03
Таджики, памирцы и ягнобцы
0,0026
38,43
70,22
Татары
0,0019
53,49
50,28
Закавказье и Крым
0,0056
42,08
41,50
Тувинцы и тофалары
0
51,37
94,59
Украинцы
0,0016
49,83
29,42
Узбеки и туркмены
0,0002
41,06
64,72
Западный Кавказ
0,9983
43,90
41,53
Максимальным значением из содержащихся в таблице является 0,9983 для этно-географической группы «Западный Кавказ». Был сделан вывод, что предки индивидуума (донора образца №1) происходят из коренных народов Западного Кавказа.
Далее, для определения территории происхождения предков индивидуума полученные вероятности происхождения из различных этно-географических групп, были картографированы. Для этого по полученным 29 значениям (Таблица 4) была вычислена сеточная цифровая модель вероятности происхождения как функция географических координат по интерполяционной формуле:
, ,
G
R
j
где  – значение вероятности происхождения в точке  с географическими координатами (), () – географические координаты точек привязки этно-географических групп (представлены в Таблице 4 в столбцах «широта» и «долгота»),  – расстояние между точками  и  на сфере радиуса =6371 км (радиус Земли), суммирование велось по всем этно-географическим группам , точки привязки которых отстоят от точки G не более чем на 1000 км.
Значения цифровой модели были вычислены в узлах квадратной сетки в проекции карты, координаты узлов сетки перед подстановкой в интерполяционную формулу пересчитаны в долготу и широту по обратным формулам проекции,
территория, занятая на карте узлами цифровой модели со значениями  выше 0,2, была признана территорией происхождения индивида.
Для отображения этой территории цифровая модель была визуализирована на географической карте способом изолиний с послойной окраской, полученная карта представлена на Фиг. 16. Территории происхождения соответствует зона на карте, ограниченная изолинией со значением 0,2. Зона в пределах этой изолинии окрашена коричневыми и желтыми оттенками. На карте (Фиг. 16) можно видеть, что территория происхождения занимает Западный Кавказ, северное Причерноморье и Крым.
Данный образец ДНК был получен от индивидуума, рожденного в г. Карачаевск (Республика Карачаево-Черкессия). Его мать и отец, как и их родители, были карачаевцами, часть из них родилась в селе Усть-Джегута, часть родилась в г. Карачаевск.
Таким образом, определенная по предложенному способу этно-географическая группа (Западный Кавказ) была определена верно, так как карачаевцы являются одним из коренных народов Западного Кавказа и территория происхождения, определенная по предложенному способу, также была определена верно, т.к. Карачаевск и Усть-Джегута входят в территорию, отмеченную на Фиг. 16.
Таким образом, пример иллюстрирует, что разработанная система позволяет правильно определять этно-географическую группу и территорию происхождения для уроженцев Кавказа.
Пример применения 2:
определение этно-географической группы населения и территории происхождения для образца №2.
Для образца №2 был известен генотип по 5229 маркерам, входящим в разработанную панель (Фиг.15). По формуле
, ,
i,
i
j
j
где  и  определены в таблице 1 для каждого полиморфизма где представляет собой целое число от 1 до 5229, и для каждой из этно-географических групп , где  представляет собой целое число от 1 до 29,
i
– генотип образца по полиморфизму , определенный на Фиг.15, в формате 0, 1 или 2, где 0 указывается для гомозиготного генотипа по референсному аллелю, 1 указывается для гетерозиготного генотипа, а 2 указывается для гомозиготного генотипа по альтернативному аллелю,
были рассчитаны значения  (вероятности происхождения) для каждой из 29 этно-географических групп. Полученные значения вероятностей представлены в Таблице 5.
Таблица 5. Значения вероятности происхождения из различных этно-географических групп, вычисленные для образца №2.
Наименования этно-географической группы (ЭГГ)
<img file="00000006.jpg" he="6" wi="4" img-content="undefined" img-format="jpg" />, вычисленное для образца №2
Широта
Долгота
Нанайцы, нивхи, орочи и ульчи
0,0000
50,97
138,72
Башкиры
0,0001
54,15
56,45
Буряты и хамнигане
0,0000
55,74
114,28
Чеченцы и ингуши
0,0001
43,33
44,95
Чукчи, коряки и ительмены
0,0000
59,38
162,26
Дагестан
0,0002
41,98
47,39
Эвенки и эвены
0,0000
56,69
141,16
Карелы и вепсы
0,0000
60,84
33,77
Казахи, каракалпаки, уйгуры и ногайцы
0,0000
45,14
57,48
Хакасы и южные алтайцы
0,0001
52,29
89,04
Ханты, манси и ненцы
0,0000
63,68
63,74
Коми и удмурты
0,0016
59,22
54,34
Киргизы
0,0000
40,54
75,46
Марийцы и чуваши
0,0000
56,05
47,47
Монголы и калмыки
0,0000
47,35
93,44
Мордва
0,0000
54,59
42,77
Осетины
0,0001
42,81
44,15
Северные русские
0,0048
59,87
40,25
Южные русские
0,9976
54,94
35,97
Северные русские Архангельской области
0,0006
65,19
43,74
Шорцы и северные алтайцы
0,0000
52,61
87,71
Сибирские татары
0,0009
57,25
70,03
Таджики, памирцы и ягнобцы
0,0003
38,43
70,22
Татары
0,0010
53,49
50,28
Закавказье и Крым
0,0000
42,08
41,50
Тувинцы и тофалары
0,0000
51,37
94,59
Украинцы
0,0002
49,83
29,42
Узбеки и туркмены
0,0001
41,06
64,72
Западный Кавказ
0,0000
43,90
41,53
Максимальным значением из содержащихся в таблице является 0,9975 для этно-географической группы «Южные русские». Был сделан вывод, что предки индивидуума (донора образца №2) происходят из южных русских популяций.
Далее, для определения территории происхождения предков индивидуума полученные вероятности происхождения из различных этно-географических групп, были картографированы. Для этого по полученным 29 значениям (Таблица 5) была вычислена сеточная цифровая модель вероятности происхождения как функция географических координат по интерполяционной формуле:
,  ,
G
R
j
где  – значение вероятности происхождения в точке  с географическими координатами (), () – географические координаты точек привязки этно-географических групп (представлены в Таблице 5 в столбцах «широта» и «долгота»),  – расстояние между точками  и  на сфере радиуса =6371 км (радиус Земли), суммирование велось по всем этно-географическим группам , точки привязки которых отстоят от точки G не более чем на 1000 км.
Значения цифровой модели были вычислены в узлах квадратной сетки в проекции карты, координаты узлов сетки перед подстановкой в интерполяционную формулу пересчитаны в долготу и широту по обратным формулам проекции,
территория, занятая на карте узлами цифровой модели со значениями  выше 0,2, была признана территорией происхождения индивида.
Для отображения этой территории цифровая модель была визуализирована на географической карте способом изолиний с послойной окраской, полученная карта представлена на Фиг. 17. Территории происхождения соответствует зона на карте, ограниченная изолинией со значением 0,2. Зона в пределах этой изолинии окрашена коричневыми и желтыми оттенками. На карте (Фиг. 17) можно видеть, что территория происхождения занимает южную часть Центральной России, а также прилегающие территории Украины, Белоруссии и Литвы.
Данный образец ДНК был получен от индивидуума, рожденного в д. Круглое Красненского района Белгородской области. Его мать и отец, как и их родители, были русскими, все они родились в Красненском районе Белгородской области.
Таким образом, определенная по предложенному способу этно-географическая группа (южные русские) была определена верно, так как русское население Белгородской области относится к южным русским, и территория происхождения, определенная по предложенному способу, также была определена верно, т.к. Белгородская области входит в территорию, отмеченную на Фиг. 17.
Таким образом, пример иллюстрирует, что разработанная система позволяет правильно определять этно-географическую группу и территорию происхождения для этнических русских.
Пример применения 3:
определение этно-географической группы населения и территории происхождения для образца №3.
Для образца №3 был известен генотип по 5229 маркерам, входящим в разработанную панель (Фиг.15). По формуле
, ,
i,
i
j
j
где  и  определены в таблице 1 для каждого полиморфизма где представляет собой целое число от 1 до 5229,и для каждой из этно-географических групп , где  представляет собой целое число от 1 до 29,
i
– генотип образца по полиморфизму , определенный на Фиг.15, в формате 0, 1 или 2, где 0 указывается для гомозиготного генотипа по референсному аллелю, 1 указывается для гетерозиготного генотипа, а 2 указывается для гомозиготного генотипа по альтернативному аллелю,
были рассчитаны значения  (вероятности происхождения) для каждой из 29 этно-географических групп. Полученные значения вероятностей представлены в Таблице 6.
Таблица 6. Значения вероятности происхождения из различных этно-географических групп, вычисленные для образца №3.
Наименования этно-географической группы (ЭГГ)
<img file="00000006.jpg" he="6" wi="4" img-content="undefined" img-format="jpg" />, вычисленное для образца №3
Широта
Долгота
Нанайцы, нивхи, орочи и ульчи
0,0015
50,97
138,72
Башкиры
0,0004
54,15
56,45
Буряты и хамнигане
0,9995
55,74
114,28
Чеченцы и ингуши
0,0000
43,33
44,95
Чукчи, коряки и ительмены
0,0004
59,38
162,26
Дагестан
0,0000
41,98
47,39
Эвенки и эвены
0,0001
56,69
141,16
Карелы и вепсы
0,0000
60,84
33,77
Казахи, каракалпаки, уйгуры и ногайцы
0,0026
45,14
57,48
Хакасы и южные алтайцы
0,0031
52,29
89,04
Ханты, манси и ненцы
0,0002
63,68
63,74
Коми и удмурты
0,0000
59,22
54,34
Киргизы
0,0038
40,54
75,46
Марийцы и чуваши
0,0004
56,05
47,47
Монголы и калмыки
0,0023
47,35
93,44
Мордва
0,0000
54,59
42,77
Осетины
0,0001
42,81
44,15
Северные русские
0,0000
59,87
40,25
Южные русские
0,0000
54,94
35,97
Северные русские Архангельской области
0,0000
65,19
43,74
Шорцы и северные алтайцы
0,0002
52,61
87,71
Сибирские татары
0,0001
57,25
70,03
Таджики, памирцы и ягнобцы
0,0000
38,43
70,22
Татары
0,0003
53,49
50,28
Закавказье и Крым
0,0000
42,08
41,50
Тувинцы и тофалары
0,0018
51,37
94,59
Украинцы
0,0000
49,83
29,42
Узбеки и туркмены
0,0013
41,06
64,72
Западный Кавказ
0,0003
43,90
41,53
Максимальным значением из содержащихся в таблице является 0,9995 для этно-географической группы «Буряты и хамнигане». Был сделан вывод, что предки индивидуума (донора образца №3) происходят из популяций бурят или хамниган.
Далее, для определения территории происхождения предков индивидуума полученные вероятности происхождения из различных этно-географических групп, были картографированы. Для этого по полученным 29 значениям (Таблица 6) была вычислена сеточная цифровая модель вероятности происхождения как функция географических координат по интерполяционной формуле:
,  ,
G
R
j
где  – значение вероятности происхождения в точке  с географическими координатами (), () – географические координаты точек привязки этно-географических групп (представлены в Таблице 6 в столбцах «широта» и «долгота»),  – расстояние между точками  и  на сфере радиуса =6371 км (радиус Земли), суммирование велось по всем этно-географическим группам , точки привязки которых отстоят от точки G не более чем на 1000 км.
Значения цифровой модели были вычислены в узлах квадратной сетки в проекции карты, координаты узлов сетки перед подстановкой в интерполяционную формулу пересчитаны в долготу и широту по обратным формулам проекции,
территория, занятая на карте узлами цифровой модели со значениями  выше 0,2, была признана территорией происхождения индивида.
Для отображения этой территории цифровая модель была визуализирована на географической карте способом изолиний с послойной окраской, полученная карта представлена на Фиг. 18. Территории происхождения соответствует зона на карте, ограниченная изолинией со значением 0,2. Зона в пределах этой изолинии окрашена коричневыми и желтыми оттенками. На карте (Фиг. 18) можно видеть, что территория происхождения занимает территорию к востоку и к северу от озера Байкал.
Данный образец ДНК был получен от индивидуума, рожденного в п.г.т. Могойтуй Агинского Бурятского округа Забайкальского края. Его мать и отец, как и их родители, были бурятами, все они родились в Забайкальском крае (Ононский, Могонтуйский, Борзенский районы).
Таким образом, определенная по предложенному способу этно-географическая группа (буряты и хамнигане) была определена верно, так как предки донора образца №3 являлись бурятами, и территория происхождения, определенная по предложенному способу, также была определена верно, т.к. Забайкальский край входит в территорию, отмеченную на Фиг. 18.
Таким образом, пример иллюстрирует, что разработанная система позволяет правильно определять этно-географическую группу и территорию происхождения для коренных народов Сибири.
Список литературы
1. База данных COSMIC (www.cancer.sanger.ac.uk/cosmic)
2. База данных Database of Genomic Variants (www.dgv.tcag.ca/dgv/app/home)
3. База данных ExSNP (www.exsnp.org)
4. База данных Genome Browser (www.genome.ucsc.edu)
5. База данных GOCVS (www.geneontology.org/page/go-cvs-repository)
6. База данных NCBI Genome Reference Consortium GRCh37 (www.ncbi.nlm.nih.gov/grc/human)
7. База данных Neandertal Genome Browser (www.projects.ensembl.org/neandertal)
8. База данных Pharma ADME (www.pharmaadme.org)
9. База данных RefSeq (www.ncbi.nlm.nih.gov/refseq)
10. Балановская Е.В., Жабагин М.К., Агджоян А.Т., Чухряева М.И., Маркина Н.В. с соавт. Популяционные биобанки: принципы организации и перспективы применения в геногеографии и персонализированной медицине // Генетика, 2016. Т. 52. № 12. С. 1371–1387. DOI: 10.7868/S001667581612002X.
11. Каталог NHGRI GWAS Catalog (www.ebi.ac.uk/gwas).
12. Песик В. Ю., Федюнин А. А., Агджоян А. Т., Чухряева М. И., Утевская О. М., Евсеева И. В., Чурносов М. И., Лепендина И., Игнашкин М. А., Богунов Ю. В., Балановская Е. В., Орехов В. А., Балановский О. П. Разнообразие региональных русских популяций по STR  маркерам, используемым при ДНК- идентификации // Генетика. - 2014. T. 50, № 6. - C. 715-723.
13. Юсупов Ю. М., Балановская Е. В., Сабитов Ж. М., Балановский О. П. Комплексные исследования этногенеза: союз геногеографии и этнологии // Вестник антропологии. Ч. - 2017. - T. 38, № 2. - C. 28-35.AlexanderD. H., Lange K. Enhancements to the ADMIXTURE algorithm for individual ancestry estimation // BMC Bioinformatics. - 2011. - T. 12. - C. 246.
14. Alhusain L, Hafez AM. Nonparametric approaches for population structure analysis. Human Genomics. 2018;12:25. doi:10.1186/s40246-018-0156-4.
15. Barbujani G., Bertorelle G. Genetics and the population history of Europe // Proc Natl Acad Sci U S A. - 2001. - T. 98, № 1. - C. 22-5.
16. Chang CC, Chow CC, Tellier LC, et al. Second-generation PLINK: rising to the challenge of larger and richer datasets. Gigascience. 2015;4:7. doi:10.1186/s13742-015-0047-8.
17. Chikhi L., Nichols R. A., Barbujani G., Beaumont M. A. Y genetic data support the Neolithic demic diffusion model // Proc Natl Acad Sci U S A. - 2002. - T. 99, № 17. - C. 11008-13.
18. Elhaik E., Tatarinova T., Chebotarev D., Piras I. S., Calo C. M., De Montis A., Atzori M., Marini M., Tofanelli S., Francalacci P., Pagani L., Tyler-Smith C., Xue Y. L., Cucca F., Schurr T. G., Gaieski J. B., Melendez C., Vilar M. G., Owings A. C., Gomez R., Fujita R., Santos F. R., Comas D., Balanovsky O., Balanovska E., Zalloua P., Soodyall H., Pitchappan R., GaneshPrasad A., Hammer M., Matisoo-Smith L., Wells R. S., Royyuru A. K., Consortium G. Geographic population structure analysis of worldwide human populations infers their biogeographical origins (vol 5, 3513, 2014) // Nature Communications. – 2016. - T. 7.
19. Hartigan JA, Wong MA. Algorithm AS 136: A K-Means Clustering Algorithm. Journal of the Royal Statistical Society. 1979;28:100-108. https://ru.scribd.com/document/293724711/A-K-Means-Clustering-Algorithm.
20. Jeong C, Balanovsky O, Lukianova E, et al. The genetic history of admixture across inner Eurasia. Nature Ecology & Evolution. 2019;3:966–976. doi:10.1038/s41559-019-0878-2.
21. Karmin M., Saag L., Vicente M., Wilson Sayres M. A., Järve M., Talas U. G., Rootsi S., Ilumäe A. M., Mägi R., Mitt M., Pagani L., Puurand T., Faltyskova Z., Clemente F., Cardona A., Metspalu E., Sahakyan H., Yunusbayev B., Hudjashov G., DeGiorgio M., Loogväli E. L., Eichstaedt C., Eelmets M., Chaubey G., Tambets K., Litvinov S., Mormina M., Xue Y., Ayub Q., Zoraqi G., Korneliussen T. S., Akhatova F., Lachance J., Tishkoff S., Momynaliev K., Ricaut F. X., Kusuma P., Razafindrazaka H., Pierron D., Cox M. P., Sultana G. N., Willerslev R., Muller C., Westaway M., Lambert D., Skaro V., Kovačevic L., Turdikulova S., Dalimova D., Khusainova R., Trofimova N., Akhmetova V., Khidiyatova I., Lichman D. V., Isakova J., Pocheshkhova E., Sabitov Z., Barashkov N. A., Nymadawa P., Mihailov E., Seng J. W., Evseeva I., Migliano A. B., Abdullah S., Andriadze G., Primorac D., Atramentova L., Utevska O., Yepiskoposyan L., Marjanovic D., Kushniarevich A., Behar D. M., Gilissen C., Vissers L., Veltman J. A., Balanovska E., Derenko M., Malyarchuk B., Metspalu A., Fedorova S., Eriksson A., Manica A., Mendez F. L., Karafet T. M., Veeramah K. R., Bradman N., Hammer M. F., Osipova L. P., Balanovsky O., Khusnutdinova E. K., Johnsen K., Remm M., Thomas M. G., Tyler-Smith C., Underhill P. A., Willerslev E., Nielsen R., Metspalu M., Villems R., Kivisild T. A recent bottleneck of Y chromosome diversity coincides with a global change in culture // Genome Res. – 2015. - T. 25, № 4. - C. 459-66.
22. Koshel SM. Geoinformation technologies in genogeography. Modern geographic cartography. 2012;158-166. (In Russ.) Кошель С.М. Геоинформационные технологии в геногеографии. Современная географическая картография. 2012; 158-166. https://www.researchgate.net/publication/294848419_Geoinformacionnye_tehnologii_v_genogeografii
23. Kosoy R, Nassir R, Tian C, et al. (2009) Ancestry informative marker sets for determining continental origin and admixture proportions in common populations in America. HumMutat 30(1) 69–78.
24. Kozlov K., Chebotarev D., Hassan M., Triska M., Triska P., Flegontov P., Tatarinova T. V. Differential Evolution approach to detect recent admixture // Bmc Genomics. – 2015. - T. 16.
25. Lee C, Abdool A, Huang C. PCA-based population structure inference with generic clustering algorithms. BMC Bioinformatics. 2009;10:73. doi:10.1186/1471-2105-10-S1-S73.
26. Liu N, Zhao H. A non-parametric approach to population structure inference using multilocus genotypes. Hum Genomics. 2006;2: 353. doi:10.1186/1479-7364-2-6-353.
27. Manichaikul A, Mychaleckyj JC, Rich SS, et al. Robust relationship inference in genome-wide association studies. Bioinformatics. 2010;26:2867-2873. doi:10.1093/bioinformatics/btq559
28. Novembre J., Johnson T., Bryc K., Kutalik Z., Boyko A. R., Auton A., Indap A., King K. S., Bergmann S., Nelson M. R., Stephens M., Bustamante C. D. Genes mirror geography within Europe // Nature. – 2008. - T. 456, № 7218. - C. 98-101.
29. Patterson N, Price AL, Reich D. Population Structure and Eigenanalysis. PLoS Genet. 2006;2:190. doi:10.1371/journal.pgen.0020190.
30. Pedregosa F, Varoquaux G, Gramfort A, et al. Scikit-learn: Machine Learning in Python. Journal of Machine Learning Research. 2011;12:2825-2830. https://www.researchgate.net/publication/51969319_Scikit-learn_Machine_Learning_in_Python
31. Rosser Z. H., Zerjal T., Hurles M. E., Adojaan M., Alavantic D., Amorim A., Amos W., Armenteros M., Arroyo E., Barbujani G., Beckman G., Beckman L., Bertranpetit J., Bosch E., Bradley D. G., Brede G., Cooper G., Corte-Real H., de Knijff P., Decorte R., Dubrova Y. E., Evgrafov O., Gilissen A., Glisic S., Golge M., Hill E. W., Jeziorowska A., Kalaydjieva L., Kayser M., Kivisild T., Kravchenko S. A., Krumina A., Kucinskas V., Lavinha J., Livshits L. A., Malaspina P., Maria S., McElreavey K., Meitinger T. A., Mikelsaar A. V., Mitchell R. J., Nafa K., Nicholson J., Norby S., Pandya A., Parik J., Patsalis P. C., Pereira L., Peterlin B., Pielberg G., Prata M. L., Previdere C., Roewer L., Rootsi S., Rubinsztein D. C., Saillard J., Santos F. R., Stefanescu G., Sykes B. C., Tolun A., Villems R., Tyler-Smith C., Jobling M. A. Y-chromosomal diversity in Europe is clinal and influenced primarily by geography, rather than by language // American Journal of Human Genetics. – 2000. - T. 67, № 6. - C. 1526-1543.
32. Semino O., Passarino G., Oefner P. J., Lin A. A., Arbuzova S., Beckman L. E., De Benedictis G., Francalacci P., Kouvatsi A., Limborska S., Marcikiae M., Mika A., Mika B., Primorac D., Santachiara-Benerecetti A. S., Cavalli-Sforza L. L., Underhill P. A. The genetic legacy of Paleolithic Homo sapiens sapiens in extant Europeans: a Y chromosome perspective // Science. – 2000. - T. 290, № 5494. - C. 1155-9.
33. Triska P, Chekanov N, Stepanov V, et al. Between Lake Baikal and the Baltic Sea: genomic history of the gateway to Europe. BMC Genet. 2017; Dec 28;18(Suppl 1):110. doi:10.1186/s12863-017-0578-3.