Assemblage et annotation d’un génome de référence humain ashkénaze

Pour la création du premier génome de référence humain assemblé à partir d’un seul individu, nous avons choisi HG002, un individu ashkénaze qui fait partie du Personal Genome Project (PGP). Le PGP utilise le modèle de consentement ouvert, la première plateforme à accès réellement ouvert pour le partage de données individuelles sur le génome humain, le phénotype et les données médicales. Le processus de consentement informe les participants potentiels des implications et des risques liés au partage des données génomiques, ainsi que de ce qu’ils peuvent attendre de leur participation. Le consentement ouvert a permis la création des premiers matériaux de référence du génome humain (HG002 est le matériau de référence NIST 8391) de Genome In A Bottle (GIAB), qui sont utilisés pour l’étalonnage, le développement de méthodes d’assemblage du génome et les mesures de performance des laboratoires. Toutes les données de séquence brutes pour ce projet ont été obtenues auprès de GIAB, où elles sont librement accessibles au public .

Nous avons assemblé le génome HG002 à partir d’une combinaison de trois ensembles de données à couverture profonde : des lectures Illumina de 249 pb, des lectures Oxford Nanopore (ONT) d’une longueur moyenne de plus de 33 kbp, et des lectures PacBio « HiFi » de haute qualité d’une longueur moyenne de 9567 pb (tableau 1).

Nous avons initialement créé deux assemblages : l’un utilisant les lectures Illumina et ONT, et le second utilisant les trois ensembles de données, y compris les lectures PacBio HiFi. L’ajout des données PacBio HiFi a donné lieu à une séquence totale légèrement plus importante dans l’assemblage (2,99 Gb contre 2,88 Gb) et à une taille N50 de contigs sensiblement plus grande (18,2 Mb contre 4,9 Mb). Cet assemblage, désigné Ash1 v0.5, a servi de base à tous les raffinements ultérieurs.

Mappage de l’assemblage sur des chromosomes

Pour créer des affectations de chromosomes pour l’assemblage Ash1 v0.5, nous avons utilisé des alignements sur GRCh38 pour mapper la plupart des échafaudages sur des chromosomes. Les étapes décrites dans la section « Méthodes » ont généré une série d’assemblages à l’échelle du chromosome progressivement améliorés, aboutissant à Ash1 v1.7. Ash1 v1.7 présente une plus grande contiguïté et des lacunes moins importantes que GRCh38, comme le montre le tableau 2. Il faut noter qu’au cours du processus de construction de ces chromosomes, une petite quantité de la séquence GRCh38 (58,3 Mb, 2 % du génome) a été utilisée pour combler les lacunes de Ash1. Ces régions comprennent certaines régions difficiles à assembler qui ont été manuellement conservées pour GRCh38. Au total, la taille estimée de toutes les lacunes dans Ash1 est de 82,9 Mbp, contre 84,7 Mbp dans GRCh38.p13.

Dans le cadre du processus d’amélioration de l’assemblage, nous avons recherché dans l’un des assemblages Ash1 préliminaires (v1.1) les 12 745 variantes structurelles isolées de haute qualité (insertions et délétions ≥ 50 pb) que Zook et al. ont identifiées en comparant les données du trio ashkénaze à GRCh37 . Cette étude a utilisé quatre technologies de séquençage différentes et plusieurs appelants de variants pour identifier les variants et filtrer les faux positifs. Sur ces 12 745 VS, 5807 sont homozygotes et 6938 sont hétérozygotes. Nous nous attendions à ce que l’assemblage Ash1 soit en accord avec presque tous les variants homozygotes. Comme Ash1 ne capture qu’un seul haplotype, nous nous attendions à ce qu’il soit en accord avec environ la moitié des VS hétérozygotes, en supposant que l’algorithme d’assemblage choisisse au hasard entre les haplotypes lorsqu’il décide de la variante à inclure dans le consensus final. Des 5807 variants homozygotes, 5284 (91%) étaient présents en utilisant nos critères de correspondance (voir la section « Méthodes »), et 3922 (56,5%) des 6938 variants hétérozygotes étaient présents. Tous les variants ont été trouvés à l’emplacement correct.

