Assemblage en annotatie van een menselijk referentiegenoom van de Asjkenazische stam

Voor het maken van het eerste menselijke referentiegenoom dat van een enkel individu wordt geassembleerd, hebben we HG002 gekozen, een Asjkenazisch individu dat deel uitmaakt van het Personal Genome Project (PGP). Het PGP maakt gebruik van het Open Toestemmingsmodel, het eerste werkelijk open-toegangsplatform voor het delen van individuele menselijke genoom-, fenotype- en medische gegevens. Bij het toestemmingsproces worden potentiële deelnemers voorgelicht over de implicaties en risico’s van het delen van genomische gegevens en over wat zij van hun deelname kunnen verwachten. Dankzij de open toestemming konden ’s werelds eerste referentiematerialen voor het menselijk genoom worden gecreëerd (HG002 is NIST-referentiemateriaal 8391) van Genome In A Bottle (GIAB), dat wordt gebruikt voor kalibratie, ontwikkeling van methoden voor genoomassemblage, en metingen van de prestaties van laboratoria . Alle ruwe sequentiegegevens voor dit project werden verkregen van GIAB, waar ze vrij beschikbaar zijn voor het publiek.

Wij assembleerden het HG002-genoom uit een combinatie van drie diepe-dekking datasets: 249-bp Illumina leest, Oxford Nanopore (ONT) leest gemiddeld meer dan 33 kbp in lengte, en hoge-kwaliteit PacBio “HiFi” leest gemiddeld 9567 bp (tabel 1).

We hebben in eerste instantie twee assemblages: een met behulp van Illumina en ONT leest, en een tweede met behulp van alle drie de datasets, met inbegrip van de PacBio HiFi leest. De toevoeging van PacBio HiFi gegevens resulteerde in iets meer totale sequentie in de montage (2,99 Gb vs 2,88 Gb) en een aanzienlijk grotere contig N50 grootte (18,2 Mb vs 4,9 Mb). Deze assemblage, aangeduid als Ash1 v0.5, was de basis voor alle latere refinements.

Het in kaart brengen van de assemblage op chromosomen

Om chromosoom toewijzingen voor de Ash1 v0.5 assemblage te maken, gebruikten we uitlijningen naar GRCh38 om de meeste van de scaffolds kaart op chromosomen. De stappen beschreven in de “Methoden” sectie genereerde een reeks van geleidelijk verbeterde chromosoom schaal assemblages, resulterend in Ash1 v1.7. Ash1 v1.7 heeft een grotere contiguïteit en kleinere lacunes dan GRCh38, zoals weergegeven in tabel 2. Merk op dat in het proces van de bouw van deze chromosomen, een kleine hoeveelheid GRCh38 sequentie (58,3 Mb, 2% van het genoom) werd gebruikt om gaten te vullen in Ash1. Deze regio’s omvatten een aantal moeilijk te monteren regio’s die handmatig zijn gecureerd voor GRCh38. In totaal is de geschatte grootte van alle hiaten in Ash1 82,9 Mbp, versus 84,7 Mbp in GRCh38.p13.

Als onderdeel van het proces ter verbetering van de assemblage hebben we een van de voorlopige Ash1-assemblages (v1.1) doorzocht op de 12.745 geïsoleerde structurele varianten van hoge kwaliteit (invoegingen en verwijderingen ≥ 50 bp) die Zook et al. identificeerden door de Ashkenazi-triogegevens te vergelijken met GRCh37 . Die studie gebruikte vier verschillende sequencing technologieën en meerdere variant callers om varianten te identificeren en vals-positieven uit te filteren. Van deze 12.745 SVs zijn 5807 homozygoot en 6938 heterozygoot. We verwachtten dat de Ash1 assemblage overeen zou komen met bijna alle homozygote varianten. Omdat Ash1 slechts één haplotype vastlegt, verwachtten we dat het met ongeveer de helft van de heterozygote SV’s overeen zou komen, ervan uitgaande dat het assemblage-algoritme willekeurig tussen de haplotypen koos bij de beslissing welke variant in de uiteindelijke consensus zou worden opgenomen. Van de 5807 homozygote varianten, waren er 5284 (91%) aanwezig volgens onze match criteria (zie de “Methoden” sectie), en 3922 (56.5%) van de 6938 heterozygote varianten waren aanwezig. Alle varianten werden op de juiste plaats gevonden.

