PEC 2: Análisis de microarrays
Aitor Moruno Cuenca
2023-05-22
1 Introducción
El ácido arquidónico, un ácido graso pertenenciente a la familia de los omega-6, ha mostrado inducir genes que están bajo control del factor de transcripción NFkappaB, según análisis de expresión diferencial con tecnología de microarrays (Affymetrix). Estudios previos ([1]) han revelado que la adición de ácido arquidónico lleva a una inducción significativa de varios genes, entre ellos: COX2, IKBA, NFKB, GMCSF, IL1B, CXCL1, TNFA, IL6, LTA, IL8, PPARG y ICAM1.
Además, la vía de señalización de PI3K se activa de forma significativa cuando se añade ácido arquidónico, esto sugiere que hay una relación entre los ácidos grasos y las vías de señalización celular implicadas en proliferación y crecimiento celular. A pesar de esto, los mecanismos e implicaciones biológicas no están claros. Estos hallazgos pueden tener un papel importante en el cáncer de próstata humano.
1.1 Objetivos
Este estudio tiene como objetivo principal replicar y profundizar en hallazgos previos, utilizando tecnología de microarrays (Affymetrix) para explorar la inducción de genes por ácidos arquidónicos. En concreto este trabajo conseguirá:
- Usar los datos del estudio previo (Hughes-Fulford M, Li C, Boonyaratanakornkit JB, Sayyah S. Omega-6 fatty acids, arachidonic acid (AA) activates PI3K signaling and induces gene expression in prostate cancer) descargándolos directamente de la base de datos Gene Expression Omnibus (GEO).
- Realizar un control de calidad de los datos crudos que proporciona el estudio, y una normalización y filtraje de estos.
- Validar la inducción de genes (mediante análisis de genes diferencialmente expresados), incluyendo COX2, IKBA, NFKB, GMCSF, IL1B, CXCL1, TNFA, IL6, LTA, IL8, PPARG y ICAM1, al usar ácido arquidónico comparado con un grupo de control.
- Llevar a cabo un análisis de significación biológica que permita replicar los resultados de estudios anteriores, en concreto determinar si los genes diferencialmente expresados están involucrados en procesos de crecimiento y proliferación celular, y activación del factor de transcripción NFkappaB.
2 Métodos
2.1 Database y softwares
Los datos con los que se trabajarán a lo largo de este estudio se han descargado directamente de GEO, descargando todos los archivos .CEL relacionados y los datos fenotípicos correspondientes. Concretamente, los datos han sido extraídos de: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE3737, el código GEO Serie es GSE3737.Se usó la plataforma de microarrays GPL96 Affymetrix Human Genome U133A Array para la obtención de tales datos.
El estudio contiene un total de 8 microarrays, es decir, trabajaremos con 8 archivos .CEL. Una vez descargados los archivos y anotaciones directamente de GEO, se leyeron los archivos mediante librerías de Bioconductor como Biobase y affy. La información del fenotipo de los arrays solo contiene una covariable que es la condición experimental determinada por el tratamiento.
A lo largo del trabajo se usará el software R (versión 4.2.3) y las librerías para análisis de datos ómicos que están integradas en Bioconductor.
2.2 Control de calidad
Para el análisis de calidad de los microarrays de Affymetrix, se utilizaron las librerías: limma, CCl4 y affyPM. Para cada array se generó un gráfico de la distribución de la señal usando la función hist y un diagrama de caja con la función boxplot. Estos gráficos permiten evaluar visualmente la distribución y variabilidad de la señal de expresión de genes en cada array.
Además, se realizó un Análisis de Componentes Principales (ACP) para visualizar la variabilidad en los datos de la expresión génica entre las diferentes muestras, este análsis se realizó con la función prcomp.
Finalmente, se llevó a cabo un análsis de agrupación jerárquica con la función hclust usando distancia euclideana y el método de agrupación por promedio.
2.3 Preprocesamiento
El preprocesamiento incluyó la normalización y el filtrado de genes. Para la normalización, se utilizó el método Robust Multi-array Average (RMA), implementado en la función rma. Este método permite ajustar la intensidad de la señal de cada array para reducir la variabilidad técnica entre los arrays.
Para el filtrado de genes, se utilizó la librería genefilter y la función nsFilter, que permite eliminar genes que no cumplen con criterios de calidad y variabilidad. El proceso finalizó guardando los datos normalizados y filtrados listos para el análisis en un objeto del tipo expressionSet
2.4 Análisis de expresión diferencial
Se detallará primero el diseño del experimento. Se obtuvieron células de cáncer de próstata PC-3, las cuales se incubaron con 5 \(\mu\)g/mL de ácido arquidónico en medio de RPMI (un tipo de cultivo celular) que contenía mg/mL de albúmina durante 2h. Por otro lado, el segundo grupo de células se trataron solo con albúmina y se mantuvieron sin activar. Es decir, trataremos en este estudio con un solo factor de tratamiento con dos niveles que diremos: PC-3 archidonic acid treated 2hr (el primer grupo de células) y PC-3 untreated (el segundo grupo). En resumen, el diseño experimental incluye la realización de microarrays tanto para muestras tratadas como no tratadas en cuatro réplicas.
Para llevar a cabo el anàlisis de expresión diferenical, se usaron las librerías de R limma y edgeR. Se definió la matriz de diseño de la siguiente forma:
\[ \mathbf{X} = \begin{bmatrix} 0 & 1 \\ 0 & 1 \\ 0 & 1 \\ 0 & 1 \\ 1 & 0 \\ 1 & 0 \\ 1 & 0 \\ 1 & 0 \end{bmatrix} \]
Donde las filas representan los arrays y las columnas el tratamiento, siendo la columna 1 el tratamiento con ácido arquidónico y la columna 2 el control. Los cuatro primeros arrays pertenecen al grupo control y los cuatro últimos al grupo de tratamiento con ácido arquidónico. Tal y como muestra la matriz de diseño, se eliminó el término de intercepción y solo se incluyó el factor tratamiento.
Se estimó un modelo lineal del tipo \(\mathbf{Y} = \beta'\mathbf{X} + \varepsilon\) a través de la función lmFit. Se definió la siguiente matriz de contraste:
\[ \textbf{C} = \begin{bmatrix} -1 \\ 1 \end{bmatrix} \]
Pues solo se realizará un test (para cada gen) con un factor de dos niveles, ya que el objetivo será contrastar la expresión diferencial del tratamiento vs control.
Posteriormente, se aplicó el método de Empirical Bayes a estos contrastes con la función eBayes para estabilizar las varianzas estimadas. Los resultados se extrajeron aplicando el ajuste de multiplicidad de Benjamini y Hochberg para contrar la tasa de descubrimiento falso (FDR) en los p-valores.
Los genes diferencialmente expresados se seleccionaron utilizando un umbral de significación estadística de p-valor corregido por FDR menor a 0.05 y un fold change absoluto mayor a 1, de manera que los genes seleccionados sean tanto estadística como biológicamente significativos en su expresión diferencial.
Por último, se generó un Volcano plot con la función volcanoplot, que muestra los genes diferencialmente expresados resaltados en el gráfico en función de sus valores de fold change y p-valores ajustados. El eje x representa los logaritmos de los fold changes, mientras que el eje y representa los logaritmos negativos de los p-valores ajustados. La línea vertical indica el umbral de fold change 1 o -1.
2.5 Anotaciones
Para realizar las anotaciones de los genes, se usaron las librerías mygene y dply. Se utilizó la base de datos hgu133a.db que proporciona anotaciones y datos de la plataforma de microarrays Affymetrix Human Genome U133A. Contiene información sobre las sondas (probesets) utilizadas en este chip. Los genes se anotaron utilizando las anotaciones disponibles en la base de datos, concretamente usando el identificador de “EntrezID”.
2.6 Análisis de significación biológica
2.6.1 Análisis de enriquecimiento de genes (Gene Enrichment Analysis)
Para realizar el análisis de enriquecimiento de genes, se utilizó la librería GOstats. Se llevó a cabo un análisis de sobre-representación utilizando la función enrichGO de la anterior librería. En este análisis, se utilizaron los genes seleccionados previamente y se ejecutó un enriquecimiento de GO (Gene Ontology) en la ontología biológica (BP). Para el ajuste de los p-valores se usó el método de Benjamini y Hochberg y un q-value máximo de 0.25 para determinar los términos enriquecidos significativos.
Los resultados del análisis de enriquecimiento GO se guardaron en un data frame y se visualizaron mediante un gráfico de puntos dotplot y un gráfico de jerarquías goplot. Además, se generó un gráfico de red de genes cnetplot para visualizar las interacciones entre los genes enriquecidos.
2.6.2 Análisis de enriquecimiento de conjunto de genes (Gene Set Enrichment Analysis)
Para el análisis de enriquecimiento de conjuntos de genes (GSEA), se utilizó la librería clusterProfiler. Primero, se ordenaron los genes ba´sandose en el valor asboluto del fold change para eliminar duplicados con el menor valor absoluto.
Luego, se reordenaron los genes basándose en el fold change para realizar el análisis GSEA utilizando la función gseKEGG. En este análisis, se utilizó el conjunto de genes reordenado, se seleccionó el organismo hsa (Homo sapiens) y se utilizó la base de datos de anotación kegg. Se estableció un umbral máximo de 0.05 para el p-valor máximo y se aplicó el método de ajuste de Benjamini y Hochberg.
Los resultados del análisis de GSEA se presentaron en un data frame utilizando la función kable de la librería kableExtra.
3 Resultado
3.1 Control de calidad y preprocesamiento
Se evaluará la calidad de los datos obtenidos de cada uno de los ocho arrays conjuntamente, de manera que podremos determinar la fiabilidad y validez de los datos obtenidos, asegurándonos que minimizamos los efectos no biológicos y los artefactos técnicos que puedan afectar a la interpretación de los resultados. Antes de empezar con el análisis, se mostrará una tabla con el array y su condición experimental:
## array treatment
## 1 GSM86079.CEL "PC-3, untreated"
## 2 GSM86080.CEL "PC-3, untreated"
## 3 GSM86081.CEL "PC-3, untreated"
## 4 GSM86082.CEL "PC-3, untreated"
## 5 GSM86083.CEL "PC-3, arachidonic acid treated 2hr"
## 6 GSM86084.CEL "PC-3, arachidonic acid treated 2hr"
## 7 GSM86085.CEL "PC-3, arachidonic acid treated 2hr"
## 8 GSM86086.CEL "PC-3, arachidonic acid treated 2hr"
Este gráfico nos muestra la distribución de intensidades de la señal de los arrays en cada experimento. El eje x representa las intensidades de señal, en escala logarítmica, y el eje y representa la frecuencia o el número de sondas que tienen tal intensidad de señal.
Observemos como hay tres arrays cuya distribución parece agruparse: GSM86084, GSM86079, GSM86083, los arrays GSM86084 y GSM86083 presentan la misma condición experimental (tratados con ácido arquidónico) mientras que GSM86079 tiene otra condición exprimental (control), de manera que parece ser que los arrays no se agrupan de forma natural basandonos únicamente en la condición experimental de tratados/no tratados. Además la distribución de las señales presentan asimetría hacia la derecha, no parecen ser una distribución normal.
A partir de los boxplots de la expresión de los arrays, vemos que tanto la variabilidad como la tendencia central no están normalizados a lo largo de los arrays, por tanto, será necesaria una normalización de los datos para corregir estos sesgos. Tampoco se observa una tendencia clara sobre la tendencia central y/o variabilidad en función de la condición experimental.
En el análisis de PCA, las dos dimensiones explican un 85,4% de la variabilidad de los datos, siendo la primera dimensión la que más variabilidad explica en un 72,7%. Como la dimensión 1 es la que más variabilidad explica, cabría esperar que esta fuese naturalmente debida a las condiciones experimentales, es decir, en función de si el array ha sido tratado con ácido arquidónico o no. En línea con el análisis de la distribución de la señal, vemos que los arrays no se agrupan naturalmente por su condición experimental pues se observan dos grupos de arrays (a partir de la dimensión 1); grupo 1: GSM86082, GSM86080, GSM86085, GSM86086; grupo 2: GSM86081, GSM86083, GSM86084, GSM86079, estos grupos contienen arrays de ambas condiciones experimentales, por tanto, los arrays no se agrupan en función de su condición experimental, esto podría producir un efecto batch causando sesgos en los resultados, debemos normalizar y filtrar los arrays.
Para confirmar que realmente existen dos grupos que separan los arrays, realizaremos un cluster jerárquico, se muestra el dendograma siguiente:
Confirmamos que los grupos que hemos podido deducir del ACP agrupan arrays con condiciones experimentales diferentes, siendo dos grupos de cuatro arrays cada uno; grupo 1: GSM86082, GSM86080, GSM86085, GSM86086; grupo 2: GSM86081, GSM86083, GSM86084, GSM86079.
Se realiza la normalización mediante Robust-Multiarray average y el filtrado, mostramos por pantalla que el objeto resultante es del tipo expressionSet.
## Background correcting
## Normalizing
## Calculating Expression## [1] "ExpressionSet"
## attr(,"package")
## [1] "Biobase"Veamos el resultado de la expresión después de la normalización y filtrado.
La señal de la distribución se ha corregido, pues ahora sigue una distribución normal y todos los arrays tienen una distribución que se asemeja, esto también queda plasmado en los boxplots pues la tendencia central y medidas de variabilidad de los arrays es suficientemente semejante. Por otro lado, del ACP y clúster jerárquico vemos que siguen habiendo dos grupos para clasificar los arrays, pero estos siguen sin estar divididos según su condición experimental pues por ejemplo los arrays GSM86082 y GSM86080 (controles) son clasificados en el mismo grupo que GSM86085 y GSM86086 (tratados), esto puede indicar que podría haber algún problema técnico con las muestras en el proceso de elaboración o algún efecto batch que no podemos identificar por falta de más covariantes e información de la muestra.
3.2 Análisis de expresión diferencial
Se usa un modelo lineal para identificar los genes diferencialmente expresados, ya sea por sobre-regulación (upregulation) o sub-regulación (downregulation), entre las diferentes condiciones experimentales (tratados con ácido arquidónico vs control). Se muestra a continuación la matriz de diseño en R:
## Treated Control
## 1 0 1
## 2 0 1
## 3 0 1
## 4 0 1
## 5 1 0
## 6 1 0
## 7 1 0
## 8 1 0
## attr(,"assign")
## [1] 1 1
## attr(,"contrasts")
## attr(,"contrasts")$group
## [1] "contr.treatment"Se ajusta un modelo lineal de mínimos cuadrados ponderados con lmFit, se muestran los coeficientes estimados de los 7 primeros genes:
## Treated Control
## 203441_s_at 7.732894 7.719993
## 212607_at 7.567524 7.425345
## 210104_at 8.486079 8.196426
## 208385_at 5.173028 5.319496
## 204212_at 8.042829 7.991760
## 209028_s_at 4.957285 4.756992El nombre de la fila es el probeset, treated el coeficiente del modelo lineal asociado a la condición de tratamiento y al control respectivamente.
Como se trata de un diseño de un factor, solo hay un contraste posible que será tratados vs control para cada uno de los genes, el objetivo a nivel biológico es identificar aquellos genes diferencialmente expresados en células de cáncer de próstata tratadas con ácido arquidónico contra aquellas que actúan como control pues solo han sido tratadas con albúmina sin activación, por tanto, el único contraste estadístico será tratados vs control. Se muestra la matriz de contrastes de R:
## Contrasts
## Levels Treated - Control
## Treated 1
## Control -1Y el estadístico de los 7 primeros contrastes realizados:
## Contrasts
## Treated - Control
## 203441_s_at 0.01290125
## 212607_at 0.14217910
## 210104_at 0.28965313
## 208385_at -0.14646761
## 204212_at 0.05106840
## 209028_s_at 0.20029284Después de aplicar ajuste por multiplicidad y Empirical Bayes, mostramos los probesets de aquellos genes que están diferencialmente expresados tanto por presentar un p-valor ajustado < 0.05 (significación estadística) como por presentar significación biológica (absolute log fold-change > 1).
## logFC AveExpr t P.Value adj.P.Val B
## 209774_x_at 5.072934 8.474089 18.543219 5.327165e-09 1.417918e-05 7.792097
## 204470_at 4.211596 8.735476 18.283695 6.100471e-09 1.417918e-05 7.741081
## 211506_s_at 3.505208 8.032616 18.097937 6.729560e-09 1.417918e-05 7.703570
## 210229_s_at 1.957622 7.717278 13.252182 1.303480e-07 2.059824e-04 6.332773
## 207850_at 4.537659 7.163907 10.833333 8.472127e-07 1.071046e-03 5.222143
## 218810_at 1.543102 7.223723 9.820219 2.070937e-06 2.181733e-03 4.626485
## 201502_s_at 2.945333 10.208619 9.474476 2.859751e-06 2.582355e-03 4.401517
## 211668_s_at 1.779871 10.169772 7.682947 1.803426e-05 1.336413e-02 3.024788
## 205476_at 1.969088 5.757488 7.634990 1.902818e-05 1.336413e-02 2.982467
## 39402_at 1.722010 7.714913 6.965201 4.134850e-05 2.613639e-02 2.357546
## 202638_s_at 1.624930 9.043722 6.709451 5.638666e-05 3.240183e-02 2.101482
## 201631_s_at 1.368266 10.150996 6.586036 6.567999e-05 3.459693e-02 1.974303
## 202644_s_at 2.268121 8.079448 6.497004 7.340835e-05 3.569340e-02 1.881067
## 210538_s_at 2.207823 5.364683 6.194624 1.079254e-04 4.872830e-02 1.554903
## ID
## 209774_x_at 209774_x_at
## 204470_at 204470_at
## 211506_s_at 211506_s_at
## 210229_s_at 210229_s_at
## 207850_at 207850_at
## 218810_at 218810_at
## 201502_s_at 201502_s_at
## 211668_s_at 211668_s_at
## 205476_at 205476_at
## 39402_at 39402_at
## 202638_s_at 202638_s_at
## 201631_s_at 201631_s_at
## 202644_s_at 202644_s_at
## 210538_s_at 210538_s_atObservemos como todos los genes diferencialmente expresados están sobre-regulados, pues su log fold change es positivo y superior a 1. Hay un total de 14 genes diferencialmente expresados. Se muestra a continuación un Volcano plot que resume gráficamente estos resultados.
En el eje x tenemos la magnitud del cambio en log2 fold change si esta es superior a 1 entonces el gen está sobre-regulado, mientras que si es menor a -1 el gen está sub-regulado. El eje y indica la significación estadística de tal magnitud, aquellos probeset en color azul resaltados en el gráfico tienen una magnitud estadísticamente significativa y, por tanto, están diferencialmente expresados en función de la condición experimental que haya. En definitiva, hay 14 probesets que indican una expresión diferente si una célula de cáncer de próstata es tratada con ácido arquidónico o con un control.
3.3 Análisis de significación biológica
Haremos un análisis de la significación biológica o análisis de enriquecimiento, para así poder entender el significado biológico de los genes diferencialmente expresados. Después de anotar los genes, mediante identificador de EntrezID, se usarán los identificadores para mapearlos a vías biológicas conocidas, que pueden ser procesos metabólicos, rutas de señalización, o cualquier otra colección de genes para realizar pruebas estadísticas y determinar si esas vías están enriquecidas en los genes de interés. Es decir, veremos si hay más genes de interés en una vía de lo que cabría esperar por puro azar.
3.4 Análsis de enriquecimiento de genes (Gene Enrichment Analysis)
En este caso analizaremos el enriquecimiento únicamente para el conjunto de genes diferencialmente expresados, es decir, los 14 genes que hemos encontrado en la sección anterior.
## logFC AveExpr t P.Value adj.P.Val B ID
## 1 5.072934 8.474089 18.543219 5.327165e-09 1.417918e-05 7.792097 209774_x_at
## 2 4.211596 8.735476 18.283695 6.100471e-09 1.417918e-05 7.741081 204470_at
## 3 3.505208 8.032616 18.097937 6.729560e-09 1.417918e-05 7.703570 211506_s_at
## 4 1.957622 7.717278 13.252182 1.303480e-07 2.059824e-04 6.332773 210229_s_at
## 5 4.537659 7.163907 10.833333 8.472127e-07 1.071046e-03 5.222143 207850_at
## 6 1.543102 7.223723 9.820219 2.070937e-06 2.181733e-03 4.626485 218810_at
## PROBEID ENTREZID
## 1 209774_x_at 2920
## 2 204470_at 2919
## 3 211506_s_at 3576
## 4 210229_s_at 1437
## 5 207850_at 2921
## 6 218810_at 80149## [1] 14## ID Description
## GO:0071222 GO:0071222 cellular response to lipopolysaccharide
## GO:0071219 GO:0071219 cellular response to molecule of bacterial origin
## GO:0071216 GO:0071216 cellular response to biotic stimulus
## GO:0019221 GO:0019221 cytokine-mediated signaling pathway
## GO:0032496 GO:0032496 response to lipopolysaccharide
## GO:0002237 GO:0002237 response to molecule of bacterial origin
## GeneRatio BgRatio pvalue p.adjust qvalue
## GO:0071222 9/14 222/18903 6.873959e-15 5.517971e-12 2.524371e-12
## GO:0071219 9/14 234/18903 1.110256e-14 5.517971e-12 2.524371e-12
## GO:0071216 9/14 261/18903 2.995683e-14 9.925695e-12 4.540824e-12
## GO:0019221 10/14 496/18903 1.288915e-13 3.202954e-11 1.465293e-11
## GO:0032496 9/14 339/18903 3.194513e-13 6.350691e-11 2.905325e-11
## GO:0002237 9/14 360/18903 5.493612e-13 9.101083e-11 4.163579e-11
## geneID Count
## GO:0071222 CXCL2/CXCL1/CXCL8/CSF2/CXCL3/ZC3H12A/NFKBIA/IL1B/TNFAIP3 9
## GO:0071219 CXCL2/CXCL1/CXCL8/CSF2/CXCL3/ZC3H12A/NFKBIA/IL1B/TNFAIP3 9
## GO:0071216 CXCL2/CXCL1/CXCL8/CSF2/CXCL3/ZC3H12A/NFKBIA/IL1B/TNFAIP3 9
## GO:0019221 CXCL2/CXCL1/CXCL8/CSF2/CXCL3/NFKBIA/CCL20/IL1B/TNFAIP3/BIRC3 10
## GO:0032496 CXCL2/CXCL1/CXCL8/CSF2/CXCL3/ZC3H12A/NFKBIA/IL1B/TNFAIP3 9
## GO:0002237 CXCL2/CXCL1/CXCL8/CSF2/CXCL3/ZC3H12A/NFKBIA/IL1B/TNFAIP3 9
Vemos como algunos de los genes diferencialmente expresados son: NFKBIA, IL1B, TNFAIP3, CXCL3, CXCL1, CXCL2, CXCL8, e involucrados en procesos de señalización mediados por citoquinas, respuestas celulares a liposacáridos y lipopolisacáridos, respuestas celulares a estímulos bióticos BIRC3 y a moléculas de origen bacteriano.
Estos procesos biológicos enriquecidos están en línea con la activación del factor de transcripción NFkappaB que menciona el estudio principal en el que se basa este trabajo ([1]), ya que algunos inductores de la actividad de NFkappaB incluyen ([2-5]): factores de necrosis tumoral como TNFAIP3, interleucina 1-beta IL1B, lipoposacáridos bacteriales, entre otros. Por tanto, nuestro enrichment analysis proporciona más evidencia de estos hallazgos.
3.5 Análisis de enriquecimiento de conjuntos de genes (Gene Set Enrichment Analysis)
Ahora en lugar de mirar los genes individuales uno por uno, GSEA analiza los conjuntos de genes predefinidos, que serán vías biológicas conodcidas, para ver si aparecen de manera significativa en la lista de genes ordenada, la ventaja será que este enfoque puede detectar cambios sutiles en la expresión de un conjunto de genes relacionados que podrían pasar desapercibidos en un análisis de genes individuales.
## ID Description setSize
## hsa04668 hsa04668 TNF signaling pathway 76
## hsa04657 hsa04657 IL-17 signaling pathway 55
## hsa04064 hsa04064 NF-kappa B signaling pathway 54
## hsa05323 hsa05323 Rheumatoid arthritis 46
## hsa04621 hsa04621 NOD-like receptor signaling pathway 85
## hsa04080 hsa04080 Neuroactive ligand-receptor interaction 91
## enrichmentScore NES pvalue p.adjust qvalue rank
## hsa04668 0.8008032 2.771255 3.874625e-18 1.201134e-15 9.054387e-16 127
## hsa04657 0.8412803 2.814030 2.367961e-17 3.670339e-15 2.766775e-15 89
## hsa04064 0.8248821 2.749018 1.332480e-15 1.376896e-13 1.037932e-13 127
## hsa05323 0.8566989 2.788488 2.195424e-15 1.701453e-13 1.282590e-13 64
## hsa04621 0.7239302 2.528972 3.131200e-14 1.941344e-12 1.463424e-12 85
## hsa04080 -0.5224231 -2.908296 4.007800e-13 2.070697e-11 1.560933e-11 799
## leading_edge
## hsa04668 tags=32%, list=2%, signal=31%
## hsa04657 tags=33%, list=1%, signal=33%
## hsa04064 tags=28%, list=2%, signal=27%
## hsa05323 tags=28%, list=1%, signal=28%
## hsa04621 tags=18%, list=1%, signal=18%
## hsa04080 tags=48%, list=13%, signal=43%
## core_enrichment
## hsa04668 CXCL2/CXCL3/CXCL1/NFKBIA/TNFAIP3/BIRC3/CXCL6/CCL20/CSF2/IL1B/ICAM1/BCL3/IRF1/LIF/JUNB/PTGS2/JUN/CXCL5/IL6/MAP3K8/CEBPB/NFKB1/BIRC2/CYLD
## hsa04657 CXCL2/CXCL3/CXCL1/CXCL8/NFKBIA/TNFAIP3/CXCL6/CCL20/CSF2/IL1B/PTGS2/JUN/CXCL5/IL6/CEBPB/MMP1/NFKB1/MMP13
## hsa04064 CXCL2/CXCL3/CXCL1/CXCL8/NFKBIA/TNFAIP3/BIRC3/PLAU/IL1B/ICAM1/BCL2A1/PTGS2/NFKB1/BIRC2/CYLD
## hsa05323 CXCL2/CXCL3/CXCL1/CXCL8/CXCL6/CCL20/CSF2/IL1B/ICAM1/JUN/CXCL5/IL6/MMP1
## hsa04621 CXCL2/CXCL3/CXCL1/CXCL8/NFKBIA/TNFAIP3/BIRC3/IL1B/JUN/IL6/PANX1/NFKB1/ITPR1/BIRC2/PLCB4
## hsa04080 GNRH2/GRM4/POMC/S1PR2/NMB/GRIN1/MLNR/OXT/P2RX5/GABRE/LTB4R/HTR1A/CALCRL/ADM2/MTNR1B/NPFF/PAQR6/MC2R/TRPV1/HTR6/SPX/CHRNA3/OPRL1/GLP1R/ADRA2C/CHRNA1/HTR4/CHRM5/CHRNB3/GABRA3/HCRT/GPR35/FSHR/SSTR3/NPFFR1/LEPR/GABRA4/CCKBR/F2/SSTR4/UCN/PTAFR/GIPR/OPRM1Vemos como el proceso de señalización NFkappaB está enriquecido para nuestro conjunto de genes, es decir, hay una mayor proporción de conjuntos de genes relacionados con la señalización de NFkappaB del que cabría esperar. También hay un enriquecimiento de señalizaciones de citoquinas proinflamatorias como TNF o IL-17, se ha demostrado que estas citoquinas pueden activar la vía de señalización PI3K ([6-7]).
4 Discusión
5 Referencias
[1] Hughes-Fulford M, Li CF, Boonyaratanakornkit J, Sayyah S. (Febrero 2006). “Arachidonic acid activates phosphatidylinositol 3-kinase signaling and induces gene expression in prostate cancer”. Cancer Res 2006 Feb 1;66(3):1427-33. PMID: 16452198.
[2] Chandel NS, Trzyna WC, McClintock DS, Schumacker PT. (Julio de 2000). “Role of oxidants in NF-kappa B activation and TNF-alpha gene transcription induced by hypoxia and endotoxin”. J Immunol, 165(2): 1013-1021. PMID 10878378.
[3] Fitzgerald DC, Meade KG, McEvoy AN, Lillis L, Murphy EP, MacHugh DE, Baird AW. (Marzo de 2007). “Tumour necrosis factor-alpha (TNF-alpha) increases nuclear factor kappaB (NFkappaB) activity in and interleukin-8 (IL-8) release from bovine mammary epithelial cells”. Vet Immunol Immunopathol, 116(1-2): 59-68. PMID 17276517. doi:10.1016/j.vetimm.2006.12.008.
[4] Renard P, Zachary MD, Bougelet C, Mirault ME, Haegeman G, Remacle J, Raes M. (Enero de 1997). “Effects of antioxidant enzyme modulations on interleukin-1-induced nuclear factor kappa B activation”. Biochem Pharmacol, 53(2): 149-160. PMID 9037247.
[5] Qin H, Wilson CA, Lee SJ, Zhao X, Benveniste EN. (Noviembre de 2005). “LPS induces CD40 gene expression through the activation of NF-kappaB and STAT-1alpha in macrophages and microglia”. Blood, 106(9): 3114-3122. PMC 1895321. PMID 16020513. doi:10.1182/blood-2005-02-0759.
[6] Yanfei Bai, Haitao Li, Rui Lv. (Septiembre 2021). “Interleukin-17 activates JAK2/STAT3, PI3K/Akt and nuclear factor-\(\kappa\)B signaling pathway to promote tumorigenesis of cervical cancer”. PMC 8461522. pmid 34630646. doi:10.3892/etm.2021.10726.
[7] Pei Lei, Yibo Gan, Yuan Xu, Lei Song. (Febrero 2017). “The inflammatory cytokina TNF-\(\alpha\) promotes the premature senescence of rat nucleus pulposus cells via the PI3K/Akt signaling pathway”. Scientific reports 7(1):42938. doi:10.1038/srep42938.
6 Apéndices
A continuación se muestra el código en R para todos los análisis del informe.
########################
## CARGA DE DATOS ##
########################
# Carga de los datos
library("Biobase")
library("affy")
# Establecemos el directorio de trabajo
setwd("C:/Users/sadel/Desktop/PEC 2/GSE3737_RAW")
mywd <- getwd()
# Cargamos primero las anotaciones al objeto target
targets <- read.AnnotatedDataFrame(file.path(mywd, "Annotations.txt"), header = TRUE,
row.names = 1, sep = "\t")
# Listamos todos los ficheros .CEL
celFiles <- list.celfiles(mywd)
# Mediante read.affybatch() leemos todos los archivos .CEL y guardamos esta informacion
# junto con las anotaciones en el objeto rawData()
rawData <- read.affybatch(filenames = file.path(mywd, celFiles), phenoData = targets)
######################
## ANALISIS CALIDAD ##
######################
library("limma")
library("CCl4")
library("affyPLM")
# Nombre de los arrays
affySampleNames <- rownames(pData(rawData))
affyColores <- 1:8
# Signal distribution
hist(rawData, main = "Signal distribution", col = affyColores, lty = 1)
legend(x="topright", legend = affySampleNames, col = affyColores, lty = 1, cex = 0.7)
# Boxplot
boxplot(rawData, main = "Distribución de las expresiones", col = affyColores, las = 2)
# ACP
plotPCA <- function ( X, labels=NULL, colors=NULL, dataDesc="", scale=FALSE)
{
pcX<-prcomp(t(X), scale=scale) # o prcomp(t(X))
loads<- round(pcX$sdev^2/sum(pcX$sdev^2)*100,1)
xlab<-c(paste("PC1",loads[1],"%"))
ylab<-c(paste("PC2",loads[2],"%"))
if (is.null(colors)) colors=1
plot(pcX$x[,1:2],xlab=xlab,ylab=ylab, col=colors,
xlim=c(min(pcX$x[,1])-100000, max(pcX$x[,1])+100000),
ylim=c(min(pcX$x[,2])-100000, max(pcX$x[,2])+100000),
)
text(pcX$x[,1],pcX$x[,2], labels, pos=3, cex=0.8)
titulo <- ifelse(dataDesc=="", "Visualización de las dos primeras componentes", dataDesc)
title(titulo, cex=0.8)
}
plotPCA(exprs(rawData), labels = affySampleNames)
clust.euclid.average <- hclust(dist(t(exprs(rawData))), method = "average")
plot(clust.euclid.average, labels=colnames(exprs(rawData)), main = "Hierarchical clustering of samples")
######################
## PRE-PROCESADO ##
######################
library("genefilter")
# Cargamos la libreria con anotaciones del microarray Affymetrix Human Genome U133A
BiocManager::install("hgu133a.db")
library("hgu133a.db")
# Normalizar mediante Robust-Multiarray Average
normalizado <- rma(rawData)
# Filtrado de genes
filtered <- nsFilter(normalizado)
eset <- filtered[["eset"]]
# Vemos como el objeto retornado es del tipo ExpressionSet
class(eset)
# Control de calidad despues de normalizar y filtrar
par(mfrow=c(1,2))
# Nombre de los arrays
affySampleNames <- rownames(pData(eset))
affyColores <- 1:8
# Signal distribution
hist(eset, main = "Signal distribution", col = affyColores, lty = 1)
legend(x="topright", legend = affySampleNames, col = affyColores, lty = 1, cex = 0.7)
boxplot(eset, main = "Distribución de las expresiones", col = affyColores, las = 2)
# ACP
plotPCA <- function ( X, labels=NULL, colors=NULL, dataDesc="", scale=FALSE)
{
pcX<-prcomp(t(X), scale=scale) # o prcomp(t(X))
loads<- round(pcX$sdev^2/sum(pcX$sdev^2)*100,1)
xlab<-c(paste("PC1",loads[1],"%"))
ylab<-c(paste("PC2",loads[2],"%"))
if (is.null(colors)) colors=1
plot(pcX$x[,1:2],xlab=xlab,ylab=ylab, col=colors,
xlim=c(min(pcX$x[,1])-20, max(pcX$x[,1])+20),
ylim=c(min(pcX$x[,2])-20, max(pcX$x[,2])+20),
)
text(pcX$x[,1],pcX$x[,2], labels, pos=3, cex=0.8)
titulo <- ifelse(dataDesc=="", "Visualización de las dos primeras componentes", dataDesc)
title(titulo, cex=0.8)
}
par(mfrow=c(1,2))
plotPCA(exprs(eset), labels = affySampleNames)
clust.euclid.average <- hclust(dist(t(exprs(eset))), method = "average")
plot(clust.euclid.average, labels=colnames(exprs(eset)), main = "Hierarchical clustering of samples")
##############################
## EXPRESIÓN DIFERENCIAL ##
##############################
library("limma")
library("edgeR")
# Convertimos a factor la variable que indica el tratamiento
group <- as.factor(pData(eset)$characteristics..original.experiment.type)
# Relabel de los levels del factor
levels(group)[1] <- "Treated"
levels(group)[2] <- "Control"
# Matriz de diseno
design <- model.matrix(~0 + group)
colnames(design) <- levels(group)
# Estimamos el modelo lineal
fit <- lmFit(eset, design)
contrast.matrix <- makeContrasts(Treated-Control, levels=design)
contrast <- contrasts.fit(fit, contrast.matrix)
contrast_ebayes <- eBayes(contrast)
# Extraemos los resultados aplicando multiplicidad de Benjamini y Hochberrg
resultados <- topTable(contrast_ebayes, adjust = "fdr", number = nrow(contrast_ebayes))
# Anadir la ID como variable y no como rowname
resultados$ID <- rownames(resultados)
# Selecionar genes diferencialmente expresados (tanto Up como Down regulated)
selected <- resultados[,"adj.P.Val"]<0.05 & abs(resultados[,"logFC"]>1)
resultados[selected,]
# VOlcano plot
coefnum = 1
opt <- par(cex.lab = 0.7)
volcanoplot(contrast_ebayes, coef=coefnum, highlight = sum(selected), names = rownames(resultados),
main =paste("Differentially expressed genes", colnames(contrast.matrix)[coefnum], sep = "\n"), xlim = c(-5,5))
abline(v=c(-1,1))
##############################
## ANOTACIONES ##
##############################
library("mygene")
library("dplyr")
# Usamos para las anotaciones la base de datos, ya cargada, hgu133a.db
anotdata <- capture.output(hgu133a())
# Anotamos los genes con ENSEMBL
annotations_ese <- AnnotationDbi::select(hgu133a.db, keys = rownames(resultados), columns = c("ENTREZID"))
# unimos las anotaciones
topanotados <- resultados %>% mutate(PROBEID=rownames(resultados)) %>% left_join(annotations_ese) %>% arrange(P.Value)
head(topanotados)
# Selecionar genes diferencialmente expresados (tanto Up como Down regulated)
selected <- topanotados[,"adj.P.Val"]<0.05 & abs(topanotados[,"logFC"]>1)
sum(selected)
##############################
## Gene Enrichment Analysis##
##############################
library("GOstats")
# Over-representation analysis
topGenes <- topanotados[selected,"ENTREZID"]
# Eliminar posibles duplicados
topGenes <- topGenes[!duplicated(topGenes)]
# Ejecutar GO enrichment analysis
ego <- enrichGO(gene = as.integer(topGenes),
keyType = "ENTREZID",
OrgDb = org.Hs.eg.db,
ont = "BP",
pAdjustMethod = "fdr",
qvalueCutoff = 0.25,
readable = TRUE)
# Resultados de GO a un data frame
ego_results <- data.frame(ego)
head(ego)
# Visualizacion de los resultados GO
dotplot(ego, showCategory=10)
# En terminos de jerarquias
goplot(ego, showCategory = 10)
# Gene network
cnetplot(ego)
################################
# GENE SET ENRICHMENT ANALYSIS##
################################
library("clusterProfiler")
# ordenar por valor absoluto de logFC para eliminar duplicados con el menor
# valor absoluto de logFC
topanotados <- topanotados[order(abs(topanotados$logFC), decreasing=T),]
topanotados <- topanotados[!duplicated(topanotados$ENTREZID),]
# re ordenar basado en logFC para hacer el GSEA
topanotados <- topanotados[order(topanotados$logFC, decreasing=T),]
genesVector <- topanotados$logFC
names(genesVector) <- topanotados$ENTREZ
set.seed(123)
gseResulti <- gseKEGG(geneList = genesVector,
organism = "hsa",
keyType = "kegg",
exponent = 1,
minGSSize = 10,maxGSSize = 500,
pvalueCutoff = 0.05,pAdjustMethod = "fdr",
# nPerm = 10000, #augmentem permutacions a 10000
verbose = TRUE,
use_internal_data = FALSE,
seed = TRUE,
eps=0,
by = "fgsea"
)
library("kableExtra")
gsea.result <- setReadable(gseResulti, OrgDb = org.Hs.eg.db, keyType ="ENTREZID")
gsea.result.df <- as.data.frame(gsea.result)
print(kable(gsea.result.df[,c("Description","setSize","NES","p.adjust")])%>% scroll_box(height = "500px"))