Nous avons également effectué des appels de petits variants (≤ 5 pb) sur Ash1 v1.1 et les avons comparés aux variants de référence HG002 v4.0 du GIAB, que nous avons utilisés pour corriger de nombreuses erreurs de substitution et d’indel (voir la section « Méthodes »), ce qui a donné Ash 1 v1.2. Nous avons ensuite réaligné l’assemblage Ash1 sur GRCh38, réappelé les variants, et comparé ces variants à l’ensemble de référence v4.1 du GIAB nouvellement développé. Parmi les variants à l’intérieur des régions de référence v4.1, les variants Ash1 correspondaient à 1 256 458 SNP homozygotes et 1 041 476 SNP hétérozygotes, et à 187 227 indels homozygotes et 193 524 indels hétérozygotes. Après avoir exclu les appels de variants à moins de 30 pb d’un vrai variant, il restait 79 269 SNP et 17 439 indels, ce qui (en supposant qu’il s’agit de toutes les erreurs dans Ash1) correspond à une valeur de qualité (QV) d’environ Q45 pour les erreurs de substitution. La plupart de ces variantes (52 191 SNP et 4629 indels) tombent dans des duplications segmentaires, représentant peut-être des duplications manquantes dans Ash1 ou un polissage imparfait par des lectures courtes. En résumé, la qualité de l’assemblage Ash1 est très élevée, avec une valeur de qualité de substitution estimée à 62 et un taux d’erreur d’indel de 2 par million de pb après exclusion des duplications segmentaires connues, des répétitions en tandem et des homopolymères.

Comparaison de l’appel de variants en utilisant Ash1 par rapport à GRCh38

L’une des motivations pour la création de nouveaux génomes de référence est qu’ils fournissent un meilleur cadre pour l’analyse des données de séquences humaines lors de la recherche de variants génétiques liés à des maladies. Par exemple, une étude sur les juifs ashkénazes qui a recueilli des données sur le génome entier (WGS) devrait utiliser un génome de référence ashkénaze plutôt que GRCh38. Comme le fond génétique est similaire, moins de variantes devraient être trouvées lors d’une recherche contre Ash1.

Pour tester cette attente, nous avons recueilli des données WGS d’un participant masculin au Personal Genome Project, PGP17 (hu34D5B9). Cet individu est estimé à 66% ashkénaze selon la base de données PGP, qui était la fraction estimée la plus élevée disponible à partir des individus PGP déjà séquencés. Nous avons téléchargé 983 220 918 lectures de 100 pb (couverture d’environ 33x) et les avons alignées sur Ash1 et GRCh38 en utilisant Bowtie2. Une fraction légèrement plus élevée de lectures (3 901 270, 0,5 %) s’est alignée sur Ash1 que sur GRCh38.

Nous avons ensuite examiné toutes les variantes mononucléotidiques (SNV, voir les « Méthodes ») entre PGP17 et chacun des deux génomes de référence. Pour simplifier l’analyse, nous n’avons pris en compte que les emplacements où PGP17 était homozygote. Dans nos comparaisons avec Ash1, nous avons d’abord identifié tous les SNV, puis examiné les données de lecture originales d’Ash1 afin de déterminer si, pour chacun de ces SNV, le génome d’Ash1 contenait un allèle différent qui correspondait à PGP17.

Au total, le nombre de sites homozygotes de PGP17 en désaccord avec Ash1 était de 1 333 345, contre 1 700 364 lorsque nous avons comparé les sites homozygotes de PGP17 à GRCh38 (fichier supplémentaire 1 : tableau S1). Nous avons ensuite examiné les données de lecture sous-jacentes de l’Ash1 pour les 1,33 million de sites SNV qui ne correspondaient pas initialement, et nous avons constaté que pour environ la moitié d’entre eux, le génome de l’Ash1 était hétérozygote, avec un allèle correspondant à PGP17. Si nous limitons les SNV aux sites où PGP17 et Ash1 sont tous deux homozygotes (plus un très petit nombre d’emplacements où Ash1 contient deux variants qui diffèrent tous deux de PGP17), cela réduit encore le nombre total de SNV à 707 756. Ainsi, nous avons trouvé un peu moins d’un million de SNV homozygotes en moins en utilisant Ash1 comme référence pour PGP17. Il est à noter qu’au lieu d’aligner les lectures sur Ash1, on pourrait les aligner sur GRCh38, puis éliminer les variants connus spécifiques à la population. Ce processus en deux étapes, bien que plus complexe, pourrait donner des résultats similaires, sauf pour les régions de Ash1 qui sont tout simplement absentes de GRCh38.

