Montering och annotering av ett Ashkenazi-referensgenom

För att skapa det första mänskliga referensgenomet som ska monteras från en enskild individ valde vi HG002, en Ashkenazi-individ som ingår i Personal Genome Project (PGP). PGP använder sig av Open Consent Model, den första plattformen med verkligt öppen tillgång för delning av individuella mänskliga genom-, fenotyp- och medicinska data . I samtyckesprocessen utbildas potentiella deltagare om konsekvenserna och riskerna med att dela genomdata och om vad de kan förvänta sig av sitt deltagande. Det öppna samtycket har gjort det möjligt att skapa världens första referensmaterial för humangenom (HG002 är NIST:s referensmaterial 8391) från Genome In A Bottle (GIAB), som används för kalibrering, utveckling av metoder för sammansättning av arvsmassor och mätningar av laboratoriernas prestanda . Alla råsekvensdata för detta projekt erhölls från GIAB, där de är fritt tillgängliga för allmänheten .

Vi sammanställde HG002-genomet från en kombination av tre djuptäckande datamängder: Illumina-avläsningar med 249 bp, Oxford Nanopore-avläsningar (ONT) med en genomsnittlig längd på över 33 kbp och högkvalitativa PacBio-”HiFi”-avläsningar med en genomsnittlig längd på 9567 bp (tabell 1).

Vi skapade inledningsvis två sammansättningar: en med hjälp av Illumina- och ONT-avläsningar och en andra med hjälp av alla tre datauppsättningarna, inklusive PacBio HiFiavläsningarna. Tillägget av PacBio HiFi-data resulterade i något mer total sekvens i sammansättningen (2,99 Gb jämfört med 2,88 Gb) och en betydligt större contig N50-storlek (18,2 Mb jämfört med 4,9 Mb). Denna samling, med beteckningen Ash1 v0.5, låg till grund för alla efterföljande förfiningar.

Mappning av samlingen på kromosomer

För att skapa kromosomtilldelningar för Ash1 v0.5-samlingen använde vi anpassningar till GRCh38 för att mappa de flesta av scaffolds på kromosomer. De steg som beskrivs i avsnittet ”Metoder” genererade en serie gradvis förbättrade kromosomskaliga sammansättningar, vilket resulterade i Ash1 v1.7. Ash1 v1.7 har större kontiguitet och mindre luckor än GRCh38, vilket framgår av tabell 2. Observera att under uppbyggnaden av dessa kromosomer användes en liten mängd GRCh38-sekvenser (58,3 Mb, 2 % av genomet) för att fylla luckor i Ash1. Dessa regioner omfattar vissa regioner som är svåra att sammanfoga och som har kurerats manuellt för GRCh38. Totalt är den uppskattade storleken på alla luckor i Ash1 82,9 Mbp, jämfört med 84,7 Mbp i GRCh38.p13.

Som en del av förbättringen av sammansättningen sökte vi i en av de preliminära Ash1-sammansättningarna (v1.1) efter de 12 745 högkvalitativa, isolerade strukturella varianter (infogningar och borttagningar ≥ 50 bp) som Zook et al. identifierade genom att jämföra data från den ashkenaziska trion med GRCh37. I den studien användes fyra olika sekvenseringstekniker och flera variant callers för att identifiera varianter och filtrera bort falska positiva varianter. Av dessa 12 745 SV är 5807 homozygota och 6938 heterozygota. Vi förväntade oss att Ash1-samlingen skulle stämma överens med nästan alla homozygota varianter. Eftersom Ash1 bara fångar en haplotyp förväntade vi oss att den skulle stämma överens med ungefär hälften av de heterozygota SV:erna, om man antar att sammansättningsalgoritmen valde slumpmässigt mellan haplotyperna när den bestämde vilken variant som skulle inkluderas i det slutliga konsensusetablissemanget. Av de 5807 homozygota varianterna fanns 5284 (91 %) med hjälp av våra matchningskriterier (se avsnittet ”Metoder”), och 3922 (56,5 %) av 6938 heterozygota varianter fanns med. Alla varianter hittades på rätt plats.

