Asamblarea și adnotarea unui genom uman de referință Ashkenazi

Pentru crearea primului genom uman de referință care să fie asamblat de la un singur individ, am ales HG002, un individ Ashkenazi care face parte din Proiectul Genomului Personal (PGP). PGP utilizează modelul Open Consent, prima platformă cu adevărat cu acces deschis pentru partajarea genomului uman individual, a fenotipului și a datelor medicale . Procesul de consimțământ îi educă pe potențialii participanți cu privire la implicațiile și riscurile partajării datelor genomice, precum și cu privire la ceea ce se pot aștepta de la participarea lor. Consimțământul deschis a permis crearea primelor materiale de referință pentru genomul uman din lume (HG002 este materialul de referință NIST 8391) de la Genome In A Bottle (GIAB), care este utilizat pentru calibrare, dezvoltarea metodelor de asamblare a genomului și măsurarea performanțelor laboratorului . Toate datele de secvență brută pentru acest proiect au fost obținute de la GIAB, unde sunt disponibile în mod gratuit pentru public .

Am asamblat genomul HG002 dintr-o combinație de trei seturi de date cu acoperire profundă: lecturi Illumina de 249 bp, lecturi Oxford Nanopore (ONT) cu o lungime medie de peste 33 kbp și lecturi PacBio „HiFi” de înaltă calitate cu o lungime medie de 9567 bp (tabelul 1).

Am creat inițial două asamblări: una folosind lecturi Illumina și ONT, iar a doua folosind toate cele trei seturi de date, inclusiv lecturile PacBio HiFi. Adăugarea datelor PacBio HiFi a avut ca rezultat o secvență totală ușor mai mare în ansamblu (2,99 Gb față de 2,88 Gb) și o dimensiune N50 a contigurilor substanțial mai mare (18,2 Mb față de 4,9 Mb). Acest ansamblu, denumit Ash1 v0.5, a stat la baza tuturor perfecționărilor ulterioare.

Corectarea ansamblului pe cromozomi

Pentru a crea atribuiri de cromozomi pentru ansamblul Ash1 v0.5, am folosit alinieri la GRCh38 pentru a cartografia majoritatea scheletelor pe cromozomi. Etapele descrise în secțiunea „Metode” au generat o serie de ansambluri la scară cromozomială îmbunătățite treptat, rezultând Ash1 v1.7. Ash1 v1.7 are o contiguitate mai mare și lacune mai mici decât GRCh38, după cum se arată în tabelul 2. Rețineți că, în procesul de construire a acestor cromozomi, o cantitate mică de secvență din GRCh38 (58,3 Mb, 2% din genom) a fost utilizată pentru a umple lacunele din Ash1. Aceste regiuni includ unele regiuni greu de asamblat care au fost curatoriate manual pentru GRCh38. În total, dimensiunea estimată a tuturor lacunelor din Ash1 este de 82,9 Mbp, față de 84,7 Mbp în GRCh38.p13.

Ca parte a procesului de îmbunătățire a ansamblului, am căutat într-unul dintre ansamblurile preliminare Ash1 (v1.1) pentru cele 12 745 de variante structurale izolate de înaltă calitate (inserții și deleții ≥ 50 pb) pe care Zook et al. le-au identificat prin compararea datelor trio-ului Ashkenazi cu GRCh37 . Studiul respectiv a utilizat patru tehnologii de secvențiere diferite și mai multe apeluri de variante pentru a identifica variantele și a filtra falsurile pozitive. Dintre aceste 12.745 de VS, 5807 sunt homozigote și 6938 sunt heterozigote. Ne așteptam ca ansamblul Ash1 să fie de acord cu aproape toate variantele homozigote. Deoarece Ash1 captează doar un singur haplotip, ne-am așteptat ca acesta să fie de acord cu aproximativ jumătate dintre VS heterozigote, presupunând că algoritmul de asamblare a ales aleatoriu între haplotipuri atunci când a decis ce variantă să includă în consensul final. Din cele 5807 variante homozigote, 5284 (91 %) au fost prezente folosind criteriile noastre de potrivire (a se vedea secțiunea „Metode”), iar 3922 (56,5 %) din 6938 de variante heterozigote au fost prezente. Toate variantele au fost găsite la locația corectă.

