Askenazi-ihmisen referenssigenomin kokoaminen ja annotaatio

Valitsimme HG002:n, joka on henkilökohtaisen genomiprojektin (Personal Genome Project, PGP) osana oleva askenazi-yksilö, luodaksemme ensimmäisen ihmisen referenssigenomin, joka kootaan yhdestä yksilöstä. PGP käyttää Open Consent Model -mallia, joka on ensimmäinen aidosti avoin foorumi yksittäisten ihmisten genomin, fenotyypin ja lääketieteellisten tietojen jakamiseen. Suostumusprosessi valistaa mahdollisia osallistujia genomitietojen jakamisen vaikutuksista ja riskeistä sekä siitä, mitä he voivat odottaa osallistumiseltaan. Avoimen suostumuksen ansiosta on voitu luoda maailman ensimmäiset ihmisen genomin vertailumateriaalit (HG002 on NIST Reference Material 8391) Genome In A Bottle (GIAB) -laitteesta, jota käytetään kalibrointiin, genomin kokoamismenetelmien kehittämiseen ja laboratorion suorituskyvyn mittaamiseen . Kaikki tämän hankkeen raakasekvenssidata saatiin GIAB:stä, jossa se on vapaasti yleisön saatavilla .

Kokosimme HG002:n genomin kolmen syvän kattavuuden datasarjan yhdistelmästä: 249-bp:n Illumina-lukemat, Oxford Nanopore (ONT) -lukemat, joiden pituus on keskimäärin yli 33 kbp, ja korkealaatuiset PacBio ”HiFi” -lukemat, joiden pituus on keskimäärin 9567 bp (taulukko 1).

Loimme aluksi kaksi kokoelmaa: toinen käyttäen Illumina- ja ONT-lukemia ja toinen käyttäen kaikkia kolmea datasarjaa, mukaan lukien PacBio HiFi -lukemat. PacBio HiFi -datan lisääminen johti hieman suurempaan kokonaissekvenssin määrään assemblaatiossa (2,99 Gb vs. 2,88 Gb) ja huomattavasti suurempaan contigin N50-kokoon (18,2 Mb vs. 4,9 Mb). Tämä kokoonpano, jota nimitettiin Ash1 v0.5:ksi, oli perustana kaikille myöhemmille tarkennuksille.

Kokoonpanon kartoittaminen kromosomeihin

Luodaksemme Ash1 v0.5 -kokoonpanon kromosomitunnukset käytimme GRCh38:aan tehtyjä linjauksia kartoittaaksemme suurimman osan telineistä kromosomeihin. Kohdassa ”Menetelmät” kuvatut vaiheet tuottivat sarjan asteittain parannettuja kromosomimittakaavan kokoonpanoja, joiden tuloksena syntyi Ash1 v1.7. Ash1 v1.7:ssä on suurempi yhtenäisyys ja pienemmät aukot kuin GRCh38:ssa, kuten taulukosta 2 käy ilmi. On huomattava, että näitä kromosomeja rakennettaessa käytettiin pieni määrä GRCh38-sekvenssiä (58,3 Mb, 2 % genomista) Ash1:n aukkojen täyttämiseen. Näihin alueisiin sisältyy joitakin vaikeasti koottavia alueita, jotka on kuratoitu manuaalisesti GRCh38:lle. Kaikkien aukkojen arvioitu koko Ash1:ssä on yhteensä 82,9 Mbp, kun se GRCh38:ssa on 84,7 Mbp.p13.

