Samling og annotation af et Ashkenazi-referencegenom

For at skabe det første menneskelige referencegenom, der skal samles fra et enkelt individ, valgte vi HG002, et Ashkenazi-individ, som er en del af Personal Genome Project (PGP). PGP anvender Open Consent Model, den første ægte open-access platform til deling af individuelle menneskelige genomer, fænotyper og medicinske data . I samtykkeprocessen informeres potentielle deltagere om konsekvenserne og risiciene ved at dele genomiske data og om, hvad de kan forvente af deres deltagelse. Det åbne samtykke har gjort det muligt at skabe verdens første referencematerialer til humane genomer (HG002 er NIST-referencemateriale 8391) fra Genome In A Bottle (GIAB), som anvendes til kalibrering, udvikling af metoder til sammensætning af genomer og målinger af laboratoriers ydeevne . Alle rå sekvensdata til dette projekt blev indhentet fra GIAB, hvor de er frit tilgængelige for offentligheden .

Vi samlede HG002-genomet fra en kombination af tre dybt dækkende datasæt: 249-bp Illumina-reads, Oxford Nanopore (ONT) reads med en gennemsnitlig længde på over 33 kbp og PacBio “HiFi”-reads af høj kvalitet med en gennemsnitlig længde på 9567 bp (tabel 1).

Vi skabte oprindeligt to samlinger: en ved hjælp af Illumina- og ONT-læsninger og en anden ved hjælp af alle tre datasæt, herunder PacBio HiFi-læsninger. Tilføjelsen af PacBio HiFi-data resulterede i lidt mere samlet sekvens i samlingen (2,99 Gb vs. 2,88 Gb) og en væsentlig større contig N50-størrelse (18,2 Mb vs. 4,9 Mb). Denne samling, der blev betegnet Ash1 v0.5, var grundlaget for alle efterfølgende forbedringer.

Mapping af samlingen på kromosomer

For at skabe kromosomtilknytninger til Ash1 v0.5-samlingen brugte vi tilpasninger til GRCh38 til at mappe de fleste af stilladserne på kromosomer. De trin, der er beskrevet i afsnittet “Metoder”, genererede en række gradvist forbedrede samlinger i kromosomskala, hvilket resulterede i Ash1 v1.7. Ash1 v1.7 har større kontiguitet og mindre huller end GRCh38, som vist i tabel 2. Bemærk, at under opbygningen af disse kromosomer blev en lille mængde GRCh38-sekvenser (58,3 Mb, 2 % af genomet) anvendt til at udfylde huller i Ash1. Disse regioner omfatter nogle regioner, der er vanskelige at samle, og som er blevet kurateret manuelt for GRCh38. I alt er den anslåede størrelse af alle huller i Ash1 82,9 Mbp, sammenlignet med 84,7 Mbp i GRCh38.p13.

Som en del af forbedringen af samlingsprocessen søgte vi i en af de foreløbige Ash1-samlinger (v1.1) efter de 12.745 isolerede strukturelle varianter af høj kvalitet (indsættelser og sletninger ≥ 50 bp), som Zook et al. identificerede ved at sammenligne Ashkenazi-triodataene med GRCh37 . I denne undersøgelse blev der anvendt fire forskellige sekventeringsteknologier og flere variant callers til at identificere varianter og filtrere falske positive varianter fra. Af disse 12 745 SV’er er 5807 homozygote og 6938 heterozygote. Vi forventede, at Ash1-samlingen ville stemme overens med næsten alle de homozygote varianter. Da Ash1 kun indfanger én haplotype, forventede vi, at den ville være enig med ca. halvdelen af de heterozygote SV’er, idet vi antog, at samlealgoritmen valgte tilfældigt mellem haplotyperne, da den besluttede, hvilken variant der skulle indgå i den endelige konsensus. Af de 5807 homozygote varianter var 5284 (91 %) til stede ved hjælp af vores matchkriterier (se afsnittet “Metoder”), og 3922 (56,5 %) af 6938 heterozygote varianter var til stede. Alle varianter blev fundet på den korrekte placering.