Am efectuat, de asemenea, apeluri de variante mici (≤ 5 pb) pe Ash1 v1.1 și le-am comparat cu variantele de referință HG002 v4.0 de la GIAB, pe care le-am folosit pentru a corecta numeroase erori de substituție și indel (a se vedea secțiunea „Metode”), obținând Ash 1 v1.2. Apoi am realiniat ansamblul Ash1 la GRCh38, am redenumit variantele și am comparat aceste variante cu setul de referință GIAB v4.1 nou dezvoltat. Dintre variantele din interiorul regiunilor de referință v4.1, variantele Ash1 au corespuns cu 1 256 458 de SNP homozigoți și 1 041 476 de SNP heterozigoți, precum și cu 187 227 de indeli homozigoți și 193 524 de indeli heterozigoți. După excluderea apelurilor de variante aflate la mai puțin de 30 pb de o variantă adevărată, au rămas 79.269 de SNP-uri și 17.439 de indeli, ceea ce (presupunând că toate acestea sunt erori în Ash1) corespunde unei valori de calitate (QV) de aproximativ Q45 pentru erorile de substituție. Cele mai multe dintre aceste variante (52 191 SNP și 4629 indels) se încadrează în duplicații segmentare, reprezentând probabil duplicații lipsă în Ash1 sau o lustruire imperfectă prin lecturi scurte. În concluzie, calitatea ansamblului Ash1 este foarte ridicată, cu o valoare estimată a calității de substituție de 62 și o rată a erorilor indel de 2 la un milion de bp după excluderea duplicațiilor segmentare cunoscute, a repetițiilor în tandem și a homopolimerilor.

Compararea apelării variantelor folosind Ash1 față de GRCh38

Una dintre motivațiile pentru crearea de noi genomuri de referință este că acestea oferă un cadru mai bun pentru analiza datelor de secvență umană atunci când se caută variante genetice legate de boli. De exemplu, un studiu al evreilor Ashkenazi care a colectat date de tip „whole-genome shotgun” (WGS) ar trebui să utilizeze un genom de referință Ashkenazi mai degrabă decât GRCh38. Deoarece fondul genetic este similar, ar trebui să se găsească mai puține variante atunci când se caută față de Ash1.

Pentru a testa această așteptare, am colectat date WGS de la un bărbat participant la Personal Genome Project, PGP17 (hu34D5B9). Acest individ este estimat a fi 66% Ashkenazi conform bazei de date PGP, care a fost cea mai mare fracție estimată disponibilă de la indivizii PGP deja secvențializați. Am descărcat 983.220.918.918.100 de citiri de 100 pb (acoperire de aproximativ 33x) și le-am aliniat atât la Ash1, cât și la GRCh38 folosind Bowtie2 . O fracțiune ușor mai mare de citiri (3 901 270, 0,5 %) s-a aliniat la Ash1 decât la GRCh38.

Apoi am examinat toate variantele cu un singur nucleotid (SNV, a se vedea „Metode”) între PGP17 și fiecare dintre cele două genomuri de referință. Pentru a simplifica analiza, am luat în considerare doar locațiile în care PGP17 era homozigot. În comparațiile noastre cu Ash1, am identificat mai întâi toate SNV-urile și apoi am examinat datele de citire originale Ash1 pentru a determina dacă, pentru fiecare dintre aceste SNV-uri, genomul Ash1 conținea o alelă diferită care se potrivea cu PGP17.

În total, numărul de situri homozigote din PGP17 care au fost în dezacord cu Ash1 a fost de 1.333.345, față de 1.700.364 atunci când am comparat siturile homozigote din PGP17 cu GRCh38 (Fișierul suplimentar 1: Tabelul S1). Am analizat apoi datele de citire Ash1 subiacente pentru cele 1,33 milioane de situri SNV care inițial nu corespundeau și am constatat că, pentru aproximativ jumătate dintre acestea, genomul Ash1 era heterozigot, cu o alelă care corespundea cu PGP17. Dacă am restrâns SNV-urile la site-urile în care PGP17 și Ash1 sunt ambele homozigote (plus un număr foarte mic de locații în care Ash1 conține două variante care diferă ambele de PGP17), acest lucru a redus și mai mult numărul total de SNV-uri, la 707 756. Astfel, am găsit cu puțin sub 1 milion mai puține SNV homozigote atunci când am folosit Ash1 ca referință pentru PGP17. Rețineți că, în loc să se alinieze la Ash1, s-ar putea alinia citirile la GRCh38 și apoi să se elimine variantele cunoscute specifice populației. Acest proces în doi pași, deși mai complex, ar putea da rezultate similare, cu excepția regiunilor din Ash1 care pur și simplu lipsesc din GRCh38.