Osana kokoonpanon parantamisprosessia etsimme yhdestä alustavasta Ash1-kokoonpanosta (v1.1) 12 745 laadukasta, eristettyä rakenteellista varianttia (insertionit ja deletionit ≥ 50 bp:n pituisia), jotka Zook ym. havaitsivat vertailemalla Ashkenazi-kolmioaineistoa GRCh:hen37 . Kyseisessä tutkimuksessa käytettiin neljää eri sekvensointitekniikkaa ja useita varianttisoittimia varianttien tunnistamiseen ja väärien positiivisten tulosten suodattamiseen. Näistä 12 745 SV:stä 5807 on homotsygoottisia ja 6938 heterotsygoottisia. Odotimme, että Ash1-kokoonpano vastaa lähes kaikkia homotsygoottisia variantteja. Koska Ash1 kattaa vain yhden haplotyypin, odotimme, että se olisi sopusoinnussa noin puolen heterotsygoottisten SV:iden kanssa, olettaen, että assosiaatioalgoritmi valitsee satunnaisesti haplotyyppien välillä päättäessään, mikä variantti sisällytetään lopulliseen konsensukseen. 5807 homotsygoottisesta variantista 5284 (91 %) esiintyi sovituskriteereillämme (ks. kohta ”Menetelmät”), ja 6938 heterotsygoottisesta variantista 3922 (56,5 %) esiintyi. Kaikki variantit löytyivät oikeasta sijainnista.

Toimme myös pieniä (≤ 5 bp) varianttipyyntöjä Ash1 v1.1:llä ja vertasimme niitä GIAB:n HG002 v4.0 -vertailuvariantteihin, joita käytimme lukuisten substituutio- ja indel-virheiden korjaamiseen (ks. ”Menetelmät”-osio), minkä tuloksena saimme tulokseksi Ash 1 v1.2:n. Tämän jälkeen sovitimme Ash1-kokoonpanon uudelleen GRCh38:aan, kutsuimme variantit uudelleen ja vertasimme näitä variantteja vastikään kehitettyyn v4.1 GIAB:n vertailujoukkoon. Vertailualueiden v4.1 sisällä olevista varianteista Ash1-variantit vastasivat 1 256 458 homotsygoottista ja 1 041 476 heterotsygoottista SNP:tä sekä 187 227 homotsygoottista ja 193 524 heterotsygoottista indeliä. Kun 30 bp:n etäisyydellä oikeasta variantista olevat varianttihuudot jätettiin pois, jäljelle jäi 79 269 SNP:tä ja 17 439 indeliä, mikä (olettaen, että nämä ovat kaikki Ash1:ssä esiintyviä virheitä) vastaa substituutiovirheiden laatuarvoa (QV), joka on noin Q45. Suurin osa näistä varianteista (52 191 SNP:tä ja 4629 indeliä) kuuluu segmentaalisiin duplikaatioihin, jotka mahdollisesti edustavat Ash1:ssä puuttuvia duplikaatioita tai lyhyiden lukukertojen epätäydellistä kiillotusta. Yhteenvetona voidaan todeta, että Ash1-kokoonpanon laatu on erittäin korkea, sillä sen arvioitu substituution laatuarvo on 62 ja indel-virheiden määrä 2 miljoonaa bp kohti sen jälkeen, kun tunnetut segmentaaliset duplikaatiot, tandemtoistot ja homopolymeerit on jätetty pois.

Vertailu varianttien kutsumisesta Ash1:n ja GRCh38:n välillä

Yksi motiiveista uusien referenssigenomien luomiselle on se, että ne tarjoavat paremmat puitteet ihmisen sekvenssidatan analysoinnille etsiessämme geneettisiä varianttimuunnoksia, jotka ovat yhteydessä sairauksiin. Esimerkiksi askenasijuutalaisia koskevassa tutkimuksessa, jossa kerättiin koko genomin shotgun (WGS) -dataa, olisi käytettävä askenasijuutalaisten referenssigenomia GRCh38:n sijasta. Koska geneettinen tausta on samankaltainen, Ash1:stä etsittäessä pitäisi löytyä vähemmän variantteja.

Testataksemme tätä odotusta keräsimme WGS-dataa Personal Genome Project -hankkeen miespuoliselta osallistujalta, PGP17:ltä (hu34D5B9). Tämän yksilön arvioidaan olevan PGP-tietokannan mukaan 66-prosenttisesti askenasialainen, mikä oli korkein arvioitu osuus, joka oli saatavilla jo sekvensoiduista PGP-yksilöistä. Latasimme 983 220 918 100-bp:n lukemat (noin 33-kertainen kattavuus) ja kohdistimme ne Bowtie2 -ohjelmalla sekä Ash1:n että GRCh38:n kanssa. Hieman suurempi osa lukemista (3 901 270, 0,5 %) kohdistui Ash1:een kuin GRCh38:aan.