Vi foretog også små (≤ 5 bp) variantopkald på Ash1 v1.1 og sammenlignede disse med HG002 v4.0-benchmarkvarianterne fra GIAB, som vi brugte til at korrigere talrige substitutions- og indelfejl (se afsnittet “Metoder”), hvilket gav Ash 1 v1.2. Vi justerede derefter Ash1-samlingen på ny til GRCh38, genkaldte varianter og sammenlignede disse varianter med det nyudviklede v4.1 GIAB-benchmark-sæt fra GIAB. Af varianterne inden for v4.1 benchmark-regionerne matchede Ash1-varianterne 1 256 458 homozygote og 1 041 476 heterozygote SNP’er og 187 227 homozygote og 193 524 heterozygote indels. Efter udelukkelse af variantopkald inden for 30 bp fra en ægte variant var der 79 269 SNP’er og 17 439 indels tilbage, hvilket (hvis man antager, at disse er alle fejl i Ash1) svarer til en kvalitetsværdi (QV) på ca. Q45 for substitutionsfejl. De fleste af disse varianter (52 191 SNP’er og 4629 indels) falder i segmentale duplikationer, hvilket muligvis repræsenterer manglende duplikationer i Ash1 eller ufuldstændig polering med korte læsninger. Sammenfattende er kvaliteten af Ash1-samlingen meget høj med en estimeret substitutions-kvalitetsværdi på 62 og en indel-fejlrate på 2 pr. million bp efter udelukkelse af kendte segmentale duplikationer, tandemrepeats og homopolymerer.

Sammenligning af variant calling ved hjælp af Ash1 versus GRCh38

En af motiverne for at skabe nye referencegenomer er, at de giver en bedre ramme for analyse af menneskelige sekvensdata, når der søges efter genetiske varianter, der er knyttet til sygdom. For eksempel bør en undersøgelse af ashkenazi-jøder, der indsamlede whole-genome shotgun-data (WGS), bruge et ashkenazi-referencegenom i stedet for GRCh38. Da den genetiske baggrund er ens, bør der findes færre varianter, når der søges mod Ash1.

For at teste denne forventning indsamlede vi WGS-data fra en mandlig deltager i Personal Genome Project, PGP17 (hu34D5B9). Denne person skønnes at være 66 % ashkenazisk i henhold til PGP-databasen, hvilket var den højeste estimerede fraktion, der var tilgængelig fra allerede sekventerede PGP-individer. Vi downloadede 983.220.918.918.100-bp-reads (ca. 33x dækning) og justerede dem til både Ash1 og GRCh38 ved hjælp af Bowtie2 . En lidt højere fraktion af læsninger (3.901.270, 0,5 %) blev justeret til Ash1 end til GRCh38.

Vi undersøgte derefter alle enkeltnukleotidvarianter (SNV’er, se “Metoder”) mellem PGP17 og hvert af de to referencegenomer. For at forenkle analysen undersøgte vi kun steder, hvor PGP17 var homozygot. I vores sammenligninger med Ash1 identificerede vi først alle SNV’er og undersøgte derefter de originale Ash1-læsedata for at afgøre, om Ash1-genomet for hver af disse SNV’er indeholdt en anden allel, der matchede PGP17.

I alt var antallet af homozygote steder i PGP17, der ikke stemte overens med Ash1, 1.333.345, mod 1.700.364, da vi sammenlignede homozygote steder i PGP17 med GRCh38 (Additional file 1: Table S1). Vi kiggede derefter på de underliggende Ash1-læsedata for de 1,33 millioner SNV-steder, der oprindeligt ikke matchede, og fandt, at for ca. halvdelen af dem var Ash1-genomet heterozygot, med en allel, der matchede PGP17. Hvis vi begrænsede SNV’erne til steder, hvor PGP17 og Ash1 begge er homozygote (plus et meget lille antal steder, hvor Ash1 indeholder to varianter, der begge adskiller sig fra PGP17), reducerede dette det samlede antal SNV’er yderligere til 707.756. Vi fandt således lige under 1 million færre homozygote SNV’er, når vi brugte Ash1 som reference for PGP17. Bemærk, at man i stedet for at tilpasse til Ash1 i stedet kunne tilpasse læsningerne til GRCh38 og derefter fjerne kendte populationsspecifikke varianter. Denne to-trins proces, selv om den er mere kompleks, kan give lignende resultater, bortset fra regioner i Ash1, der simpelthen mangler i GRCh38.

