Computational assessment of methods for inferring intratumoral genomic heterogeneity

Resumen

1. Evolución del cáncer y genómica El cáncer es el resultado de la acumulación de mutaciones genéticas y/o epigenéticas en células somáticas sanas que conducen al desarrollo y crecimiento de tumores. Si bien las fuerzas evolutivas como la mutación, la selección, la deriva genética y la migración dan forma al crecimiento y desarrollo del tumor, la evolución del cáncer somático difiere de la evolución de la línea germinal en varios aspectos. La diferencia más importante es que la evolución somática es asexual, se produce a través de la mitosis celular y, por lo tanto, sin recombinación genética. Como explico en el Capítulo 1, las células cancerosas evolucionan a lo largo de linajes clonales y la población celular que porta el mismo conjunto de mutaciones –también llamadas variantes– se llama clon. Sin embargo, debemos tener en cuenta que después de la mitosis las copias del genoma resultantes no son idénticas debido a errores durante la replicación del ADN, y que carecemos de la capacidad para detectar todas las mutaciones en un tumor. Por tanto, de ahora en adelante utilizaré una definición operativa en la que un clon es una población de células que comparte el mismo genotipo según las variantes somáticas detectadas en las muestras tumorales. Históricamente, los tumores individuales se han clasificado en tipos histológicos de cáncer, según el órgano en el que se originan. Sin embargo, dentro de un tipo de cáncer existe una gran variabilidad en la composición genómica de los tumores entre distintos pacientes. Aún más, un tumor dado contiene típicamente una multitud de poblaciones de células tumorales con distinta composición genómica, es decir, múltiples clones. Esta diversidad genética a menudo se denomina heterogeneidad genómica intratumoral (ITH). En los últimos años, los avances en las tecnologías de secuenciación de ADN de alto rendimiento (HTS) han permitido caracterizar la ITH con un detalle sin precedentes y revelar su profundo impacto en la progresión del cáncer, su tratamiento y la resistencia a los medicamentos. Las variantes genómicas son una característica dominante y definitoria de las células cancerosas. Dentro de un solo tumor existen diferentes tipos de variantes genómicas somáticas que son el resultado de una variedad de procesos mutagénicos, como las variantes de un solo nucleótido (SNV), inserciones o deleciones cortas (indeles) y variantes estructurales (SV), que incluyen las variantes numéricas de copia (CNV). Además, los polimorfismos de un solo nucleótido (SNP) de la línea germinal se encuentran en todas las células somáticas, ya estaban presentes en el cigoto del que se originan todas las células somáticas. Por el contrario, las mutaciones somáticas se encuentran solo en un subconjunto de células somáticas. Un desafío particular en la detección de variantes tumorales surge del hecho de que muchas variantes somáticas están presentes en una fracción relativamente pequeña de células tumorales. La proporción de células que porta una variante determinada se denomina prevalencia celular. Es importante destacar que el descubrimiento de que la ITH es prevalente en la mayoría de los tipos de cáncer implica que una sola muestra localizada, o biopsia, no es representativa de todo el tumor. De hecho, se ha demostrado que la secuenciación de una sola biopsia tiene un sesgo de muestreo significativo que puede mitigarse secuenciando por separado varias muestras. Como resultado, la secuenciación masiva multirregional (M-seq) de tumores se ha convertido en un procedimiento común en la investigación del cáncer. De este modo, se han desarrollado estrategias de llamada de variantes de múltiples muestras ad hoc para aprovechar explícitamente la heterogeneidad genómica contenida en múltiples muestras, en el tiempo y/o en el espacio. Los estudios de cáncer típicos secuencian millones de células a la vez (es decir, secuenciación conjunta de ADN en masa; ver más abajo), lo que implica la ausencia de secuencias individuales. La deconvolución clonal es el proceso de inferir clones y sus proporciones a partir de una o más muestras a granel (“bulk”) de un tumor. Para hacer esto, las mutaciones tumorales se agrupan en clones de acuerdo con sus frecuencias de alelos variantes, y se determina la prevalencia celular de cada clon en cada muestra. El problema de la deconvolución clonal se complica por el ruido inherente al proceso de secuenciación del ADN, así como por la complejidad del genoma tumoral en sí, como resultado de su inestabilidad genómica. Por lo tanto, se han diseñado una gran cantidad de herramientas para abordar la descomposición de muestras de tumores mixtos en los subclones constituyentes. 2. Motivación y objetivos de la tesis Tanto la detección de variantes como las herramientas de deconvolución clonal constituyen instrumentos cruciales en la caracterización de ITH. Sin embargo, no está claro qué tan bien se desempeñan las diversas herramientas en su tarea en diferentes condiciones. En esta tesis, mi objetivo principal es evaluar el rendimiento de los métodos para caracterizar la ITH en datos de tumores M-seq. Para lograr este objetivo, me propuse asumir los siguientes tres objetivos principales. Objetivo 1. Desarrollar una herramienta de simulación para la generación de datos M-seq que incorpore los atributos particulares de los genomas del cáncer. Estos atributos comprenden la presencia de múltiples clones relacionados evolutivamente, la inestabilidad genómica reflejada por las variantes del número de copias, la heterogeneidad espacial y la mezcla de células sanas, entre otros. Objetivo 2. Evaluar el rendimiento de los métodos de identificación de variantes somáticas basados ​​en datos M-seq de tumores simulados, para comprender con cuánta precisión se pueden capturar la carga mutacional y la estructura de las subpoblaciones tumorales con las estrategias actuales. Objetivo 3. Evaluar la precisión de los métodos de deconvolución clonal utilizando conjuntos de datos de tumores M-seq simulados, con respecto a la correcta identificación de clones, de su prevalencia celular y de sus relaciones ancestrales. 3. Simulación de heterogeneidad genómica intratumoral Para evaluar el desempeño de cualquier método de inferencia sobre cualquier tipo de datos, es necesario conocer las verdaderas características de los elementos que se infieren. En el caso de la ITH, los elementos principales son las variantes genómicas, los genotipos de los clones con respecto a estas variantes y la prevalencia celular de cada clon en cada muestra. Además de los elementos que se van a inferir, hay una multitud de factores que potencialmente influyen en el rendimiento de la inferencia, como el número de regiones tumorales muestreadas y la profundidad de secuenciación. En ausencia de datos con valores conocidos para los parámetros que se quieren inferir, las simulaciones mediante ordenador pueden ser una herramienta tremendamente útil para distintos propósitos. En el Capítulo 2 describo un marco computacional, llamado TumGenomSim, que desarrollé durante mi tesis para simular datos tumorales M-seq con el propósito de evaluar las estrategias de detección y descripción de la ITH. TumGenomSim cumple una serie de requisitos funcionales importantes cuya combinación no se encuentra en ninguna otra herramienta de simulación de tumores. Está guiada por la filogenia, ya que las mutaciones del tumor simulado son consistentes con una filogenia de clones subyacente que fija las relaciones ancestrales entre los clones. También considera el muestreo multirregional, generando múltiples muestras contemporáneas del mismo tumor. Además, es capaz de simular lecturas de secuenciación, haciendo coincidir los perfiles de error de diferentes plataformas HTS, pero siendo capaz al mismo tiempo de generar directamente recuentos de lecturas de HTS, lo que es mucho más eficiente desde el punto de vista computacional. Además, esta herramienta es capaz de generar variantes del número de copias, las cuales ocurren con frecuencia en muchos tumores, y que pueden distorsionar las frecuencias alélicas variantes de los SNV que se encuentran en estas regiones. Finalmente, TumGenomSim incluye un modelo explícito de evolución espacial para generar la distribución de clones entre las diferentes regiones tumorales. Además de TumGenomSim, que implementé en C++, también desarrollé un conjunto de funciones de utilidad en Python que he llamado simtools . La funcionalidad de las simtools incluye la generación de árboles de clones, la simulación de la evolución del rango geográfico y el muestreo repetido de los parámetros de simulación para generar conjuntos de réplicas para experimentos de simulación a gran escala. Una simulación completa utilizando simtools y TumGenomSim sigue una secuencia de pasos: (1) Se simula un árbol filogenético que describe las relaciones ancestrales entre los clones en un tumor. (2) La distribución espacial de los clones se define mediante una matriz de prevalencia que describe la frecuencia celular en el que cada clon está presente en cada muestra, generado siguiendo un modelo de dispersión, extinción y cladogénesis (DEC). (3) El usuario puede proporcionar un genoma de referencia o puede generarse in silico. (4) Los SNPs de la línea germinal se añaden al genoma de referencia para generar un genoma diploide sano. (5) Las variantes tumorales se generan como un conjunto de SNV somáticos, según las firmas mutacionales de COSMIC. (6) Se generan eventos de CNV que copian o eliminan una determinada porción del genoma. (7) Dependiendo del modo de ejecución, las lecturas de secuenciación o los recuentos de lectura se simulan en función de los genomas de cada clon. 4. Evaluación de estrategias de detección de variantes para datos secuenciación multirregional de tumores En el Capítulo 3 comparo el rendimiento de 16 herramientas de detección SNV en datos M-seq simulados y empíricos, en diferentes escenarios, incluida una cobertura de secuenciación variable y niveles crecientes de estructura espacial. De esta manera, evalué el rendimiento de estas herramientas para contestar distintas cuestiones específicas, como (1) la cuantificación de la carga total de mutaciones a nivel del tumor, (2) la detección de variantes en biopsias individuales y (3) la detección de variantes raras. Para generar los datos M-seq, utilicé mi herramienta TumGenomSim (descrita en el Capítulo 2) para evolucionar in silico genomas con SNPs y SNVs a lo largo de filogenias de clones determinadas. Generé un total de 120 réplicas de tumores con muestras multirregionales. Para cada tumor, construí un árbol filogenético aleatorio con seis clones. Fijé el número de regiones muestreadas en seis –una sana más cinco tumorales– y distribuí los clones entre ellas de acuerdo con tres escenarios de mezcla arbitrarios: baja (pocos clones por muestra, poca superposición de clones entre muestras); moderado (número moderado de clones por muestra; superposición moderada entre muestras); y alto (gran número de clones por muestra; alto solapamiento entre muestras). Basé cada genoma tumoral en una referencia de 3 Mb y 3000 SNPs de la línea germinal, e introduje 100 SNV somáticos en cada tumor, distribuidos entre los clones. Basándome en los genomas de cada clon, simulé lecturas de secuenciación de extremo emparejado (PE) de Illumina a diferentes niveles de profundidad de 30x, 50x, 100x y 300x. En este estudio consideré 16 estrategias de identificación de variantes que agrupé en tres categorías: “marginal”, “conjunta” y “de dos pasos”. En la estrategia “marginal”, incluí herramientas de detección de de variantes somáticas que se utilizan habitualmente para identificar SNVs en muestras emparejadas tumor/normales. En este caso utilicé estas herramientas en un escenario multirregional simplemente aplicándolas de forma independiente para cada una de las regiones simuladas. Bajo la estrategia “de dos pasos”, incluí MuClone y SNV-PPILP. MuClone incorpora información previa sobre la composición clonal del tumor para clasificar el conjunto de variantes candidatas; SNV-PPILP asume que todas las muestras son monofiléticas y aplica la suposición de una “filogenia perfecta” para filtrar un el conjunto de variantes. Bajo la estrategia “conjunta”, consideré MultiSNV, que emplea un modelo estadístico que considera simultáneamente todas las muestras del mismo tumor, Mutect2-multi, que agrega las lecturas de todas las muestras de tumores, y HaplotypeCaller, que aunque se desarrolló en principio para identificar variantes de la línea germinal entre individuos, también parece aplicable a poblaciones de células somáticas dentro de individuos. Para saber si las tendencias identificadas en las simulaciones también se observan al analizar datos reales, también identifiqué variantes usando las mismas estrategias en cuatro conjuntos de datos empíricos M-seq de secuenciación de exoma completo, con tres a cinco muestras regionales cada uno. Finalmente, evalué los requisitos computacionales de cada uno de los métodos. Mis resultados indican que para conjuntos de datos con baja profundidad de secuenciación o alta mezcla entre clones, Mutect2 en modo multimuestra es la herramienta más precisa. Sin embargo, en condiciones menos complejas, los métodos que utilizan simultáneamente la información de todas las muestras no funcionaron mejor que algunos de los métodos estándar que analizan muestra por muestra. Para describir la carga mutacional total o detectar variantes regionales, Mutect2 y MuTect1 parecen ser las mejores herramientas. Mutect2 y Neusomatic fueron superiores para caracterizar diferencias entre regiones, mientras que Mutect2 y Mutect1 proporcionan las distribuciones de frecuencia de alelos variantes más precisas. Como tendencia general, una mayor mezcla clonal parece implicar una menor sensibilidad (proporción de las variantes verdaderas que se capturan) a nivel de muestras individuales, aunque Mutect2-multi-F y MuClone se vieron menos afectados en este sentido. Para la precisión (proporción de las variantes capturadas que son verdaderas) observé la tendencia opuesta, de manera que la mayoría de los métodos lograron mejores resultados en el escenario de alta mezcla clonal. Con respecto a la profundidad de secuenciación, en general, la mayoría de las herramientas lograron una sensibilidad significativamente mejor a nivel de muestra individual a profundidades más altas. La precisión fue generalmente muy alta en todos los niveles de profundidad. Como era de esperar, se puede observar una clara compensación de rendimiento entre la sensibilidad y la precisión estadística al comparar el número de falsos positivos (FP) y falsos negativos (FN) para las diferentes frecuencias de alelos variantes (VAF) observadas. Las estrategias que tienden a evitar los FP a bajas frecuencias tienden a producir un mayor número de FN. Otros métodos (por ejemplo, MultiSNV, Mutect1, Mutect2-single, Strelka2) logran impulsar el umbral de detección hacia frecuencias más bajas. El número creciente de llamadas FP por Mutect2-multi-F hacia los VAF más bajos observados indica que en un cierto punto, una mayor recuperación de variantes verdaderas tiene el costo de más FP. En general, es importante tener en cuenta que las estrategias de identificación de variantes evaluadas se comportan de manera muy diferente entre sí, y conocer sus errores específicos puede servir para mejorarlos. Por ejemplo, VarDict genera muchos FP cuando las variantes están a baja frecuencia, lo que podrían evitarse aplicando un simple filtro. Por otro lado, SNooPer infiere erróneamente que muchos SNPs de la línea germinal son somáticos, lo que puede mitigarse mediante la comparación con SNP conocidos disponibles en bases de datos. 4. Comparativa de herramientas de deconvolución clonal para datos de secuenciación multirregional de tumores Con el fin de proporcionar una visión más completa de nuestra capacidad actual para describir con precisión la ITH, en el Capítulo 4 realicé simulaciones extensivas de conjuntos de datos M-seq de tumores con TumGenomSim con la intención de comparar el rendimiento de siete herramientas de deconvolución clonal (CD): Cloe, CloneFinder, LICHeE, PhyloWGS, PyClone, PyClone-VI y SciClone. Seleccioné estas herramientas en particular porque se utilizan habitualmente (con muchas citas, por ejemplo, LICHeE, PhyloWGS, PyClone, SciClone), o implemtaban métodos nuevos (por ejemplo, CloneFinder, PyClone-VI), tenían características interesantes (por ejemplo, relajar el modelo de mutación de sitios infinitos, como con Cloe), y, en todos los casos, porque era relativamente fáciles de ejecutar con cientos de conjuntos de datos en un tiempo razonable. En total, simulé 500 réplicas tumorales en varios escenarios biológicos que incluían distintos números de clones, números de muestras, pureza del tumor, profundidad de secuenciación, número de SNV y carga de CNV. Para cada réplica, simulé un genoma de 10 Mb, con 10 cromosomas, con 100,000 SNPs germinales. Incluí entre dos y diez clones, y analicé tres y diez muestras regionales por tumor. A cada tumor le asigné un grado variable de contaminación con células sanas, de modo que la pureza del tumor fuese de entre el 20 y el 100%. Además, simulé entre 10 y 10.000 mutaciones somáticas, una fracción variable de las cuales eran CNVs, siendo el resto SNVs. Para imitar un análisis típico, calculé el estado del número de copia local (CN) de los SNV a partir de recuentos de lecturas de HTS utilizando una herramienta de inferencia de CN, y usé solo SNV en regiones diploides. Para comparar el rendimiento de las herramientas de CD utilicé métricas que cuantifican la precisión con respecto al número inferido de clones, los genotipos de clones, la prevalencia de clones y la estructura del árbol de clones. Además, registré los requisitos computacionales de las diferentes herramientas. Idealmente, el número inferido de clones debería coincidir con el número simulado, pero en general estas dos variables se correlacionaron débilmente en el mejor de los casos. LICHeE mostró la correlación más alta, seguida por SciClone y PyClone-VI. Para algunas herramientas de CD, el error del número de clones inferidos se correlacionó positivamente con el número de muestras regionales, como CloneFinder, LICHeE y SciClone. Con una mayor pureza del tumor, los métodos de CD parecen inferir un número de clones más cercano al número real, pero esta tendencia no es muy fuerte. La profundidad de secuenciación no afectó significativamente el error en el número de clones inferidos, a excepción de CloneFinder, que se benefició de una mayor profundidad. Para todas las herramientas, aparte de LICHeE y PyClone, el error del número de clones se correlacionó positivamente con la proporción del genoma afectado por CNVs, es decir, por la carga de CNVs. Las diferencias en la precisión del genotipado entre Pyclone-VI, SciClone, LICHeE, CloneFinder y PyClone no fueron estadísticamente significativas, pero Cloe y PhyloWGS funcionaron peor que el resto. En general, la precisión del genotipo fue mejor con más de cuatro muestras, pero observé una clara correlación con el número de muestras regionales, con la pureza del tumor o con la profundidad de secuenciación. La mayoría de los métodos fueron en gran medida independientes del número de SNVs, aparte de PyClone y Cloe, que fueron menos precisos con conjuntos de SNVs más grandes. LICHeE, PyClone, PyClone-VI y SciClone se comportaron significativamente peor con una proporción cada vez mayor del genoma afectado por CNVs. Por el contrario, PhyloWGS mejoró con una mayor carga de CNVs. Aparte de PhyloWGS, todas las herramientas infirieron la prevalencia celular de los clones con bastante precisión, aunque el error fue menor cuando el número de clones simulados fue mayor en el caso de CloneFinder, LICHeE y SciClone. El número de muestras, el número de SNVs o la profundidad de secuenciación no influyó demasiado en la inferencia de la prevalencia clonal. La pureza del tumor se correlacionó positivamente con el error de prevalencia para todas las herramientas, excepto PhyloWGS. Cloe, LICHeE y SciClone arrojaron un error de prevalencia más alto para los conjuntos de datos con más carga de CNV, mientras que se redujo para PhyloWGS. LICHeE obtuvo los árboles de clones más precisos. Como tendencia general, los árboles inferidos fueron peores cuando se simuló un mayor número de clones, como era de esperar. El número de muestras regionales por tumor, la pureza del tumor, el número de SNVs o la profundidad de secuenciación no tuvieron una gran influencia en la precisión del árbol de clones. CloneFinder y PhyloWGS se beneficiaron de una mayor carga de CNV, mientras que LICHeE mostró efectos mixtos (beneficiosos para DISC, perjudiciales para CASet). Cloe no se vio afectado en gran medida por la carga de CNV. En total, dado mi modelo de simulación y los escenarios estudiados, un diseño experimental M-seq eficiente podría consistir en cuatro a nueve muestras regionales, una profundidad de secuenciación de 100 a 300x y 100 a 500 SNV. En estas condiciones, la pureza del tumor no debería ser una preocupación importante. Para analizar este tipo de experimentos, LICHeE parece ser la herramienta más competitiva de entre las estudiadas en este capítulo. 5. Conclusiones principales El trabajo presentado en esta tesis refleja mis esfuerzos para evaluar la precisión de los métodos actuales para la detección de variantes somáticas y la deconvolución clonal a partir de datos de secuenciación multiregional de tumores. Mis principales conclusiones son: La simulación de datos de secuenciación multiregional que reflejan parámetros tumorales importantes requiere una cierta complejidad del modelo generador, debiendo incorporar las relaciones ancestrales entre clones, un almacenamiento no trivial del genoma del clon que permita rastrear fácilmente las copias alélicas individuales de las distintas regiones genómicas y los correspondientes SNV, y un componente espacial que asegure la coherencia entre la distribución geográfica de los clones y la genealogía. En mi trabajo, hice un esfuerzo por incluir todos estos aspectos en un marco de simulación integral, denominado TumGenomSim, que formó la base para mis estudios de evaluación comparativa sobre la detección de variantes multirregionales, utilizando lecturas de secuenciación simuladas, y sobre deconvolución clonal, utilizando datos de recuento de lecturas simuladas. Espero que también pueda ser útil para otros investigadores. Las herramientas de detección de variantes somáticas que estudié utilizando datos simulados representaron diferentes opciones para equilibrar la sensibilidad y la precisión estadística. Algunas estrategias parecen poner más énfasis en la sensibilidad al precio de una menor precisión, y viceversa. La elección de la mejor herramienta dependerá finalmente de la pregunta de investigación. Según mi análisis, Mutect2 logró el mejor equilibrio, pero ofrece dos modos diferentes; por lo tanto, si la sensibilidad es más importante, Mutect2 en modo de muestra múltiple podría ser la mejor opción, y cuando la precisión sea más importante, Mutect2 en modo de muestra única puede ser más apropiado. Las variantes identificadas por las herramientas de detección en los datos empíricos de M-seq no se solapan demasiado entre distintas herramientas. En algunos casos, puede ser aconsejable confiar en más de una herramienta para identificar las distintas variantes somáticas con mayor precisión. Los enfoques de “dos pasos”, que tienen como objetivo aprovechar la información multirregional,, como la composición clonal o las restricciones filogenéticas entre muestras, dependen en gran medida de la precisión de la información externa. En mis simulaciones estos métodos son menos precisos que los que identifican variantes basadas exclusivamente en los datos de secuenciación. La mayoría de las herramientas de deconvolución clonal estudiadas funcionaron comparativamente bien con un número intermedio de muestras, una profundidad de secuenciación moderada y un número relativamente bajo de mutaciones. Sin embargo, la precisión de las herramientas de deconvolución clonal finalmente puede depender más de las propiedades inherentes al tumor que se analiza, que del diseño experimental. Finalmente, el problema de la deconvolución clonal a partir de datos de tumores M-seq todavía parece estar lejos de resolverse con un alto grado de precisión, al menos en las condiciones que yo he estudiado en esta tesis. Algunas de las herramientas de deconvolución clonal desarrolladas hace algunos años pueden no ser adecuadas para los conjuntos de datos HTS actuales con muchos miles de SNV y / o muchas muestras regionales debido a limitaciones de tiempo de ejecución. Se requieren algoritmos más eficientes y las mejoras de métodos más antiguos, como PyClone-VI, son una contribución muy útil al campo.