Asamblaje y anotación de un genoma humano de referencia asquenazí

Para la creación del primer genoma humano de referencia ensamblado a partir de un solo individuo, elegimos a HG002, un individuo asquenazí que forma parte del Proyecto Genoma Personal (PGP). El PGP utiliza el Modelo de Consentimiento Abierto, la primera plataforma de acceso verdaderamente abierto para compartir el genoma humano individual, el fenotipo y los datos médicos. El proceso de consentimiento educa a los participantes potenciales sobre las implicaciones y los riesgos de compartir datos genómicos, y sobre lo que pueden esperar de su participación. El consentimiento abierto ha permitido la creación de los primeros materiales de referencia del mundo sobre el genoma humano (el HG002 es el material de referencia 8391 del NIST) de Genome In A Bottle (GIAB), que se está utilizando para la calibración, el desarrollo de métodos de ensamblaje del genoma y la medición del rendimiento del laboratorio. Todos los datos de las secuencias en bruto para este proyecto se obtuvieron de GIAB, donde están disponibles de forma gratuita para el público.

Ensamblamos el genoma HG002 a partir de una combinación de tres conjuntos de datos de cobertura profunda: lecturas de Illumina de 249 pb, lecturas de Oxford Nanopore (ONT) con una longitud media de más de 33 kbp y lecturas de PacBio «HiFi» de alta calidad con una media de 9567 pb (Tabla 1).

Inicialmente creamos dos ensamblajes: uno utilizando lecturas de Illumina y ONT, y un segundo utilizando los tres conjuntos de datos, incluyendo las lecturas de PacBio HiFi. La adición de los datos de PacBio HiFi dio como resultado una secuencia total ligeramente mayor en el ensamblaje (2,99 Gb frente a 2,88 Gb) y un tamaño N50 de contig sustancialmente mayor (18,2 Mb frente a 4,9 Mb). Este ensamblaje, designado Ash1 v0.5, fue la base para todos los refinamientos posteriores.

Mapeo del ensamblaje en los cromosomas

Para crear asignaciones de cromosomas para el ensamblaje Ash1 v0.5, utilizamos alineaciones con GRCh38 para mapear la mayoría de los andamios en los cromosomas. Los pasos descritos en la sección «Métodos» generaron una serie de ensamblajes a escala cromosómica gradualmente mejorados, dando como resultado Ash1 v1.7. Ash1 v1.7 tiene mayor contigüidad y menores huecos que GRCh38, como se muestra en la Tabla 2. Obsérvese que en el proceso de construcción de estos cromosomas se utilizó una pequeña cantidad de secuencia de GRCh38 (58,3 Mb, el 2% del genoma) para rellenar huecos en Ash1. Estas regiones incluyen algunas regiones difíciles de ensamblar que han sido curadas manualmente para GRCh38. En total, el tamaño estimado de todos los huecos en Ash1 es de 82,9 Mbp, frente a 84,7 Mbp en GRCh38.p13.

Como parte del proceso de mejora del ensamblaje, buscamos en uno de los ensamblajes preliminares de Ash1 (v1.1) las 12.745 variantes estructurales aisladas de alta calidad (inserciones y deleciones ≥ 50 pb) que Zook et al. identificaron comparando los datos del trío Ashkenazi con GRCh37. Ese estudio utilizó cuatro tecnologías de secuenciación diferentes y múltiples llamadores de variantes para identificar las variantes y filtrar los falsos positivos. De estos 12.745 SVs, 5807 son homocigotos y 6938 son heterocigotos. Esperábamos que el ensamblaje de Ash1 coincidiera con casi todas las variantes homocigóticas. Dado que Ash1 captura sólo un haplotipo, esperábamos que coincidiera con aproximadamente la mitad de los SV heterocigotos, suponiendo que el algoritmo de ensamblaje eligiera al azar entre los haplotipos al decidir qué variante incluir en el consenso final. De las 5807 variantes homocigóticas, 5284 (91%) estaban presentes utilizando nuestros criterios de coincidencia (véase la sección «Métodos»), y 3922 (56,5%) de las 6938 variantes heterocigóticas estaban presentes. Todas las variantes se encontraron en la localización correcta.