Sammenligning med almindelige Ashkenazi-varianter

For at undersøge, i hvilket omfang Ash1 indeholder kendte, almindelige Ashkenazi-varianter (i forhold til GRCh38), undersøgte vi SNV’er, der er rapporteret med høj frekvens i en Ashkenazi-population fra Genome Aggregation Database (gnomAD) . GnomAD v3.0 indeholder SNV-opkald fra short-read helgenomdata fra 1662 Ashkenazi-individer. Da nogle varianter kun blev kaldt i en delmængde af disse individer, overvejede vi kun variantsteder, der blev rapporteret hos mindst 200 personer. Derefter indsamlede vi større allel-SNV’er, idet vi krævede, at allelfrekvensen skulle være over 0,5 i den udtagne population. Vi begrænsede yderligere vores analyse til enkeltbasesubstitutioner. Dette gav os 2 008 397 gnomAD SNV-steder, hvor den ashkenaziske major allel adskilte sig fra GRCh38.

Vi var i stand til præcist at kortlægge 1 790 688 af de 2 008 397 gnomAD-steder fra GRCh38 på Ash1 (se afsnittet “Metoder”). Vi sammenlignede derefter GRCh38-basen med den ashkenaziske hovedallelbase på hver position, og vi undersøgte også de alternative alleler i Ash1 på steder, hvor den er heterozygot. For steder, hvor læsedataene viste, at HG002 var heterozygot og havde både Ashkenazi-majorallelen og GRCh38-allelen, udskiftede vi Ash1-basen, hvis det var nødvendigt, for at sikre, at den stemte overens med majorallelen. Efter disse udskiftninger indeholdt Ash1 den ashkenaziske majorallel på 88 % (1 580 866) af de 1,79 mio. steder. På de resterende steder matchede Ash1 enten GRCh38-allelen, fordi HG002 er homozygot for referenceallelen (204 729 steder), eller også indeholdt den en tredje allel, der hverken matchede GRCh38 eller gnomAD-majorallelen (5093 steder). Disse ændringer bør yderligere reducere antallet af rapporterede SNV’er, når et ashkenazisk individ tilpasses til Ash1.

Værd at bemærke er, at efterhånden som hyppigheden af major-allelen i gnomAD Ashkenazi-populationen stiger, stiger også andelen af steder, hvor Ash1 matchede major-allelen. For eksempel er over 98 % af de SNV’er, der har en allelfrekvens > 0,9 i gnomAD Ashkenazi-populationen, repræsenteret i Ash1 (tabel 3).

Annotation

En vigtig del af ethvert referencegenom er annotation: en samling af alle gener og andre træk, der findes på genomet. For at Ash1 kan fungere som et ægte referencegenom, har vi kortlagt alle de kendte gener, der anvendes af det videnskabelige samfund, på dets koordinatsystem, idet vi har anvendt de samme gennavne og identifikatorer. Der er blevet oprettet flere forskellige annotationsdatabaser for GRCh38, og for Ash1 valgte vi at bruge CHESS-databasen for menneskelige gener, fordi den er omfattende og omfatter alle de proteinkodende gener i både GENCODE og RefSeq, de to andre almindeligt anvendte gendatabaser, og fordi den bevarer de identifikatorer, der anvendes i disse kataloger. De ikke-kodende gener er forskellige i de tre databaser, men CHESS har det største antal genloci og isoformer. Vi brugte CHESS version 2.2, som har 42.167 gener på de primære kromosomer (eksklusive GRCh38 alternative stilladser), hvoraf 20.197 er proteinkodende.