Tarkastelimme tämän jälkeen kaikki yksittäisten nukleotidivarianttien (SNV:t, ks. ”Menetelmät”) väliset erot PGP17:n ja kummankin referenssigenomin välillä. Analyysin yksinkertaistamiseksi otimme huomioon vain paikat, joissa PGP17 oli homotsygoottinen. Vertailuissamme Ash1:een tunnistimme ensin kaikki SNV:t ja tarkastelimme sitten alkuperäistä Ash1:n lukutietoa määrittääksemme, oliko Ash1-genomissa kunkin tällaisen SNV:n osalta eri alleeli, joka vastasi PGP17:ää.

Kokonaisuudessaan niiden homotsygoottisten kohtien määrä PGP17:ssä, jotka olivat erimielisiä Ash1:n kanssa, oli 1333345, kun taas niiden määrä oli 1700364 verrattaessa homotsygoottisia kohtia PGP17:ssä GRCh38:aan (Lisätiedosto 1: Taulukko S1). Tämän jälkeen tarkastelimme taustalla olevia Ash1-lukutietoja niiden 1,33 miljoonan SNV-kohdan osalta, jotka alun perin eivät vastanneet toisiaan, ja havaitsimme, että noin puolessa niistä Ash1-genomi oli heterotsygoottinen, ja yksi alleeli vastasi PGP17:ää. Jos rajoitimme SNV:t kohtiin, joissa PGP17 ja Ash1 ovat molemmat homotsygoottisia (sekä hyvin pieneen määrään kohtia, joissa Ash1:ssä on kaksi varianttia, jotka molemmat eroavat PGP17:stä), tämä vähensi SNV:iden kokonaismäärää entisestään 707 756:een. Löysimme siis hieman alle miljoona homotsygoottista SNV:tä vähemmän, kun käytimme Ash1:tä PGP17:n vertailukohtana. Huomattakoon, että sen sijaan, että lukemat olisi kohdistettu Ash1:een, ne voitaisiin sen sijaan kohdistaa GRCh38:aan ja sitten poistaa tunnetut populaatiokohtaiset variantit. Tämä kaksivaiheinen prosessi, vaikkakin monimutkaisempi, saattaa tuottaa samanlaisia tuloksia, lukuun ottamatta Ash1:n alueita, jotka yksinkertaisesti puuttuvat GRCh38:sta.

Vertailu yleisiin askenasivariantteihin

Tarkastellaksemme, missä määrin Ash1 sisältää tunnettuja, yleisiä askenasivariantteja (suhteessa GRCh38:aan), tarkastelimme SNV:itä, jotka on raportoitu korkealla frekvenssillä askenasipopulaatiossa Genome Aggregation Database -tietokannasta (gnomAD) . GnomAD v3.0 sisältää SNV-kutsuja 1662 ashkenazi-yksilön lyhyen lukukerran kokogenomitiedoista. Koska joitakin variantteja kutsuttiin vain osajoukossa näistä henkilöistä, otimme huomioon vain varianttipaikat, jotka oli ilmoitettu vähintään 200 henkilöllä. Tämän jälkeen keräsimme suurten alleelien SNV:t, ja edellytimme, että alleelien frekvenssi on yli 0,5 näytteeksi otetussa populaatiossa. Rajoitimme analyysimme lisäksi yhden emäksen substituutioihin. Näin saimme 2 008 397 gnomAD SNV-kohtaa, joissa askenasien pääalleeli erosi GRCh38:sta.