También realizamos llamadas de variantes pequeñas (≤ 5 pb) en Ash1 v1.1 y las comparamos con las variantes de referencia HG002 v4.0 de GIAB, que utilizamos para corregir numerosos errores de sustitución e indel (véase la sección «Métodos»), obteniendo Ash 1 v1.2. A continuación, volvimos a alinear el ensamblaje de Ash1 con GRCh38, volvimos a llamar a las variantes y las comparamos con el nuevo conjunto de referencia v4.1 de GIAB. De las variantes dentro de las regiones de referencia v4.1, las variantes de Ash1 coincidieron con 1.256.458 SNPs homocigóticos y 1.041.476 heterocigóticos, y 187.227 indels homocigóticos y 193.524 heterocigóticos. Tras excluir las llamadas de variantes a menos de 30 pb de una variante verdadera, quedaron 79.269 SNP y 17.439 indels, lo que (asumiendo que todos son errores en Ash1) corresponde a un valor de calidad (QV) de aproximadamente Q45 para los errores de sustitución. La mayoría de estas variantes (52.191 SNPs y 4629 indels) caen en duplicaciones segmentarias, posiblemente representando duplicaciones faltantes en Ash1 o un pulido imperfecto por lecturas cortas. En resumen, la calidad del ensamblaje de Ash1 es muy alta, con un valor estimado de calidad de sustitución de 62 y una tasa de error de indel de 2 por millón de pb después de excluir las duplicaciones segmentarias conocidas, las repeticiones en tándem y los homopolímeros.

Comparación de la llamada de variantes utilizando Ash1 frente a GRCh38

Una de las motivaciones para crear nuevos genomas de referencia es que proporcionan un mejor marco para analizar los datos de secuencias humanas cuando se buscan variantes genéticas vinculadas a enfermedades. Por ejemplo, un estudio sobre los judíos asquenazíes que recogiera datos de escopeta de genoma completo (WGS) debería utilizar un genoma de referencia asquenazí en lugar de GRCh38. Debido a que los antecedentes genéticos son similares, se deberían encontrar menos variantes cuando se busque contra Ash1.

Para probar esta expectativa, recogimos datos WGS de un participante masculino en el Proyecto Genoma Personal, PGP17 (hu34D5B9). Se estima que este individuo es un 66% asquenazí según la base de datos del PGP, que era la fracción estimada más alta disponible de los individuos del PGP ya secuenciados. Descargamos 983.220.918.100 pb de lecturas (aproximadamente 33 veces la cobertura) y las alineamos con Ash1 y GRCh38 utilizando Bowtie2 . Una fracción ligeramente mayor de lecturas (3.901.270, 0,5%) se alineó con Ash1 que con GRCh38.

A continuación, examinamos todas las variantes de un solo nucleótido (SNVs, ver «Métodos») entre PGP17 y cada uno de los dos genomas de referencia. Para simplificar el análisis, sólo consideramos las localizaciones en las que PGP17 era homocigota. En nuestras comparaciones con Ash1, primero identificamos todos los SNV y luego examinamos los datos de lectura originales de Ash1 para determinar si, para cada uno de esos SNV, el genoma de Ash1 contenía un alelo diferente que coincidía con el de PGP17.

