Capítulo 4 Análisis masivo e integrador de conjuntos de genes

4.1 Motivación

El surgimiento y rápido avance de las tecnologías de secuenciación de alto rendimiento han llevado a la disponibilidad de miles de experimentos en repositorios públicos como son el Gene Expression Omnibus (Barrett et al., 2012) y Array Express (Kolesnikov et al., 2014). Más aún, se han llevado a cabo varios proyectos de gran escala, centrados en la comprensión integrada de la complejidad de las enfermedades humanas, como The Cancer Genome Atlas (TCGA) (Weinstein et al., 2013) y el US-LACRN (Llera et al., 2015) para el estudio del cáncer. Comúnmente, los conjuntos de datos disponibles incluyen diferentes tipos de datos provenientes de múltiples tecnologías ómicas para los mismos sujetos, proporcionando información a diferentes niveles moleculares y creando oportunidades sin precedentes para estudiar enfermedades humanas. La integración de estos datos genéticos, genómicos y proteómicos puede conducir a una mejor identificación de grupos homogéneos de sujetos, así como de patrones biológicos comunes y distintivos entre los grupos. La caracterización de grupos de sujetos que presentan una misma condición patológica, contribuye al desarrollo de estrategias diagnósticas y de tratamiento más adecuadas (Chen & Snyder, 2013; Kedaigle & Fraenkel, 2018). En este sentido, el auge de la medicina personalizada y la disponibilidad de repositorios de datos moleculares de alto rendimiento han aumentado la necesidad de herramientas adecuadas que permitan a los investigadores traslacionales gestionar y explorar estas fuentes de datos (Canuel, Rance, Avillach, Degoulet, & Burgun, 2014; Fernandez & Casares, 2018).

La integración y el análisis de grandes conjuntos de datos provenientes de múltiples fuentes ómicas se ha convertido en una tarea cada vez más desafiante. Los primeros esfuerzos para identificar patrones distintivos relacionados con fenotipos o contrastes (comparaciones entre pares de fenotipos) específicos se han basado en el análisis de expresión diferencial. Estos esfuerzos condujeron a la identificación de diferentes listas de genes, proteínas u otras características genómicas, con una ligera superposición entre los conjuntos de datos (Creighton et al., 2018; Hoadley et al., 2014; Korkola et al., 2007; Meng, Kuster, Culhane, & Gholami, 2014; Metzger Filho, Ignatiadis, & Sotiriou, 2011; Rodriguez et al., 2016a; Verhaak et al., 2010; Weinstein et al., 2013). Por el contrario, cuando la estrategia de análisis se basaba en el enriquecimiento funcional, se identificaron grandes cantidades de vías biológicas, funciones moleculares y procesos biológicos en común entre diferentes conjuntos de datos (Creighton et al., 2018; Meng et al., 2014; Rodriguez et al., 2016a). Estos resultados denotan que, a pesar de las posibles diferencias en las características ómicas probadas en varios ensayos, la biología subyacente presenta características similares, lo que lleva a resultados fenotípicos o pronósticos similares. Por lo tanto, el análisis funcional de Conjuntos de Genes (CG) resulta en un mejor enfoque para la integración de múltiples conjuntos de datos.

El principal inconveniete de los análisis de expresión diferencial o de enriquecimiento funcional es la falta de integración de datos de múltiples ómicas. Esto ha motivado al surgimiento de un campo de investigación, conocido como Big Omics Analysis, en el que la integración de tipos de datos dispares se ha vuelto cada vez más importante para capturar la heterogeneidad de los procesos biológicos (Fernandez & Casares, 2018; Liu, Shen, & Pan, 2016). Los enfoques recientes se basan en aprendizaje automático (Shen, Olshen, & Ladanyi, 2009), análisis de redes (Wang et al., 2014), o fusión de patrones (Shi et al., 2017). Estas herramientas se basan principalmente en el análisis de expresión o patrones de coexpresión de un único conjunto de datos, lo que supone una gran limitación, ya que se debe contar con datos de cada ómica para cada sujeto. Además, ninguna de estas herramientas es capaz de analizar el comportamiento funcional de la enfermedad o fenotipo; y mucho menos analizar multiples cohortes (poblaciones). Si bien una herramienta ampliamente utilizada llamada PARADIGM (Vaske et al., 2010) ha sido desarrollada para analizar los aspectos funcionales de múltiples fuentes ómicas, sólo puede proporcionar información a nivel de sujetos individuales, y por lo tanto, omite el análisis de grupos de interés, fenotipos, conjuntos de datos o tecnologías ómicas.

El análisis de grandes fuentes de datos ómicas resulta particularmente útil para el estudio de enfermedades heterogéneas, como el Cáncer de Mama (CM) (Korkola et al., 2007; Kwa, Makris, & Esteva, 2017; Reis-Filho & Pusztai, 2011). En este sentido, Perou et al. (Perou et al., 2000) definió cuatro subtipos intrínsecos de CM -introducidos en la Sección 2.3.1- que presentan diferentes perfiles. Estos subtipos han sido ampliamente evaluados a través del análisis independiente de diferentes fuentes de datos incluyendo genes, expresión de microARN, variación del número de copias, y proteínas (Network & others, 2012). Aunque estos análisis condujeron a la validación del esquema de clasificación de Perou et al., también han dado lugar a más preguntas sobre los propios subtipos. Por ejemplo, se ha descubierto que los sujetos con CM están agrupados en aún más grupos de subtipos (Curtis et al., 2012). Más aún, se han identificado nuevos grupos de CM dentro de los subtipos previamente definidos (Aure et al., 2017). Asimismo, algunos sujetos con CM no pudieron ser asignados a ningún subtipo existente (Fresno et al., 2016). Por lo tanto, el éxito de la medicina traslacional en CM aún se ve limitado por la capacidad de las herramientas de análisis integradoras y masivas para descifrar los procesos biológicos aún desconocidos que subyacen a cada uno de estos subtipos de CM. Para lograr este éxito, resulta fundamental una herramienta que permita al investigador comparar bases de datos tanto de diversas fuentes ómicas como de diversas poblaciones o incluso enfermedades.