Pystyimme kartoittamaan täsmällisesti 1 790 688 GRCh38:n 2 008 397 gnomAD-kohdasta 1 790 688 kohdetta Ash1:een (ks. kohta ”Menetelmät”). Tämän jälkeen vertasimme GRCh38:n emästä Ashkenazin pääalleelin emäkseen kussakin paikassa, ja tutkimme myös Ash1:n vaihtoehtoisia alleeleja paikoissa, joissa se on heterotsygoottinen. Niissä kohdissa, joissa lukutiedot osoittivat, että HG002 oli heterotsygoottinen ja että sillä oli sekä Ashkenazi major -alleeli että GRCh38-alleeli, korvasimme tarvittaessa Ash1:n emäksen varmistaaksemme, että se vastasi major-alleelia. Näiden korvausten jälkeen Ash1 sisälsi Ashkenazi major -alleelin 88 prosentissa (1 580 866) 1,79 miljoonasta kohdasta. Jäljelle jääneissä paikoissa Ash1 joko vastasi GRCh38-alleelia, koska HG002 on homotsygootti referenssialleelin suhteen (204 729 paikkaa), tai se sisälsi kolmannen alleelin, joka ei vastannut GRCh38- eikä gnomAD-pääalleelia (5093 paikkaa). Näiden muutosten pitäisi edelleen vähentää raportoitujen SNV:iden määrää, kun askenasi-yksilöä kohdistetaan Ash1:een.

Huomionarvoista on, että kun pääalleelin frekvenssi gnomAD-askenasipopulaatiossa kasvaa, kasvaa myös niiden paikkojen osuus, joissa Ash1 vastasi pääalleelia. Esimerkiksi niistä SNV:istä, joiden alleelifrekvenssi on > 0,9 gnomAD-askenasipopulaatiossa, yli 98 % on edustettuna Ash1:ssä (taulukko 3).

Annotaatio

Tärkeä osa mitä tahansa referenssigenomia on annotaatio: kokoelma kaikista genomista löytyvistä geeneistä ja muista ominaisuuksista. Jotta Ash1 voisi toimia todellisena referenssigenomina, olemme kartoittaneet kaikki tiedeyhteisön käyttämät tunnetut geenit sen koordinaatistoon käyttäen samoja geenien nimiä ja tunnuksia. GRCh38:lle on luotu useita erilaisia annotaatiotietokantoja, ja Ash1:n osalta päätimme käyttää CHESS-ihmisen geenitietokantaa, koska se on kattava ja sisältää kaikki proteiineja koodaavat geenit sekä GENCODE- että RefSeq-tietokannoista, kahdesta muusta laajalti käytetystä geenitietokannasta, ja koska siinä on säilytetty mainituissa tietokannoissa käytetyt tunnisteet. Koodaamattomat geenit eroavat toisistaan kolmessa tietokannassa, mutta CHESS sisältää eniten geenilokuksia ja isomuotoja. Käytimme CHESSin versiota 2.2, jossa on 42 167 geeniä primaarikromosomeissa (pois lukien GRCh38-vaihtoehtoiset telineet), joista 20 197 on proteiineja koodaavia.

Geenien kartoittaminen yhdestä kokoonpanosta toiseen on monimutkainen tehtävä erityisesti sellaisten geenien kohdalla, jotka esiintyvät hyvin samankaltaisissa, monikopioisissa geeniperheissä. Tehtävä on vieläkin monimutkaisempi, kun kaksi assemblaatiota edustavat eri yksilöitä (eikä vain saman yksilön eri assemblaatioita), koska yksilöiden välillä on yhden nukleotidin eroja, insertioita, deletioita, uudelleenjärjestelyjä ja todellisia kopioluvun vaihteluita. Tarvitsimme menetelmän, joka olisi kestävä kaikkien näiden potentiaalisten erojen suhteen.

Tämän ongelman ratkaisemiseksi käytimme hiljattain kehitettyä Liftoff-kartoitustyökalua, joka oli kokeissamme ainoa työkalu, jolla pystyttiin kartoittamaan lähes kaikki ihmisen geenit yksilöstä toiseen. Liftoff ottaa kaikki geenit ja transkriptit genomista ja kartoittaa ne, kromosomi kromosomilta, toiseen genomiin. Kaikkia geenejä, joita ei pystytä kartoittamaan samaan kromosomiin, Liftoff yrittää kartoittaa eri kromosomien välillä. Toisin kuin muut työkalut, se ei luota yksityiskohtaiseen karttaan, joka on rakennettu koko genomin kohdistuksesta, vaan se kartoittaa jokaisen geenin erikseen. Geenit kohdistetaan transkriptiotasolla, mukaan lukien intronit, jotta prosessoituja pseudogeenejä ei erehdyksessä tunnisteta geeneiksi.