En total, el número de sitios homocigóticos en PGP17 que no coincidían con Ash1 era de 1.333.345, frente a 1.700.364 cuando comparamos los sitios homocigóticos en PGP17 con GRCh38 (Archivo adicional 1: Tabla S1). A continuación, examinamos los datos de lectura subyacentes de Ash1 para los 1,33 millones de sitios SNV que inicialmente no coincidían, y descubrimos que para aproximadamente la mitad de ellos, el genoma de Ash1 era heterocigoto, con un alelo que coincidía con PGP17. Si restringimos los SNVs a los sitios en los que PGP17 y Ash1 son ambos homocigotos (además de un número muy pequeño de lugares en los que Ash1 contiene dos variantes que difieren ambas de PGP17), esto redujo el número total de SNVs aún más, a 707.756. Por lo tanto, encontramos algo menos de 1 millón de SNVs homocigotos al utilizar Ash1 como referencia para PGP17. Nótese que en lugar de alinear con Ash1, se podrían alinear las lecturas con GRCh38 y luego eliminar las variantes específicas de la población conocida. Este proceso de dos pasos, aunque más complejo, podría producir resultados similares, excepto para las regiones de Ash1 que simplemente faltan en GRCh38.

Comparación con variantes comunes asquenazíes

Para examinar hasta qué punto Ash1 contiene variantes asquenazíes conocidas y comunes (en relación con GRCh38), examinamos los SNVs reportados con alta frecuencia en una población asquenazí de la Base de Datos de Agregación del Genoma (gnomAD) . GnomAD v3.0 contiene llamadas de SNV de datos de genoma completo de lectura corta de 1662 individuos asquenazíes. Debido a que algunas variantes sólo fueron llamadas en un subconjunto de estos individuos, consideramos sólo los sitios de variantes que fueron reportados en un mínimo de 200 personas. A continuación, recopilamos las SNV de los alelos principales, requiriendo que la frecuencia del alelo fuera superior a 0,5 en la población muestreada. Además, restringimos nuestro análisis a las sustituciones de una sola base. Esto nos proporcionó 2.008.397 sitios SNV gnomAD en los que el alelo mayor asquenazí difería de GRCh38.

Pudimos mapear con precisión 1.790.688 de los 2.008.397 sitios gnomAD de GRCh38 en Ash1 (véase la sección «Métodos»). A continuación, comparamos la base de GRCh38 con la base del alelo mayor de Ashkenazi en cada posición, y también examinamos los alelos alternativos de Ash1 en los sitios donde es heterocigoto. En los sitios en los que los datos de lectura mostraban que HG002 era heterocigota y tenía tanto el alelo mayor asquenazí como el alelo GRCh38, sustituimos la base de Ash1, si era necesario, para asegurarnos de que coincidía con el alelo mayor. Tras estas sustituciones, Ash1 contenía el alelo mayor asquenazí en el 88% (1.580.866) de los 1,79 millones de sitios. En los sitios restantes, Ash1 coincidía con el alelo GRCh38 porque HG002 es homocigoto para el alelo de referencia (204.729 sitios), o contenía un tercer alelo que no coincidía ni con GRCh38 ni con el alelo mayor gnomAD (5093 sitios). Estas modificaciones deberían reducir aún más el número de SNVs reportados cuando se alinea un individuo Ashkenazi con Ash1.

Digno de mención es que, a medida que la frecuencia del alelo mayor en la población Ashkenazi de gnomAD aumenta, la proporción de sitios donde Ash1 coincidió con el alelo mayor también aumenta. Por ejemplo, de los SNVs que tienen una frecuencia alélica > 0.9 en la población Ashkenazi gnomAD, más del 98% están representados en Ash1 (Tabla 3).

Anotación

Una parte vital de cualquier genoma de referencia es la anotación: la colección de todos los genes y otras características encontradas en el genoma. Para que Ash1 funcione como un verdadero genoma de referencia, hemos mapeado todos los genes conocidos utilizados por la comunidad científica en su sistema de coordenadas, utilizando los mismos nombres e identificadores de genes. Se han creado varias bases de datos de anotaciones diferentes para GRCh38, y para Ash1, elegimos utilizar la base de datos de genes humanos CHESS porque es exhaustiva, ya que incluye todos los genes codificadores de proteínas tanto en GENCODE como en RefSeq, las otras dos bases de datos de genes ampliamente utilizadas, y porque conserva los identificadores utilizados en esos catálogos. Los genes no codificantes difieren entre las tres bases de datos, pero CHESS es la que tiene el mayor número de loci génicos e isoformas. Usamos la versión 2.2 de CHESS, que tiene 42.167 genes en los cromosomas primarios (excluyendo los andamios alternativos de GRCh38), de los cuales 20.197 son codificantes de proteínas.