Comparare cu variantele comune Ashkenazi

Pentru a examina măsura în care Ash1 conține variante Ashkenazi cunoscute și comune (în raport cu GRCh38), am examinat SNV-urile raportate cu frecvență ridicată într-o populație Ashkenazi din Genome Aggregation Database (gnomAD) . GnomAD v3.0 conține apeluri SNV din datele de genom întreg cu citire scurtă de la 1662 de indivizi Ashkenazi. Deoarece unele variante au fost apelate doar la un subset al acestor indivizi, am luat în considerare doar site-urile de variante care au fost raportate la un minim de 200 de persoane. Am colectat apoi SNV-urile alelelor majore, necesitând ca frecvența alelei să fie mai mare de 0,5 în populația eșantionată. Am restrâns în continuare analiza noastră la substituțiile cu o singură bază. Acest lucru ne-a oferit 2.008.397 de situri SNV gnomAD în care alela majoră Ashkenazi diferă de GRCh38.

Am reușit să cartografiem cu precizie 1.790.688 din cele 2.008.397 de situri gnomAD din GRCh38 pe Ash1 (a se vedea secțiunea „Metode”). Am comparat apoi baza GRCh38 cu baza alelei majore Ashkenazi la fiecare poziție și am examinat, de asemenea, alelele alternative din Ash1 în siturile în care este heterozigotă. Pentru locurile în care datele de citire au arătat că HG002 era heterozigotă și avea atât alela majoră Ashkenazi, cât și alela GRCh38, am înlocuit baza Ash1, dacă a fost necesar, pentru a ne asigura că aceasta corespundea cu alela majoră. După aceste înlocuiri, Ash1 conținea alela majoră Ashkenazi în 88% (1.580.866) din cele 1,79 milioane de situri. La siturile rămase, Ash1 fie corespundea alelei GRCh38, deoarece HG002 este homozigotă pentru alela de referință (204.729 situri), fie conținea o a treia alelă care nu corespundea nici cu GRCh38, nici cu alela majoră gnomAD (5093 situri). Aceste modificări ar trebui să reducă și mai mult numărul de SNV-uri raportate atunci când se aliniază un individ Ashkenazi la Ash1.

Merită remarcat faptul că, pe măsură ce crește frecvența alelei majore în populația Ashkenazi gnomAD, crește și proporția de situri în care Ash1 se potrivește cu alela majoră. De exemplu, dintre SNV-urile care au o frecvență a alelei > 0,9 în populația Ashkenazi gnomAD, peste 98% sunt reprezentate în Ash1 (Tabelul 3).

Anotație

O parte vitală a oricărui genom de referință este adnotarea: colectarea tuturor genelor și a altor caracteristici găsite în genom. Pentru a permite ca Ash1 să funcționeze ca un adevărat genom de referință, am cartografiat toate genele cunoscute utilizate de comunitatea științifică pe sistemul său de coordonate, folosind aceleași nume și identificatori de gene. Au fost create mai multe baze de date de adnotare diferite pentru GRCh38, iar pentru Ash1, am ales să folosim baza de date de gene umane CHESS, deoarece este cuprinzătoare, incluzând toate genele codificatoare de proteine atât în GENCODE, cât și în RefSeq, celelalte două baze de date de gene utilizate pe scară largă, și deoarece păstrează identificatorii utilizați în aceste cataloage. Genele necodificatoare diferă între cele trei baze de date, dar CHESS are cel mai mare număr de loci și izoforme de gene. Am folosit versiunea 2.2 a CHESS, care are 42 167 de gene pe cromozomii primari (excluzând schelele alternative GRCh38), dintre care 20 197 sunt codificatoare de proteine.

