Az askenázi humán referencia genom összeszerelése és annotációja

A személyi genomprojekt (PGP) részét képező HG002 askenázi személyt választottuk az első olyan humán referencia genom létrehozásához, amelyet egyetlen egyedből állítottunk össze. A PGP az Open Consent Model-t használja, az első valóban nyílt hozzáférésű platformot az egyéni emberi genom, fenotípus és orvosi adatok megosztására . A hozzájárulási folyamat felvilágosítja a potenciális résztvevőket a genomikai adatok megosztásának következményeiről és kockázatairól, valamint arról, hogy mire számíthatnak a részvételtől. A nyílt hozzájárulás lehetővé tette a világ első humán genom referenciaanyagainak létrehozását (a HG002 a NIST 8391 referenciaanyag) a Genome In A Bottle (GIAB) által, amelyet kalibrálásra, genom-összeszerelési módszerek fejlesztésére és laboratóriumi teljesítménymérésekre használnak . A projekthez szükséges összes nyers szekvenciaadatot a GIAB-től szereztük be, ahol a nyilvánosság számára szabadon hozzáférhetőek .

A HG002 genomot három mély lefedettségű adatkészlet kombinációjából állítottuk össze: 249 bp hosszúságú Illumina olvasatok, átlagosan több mint 33 kbp hosszúságú Oxford Nanopore (ONT) olvasatok és kiváló minőségű, átlagosan 9567 bp hosszúságú PacBio “HiFi” olvasatok (1. táblázat).

Először két összeállítást készítettünk: az egyiket az Illumina és az ONT olvasatok felhasználásával, a másodikat pedig mindhárom adatkészlet felhasználásával, beleértve a PacBio HiFi olvasatait is. A PacBio HiFi adatok hozzáadása valamivel több teljes szekvenciát eredményezett az összeállításban (2,99 Gb vs. 2,88 Gb) és lényegesen nagyobb contig N50 méretet (18,2 Mb vs. 4,9 Mb). Ez az Ash1 v0.5-nek nevezett összeállítás volt az alapja minden későbbi finomításnak.

Az összeállítás kromoszómákra való leképezése

Az Ash1 v0.5 összeállítás kromoszóma-hozzárendelésének létrehozásához a GRCh38-hoz való igazításokat használtuk a legtöbb scaffold kromoszómára való leképezéséhez. A “Módszerek” fejezetben leírt lépésekkel fokozatosan javított kromoszómaszintű összeállítások sorozatát hoztuk létre, amelynek eredménye az Ash1 v1.7 lett. Az Ash1 v1.7 nagyobb egybefüggőséggel és kisebb hézagokkal rendelkezik, mint a GRCh38, amint azt a 2. táblázat mutatja. Megjegyzendő, hogy e kromoszómák összeállítása során a GRCh38 szekvenciájának egy kis részét (58,3 Mb, a genom 2%-a) felhasználtuk az Ash1-ben lévő hézagok kitöltésére. Ezek a régiók tartalmaznak néhány nehezen összeilleszthető régiót, amelyeket kézzel kuratáltak a GRCh38 számára. Összességében az Ash1 összes hézagának becsült mérete 82,9 Mbp, míg a GRCh38-ban 84,7 Mbp.p13.

Az összeállítás javítási folyamat részeként az egyik előzetes Ash1-összeállítást (v1.1) átkutattuk a 12 745 jó minőségű, izolált szerkezeti variáns (beillesztések és törlések ≥ 50 bp) után, amelyeket Zook et al. az askenázi trió adatainak a GRCh37 -vel való összehasonlításával azonosított. Az említett tanulmány négy különböző szekvenálási technológiát és több variánshívót használt a variánsok azonosításához és a hamis pozitív eredmények kiszűréséhez. E 12 745 SV közül 5807 homozigóta és 6938 heterozigóta. Azt vártuk, hogy az Ash1 összeállítása majdnem az összes homozigóta variánssal megegyezik. Mivel az Ash1 csak egy haplotípust rögzít, arra számítottunk, hogy a heterozigóta SV-k körülbelül felével fog egyetérteni, feltételezve, hogy az összeszerelő algoritmus véletlenszerűen választ a haplotípusok között, amikor eldönti, hogy melyik variáns kerüljön be a végső konszenzusba. Az 5807 homozigóta variánsból 5284 (91%) volt jelen az egyezés kritériumai alapján (lásd a “Módszerek” fejezetet), és 6938 heterozigóta variánsból 3922 (56,5%) volt jelen. Minden változatot a megfelelő helyen találtunk.