Mapear los genes de un ensamblaje a otro es una tarea compleja, particularmente para los genes que ocurren en familias de genes altamente similares y multicopiados. La tarea es aún más compleja cuando los dos conjuntos representan a diferentes individuos (en lugar de simplemente diferentes conjuntos del mismo individuo), debido a la presencia de diferencias de un solo nucleótido, inserciones, deleciones, reordenamientos y variaciones genuinas en el número de copias entre los individuos. Necesitábamos un método que fuera robusto frente a todas estas diferencias potenciales.

Para abordar este problema, utilizamos la herramienta de mapeo Liftoff, recientemente desarrollada, que en nuestros experimentos fue la única herramienta que pudo mapear casi todos los genes humanos de un individuo a otro. Liftoff toma todos los genes y transcripciones de un genoma y los mapea, cromosoma a cromosoma, a un genoma diferente. En el caso de los genes que no se asignan al mismo cromosoma, Liftoff intenta asignarlos a otros cromosomas. A diferencia de otras herramientas, no se basa en un mapa detallado construido a partir de un alineamiento de todo el genoma, sino que mapea cada gen individualmente. Los genes se alinean a nivel de transcripción, incluyendo los intrones, para que los pseudogenes procesados no se identifiquen erróneamente como genes.

Intentamos mapear los 310.901 transcritos de 42.167 loci de genes en los cromosomas primarios de GRCh38 a Ash1. En total, mapeamos con éxito 309.900 (99,7%) transcritos de 42.070 loci génicos en los cromosomas principales (Archivo adicional 1: Tabla S2). Consideramos que un transcrito se ha mapeado con éxito si la secuencia de ARNm en Ash1 es al menos un 50% más larga que la secuencia de ARNm en GRCh38. Sin embargo, la gran mayoría de los transcritos superan ampliamente este umbral, con el 99% de los transcritos mapeados con una cobertura mayor o igual al 95% (Archivo adicional 2: Figura S2). La identidad de secuencia de los transcritos mapeados es igualmente alta, con un 99% de transcritos mapeados con una identidad de secuencia mayor o igual al 94% (Archivo adicional 2: Figura S3).

Genes translocalizados

De aquellos genes con al menos una isoforma mapeada con éxito, 42.059 (99,7%) se mapearon en las localizaciones correspondientes del mismo cromosoma en Ash1. De los 108 genes que inicialmente no se mapearon, 11 genes se mapearon en un cromosoma diferente en 7 bloques distintos (mostrados en la Tabla 4), sugiriendo una translocación entre los dos genomas. Curiosamente, 16 de las 22 localizaciones implicadas en las translocaciones se encontraban en regiones subteloméricas, lo que ocurrió en 8 pares en los que ambas localizaciones estaban cerca de los telómeros. Esto es coherente con estudios anteriores que informan de que los reordenamientos que implican telómeros y subtelómeros pueden ser una forma común de translocación en los seres humanos.

Examinamos la translocación entre los cromosomas 15 y 20, que contiene tres de los genes de la Tabla 4, observando más de cerca la alineación entre GRCh38 y Ash1. La translocación se encuentra en el telómero de ambos cromosomas, desde la posición 65.079.275 a 65.109.824 (30.549 pb) de Ash1 chr20 y 101.950.338 a 101.980.928 (30.590 pb) de GRCh39 chr15. Para confirmar la translocación, alineamos un conjunto independiente de lecturas muy largas de PacBio, todas procedentes de HG002, con el ensamblaje de Ash1 v1.7 (véase la sección «Métodos») y evaluamos la región alrededor del punto de rotura en la cr20. Estos alineamientos muestran una cobertura profunda y consistente que se extiende muchos kilobases a ambos lados del punto de ruptura, apoyando la corrección del ensamblaje de Ash1 (Fig. 1).

