Assemblaggio e annotazione di un genoma umano di riferimento Ashkenazi

Per la creazione del primo genoma umano di riferimento assemblato da un singolo individuo, abbiamo scelto HG002, un individuo Ashkenazi che fa parte del Personal Genome Project (PGP). Il PGP utilizza il modello di consenso aperto, la prima vera piattaforma ad accesso aperto per la condivisione di genoma umano individuale, fenotipo e dati medici. Il processo di consenso educa i potenziali partecipanti sulle implicazioni e i rischi della condivisione dei dati genomici e su ciò che possono aspettarsi dalla loro partecipazione. Il consenso aperto ha permesso la creazione dei primi materiali di riferimento del genoma umano al mondo (HG002 è il materiale di riferimento NIST 8391) da Genome In A Bottle (GIAB), che viene utilizzato per la calibrazione, lo sviluppo di metodi di assemblaggio del genoma e le misurazioni delle prestazioni del laboratorio. Tutti i dati di sequenza grezzi per questo progetto sono stati ottenuti da GIAB, dove è liberamente disponibile al pubblico.

Abbiamo assemblato il genoma HG002 da una combinazione di tre set di dati a copertura profonda: 249 bp Illumina legge, Oxford Nanopore (ONT) legge in media oltre 33 kbp di lunghezza, e PacBio di alta qualità “HiFi” legge in media 9567 bp (Tabella 1).

All’inizio abbiamo creato due assemblaggi: uno utilizzando Illumina e ONT letture, e un secondo utilizzando tutti e tre i set di dati, comprese le letture PacBio HiFi. L’aggiunta dei dati PacBio HiFi ha portato a una sequenza totale leggermente superiore nell’assemblaggio (2,99 Gb contro 2,88 Gb) e una dimensione N50 dei contig sostanzialmente più grande (18,2 Mb contro 4,9 Mb). Questo assemblaggio, designato Ash1 v0.5, è stato la base per tutti i successivi perfezionamenti.

Mappatura dell’assemblaggio sui cromosomi

Per creare assegnazioni cromosomiche per l’assemblaggio Ash1 v0.5, abbiamo usato allineamenti a GRCh38 per mappare la maggior parte degli scaffold sui cromosomi. I passaggi descritti nella sezione “Metodi” hanno generato una serie di assemblaggi su scala cromosomica gradualmente migliorati, con il risultato di Ash1 v1.7. Ash1 v1.7 ha una maggiore contiguità e minori lacune rispetto a GRCh38, come mostrato nella tabella 2. Si noti che nel processo di costruzione di questi cromosomi, una piccola quantità di sequenza di GRCh38 (58,3 Mb, 2% del genoma) è stata utilizzata per riempire le lacune in Ash1. Queste regioni includono alcune regioni difficili da assemblare che sono state curate manualmente per GRCh38. In totale, la dimensione stimata di tutte le lacune in Ash1 è 82,9 Mbp, contro 84,7 Mbp in GRCh38.p13.

Come parte del processo di miglioramento dell’assemblaggio, abbiamo cercato in uno degli assemblaggi preliminari di Ash1 (v1.1) le 12.745 varianti strutturali isolate di alta qualità (inserzioni e delezioni ≥ 50 bp) che Zook et al. hanno identificato confrontando i dati del trio Ashkenazi con GRCh37 . Questo studio ha utilizzato quattro diverse tecnologie di sequenziamento e molteplici chiamanti di varianti per identificare le varianti e filtrare i falsi positivi. Di queste 12.745 SV, 5807 sono omozigoti e 6938 sono eterozigoti. Ci aspettavamo che l’assemblaggio Ash1 fosse d’accordo con quasi tutte le varianti omozigoti. Poiché Ash1 cattura solo un aplotipo, ci aspettavamo che sarebbe stato d’accordo con circa la metà delle SV eterozigote, assumendo che l’algoritmo di assemblaggio scegliesse casualmente tra gli aplotipi quando decideva quale variante includere nel consenso finale. Delle 5807 varianti omozigote, 5284 (91%) erano presenti utilizzando i nostri criteri di corrispondenza (vedi la sezione “Metodi”), e 3922 (56,5%) delle 6938 varianti eterozigote erano presenti. Tutte le varianti sono state trovate nella posizione corretta.