Az Ash1 v1.1-en kis (≤ 5 bp) variánshívásokat is végeztünk, és ezeket összehasonlítottuk a GIAB HG002 v4.0 benchmark variánsokkal, amelyekkel számos szubsztitúciós és indel hibát korrigáltunk (lásd a “Módszerek” szakaszt), így kaptuk meg az Ash 1 v1.2-t. Ezután az Ash1 összeállítását újra összehangoltuk a GRCh38-hoz, újrahívtuk a variánsokat, és ezeket a variánsokat összehasonlítottuk az újonnan kifejlesztett v4.1 GIAB benchmark készlettel. A v4.1 benchmark régiókon belüli variánsok közül az Ash1 variánsok 1 256 458 homozigóta és 1 041 476 heterozigóta SNP-vel, valamint 187 227 homozigóta és 193 524 heterozigóta indellel egyeztek meg. A valódi variánstól 30 bp-en belüli variánshívások kizárása után 79 269 SNP és 17 439 indel maradt, ami (feltételezve, hogy ezek mind hibák az Ash1-ben) körülbelül Q45-ös minőségi értéknek (QV) felel meg a szubsztitúciós hibák esetében. E variánsok többsége (52 191 SNP és 4629 indel) szegmentális duplikációkba esik, ami valószínűleg az Ash1-ben hiányzó duplikációkat vagy a rövid leolvasások általi tökéletlen polírozást jelenti. Összefoglalva, az Ash1 összeállítás minősége nagyon magas, a szubsztitúciós minőség becsült értéke 62, az indel hibaarány pedig 2 per millió bp az ismert szegmentális duplikációk, tandem ismétlődések és homopolimerek kizárása után.

A variánshívás összehasonlítása az Ash1 és a GRCh38 használatával

Az új referencia genomok létrehozásának egyik motivációja az, hogy jobb keretet biztosítanak az emberi szekvenciaadatok elemzéséhez a betegségekkel összefüggő genetikai variánsok keresése során. Például az askenázi zsidókról szóló, teljes genom shotgun (WGS) adatokat gyűjtő tanulmánynak a GRCh38 helyett az askenázi referenciagenomot kell használnia. Mivel a genetikai háttér hasonló, kevesebb variánst kellene találni, ha az Ash1 ellenében keresünk.

Az elvárás tesztelésére a Personal Genome Project egyik férfi résztvevőjétől, a PGP17-től (hu34D5B9) gyűjtöttünk WGS-adatokat. Ez az egyén a PGP-adatbázis szerint becslések szerint 66%-ban askenázi, ami a legmagasabb becsült frakció volt, amely a már szekvenált PGP egyénekből rendelkezésre állt. Letöltöttünk 983 220 918 100 bp hosszúságú olvasatokat (körülbelül 33x-os lefedettség), és Bowtie2 segítségével összehangoltuk őket mind az Ash1-hez, mind a GRCh38-hoz. A leolvasások valamivel nagyobb hányada (3 901 270, 0,5%) igazodott az Ash1-hez, mint a GRCh38-hoz.

Ezután megvizsgáltuk a PGP17 és a két referencia genom közötti összes egynukleotid variánst (SNV-k, lásd a “Módszerek”-t). Az elemzés egyszerűsítése érdekében csak azokat a helyeket vettük figyelembe, ahol a PGP17 homozigóta volt. Az Ash1-gyel való összehasonlításaink során először azonosítottuk az összes SNV-t, majd megvizsgáltuk az eredeti Ash1 leolvasási adatokat annak megállapítására, hogy minden egyes ilyen SNV esetében az Ash1 genom tartalmaz-e olyan eltérő allélt, amely megegyezik a PGP17-gyel.

Az Ash1-gyel nem egyező homozigóta helyek száma a PGP17-ben összesen 1 333 345 volt, szemben az 1 700 364-gyel, amikor a PGP17 homozigóta helyeit a GRCh38-hoz hasonlítottuk (Additional file 1: S1 táblázat). Ezután megnéztük az alapul szolgáló Ash1 leolvasási adatokat az 1,33 millió SNV-helyre vonatkozóan, amelyek eredetileg nem egyeztek, és azt találtuk, hogy ezek körülbelül felénél az Ash1 genom heterozigóta volt, és az egyik allél megegyezett a PGP17-gyel. Ha az SNV-ket azokra a helyekre korlátoztuk, ahol a PGP17 és az Ash1 egyaránt homozigóta (valamint egy nagyon kis számú olyan helyre, ahol az Ash1 két olyan változatot tartalmaz, amelyek mindkettő eltér a PGP17-től), ez még tovább csökkentette az SNV-k teljes számát, 707 756-ra. Így valamivel kevesebb, mint 1 millióval kevesebb homozigóta SNV-t találtunk, amikor az Ash1-et használtuk a PGP17 referenciájaként. Megjegyzendő, hogy az Ash1-hez való igazítás helyett a GRCh38-hoz is igazíthatjuk a leolvasásokat, majd eltávolíthatjuk az ismert populációspecifikus variánsokat. Ez a kétlépéses folyamat, bár bonyolultabb, hasonló eredményeket adhat, kivéve az Ash1 azon régióit, amelyek egyszerűen hiányoznak a GRCh38-ból.