4.2 Adaptación del IFA a otras fuentes de datos

Como se mencionó en la Sección 3.5.6, el pipeline del IFA lleva a cabo sus análisis utilizando los métodos dE_F_none y mGSZ. Para llevar a cabo dichos análisis, la lista de genes Diferencialmente Expresados (DE) y el rankeo son calculados mediante un modelo lineal obtenido por la función eBayes de la librería limma de R. La obtención de los genes DE se lleva a cabo como un paso anterior a utilizar dE_F_none, sin embargo, para el caso de mGSZ, el rankeo de los genes se realiza internamente por dicha librería.

El ajuste del modelo lineal de eBayes asume que los datos de entrada siguen una distribución Normal, lo cual, como se vió en la Sección 2.1, se cumple para datos provenientes de microarreglos de ADN y de iTRAQ. Sin embargo, este supuesto no se cumple para datos de Secuenciación de ARN, donde vimos que, al ser datos de conteos, se asemejan a una distribución de Poisson o Binomial Negativa.

Dado el interés de desarrollar una herramienta que permita comparar, desde un punto de vista funcional, una cantidad masiva de bases de datos, tanto de diversas poblaciones como de una variedad de fuentes ómicas. Es por esto que resultó necesario adaptar el pipeline previo del IFA de modo que permita llevar a cabo el análisis a partir de datos de conteos. En particular, se debió inspeccionar el paquete R mGSZ y realizarle las modificaciones pertinentes sobre la metodología de rankeo de genes.

El IFA utiliza una misma herramienta (y modelo estadístico), tanto para el ASR como para el análisis de PFC, para lograr un análisis unificado. La adaptación del IFA a datos de conteos se realizó manteniendo esta idea de unificación. En este sentido, tanto para la obtención de los genes DE como para rankearlos, se implementó una modificación del pipeline IFA que dependiendo del tipo de datos de entrada utiliza la función eBayes (para datos Normales) ó voom (para conteos). La función voom también se encuentra desarrollada en la librería limma, y se desarrolló principalmente con el fin de detectar genes DE para datos de conteos, la idea de esta propuesta es poder modelar los datos como si fueran Normales mediante mínimos cuadrados generalizados (Law, Chen, Shi, & Smyth, 2014). Para ello se realiza una transformación sobre los datos de conteos de manera de poder corregir la relación media-varianza de los datos, a partir de un ajuste loess. De esta manera, a partir de la matriz de conteos se obtiene una matriz de pesos asociada a la corrección necesaria y por consiguiente poder ajustar el modelo lineal mediante mínimos cuadrados generalizados. Como resultado de esta reimplementación, se obtuvo una alternativa del IFA que, de un modo transparente al usuario, permite llevar a cabo el análisis a partir tanto de datos Normales como de conteos.

4.3 Herramienta desarrollada

En el presente capítulo se introduce la librería R desarrollada: MIGSA -del inglés Massive and Integrative Gene Set Analysis-. MIGSA tiene como objetivo identificar la presencia de patrones funcionales comunes o distintivos de diversos fenotipos. Para el análisis, los términos biológicos, clínicos, o definidos por el usuario (de interés específico) se definen mediante CG. MIGSA permite a los usuarios realizar un análisis comparativo integrador de los datasets combinando diferentes experimentos de diferentes plataformas ómicas. Nuestra librería supera la pérdida de información producida por la fusión de diferentes fuentes de datos, explota la complementariedad de los diferentes niveles ómicos y permite realizar consultas directas y supervisadas basadas en cada fenotipo.

Como se explica esquemáticamente en la Figura 4.1, MIGSA toma como entrada una lista de matrices de expresión (\(L=\{M_d\}; d \in \{1,\ldots,D\}\)). Cada matriz contiene datos del \(d\)-ésimo dataset, con \(S_d\) sujetos medidos a través de cualquier tecnología ómica (por ejemplo, genes, proteínas, etc.) que pueda ser anotada con identificadores Entrez. Para el \(d\)-ésimo conjunto de datos, MIGSA requiere la identificación del fenotipo de cada sujeto \(F_d=[f_{d1},...,f_{dS_d}]\), con \(f_{di} \in \{a,b,...\}\) \((i\in \{1,...,S_d\})\), para identificar cada experimento \(E_{d}^{a,b}\) que resulta de la matriz de expresión \(M_d\) y el contraste de fenotipo a vs. b. Luego, en paralelo, MIGSA aplica individualmente a cada experimento el Análisis Funcional Integrador (IFA) descrito en el Cápitulo 3. Vale la pena mencionar que el IFA aplicado por MIGSA presenta dos novedades respecto a su publicación original (Rodriguez et al., 2016a): ahora permite analizar datos de conteos (como los de secuenciación de ARN), microarreglos de ADN, y proteómicos de iTRAQ; e identifica los genes específicos que contribuyeron al enriquecimiento de cada CG.

