Zusammenbau und Annotation eines Ashkenazi-Human-Referenzgenoms

Für die Erstellung des ersten menschlichen Referenzgenoms, das aus einem einzelnen Individuum zusammengestellt wurde, wählten wir HG002, ein Ashkenazi-Individuum, das Teil des Personal Genome Project (PGP) ist. Das PGP verwendet das Open Consent Model, die erste wirklich frei zugängliche Plattform für die gemeinsame Nutzung von individuellen menschlichen Genom-, Phänotyp- und medizinischen Daten. Der Einwilligungsprozess klärt potenzielle Teilnehmer über die Auswirkungen und Risiken der Weitergabe von Genomdaten auf und darüber, was sie von ihrer Teilnahme erwarten können. Die offene Zustimmung hat die Erstellung des weltweit ersten Referenzmaterials für das menschliche Genom (HG002 ist das NIST-Referenzmaterial 8391) von Genome In A Bottle (GIAB) ermöglicht, das für die Kalibrierung, die Entwicklung von Genomassemblierungsmethoden und die Messung der Laborleistung verwendet wird. Alle Rohsequenzdaten für dieses Projekt wurden von GIAB bezogen, wo sie für die Öffentlichkeit frei zugänglich sind.

Wir assemblierten das HG002-Genom aus einer Kombination von drei Datensätzen mit hoher Abdeckung: 249-bp Illumina-Reads, Oxford Nanopore (ONT) Reads mit einer durchschnittlichen Länge von über 33 kbp und hochwertige PacBio „HiFi“-Reads mit einer durchschnittlichen Länge von 9567 bp (Tabelle 1).

Wir haben zunächst zwei Assemblierungen erstellt: eine unter Verwendung von Illumina- und ONT-Reads und eine zweite unter Verwendung aller drei Datensätze, einschließlich der PacBio HiFi-Reads. Die Hinzufügung der PacBio HiFi-Daten führte zu einer geringfügig größeren Gesamtsequenz in der Assemblierung (2,99 Gb vs. 2,88 Gb) und einer wesentlich größeren N50-Kontig-Größe (18,2 Mb vs. 4,9 Mb). Diese Assemblierung, die als Ash1 v0.5 bezeichnet wurde, war die Grundlage für alle nachfolgenden Verfeinerungen.

Mapping der Assemblierung auf Chromosomen

Um Chromosomenzuordnungen für die Ash1 v0.5-Assemblierung zu erstellen, verwendeten wir Alignments zu GRCh38, um die meisten Gerüste auf Chromosomen zuzuordnen. Die im Abschnitt „Methoden“ beschriebenen Schritte führten zu einer Reihe von schrittweise verbesserten Chromosomenzuordnungen, die in Ash1 v1.7 mündeten. Ash1 v1.7 weist eine größere Kontiguität und kleinere Lücken als GRCh38 auf, wie in Tabelle 2 dargestellt. Es ist zu beachten, dass bei der Erstellung dieser Chromosomen ein kleiner Teil der GRCh38-Sequenz (58,3 Mb, 2 % des Genoms) verwendet wurde, um Lücken in Ash1 zu füllen. Diese Regionen umfassen einige schwer zu assemblierende Regionen, die manuell für GRCh38 kuratiert wurden. Insgesamt beträgt die geschätzte Größe aller Lücken in Ash1 82,9 Mbp, gegenüber 84,7 Mbp in GRCh38.p13.