Egyeztetés a közös askenázi változatokkal

Az Ash1 ismert, közös askenázi változatok (GRCh38-hoz viszonyított) tartalmának vizsgálatához megvizsgáltuk a Genome Aggregation Database (gnomAD) alapján egy askenázi populációban magas gyakorisággal jelentett SNV-ket. A GnomAD v3.0 1662 askenázi egyéntől származó rövid leolvasású teljes genom adatokból származó SNV hívásokat tartalmaz. Mivel egyes variánsokat csak ezeknek az egyéneknek egy részhalmazában hívtak meg, csak olyan variánshelyeket vettünk figyelembe, amelyeket legalább 200 személynél jelentettek. Ezután összegyűjtöttük a fő allél SNV-ket, amihez az allélfrekvenciának 0,5 fölött kellett lennie a mintavételezett populációban. Elemzésünket továbbá az egybázisú szubsztitúciókra korlátoztuk. Így 2 008 397 gnomAD SNV-helyet kaptunk, ahol az askenázi főallél különbözött a GRCh38-tól.

A GRCh38 2 008 397 gnomAD-helyéből 1 790 688-at pontosan le tudtunk térképezni az Ash1-re (lásd a “Módszerek” részt). Ezután minden egyes pozícióban összehasonlítottuk a GRCh38 bázisát az askenázi főallél bázisával, és megvizsgáltuk az Ash1 alternatív alléljait is azokon a helyeken, ahol heterozigóta. Azokon a helyeken, ahol a leolvasási adatok azt mutatták, hogy a HG002 heterozigóta, és mind az askenázi fő alléllel, mind a GRCh38 alléllal rendelkezik, szükség esetén kicseréltük az Ash1 bázisát, hogy az megegyezzen a fő alléllal. E cserék után az Ash1 az 1,79 millió hely 88%-án (1 580 866) tartalmazta az askenázi fő allélt. A fennmaradó helyeken az Ash1 vagy a GRCh38 alléllal egyezett meg, mivel a HG002 homozigóta a referenciaallélre (204 729 hely), vagy egy harmadik allélt tartalmazott, amely sem a GRCh38, sem a gnomAD főalléllel nem egyezett meg (5093 hely). Ezeknek a módosításoknak tovább kell csökkenteniük a bejelentett SNV-k számát, amikor egy askenázi egyedet az Ash1-hez igazítunk.

Érdemes megjegyezni, hogy ahogy nő a fő allél gyakorisága a gnomAD askenázi populációban, úgy nő azon helyek aránya is, ahol az Ash1 megfelelt a fő allélnak. Például a gnomAD askenázi populációban > 0,9 allélfrekvenciával rendelkező SNV-k több mint 98%-a szerepel az Ash1-ben (3. táblázat).

Annotáció

Minden referencia genom létfontosságú része az annotáció: a genomban található összes gén és egyéb jellemző összegyűjtése. Annak érdekében, hogy az Ash1 valódi referencia genomként működhessen, a tudományos közösség által használt összes ismert gént feltérképeztük a koordinátarendszerére, azonos génnevek és azonosítók használatával. A GRCh38-hoz számos különböző annotációs adatbázist hoztak létre, az Ash1 esetében pedig a CHESS humán génadatbázist választottuk, mivel az átfogó, és tartalmazza mind a GENCODE, mind a RefSeq, a két másik széles körben használt génadatbázisban szereplő összes fehérjekódoló gént, és mivel megtartja az említett katalógusokban használt azonosítókat. A nem kódoló gének a három adatbázisban különböznek, de a CHESS tartalmazza a legtöbb génlokuszt és izoformát. A CHESS 2.2-es verzióját használtuk, amely 42 167 gént tartalmaz az elsődleges kromoszómákon (kivéve a GRCh38 alternatív scaffoldokat), amelyek közül 20 197 fehérjekódoló.