Figure 4.1: Flujo de trabajo del Análisis Masivo e Integrador de Conjuntos de Genes (MIGSA). De una lista de matrices de expresión \(L=\{M_1,....,M_D\}\), cada matriz con la identificación del fenotipo (por ejemplo \(a\), \(b\), \(c\), \(d\)) para cada muestra, se realiza la lista de experimentos (Es), donde cada experimento compara un par de fenotipos de una matriz (por ejemplo \(E_{1}^{a,b}\) contrasta los fenotipos \(a\) y \(b\) para la matriz de expresión \(M_1\)). MIGSA toma como entrada esta lista de Es y aplica, individualmente a cada experimento, una versión mejorada del IFA previamente presentado, de manera paralela. Esta versión mejorada del IFA toma una matriz de expresión y, de acuerdo con su tipo de datos, utiliza una función \(f\) diferente para obtener tanto la expresión diferencial como el rankeo de los genes, y luego realiza ambos Análisis de Sobre-Representación (ASR) y de Puntuación Funcional de Clase (FCS). Para este paso, se puede utilizar una colección generada por el usuario de conjuntos de genes ó seleccionar entre más de 130 colecciones proporcionadas por MIGSA. Los resultados de cada experimento se almacenan en un cubo de datos Q. Como salida, MIGSA proporciona diferentes alternativas de exploración y visualización, permitiendo identificar fácilmente los Conjuntos de Genes Enriquecidos (CGE) para cada experimento y, para cada conjunto de genes, los Genes Importantes para su Enriquecimiento (GIE).

Los resultados de MIGSA se resumen en un cubo de datos tridimensional, \(Q^{Es,CGE,GIE}\), con tres ejes o niveles de abstracción: Experimentos (Es; es decir, contrastes de fenotipos), Conjuntos de Genes Enriquecidos (CGE; CG estadísticamente desregulados para un umbral de significancia), y la lista de Genes Importantes para el Enriquecimiento (GIE; es decir, para cada CG, los genes que fueron los principales responsables de su desregulación en un experimento específico). Así, para el experimento \(E_{d}^{a,b}\), MIGSA proporciona la lista de conjuntos de genes enriquecidos \(CGE_{d}^{a,b}\), y la lista de genes importantes de enriquecimiento \(GIE_{d}^{a,b}\) (ver Figura 4.1). Nótese que en el cubo \(Q\) los ejes son: \(Es=\cup_d\cup_{a,b}E_{d}^{a,b}\), \(CGE=\cup_d\cup_{a,b}CGE_{d}^{a,b}\), y \(GIE=\cup_d\cup_{a,b}GIE_{d}^{a,b}\).

Un gen \(G\) será considerado un GIE para un determinado conjunto de genes \(CG\) y experimento \(E_{d}^{a,b}\), si \(G \in CG\), y al menos uno de:

• \(CG\) fue enriquecido por el análisis de sobre-representación, y el gen \(G\) se encontró diferencialmente expresado en \(E_{d}^{a,b}\).
• \(CG\) fue enriquecido por puntuación funcional de clase, y el gen \(G\) conformó el leading-edge (Subramanian et al., 2005) de \(CG\), es decir, la posición de \(G\), en el rankeo, se encuentra a la izquierda/derecha de la gráfica positiva/negativa del Enrichment Score de \(CG\).

MIGSA proporciona varias herramientas gráficas, que permiten a los usuarios explorar y visualizar fácilmente los resultados de sus datos. Por ejemplo, el primer eje del cubo MIGSA (\(Q^{CGE,Es}\)) puede representarse fácilmente mediante un heatmap. En este caso, los CG se distribuyen en filas y los experimentos en columnas. En particular, las filas y columnas se encuentran ubicadas de acuerdo a un ordenamiento dado por la distancia Jaccard y enlace promedio. El heatmap de MIGSA resume eficazmente los resultados de múltiples IFA y ofrece a los usuarios la capacidad de visualizar patrones de enriquecimiento funcional relacionados con contrastes, fenotipos o incluso con fuentes de datos específicas. El heatmap incluye un dendrograma superior que muestra los grupos formados por los resultados de todos los experimentos. Este dendrograma permite una rápida identificación de clusters y sub-clusters de experimentos, y su asociación con los fenotipos estudiados.

4.4 Validación de la herramienta

Para demostrar la utilidad de MIGSA como herramienta de descubrimiento, se utilizó con el fin de lograr una caracterización funcional de cada subtipo molecular PAM50 de CM. Se aplicó MIGSA para analizar una gran cantidad de bases de datos de CM donde se evaluaron todas las combinaciones posibles de pares de subtipos. La eficiencia de MIGSA se vería reflajada en encontrar genes y patrones funcionales que caractericen los diversos fenotipos analizados. En este sentido, dado que se analizaron experimentos pertenecientes a una misma enfermedad, será esperable encontrar patrones en común compartidos entre todos los contrastes. Pero también, y quizás de mayor interés biológico, resultará fundamental detectar patrones específicos que caractericen cada contraste, y más particularmente, cada subtipo.

4.5 Datos de entrada

4.5.1 Matrices de expresión