Als Teil des Prozesses zur Verbesserung der Assemblierung durchsuchten wir eine der vorläufigen Ash1-Assemblierungen (v1.1) nach den 12.745 qualitativ hochwertigen, isolierten strukturellen Varianten (Insertionen und Deletionen ≥ 50 bp), die Zook et al. durch den Vergleich der Ashkenazi-Trio-Daten mit GRCh37 identifiziert hatten. In dieser Studie wurden vier verschiedene Sequenziertechnologien und mehrere Varianten-Caller verwendet, um Varianten zu identifizieren und falsch-positive Ergebnisse herauszufiltern. Von diesen 12.745 SVs sind 5807 homozygot und 6938 heterozygot. Wir erwarteten, dass die Ash1-Zusammenstellung mit fast allen homozygoten Varianten übereinstimmen würde. Da Ash1 nur einen Haplotyp erfasst, erwarteten wir, dass es mit etwa der Hälfte der heterozygoten SVs übereinstimmen würde, wenn man davon ausgeht, dass der Assembler-Algorithmus bei der Entscheidung, welche Variante in den endgültigen Konsens aufgenommen wird, zufällig zwischen den Haplotypen wählt. Von den 5807 homozygoten Varianten waren 5284 (91 %) nach unseren Übereinstimmungskriterien (siehe Abschnitt „Methoden“) vorhanden, und 3922 (56,5 %) von 6938 heterozygoten Varianten waren vorhanden. Alle Varianten wurden an der richtigen Stelle gefunden.

Wir haben auch kleine (≤ 5 bp) Variantenaufrufe auf Ash1 v1.1 durchgeführt und diese mit den HG002 v4.0 Benchmark-Varianten von GIAB verglichen, mit denen wir zahlreiche Substitutions- und Indel-Fehler korrigiert haben (siehe Abschnitt „Methoden“), was zu Ash 1 v1.2 führte. Anschließend haben wir die Ash1-Assembly an GRCh38 neu ausgerichtet, Varianten neu aufgerufen und diese Varianten mit dem neu entwickelten v4.1 GIAB-Benchmark-Set verglichen. Von den Varianten innerhalb der v4.1 Benchmark-Regionen stimmten die Ash1-Varianten mit 1.256.458 homozygoten und 1.041.476 heterozygoten SNPs sowie 187.227 homozygoten und 193.524 heterozygoten Indels überein. Nach Ausschluss von Variantenaufrufen innerhalb von 30 bp einer echten Variante blieben 79.269 SNPs und 17.439 Indels übrig, was (unter der Annahme, dass es sich dabei ausschließlich um Fehler in Ash1 handelt) einem Qualitätswert (QV) von etwa Q45 für Substitutionsfehler entspricht. Die meisten dieser Varianten (52.191 SNPs und 4629 Indels) fallen in segmentale Duplikationen, die möglicherweise fehlende Duplikationen in Ash1 oder unvollkommenes Polishing durch kurze Reads darstellen. Zusammenfassend lässt sich sagen, dass die Qualität des Ash1-Assemblies sehr hoch ist, mit einem geschätzten Substitutionsqualitätswert von 62 und einer Indel-Fehlerrate von 2 pro Million bp, nachdem bekannte segmentale Duplikationen, Tandemwiederholungen und Homopolymere ausgeschlossen wurden.

Vergleich des Variantenaufrufs unter Verwendung von Ash1 im Vergleich zu GRCh38

Eine der Motivationen für die Schaffung neuer Referenzgenome besteht darin, dass sie einen besseren Rahmen für die Analyse menschlicher Sequenzdaten bei der Suche nach genetischen Varianten im Zusammenhang mit Krankheiten bieten. Zum Beispiel sollte eine Studie über Ashkenazi-Juden, die Whole-Genome-Shotgun-Daten (WGS) gesammelt hat, ein Ashkenazi-Referenzgenom und nicht GRCh38 verwenden. Da der genetische Hintergrund ähnlich ist, sollten bei der Suche nach Ash1 weniger Varianten gefunden werden.

Um diese Erwartung zu testen, haben wir WGS-Daten eines männlichen Teilnehmers am Personal Genome Project, PGP17 (hu34D5B9), gesammelt. Dieses Individuum ist laut der PGP-Datenbank schätzungsweise zu 66 % aschkenasisch, was der höchste geschätzte Anteil war, der von bereits sequenzierten PGP-Individuen verfügbar war. Wir luden 983.220.918.100-bp-Reads herunter (ca. 33x Abdeckung) und richteten sie mit Bowtie2 sowohl an Ash1 als auch an GRCh38 aus. Ein etwas höherer Anteil der Reads (3.901.270, 0,5 %) wurde an Ash1 als an GRCh38 ausgerichtet.