We hebben ook kleine (≤ 5 bp) variant calls gedaan op Ash1 v1.1 en deze vergeleken met de HG002 v4.0 benchmark varianten van GIAB, die we gebruikt hebben om talrijke substitutie en indel fouten te corrigeren (zie de “Methoden” sectie), wat Ash1 v1.2 opleverde. Vervolgens hebben we de Ash1 assemblage opnieuw uitgelijnd met GRCh38, de varianten hernoemd, en deze varianten gebenchmarkt tegen de nieuw ontwikkelde v4.1 GIAB benchmark set. Van de varianten binnen de v4.1 benchmark regio’s, kwamen de Ash1 varianten overeen met 1.256.458 homozygote en 1.041.476 heterozygote SNPs, en 187.227 homozygote en 193.524 heterozygote indels. Na het uitsluiten van variant calls binnen 30 bp van een echte variant, bleven 79.269 SNPs en 17.439 indels over, wat (ervan uitgaande dat dit allemaal fouten in Ash1 zijn) overeenkomt met een kwaliteitswaarde (QV) van ongeveer Q45 voor substitutiefouten. De meeste van deze varianten (52.191 SNPs en 4629 indels) vallen in segmentale duplicaties, die mogelijk ontbrekende duplicaties in Ash1 vertegenwoordigen of onvolmaakte polijsting door korte lezingen. Samenvattend is de kwaliteit van de Ash1-assemblage zeer hoog, met een geschatte substitutiekwaliteitswaarde van 62 en een indelfoutenpercentage van 2 per miljoen bp na het uitsluiten van bekende segmentale duplicaties, tandemherhalingen en homopolymeren.

Vergelijking van variantenoproeping met behulp van Ash1 versus GRCh38

Eén van de motivaties voor het creëren van nieuwe referentie-genomen is dat ze een beter kader bieden voor het analyseren van menselijke sequentiegegevens bij het zoeken naar genetische varianten die verband houden met ziekte. Bijvoorbeeld, een studie van Ashkenazi Joden die whole-genome shotgun (WGS) gegevens verzamelde zou een Ashkenazi referentiegenoom moeten gebruiken in plaats van GRCh38. Omdat de genetische achtergrond vergelijkbaar is, zouden er minder varianten moeten worden gevonden bij het zoeken tegen Ash1.

Om deze verwachting te testen, hebben we WGS-gegevens verzameld van een mannelijke deelnemer aan het Personal Genome Project, PGP17 (hu34D5B9). Dit individu wordt geschat op 66% Ashkenazi volgens de PGP-database, wat de hoogste geschatte fractie was die beschikbaar was van reeds gesequenteerde PGP-individuen. We downloadden 983.220.918.100-bp leest (ongeveer 33x dekking) en uitgelijnd ze aan zowel Ash1 en GRCh38 met behulp van Bowtie2 . Een iets hogere fractie van gelezen (3.901.270, 0,5%) uitgelijnd op Ash1 dan op GRCh38.

We onderzochten vervolgens alle single-nucleotide varianten (SNV’s, zie de “Methoden”) tussen PGP17 en elk van de twee referentie-genomen. Om de analyse te vereenvoudigen, hebben we alleen gekeken naar locaties waar PGP17 homozygoot was. In onze vergelijkingen met Ash1 identificeerden we eerst alle SNV’s en onderzochten vervolgens de oorspronkelijke Ash1 leesgegevens om te bepalen of, voor elk van die SNV’s, het Ash1 genoom een ander allel bevatte dat overeenkwam met PGP17.

In totaal was het aantal homozygote locaties in PGP17 dat niet overeenkwam met Ash1 1.333.345, versus 1.700.364 toen we homozygote locaties in PGP17 vergeleken met GRCh38 (Additional file 1: Tabel S1). Vervolgens keken we naar de onderliggende gegevens van de 1,33 miljoen SNV’s die aanvankelijk niet overeenkwamen, en ontdekten dat voor ongeveer de helft daarvan het genoom van Ash1 heterozygoot was, met één allel dat overeenkwam met PGP17. Als we de SNV’s beperken tot locaties waar PGP17 en Ash1 beide homozygoot zijn (plus een zeer klein aantal locaties waar Ash1 twee varianten bevat die beide verschillen van PGP17), reduceert dit het totale aantal SNV’s nog verder, tot 707.756. We vonden dus iets minder dan 1 miljoen homozygote SNV’s wanneer we Ash1 als referentie voor PGP17 gebruikten. Merk op dat in plaats van uitlijning op Ash1, men in plaats daarvan zou kunnen de leest uit te lijnen naar GRCh38 en verwijder vervolgens bekende populatie-specifieke varianten. Deze twee-staps proces, hoewel complexer, zou kunnen opleveren vergelijkbare resultaten, behalve voor regio’s van Ash1 die gewoon ontbreken in GRCh38.