Yritimme kartoittaa kaikki 310 901 transkriptiä 42 167 geenilokuksen 42 167 geenilokista GRCh38:n primaarikromosomeista Ash1:een. Kaiken kaikkiaan onnistuimme kartoittamaan 309 900 (99,7 %) transkriptiä 42 070 geenilokista pääkromosomeihin (lisätiedosto 1: taulukko S2). Katsoimme, että transkripti on onnistuneesti kartoitettu, jos mRNA-sekvenssi Ash1:ssä on vähintään 50 prosenttia yhtä pitkä kuin mRNA-sekvenssi GRCh38:ssa. Valtaosa transkripteistä kuitenkin ylitti tämän kynnysarvon huomattavasti, sillä 99 prosenttia transkripteistä kartoitettiin vähintään 95 prosentin kattavuudella (lisätiedosto 2: kuva S2). Kartoitettujen transkriptien sekvenssi-identiteetti on yhtä korkea, sillä 99 % transkripteistä kartoitettiin vähintään 94 %:n sekvenssi-identiteetillä (Additional file 2: Figure S3).

Transkriptoidut geenit

Niistä geeneistä, joilla oli vähintään yksi onnistuneesti kartoitettu isoformi, 42 059 (99,7 %) kartoitettiin vastaaville paikoille samassa kromosomissa Ash1:ssä. Niistä 108 geenistä, joita ei alun perin onnistuttu kartoittamaan, 11 geeniä karttui eri kromosomiin seitsemässä eri lohkossa (taulukossa 4), mikä viittaa translokaatioon näiden kahden genomin välillä. Mielenkiintoista on, että 22:sta translokaatioon osallistuneesta paikasta 16 sijaitsi subtelomeerisilla alueilla, jotka esiintyivät kahdeksassa parissa, joissa molemmat sijainnit olivat lähellä telomeerejä. Tämä on yhdenmukaista aiempien tutkimusten kanssa, joissa on raportoitu, että telomeereihin ja subtelomeereihin liittyvät uudelleenjärjestelyt voivat olla yleinen translokaation muoto ihmisillä.

Tarkastelimme kromosomien 15 ja 20 välistä translokaatiota, joka sisältää kolme taulukossa 4 esitetyistä geeneistä, tarkastelemalla tarkemmin GRCh38:n ja Ash1:n välistä linjausta. Translokaatio on molempien kromosomien telomeerissa, Ash1:n chr20:n kohdasta 65 079 275 – 65 109 824 (30 549 bp) ja GRCh39:n chr15:n kohdasta 101 950 338 – 101 980 928 (30 590 bp). Translokaation vahvistamiseksi linjasimme riippumattoman joukon erittäin pitkiä PacBio-lukemia, kaikki HG002:sta, Ash1 v1.7 -kokoonpanoon (ks. jakso ”Menetelmät”) ja arvioimme katkaisukohdan ympärillä olevan alueen chr20:ssä. Nämä linjaukset osoittavat syvää, johdonmukaista kattavuutta, joka ulottuu useita kilobaaseja katkeamiskohdan molemmin puolin, mikä tukee Ash1-kokoonpanon oikeellisuutta (kuva 1).

Snapshot, jossa näkyy pitkien PacBio-lukujen kohdistukset Ash1-genomiin, keskitetty kromosomissa 20 olevan paikan vasempaan päähän (sijainti 65,079,275), jossa tapahtui translokaatio kromosomien 15 (GRCh38) ja 20 (Ash1) välillä. Kuvan yläosassa on koordinaatit chr20:ssä. Sen alapuolella on lukupeiton histogrammi, ja yksittäiset lukemat täyttävät kuvan alaosan. Lukemissa olevat indelit, jotka näkyvät värillisinä palkkeina kussakin lukemassa, johtuvat pitkien lukemien suhteellisen suuresta virheprosentista

Puuttuvat geenit