Vi gjorde också små (≤ 5 bp) variantutrop på Ash1 v1.1 och jämförde dessa med HG002 v4.0 benchmarkvarianterna från GIAB, som vi använde för att korrigera ett stort antal substitutions- och indelfel (se avsnittet ”Metoder”), vilket resulterade i Ash 1 v1.2. Vi återjusterade sedan Ash1-samlingen till GRCh38, kallade om varianter och jämförde dessa varianter med den nyutvecklade v4.1 benchmarkuppsättningen från GIAB. Av varianterna inom v4.1:s riktmärkesregioner matchade Ash1-varianterna 1 256 458 homozygota och 1 041 476 heterozygota SNP:er samt 187 227 homozygota och 193 524 heterozygota indels. Efter att ha uteslutit varianter inom 30 bp från en sann variant återstod 79 269 SNP och 17 439 indels, vilket (om man antar att alla dessa är fel i Ash1) motsvarar ett kvalitetsvärde (QV) på ungefär Q45 för substitutionsfel. De flesta av dessa varianter (52 191 SNP och 4629 indels) faller i segmentala duplikationer, vilket möjligen representerar saknade duplikationer i Ash1 eller ofullständig polering med korta läsningar. Sammanfattningsvis är kvaliteten på Ash1-samlingen mycket hög, med ett uppskattat substitutions-kvalitetsvärde på 62 och en indel-felfrekvens på 2 per miljon bp efter att kända segmentala duplikationer, tandemrepetitioner och homopolymerer har uteslutits.

Komparation av variant calling med hjälp av Ash1 jämfört med GRCh38

Ett av motiven för att skapa nya referensgenom är att de ger ett bättre ramverk för att analysera humana sekvensdata när man letar efter genetiska varianter som är kopplade till sjukdom. Till exempel bör en studie av ashkenaziska judar som samlat in WGS-data (whole-genome shotgun) använda ett ashkenaziskt referensgenom i stället för GRCh38. Eftersom den genetiska bakgrunden är likartad bör färre varianter hittas när man söker mot Ash1.

För att testa denna förväntan samlade vi in WGS-data från en manlig deltagare i Personal Genome Project, PGP17 (hu34D5B9). Denna individ uppskattas vara 66 % ashkenazisk enligt PGP-databasen, vilket var den högsta uppskattade fraktion som fanns tillgänglig från redan sekvenserade PGP-individer. Vi hämtade 983 220 918 100-bp reads (ungefär 33x täckning) och anpassade dem till både Ash1 och GRCh38 med Bowtie2 . En något högre andel läsningar (3 901 270, 0,5 %) anpassades till Ash1 än till GRCh38.

Vi undersökte sedan alla enskilda nukleotidvarianter (SNV, se ”Metoder”) mellan PGP17 och vart och ett av de två referensgenomerna. För att förenkla analysen tog vi endast hänsyn till platser där PGP17 var homozygot. I våra jämförelser med Ash1 identifierade vi först alla SNV:er och undersökte sedan de ursprungliga avlästa Ash1-data för att avgöra om Ash1-genomet för var och en av dessa SNV:er innehöll en annan allel som stämde överens med PGP17.