Wir untersuchten dann alle Einzelnukleotidvarianten (SNVs, siehe „Methoden“) zwischen PGP17 und jedem der beiden Referenzgenome. Um die Analyse zu vereinfachen, haben wir nur die Stellen berücksichtigt, an denen PGP17 homozygot war. Bei unseren Vergleichen mit Ash1 identifizierten wir zunächst alle SNVs und untersuchten dann die ursprünglichen Ash1-Lesedaten, um festzustellen, ob das Ash1-Genom für jede dieser SNVs ein anderes Allel enthielt, das mit PGP17 übereinstimmte.

Insgesamt betrug die Anzahl der homozygoten Stellen in PGP17, die nicht mit Ash1 übereinstimmten, 1.333.345, gegenüber 1.700.364, als wir homozygote Stellen in PGP17 mit GRCh38 verglichen (Additional file 1: Table S1). Wir untersuchten dann die zugrundeliegenden Ash1-Lesedaten für die 1,33 Millionen SNV-Stellen, die anfänglich nicht übereinstimmten, und stellten fest, dass das Ash1-Genom bei etwa der Hälfte von ihnen heterozygot war, wobei ein Allel mit PGP17 übereinstimmte. Wenn wir die SNVs auf Stellen beschränkten, an denen PGP17 und Ash1 beide homozygot sind (plus eine sehr kleine Anzahl von Stellen, an denen Ash1 zwei Varianten enthält, die sich beide von PGP17 unterscheiden), reduzierte sich die Gesamtzahl der SNVs noch weiter auf 707.756. Wir fanden also knapp 1 Million weniger homozygote SNVs, wenn wir Ash1 als Referenz für PGP17 verwendeten. Anstelle eines Alignments an Ash1 könnte man die Reads auch an GRCh38 alignieren und dann bekannte populations-spezifische Varianten entfernen. Dieser zweistufige Prozess ist zwar komplexer, könnte aber zu ähnlichen Ergebnissen führen, mit Ausnahme von Regionen von Ash1, die in GRCh38 einfach fehlen.

Vergleich mit häufigen Ashkenazi-Varianten

Um zu prüfen, inwieweit Ash1 bekannte, häufige Ashkenazi-Varianten (im Vergleich zu GRCh38) enthält, haben wir SNVs untersucht, die mit hoher Häufigkeit in einer Ashkenazi-Population aus der Genome Aggregation Database (gnomAD) gemeldet wurden. GnomAD v3.0 enthält SNV-Calls aus Short-Read-Ganzgenomdaten von 1662 Ashkenazi-Individuen. Da einige Varianten nur bei einer Teilmenge dieser Personen genannt wurden, haben wir nur Varianten berücksichtigt, die bei mindestens 200 Personen gemeldet wurden. Anschließend sammelten wir SNVs mit Hauptallelen, wobei die Allelhäufigkeit in der untersuchten Population über 0,5 liegen musste. Außerdem beschränkten wir unsere Analyse auf Einzelbasensubstitutionen. Auf diese Weise erhielten wir 2.008.397 gnomAD-SNV-Stellen, an denen sich das Ashkenazi-Hauptallel von GRCh38 unterscheidet.