Comparaison avec les variants ashkénazes communs

Pour examiner dans quelle mesure Ash1 contient des variants ashkénazes connus et communs (par rapport à GRCh38), nous avons examiné les SNV signalés à haute fréquence dans une population ashkénaze à partir de la base de données d’agrégation du génome (gnomAD) . GnomAD v3.0 contient des appels de SNV à partir de données de génome entier à lecture courte provenant de 1662 individus ashkénazes. Comme certains variants n’ont été appelés que dans un sous-ensemble de ces individus, nous n’avons pris en compte que les sites de variants qui ont été signalés chez un minimum de 200 personnes. Nous avons ensuite recueilli les SNV d’allèles majeurs, en exigeant que la fréquence de l’allèle soit supérieure à 0,5 dans la population échantillonnée. Nous avons ensuite restreint notre analyse aux substitutions d’une seule base. Cela nous a donné 2 008 397 sites SNV gnomAD où l’allèle majeur ashkénaze différait de celui de GRCh38.

Nous avons pu cartographier précisément 1 790 688 des 2 008 397 sites gnomAD de GRCh38 sur Ash1 (voir la section « Méthodes »). Nous avons ensuite comparé la base de GRCh38 à la base de l’allèle majeur Ashkénaze à chaque position, et nous avons également examiné les allèles alternatifs dans Ash1 aux sites où il est hétérozygote. Pour les sites où les données de lecture ont montré que HG002 était hétérozygote et possédait à la fois l’allèle majeur ashkénaze et l’allèle GRCh38, nous avons remplacé la base de l’Ash1, si nécessaire, pour nous assurer qu’elle correspondait à l’allèle majeur. Après ces remplacements, Ash1 contenait l’allèle majeur ashkénaze sur 88 % (1 580 866) des 1,79 million de sites. Sur les sites restants, Ash1 correspondait soit à l’allèle GRCh38 parce que HG002 est homozygote pour l’allèle de référence (204 729 sites), soit il contenait un troisième allèle ne correspondant ni à GRCh38 ni à l’allèle majeur gnomAD (5093 sites). Ces modifications devraient réduire davantage le nombre de SNV signalés lors de l’alignement d’un individu ashkénaze sur Ash1.

Il convient de noter que, à mesure que la fréquence de l’allèle majeur dans la population ashkénaze gnomAD augmente, la proportion de sites où Ash1 correspond à l’allèle majeur augmente également. Par exemple, parmi les SNV qui ont une fréquence d’allèle > 0,9 dans la population ashkénaze gnomAD, plus de 98% sont représentés dans Ash1 (tableau 3).

Annotation

Une partie vitale de tout génome de référence est l’annotation : la collection de tous les gènes et autres caractéristiques trouvés sur le génome. Pour permettre à Ash1 de fonctionner comme un véritable génome de référence, nous avons cartographié tous les gènes connus utilisés par la communauté scientifique sur son système de coordonnées, en utilisant les mêmes noms et identifiants de gènes. Plusieurs bases de données d’annotation différentes ont été créées pour GRCh38, et pour Ash1, nous avons choisi d’utiliser la base de données de gènes humains CHESS parce qu’elle est complète, incluant tous les gènes codant pour des protéines à la fois dans GENCODE et RefSeq, les deux autres bases de données de gènes largement utilisées, et parce qu’elle conserve les identifiants utilisés dans ces catalogues. Les gènes non codants diffèrent entre les trois bases de données, mais CHESS possède le plus grand nombre de loci et d’isoformes de gènes. Nous avons utilisé la version 2.2 de CHESS, qui compte 42 167 gènes sur les chromosomes primaires (à l’exclusion des échafaudages alternatifs GRCh38), dont 20 197 codent pour des protéines.