Totalt var antalet homozygota platser i PGP17 som inte stämde överens med Ash1 1 333 345, jämfört med 1 700 364 när vi jämförde homozygota platser i PGP17 med GRCh38 (Additional file 1: Table S1). Vi tittade sedan på de underliggande Ash1-avläsningsdata för de 1,33 miljoner SNV-platser som ursprungligen inte stämde överens, och fann att för ungefär hälften av dem var Ash1-genomet heterozygot, med en allel som stämde överens med PGP17. Om vi begränsade SNV:erna till platser där PGP17 och Ash1 båda är homozygota (plus ett mycket litet antal platser där Ash1 innehåller två varianter som båda skiljer sig från PGP17), minskade det totala antalet SNV:er ytterligare, till 707 756. Vi hittade alltså knappt 1 miljon färre homozygota SNV:er när vi använde Ash1 som referens för PGP17. Observera att man i stället för att anpassa till Ash1 kan man i stället anpassa läsningarna till GRCh38 och sedan ta bort kända populationsspecifika varianter. Denna tvåstegsprocess, även om den är mer komplex, kan ge liknande resultat, med undantag för regioner i Ash1 som helt enkelt saknas i GRCh38.

Bjämförelse mot vanliga ashkenaziska varianter

För att undersöka i vilken utsträckning Ash1 innehåller kända, vanliga ashkenaziska varianter (i förhållande till GRCh38) undersökte vi SNV:er som rapporterats högfrekventa i en ashkenazisk population från Genome Aggregation Database (gnomAD) . GnomAD v3.0 innehåller SNV-anrop från kortlästa helgenomdata från 1662 ashkenaziska individer. Eftersom vissa varianter endast kallades i en delmängd av dessa individer tog vi endast hänsyn till variantplatser som rapporterades hos minst 200 personer. Vi samlade sedan in SNV:er med större allel, vilket kräver att allelfrekvensen är över 0,5 i den provtagna populationen. Vi begränsade vidare vår analys till single-base-substitutioner. Detta gav oss 2 008 397 gnomAD SNV-platser där den ashkenaziska huvudallelen skilde sig från GRCh38.

Vi kunde exakt kartlägga 1 790 688 av de 2 008 397 gnomAD-platserna från GRCh38 på Ash1 (se avsnittet ”Metoder”). Vi jämförde sedan GRCh38-basen med basen för den ashkenaziska huvudallelen vid varje position, och vi undersökte också de alternativa allelerna i Ash1 på platser där den är heterozygot. För platser där avläsningsdata visade att HG002 var heterozygot och hade både Ashkenazi-majorallelen och GRCh38-allelen, bytte vi ut Ash1-basen, om det var nödvändigt, för att se till att den stämde överens med majorallelen. Efter dessa ersättningar innehöll Ash1 Ashkenazi major allelen på 88 % (1 580 866) av de 1,79 miljoner platserna. På de återstående platserna motsvarade Ash1 antingen GRCh38-allelen eftersom HG002 är homozygot för referensallelen (204 729 platser), eller så innehöll Ash1 en tredje allel som varken motsvarade GRCh38 eller gnomAD-majorallelen (5093 platser). Dessa ändringar bör ytterligare minska antalet rapporterade SNV:er när man anpassar en ashkenazisk individ till Ash1.

Värt att notera är att i takt med att frekvensen av huvudallelen i gnomAD Ashkenazipopulationen ökar, ökar också andelen platser där Ash1 matchade huvudallelen. Till exempel, av SNV:er som har en allelfrekvens > 0,9 i gnomAD Ashkenazi-populationen är över 98 % representerade i Ash1 (tabell 3).

Annotering

En viktig del av varje referensgenom är annotering: samlingen av alla gener och andra egenskaper som finns på genomet. För att Ash1 ska kunna fungera som ett verkligt referensgenom har vi kartlagt alla kända gener som används av forskarsamhället på dess koordinatsystem och använt samma gennamn och identifierare. Flera olika annoteringsdatabaser har skapats för GRCh38, och för Ash1 valde vi att använda CHESS-databasen för mänskliga gener eftersom den är heltäckande och omfattar alla proteinkodande gener i både GENCODE och RefSeq, de två andra allmänt använda gendatabaserna, och eftersom den behåller de identifierare som används i dessa kataloger. De icke-kodande generna skiljer sig åt mellan de tre databaserna, men CHESS har det största antalet genloci och isoformer. Vi använde CHESS version 2.2, som har 42 167 gener på de primära kromosomerna (exklusive GRCh38 alternativa scaffolds), varav 20 197 är proteinkodande.