Wir konnten 1.790.688 der 2.008.397 gnomAD-Stellen von GRCh38 präzise auf Ash1 abbilden (siehe Abschnitt „Methoden“). Anschließend verglichen wir die GRCh38-Basis mit der Basis des Ashkenazi-Hauptallels an jeder Position, und wir untersuchten auch die alternativen Allele in Ash1 an Stellen, an denen es heterozygot ist. An Stellen, an denen die Lesedaten zeigten, dass HG002 heterozygot war und sowohl das Ashkenazi-Hauptallel als auch das GRCh38-Allel aufwies, ersetzten wir die Ash1-Base, falls erforderlich, um sicherzustellen, dass sie mit dem Hauptallel übereinstimmt. Nach diesen Ersetzungen enthielt Ash1 das Ashkenazi-Hauptallel an 88 % (1.580.866) der 1,79 Millionen Stellen. An den verbleibenden Stellen entsprach Ash1 entweder dem GRCh38-Allel, da HG002 homozygot für das Referenzallel ist (204.729 Stellen), oder es enthielt ein drittes Allel, das weder mit GRCh38 noch mit dem gnomAD-Hauptallel übereinstimmte (5093 Stellen). Diese Änderungen sollten die Anzahl der gemeldeten SNVs beim Alignment eines Ashkenazi-Individuums an Ash1 weiter reduzieren.

Bemerkenswert ist, dass mit zunehmender Häufigkeit des Hauptallels in der gnomAD-Ashkenazi-Population auch der Anteil der Stellen zunimmt, an denen Ash1 mit dem Hauptallel übereinstimmt. Beispielsweise sind von den SNVs, die in der gnomAD Ashkenazi-Population eine Allelhäufigkeit > 0,9 aufweisen, über 98 % in Ash1 vertreten (Tabelle 3).

Annotation

Ein wesentlicher Bestandteil eines jeden Referenzgenoms ist die Annotation: die Sammlung aller Gene und anderer Merkmale, die auf dem Genom gefunden wurden. Damit Ash1 als echtes Referenzgenom fungieren kann, haben wir alle bekannten Gene, die von der wissenschaftlichen Gemeinschaft verwendet werden, auf sein Koordinatensystem abgebildet und dieselben Gennamen und Bezeichnungen verwendet. Für GRCh38 wurden verschiedene Annotationsdatenbanken erstellt, und für Ash1 haben wir uns für die menschliche Gendatenbank CHESS entschieden, weil sie umfassend ist und alle proteinkodierenden Gene sowohl in GENCODE als auch in RefSeq, den beiden anderen weit verbreiteten Gendatenbanken, enthält und weil sie die in diesen Katalogen verwendeten Bezeichner beibehält. Die nicht kodierenden Gene sind in den drei Datenbanken unterschiedlich, aber CHESS enthält die größte Anzahl von Genorten und Isoformen. Wir haben CHESS Version 2.2 verwendet, die 42.167 Gene auf den primären Chromosomen enthält (ohne die alternativen Gerüste von GRCh38), von denen 20.197 für Proteine kodieren.

Die Zuordnung von Genen von einer Baugruppe zu einer anderen ist eine komplexe Aufgabe, insbesondere bei Genen, die in sehr ähnlichen Genfamilien mit mehreren Kopien auftreten. Die Aufgabe ist sogar noch komplexer, wenn die beiden Assemblies verschiedene Individuen repräsentieren (und nicht einfach verschiedene Assemblies desselben Individuums), da es Unterschiede zwischen den einzelnen Nukleotiden, Insertionen, Deletionen, Rearrangements und echte Variationen der Kopienzahl zwischen den Individuen gibt. Wir brauchten eine Methode, die angesichts all dieser potenziellen Unterschiede robust ist.

Um dieses Problem zu lösen, verwendeten wir das kürzlich entwickelte Liftoff-Mapping-Tool, das in unseren Experimenten das einzige Tool war, das fast alle menschlichen Gene von einem Individuum auf ein anderes übertragen konnte. Liftoff nimmt alle Gene und Transkripte aus einem Genom und ordnet sie Chromosom für Chromosom einem anderen Genom zu. Bei allen Genen, die sich nicht auf demselben Chromosom abbilden lassen, versucht Liftoff, sie chromosomenübergreifend abzubilden. Im Gegensatz zu anderen Tools stützt es sich nicht auf eine detaillierte Karte, die aus einem Alignment des gesamten Genoms erstellt wurde, sondern kartiert jedes Gen einzeln. Die Gene werden auf der Transkript-Ebene ausgerichtet, einschließlich der Introns, so dass prozessierte Pseudogene nicht fälschlicherweise als Gene identifiziert werden.