Corectarea genelor de la un ansamblu la altul este o sarcină complexă, în special pentru genele care apar în familii de gene foarte asemănătoare, cu mai multe copii. Sarcina este și mai complexă atunci când cele două ansambluri reprezintă indivizi diferiți (mai degrabă decât pur și simplu ansambluri diferite ale aceluiași individ), din cauza prezenței diferențelor de un singur nucleotid, a inserțiilor, a delețiilor, a rearanjamentelor și a variațiilor autentice ale numărului de copii între indivizi. Aveam nevoie de o metodă care să fie robustă în fața tuturor acestor potențiale diferențe.

Pentru a aborda această problemă, am folosit instrumentul de cartografiere Liftoff, dezvoltat recent, care, în experimentele noastre, a fost singurul instrument care a putut cartografia aproape toate genele umane de la un individ la altul. Liftoff ia toate genele și transcriptele dintr-un genom și le mapează, cromozom cu cromozom, la un alt genom. Pentru toate genele care nu reușesc să fie cartografiate pe același cromozom, Liftoff încearcă să le cartografieze între cromozomi. Spre deosebire de alte instrumente, acesta nu se bazează pe o hartă detaliată construită dintr-o aliniere a întregului genom, ci, în schimb, cartografiază fiecare genă în parte. Genele sunt aliniate la nivel de transcripție, inclusiv intronii, astfel încât pseudogenele procesate să nu fie identificate în mod eronat ca gene.

Am încercat să cartografiem toate cele 310.901 transcripții de la 42.167 de loci genetici de pe cromozomii primari din GRCh38 la Ash1. În total, am reușit să cartografiem cu succes 309.900 (99,7%) transcripte din 42.070 loci genetici pe cromozomii principali (Fișier suplimentar 1: Tabelul S2). Am considerat că un transcript a fost cartografiat cu succes dacă secvența de ARNm din Ash1 este cel puțin 50% mai lungă decât secvența de ARNm din GRCh38. Cu toate acestea, marea majoritate a transcriptelor depășesc cu mult acest prag, 99% dintre transcripte fiind cartografiate la o acoperire mai mare sau egală cu 95% (Fișier suplimentar 2: Figura S2). Identitatea de secvență a transcriptelor cartografiate este la fel de ridicată, cu 99% dintre transcripte cartografiate cu o identitate de secvență mai mare sau egală cu 94% (Fișier suplimentar 2: Figura S3).

Gene translocate

Dintre genele cu cel puțin o izoformă cartografiată cu succes, 42 059 (99,7%) au fost cartografiate la locațiile corespunzătoare pe același cromozom în Ash1. Din cele 108 gene care inițial nu au reușit să se cartografieze, 11 gene s-au cartografiat pe un cromozom diferit în 7 blocuri distincte (prezentate în tabelul 4), sugerând o translocație între cele două genomuri. Este interesant faptul că 16 din cele 22 de locații implicate în translocații se aflau în regiuni subtelomerice, care au apărut în 8 perechi în care ambele locații se aflau în apropierea telomerilor. Acest lucru este în concordanță cu studiile anterioare care au raportat că rearanjamentele care implică telomeri și subtelomeri pot fi o formă comună de translocație la om .

Am examinat translocația dintre cromozomii 15 și 20, care conține trei dintre genele din tabelul 4, analizând mai atent alinierea dintre GRCh38 și Ash1. Translocația se află la nivelul telomerului ambilor cromozomi, de la poziția 65,079,275 la 65,109,824 (30,549 pb) din Ash1 chr20 și de la 101,950,338 la 101,980,928 (30,590 pb) din GRCh39 chr15. Pentru a confirma translocația, am aliniat un set independent de citiri PacBio foarte lungi, toate din HG002, la ansamblul Ash1 v1.7 (a se vedea secțiunea „Metode”) și am evaluat regiunea din jurul punctului de ruptură de pe chr20. Aceste alinieri arată o acoperire profundă și consistentă care se extinde pe mai multe kilobaze de ambele părți ale punctului de ruptură, susținând corectitudinea ansamblului Ash1 (Fig. 1).