Att kartlägga gener från en sammansättning till en annan är en komplicerad uppgift, särskilt för gener som förekommer i mycket likartade genfamiljer med flera kopior. Uppgiften är ännu mer komplicerad när de två sammansättningarna representerar olika individer (snarare än bara olika sammansättningar av samma individ), på grund av förekomsten av skillnader mellan enskilda nukleotider, insättningar, borttagningar, omarrangemang och genuina variationer i antalet kopior mellan individerna. Vi behövde en metod som skulle vara robust mot alla dessa potentiella skillnader.

För att lösa detta problem använde vi det nyligen utvecklade kartläggningsverktyget Liftoff, som i våra experiment var det enda verktyg som kunde kartlägga nästan alla mänskliga gener från en individ till en annan. Liftoff tar alla gener och transkript från ett genom och kartlägger dem, kromosom för kromosom, till ett annat genom. För alla gener som inte kan kartläggas till samma kromosom försöker Liftoff kartlägga dem över kromosomerna. Till skillnad från andra verktyg förlitar sig Liftoff inte på en detaljerad karta som byggs upp från en anpassning av hela genomet, utan kartlägger varje gen individuellt. Generna anpassas på transkriptnivå, inklusive introner, så att bearbetade pseudogener inte felaktigt identifieras som gener.

Vi försökte kartlägga alla 310 901 transkript från 42 167 genloci på de primära kromosomerna i GRCh38 till Ash1. Totalt lyckades vi kartlägga 309 900 (99,7 %) transkript från 42 070 genloci på huvudkromosomerna (Additional file 1: Table S2). Vi ansåg att en transkript har kartlagts framgångsrikt om mRNA-sekvensen i Ash1 är minst 50 % så lång som mRNA-sekvensen i GRCh38. Den stora majoriteten av transkript överskrider dock vida detta tröskelvärde, med 99 % av transkript som kartläggs med en täckning som är större än eller lika med 95 % (Additional file 2: Figur S2). Sekvensidentiteten hos de kartlagda transkriptionerna är lika hög, med 99 % av transkriptionerna som kartläggs med en sekvensidentitet som är större än eller lika med 94 % (Additional file 2: Figure S3).

Translokaliserade gener

Av de gener med minst en isoform som kartlagts med framgång kartlades 42 059 (99,7 %) till motsvarande platser på samma kromosom i Ash1. Av de 108 gener som inledningsvis misslyckades med att kartlägga kartlades 11 gener till en annan kromosom i 7 olika block (visas i tabell 4), vilket tyder på en translokation mellan de två genomerna. Intressant nog fanns 16 av de 22 platser som var inblandade i translokationerna i subtelomeriska regioner, vilket inträffade i 8 par där båda platserna låg nära telomerer. Detta stämmer överens med tidigare studier som rapporterar att omläggningar som involverar telomerer och subtelomerer kan vara en vanlig form av translokation hos människor .

Vi undersökte translokationen mellan kromosomerna 15 och 20, som innehåller tre av generna i tabell 4, genom att titta närmare på anpassningen mellan GRCh38 och Ash1. Translokationen ligger vid telomeren på båda kromosomerna, från position 65 079 275 till 65 109 824 (30 549 bp) på Ash1 chr20 och 101 950 338 till 101 980 928 (30 590 bp) på GRCh39 chr15. För att bekräfta translokationen anpassade vi en oberoende uppsättning mycket långa PacBio-avläsningar, alla från HG002, till Ash1 v1.7-samlingen (se avsnittet ”Metoder”) och utvärderade området runt brytpunkten på chr20. Dessa anpassningar visar en djup, konsekvent täckning som sträcker sig över många kilobaser på båda sidor av brytpunkten, vilket stödjer att Ash1-samlingen är korrekt (fig. 1).