Wir haben versucht, alle 310.901 Transkripte von 42.167 Genorten auf den primären Chromosomen in GRCh38 auf Ash1 zu übertragen. Insgesamt konnten wir 309.900 (99,7 %) Transkripte von 42.070 Genloci auf den Hauptchromosomen erfolgreich kartieren (Zusatzdatei 1: Tabelle S2). Wir betrachteten ein Transkript als erfolgreich kartiert, wenn die mRNA-Sequenz in Ash1 mindestens 50% so lang ist wie die mRNA-Sequenz auf GRCh38. Die überwiegende Mehrheit der Transkripte liegt jedoch weit über diesem Schwellenwert, wobei 99 % der Transkripte mit einer Abdeckung von 95 % oder mehr kartiert wurden (Zusätzliche Datei 2: Abbildung S2). Die Sequenzidentität der kartierten Transkripte ist ähnlich hoch, wobei 99 % der Transkripte mit einer Sequenzidentität von mehr als oder gleich 94 % kartiert wurden (Zusatzdatei 2: Abbildung S3).

Translozierte Gene

Von den Genen mit mindestens einer erfolgreich kartierten Isoform wurden 42 059 (99,7 %) an den entsprechenden Stellen auf demselben Chromosom in Ash1 kartiert. Von den 108 Genen, die zunächst nicht kartiert werden konnten, wurden 11 Gene in 7 verschiedenen Blöcken auf einem anderen Chromosom kartiert (siehe Tabelle 4), was auf eine Translokation zwischen den beiden Genomen schließen lässt. Interessanterweise befanden sich 16 der 22 an den Translokationen beteiligten Stellen in subtelomeren Regionen, die in 8 Paaren auftraten, in denen beide Stellen in der Nähe von Telomeren lagen. Dies stimmt mit früheren Studien überein, in denen berichtet wurde, dass Umlagerungen, an denen Telomere und Subtelomere beteiligt sind, eine häufige Form der Translokation beim Menschen sein können.

Wir haben die Translokation zwischen den Chromosomen 15 und 20, die drei der Gene in Tabelle 4 enthält, genauer untersucht, indem wir uns das Alignment zwischen GRCh38 und Ash1 angesehen haben. Die Translokation befindet sich am Telomer beider Chromosomen, von Position 65.079.275 bis 65.109.824 (30.549 bp) von Ash1 chr20 und 101.950.338 bis 101.980.928 (30.590 bp) von GRCh39 chr15. Um die Translokation zu bestätigen, haben wir einen unabhängigen Satz sehr langer PacBio-Reads, alle von HG002, an die Ash1 v1.7-Assembly (siehe Abschnitt „Methoden“) angeglichen und die Region um den Bruchpunkt auf chr20 ausgewertet. Diese Alignments zeigen eine tiefe, konsistente Abdeckung, die sich über viele Kilobasen auf beiden Seiten des Bruchpunkts erstreckt, was die Korrektheit der Ash1-Assemblierung unterstützt (Abb. 1).

Schnappschuss, der Alignments von langen PacBio-Reads zum Ash1-Genom zeigt, zentriert auf das linke Ende der Stelle in Chromosom 20 (Position 65.079.275), an der eine Translokation zwischen Chromosom 15 (GRCh38) und 20 (Ash1) stattfand. Der obere Teil der Abbildung zeigt die Koordinaten auf Chromr20. Darunter befindet sich ein Histogramm der Leseabdeckung, und die einzelnen Reads füllen den unteren Teil der Abbildung. Die Indels in den Reads, die als farbige Balken auf jedem Read dargestellt sind, sind auf die relativ hohe Fehlerrate der langen Reads zurückzuführen

Fehlende Gene

