-
Notifications
You must be signed in to change notification settings - Fork 1
/
SPIA - español.Rmd
149 lines (85 loc) · 5.21 KB
/
SPIA - español.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
---
title: "Pathway Enrichment Analysis"
subtitle: "(Análisis de enriquecimiento de vías)"
author: "Karina y Jesica Formoso"
output: html_document
---
Las tecnologías de genómica, transcriptómica y proteómica generan grandes cantidades de datos que luego deben ser analizados de la manera más eficaz posible. Una tendencia cada vez más marcada es la de analizar los genes obtenidos en grupos funcionalmente relacionados. Esto se logra mediante el **Pathway enrichment analysis**. El objetivo de este método es identificar grupos de genes con cambios de expresión posiblemente moderados pero coordinados en diferentes condiciones biológicas. Hay distintos tipos de análisis que se pueden realizar y se clasifican en:
- Competitivo vs. Autónomo (Competitive vs. Self-contained).
- Topológico vs. no topológico (Topological vs Non-Tophological).
Aquí vamos a utilizar SPIA que es un análisis mixto y topológica. Este método calcula un valor teniendo en cuenta el fold change de los genes, el score para el pathway enrichment y la topología de las vías de señalización.
Para trabajar con SPIA es necesario contar con la lista de genes diferencialmente expresados con sus log fold changes y la lista completa de genes en la plataforma.
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
```
Para este ejemplo vamos a utilizar la base de datos **colorectal cancer [Affymetrix geneChip technology (GEE GSE4107)]** y los paquetes SPIA y tidyverse:
#### Instalar SPIA
source("https://bioconductor.org/biocLite.R")
biocLite("SPIA")
#### Instalar tidyverse
install.packages(tidyverse)
```{r librerias}
library(tidyverse)
library(SPIA)
select <- dplyr::select
```
## Análisis topológico de enriquecimiento de vías con SPIA.
Vamos a cargar el set de datos **colorectalcancer**, el cual incluye un dataset llamado "top". Con la función head() podemos ver las primeras 6 observaciones o filas.
```{r data}
data(colorectalcancer)
head(top)
```
Usamos el paquete hgu133plus2.db que contiene Affymetrix Human Genome U133 Plus 2.0 Array annotation data para asignarle a cada ID de la base **top** un **ENTREZ ID**.
```{r nombresgenes}
library(hgu133plus2.db)
x <- hgu133plus2ENTREZID
top$ENTREZ <- unlist(as.list(x[top$ID]))
```
Vamos a seleccionar las observaciones que no tengan datos faltantes y eliminar los datos duplicados en ENTREZ. Como las observaciones están ordenadas según el log2FoldChange (logFC en la base), al quedarnos con la primera aparición de cada ENTREZ, retenemos la más significativa de cada probset.
```{r limpiarTOP}
top <- top %>%
filter(!is.na(ENTREZ),
!duplicated(ENTREZ))
```
Además, vamos a seleccionar las observaciones que tienen un p-valor ajustado menor a 0.1.
```{r padj}
tg1 <- filter(top, adj.P.Val < 0.1)
```
Luego, creamos un vector con los valores de logFC de la nueva base y utilizamos la variable ENTREZ para nombrar los valores del vector.
```{r col}
DE_Colorectal <- tg1$logFC
names(DE_Colorectal) <- tg1$ENTREZ
```
Creamos un segundo vector con los ENTREZ de la base TOP. Esta es la base original pero sin los valores faltantes o repetidos.
```{r tot}
ALL_Colorectal <- top$ENTREZ
```
Corremos el análisis utilizando la función SPIA(). Utilizamos el método de "fisher" para estudiar la significancia de la representación de nuestros genes en la vía de señalización.
```{r SPIA, results=FALSE}
resultados <- spia(de = DE_Colorectal, all = ALL_Colorectal, organism = "hsa", nB = 2000, plots = TRUE, verbose = TRUE, combine = "fisher")
```
Eliminemos la columna "KEGGLINK" y veamos las primeras líneas de nuestros resultados:
```{r resultados}
resultados %>%
select(-KEGGLINK) %>%
head()
```
El output obtenido es el siguiente:
### Columnas:
- pSize: el número de genes en la vía de señalización.
- NDE: el número de genes diferencialmente expresados en esta vía.
- tA: la acumulación total de perturbaciones observadas en la vía.
- pNDE: la probabilidad (p-valor) de observar al menos genes NDE en la vía utilizando un modelo hipergeométrico (similar a ORA)
- pPERT: la probabilidad (p-valor) de observar una acumulación total más extrema que tA solo por casualidad
- pG: el p-valor obtenido al combinar pNDE y pPERT
- pGFdr y pGFWER: los valores pG ajustados por la tasa de descubrimiento falso (FDR) y el método de Bonferroni, respectivamente.
- State: indica si la vía está inhibida o activada.
- KEGGLINK: proporciona un enlace web al sitio web de KEGG que muestra la imagen de la ruta con los genes expresados diferencialmente resaltados en rojo.
### Gráficos:
Como indicamos a la función SPIA que genere los gráficos asociados a estos resultados (plots = TRUE), estos se guardan como un archivo pdf en el directorio en el que estamos trabajando.
Podemos ver las vías significativamente desreguladas al ver la sobrerrepresentación y las perturbaciones de cada vía.
```{r plotP}
plotP(resultados, threshold = 0.05)
```
En este gráfico, cada vía es un punto y las coordenadas son el logaritmo de pNDE (usando un modelo hipergeométrico) y el p-valor de las perturbaciones (pPERT). Las líneas oblicuas muestran las regiones de importancia según la evidencia combinada.
link al código completo: ....