Mapper les gènes d’un assemblage à un autre est une tâche complexe, en particulier pour les gènes qui apparaissent dans des familles de gènes multicopies très similaires. La tâche est encore plus complexe lorsque les deux assemblages représentent des individus différents (plutôt que simplement des assemblages différents du même individu), en raison de la présence de différences mononucléotidiques, d’insertions, de délétions, de réarrangements et de véritables variations du nombre de copies entre les individus. Nous avions besoin d’une méthode qui soit robuste face à toutes ces différences potentielles.

Pour résoudre ce problème, nous avons utilisé l’outil de cartographie Liftoff récemment développé, qui, dans nos expériences, était le seul outil capable de cartographier presque tous les gènes humains d’un individu à l’autre. Liftoff prend tous les gènes et les transcrits d’un génome et les fait correspondre, chromosome par chromosome, à un autre génome. Pour tous les gènes qui ne parviennent pas à être cartographiés sur le même chromosome, Liftoff tente de les cartographier sur plusieurs chromosomes. Contrairement à d’autres outils, il ne s’appuie pas sur une carte détaillée construite à partir d’un alignement du génome entier, mais cartographie chaque gène individuellement. Les gènes sont alignés au niveau des transcriptions, y compris les introns, de sorte que les pseudogènes transformés ne seront pas identifiés par erreur comme des gènes.

Nous avons tenté de cartographier les 310 901 transcriptions de 42 167 loci de gènes sur les chromosomes primaires dans GRCh38 à Ash1. Au total, nous avons réussi à cartographier 309 900 (99,7 %) transcrits de 42 070 loci génétiques sur les chromosomes principaux (fichier supplémentaire 1 : tableau S2). Nous avons considéré qu’une transcription était cartographiée avec succès si la séquence d’ARNm dans Ash1 est au moins 50% plus longue que la séquence d’ARNm sur GRCh38. Cependant, la grande majorité des transcrits dépassent largement ce seuil, 99 % des transcrits étant cartographiés avec une couverture supérieure ou égale à 95 % (fichier supplémentaire 2 : figure S2). L’identité de séquence des transcrits cartographiés est également élevée, avec 99 % des transcrits cartographiés avec une identité de séquence supérieure ou égale à 94 % (fichier supplémentaire 2 : figure S3).

Gènes translocalisés

Parmi les gènes dont au moins une isoforme a été cartographiée avec succès, 42 059 (99,7 %) ont été cartographiés aux emplacements correspondants sur le même chromosome dans Ash1. Parmi les 108 gènes qui n’ont pas réussi à s’associer, 11 gènes se sont associés à un chromosome différent dans 7 blocs distincts (indiqués dans le tableau 4), ce qui suggère une translocation entre les deux génomes. Il est intéressant de noter que 16 des 22 emplacements impliqués dans les translocations se trouvaient dans des régions subtélomériques, ce qui s’est produit dans 8 paires où les deux emplacements étaient proches des télomères. Ceci est cohérent avec des études précédentes rapportant que les réarrangements impliquant les télomères et les subtélomères peuvent être une forme commune de translocation chez les humains .

Nous avons examiné la translocation entre les chromosomes 15 et 20, qui contient trois des gènes du tableau 4, en regardant de plus près l’alignement entre GRCh38 et Ash1. La translocation se trouve au niveau du télomère des deux chromosomes, de la position 65 079 275 à 65 109 824 (30 549 pb) du chr20 de Ash1 et 101 950 338 à 101 980 928 (30 590 pb) du chr15 de GRCh39. Pour confirmer la translocation, nous avons aligné un ensemble indépendant de lectures PacBio très longues, toutes provenant de HG002, sur l’assemblage Ash1 v1.7 (voir la section « Méthodes ») et évalué la région autour du point de rupture sur le chr20. Ces alignements montrent une couverture profonde et cohérente s’étendant sur plusieurs kilobases de part et d’autre du point de rupture, soutenant l’exactitude de l’assemblage Ash1 (Fig. 1).