Imagine care arată alinierile citirilor lungi PacBio la genomul Ash1, centrate pe capătul stâng al locației din cromozomul 20 (poziția 65,079,275) unde a avut loc o translocație între cromozomul 15 (GRCh38) și 20 (Ash1). Partea de sus a figurii arată coordonatele de pe chr20. Sub aceasta se află o histogramă a acoperirii citirilor, iar citirile individuale completează partea de jos a figurii. Indelurile din lecturi, prezentate sub formă de bare colorate pe fiecare citire, se datorează ratei relativ ridicate de eroare a citirilor lungi

Gene lipsă

Șaizeci și două de gene nu au reușit să se cartografieze în întregime de pe GRCh38 pe Ash1, iar alte 32 de gene s-au cartografiat doar parțial (sub pragul de acoperire de 50 %), după cum se arată în tabelul 5. Toate genele care nu au reușit să se cartografieze sau care s-au cartografiat parțial erau membre ale unor familii de mai multe gene și, în fiecare caz, a existat cel puțin o altă copie a genei lipsă prezentă în Ash1, cu o identitate medie de 98,5 %. Astfel, nu există niciun caz în care o genă prezentă în GRCh38 să fie complet absentă din Ash1; genele prezentate în tabelul 5 reprezintă cazuri în care Ash1 are mai puțini membri ai unei familii multigene. Trei gene suplimentare (2 codificatoare de proteine, 1 lncRNA) au fost cartografiate la două contig-uri neplasate, ceea ce va oferi un ghid pentru plasarea acestor contig-uri în versiunile viitoare ale ansamblului Ash1.

După cartografierea genelor pe Ash1, am extras secvențele codificatoare din transcriptele care au fost cartografiate complet (acoperire egală cu 100%), le-am aliniat la secvențele codificatoare din GRCh38 și am numit variantele în raport cu GRCh38 (a se vedea secțiunea „Metode”). În cadrul celor 35.513.365 pb din aceste transcripte codificatoare de proteine, am găsit 20.864 de variante și indeluri de un singur nucleotid. Paisprezece mii patru sute nouăzeci și nouă dintre aceste variante au fost încadrate în regiunile GIAB „apelabile” pentru variantele cu grad ridicat de încredere, deși 3963 dintre acestea se aflau în regiunile repetitive GIAB „dificile”, pentru care alinierile sunt adesea ambigue. Din cele 10 536 de variante care nu se aflau în aceste regiuni dificile, 10 456 (99,2 %) au fost în concordanță cu setul de variante GIAB cu grad ridicat de încredere. În regiunile dificile, 3804/3963 (96,0%) au fost de acord cu setul GIAB.

Apoi am notat modificările de aminoacizi cauzate de variante și de cartografierea incompletă pentru toate secvențele codificatoare de proteine. Din 124.238 de transcripții codificatoare de proteine din 20.197 de gene, 92.600 (74,5%) au secvențe proteice 100% identice. Alte 26.566 (21,4%) au cel puțin o modificare de aminoacizi, dar produc proteine cu lungime identică, iar 1561 (1,3%) au mutații de păstrare a cadrului care inserează sau elimină unul sau mai mulți aminoacizi, lăsând restul proteinei neschimbate. Tabelul 6 prezintă statistici cu privire la toate modificările secvențelor de proteine. Dacă o proteină a avut mai mult de 1 variantă, am numărat-o la varianta cea mai consecventă, de exemplu, dacă o proteină a avut o variantă missense și un codon de oprire prematură, o includem în grupul „stop câștigat”.

De un interes deosebit sunt acele transcripte cu variante care perturbă semnificativ secvența proteică și pot duce la pierderea funcției. Printre acestea se numără transcripții afectate de o deplasare de cadru (2158), pierdere de stop (58), câștig de stop (416), pierdere de start (58) sau trunchiere datorată unei cartografieri incomplete (564). Aceste izoforme perturbate reprezintă 885 de loci genetici; cu toate acestea, 505 dintre aceste gene au cel puțin o altă izoformă care nu este afectată de o variantă perturbatoare. Rămân 380 de gene în care toate izoformele au cel puțin o perturbare; lista completă este furnizată în Fișierul suplimentar 1: Tabelul S1.

.