Sampling
Piśmienną świadomą zgodę uzyskano od osobnika Saami, którego genom analizowano w tym badaniu, które zostało zatwierdzone przez komisję etyczną szpitala w Helsinkach i Uusimaa (decyzja 329/13/03/00/2013) oraz komisję etyczną Uniwersytetu w Lipsku (numer referencyjny zatwierdzenia 398-13-16122013).
Pobieranie i ekstrakcja starożytnego DNA wymaga ścisłej procedury w celu uniknięcia kontaminacji wprowadzonej przez współczesny materiał genetyczny. W przypadku 13 osobników z epoki żelaza z Finlandii, które mieliśmy do dyspozycji, pobieranie próbek odbywało się w sterylnym pomieszczeniu przeznaczonym do pracy ze starożytnym DNA, w Instytucie Nauk Archeologicznych w Tybindze. Wstępny proces pracy obejmował dokumentowanie, fotografowanie i przechowywanie próbek w indywidualnych, oznaczonych kodami identyfikacyjnymi plastikowych probówkach i plastikowych torebkach. W wyniku wczesnych badań pilotażowych, próbki zębów, których używaliśmy, były rozdrobnione, a część zębiny została usunięta. Pozostała zębina została pobrana poprzez ostrożne oddzielenie jej od szkliwa za pomocą wiertła dentystycznego i schłodzonych diamentowych głowic wiertarskich, obracanych z prędkością poniżej 15 obr/min, w celu uniknięcia możliwego uszkodzenia starożytnego DNA spowodowanego wysoką temperaturą.
Dla próbek z miejsc Bolshoy i Chalmny Varre, użyliśmy resztek proszku zębowego, który został pierwotnie przetworzony w Instytucie Antropologii na Uniwersytecie w Moguncji dla celów replikacji, jak opisano w ref. 24. W skrócie, etapy przygotowania próbki obejmowały naświetlanie promieniami UV przez 30-45 min, a następnie delikatne przetarcie powierzchni rozcieńczonymi komercyjnymi wybielaczami. Zęby zostały następnie wypiaskowane przy użyciu ścierniwa z tlenku glinu (Harnisch & Rieth) i zmielone na drobny proszek przy użyciu młynka z mieszadłem (Retsch).
Kalibracja daty radiowęglowej
Wykalibrowaliśmy datę radiowęglową Bolszoj, podaną w ref. 24,55 jako 3473 ± 42 lata BP, używając Intcal1356 jako krzywej kalibracyjnej, stosując OxCal 4.357.
Ekstrakcja DNA i przygotowanie biblioteki
DNA z sześciu próbek z Bolszoj i dwóch z Chalmny Varre zostało wyekstrahowane w obiektach starożytnego DNA Instytutu Maxa Plancka dla Nauki o Historii Człowieka (MPI-SHH) w Jenie, Niemcy. Ekstrakcja dla próbek z Levänluhta była podobnie przeprowadzona w pomieszczeniach czystych Instytutu Nauk Archeologicznych w Tybindze. Dla każdej próbki, ~ 50 mg sproszkowanej zębiny zostało użyte do procedury ekstrakcji specjalnie zaprojektowanej do pobierania starożytnego DNA58. Bufor ekstrakcyjny zawierający 0.45 M EDTA, pH 8.0 (Life Technologies) i 0.25 mg/ml Proteinazy K (Sigma-Aldrich) został dodany do proszku kostnego i inkubowany w temperaturze 37 °C z rotacją przez noc. Supernatant oddzielono od granulatu proszku kostnego przez wirowanie (14 000 obr./min). Do supernatantu dodano bufor wiążący składający się z 5 M GuHCL (Sigma Aldrich) i 40% izopropanolu (Merck), wraz z 400 μl 1 M octanu sodu (pH 5,5), a roztwór oczyszczono, wirując go przez kolumnę oczyszczającą podłączoną do lejka montażowego High Pure Extender (8 min. w 1500 obr./min, z powolnym przyspieszaniem). Następnie kolumnę odwirowano w probówce zbierającej (1 min 14 000 obr./min) 1-2 razy, aby zmaksymalizować wydajność. Następnie przeprowadzono dwa kolejne etapy płukania z użyciem 450 μl buforu płuczącego (High Pure Viral Nucleic Acid Large Volume Kit) oraz dwa etapy wirowania na sucho przez 1 min przy 14 000 obr. Ostateczną całkowitą objętość 100 μl eluatu uzyskano w dwóch oddzielnych rundach elucji 50 μl TET (10 mM Tris-HCL, 1 mM EDTA pH 8,0, 0,1% Tween20), z których każda wirowała przez 1 min przy 14 000 obr/min w świeżej probówce Eppendorf o pojemności 1,5 ml. Kontrole negatywne (bufor zamiast próbki) były przetwarzane równolegle w stosunku 1 kontrola na 7 próbek.
Z ekstraktu 100 μl, 20 μl zostało użyte do immortalizacji DNA próbki jako dwuniciowa biblioteka. Procedura obejmowała naprawę tępego końca, ligację adaptera i etapy wypełniania adaptera, jak opisali Meyer i Kircher59. Podczas etapu naprawy tępego końca, mieszanina 0,4 U/μl T4 PNK (kinaza polinukleotydowa) i 0,024 U/μl T4 polimerazy DNA, 1× NEB buffer 2 (NEB), 100 μM dNTP mix (Thermo Scientific), 1 mM ATP (NEB) i 0.Do szablonowego DNA dodawano 8 mg/ml BSA (NEB), a następnie inkubowano w termocyklerze (15 min 15 °C, 15 min 25 °C) i oczyszczano za pomocą zestawu MinElute (QIAGEN). Produkt był eluowany w 18 μl buforu TET. Etap ligacji adaptera obejmował mieszaninę 1× Quick Ligase Buffer (NEB), 250 nM adapterów Illumina (Sigma-Aldrich) i 0,125 U/μl Quick Ligase (NEB), dodaną do 18 μl eluatu, po czym następowała 20-minutowa inkubacja i drugi etap oczyszczania za pomocą kolumn MinElute, tym razem w 20 μl eluatu. Do etapu wypełniania dodawano mieszaninę 0,4 U/μl polimerazy Bst i 125 μM mieszaniny dNTP, a następnie mieszaninę inkubowano w termocyklerze (30 min 37 °C, 10 min 80 °C). Dla wszystkich 13 ekstraktów z Levänluhta wytworzono biblioteki bez traktowania Uracil-DNA-glikozylazą (UDG). Dodatkowo, biblioteki poddane obróbce UDG w połowie zostały wytworzone dla siedmiu z oryginalnych 13 ekstraktów z Levänluhta, oraz dla wszystkich ekstraktów z Bolshoy i Chalmny Varre. Aby wprowadzić obróbkę UDG, do przygotowania biblioteki włączono etap wstępny, w którym do 20 μl ekstraktu dodawano 250 U enzymu USER (NEB), po czym następowała inkubacja w 37 °C przez 30 min, a następnie w 12 °C przez 1 min. Po tym ponownie dodawano 200 UGI (inhibitor glikozylazy Uracylu, NEB) i następowała kolejna identyczna inkubacja w celu zatrzymania enzymatycznego wydzielania miejsc deaminowanych, jak opisano w60. Dla każdej biblioteki włączono unikalną parę indeksów o długości ośmiu pb, stosując polimerazę Pfu Turbo Cx Hotstart DNA i program termocyklingu o następującym profilu temperaturowym: wstępna denaturacja (98 °C przez 30 s), cykl denaturacji/ wyżarzania/ wydłużania (98 °C przez 10 s/ 60 °C przez 20 s/ 72 °C przez 20 s) i końcowe wydłużanie w 72 °C przez 10 min61. Proszek kostny z niedźwiedzia jaskiniowego był przetwarzany równolegle, służąc jako kontrola pozytywna. Negatywne kontrole dla obu etapów ekstrakcji i przygotowania biblioteki były przechowywane obok próbek przez cały czas trwania pracy.
Wydajność eksperymentu zapewniono przez ilościowe określenie stężenia bibliotek na qPCR (Roche) z wykorzystaniem podwielokrotnień z bibliotek przed i po indeksowaniu. Liczba kopii molekularnych w preindeksowanych bibliotekach wahała się od ~ 10E8 do ~ 10E9 kopii/μl, wskazując na udane przygotowanie biblioteki, podczas gdy indeksowane biblioteki wahały się od ~ 10E10 do ~ 10E12 kopii/μl, stwierdzając dopuszczalną wydajność indeksowania. Kontrole negatywne wykazywały 4-5 rzędów wielkości niższe stężenie niż próbki, wskazując na niski poziom zanieczyszczeń z etapów obróbki laboratoryjnej.
Biblioteki amplifikowano metodą PCR, dla ilości cykli odpowiadającej stężeniom indeksowanych bibliotek, przy użyciu polimerazy AccuPrime Pfx (5 μl szablonu biblioteki, 2 U polimerazy DNA AccuPrime Pfx firmy Invitrogen, 1 U gotowego 10× PCR mastermix i 0.3 μM primerów IS5 i IS6, na każde 100 μl reakcji) z profilem termicznym 2 min denaturacji w 95 °C, 3-9 cykli składających się z 15 s denaturacji w 95 °C, 30 s annealingu w 60 °C, 2 min elongacji w 68 °C i 5 min elongacji w 68 °C. Wzmocnione biblioteki były oczyszczane przy użyciu kolumn MinElute spin z zastosowaniem standardowego protokołu dostarczonego przez producenta (Qiagen), i kwantyfikowane do sekwencjonowania przy użyciu Agilent 2100 Bioanalyzer DNA 1000 chip.
Dla współczesnego osobnika Saami, całkowite DNA było ekstrahowane fenolowo-chloroformowo i fizycznie ścinane przy użyciu fragmentacji COVARIS. Zmodyfikowana biblioteka Illumina została przygotowana przy użyciu naprawy tępych końców, a następnie A-tailing 3′-końca i ligacji rozwidlonych adapterów. Po PCR indeksującym następowało wycięcie fragmentów o długości od 500 do 600 bp z 2% żelu agarozowego.
Sekwencjonowanie metodą capture &
Użyliśmy procedury in-solution capture z ref. 62 do wzbogacenia naszych bibliotek o fragmenty DNA pokrywające się z 1,237,207 zmiennymi pozycjami w ludzkim genomie4. Sekwencje użyte jako przynęta, dołączone do kulek magnetycznych, zostały zmieszane z szablonem próbki DNA w roztworze i pozostawione do hybrydyzacji z docelowym DNA podczas 24-godzinnej inkubacji w temperaturze 60 °C w piecu rotacyjnym. 4-6 próbek było łączonych w równych proporcjach masowych w łącznej ilości 2000 ng DNA. Pobrane biblioteki były sekwencjonowane (75 bp single-end, plus dodatkowe sparowane-end dla trzech nie-UDG bibliotek osobników z Levänluhta) na platformie Illumina HiSeq 4000 w Max Planck Institute for the Science of Human History w Jenie. Spośród 13 oryginalnie przetworzonych próbek z epoki żelaza z Finlandii, siedem okazało się być odpowiedniej jakości do wykorzystania w dalszych analizach. Współczesny genom Saamów został zsekwencjonowany w Genome Analyser II (8 pasów, 125 bp paired-end) w Max Planck Institute for Evolutionary Anthropology w Lipsku.
Przetwarzanie zsekwencjonowanych odczytów
Użyliśmy EAGER63 (wersja 1.92.50) do przetwarzania sekwencjonowanych odczytów, stosując domyślne parametry (patrz poniżej) dla danych sekwencjonowania pojedynczego końca pochodzących od człowieka, poddanych obróbce UDG-połowy, podczas przetwarzania bibliotek UDG-połowy dla wszystkich osobników. W szczególności, AdapterRemoval został użyty do przycięcia adapterów sekwencjonowania z naszych odczytów, z minimalnym nakładaniem się 1 bp, przy minimalnej jakości baz 20 i minimalnej długości sekwencji 30 bp. BWA aln (wersja 0.7.12-r1039, https://sourceforge.net/projects/bio-bwa/files)64 wykorzystano do mapowania odczytów do ludzkiej sekwencji referencyjnej hs37d5, z długością ziarna (-l) równą 32, maksymalną liczbą różnic (-n) równą 0,01, bez filtrowania jakości. Usuwanie duplikatów przeprowadzono przy użyciu programu DeDup v0.12.1. Obliczenia uszkodzeń deaminacji zasad końcowych wykonano przy użyciu mapDamage, określając długość (-l) 100 bp (Tabela 1).
Do analiz downstreamowych użyto bamutils (wersja 1.0.13, https://github.com/statgen/bamUtil.git) TrimBam do obcięcia dwóch zasad na początku i końcu wszystkich odczytów. Procedura ta eliminuje pozycje, które są dotknięte deaminacją, usuwając w ten sposób błędy genotypowania, które mogłyby powstać z powodu starożytnych uszkodzeń DNA.
Dla trzech osobników z Levänluhta przekraczających progowe pokrycie 1% we wstępnym przesiewaniu, użyliśmy bibliotek nie poddanych obróbce UUDG, aby potwierdzić autentyczność starożytnych danych. Dla tych nieprzetworzonych bibliotek przeprowadzono dwie rundy sekwencjonowania, które przetworzono przy użyciu programu EAGER z powyższymi parametrami, ale z określeniem trybu obróbki non-UDG i ustawieniem właściwego typu sekwencjonowania między bibliotekami. Połączone odczyty wyodrębniono z wynikowych plików bam i połączono z plikiem bam zawierającym odczyty z przebiegu sekwencjonowania pojedynczego końca za pomocą programu samtools merge (wersja 1.3)65.
Nowoczesny genom Saami wygenerowano przy użyciu programu Ibis do wywoływania zasad i wewnętrznego skryptu do przycinania adapterów. Uzyskane odczyty zostały następnie wyrównane do ludzkiego genomu referencyjnego hs37d5 przy użyciu bwa 0.5.9-r16 (parametry -e 20 -o 2 -n 0.01).
Genotypowanie
Do genotypowania 15 starożytnych osobników użyliśmy niestandardowego programu (pileupCaller). Plik pileup został wygenerowany przy użyciu samtools mpileup z parametrami -q 30 -Q 30 -B zawierającymi tylko miejsca pokrywające się z naszym panelem przechwytującym. Z tego pliku, dla każdego osobnika i każdego SNP w panelu 1240K, losowano jeden odczyt pokrywający SNP i dokonywano pseudo-haploidalnego wywołania, tzn. zakładano, że starożytny osobnik jest homozygotyczny dla allelu w losowo wylosowanym odczycie dla danego SNP. PileupCaller jest dostępny na stronie https://github.com/stschiff/sequenceTools.git.
Dla trzech bibliotek z Levänluhta, które nie zostały poddane obróbce UDG, genotypowaliśmy tylko transwersje, eliminując w ten sposób artefakty pośmiertnych uszkodzeń C- > T z dalszych analiz.
Strzelany genom współczesnego osobnika Saami został poddany genotypowaniu przy użyciu GATK (wersja 1.3-25-g32cdef9) Unified Genotyper po uprzednim wyrównaniu indeli. Wywołania wariantów były filtrowane dla wariantów z wynikiem jakości powyżej 30, a skrypt został użyty do konwersji wariantów do formatu EigenStrat.
Dane zostały połączone z dużym zbiorem danych składającym się z 3871 starożytnych i współczesnych osobników genotypowanych na tablicach Human Origins i/lub 1240K SNP, przy użyciu mergeit.
Określanie płci
Aby określić płeć genetyczną każdego starożytnego osobnika, obliczyliśmy pokrycie na autosomach, jak również na każdym chromosomie płci. Użyliśmy niestandardowego skryptu (https://github.com/TCLamnidis/Sex.DetERRmine) do obliczenia każdego względnego pokrycia, jak również związanych z nimi pasków błędów (Rysunek uzupełniający 1, Uwaga uzupełniająca 3 dla więcej informacji na temat obliczania błędów). Oczekuje się, że samice będą miały wskaźnik x równy 1 i wskaźnik y równy 0, podczas gdy samce będą miały zarówno wskaźnik x, jak i wskaźnik y równy 0,5 (ref. 49).
Autoryzacja
Pierwej potwierdziliśmy, że wzór deaminacji na terminalnych bazach odczytów DNA był charakterystyczny dla starożytnego DNA (1.04-4,5% dla bibliotek nie-UDG i 4,7-9,5% dla bibliotek nie-UDG, patrz Tabela uzupełniająca 1), używając programu mapDamage (wersja 2.0.6)66. Przeprowadziliśmy szereg różnych testów, aby zapewnić autentyczność naszych starożytnych danych. Dla osobników męskich zbadaliśmy polimorfizmy na chromosomie X27 za pomocą pakietu oprogramowania ANGSD (wersja 0.910)67. Ujawniło to solidne szacunki kontaminacji dla 2 osobników płci męskiej z Bolshoy i 1 osobnika płci męskiej z Chalmny-Varre. Wszystkie one były poniżej 1,6% zanieczyszczenia (Tabela 1). Dla osobników żeńskich z tych dwóch miejsc, zauważamy, że są one projektowane blisko samców w przestrzeni PCA (Fig. 2a, Supplementary Figure 3), sugerując ograniczone efekty potencjalnego zanieczyszczenia. Dodatkowo, wygenerowaliśmy zbiór danych przefiltrowanych przez PMD dla wszystkich osobników używając pmdtools (wersja 0.60)30. Filtrowanie PMD zostało wykonane przy użyciu genomu referencyjnego zamaskowanego dla wszystkich pozycji na panelu przechwytującym 1240K, aby uniknąć systematycznych błędów allelicznych na analizowanych pozycjach SNP. Ustawiliśmy próg pmd na 3, który zgodnie z oryginalną publikacją30 skutecznie eliminuje potencjalne współczesne zanieczyszczenia na podstawie braku modyfikacji zasad zgodnych z deaminacją.
Aby zapewnić bardziej ilościowe oszacowanie możliwego zanieczyszczenia u samic, użyliśmy programu ContamMix (wersja 1.0-10)29 do oszacowania zanieczyszczenia mitochondrialnego. Wyodrębniliśmy odczyty mapujące do referencji mitochondrialnej dla każdego ze starożytnych osobników za pomocą samtools (wersja 1.3)65. Następnie wygenerowaliśmy mitochondrialną sekwencję konsensusu dla każdego ze starożytnych osobników przy użyciu Geneious (wersja 10.0.9, http://www.geneious.com,68), i wywołując N dla wszystkich miejsc o pokryciu niższym niż 5. Na koniec wszystkie odczyty mitochondrialne wyrównano do odpowiednich sekwencji konsensusu, używając bwa aln (wersja 0.7.12-r1039)64 z maksymalną liczbą różnic w nasionach (-k) ustawioną na 5 i maksymalną liczbą różnic (-n) na 10, oraz bwa samse (wersja 0.7.12-r1039)64. Wygenerowano wielokrotne wyrównanie sekwencji konsensusu i referencyjnego zestawu 311 genomów mitochondrialnych69 , używając programu mafft (wersja v7.305b)70,71,72 z parametrem –auto. Wyrównanie odczytów, jak również wielokrotne wyrównanie konsensusu i 311 referencyjnych genomów mitochondrialnych zostały następnie dostarczone do ContamMix. Zgłaszamy tutaj tryb a posteriori zanieczyszczenia, wraz z górnymi i dolnymi granicami 95% przedziału posterior (Tabela 1).
Dla dodatkowego uwierzytelnienia, przeprowadziliśmy nadzorowane ADMIXTURE28 (wersja 1.3.0) dla wszystkich próbek, używając sześciu współczesnych populacji (Atayal, French, Kalash, Karitiana, Mbuti i Papuan) jako zdefiniowanych klastrów genetycznych, aby zlokalizować wszelkie duże różnice w klastrach genetycznych wśród osobników z tego samego miejsca (Supplementary Figure 2). Przetestowaliśmy moc tej metody do wykrywania zanieczyszczeń i stwierdziliśmy, że może ona wykryć zanieczyszczenia, które są daleko spokrewnione z przodkami obecnymi w badanych osobnikach już przy wskaźnikach 5-8%, ale brakuje jej mocy do identyfikacji zanieczyszczeń blisko spokrewnionych z badanymi osobnikami (patrz Dodatkowa uwaga 2). Nie zaobserwowaliśmy znaczących różnic (w ramach naszej rozdzielczości) we wzorcach rodowych pomiędzy starożytnymi osobnikami z tego samego miejsca, z wyjątkiem Levänluhta, gdzie indywidualna próbka JK2065 wydaje się pochodzić z innego rodowodu. Przypisaliśmy jej zatem w tym badaniu osobną etykietę populacji, Levänluhta_B.
Wreszcie, używając smartpca, dokonaliśmy projekcji zbiorów danych przefiltrowanych i nie przefiltrowanych przez PMD na ten sam zestaw głównych składowych skonstruowanych na współczesnych populacjach europejskich, aby zapewnić, że starożytne osobniki pozostaną w przybliżeniu na równoważnych pozycjach niezależnie od filtrowania PMD. Było to możliwe dla wszystkich próbek z obróbką UDG-połowy, z wyjątkiem osobników z Levänluhta, które reprezentowały zbyt małe uszkodzenia, by można było zastosować filtrowanie PMD. W odniesieniu do tego miejsca, oparliśmy się na bibliotekach nie-UDG (używając tylko transwersji), które zostały wygenerowane dla trzech osobników użytych w głównej analizie. Stwierdziliśmy, że w ramach oczekiwanego szumu spowodowanego niską liczbą SNP, wszystkie próbki wykazują spójność między przefiltrowanymi i niefiltrowanymi zestawami danych, sugerując niską ilość zanieczyszczeń we wszystkich próbkach (Dodatkowa Figura 3a, b). Cztery dodatkowe osobniki z Levänluhta zostały wykluczone z głównej analizy i z tego testu uwierzytelniającego z powodu niskiego pokrycia (< 15 000 pokrytych SNPs) i braku bibliotek nie-UDG.
Statystyki F
Wszystkie programy wykorzystywane do obliczania statystyk F w tym badaniu można znaleźć jako część pakietu Admixtools (https://github.com/DReichLab/AdmixTools)2,42.
Do wszystkich obliczeń F3 użyliśmy qp3Pop (wersja 412).
Statystyki F4 zostały obliczone przy użyciu qpDstat (wersja 711), a qpAdm (wersja 632)2 został użyty do oszacowania proporcji mieszaniny przy użyciu: Źródła (Lewe Populacje): Nganasan; WHG; EHG; Yamnaya_Samara; LBK_PL. Outgroups (Right Populations): OG1: Mbuti; CHG; Israel_Natufian; Onge; Villabruna; Ami; Mixe. OG2: Mbuti; CHG; Onge; Villabruna; Ami; Mixe. OG3: Mixe; CHG; Israel_Natufian; Villabruna; Onge; Ami. OG4: Mbuti; Israel_Natufian; Onge; Villabruna; Ami; Mixe. OG5: Mbuti; Samara_HG; CHG; Israel_Natufian; Villabruna; Ami.
Aby upewnić się, że zestawy outgroup mają wystarczającą moc do rozróżnienia rodowodów obecnych w źródłach, uruchomiliśmy qpWave (wersja 410) używając tylko źródeł jako lewych populacji i każdego zestawu outgroup jako praw. Wszystkie takie przebiegi qpWave były zgodne tylko z maksymalną rangą, co oznacza, że wszystkie zestawy outgroup miały wystarczającą moc do rozróżnienia pięciu różnych źródeł. Wszystkie modele qpWave i qpAdm były uruchamiane przy użyciu opcji allsnps: YES. Gdy pięciokierunkowe modele domieszkowe dostarczone przez qpAdm miały wartości p powyżej 0.05, ale zawierały nierealne proporcje mieszanki i jednemu ze źródeł przypisano ujemną proporcję mieszanki, uruchomiliśmy model ponownie z tym źródłem, które zostało wykluczone. Dla każdej badanej populacji, jeśli zestaw grup zewnętrznych OG1 nie dał działającego pełnego modelu (p < 0.05), próbowaliśmy alternatywnych zestawów grup zewnętrznych z usuniętą jedną prawą populacją. W ten sposób otrzymaliśmy zestawy outgroup OG2-4. W przypadku Levänluhta, wiele zestawów outgroup wyprodukowało działające modele, które są wymienione w Supplementary Data 4. Model domieszki wymagający minimalnej liczby źródeł, a jednocześnie zapewniający wykonalne proporcje domieszki jest zawsze pokazany. W przypadku PWC ze Szwecji, gdzie żaden z zestawów grup zewnętrznych OG1-4 nie wytworzył modelu roboczego, użyto poprawionego zestawu właściwych populacji (OG5), który zawiera Samara_HG, aby zapewnić większą moc odróżniania przodków łowców-zbieraczy. Preferowaliśmy modele z OG1-4 zamiast OG5, ponieważ OG5 zawiera więcej starożytnych genomów z potencjalnymi błędami w odpowiednich populacjach, co częściej powoduje nieudane modele dla współczesnych populacji testowych. Wykluczone źródła w minimalnych modelach zostały określone jako N/A (Supplementary Data 4). Jeśli można było zrezygnować z Yamnaya lub EHG (jak w przypadku Levänluhta), pokazujemy model, który jest bardziej zgodny z poprzednimi publikacjami3,7,8,45 na Rys. 4, ale pokazujemy oba modele w Danych uzupełniających 4.
Analiza składowych głównych
Użyliśmy smartpca (wersja #16000)73 (https://github.com/DReichLab/EIG) do przeprowadzenia analizy składowych głównych (PCA), używając lsqproject: YES i shrinkmode: YES parameters.
Dla euroazjatyckiej PCA (ryc. 2a), do skonstruowania głównych składowych wykorzystano następujące populacje: Abchazi, Adygei, Albańczycy, Ałtajczycy, Ami, Ormianie, Atayal, Avar.SG, Azeri_WGA, Balkar, Balochi, Baskowie, BeduinA, BeduinB, Białorusini, Borneo, Brahui, Bułgarzy, Buriat.SG, kambodżański, Kanaryjczycy, Czeczenia, Czuwasz, chorwacki, Cypr, czeski, Dai, Daur, Dołgan, Druze, angielski, estoński, parzysty, fiński, francuski, gruziński, grecki, Han, Hazara, Hezhen, węgierski, islandzki, irański, włoski_północny, włoski/południowy, japoński, Żyd_Aszkenazi, Żyd_Georgijski, Żyd_Irański, Żyd_Irakijski, Żyd_Libijski, Żyd_Marokański, Żyd_Tunisjanin, Żyd_Turecki, Żyd_Jemeński, Jordański, Kałasz, Kałmuk, Kinh, Koreański, Kumyk, Kurd_WGA, kirgiski, Lahu, libański, Lezgin, litewski, Makrani, Mala, maltański, Mansi, Miao, Mongola, mordowski, Naxi, Nganasan, Nogai, północnoosetyjski.DG, norweski, orkadyjski, Oroqen, palestyński, pathan, rosyjski, Saami.DG, Saami_WGA, sardyński, saudyjski, szkocki, Selkup, Semende, She, Sherpa.DG, sycylijski, hiszpański, hiszpański_północny, syryjski, tadżycki, tajski, tybetański.DG, Tu, Tubalar, Tujia, Turkish, Turkmen, Tuvinian, Ukrainian, Ulchi, Uygur, Uzbek, Xibo, Yakut, Yi, Yukagir.
Dla zachodnioeuroazjatyckiej PCA (Supplementary Figure 3a, b), następujące populacje zostały użyte do skonstruowania głównych komponentów: Abchazi, Adygei, Albańczycy, Ormianie, Bałkarzy, Baskowie, BeduiniA, BeduiniB, Białorusini, Bułgarzy, Kanaryjczycy_Islandczycy, Czeczeni, Czuwasze, Chorwaci, Cypryjczycy, Czesi, Druzowie, Anglicy, Estończycy, Finowie, Francuzi, Gruzini, Grecy, Węgrzy, Islandczycy, Irańczycy, Włosi_Północni, Włosi_Południowi, Żydzi_Aszkenazyjscy, Żydzi_Georgijscy, Żyd_Irański, Żyd_Irakijski, Żyd_Libijski, Żyd_Marokański, Żyd_Tuński, Żyd_Turecki, Żyd_Jemeński, Jordański, Kumyk, Libański, Leżajski, Litewski, Maltański, Mordowski, Północno_Osetyjski, Norweski, Orkadyjski, Palestyński, Polski, Rosyjski, Sardyński, Saudyjski, Szkocki, Sycylijski, Hiszpański, Hiszpańsko_Północny, Syryjski, Turecki, Ukraiński.
AnalizaADMIXTURE
ADMIXTURE28 przeprowadzono w wersji 1.3.0, po wykluczeniu wariantów o częstości alleli mniejszej niż 0,01 i po przycięciu LD przy użyciu plink (wersja 1.90b3.29)74 z rozmiarem okna 200, rozmiarem kroku 5 i progiem R2 0,5 (https://www.genetics.ucla.edu/software/admixture/download.html). Dla każdej wartości K przeprowadzono pięć powtórzeń, przy czym wartości K mieściły się w zakresie od 2 do 15. Wykorzystano populacje: Ami, Ami.DG, ormiańska, Atayal, Atayal.DG, Bałochi, baskijska, BedouinB, białoruska, Brahmin_Tiwari, Brahui, Chuvash, chorwacka, cypryjska, czeska, angielska, estońska, parzysta, fińska, fińska.DG, francuski, grecki, GujaratiB, Hadza, Han, węgierski, islandzki, Kalash, Karitiana, litewski, Makrani, Mala, Mansi, Mansi.DG, Mari.SG, Mbuti, Mbuti.DG, Mixe, mordowski, Nganasan, norweski, Onge, orkadyjski, papuaski, Pima, rosyjski, Saami.DG, ModernSaami, sardyński, szkocki, Selkup, hiszpański, ukraiński, Ulchi, Yoruba, ALPC_Hungary_MN, Baalberge_MN, Baltic_BA, Baltic_CCC, Baltic_CWC, Baltic_LN, BolshoyOleniOstrov, Bu_kk_Culture_Hungary_MN, ChalmnyVarre, CHG, EHG, Esperstedt_MN, Ganj_Dareh_Iran_Neolit, Hungary_MN, Hungary_Neolithic, Iran_Chalcolithic, JK2065, Koros_Hungary_PL, Kunda, Latvia_HG3, Latvia_MN1, LBK_PL, LBK_Hungary_PL, Levanluhta, Narva, PWC_Sweden_NHG.SG, Scandinavia_LNBA, SHG, Sweden_HG.SG, TRB, Ukraine_HG1, Ukraine_N1, WHG, Yamnaya_Samara.
Odkryliśmy, że K = 11 skutkuje najniższym błędem walidacji krzyżowej, jak pokazano na rysunku 4b.
Haplotypowanie chromosomalne Y
Przypisaliśmy starożytnych mężczyzn do haplogrup Y za pomocą programu yHaplo (https://github.com/23andMe/yhaplo)75. W skrócie, program ten zapewnia automatyczne przeszukiwanie drzewa haplogrup Y (jak podano w yHaplo, dostęp z ISOGG w dniu 04 stycznia 2016 r.) od korzenia do gałęzi downstream w oparciu o obecność alleli pochodnych i przypisuje najbardziej downstreamową haplogrupę z allelami pochodnymi. Dla około 15 000 Y-chromosomalnych SNPs obecnych zarówno w naszym panelu przechwytującym, jak i w dwóch opublikowanych zestawach danych76,77, losowo próbkowaliśmy pojedynczą bazę i używaliśmy jej jako genotypu haploidalnego. Użyliśmy niestandardowego skryptu do konwersji genotypów EigenStrat do formatu yHaplo. Zgłaszamy haplogrupę przypisaną najdalej w dół, dostarczoną przez program (Tabela 1). Sprawdziliśmy również ręcznie status pochodnych i brak mutacji definiujących wyznaczoną haplogrupę, ponieważ brakujące informacje mogą prowadzić do przedwczesnego zatrzymania automatycznego wyszukiwania.
Mitochondrialne haplotypowanie
Importowaliśmy przycięte odczyty mitochondrialne dla każdego osobnika z jakością mapowania >30 do Geneious (wersja 10.0.9, https://www.geneious.com)68 i dokonaliśmy reasemblacji tych odczytów do genomu referencyjnego RSRS78, używając mappera Geneious, ze średnią czułością i 5 iteracjami. Użyliśmy wbudowanego automatycznego wywoływacza wariantów w Geneious, aby znaleźć polimorfizmy mitochondrialne z minimalnym pokryciem 3 i minimalną częstotliwością wariantów 0,67. Otrzymane warianty zostały wyeksportowane do Excela i ręcznie porównane z SNPs zgłoszonymi w filogenezie mtDNA online (mtDNA tree Build 17, 18 Feb 2016, http://www.phylotree.org/). Pozycje nukleotydów 309.1 C(C), 315.1C, indele AC na 515-522, 16182C, 16183C, 16193.1C(C) i 16519 zostały zamaskowane i nie zostały uwzględnione w naszych wywołaniach haplotypów.
Fenotypowe SNPs
Użyliśmy samtools mpileup (wersja 1.3)65, filtrując dla jakości map (-Q) i baz (-q) 30, wyłączając per-Base Alignment Quality (-B), na przyciętych plikach bam, aby wygenerować pileup odczytów mapujących do zestawu 43 fenotypowych SNPs4,40,41,79, które są częścią naszego panelu przechwytywania genomu. Niestandardowy skrypt Pythona został użyty do przekształcenia pileupu w tabelę zawierającą liczbę odczytów wspierających każdy allel (Supplementary Data 2).
Dostępność kodu
Wszystkie programy opisane po raz pierwszy w tym badaniu są swobodnie dostępne w repozytoriach online. Sex.DetERRmine: https://github.com/TCLamnidis/Sex.DetERRmine
ContaminateGenotypes: https://github.com/TCLamnidis/ContaminateGenotypes
https://github.com/TCLamnidis/ContaminateGenotypes