Vergelijking met gemeenschappelijke Ashkenazi varianten

Om te onderzoeken de mate waarin Ash1 bekende, gemeenschappelijke Ashkenazi varianten (ten opzichte van GRCh38) bevat, onderzochten we SNV’s gemeld bij hoge frequentie in een Ashkenazi bevolking van de Genome Aggregation Database (gnomAD) . GnomAD v3.0 bevat SNV-oproepen uit short-read genoomwijde gegevens van 1662 Ashkenazi individuen. Omdat sommige varianten slechts in een subset van deze individuen werden opgeroepen, hebben we alleen gekeken naar variant sites die in een minimum van 200 personen werden gerapporteerd. Vervolgens verzamelden we de belangrijkste allel SNV’s, waarbij de allelfrequentie hoger dan 0,5 moest zijn in de bemonsterde populatie. We beperkten onze analyse verder tot single-base substituties. Dit gaf ons 2.008.397 gnomAD SNV sites waar het Ashkenazi major allel verschilde van GRCh38.

We waren in staat om 1.790.688 van de 2.008.397 gnomAD sites van GRCh38 nauwkeurig in kaart te brengen op Ash1 (zie de sectie “Methoden”). Vervolgens vergeleken we de GRCh38 basis met de Ashkenazi major allel basis op elke positie, en we onderzochten ook de alternatieve allelen in Ash1 op plaatsen waar het heterozygoot is. Voor plaatsen waar de afgelezen gegevens aantoonden dat HG002 heterozygoot was en zowel het Ashkenazi major allel als het GRCh38 allel had, hebben we, indien nodig, de Ash1 base vervangen, om er zeker van te zijn dat deze overeenkwam met het major allel. Na deze vervangingen bevatte Ash1 het belangrijkste Ashkenazi allel op 88% (1.580.866) van de 1,79 miljoen locaties. Op de resterende plaatsen kwam Ash1 overeen met het GRCh38 allel omdat HG002 homozygoot is voor het referentie allel (204.729 plaatsen), of het bevatte een derde allel dat noch met GRCh38 noch met het gnomAD major allel overeenkwam (5093 plaatsen). Deze modificaties zouden het aantal gerapporteerde SNV’s verder moeten verminderen wanneer een Ashkenazi-individu wordt uitgelijnd met Ash1.

Opgemerkt moet worden dat, naarmate de frequentie van het hoofdallel in de gnomAD Ashkenazi-populatie toeneemt, het aandeel locaties waar Ash1 overeenkwam met het hoofdallel ook toeneemt. Van de SNV’s met een allelfrequentie > 0,9 in de gnomAD Ashkenazi-populatie is bijvoorbeeld meer dan 98% vertegenwoordigd in Ash1 (tabel 3).

Annotatie

Een vitaal onderdeel van elk referentiegenoom is annotatie: de verzameling van alle genen en andere kenmerken die in het genoom worden aangetroffen. Om Ash1 als een echt referentiegenoom te laten fungeren, hebben wij alle bekende genen die door de wetenschappelijke gemeenschap worden gebruikt, in kaart gebracht op zijn coördinatensysteem, met gebruikmaking van dezelfde gennamen en identificatoren. Er zijn verschillende annotatie databases gemaakt voor GRCh38, en voor Ash1 hebben we gekozen voor de CHESS menselijke gen-database, omdat deze veelomvattend is en alle eiwit-coderende genen bevat in zowel GENCODE als RefSeq, de twee andere veel gebruikte gen-databases, en omdat het de identifiers behoudt die in die catalogi worden gebruikt. De niet-coderende genen verschillen tussen de drie databases, maar CHESS heeft het grootste aantal genloci en isovormen. Wij gebruikten CHESS versie 2.2, die 42.167 genen op de primaire chromosomen bevat (exclusief de GRCh38 alternatieve scaffolds), waarvan 20.197 eiwitcoderend zijn.