Abbiamo anche fatto piccole (≤ 5 bp) chiamate di varianti su Ash1 v1.1 e le abbiamo confrontate con le varianti di riferimento HG002 v4.0 di GIAB, che abbiamo usato per correggere numerosi errori di sostituzione e indel (vedi la sezione “Metodi”), ottenendo Ash 1 v1.2. Abbiamo poi ri-allineato l’assemblaggio Ash1 a GRCh38, abbiamo ri-chiamato le varianti e abbiamo confrontato queste varianti con il nuovo set di riferimento GIAB v4.1. Delle varianti all’interno delle regioni di riferimento v4.1, le varianti Ash1 hanno abbinato 1.256.458 omozigoti e 1.041.476 eterozigoti SNPs, e 187.227 omozigoti e 193.524 eterozigoti indel. Dopo aver escluso le chiamate di varianti entro 30 bp di una vera variante, sono rimasti 79.269 SNPs e 17.439 indel, che (assumendo che questi siano tutti errori in Ash1) corrisponde a un valore di qualità (QV) di circa Q45 per gli errori di sostituzione. La maggior parte di queste varianti (52.191 SNPs e 4629 indel) cadono in duplicazioni segmentali, forse rappresentando duplicazioni mancanti in Ash1 o lucidatura imperfetta da letture brevi. In sintesi, la qualità dell’assemblaggio Ash1 è molto alta, con un valore di qualità di sostituzione stimato di 62 e un tasso di errore indel di 2 per milione di bp dopo aver escluso duplicazioni segmentali note, ripetizioni tandem e omopolimeri.

Confronto della chiamata delle varianti usando Ash1 contro GRCh38

Una delle motivazioni per la creazione di nuovi genomi di riferimento è che essi forniscono un quadro migliore per analizzare i dati di sequenza umana nella ricerca di varianti genetiche legate alla malattia. Per esempio, uno studio sugli ebrei ashkenaziti che ha raccolto dati WGS (whole-genome shotgun) dovrebbe usare un genoma di riferimento ashkenazita piuttosto che GRCh38. Poiché il background genetico è simile, meno varianti dovrebbero essere trovate quando si cerca contro Ash1.

Per testare questa aspettativa, abbiamo raccolto dati WGS da un partecipante maschio al Personal Genome Project, PGP17 (hu34D5B9). Questo individuo è stimato essere il 66% Ashkenazi secondo il database PGP, che era la più alta frazione stimata disponibile da individui PGP già sequenziati. Abbiamo scaricato 983.220.918.100 bp di letture (circa 33x di copertura) e le abbiamo allineate sia ad Ash1 che a GRCh38 usando Bowtie2 . Una frazione leggermente più alta di letture (3,901,270, 0.5%) si è allineata ad Ash1 che a GRCh38.

Abbiamo poi esaminato tutte le varianti a singolo nucleotide (SNVs, vedi “Metodi”) tra PGP17 e ciascuno dei due genomi di riferimento. Per semplificare l’analisi, abbiamo considerato solo le posizioni in cui PGP17 era omozigote. Nei nostri confronti con Ash1, abbiamo prima identificato tutti gli SNV e poi esaminato i dati originali di Ash1 per determinare se, per ciascuno di questi SNV, il genoma Ash1 conteneva un allele diverso che corrispondeva a PGP17.

In totale, il numero di siti omozigoti in PGP17 che non erano in accordo con Ash1 era 1.333.345, contro 1.700.364 quando abbiamo confrontato i siti omozigoti in PGP17 con GRCh38 (file aggiuntivo 1: Tabella S1). Abbiamo poi guardato i dati sottostanti Ash1 leggere per i 1,33 milioni di siti SNV che inizialmente non corrispondevano, e ha trovato che per circa la metà di loro, il genoma Ash1 era eterozigote, con un allele corrispondente PGP17. Se abbiamo ristretto gli SNV ai siti in cui PGP17 e Ash1 sono entrambi omozigoti (più un numero molto piccolo di posizioni in cui Ash1 contiene due varianti che differiscono entrambe da PGP17), questo ha ridotto il numero totale di SNV ancora di più, a 707.756. Così, abbiamo trovato poco meno di 1 milione di SNV omozigoti in meno quando si usa Ash1 come riferimento per PGP17. Si noti che piuttosto che allineare ad Ash1, si potrebbe invece allineare le letture a GRCh38 e poi rimuovere le varianti note specifiche della popolazione. Questo processo in due fasi, anche se più complesso, potrebbe dare risultati simili, tranne che per le regioni di Ash1 che sono semplicemente mancanti da GRCh38.