A gének feltérképezése egyik összeállításból a másikba összetett feladat, különösen a nagyon hasonló, több kópiát tartalmazó géncsaládokban előforduló gének esetében. A feladat még összetettebb, ha a két összeállítás különböző egyedeket képvisel (és nem egyszerűen ugyanazon egyed különböző összeállításait), mivel az egyedek között egynukleotid eltérések, inszerciók, deléciók, átrendeződések és valódi eltérések vannak a kópiaszámban. Olyan módszerre volt szükségünk, amely mindezekkel a lehetséges különbségekkel szemben robusztus.

A probléma megoldására a nemrégiben kifejlesztett Liftoff térképező eszközt használtuk, amely kísérleteink során az egyetlen olyan eszköz volt, amely szinte az összes emberi gént le tudta térképezni egyik egyedről a másikra. A Liftoff az összes gént és transzkriptumot egy genomból veszi, és kromoszómáról kromoszómára feltérképezi őket egy másik genomra. Azokat a géneket, amelyeket nem sikerül ugyanarra a kromoszómára leképezni, a Liftoff megpróbálja kromoszómák között leképezni. Más eszközökkel ellentétben nem támaszkodik a teljes genom összehangolásából felépített részletes térképre, hanem minden egyes gént külön-külön térképez le. A géneket a transzkriptumok szintjén igazítja, beleértve az intronokat is, így a feldolgozott pszeudogének nem lesznek tévesen génekként azonosítva.

A GRCh38 elsődleges kromoszómáin lévő 42 167 génloci mind a 310 901 transzkriptumát megpróbáltuk feltérképezni az Ash1-re. Összesen 42 070 génlókuszból 309 900 (99,7%) transzkriptet sikerült leképeznünk a fő kromoszómákra (Additional file 1: S2 táblázat). Egy transzkriptet akkor tekintettünk sikeresen leképezettnek, ha az Ash1-ben található mRNS-szekvencia legalább 50%-ban olyan hosszú, mint a GRCh38-on található mRNS-szekvencia. A transzkriptek túlnyomó többsége azonban jóval meghaladta ezt a küszöbértéket, a transzkriptek 99%-a 95%-nál nagyobb vagy azzal egyenlő lefedettséggel térképeződött le (Additional file 2: S2 ábra). A leképezett transzkriptek szekvenciaazonossága hasonlóan magas, a transzkriptek 99%-a 94%-nál nagyobb vagy azzal egyenlő szekvenciaazonossággal térképeződött le (Additional file 2: S3 ábra).

Translokált gének

A legalább egy sikeresen leképezett izoformával rendelkező gének közül 42 059 (99,7%) az Ash1 azonos kromoszómájának megfelelő helyére térképeződött. Az eredetileg sikertelenül térképezett 108 gén közül 11 gén 7 különböző kromoszómára térképeződött 7 különböző blokkban (lásd a 4. táblázatot), ami a két genom közötti transzlokációra utal. Érdekes módon a transzlokációban érintett 22 helyből 16 szubtelomerikus régiókban volt, ami 8 olyan párban fordult elő, ahol mindkét hely telomerek közelében volt. Ez összhangban van korábbi tanulmányokkal, amelyek arról számoltak be, hogy a telomereket és szubtelomereket érintő átrendeződések a transzlokáció gyakori formája lehetnek az emberben .

A 15. és 20. kromoszóma közötti transzlokációt, amely a 4. táblázatban szereplő három gént tartalmazza, a GRCh38 és az Ash1 közötti összehangolást közelebbről megvizsgálva vizsgáltuk. A transzlokáció mindkét kromoszóma telomerjénél található, az Ash1 chr20 65,079,275 és 65,109,824 (30,549 bp), valamint a GRCh39 chr15 101,950,338 és 101,980,928 (30,590 bp) közötti pozíciójában. A transzlokáció megerősítése érdekében összehangoltunk egy független, nagyon hosszú PacBio leolvasásokból álló sorozatot, amelyek mind a HG002-ből származnak, az Ash1 v1.7 összeállítással (lásd a “Módszerek” részt), és kiértékeltük a töréspont körüli régiót a chr20-on. Ezek az illesztések mély, konzisztens, több kilobázist átfogó lefedettséget mutatnak a töréspont mindkét oldalán, ami alátámasztja az Ash1 összeállítás helyességét (1. ábra).