Kuudellakymmenelläkahdella geenillä ei onnistuttu kartoittamaan kokonaan GRCh38:sta Ash1:een, ja toiset 32 geeniä kartoitettiin vain osittain (alle 50 prosentin kattavuuden kynnyksen), kuten taulukosta 5 näkyy. Kaikki geenit, jotka eivät kartoittuneet tai jotka kartoittuivat osittain, kuuluivat monigeeniperheisiin, ja jokaisessa tapauksessa Ash1:ssä oli vähintään yksi toinen kopio puuttuvasta geenistä, jonka keskimääräinen identiteetti oli 98,5 prosenttia. Näin ollen ei ole yhtään tapausta, jossa GRCh38:ssa esiintyvä geeni puuttuisi kokonaan Ash1:stä; taulukossa 5 esitetyt geenit edustavat tapauksia, joissa Ash1:ssä on vähemmän monigeeniperheen jäseniä. Kolme ylimääräistä geeniä (2 proteiinia koodaavaa, 1 lncRNA) karttui kahteen sijoittamattomaan kontigiin, mikä toimii ohjenuorana näiden kontigien sijoittamisessa Ash1-kokoonpanon tulevissa versioissa.

Kartoitettuamme geenit Ash1:een poimimme koodaavat sekvenssit transkripteistä, jotka karttuivat täydellisesti (kattavuus 100 %), kohdistimme ne GRCh38:n koodaaviin sekvensseihin ja kutsuimme variantteja suhteessa GRCh38:aan (ks. jakso ”Menetelmät”). Näiden proteiineja koodaavien transkriptien 35 513 365 bp:n sisältä löytyi 20 864 yhden nukleotidin varianttia ja indeliä. Neljätoista tuhatta neljäsataa yhdeksänkymmentäyhdeksän näistä varianteista sijoittui GIAB:n ”kutsuttaville” alueille korkean luotettavuuden varianteille, joskin 3963 näistä sijaitsi GIAB:n ”vaikeilla” toistuvilla alueilla, joiden kohdistus on usein epäselvä. Niistä 10 536 variantista, jotka eivät sijainneet näillä vaikeilla alueilla, 10 456 (99,2 %) oli yhtäpitäviä GIAB:n korkean varmuuden varianttijoukon kanssa. Vaikeilla alueilla 3804/3963 (96,0 %) oli yhtäpitäviä GIAB-joukon kanssa.

Tämän jälkeen annotoimme varianttien ja epätäydellisen kartoituksen aiheuttamat muutokset aminohapoissa kaikkien proteiineja koodaavien sekvenssien osalta. 20 197 geenin 124 238:sta proteiinia koodaavasta transkriptiosta 92 600:lla (74,5 %) oli 100-prosenttisesti identtiset proteiinisekvenssit. Toisissa 26 566:ssa (21,4 %) on vähintään yksi aminohappomuutos, mutta ne tuottavat pituudeltaan identtisiä proteiineja, ja 1561:ssä (1,3 %) on runkoa säilyttäviä mutaatioita, jotka lisäävät tai poistavat yhden tai useampia aminohappoja jättäen loput proteiinista ennalleen. Taulukossa 6 esitetään tilastoja kaikista proteiinisekvenssien muutoksista. Jos proteiinissa oli useampi kuin yksi muunnos, laskimme sen seurauksiltaan merkittävimmän muunnoksen alle, eli jos proteiinissa oli missense-muunnos ja ennenaikainen stop-kodoni, laskimme sen ”stop gained”-ryhmään.

Erityisen kiinnostuksen kohteena ovat transkriptit, joiden variantit häiritsevät merkittävästi proteiinisekvenssiä ja saattavat johtaa toiminnan menetykseen. Näihin kuuluvat transkriptit, joihin vaikuttaa frameshift (2158), stop loss (58), stop gain (416), start loss (58) tai epätäydellisestä kartoituksesta johtuva typistyminen (564). Nämä häiriintyneet isomuodot edustavat 885 geenilokusta; 505:llä näistä geeneistä on kuitenkin vähintään yksi toinen isomuoto, johon häiriintynyt variantti ei vaikuta. Jäljelle jää 380 geeniä, joissa kaikissa isoformeissa on vähintään yksi häiriö; täydellinen luettelo on lisätiedostossa 1: taulukko S1.