Se analizaron veinticinco conjuntos de datos de expresión génica de cáncer de mama basados en microarreglos de ADN (4.145 sujetos) previamente reportados por Haibe-Kains et al. (Haibe-Kains et al., 2012) e identificados como microarreglos comerciales por Fresno et al. (Fresno et al., 2016), estos conjuntos de datos son los de CM presentes en la Tabla ??, excepto por los de TCGA. Para cada sujeto de los datasets de Haibe-Kains (HKds), se tomó la clasificación PAM50 basada en el algoritmo pbcmc, según se detalla en Fresno et al. El algoritmo pbcmc devuelve un nivel de confianza de que un sujeto efectivamente pertenezca a un subtipo o no, de este modo se evita analizar sujetos que posiblemente estén mal asignados por genefu, y que por ende, provoquen un sesgo en los resultados. Se conservaron aquellos sujetos clasificados efectivamente como Basal, Luminal A, Luminal B y Her2+, mientras que los sujetos clasificados como Normal o No Asignados fueron descartados, lo que resultó en un total de 2.756 sujetos analizados.

Además de los HKds, se analizó una población de 97 sujetos que cuentan tanto con datos transcriptómicos (secuenciación de ARN y microarreglos de ADN) y proteómicos (iTRAQ) de TCGA; se seleccionaron los mismos sujetos para cada una de las bases de datos para evitar posibles sesgos. La clasificación PAM50 para estos sujetos se obtuvo a través de la librería pbcmc utilizando los datos de microarreglos y la opción de estandarización de escala. Mediante pbcmc se obtuvieron 25, 25, 12, 20, 0, y 15 sujetos identificados como Basal, Her2+, Luminal A, Luminal B, Normal, y No Asignados, respectivamente.

Para cada conjunto de datos, se promediaron los genes con múltiples sondas o detecciones. Para cada contraste de dos subtipos, sólo se incluyeron los genes detectados en al menos el 70% de los sujetos de cada subtipo. Además, para los datos del secuenciación de ARN, sólo se incluyeron los genes con más de 10 conteos en promedio por subtipo. Para un conjunto de datos dado y para cada subtipo dado, si había menos de ocho sujetos presentes, entonces estos sujetos fueron descartados para evitar tener una población no representativa, y en consecuencia no se analizaron los contrastes con ese subtipo para ese conjunto de datos.

Luego de los filtros mencionados, se obtuvieron 118 experimentos que contrastan los subtipos de CM para los HKds y 18 para TCGA. Para el análisis de sobre-representación del IFA de MIGSA, se eligieron parámetros de selección de expresión diferencial que incluyan entre el 4 y 6% de los genes como diferencialmente expresados.

4.5.2 Conjuntos de genes

MIGSA permite realizar el IFA analizando más de 130 bases de datos disponibles de CG conocidos (Kuleshov et al., 2016), como así también utilizar cualquier colección de CG especificada por el usuario. En este trabajo se obtuvieron (org.Hs.eg.db v3.5.0) y analizaron las tres categorías de CG de GO (Consortium, 2016), alcanzando 15,796, 1,902, y 4,604 CG de las ontologías de Procesos Biológicos (PB), Componentes Celulares (CC), y Funciones Moleculares (FM), respectivamente. Un CG se consideró CGE si obtuvo un p-valor \(\le 0,01\) por el IFA.

4.6 Estrategias para el análisis de eficiencia

Como primer instancia, se aplicó MIGSA sobre los experimentos de los HKds. Con el cubo resultante de MIGSA, se combinaron los resultados para definir los conceptos:

• Términos Enriquecidos en Consenso por Contraste (TECC): Para cada contraste, esos CG (y GIE) enriquecidos en al menos K% de los experimentos del contraste (donde K es un nivel definido por el usuario, aquí usamos \(K=50\)), es decir, \(\boldsymbol{TECC}^{a,b}=\{CG:\sum_{j=1}^D I(CG \in CGE_j^{a,b}) \ge \frac{K}{100}N^{a,b}\}\) donde I es la función indicadora, D es el número de datasets (aquí \(D=25\) HKds), y \(N^{a,b}\) el es el número de experimentos que contrastan el subtipo a con b.
• Perfil Transcriptómico Funcional (PTF): Para cada subtipo, aquellos TECC (y sus respectivos GIE) presentes en común al contrastar un subtipo contra cada uno de los restantes, es decir, \(\boldsymbol{PTF}_a=\cap_{b \neq a} \boldsymbol{TECC}^{a,b}\).

Luego, se analizaron los TECC para obtener un perfil funcional para cada contraste de los subtipos de CM. Adicionalmente, se obtuvieron los PTF basados en los experimentos de los HKds y se exploraron conjuntamente con los GIE correspondientes, para obtener una caracterización funcional de cada subtipo de CM.

Posteriormente, para demostrar la capacidad del paquete MIGSA de integrar múltiples conjuntos de datos ómicos, se aplicó MIGSA a los conjuntos de datos de TCGA de CM, y los resultados se compararon con los TECC definidos con los experimentos de los HKds. Para evaluar la concordancia entre los resultados de MIGSA sobre los conjuntos de datos de transcriptómica y proteómica del TCGA, se analizaron los CGE (y los GIE) para definir:

• Términos Enriquecidos por Transcriptómica (TET): Para cada subtipo, aquellos CG enriquecidos consistentemente en experimentos transcriptómicos. Para identificar los TET, primero, para cada contraste, se obtuvo la intersección de CGE tanto en los datos de secuenciación de ARN como en los de microarreglos de TCGA. Luego, para cada subtipo, el conjunto de términos presentes en los tres contrastes del subtipo dado se denominó Términos Enriquecidos por Transcriptómica.
• Términos Enriquecidos por Proteómica (TEP): Para cada subtipo, aquellos CG enriquecidos consistentemente en experimentos proteómicos. Para obtener los TEP, primero, se identificaron los CGE encontrados en los experimentos de iTRAQ del TCGA para cada contraste. Luego, se seleccionó el conjunto de términos presentes en todos los contrastes de un subtipo dado.