Imagen que muestra las alineaciones de las lecturas largas de PacBio con el genoma de Ash1, centradas en el extremo izquierdo de la localización en el cromosoma 20 (posición 65.079.275) donde se produjo una translocación entre el cromosoma 15 (GRCh38) y el 20 (Ash1). La parte superior de la figura muestra las coordenadas en el cromosoma 20. Debajo hay un histograma de la cobertura de las lecturas, y las lecturas individuales llenan la parte inferior de la figura. Los indels en las lecturas, mostrados como barras de color en cada lectura, se deben a la tasa de error relativamente alta de las lecturas largas

Genes perdidos

Sesenta y dos genes no pudieron mapearse completamente desde GRCh38 a Ash1, y otros 32 genes se mapearon sólo parcialmente (por debajo del umbral de cobertura del 50%), como se muestra en la Tabla 5. Todos los genes que no se mapearon o que se mapearon parcialmente eran miembros de familias multigénicas, y en todos los casos, había al menos otra copia del gen que faltaba presente en Ash1, con una identidad media del 98,5%. Por lo tanto, no hay ningún caso de un gen que esté presente en GRCh38 y que esté totalmente ausente en Ash1; los genes mostrados en la Tabla 5 representan casos en los que Ash1 tiene menos miembros de una familia multigénica. Tres genes adicionales (2 que codifican proteínas, 1 lncRNA) se mapearon a dos contigs no colocados, lo que proporcionará una guía para colocar esos contigs en futuras versiones del ensamblaje de Ash1.

Después de mapear los genes en Ash1, extrajimos las secuencias codificantes de las transcripciones que se mapearon completamente (cobertura igual al 100%), las alineamos con las secuencias codificantes de GRCh38, y llamamos a las variantes relativas a GRCh38 (ver la sección «Métodos»). Dentro de los 35.513.365 pb de estos transcritos codificantes de proteínas, encontramos 20.864 variantes de un solo nucleótido e indels. Catorce mil cuatrocientos noventa y nueve de estas variantes se encontraban dentro de las regiones «seleccionables» de GIAB para variantes de alta confianza, aunque 3963 de ellas se encontraban en regiones repetitivas «difíciles» de GIAB, para las que los alineamientos suelen ser ambiguos. De las 10.536 variantes que no se encontraban en estas regiones difíciles, 10.456 (99,2%) coincidían con el conjunto de variantes de alta confianza del GIAB. En las regiones difíciles, 3804/3963 (96,0%) coincidieron con el conjunto GIAB.

A continuación, anotamos los cambios en los aminoácidos causados por las variantes y el mapeo incompleto para todas las secuencias que codifican proteínas. De los 124.238 transcritos que codifican proteínas de 20.197 genes, 92.600 (74,5%) tienen secuencias proteicas 100% idénticas. Otros 26.566 (21,4%) tienen al menos un cambio de aminoácido pero dan lugar a proteínas con la misma longitud, y 1.561 (1,3%) tienen mutaciones que preservan el marco y que insertan o eliminan uno o más aminoácidos, dejando el resto de la proteína sin cambios. La tabla 6 muestra las estadísticas de todos los cambios en las secuencias de las proteínas. Si una proteína tenía más de una variante, la contamos bajo la variante más consecuente, es decir, si una proteína tenía una variante de sentido erróneo y un codón de parada prematuro, la incluimos en el grupo de «parada ganada».

De particular interés son aquellos transcritos con variantes que alteran significativamente la secuencia de la proteína y pueden resultar en una pérdida de función. Entre ellos se incluyen los transcritos afectados por un desplazamiento de marco (2158), pérdida de parada (58), ganancia de parada (416), pérdida de inicio (58) o truncamiento debido a un mapeo incompleto (564). Estas isoformas alteradas representan 885 loci de genes; sin embargo, 505 de estos genes tienen al menos otra isoforma que no está afectada por una variante alteradora. Esto deja 380 genes en los que todas las isoformas tienen al menos una interrupción; la lista completa se proporciona en el archivo adicional 1: Tabla S1.