A hosszú PacBio leolvasásoknak az Ash1 genomhoz való illesztését bemutató pillanatfelvétel, a 20-as kromoszómán lévő hely bal végére (65,079,275 pozíció) központosítva, ahol a 15-ös (GRCh38) és a 20-as (Ash1) kromoszóma között transzlokáció történt. Az ábra felső része a chr20 koordinátáit mutatja. Alatta a leolvasási lefedettség hisztogramja látható, az egyes leolvasások pedig az ábra alsó részét töltik ki. A leolvasásokban található indelek, amelyeket az egyes leolvasásokon színes sávok jeleznek, a hosszú leolvasások viszonylag magas hibaarányából adódnak

Hiányos gének

Százkilencvenkét génnek nem sikerült teljesen leképezni a GRCh38-ról az Ash1-re, és további 32 gén csak részben (az 50%-os lefedettségi küszöb alatt) térképeződött le, amint azt az 5. táblázat mutatja. A nem vagy csak részben leképeződő gének mindegyike több génből álló géncsalád tagja volt, és minden esetben a hiányzó génnek legalább egy másik példánya is jelen volt az Ash1-ben, átlagosan 98,5%-os azonosság mellett. Tehát egyáltalán nincs olyan eset, amikor egy gén jelen van a GRCh38-ban, és teljesen hiányzik az Ash1-ből; az 5. táblázatban bemutatott gének olyan eseteket jelentenek, amikor az Ash1-ben kevesebb tagja van egy több génből álló családnak. Három további gén (2 fehérjekódoló, 1 lncRNS) két nem elhelyezett kontighoz térképeződött, ami útmutatást nyújt majd e kontigok elhelyezéséhez az Ash1 assembly jövőbeli kiadásaiban.

A gének Ash1-re történő leképezése után kivontuk a kódoló szekvenciákat a teljesen leképezett (100%-os lefedettségű) transzkriptekből, összehangoltuk őket a GRCh38 kódoló szekvenciáival, és a GRCh38-hoz viszonyított variánsokat neveztük meg (lásd a “Módszerek” fejezetet). Ezekben a fehérjéket kódoló transzkriptekben lévő 35 513 365 bp-n belül 20 864 egynukleotid variánst és indelt találtunk. E variánsok közül tizennégyezer-négyszázkilencvenkilenc a GIAB “hívható” régióiba esett a nagy megbízhatóságú variánsok számára, bár ezek közül 3963 a GIAB “nehéz” ismétlődő régióiban volt, amelyek esetében az összehangolások gyakran nem egyértelműek. A 10 536 olyan variánsból, amely nem ezekben a nehéz régiókban volt, 10 456 (99,2%) megegyezett a GIAB magas megbízhatóságú variánskészletével. A nehéz régiókban 3804/3963 (96,0%) megegyezett a GIAB-készlettel.

Ezután az összes fehérjekódoló szekvenciára vonatkozóan megjegyeztük a variánsok és a hiányos leképezés által okozott aminosavváltozásokat. A 20 197 génből származó 124 238 fehérjekódoló transzkriptumból 92 600 (74,5%) 100%-ban azonos fehérjeszekvenciával rendelkezik. További 26 566 (21,4%) legalább egy aminosavváltozást tartalmaz, de azonos hosszúságú fehérjéket eredményez, és 1561 (1,3%) olyan keretmegőrző mutációval rendelkezik, amely egy vagy több aminosavat illeszt be vagy töröl, a fehérje többi részét változatlanul hagyva. A 6. táblázat a fehérjeszekvenciák összes változására vonatkozó statisztikákat mutatja be. Ha egy fehérje 1-nél több variánssal rendelkezett, akkor azt a legkövetkezetesebb variáns alá számoltuk, azaz ha egy fehérje missense variánssal és korai stop kodonnal rendelkezett, akkor a “stop gained” csoportba soroltuk.

által okozott korai stop kodonokra utal

Kiemelten érdekesek azok a transzkriptek, amelyek variánsai jelentősen megzavarják a fehérjeszekvenciát, és funkcióvesztést eredményezhetnek. Ezek közé tartoznak a frameshift (2158), stopvesztés (58), stopnyereség (416), startvesztés (58) vagy hiányos leképezés miatti csonkulás (564) által érintett transzkriptek. Ezek a megszakadt izoformák 885 génlokuszt képviselnek; azonban e gének közül 505-nek van legalább 1 másik izoformája, amelyet nem érint a megszakadó variáns. Így 380 olyan gén marad, amelynek minden izoformája legalább egy megszakítással rendelkezik; a teljes lista az 1. kiegészítő fájlban található: S1. táblázat.