Confronto con le varianti comuni Ashkenazi

Per esaminare la misura in cui Ash1 contiene noti, varianti comuni Ashkenazi (rispetto a GRCh38), abbiamo esaminato SNVs riportati ad alta frequenza in una popolazione Ashkenazi dal Genome Aggregation Database (gnomAD) . GnomAD v3.0 contiene le chiamate SNV dai dati del genoma intero a lettura breve di 1662 individui ashkenaziti. Poiché alcune varianti sono state chiamate solo in un sottoinsieme di questi individui, abbiamo considerato solo i siti di varianti che sono stati riportati in un minimo di 200 persone. Abbiamo poi raccolto gli SNV degli alleli maggiori, richiedendo che la frequenza allelica fosse superiore a 0,5 nella popolazione campionata. Abbiamo ulteriormente limitato la nostra analisi alle sostituzioni di una singola base. Questo ci ha dato 2.008.397 siti gnomAD SNV in cui l’allele maggiore ashkenazita differiva da GRCh38.

Siamo stati in grado di mappare con precisione 1.790.688 dei 2.008.397 siti gnomAD di GRCh38 su Ash1 (vedi la sezione “Metodi”). Abbiamo poi confrontato la base di GRCh38 con la base dell’allele maggiore Ashkenazi in ogni posizione, e abbiamo anche esaminato gli alleli alternativi in Ash1 nei siti in cui è eterozigote. Per i siti in cui i dati di lettura hanno mostrato che HG002 era eterozigote e aveva sia l’allele maggiore Ashkenazi che l’allele GRCh38, abbiamo sostituito la base Ash1, se necessario, per garantire che corrispondesse all’allele maggiore. Dopo queste sostituzioni, Ash1 conteneva l’allele maggiore ashkenazita nell’88% (1.580.866) dei 1,79 milioni di siti. Nei siti rimanenti, Ash1 o corrispondeva all’allele GRCh38 perché HG002 è omozigote per l’allele di riferimento (204.729 siti), o conteneva un terzo allele che non corrispondeva né a GRCh38 né all’allele maggiore gnomAD (5093 siti). Queste modifiche dovrebbero ridurre ulteriormente il numero di SNV riportati quando si allinea un individuo ashkenazita ad Ash1.

Degno di nota è che, con l’aumentare della frequenza dell’allele maggiore nella popolazione ashkenazita gnomAD, aumenta anche la proporzione di siti in cui Ash1 corrisponde all’allele maggiore. Per esempio, di SNVs che hanno una frequenza allelica > 0,9 nella popolazione gnomAD Ashkenazi, oltre il 98% sono rappresentati in Ash1 (Tabella 3).

Annotazione

Una parte vitale di ogni genoma di riferimento è l’annotazione: la raccolta di tutti i geni e altre caratteristiche trovate sul genoma. Per permettere ad Ash1 di funzionare come un vero genoma di riferimento, abbiamo mappato tutti i geni conosciuti usati dalla comunità scientifica sul suo sistema di coordinate, usando gli stessi nomi di geni e identificatori. Diversi database di annotazioni sono stati creati per GRCh38, e per Ash1, abbiamo scelto di utilizzare il database dei geni umani CHESS perché è completo, includendo tutti i geni codificanti le proteine sia in GENCODE che in RefSeq, gli altri due database di geni ampiamente utilizzati, e perché mantiene gli identificatori utilizzati in questi cataloghi. I geni non codificanti differiscono tra i tre database, ma CHESS ha il maggior numero di loci genici e isoforme. Abbiamo usato CHESS versione 2.2, che ha 42.167 geni sui cromosomi primari (esclusi gli scaffold alternativi GRCh38), di cui 20.197 sono codificanti per le proteine.

