El esquema conceptual de GREIN se muestra en la Fig. 1. Los conjuntos de datos individuales de RNA-seq son procesados por el pipeline GREP2 y almacenados localmente como R Expression Sets. El usuario puede acceder y analizar conjuntos de datos preprocesados a través de la interfaz gráfica de usuario (GUI) de GREIN o enviar para su procesamiento conjuntos de datos que aún no han sido procesados. Los flujos de trabajo impulsados por la GUI facilitan el examen y la visualización de los datos, el análisis estadístico, la construcción de firmas transcripcionales y la interpretación de biología de sistemas de los genes expresados diferencialmente (DE). Tanto GREIN como el back-end pipeline (GREP2) están escritos en R y publicados como contenedor Docker y paquete R respectivamente. Las interfaces gráficas de usuario para GREIN se implementan en Shiny16, un marco web para construir aplicaciones web dinámicas en R. La instancia web en https://shiny.ilincs.org/grein se despliega a través de un robusto enjambre Docker de servidores Shiny de carga equilibrada. La infraestructura completa de GREIN, incluyendo la tubería de procesamiento se despliega a través de contenedores Docker.
Flujo de trabajo esquemático de GREP2, interfaz web y salidas de GREIN. Los conjuntos de datos GEO se procesan sistemáticamente utilizando el pipeline de GREP2 y se almacenan en la biblioteca de conjuntos de datos back-end. Los flujos de trabajo de GREIN impulsados por la GUI facilitan el análisis y la visualización de los conjuntos de datos procesados.
Los flujos de trabajo impulsados por la GUI en GREIN facilitan los escenarios típicos de reutilización de los datos de RNA-seq, como el examen de las medidas de control de calidad y la visualización de los patrones de expresión en todo el conjunto de datos, el tamaño de la muestra y el análisis de potencia con el fin de informar el diseño experimental de futuros estudios, la expresión genética diferencial estadística, el enriquecimiento de la lista de genes y el análisis de redes. Además de la comparación estándar de dos grupos, el módulo de análisis de expresión génica diferencial también admite el ajuste de un modelo lineal generalizado que tiene en cuenta las covariables o los efectos de lote. Las herramientas interactivas de visualización y exploración implementadas incluyen el análisis de conglomerados, los mapas térmicos interactivos, el análisis de componentes principales (PCA), la incrustación estocástica de vecinos t-distribuida (t-SNE), etc. (Tabla suplementaria S1). El usuario también puede buscar anotaciones ontológicas de muestras de RNA-seq humano y conjuntos de datos proporcionados por el proyecto MetaSRA13. Cada muestra de RNA-seq humana procesada está etiquetada con el mapeo MetaSRA de ontologías biomédicas, incluyendo Disease Ontology, Cell Ontology, Experimental Factor Ontology, Cellosaurus y Uberon. La interpretación biológica de las expresiones genéticas diferenciales se ve facilitada por los enlaces directos a otras herramientas en línea para realizar los típicos análisis post-hoc, como la lista de genes y el análisis de enriquecimiento de vías y el análisis de redes de genes expresados diferencialmente (DE). La conexión con estos servicios web analíticos se implementa enviando la firma de expresión génica diferencial (es decir, la lista de cambios medios en la expresión génica y los valores p asociados para todos los genes regulados arriba/abajo analizados) a iLINCS17 (Integrative LINCS). iLINCS también proporciona el análisis de conectividad de las firmas para las firmas del Mapa de Conectividad L1000 recientemente publicado18. Las instrucciones detalladas paso a paso sobre los flujos de trabajo de análisis de GREIN se proporcionan en el material suplementario y en la sección de «Ayuda» de GREIN.
Funcionalidades clave
Búsqueda o envío para su procesamiento
El usuario puede buscar un conjunto de datos GEO ya procesado en la casilla «Search for GEO series (GSE) accession» o enviar un conjunto de datos para su procesamiento si el conjunto de datos no está ya procesado (Supplementary Fig. S2). En este momento, la gran mayoría de los conjuntos de datos GEO de ARN-seq de humanos, ratones y ratas han sido preprocesados y el envío por parte del usuario de conjuntos de datos GEO para su procesamiento sólo será necesario ocasionalmente. El usuario puede comprobar el estado de procesamiento del conjunto de datos solicitado en la pestaña «Consola de procesamiento» (Supplementary Fig. S3). Otras opciones de búsqueda incluyen la búsqueda de palabras clave a través de los metadatos de los conjuntos de datos y la búsqueda de muestras a través de ontologías biomédicas a través de las anotaciones ontológicas de MetaSRA.
Explore dataset
GREIN permite el acceso tanto a los datos brutos como a los normalizados (recuentos por millón y transcripción por millón) a nivel de genes y transcripciones. GREIN viene con varias herramientas interactivas y personalizables para visualizar los patrones de expresión, tales como mapas de calor interactivos de genes y muestras agrupadas, gráficos de densidad para todas las muestras o un subconjunto de ellas, análisis de variabilidad entre grupos y dentro de ellos mediante análisis de reducción de la dimensionalidad en 2D y 3D y visualizaciones como PCA y t-SNE (Fig. 2). El usuario también puede visualizar el perfil de expresión de cada gen por separado (Fig. Suplementaria S6).
Tramas de análisis exploratorio en GREIN. (A) El mapa de calor de correlación muestra una mayor correlación dentro de las líneas celulares y una baja correlación entre las líneas celulares. En general, las correlaciones altas dentro de cada línea celular indican una alta calidad de los perfiles transcripcionales. (B) Agrupación jerárquica basada en la correlación de Pearson de los 500 genes más variables según la desviación absoluta media como medida de variabilidad. Los datos están normalizados y centrados en la media. (C) Gráfico de análisis de componentes principales tridimensional de las líneas celulares. (D) El gráfico t-SNE bidimensional de la condición de tratamiento y la línea celular muestra una clara separación de las líneas celulares, y luego las fracciones de ARN que indican dos fuentes dominantes de la variabilidad entre los perfiles de ARN-seq.
Control de calidad
La calidad de los datos de ARN-seq en los repositorios públicos sigue siendo un problema importante. En un estudio reciente de Deelen et al.19, más de la mitad de las 65.000 muestras públicas de RNA-seq procesadas tuvieron que ser eliminadas debido a problemas de control de calidad. En lugar de eliminar las muestras, GREIN proporciona un informe completo de control de calidad (QC) de los datos de la secuencia en bruto y el mapeo de la secuencia para cada muestra (Supplementary Fig. S7), y permite al usuario tomar una decisión sobre qué muestras deben ser excluidas de los análisis posteriores.
Análisis de potencia estadística
El módulo de análisis de potencia en GREIN facilita el cálculo y la visualización de la potencia estadística de la detección de genes expresados diferencialmente en futuros estudios que utilizan muestras biológicas similares. La estimación del tamaño adecuado de la muestra para futuros estudios con muestras biológicas similares es a menudo el factor clave que motiva el reanálisis de los datos de RNA-seq. El análisis de potencia también facilita el análisis post-hoc de las tasas de falsos negativos en el conjunto de datos actual. La falta de potencia estadística y las diferencias de potencia estadística entre los genes pueden producir resultados falsos negativos que lleven a conclusiones erróneas20. El segmento «Curva de potencia» proporciona estimaciones de potencia para diferentes números de muestras basadas en un solo gen (Fig. 3A). El usuario puede modificar los valores por defecto de los parámetros. El gráfico «Detectabilidad de los genes» visualiza la estimación de potencia de cada uno de los genes basándose en los grupos seleccionados y en la dispersión de los genes (Fig. 3B). La cobertura media de los genes se traza contra su variabilidad biológica y se muestran en dos conjuntos basados en su estado de detectabilidad (potencia ≥0,8 y potencia < 0,8).
Análisis de potencia para evaluar los cambios transcripcionales en la línea celular MCF10A no maligna. (A) Estimaciones de potencia basadas en un solo gen para diferentes números de muestras en cada grupo con un cambio mínimo de 2 pliegues y una significación estadística α = 0,01. (B) Detectabilidad por genes en el plano log2CPM-BCOV con FDR ≤0,1 y dos muestras en cada grupo.
Expresión génica diferencial
Crear e interpretar la firma de expresión génica diferencial es un escenario de análisis típico en los experimentos de RNA-seq. Con GREIN, el usuario puede crear una firma comparando la expresión génica entre dos grupos de muestras con o sin ajustes por covariables experimentales o efectos de lote. GREIN puede manejar diseños experimentales complejos proporcionando la flexibilidad de reorganizar grupos y subgrupos o seleccionar muestras específicas. La firma de expresión diferencial puede visualizarse a través de gráficos interactivos que incluyen un mapa de calor de los principales genes desregulados diferencialmente (Fig. suplementaria S15) clasificados por la tasa de falsos descubrimientos (FDR), un gráfico de cambio de pliegue logarítmico frente a la expresión media logarítmica (MA) (Fig. suplementaria S16) y un gráfico de detectabilidad de genes (Fig. suplementaria S17). La firma de expresión diferencial, con o sin tener en cuenta los resultados potencialmente falsos negativos, puede exportarse directamente a iLINCS para el análisis de enriquecimiento y conectividad.
Caso de uso: Análisis de la regulación transcripcional y traslacional de la hipoxia en líneas celulares de cáncer de mama no maligno y triple negativo
Demostramos el uso de GREIN volviendo a analizar unos datos de RNA-seq publicados recientemente en GEO (GSE104193). Sesé et al.21 examinaron la regulación transcripcional y traslacional del subtipo de cáncer de mama triple negativo (TNBC) refractario a las hormonas bajo una combinación de tratamiento con hipoxia e inhibidores de mTOR (diana mecánica de la rapamicina). En concreto, los autores analizaron los perfiles de expresión de células TNBC (MDA-MB-231) y epiteliales de mama no malignas (MCF10A) expuestas a condiciones de normoxia (21% de O2) e hipoxia (0,5% de O2) y/o tratadas con un inhibidor de mTORC1 y -2 PP242. Cada una de las muestras fue secuenciada para el ARNm total (T) y unido a polisomas (P). El conjunto de datos contiene 32 muestras, que representan dos réplicas biológicas para cada combinación de línea celular, nivel de oxígeno, estado de tratamiento y fracción de ARNm.
El análisis exploratorio del conjunto de datos procesados en GREIN (Fig. 2) muestra que la fuente más fuerte de variación entre las muestras proviene de las diferencias entre las dos líneas celulares. Esto se ve reforzado por el análisis de correlación de los perfiles de expresión completos (Fig. 2A), la agrupación jerárquica de los 500 genes más variables basada en la desviación absoluta media (Fig. 2B), el gráfico PCA 3D de las muestras (Fig. 2C) y el gráfico t-SNE 2D (Fig. 2D). Además, las altas correlaciones entre los perfiles de expresión para la misma línea celular (Fig. 2A) indican una buena relación señal-ruido en las mediciones de la expresión génica. La subestructura adicional de los datos indicada por el gráfico 2D t-SNE se ha examinado pintando las muestras según diferentes atributos (Fig. Suplementaria S1). Este análisis reveló que las separaciones dentro de cada línea celular son inducidas por diferentes fracciones de ARNm y luego por diferencias entre las condiciones experimentales.
A continuación, utilizamos GREIN para realizar un análisis de potencia estadística basado en el patrón de variabilidad biológica observado en este conjunto de datos. Consideramos los perfiles transcripcionales de cada línea celular expuesta a la hipoxia y tratada con o sin PP242, lo que conduce a cuatro comparaciones. Asumiendo una diferencia de expresión de al menos dos veces entre los grupos, a la significación estadística de α = 0,01, y con sólo dos réplicas en cada grupo, el poder estadístico de un gen para ser detectado como diferencialmente expresado es inferior a 0,55 en todas las comparaciones (Tabla 2). Nuestro análisis indica que se necesitarían cuatro réplicas por grupo para alcanzar un 80% de potencia detectando un cambio de dos veces en la expresión (Tabla 2 y Fig. 3A). En un experimento típico de RNA-seq, una profundidad de secuenciación de 20-30 millones es suficiente para cuantificar la expresión de casi todos los genes4,22 , lo que también es evidente en este conjunto de datos. También evaluamos el poder estadístico de cada gen para ser detectado como diferencialmente expresado a partir del gráfico de «Detectabilidad de genes». Los valores medios del logaritmo de los recuentos por millón (CPM) de los genes se trazaron frente al coeficiente de variación biológica (BCOV) del gen y se calculó la potencia para los genes correspondientes (Fig. 3B). Se utilizó una tasa de falsos descubrimientos controlada de 0,05 y un porcentaje esperado de verdaderos positivos del 10% para estimar la significación estadística. Definimos que un gen es detectable como diferencialmente expresado en condiciones de hipoxia si su potencia es de 0,8 o superior. Como se esperaba, existe una relación inversa entre el BCOV y la potencia (Fig. 3B). Además, la potencia para detectar la expresión diferencial de un gen aumenta con un mayor log CPM o tamaño del efecto.
Uno de los objetivos del estudio era analizar los cambios transcripcionales en condiciones hipóxicas y normoxicas con y sin tratamiento con PP242 en las líneas celulares MCF10A y MDA-MB-231. Creamos firmas transcripcionales de hipoxia e hipoxia + PP242 en el ARNm total mediante el análisis de expresión diferencial entre las muestras de hipoxia e hipoxia + PP242, respectivamente, frente a las muestras de control, mientras se ajustaba el efecto de lote tratando la «réplica» como una covariable, para cada línea celular por separado. Encontramos un mayor número de genes expresados diferencialmente (DE) en las líneas celulares MCF10A en comparación con MDA-MB-231 tanto en hipoxia como en hipoxia + PP242 (Fig. 4A), lo que indica que quizás la línea celular tumoral está mejor equipada para hacer frente a la hipoxia. Este análisis también mostró que la mayoría de los genes no expresados diferencialmente tampoco son detectables, lo que indica que pueden representar resultados falsos negativos. Esto está en consonancia con el análisis de potencia que muestra que se necesitarían 4 muestras por grupo para identificar de forma consistente los genes expresados diferencialmente con un BCOV medio. Para identificar genes de menor expresión se necesitaría un tamaño de muestra aún mayor.
Expresión diferencial y detectabilidad de los genes. (A) El número de genes (escala log10) no expresados diferencialmente y no detectables (NDE&NDT), expresados diferencialmente (DE), y no expresados diferencialmente pero detectables (NDE&DT) en las comparaciones con normoxia para la fracción de ARNm total. Llamamos a un gen detectable (DT) si su potencia ≥0,8 y diferenciable si FDR < 0,05. (B) El gráfico de detectabilidad de genes para la primera comparación (MCF10A e hipoxia) que visualiza la lista de genes mencionada anteriormente junto con sus respectivos cambios de pliegues (FC).
Para interpretar los genes expresados diferencialmente en términos de vías biológicas afectadas, enviamos las firmas de expresión génica diferencial de la hipoxia a herramientas de enriquecimiento en línea (DAVID23, ToppGene24, Enrichr25 y Reactome26) a través de iLINCS. Las firmas enviadas incluían una lista combinada de genes DE y NDE&DT que representaban probables verdaderos positivos y verdaderos negativos. Los genes se seleccionaron basándose en un límite de 0,7 y 0,01 para la potencia estadística y el FDR respectivamente. La figura 5 ilustra los resultados de enriquecimiento obtenidos de ToppGene para la firma de hipoxia MCF10. Los genes enriquecidos significativamente (FDR < 0.05) Las 10 principales categorías de ontología génica (GO) de ToppGene y de la herramienta de anotación funcional DAVID incluyen la respuesta a la hipoxia, la respuesta a la disminución de los niveles de oxígeno, la angiogénesis, la regulación de la proliferación celular, el proceso de oxidación-reducción y la respuesta al estímulo abiótico que son comunes en ambas líneas celulares (Tabla suplementaria S2 y Tabla suplementaria S3). La mayoría de estas categorías son consistentes con el estudio original. Además, la suite ToppGene identificó la red del factor de transcripción inducido por la hipoxia (HIF-1-alfa) que se activó en ambas líneas celulares (Tabla Suplementaria S4 y Tabla Suplementaria S5).
Reducción de algunas de las categorías significativas de vías y ontología genética (GO) de ToppGene a través de iLINCS. Estas categorías se encuentran en la comparación entre hipoxia y normoxia en la línea celular MCF10A utilizando una lista combinada de genes DE y NDE&DT. La línea vertical roja es el punto de corte seleccionado de 0,05.
Finalmente, utilizamos la conexión GREIN con iLINCS para «conectar» la firma cargada con las firmas de consenso LINCS27 (CGS) de derribo de genes18. Encontramos 3.727 firmas de consenso LINCS que estaban significativamente (pValue < 0,05) conectadas con nuestra firma cargada. Los genes objetivo de las 100 firmas más conectadas se seleccionaron para un análisis de enriquecimiento adicional. Encontramos la respuesta celular a la hipoxia y la regulación del Factor Inducible por Hipoxia (HIF) por el oxígeno en la lista de las 10 principales vías activadas en ambas líneas celulares (Tabla Suplementaria S6 y Tabla Suplementaria S7). Si bien este análisis arroja categorías funcionales enriquecidas similares a las del análisis de enriquecimiento inicial, complementa el análisis original al implicar a varios genes diana que no se expresan de forma diferencial, aunque están lo suficientemente expresados como para ser detectables según nuestros análisis de potencia. La unión de estos dos resultados implica a estos genes como potenciales reguladores de alto nivel de la respuesta a la hipoxia.