Mapping af gener fra en samling til en anden er en kompleks opgave, især for gener, der forekommer i meget ensartede, multikopierede genfamilier. Opgaven er endnu mere kompleks, når de to samlinger repræsenterer forskellige individer (snarere end blot forskellige samlinger af det samme individ), på grund af tilstedeværelsen af enkelt-nucleotidforskelle, insertioner, deletioner, rearrangementer og ægte variationer i antallet af kopier mellem individerne. Vi havde brug for en metode, der ville være robust over for alle disse potentielle forskelle.

For at løse dette problem brugte vi det nyligt udviklede Liftoff-kortlægningsværktøj, som i vores eksperimenter var det eneste værktøj, der kunne kortlægge næsten alle menneskelige gener fra et individ til et andet. Liftoff tager alle gener og transskriptioner fra et genom og kortlægger dem, kromosom for kromosom, til et andet genom. For alle gener, som ikke kan kortlægges på det samme kromosom, forsøger Liftoff at kortlægge dem på tværs af kromosomer. I modsætning til andre værktøjer er det ikke afhængig af et detaljeret kort, der er opbygget ud fra en alignment af hele genomet, men i stedet kortlægger det hvert enkelt gen individuelt. Generne tilpasses på transkriptniveau, herunder introner, så forarbejdede pseudogener ikke fejlagtigt identificeres som gener.

Vi forsøgte at kortlægge alle 310.901 transkripter fra 42.167 genloci på de primære kromosomer i GRCh38 til Ash1. I alt lykkedes det os at kortlægge 309 900 (99,7 %) transkripter fra 42 070 genloci på hovedkromosomerne (Yderligere fil 1: Tabel S2). Vi anså et transkript for at være kortlagt med succes, hvis mRNA-sekvensen i Ash1 er mindst 50 % så lang som mRNA-sekvensen på GRCh38. Langt størstedelen af transskriptionerne overskrider imidlertid i høj grad denne tærskel, idet 99 % af transskriptionerne kortlægges med en dækning større end eller lig med 95 % (Yderligere fil 2: Figur S2). Sekvensidentiteten af de kortlagte transkripter er ligeledes høj, idet 99 % af transkripterne kortlægges med en sekvensidentitet større end eller lig med 94 % (Yderligere fil 2: Figur S3).

Translokaliserede gener

Af de gener med mindst én succesfuldt kortlagt isoform kortlægges 42.059 (99,7 %) til de tilsvarende steder på det samme kromosom i Ash1. Af de 108 gener, der oprindeligt ikke kunne kortlægges, blev 11 gener kortlagt til et andet kromosom i 7 forskellige blokke (vist i tabel 4), hvilket tyder på en translokation mellem de to genomer. Det er interessant, at 16 af de 22 steder, der var involveret i translokationerne, lå i subtelomeriske regioner, hvilket skete i 8 par, hvor begge steder var nær telomerer. Dette er i overensstemmelse med tidligere undersøgelser, der rapporterer, at rearrangementer, der involverer telomerer og subtelomerer, kan være en almindelig form for translokation hos mennesker .

Vi undersøgte translokationen mellem kromosom 15 og 20, som indeholder tre af generne i tabel 4, ved at se nærmere på tilpasningen mellem GRCh38 og Ash1. Translokationen befinder sig ved telomeren på begge kromosomer, fra position 65.079.275 til 65.109.824 (30.549 bp) på Ash1 chr20 og 101.950.338 til 101.980.928 (30.590 bp) på GRCh39 chr15. For at bekræfte translokationen tilpassede vi et uafhængigt sæt meget lange PacBio-reads, alle fra HG002, til Ash1 v1.7-samlingen (se afsnittet “Metoder”) og evaluerede området omkring brudpunktet på chr20. Disse tilpasninger viser en dyb, konsistent dækning, der strækker sig over mange kilobaser på begge sider af brudpunktet, hvilket understøtter, at Ash1-samlingen er korrekt (Fig. 1).