Mappare i geni da un assemblaggio all’altro è un compito complesso, in particolare per i geni che si presentano in famiglie di geni altamente simili e con più copie. Il compito è ancora più complesso quando i due assemblaggi rappresentano individui diversi (piuttosto che semplicemente diversi assemblaggi dello stesso individuo), a causa della presenza di differenze di singolo nucleotide, inserzioni, delezioni, riarrangiamenti, e genuine variazioni nel numero di copie tra gli individui. Avevamo bisogno di un metodo che fosse robusto di fronte a tutte queste potenziali differenze.

Per affrontare questo problema, abbiamo usato lo strumento di mappatura Liftoff recentemente sviluppato, che nei nostri esperimenti è stato l’unico strumento che poteva mappare quasi tutti i geni umani da un individuo all’altro. Liftoff prende tutti i geni e le trascrizioni da un genoma e li mappa, cromosoma per cromosoma, su un altro genoma. Per tutti i geni che non riescono a mappare sullo stesso cromosoma, Liftoff cerca di mapparli attraverso i cromosomi. A differenza di altri strumenti, non si basa su una mappa dettagliata costruita da un allineamento dell’intero genoma, ma invece mappa ogni gene individualmente. I geni sono allineati a livello di trascrizione, compresi gli introni, in modo che gli pseudogeni elaborati non vengano erroneamente identificati come geni.

Abbiamo cercato di mappare tutti i 310.901 trascritti da 42.167 loci genici sui cromosomi primari in GRCh38 ad Ash1. In totale, abbiamo mappato con successo 309.900 (99,7%) trascrizioni da 42.070 loci genici sui cromosomi principali (file aggiuntivo 1: tabella S2). Abbiamo considerato una trascrizione mappata con successo se la sequenza di mRNA in Ash1 è almeno il 50% più lunga della sequenza di mRNA su GRCh38. Tuttavia, la stragrande maggioranza delle trascrizioni supera ampiamente questa soglia, con il 99% delle trascrizioni mappate con una copertura maggiore o uguale al 95% (Additional file 2: Figura S2). L’identità di sequenza delle trascrizioni mappate è ugualmente alta, con il 99% delle trascrizioni mappate con un’identità di sequenza maggiore o uguale al 94% (Additional file 2: Figure S3).

Geni traslocati

Di quei geni con almeno una isoforma mappata con successo, 42.059 (99,7%) hanno mappato nelle posizioni corrispondenti sullo stesso cromosoma in Ash1. Dei 108 geni che inizialmente non sono riusciti a mappare, 11 geni hanno mappato su un cromosoma diverso in 7 blocchi distinti (indicati nella tabella 4), suggerendo una traslocazione tra i due genomi. È interessante notare che 16 delle 22 posizioni coinvolte nelle traslocazioni erano in regioni subtelomeriche, che si sono verificate in 8 coppie in cui entrambe le posizioni erano vicine ai telomeri. Questo è coerente con studi precedenti che riportano che i riarrangiamenti che coinvolgono telomeri e subtelomeri possono essere una forma comune di traslocazione negli esseri umani.

Abbiamo esaminato la traslocazione tra i cromosomi 15 e 20, che contiene tre dei geni nella tabella 4, guardando più da vicino l’allineamento tra GRCh38 e Ash1. La traslocazione è al telomero di entrambi i cromosomi, dalla posizione 65.079.275 a 65.109.824 (30.549 bp) di Ash1 chr20 e da 101.950.338 a 101.980.928 (30.590 bp) di GRCh39 chr15. Per confermare la traslocazione, abbiamo allineato un set indipendente di letture PacBio molto lunghe, tutte da HG002, all’assembly Ash1 v1.7 (vedi la sezione “Metodi”) e abbiamo valutato la regione intorno al breakpoint su chr20. Questi allineamenti mostrano una copertura profonda e coerente che si estende per molti kilobasi su entrambi i lati del breakpoint, sostenendo la correttezza dell’assemblaggio Ash1 (Fig. 1).