Zweiundsechzig Gene konnten nicht vollständig von GRCh38 auf Ash1 abgebildet werden, und weitere 32 Gene konnten nur teilweise abgebildet werden (unter dem Schwellenwert von 50 % Abdeckung), wie in Tabelle 5 dargestellt. Alle Gene, die nicht oder nur teilweise kartiert werden konnten, waren Mitglieder von Multi-Gen-Familien, und in jedem Fall gab es mindestens eine weitere Kopie des fehlenden Gens in Ash1, mit einer durchschnittlichen Identität von 98,5 %. Es gibt also keine Fälle, in denen ein Gen, das in GRCh38 vorhanden ist, in Ash1 gänzlich fehlt; die in Tabelle 5 aufgeführten Gene repräsentieren Fälle, in denen Ash1 weniger Mitglieder einer Multi-Gen-Familie aufweist. Drei zusätzliche Gene (2 proteinkodierende, 1 lncRNA) wurden zwei nicht platzierten Contigs zugeordnet, was einen Anhaltspunkt für die Platzierung dieser Contigs in zukünftigen Versionen der Ash1-Assembly bieten wird.

Nach der Zuordnung der Gene zu Ash1 haben wir die kodierenden Sequenzen aus den Transkripten extrahiert, die vollständig zugeordnet wurden (Abdeckung gleich 100 %), sie an die kodierenden Sequenzen von GRCh38 angeglichen und die Varianten relativ zu GRCh38 genannt (siehe Abschnitt „Methoden“). Innerhalb der 35.513.365 bp in diesen proteinkodierenden Transkripten fanden wir 20.864 Einzelnukleotidvarianten und Indels. Vierzehntausendvierhundertneunundneunzig dieser Varianten fielen in die GIAB-„Callable“-Regionen für High-Confidence-Varianten, obwohl 3963 von ihnen in GIAB-„schwierigen“ repetitiven Regionen lagen, für die Alignments oft mehrdeutig sind. Von den 10.536 Varianten, die sich nicht in diesen schwierigen Regionen befanden, stimmten 10.456 (99,2 %) mit dem GIAB-Variantensatz mit hoher Zuverlässigkeit überein. In den schwierigen Regionen stimmten 3804/3963 (96,0 %) mit dem GIAB-Satz überein.

Wir haben dann die durch Varianten und unvollständiges Mapping verursachten Änderungen der Aminosäuren für alle proteinkodierenden Sequenzen annotiert. Von 124.238 proteincodierenden Transkripten aus 20.197 Genen haben 92.600 (74,5 %) 100 % identische Proteinsequenzen. Weitere 26 566 (21,4 %) weisen mindestens eine Aminosäureänderung auf, ergeben aber Proteine mit identischer Länge, und 1561 (1,3 %) weisen rahmenerhaltende Mutationen auf, bei denen eine oder mehrere Aminosäuren eingefügt oder entfernt werden, während der Rest des Proteins unverändert bleibt. Tabelle 6 enthält statistische Angaben zu allen Veränderungen in den Proteinsequenzen. Wenn ein Protein mehr als eine Variante aufwies, zählten wir es unter der folgenreichsten Variante, d. h. wenn ein Protein eine Missense-Variante und ein vorzeitiges Stoppcodon aufwies, zählten wir es zur Gruppe „Stop gained“.

Von besonderem Interesse sind jene Transkripte mit Varianten, die die Proteinsequenz erheblich stören und zu einem Funktionsverlust führen können. Dazu gehören Transkripte, die von einem Frameshift (2158), einem Stop-Loss (58), einem Stop-Gain (416), einem Start-Loss (58) oder einer Trunkierung aufgrund eines unvollständigen Mappings (564) betroffen sind. Diese gestörten Isoformen repräsentieren 885 Genorte; 505 dieser Gene haben jedoch mindestens eine weitere Isoform, die nicht von einer störenden Variante betroffen ist. Somit verbleiben 380 Gene, bei denen alle Isoformen mindestens eine Störung aufweisen; die vollständige Liste ist in Zusatzdatei 1: Tabelle S1 enthalten.