Snapshot som visar anpassningar av långa PacBio-reads till Ash1-genomet, centrerat på den vänstra änden av den plats på kromosom 20 (position 65 079 275) där en translokation inträffade mellan kromosom 15 (GRCh38) och 20 (Ash1). Den övre delen av figuren visar koordinaterna på kr20. Nedanför visas ett histogram över lästäckningen, och de enskilda läsningarna fyller den nedre delen av figuren. Indels i avläsningarna, som visas som färgade staplar på varje avläsning, beror på den relativt höga felprocenten för de långa avläsningarna

Missade gener

Sextiotvå gener misslyckades helt med att kartlägga från GRCh38 till Ash1, och ytterligare 32 gener kartlades endast delvis (under tröskelvärdet på 50 % täckning), vilket framgår av tabell 5. Alla gener som inte lyckades kartläggas eller som kartlades delvis var medlemmar av flergenfamiljer, och i varje fall fanns det minst en annan kopia av den saknade genen i Ash1, med en genomsnittlig identitet på 98,5 %. Det finns alltså inga fall alls av en gen som finns i GRCh38 och som helt saknas i Ash1; de gener som visas i tabell 5 representerar fall där Ash1 har färre medlemmar av en multigenfamilj. Ytterligare tre gener (2 proteinkodande, 1 lncRNA) kartlades till två oplacerade contigs, vilket kommer att ge en vägledning för att placera dessa contigs i framtida utgåvor av Ash1-samlingen.

Efter att ha kartlagt generna på Ash1 extraherade vi de kodande sekvenserna från transkript som kartlades fullständigt (täckning lika med 100 %), anpassade dem till de kodande sekvenserna från GRCh38 och kallade varianter i förhållande till GRCh38 (se avsnittet ”Metoder”). Inom de 35 513 365 bp i dessa proteinkodande transkriptioner hittade vi 20 864 enskilda nukleotidvarianter och indels. Fjorton tusen fyra hundra nittionio av dessa varianter föll inom GIAB:s ”callable”-regioner för varianter med hög konfidens, även om 3963 av dessa fanns i GIAB:s ”svåra” repetitiva regioner, för vilka anpassningar ofta är tvetydiga. Av de 10 536 varianter som inte fanns i dessa svåra regioner stämde 10 456 (99,2 %) överens med GIAB:s uppsättning varianter med hög konfidens. I de svåra regionerna stämde 3804/3963 (96,0 %) överens med GIAB-uppsättningen.

Vi annoterade sedan de förändringar i aminosyror som orsakats av varianter och ofullständig kartläggning för alla proteinkodande sekvenser. Av 124 238 proteinkodande transkript från 20 197 gener har 92 600 (74,5 %) 100 % identiska proteinsekvenser. Ytterligare 26 566 (21,4 %) har minst en aminosyraförändring men ger proteiner med identisk längd, och 1561 (1,3 %) har rambevarande mutationer som infogar eller tar bort en eller flera aminosyror och lämnar resten av proteinet oförändrat. Tabell 6 visar statistik över alla förändringar i proteinsekvenser. Om ett protein hade mer än 1 variant räknade vi det under den mest följdriktiga varianten, dvs. om ett protein hade en missense-variant och en för tidig stoppkodon inkluderade vi det i gruppen ”stop gained”.

Av särskilt intresse är de transkript med varianter som påtagligt stör proteinsekvensen och som kan leda till förlust av funktion. Dessa inkluderar transkript som påverkas av en ramförskjutning (2158), stoppförlust (58), stoppförstärkning (416), startförlust (58) eller trunkering på grund av ofullständig kartläggning (564). Dessa störda isoformer representerar 885 genloci. 505 av dessa gener har dock minst en annan isoform som inte påverkas av en störande variant. Detta lämnar 380 gener där alla isoformer har minst en störning. En fullständig lista finns i Additional file 1: Table S1.