Het in kaart brengen van genen van de ene assemblage naar de andere is een complexe taak, vooral voor genen die voorkomen in zeer vergelijkbare, multi-kopieën genfamilies. De taak is nog complexer wanneer de twee assemblages verschillende individuen vertegenwoordigen (in plaats van eenvoudigweg verschillende assemblages van hetzelfde individu), als gevolg van de aanwezigheid van single-nucleotide verschillen, inserties, deleties, herschikkingen, en echte variaties in kopie-aantal tussen de individuen. We hadden een methode nodig die robuust zou zijn in het gezicht van al deze potentiële verschillen.

Om dit probleem aan te pakken, gebruikten we het recent ontwikkelde Liftoff mapping tool, dat in onze experimenten het enige instrument was dat bijna alle menselijke genen van het ene individu naar het andere in kaart kon brengen. Liftoff neemt alle genen en transcripten van een genoom en brengt ze, chromosoom voor chromosoom, in kaart in een ander genoom. Voor alle genen die niet op hetzelfde chromosoom in kaart kunnen worden gebracht, probeert Liftoff ze over chromosomen heen in kaart te brengen. In tegenstelling tot andere hulpmiddelen vertrouwt het niet op een gedetailleerde kaart die is opgebouwd uit een uitlijning van het hele genoom, maar brengt het elk gen afzonderlijk in kaart. Genen worden uitgelijnd op transcriptniveau, inclusief introns, zodat verwerkte pseudogenen niet ten onrechte als genen zullen worden geïdentificeerd.

We hebben geprobeerd om alle 310.901 transcripten van 42.167 genloci op de primaire chromosomen in GRCh38 naar Ash1 in kaart te brengen. In totaal hebben we met succes in kaart gebracht 309.900 (99,7%) transcripten van 42.070 gen loci op de belangrijkste chromosomen (Additional file 1: Tabel S2). We beschouwden een transcript als succesvol in kaart gebracht als de mRNA-sequentie in Ash1 is ten minste 50% zo lang als de mRNA-sequentie op GRCh38. Echter, de overgrote meerderheid van de transcripten veel meer dan deze drempel, met 99% van de transcripten in kaart gebracht op een dekking groter dan of gelijk aan 95% (Additional file 2: Figuur S2). De sequentie-identiteit van de in kaart gebrachte transcripten is vergelijkbaar hoog, met 99% van de transcripten in kaart brengen met een sequentie-identiteit groter dan of gelijk aan 94% (Additional file 2: Figuur S3).

Getranslokeerde genen

Van die genen met ten minste één met succes in kaart gebracht isovorm, 42.059 (99,7%) in kaart gebracht op de overeenkomstige locaties op hetzelfde chromosoom in Ash1. Van de 108 genen die aanvankelijk niet in kaart werden gebracht, werden 11 genen in kaart gebracht op een ander chromosoom in 7 verschillende blokken (weergegeven in tabel 4), wat een translocatie tussen de twee genomen suggereert. Interessant is dat 16 van de 22 locaties die betrokken waren bij de translocaties zich in subtelomere regio’s bevonden, die zich voordeden in 8 paren waar beide locaties zich in de buurt van telomeren bevonden. Dit is consistent met eerdere studies die meldden dat herschikkingen waarbij telomeren en subtelomeren betrokken zijn, een veel voorkomende vorm van translocatie bij de mens kunnen zijn.

We onderzochten de translocatie tussen chromosomen 15 en 20, die drie van de genen in Tabel 4 bevat, door de uitlijning tussen GRCh38 en Ash1 nauwkeuriger te bekijken. De translocatie bevindt zich op het telomeer van beide chromosomen, van positie 65.079.275 tot 65.109.824 (30.549 bp) van Ash1 chr20 en 101.950.338 tot 101.980.928 (30.590 bp) van GRCh39 chr15. Om de translocatie te bevestigen, hebben wij een onafhankelijke reeks zeer lange PacBio-reads, alle van HG002, uitgelijnd met de Ash1 v1.7 assemblage (zie de sectie “Methoden”) en de regio rond het breekpunt op chr20 geëvalueerd. Deze uitlijningen laten een diepe, consistente dekking zien die zich uitstrekt over vele kilobases aan beide zijden van het breekpunt, wat de juistheid van de Ash1-assemblage ondersteunt (fig. 1).