Finalmente, para cada subtipo, se evaluó la concordancia en la detección de CGE a través de diferentes plataformas y conjuntos de datos ómicos comparando los TET y TEP del TCGA, con los PTF encontrados analizando los experimentos de los HKds. Utilizando MIGSA, se llevo a cabo una inspección de los árboles de GO para encontrar intersecciones entre PTF, TET y TEP para las categorías de PB, FM y CC de cada subtipo.

4.7 Resultados

4.7.1 MIGSA sobre datasets de microarreglos

Se analizaron, con MIGSA, los 118 experimentos de HKds. En la Figura 4.2 se observa el heatmap que ilustra el primer eje del cubo de resultados (\(Q^{CGE,Es}\)), donde el dendrograma superior revela los clusters formados por los resultados de todos los experimentos que contrastan los subtipos de CM. Al analizar los dendrogramas se nota la conformación de clusters generados principalmente por los subtipos contrastados (no se encontró efecto de lote o plataforma). En particular, el cluster más grande contiene todos los experimentos que contrastan el subtipo Luminal A (marca LA en la Figura 4.2; 63 experimentos/columnas). Dentro de este cluster, se identificaron dos sub-cluster distintos. Uno de estos sub-clusters consiste en 23 de un total de 25 (92%) experimentos evaluando Luminal A contra Basal (marca LA-Ba). El otro sub-cluster incluyó 18 de 20 (90%) experimentos contrastando los dos subtipos Luminales (marca LA-LB). Los experimentos que comparan Luminal A contra Her2+ (LA-H) fueron agrupados principalmente en dos clusters que contienen seis y nueve de 18 (83%) experimentos respectivamente.

El segundo gran cluster divisable en la Figura 4.2 incluye todos los experimentos que contrastan los subtipos Basal con Luminal B ó Her2+ (marca Ba). En este cluster, se identificaron dos sub-clusters, uno con 15 de 20 (75%) experimentos contrastando Basal con Luminal B (marca Ba-LB), y otro con 12 de 18 (67%) experimentos contrastando Basal con Her2+ (marca Ba-H). Además de estos dos grandes sub-clusters, se diferencia uno con experimentos que contrastan Luminal B con Her2+ (marca LB-H), con 14 de 17 (82%) experimentos. En particular, este último cluster tiene la menor cantidad de términos enriquecidos, lo que sugiere una gran similitud entre estos dos subtipos.

Figure 4.2: Conjuntos de genes enriquecidos por análisis funcional para cada experimento. Heatmap: conjuntos de genes en filas, experimentos en columnas; ordenados por distancia Jaccard y enlace promedio. Conjuntos de genes enriquecidos en rojo. Las casillas S1 y S2 contienen conjuntos de genes enriquecidos en común para subtipos, y las casillas C1, C2, C3 y C4 para contrastes (comparación entre pares de subtipos). Barras de columnas (de arriba abajo): datasets (cada color representa uno de los 25 datasets analizados); contraste; subtipo1 y subtipo2, salmón es Luminal A (LA), turquesa es Basal (Ba), púrpura es Her2+ (H), y el verde es Luminal B (LB). Clusters de dendrogramas: Se demarcan para los subtipos Basal (Ba) y Luminal A (LA), y uno por cada contraste.

El heatmap de la Figura 4.2 también reveló patrones de enriquecimiento comunes específicos de subtipos para todos los subtipos contrastados. Por ejemplo, las cajas S1 muestran CGE que se encuentran comúnmente en los contrastes con el subtipo Basal, mientras que la caja S2 indica patrones de enriquecimiento de Luminal A. Además, también se observan patrones de enriquecimiento específicos de contrastes, es decir, conjuntos de términos que denotan diferencias entre dos subtipos determinados. Por ejemplo, la caja C1 contiene aquellos CG que diferencian Luminal A de Luminal B, la caja C2 para Luminal A vs. Her2+, C3 para Basal vs. Luminal B, y C4 para Luminal B vs. Her2+.

Las capacidades exploratorias de MIGSA se utilizaron para sub-dividir el heatmap mostrado en la Figura 4.2 en sub-heatmaps específicos de cada subtipo (es decir, Basal vs. todos los demás, Her2+ vs. todos los demás, Luminal B vs. todos los demás, y Luminal A vs. todos los demás). Estos sub-heatmap proporcionaron una visión detallada de los patrones de enriquecimiento de cada subtipo, y pueden observarse en (Rodriguez, Merino, Llera, & Fernández, 2019). En particular, los dendrogramas de estas figuras mostraron una mejor agrupación por contrastes, es decir, clusters para cada contraste que agruparon más experimentos que en la Figura 4.2.

El número de TECC y de términos en PTF encontrados se presenta en la Tabla ??. Como puede observarse, el número más bajo de TECC (184) se encuentra para el contraste Luminal B vs. Her2+, mientras que el contraste Luminal A vs. Luminal B muestra el valor más alto, con 983 CG enriquecidos. En términos de los PTF, Luminal B y Her2+ fueron los que presentaron el menor número de términos, mientras que Luminal A alcanzó el mayor número de CG en su PTF, seguido por el subtipo Basal.