Snapshot, der viser alignments af lange PacBio-reads til Ash1-genomet, centreret på den venstre ende af det sted på kromosom 20 (position 65,079,275), hvor der skete en translokation mellem kromosom 15 (GRCh38) og 20 (Ash1). Den øverste del af figuren viser koordinaterne på chr20. Herunder ses et histogram af læsedækningen, og de enkelte læsninger fylder den nederste del af figuren. Indels i læsningerne, der er vist som farvede bjælker på hver læsning, skyldes den relativt høje fejlfrekvens for de lange læsninger

Missing genees

Det lykkedes ikke at kortlægge 62 gener helt fra GRCh38 til Ash1, og yderligere 32 gener kortlagde kun delvist (under 50 %-dækningstærsklen), som det fremgår af tabel 5. Alle de gener, der ikke kunne kortlægges eller som kortlagdes delvist, var medlemmer af flergenfamilier, og i alle tilfælde var der mindst én anden kopi af det manglende gen til stede i Ash1 med en gennemsnitlig identitet på 98,5 %. Der er således overhovedet ingen tilfælde af et gen, der er til stede i GRCh38, og som er helt fraværende i Ash1; de gener, der er vist i tabel 5, repræsenterer tilfælde, hvor Ash1 har færre medlemmer af en flergenfamilie. Tre yderligere gener (2 proteinkoderende, 1 lncRNA) er kortlagt til to ikke-placerede contigs, hvilket vil give en vejledning til placering af disse contigs i fremtidige udgivelser af Ash1-samlingen.

Når vi havde kortlagt generne på Ash1, ekstraherede vi de kodende sekvenser fra transkriptioner, der blev kortlagt fuldstændigt (dækning lig med 100 %), tilpassede dem til de kodende sekvenser fra GRCh38 og kaldte varianter i forhold til GRCh38 (se afsnittet “Metoder”). Inden for de 35.513.365 bp i disse proteinkodende transskriptioner fandt vi 20.864 enkelt-nucleotidvarianter og indels. Fjorten tusinde fire hundrede nioghalvfems af disse varianter faldt inden for GIAB’s “callable”-regioner for varianter med høj tillid, selv om 3963 af disse var i GIAB’s “vanskelige” gentagne regioner, for hvilke tilpasninger ofte er tvetydige. Af de 10 536 varianter, der ikke befandt sig i disse vanskelige regioner, stemte 10 456 (99,2 %) overens med GIAB’s sæt af varianter med høj konfidens. I de vanskelige regioner stemte 3804/3963 (96,0 %) overens med GIAB-sættet.

Vi annoterede derefter ændringerne i aminosyrer forårsaget af varianter og ufuldstændig kortlægning for alle protein-kodende sekvenser. Ud af 124.238 proteinkodende transkript fra 20.197 gener har 92.600 (74,5 %) 100 % identiske proteinsekvenser. Andre 26 566 (21,4 %) har mindst én aminosyreændring, men giver proteiner med identisk længde, og 1561 (1,3 %) har rammebevarende mutationer, der indsætter eller sletter en eller flere aminosyrer, men lader resten af proteinet forblive uændret. Tabel 6 viser statistikker over alle ændringer i proteinsekvenser. Hvis et protein havde mere end 1 variant, talte vi det under den mest konsekvensfulde variant, dvs. hvis et protein havde en missense-variant og et for tidligt stopkodon, inkluderede vi det i gruppen “stop gained”.

Af særlig interesse er de transkripter med varianter, der i væsentlig grad forstyrrer proteinsekvensen og kan resultere i tab af funktion. Disse omfatter transkriptioner, der er påvirket af et frameshift (2158), stoptab (58), stopforøgelse (416), starttab (58) eller trunkering på grund af ufuldstændig kortlægning (564). Disse forstyrrede isoformer repræsenterer 885 genloci; 505 af disse gener har dog mindst én anden isoform, som ikke er påvirket af en forstyrrende variant. Dette efterlader 380 gener, hvor alle isoformer har mindst én afbrydelse; den fuldstændige liste findes i Additional file 1: Table S1.