Snapshot met uitlijningen van lange PacBio-reads aan het Ash1-genoom, gecentreerd op het linkeruiteinde van de locatie in chromosoom 20 (positie 65.079.275) waar een translocatie optrad tussen chromosoom 15 (GRCh38) en 20 (Ash1). Het bovenste deel van de figuur toont de coördinaten op chr20. Daaronder is een histogram van de leesdekking, en de individuele lezingen vullen het onderste deel van de figuur. De indels in de leest, weergegeven als gekleurde balken op elke leest, zijn te wijten aan de relatief hoge foutenmarge van de lange leest

Missing genen

Tweeënzestig genen niet volledig in kaart van GRCh38 op Ash1, en nog eens 32 genen in kaart slechts gedeeltelijk (onder de 50% dekking drempel), zoals weergegeven in tabel 5. Alle genen die niet in kaart gebracht van GRCh38 op Ash1, en nog eens 32 genen in kaart slechts gedeeltelijk (onder de 50% dekking drempel), zoals weergegeven in tabel 5. Alle genen die niet of slechts gedeeltelijk gekarteerd werden, waren leden van multi-gen families, en in elk geval was er minstens één andere kopie van het ontbrekende gen aanwezig in Ash1, met een gemiddelde identiteit van 98.5%. Er zijn dus helemaal geen gevallen van een gen dat in GRCh38 aanwezig is en dat in Ash1 geheel ontbreekt; de genen in Tabel 5 zijn gevallen waarin Ash1 minder leden van een multi-gen familie heeft. Drie extra genen (2 eiwitcoderend, 1 lncRNA) zijn in kaart gebracht met twee niet-geplaatste contigs, die een leidraad zullen vormen voor het plaatsen van deze contigs in toekomstige versies van de Ash1-assemblage.

Na het in kaart brengen van de genen op Ash1, hebben we de coderende sequenties geëxtraheerd van transcripten die volledig in kaart waren gebracht (dekking gelijk aan 100%), deze uitgelijnd met de coderende sequenties van GRCh38, en varianten genoemd ten opzichte van GRCh38 (zie de sectie “Methoden”). Binnen de 35.513.365 bp in deze eiwit-coderende transcripten, vonden we 20.864 single-nucleotide varianten en indels. Veertienduizend vierhonderd negenennegentig van deze varianten vielen binnen de GIAB “opvraagbare” regio’s voor varianten met hoge betrouwbaarheid, hoewel 3963 van deze varianten in GIAB “moeilijke” repetitieve regio’s vielen, waarvoor uitlijningen vaak dubbelzinnig zijn. Van de 10.536 varianten die zich niet in deze moeilijke regio’s bevonden, kwamen er 10.456 (99,2%) overeen met de GIAB “highconfidence” variantenreeks. In de moeilijke regio’s, kwamen 3804/3963 (96,0%) overeen met de GIAB set.

We annoteerden vervolgens de veranderingen in aminozuren veroorzaakt door varianten en onvolledige mapping voor alle eiwit-coderende sequenties. Van de 124.238 eiwit-coderende transcripten van 20.197 genen, hebben 92.600 (74,5%) 100% identieke eiwit-sequenties. Nog eens 26.566 (21,4%) hebben ten minste één aminozuurverandering, maar leveren eiwitten op met dezelfde lengte, en 1561 (1,3%) hebben frame-behoudende mutaties die één of meer aminozuren invoegen of verwijderen, waarbij de rest van het eiwit ongewijzigd blijft. Tabel 6 toont statistieken over alle veranderingen in eiwitsequenties. Als een eiwit meer dan 1 variant had, telden we het onder de meest ingrijpende variant, d.w.z. als een eiwit een missense-variant en een voortijdig stopcodon had, namen we het op in de “stop gained”-groep.

Van bijzonder belang zijn de transcripten met varianten die de eiwitsequentie aanzienlijk verstoren en kunnen resulteren in functieverlies. Hiertoe behoren transcripten met een frameshift (2158), stop loss (58), stop gain (416), start loss (58), of truncatie als gevolg van onvolledige mapping (564). Deze verstoorde isovormen vertegenwoordigen 885 gen loci; echter, 505 van deze genen hebben tenminste 1 andere isovorm die niet beïnvloed is door een verstorende variant. Dit laat 380 genen over waarin alle isovormen ten minste één verstoring hebben; de volledige lijst is te vinden in Additional file 1: Table S1.