Cabe señalar que los PTF no son mutuamente excluyentes, lo cuál era de esperar ya que todos los fenotipos contrastados son subtipos de la misma enfermedad. Veintisiete CG se encontraron presentes, en común, en los PTF de los cuatro subtipos de CM, mientras que 506 CG fueron únicos (exclusivos) para un solo subtipo. El número de CG exclusivos observados para cada PTF fue de 103 para Basal, uno para Luminal B, siete para Her2+ y 395 para Luminal A. Además, se encontraron 33 CG en común entre dos PTF diferentes: diez en \(\boldsymbol{PTF}_{Basal} \cap \boldsymbol{PTF}_{Her2+}\), ocho en \(\boldsymbol{PTF}_{Basal} \cap \boldsymbol{PTF}_{Luminal B}\), siete en \(\boldsymbol{PTF}_{Basal} \cap \boldsymbol{PTF}_{Luminal A}\), y ocho en \(\boldsymbol{PTF}_{Luminal B} \cap \boldsymbol{PTF}_{Luminal A}\).

Conjuntos de genes en consenso. Número de Términos Enriquecidos en Consenso por Contraste (TECC), definidos como aquellos conjuntos de genes enriquecidos en al menos el 50% de los experimentos evaluados para cada contraste (fila vs. columna). La columna “Intersección” indica el número de conjuntos de genes encontrados en todos los TECC de cada subtipo dado, lo cual define el Perfil Transcriptómico Funcional (PTF) de cada subtipo de cáncer de mama.

Contraste	Basal	Luminal B	Her2+	Luminal A	Intersección
Basal	-	488	437	768	155
Luminal B		-	184	983	44
Her2+			-	707	44
Luminal A				-	437

4.7.2 Exploración de los PTF de los subtipos de cáncer de mama

La exploración de los PTF reveló varios CG de GO previamente relacionados con los subtipos de CM. La visualización basada en árboles de la estructura GO, proporcionada por MIGSA, se utilizó para visualizar el PTF de cada subtipo en las categorías PB, CC y FM de GO. Los árboles de GO de los PTF de cada subtipo pueden encontrarse en (Rodriguez et al., 2019), donde cada término fue codificado por colores de acuerdo al PTF al que pertenece. La visualización en árbol nos permitió explorar más a fondo las relaciones entre los PTF de cada subtipo. La exploración de los árboles reveló que los términos no estaban dispersos al azar, sino que se encontraban conformando ramas a diversas profundidades, dominadas principalmente por subtipos particulares. Además, los CG presentes en común en los PTF de todos los subtipos fueron encontrados ubicados en nodos poco profundos que involucran términos generales y comúnmente encontrados, tales como citoplasma, desarrollo de sistemas, muerte y proliferación celular.

El análisis del árbol de PB reveló que las ramas que involucran términos relacionados con la regulación de la vía de señalización de los receptores de estrógeno sólo aparecieron en el PTF del subtipo Basal. Esto está motivado por el hecho de que, genes implicados en estos CG, incluyendo ESR1, se han encontrado previamente sobre-expresados en los subtipos Luminales y en algunos tumores Her2+ (Network & others, 2012), y sub-expresados en los sujetos Basal (Valentin, Da Silva, Privat, Alaoui-Jamali, & Bignon, 2012). Por esta razon, estos términos no se encontraron entre los TECC de los contrastes que involucran a los subtipos Luminales y HER2+ y, en consecuencia, no estaban presentes en el PTF de estos subtipos. Además, los términos asociados con el filamento de actina, la transición epitelio-mesénquima y el estrés del retículo endoplásmico aparecieron en el árbol de PB como PTF de Basal. Estas ramas son procesos bien conocidos relacionados con los tumores Basales, que los diferencian de otros subtipos (Guen et al., 2017). Los genes MLPH y FOXC1 se encontraron entre los GIE presentes en estas ramas, y fueron consistentemente sub y sobre-expresados, respectivamente, en todos los experimentos en los que se contrastó Basal. Cabe mencionar que estos dos genes están involucrados en una firma molecular recientemente desarrollada para cánceres de mama triple negativos, altamente relacionados con tumores Basales (Santuario-Facio et al., 2017). Además, los GIE presentes en términos relacionados con la diferenciación celular, como GATA3 y FOXA1, se encontraron desregulados en el subtipo Basal para todos los contrastes evaluados, en concordancia con los hallazgos del “Grupo del Atlas del Genoma del Cáncer” (Network & others, 2012). A su vez, ERBB4, un oncogén drogable (un gen cuya expresión puede afectarse con drogas), se encontró desregulado en el 95% de los contrastes con Basal evaluados, lo que sugiere que los sujetos Basales no se verían afectados por los fármacos que ataquen el gen ERBB4.

El subtipo Luminal A resultó ser muy diferente de los demás en lo que se refiere a los procesos de regulación del ciclo celular y replicación del ADN, en donde los términos como la transición G1/S del ciclo celular y la unión del origen de la replicación del ADN aparecieron como su PTF. Es bien sabido que los procesos de progresión a través del ciclo celular y la replicación del ADN están profundamente involucrados en la proliferación (Hanahan & Weinberg, 2000), que es la principal diferencia entre Luminal A y los subtipos más agresivos como Luminal B, Her2+, y Basal (Prat & Perou, 2011; Yersal & Barutca, 2014). En particular, estos términos muestran la principal diferencia entre los dos subtipos Luminales, ya que los tumores Luminal B expresan clásicamente niveles más altos de Ki67 (Cheang et al., 2009). Adicionalmente, se encontró que los GIE que provocan la proliferación celular, tales como AURKA, PCNA, CHEK2, y RAD51, estaban sub-expresados en Luminal A con respecto a los otros subtipos. Mientras que, genes como ERBB2, GRB7, y FGFR4 fueron encontrados sobre-expresados en todos los experimentos en los que Her2+ fue contrastado, en concordancia con lo informado por el “Grupo del Atlas del Genoma del Cáncer” (Network & others, 2012).