Snapshot che mostra gli allineamenti di lunghe letture PacBio al genoma Ash1, centrato sull’estremità sinistra della posizione nel cromosoma 20 (posizione 65.079.275) dove si è verificata una traslocazione tra il cromosoma 15 (GRCh38) e 20 (Ash1). La parte superiore della figura mostra le coordinate sul chr20. Sotto c’è un istogramma della copertura delle letture, e le singole letture riempiono la parte inferiore della figura. Gli indel nelle letture, mostrati come barre colorate su ogni lettura, sono dovuti al tasso di errore relativamente alto delle letture lunghe

Geni mancanti

Sessantadue geni non sono riusciti interamente a mappare da GRCh38 su Ash1, e altri 32 geni hanno mappato solo parzialmente (sotto la soglia di copertura del 50%), come mostrato nella tabella 5. Tutti i geni che non sono riusciti a mappare o che hanno mappato parzialmente erano membri di famiglie multi-gene, e in ogni caso, c’era almeno un’altra copia del gene mancante presente in Ash1, con un’identità media del 98,5%. Quindi, non ci sono affatto casi di un gene che è presente in GRCh38 e che è interamente assente in Ash1; i geni mostrati nella tabella 5 rappresentano casi in cui Ash1 ha meno membri di una famiglia multigenica. Tre geni aggiuntivi (2 codificanti proteine, 1 lncRNA) mappati a due contig non posizionati, che forniranno una guida per posizionare quei contig nelle future versioni dell’assembly Ash1.

Dopo aver mappato i geni su Ash1, abbiamo estratto le sequenze codificanti dai trascritti che hanno mappato completamente (copertura pari al 100%), li abbiamo allineati alle sequenze codificanti di GRCh38 e abbiamo chiamato le varianti relative a GRCh38 (vedi la sezione “Metodi”). All’interno dei 35.513.365 bp di queste trascrizioni codificanti le proteine, abbiamo trovato 20.864 varianti a singolo nucleotide e indel. Quattordicimila quattrocento novantanove di queste varianti sono caduti all’interno del GIAB “callable” regioni per le varianti ad alta fiducia, anche se 3963 di questi erano in GIAB “difficile” regioni ripetitive, per i quali gli allineamenti sono spesso ambigui. Delle 10.536 varianti non presenti in queste regioni difficili, 10.456 (99,2%) hanno concordato con il set di varianti ad alta confidenza GIAB. Nelle regioni difficili, 3804/3963 (96,0%) hanno concordato con il set GIAB.

Abbiamo poi annotato i cambiamenti negli aminoacidi causati dalle varianti e dalla mappatura incompleta per tutte le sequenze codificanti le proteine. Su 124.238 trascrizioni codificanti proteine da 20.197 geni, 92.600 (74,5%) hanno sequenze proteiche identiche al 100%. Altri 26.566 (21,4%) hanno almeno un cambiamento di aminoacido ma producono proteine con la stessa lunghezza, e 1561 (1,3%) hanno mutazioni che conservano il frame che inseriscono o eliminano uno o più aminoacidi, lasciando il resto della proteina invariato. La tabella 6 mostra le statistiche su tutti i cambiamenti nelle sequenze delle proteine. Se una proteina aveva più di 1 variante, l’abbiamo contata sotto la variante più conseguente, cioè, se una proteina aveva una variante missenso e un codone di stop prematuro, l’abbiamo inclusa nel gruppo “stop guadagnato”.

Di particolare interesse sono quelle trascrizioni con varianti che interrompono significativamente la sequenza proteica e possono risultare in perdita di funzione. Questi includono trascrizioni affette da un frameshift (2158), perdita di stop (58), guadagno di stop (416), perdita di inizio (58), o troncamento dovuto alla mappatura incompleta (564). Queste isoforme interrotte rappresentano 885 loci genici; tuttavia, 505 di questi geni hanno almeno 1 altra isoforma che non è influenzata da una variante di interruzione. Questo lascia 380 geni in cui tutte le isoforme hanno almeno una interruzione; l’elenco completo è fornito nel file aggiuntivo 1: Tabella S1.