Snapshot montrant les alignements de longues lectures PacBio au génome Ash1, centré sur l’extrémité gauche de l’emplacement dans le chromosome 20 (position 65,079,275) où une translocation s’est produite entre le chromosome 15 (GRCh38) et 20 (Ash1). La partie supérieure de la figure montre les coordonnées sur le chr20. En dessous se trouve un histogramme de la couverture des lectures, et les lectures individuelles remplissent la partie inférieure de la figure. Les indels dans les lectures, représentés par des barres colorées sur chaque lecture, sont dus au taux d’erreur relativement élevé des longues lectures

Gènes manquants

Sixante-deux gènes n’ont pas réussi à se mapper entièrement de GRCh38 sur Ash1, et 32 autres gènes ne se sont mappés que partiellement (sous le seuil de couverture de 50%), comme le montre le tableau 5. Tous les gènes qui n’ont pas pu être cartographiés ou qui l’ont été partiellement étaient membres de familles multigéniques et, dans tous les cas, il y avait au moins une autre copie du gène manquant présent dans Ash1, avec une identité moyenne de 98,5 %. Il n’y a donc aucun cas où un gène présent dans GRCh38 est entièrement absent de Ash1 ; les gènes indiqués dans le tableau 5 représentent des cas où Ash1 a moins de membres d’une famille multigénique. Trois gènes supplémentaires (2 codant pour des protéines, 1 lncRNA) ont cartographié à deux contigs non placés, ce qui fournira un guide pour placer ces contigs dans les futures versions de l’assemblage Ash1.

Après avoir cartographié les gènes sur Ash1, nous avons extrait les séquences codantes des transcrits qui ont été cartographiés complètement (couverture égale à 100%), nous les avons alignés sur les séquences codantes de GRCh38, et nous avons appelé les variants par rapport à GRCh38 (voir la section « Méthodes »). Dans les 35 513 365 pb de ces transcriptions codant pour des protéines, nous avons trouvé 20 864 variants et indels mononucléotidiques. Quatorze mille quatre cent quatre-vingt-dix-neuf de ces variants se trouvaient dans les régions  » appelables  » du GIAB pour les variants à haute confiance, bien que 3963 d’entre eux se trouvaient dans les régions répétitives  » difficiles  » du GIAB, pour lesquelles les alignements sont souvent ambigus. Sur les 10 536 variants qui ne se trouvaient pas dans ces régions difficiles, 10 456 (99,2 %) étaient en accord avec l’ensemble de variants de haute confiance de la GIAB. Dans les régions difficiles, 3804/3963 (96,0%) étaient en accord avec l’ensemble GIAB.

Nous avons ensuite annoté les changements d’acides aminés causés par les variants et les cartographies incomplètes pour toutes les séquences codant pour des protéines. Sur les 124 238 transcrits codant pour des protéines provenant de 20 197 gènes, 92 600 (74,5%) ont des séquences protéiques 100% identiques. Par ailleurs, 26 566 (21,4 %) présentent au moins un changement d’acide aminé mais produisent des protéines de longueur identique, et 1 561 (1,3 %) présentent des mutations préservant le cadre qui insèrent ou suppriment un ou plusieurs acides aminés, laissant le reste de la protéine inchangé. Le tableau 6 présente des statistiques sur toutes les modifications des séquences de protéines. Si une protéine avait plus d’une variante, nous l’avons comptée sous la variante la plus conséquente, c’est-à-dire que si une protéine avait une variante faux-sens et un codon stop prématuré, nous l’incluons dans le groupe « stop gagné ».

Les transcrits présentant des variantes qui perturbent de manière significative la séquence protéique et peuvent entraîner une perte de fonction sont particulièrement intéressants. Il s’agit notamment des transcrits affectés par un décalage de cadre (2158), une perte d’arrêt (58), un gain d’arrêt (416), une perte de début (58) ou une troncature due à une cartographie incomplète (564). Ces isoformes perturbées représentent 885 loci génétiques ; cependant, 505 de ces gènes ont au moins une autre isoforme qui n’est pas affectée par une variante perturbatrice. Cela laisse 380 gènes dans lesquels toutes les isoformes ont au moins une perturbation ; la liste complète est fournie dans le fichier additionnel 1 : Tableau S1.