4.7.3 Integración de resultados de MIGSA de TCGA y Haibe-Kains

Los resultados de MIGSA para cada tecnología de expresión del TCGA se combinaron con los resultados de los HKds. La Figura 4.3 muestra el heatmap del eje \(Q^{CGE,Es}\) del cubo de resultados de MIGSA. El análisis del dendrograma superior mostró que los contrastes evaluados a partir de los datos del TCGA se asemejan a los respectivos resultados de HKds.

Como en la Figura 4.2, al analizar la Figura 4.3 se encontró un gran cluster con 11 de 12 experimentos contrastando Luminal A (marca LA), y otro cluster incluyendo los ocho experimentos contrastando el subtipo Basal con Luminal B ó Her2+ (marca Ba). Curiosamente, ambos clusters identificados contienen un cluster más específico, que agrupa únicamente experimentos transcriptómicos (marcas LA\(^T\) y Ba\(^T\)). Además, para cada contraste, se encontró un cluster que agrupaba solamente datos transcriptómicos, es decir, de microarreglos y secuenciación de ARN del TCGA, junto con los TECC de HKds (marcas T). Este resultado se encuentra en acuerdo con lo reportado por el “Grupo del Atlas del Genoma del Cáncer” (Network & others, 2012) para un subgrupo de sujetos Luminal B y Her2+. En concordancia con los resultados del “Grupo del Atlas del Genoma del Cáncer” para arrays de proteínas en fase inversa (Network & others, 2012), los experimentos proteómicos que contrastan Basal con Luminal B ó Her2+ se agruparon cerca de sus contrapartes transcriptómicas. Por otro lado, para los demás contrastes, la expresión de proteínas parece proporcionar información funcional complementaria a su contraparte transcriptómica.

Para cada subtipo, se obtuvieron los TET y los TEP, y los correspondientes GIE, y se analizaron en conjunto con los correspondientes PTF encontrados previamente. En la Tabla ??, para cada subtipo, se presenta el número de términos que se encuentran exclusivamente en PTF, TET y TEP, y los que se encuentran exclusivamente en cada intersección. Analizando estos los resultados, encontramos que la concordancia entre los PTF y su contraparte del TCGA es notable. Por ejemplo, de los 680 términos que pertenecen a los PTF (suma de celdas de la fila Total para columnas que involucran los PTF), el 38% también se encontró en TEP o TET de los datos del TCGA. En particular, los subtipos Basal y Luminal A fueron los únicos que exhibieron términos en común entre transcriptómica (PTF y TET) y proteómica (TEP). También se encontró que casi el 80% de esos términos en los TET también se encontraron en los PTF, evidenciando la alta consistencia de MIGSA para analizar los conjuntos de datos de transcriptómica de los HKds y TCGA. Por el contrario, y para todos los subtipos de CM, esos términos enriquecidos con datos proteómicos no se encontraron frecuentemente entre los TET. Este hallazgo sugiere que los TEP proporcionan información complementaria a la contraparte transcriptómica, como se reportó previamente (Meng et al., 2014). Sin embargo, será necesaria una investigación más profunda que incluya un mayor número de datasets proteómicos para poder definir un Perfil Proteómico Funcional.

Figure 4.3: Conjuntos de genes enriquecidos por análisis funcional para TCGA. Heatmap: conjuntos de genes en filas, experimentos en columnas; ordenados por distancia Jaccard y enlace promedio. Conjuntos de genes enriquecidos en rojo. Barras de columnas (de arriba abajo): salmón es “términos enriquecidos por consenso para los conjuntos de datos de Haibe-Kains”, en turquesa datos de TCGA; azul datos de secuenciación de ARN, verde de microarreglos, en salmón de iTRAQ; contraste (comparación entre pares de subtipos); subtipo1 y subtipo2, salmón es Luminal A, turquesa es Basal, púrpura es Her2+, y el verde es Luminal B. Clusters de dendrogramas: Se demarcan para los subtipos Basal (Ba) y Luminal A (LA), para experimentos transcriptómicos de Basal (Ba\(^T\)) y de Luminal A (LA\(^T\)), y para experimentos transcriptómicos específicos por contraste (T).

4.7.4 Exploración de los resultados de TCGA en conjunto con los PTF

Para cada subtipo, se analizaron los árboles de GO incluyendo los términos enriquecidos por cada categoría presente en las columnas de la Tabla ??. Como se señaló anteriormente, los términos enriquecidos únicamente en datos proteómicos se encontraron principalmente en los subtipos Basal y Her2+. Por otra parte, para Luminal A y Luminal B los árboles están dominados principalmente por ramas enriquecidas en transcriptómica.

Un hallazgo interesante al analizar los TEP fue la identificación del gen INPP4B como un GIE al contrastar los subtipos Basal y Her2+ con los Luminales. En particular, la sobre-expresión de INPP4B detectada en los tumores Luminales es consistente con su perfil de expresión reportado por el “Grupo del Atlas del Genoma del Cáncer” (Network & others, 2012) para datos de arrays de proteínas en fase inversa.

Para los sujetos Basales, los GIE tales como DHRS3, LDHB, HSD17B8, CBR4, CYB5A, y PHGDH fueron encontrados en ramas relacionadas con la actividad de la oxidoreductasa. De estos, los GIE de deshidrogenasa, DHRS3, LDHB, y el PHGDH fueron encontrados sobre-expresados en experimentos proteómicos en los cuales el subtipo Basal fue contrastado. Apoyando este hallazgo, la sobre-expresión de la proteína PHGDH ha sido previamente relacionada con tumores Basales (Gromova et al., 2015). Además, este gen se ha asociado con la proliferación celular, la migración y los procesos de invasión, y recientemente se ha identificado como un biomarcador inmunohistoquímico de cáncer (Song, Feng, Lu, Lin, & Dong, 2018).

El análisis de los árboles de Luminal B incluye una rama de TEP relacionada con la secreción hormonal. En esta rama, se encontraron genes tales como GATA3, IRS1, BAD, SYTL4, NDUFAF2, GPLD1, EGFR, FGB, FGG, y LYN. En particular, para el contraste proteómico Luminal B vs. Her2+, EGFR se encontró des-regulado para Luminal B, mientras que GATA3 estaba sobre-expresado, tal como se informa para arrays de proteínas en fase inversa por el “Grupo del Atlas del Genoma del Cáncer” (Network & others, 2012). Otra rama, relacionada con los ribosomas, fue encontrada para Luminal B, en la cual varios genes ribosómicos mitocondriales fueron encontrados como GIE: MRPL18, MRPL28, MRPL37, MRPL45, MRPL48, MRPL50, MRPL53, y MRPL54; todos ellos fueron encontrados sobre-expresados en Luminal B contrastándolo con Luminal A.

Como se muestra en la Tabla ??, 29 términos previamente identificados en los PTF de los subtipos de CM también fueron detectados en los respectivos TEP y TET. Entre estos, se encontraron términos que diferencian el subtipo Basal, como el término de FM de unión al distroglicano; y el de PB de regulación positiva de la vía de señalización del receptor de estrógeno intracelular. Para el subtipo Luminal A, se detectó una rama de CC relacionada con el cromosoma condensado, y términos de PB asociados a ciclos celulares mitóticos y meióticos. Veintidós GIE fueron encontrados contribuyendo al enriquecimiento de términos para PTF, TEP y TET. De ellos, 11 fueron identificados para el subtipo Basal (por ejemplo, AR, GATA3, FOXA1, NFIB, AGR2, CA12, y DSC2) y 11 para los tumores Luminal A (por ejemplo, TOP2A, CCNB1, CDC20, UBE2C, TTK, y KIF2C). Cabe destacar que todos estos genes se comportaron de acuerdo con las observaciones realizadas por el “Grupo del Atlas del Genoma del Cáncer” (Network & others, 2012).

4.8 Conclusiones

En este capítulo se presentó el paquete R MIGSA para el Análisis Masivo e Integrador de Conjuntos de Genes. Esta herramienta realiza un análisis funcional completo/integrador, masivamente sobre varias matrices de expresión de diferentes repositorios y fuentes ómicas. MIGSA facilita la búsqueda de patrones funcionales -que son denotados por conjuntos de genes y que explican fenotipos particulares- por medio de herramientas exploratorias y gráficas. MIGSA también permite a los usuarios indagar en conjuntos de genes particulares para buscar los genes específicos responsables del enriquecimiento funcional. Las capacidades analíticas y exploratorias de MIGSA como herramienta de análisis funcional masivo, en conjunto con el hecho de que a nuestro entender, es la única herramienta existente que permite la comparación de múltiples experimentos desde un punto de vista funcional, la convierten en un enfoque novedoso.

La utilidad de MIGSA quedó demostrada debido a la caracterización funcional de los subtipos de Cáncer de Mama (CM) definidos por Perou et al. La aplicación de MIGSA sobre 118 experimentos de microarreglos de CM resultó en la identificación de Perfiles Transcriptómicos Funcionales (PTF) y Genes Importantes de Enriquecimiento (GIE) para cada subtipo. La fiabilidad de nuestros hallazgos fue confirmada usando resultados previamente reportados en la literatura. La identificación de GIE con un comportamiento de expresión consistente a través de conjuntos de datos resulta útil para el descubrimiento de genes provocadores de enfermedades, y por ende, buenos blancos a atacar por drogas.

La capacidad de MIGSA de integrar múltiples conjuntos de datos ómicos se demostró aplicándolo a los datasets de CM del TCGA compuestos tanto de datos transcriptómicos (secuenciación de ARN y microarreglos) como proteómicos (iTRAQ). Los resultados obtenidos se utilizaron para identificar, para cada subtipo, el conjunto de Términos Enriquecidos por Transcriptómica (TET) y de Términos Enriquecidos por Proteómica (TEP). Al integrar y comparar estos resultados con los PTF encontrados previamente, se descubrió una alta consistencia de patrones funcionales entre los términos enriquecidos en datos de TCGA y en los PTF. Dado que MIGSA es una herramienta tan flexible, el procedimiento presentado puede ser replicado tanto en otros subtipos de CM como en otras enfermedades con el fin de explorar perfiles funcionales. El paquete MIGSA se encuentra libremente disponible -bajo una licencia GPLv2- en el repositorio Bioconductor https://bioconductor.org/packages/MIGSA/, y cuenta con más de 3.000 descargas según las estadísticas del repositorio https://bioconductor.org/packages/stats/bioc/MIGSA/ desde su primera versión, en Marzo de 2017.