Download Diapositiva 1 - Departamento de Matematica

¿Cómo analizar los datos crudos de
microarrays?
Guía práctica de análisis de datos de microarrays
Juan Pablo Fededa y Carlos Rocco
Análisis Exploratorio y Confirmatorio de Datos de Experimentos de Microarrays
Primer Cuatrimestre 2006
Instituto de Cálculo y Departamento de Matemática.
Facultad de Ciencias Exactas y Naturales
Universidad de Buenos Aires
1
IMPORTANTE: Esta es una guía acerca de cómo analizar los datos
que se obtienen del scanneado de un microarray de dos colores. No
contiene información acerca de cómo construir, ensayar, diseñar un
experimento de y/ó scannear microarrays. El objetivo final es ayudar a
un biólogo a que, a partir del dato original de scanneado de un
microarray, obtenga una lista de genes diferencialmente expresados en
las dos condiciones confrontadas. Tampoco se hablará de data mining y
generación de patrones de comportamiento génico.
El análisis de datos que promueve esta guía involucra el uso de R, un
software libre de estadística sobre el cuál se fueron construyendo
múltiples herramientas para analizar microarrays.
Se asume que el lector tiene un background matemático básico y que
entiende los conceptos fundamentales de la tecnología de microarrays.
2
Microarrays……………
Fabricación de un microarray
FLUJO DE TRABAJO
u$s
Diseño del experimento
Preparación de la muestra - hibridación
Scanning
Análisis de datos - lista genes diferencialmente
expresados
Data mining – clustering
3
Microarrays……………
FLUJO DE TRABAJO
Fabricación de un microarray
Diseño del experimento
Preparación de la muestra - hibridación
Scanning
Análisis de datos - lista genes diferencialmente
expresados
Data mining – clustering (……..2007!!!!)
4
FLUJO DE TRABAJO
Análisis de Datos de Microarrays
Obtención de los datos crudos – Ingreso de
los datos al R – carga de paquetes de
análisis
Gráficos de estado del microarray – toma
de decisión ( sigo con el análisis o dejo? )
Normalizaciónes
Gráficos de normalización – toma de
decisión ( sigo con el análisis, realizo otras
normalizaciones, dejo? )
Identificación de genes diferencialmente
expresados
Gráficos y tablas de genes diferencialmente
expresados
5
FLUJO DE TRABAJO
Análisis de Datos de Microarrays
Obtención de los datos crudos – Ingreso de
los datos al R – carga de paquetes de
análisis
Gráficos de estado del microarray – toma
de decisión ( sigo con el análisis o dejo? )
Normalizaciónes
R
Gráficos de normalización – toma de
decisión ( sigo con el análisis, realizo otras
normalizaciones, dejo? )
Identificación de genes diferencialmente
expresados
Gráficos y tablas de genes diferencialmente
expresados
6
Obtención de R y paquetes de
análisis
Carga de paquetes de análisis
Obtención de los datos crudos
Ingreso de los datos al R
7
-1: Obtener el software R !!!! de http://cran.r-project.org/
8
-1: Obtener el software R !!!! de http://cran.r-project.org/
9
-1: Obtener el software R !!!! de http://cran.r-project.org/
10
0: Obtener ( bajar) los siguientes paquetes de trabajo de
R en http://www.bioconductor.org/packages/bioc/1.8/index.html :
Affyaffyio, arrayQuality, Biobase, colorspace, convert,
gridBase, hexbin, limma, marray, RColorBrewer, vsn
11
Una vez hecho esto instalar el R y dentro de la consola
del R instalar los paquetes bajados en archivos zip :
NOTA: Se pueden actualizar e instalar desde R (en otra de las opciones de12
package).
Una vez hecho esto instalar el R y dentro de la consola
del R instalar los paquetes bajados en archivos zip :
13
Una vez hecho esto instalar el R y dentro de la consola
del R instalar los paquetes bajados en archivos zip :
14
Antes que nada: vamos a trabajar escribiendo nuestras
ordenes en un script, para lo cuál primero hay que abrir
un script:
15
Es más fácil trabajar con las ventanas de la consola y del
script (Editor) en paralelo ( tile):
16
Entonces…..esto funciona escribiendo la orden en el script y
luego se utiliza el botón run line para correrla en la consola:
SCRIPT
CONSOLA
NOTA: run line se puede hacer con F5.
17
La orden se corre en la consola:
18
Cargamos el paquete marray para empezar a trabajar (al
cargar marray se carga automaticamente limma )
19
Obtención de R y paquetes de
análisis
Carga de paquetes de análisis
Obtención de los datos crudos
Ingreso de los datos al R
20
OK; los microarrays sobre los cuáles vamos a trabajar como
ejemplo corresponden al siguiente paper:
21
A los datos originales los bajamos de la base de datos de
microarrays de stanford: http://genome-www5.stanford.edu/
22
Vamos a la carpeta SMD, en donde están los datos crudos:
23
Particularmente vamos a analizar los datos de los
microarrays correspondientes a los experimentos de
inmunoprecipitación de PUMILIO en moscas adultas:
24
Vamos a utilizar los datos originales del scanner (ORI DATA),
en este caso son 4 experimentos de hibridación: exptID 52987,
53530, 54253 y 54596
25
Utilizamos (bajamos) los archivos .gpr (genepix):
26
NOTA: Diferentes scanners generan
diferentes tipos de archivos con
diferentes extensiones; los paquetes de
análisis de datos que utilizaremos leen la
mayoría de estas extensiones.
Agilent
Arrayvision
Bluefuse
Genepix
spot.close.open
genepix.custom
genepix.median
Imagene
Smd
Quantarray
Scanarrayexpress
smd.old
Spot
El procedimiento para ingresar estas files en R es el mismo
que el que veremos a continuación.
27
Obtención de R y paquetes de
análisis
Carga de paquetes de análisis
Obtención de los datos crudos
Ingreso de los datos al R
28
Volviendo a R, cambiamos de directorio en el cuál están
nuestros archivos originales salidos del scanner (ej.: .gpr
para genepix):
29
Cambiamos de directorio en el cuál están nuestros archivos
originales salidos del scanner (ej.: .gpr para genepix):
30
Cambiamos de directorio en el cuál están nuestros archivos
originales salidos del scanner (ej.: .gpr para genepix):
31
Empezamos a trabajar en R:
dirección.experimentospumilio <-"D://juan//análisis exploratorio y confirmatorio de
experimentos de microarrays//curso microarrays//pumilio//Experiment Sets Name
Organism Type Created By Options“
Esta es la otra opción para indicar la carpeta en donde están nuestros datos (Por
defecto lee en la carpeta de trabajo).
setwd(dirección.experimentospumilio)
R va a leer los datos desde esta dirección seteada.
archivosexperimentospumilio <- c("52987.gpr", "53530.gpr", "54253.gpr", "54596.gpr")
Generamos un objeto marrayRaw con los archivos de los arrays que vamos a
analizar.
crudosexperimentospumilio<-read.GenePix(archivosexperimentospumilio)
Generamos un objeto que NO ignore las estimaciones de las intensidades del
background para los canales verde ( G ) y rojo ( R ).
crudosexperimentospumilio2<-read.GenePix(archivosexperimentospumilio, name.Gb=
NULL, name.Rb= NULL)
Generamos un objeto que ignore las estimaciones de las intensidades del background
para los canales verde ( G ) y rojo ( R ) fijando el valor NULL para los argumentos
32
name.Gb y name.Rb.
33
Estructura de los datos que estamos analizando en el
paquete marray:
etceteras
Gb (background Cy3)
Rb (background Cy5)
Gf (foreground Cy3)
microarreglos
genes
Rf (foreground Cy5)
34
summary(crudosexperimentospumilio)
Nos muestra un resumen de las características del objeto
crudosexperimentospumilio que generamos:
Son datos sin normalizar, es
clase marray, contiene los
datos de los 4 microarrays
17664 spots organizados en 8
filas x 4 columnas de
subconjuntos ( grids ), cada grid
contiene 23 filas x 24 columnas
de spots
35
slotNames (crudosexperimentospumilio)
Nos muestra la estructura de slots del objeto crudosexperimentopumilio.
slotNames (crudosexperimentospumilio2)
Nos muestra la estructura de slots del objeto crudosexperimentopumilio2.
Ambos objetos (crudosexperimentospumilio, crudosexperimentospumilio2 )
tienen la misma estructura de slots. Los primeros 5 ("maRf" "maGf"
"maRb" "maGb" "maW“) son las matrices que contienen la información
cuantitativa básica extraída de los archivos .gpr. Los demás están
asociados con los archivos de la estructura de los arreglos (layout) y las
36
anotaciones.
crudosexperimentospumilio@maRf
Objeto con las intensidades del foreground para el canal rojo.
crudosexperimentospumilio@maRb
Objeto con las intensidades del background para el canal rojo.
crudosexperimentospumilio@maGf
Objeto con las intensidades del foreground para el canal verde.
crudosexperimentospumilio@maGb
Objeto con las intensidades del background para el canal verde.
crudosexperimentospumilio@maLayout
Objeto con la estructura o geometría del arreglo.
crudosexperimentospumilio@maGnames
Objeto con los nombres de los genes.
crudosexperimentospumilio@maTargets
Objeto con información acerca de que muestras fueron hibridadas en cada
37
arreglo.
crudosexperimentospumilio@maGf[1:5,]
Generamos un objeto con las intensidades del foreground para el canal
verde para las 5 primeras filas (genes) para todos los microarrays
En este caso el microarray de
la columna 1 es 52987
[1:5, 1 ]
microarrays
genes
Recuerden la estructura
de datos de marray!!!
38
Existen métodos específicos para objetos de la clase marrayRaw
(como crudosexperimentospumilio), que luego se utilizaran para tomar
decisiones acerca de los datos analizados:
maA : matriz de log2 de intensidades (con corrección por
background)
maM : matriz de log2 de cocientes (con corrección por background)
maLR : matriz de log2 de ( intensidades - background) para el
canal rojo
maLG: matriz de log2 de ( intensidades - background) para el canal
verde
LR = log2 ( Rf -Rb )
LG =log2 ( Gf - Gb )
A = 0.5 x ( LR + LG )
M= LR - LG
39
Los datos se transforman aplicando log2
para que las intensidades se distribuyan
en forma aproximadamente simétrica.
Esto mejora la visualización gráfica de
los datos. También se intenta que se
reduzca la relación entre la intensidad y
la varianza que aparece cuando se
realizan experimentos con
replicaciones. La mayoría de los
paquetes de análisis de microarrays
utilizan log2 .
40
Gráficos espaciales del estado del
microarray
41
image(crudosexperimentospumilio[,1], xvar="maM“)
Generamos un gráfico de los valores de M del primer microarray ( [ ,1 ] ) de
crudosexperimentospumilio:
valores altos de M
(cy5 > cy3)
valores bajos de M
(cy5 < cy3)
42
image(crudosexperimentospumilio[,1], xvar="maRb“)
Generamos un gráfico de los valores del background para el canal rojo del
primer microarray ( [ ,1 ] ) de crudosexperimentospumilio:
valores altos de
background en el
canal rojo
valores bajos de
background en el
canal rojo
43
par(mfrow=c(2,2))
Generamos un espacio gráfico de 2 x 2 gráficos.
image (crudosexperimentospumilio[,1], xvar="maSpotCol", bar=FALSE)
Generamos un gráfico del vector columna para cada spot en el primer microarray (
[ ,1 ] ) de crudosexperimentospumilio ( bar=FALSE : para que no dibuje la barra de
colores).
image (crudosexperimentospumilio[,1], xvar="maPrintTip", bar=FALSE)
Generamos un gráfico del vector print-tip para cada spot en el primer microarray (
[ ,1 ] ) de crudosexperimentospumilio
image (crudosexperimentospumilio[,1], xvar="maControls", col=heat.colors(10),
bar=FALSE)
Generamos un gráfico espacial que indica en donde se encuentran los spots
control en el primer array de ( [ ,1 ] ) de crudosexperimentospumilio, indicados en
rojo
image (crudosexperimentospumilio[,1], xvar="maPlate", bar=FALSE)
Generamos un gráfico de los valores de M del primer microarray ( [ ,1 ] ) de
crudosexperimentospumilio
44
45
boxplot(crudosexperimentospumilio[,1], xvar = "maPrintTip", yvar = "maM",main="arreglo
52987.gpr")
Podemos generar un gráfico boxplot que nos muestre cómo varía M en función del printtip para el primer microarray ( [ ,1 ] ) de crudosexperimentospumilio:
46
boxplot(crudosexperimentospumilio,yvar="maM")
Podemos generar un gráfico boxplot que nos muestre cómo varía M en función del array
para crudosexperimentospumilio:
47
library (arrayQuality)
Cargamos el paquete arrayQuality .
maQualityPlots(crudosexperimentospumilio [,1])
Generamos un gráfico maQualityPlots para evaluar a priori la calidad delmicroarray:
48
1. MA-plot de los datos crudos sin sustracción de background. Cada línea
coloreada representa la curva loess para cada grupo de print-tip. Estas
curvas representan que cantidad de normalización se deberá realizar a cada
grupo de datos. Mientras más separadas las curvas de la línea 0 de M, mayor
normalización será necesaria. Mientras mayor la normalización necesaria,
menor la calidad del experimento.
49
2. MA-plot usando el paquete hexbin que destaca la densidad de puntos,
de los datos normalizados; amarillo = alta densidad de puntos, azul = baja
densidad. Este diagrama nos indica como quedó la normalización por
grupo de print-tip (default).
50
3.Gráfico espacial de los valores de M
crudos sin sustracción del fondo.
Amarillo = valores bajos de M, azul =
valores altos de M, blanco = valores
faltantes . Esto nos permite visualizar
una hibridación despareja.
51
4. Gráfico espacial de los rangos de
los valores de M normalizados. Por
defecto, se utiliza la normalización
loess por print-tip-group. Además,
los puntos marcados (flagged) son
señalados por medio de un
recuadro negro. Este tipo de
representación gráfica permite
verificar que la normalización ha
quitado efectos espaciales.
52
5. Gráfico espacial de los valores de A
crudos sin sustracción del fondo. El
color representa la intensidad de la
señal.
53
6. Histogramas del cociente señal ruido
(signal-to-noise ratio) (S2N) para cada canal
Cy5 y Cy3.
S2N = log2(intensidad del spot / intensidad del
fondo)
En la parte superior de cada histograma se
imprimen la media y varianza de la señal.
Además, se superponen curvas de densidad
de SNR estratificado por diversos tipos de
control (estado). Los esquemas de color de las
curvas de densidad están dados la tabla 1. El
SNR es un buen indicador para los problemas
del tinte. Las líneas de densidad de controles
negativos (por ej. secuencias artificiales que
han sido diseñadas para no presentar
homología con ningún genoma) y vacíos
deben estar cercanas, casi superpuestas.
54
7. Diagrama de punto de los valores M
normalizados para los controles. Los
controles con más de 3 réplicas se
representan en el eje Y, el esquema de
colores se representa en la tabla 1. Los
valores de los controles M deben estar
concentrados y cerca de 0.
55
8. Diagrama del punto de A, sin la
sustracción del fondo. Los controles
con más de 3 réplicas se representan
en el eje Y, el esquema de colores se
representa en la tabla 1. La intensidad
de los controles positivos debe estar en
la región de alta intensidad, los
controles negativos y vacíos deben
estar en la región de intensidad más
baja. Los rangos de los controles
positivos y los negativos/vacíos deben
estar separados.
56
Normalizaciónes
Gráficos de normalización
57
Normalizaciones de datos : ¿ POR QUE NORMALIZAMOS?
Para analizar los datos es necesario reducir el error sistemático. Como
suponemos que sólo una pequeña cantidad de genes se expresa
diferencialmente, las distribuciones de intensidades deben ser
uniformes espacialmente dentro del arreglo y entre arreglos.
Entonces… normalizaremos primero dentro de un microarray utilizando
la función maNorm y luego, si fuese necesario, entre arreglos con la
función maNormScale.
¿Cómo nos damos cuenta si es necesario?
Luego de la normalización dentro de cada microarray, vamos a realizar
diferentes gráficos que nos van a indicar si la escala de los datos entre
microarrays es diferente o no. Si es muy diferente, tendremos que
normalizar por escala.
Recuerden que por cada paso de normalización estamos perdiendo
información.
58
Normalización dentro de cada microarray:
experimentospumilionorm <- maNorm(crudosexperimentospumilio, norm= "p")
Generamos un objeto que contiene los datos normalizados por grupo de printtip ( norm = “p” : normaliza por intensidad DENTRO de cada aguja ) de
crudosexperimentospumilio.
summary (experimentospumilionorm)
Nos muestra un resumen de las características del objeto
experimentospumilionorm que generamos.
summary (crudosexperimentospumilio)
Nos muestra nuevamente un resumen de las características del objeto
crudosexperimentospumilio.
59
Normalización dentro de cada microarray:
Podemos comparar la parte estadística de crudosexperimentospumilio (sin
normalizar) con la misma de experimentospumilionorm ( mismos datos
normalizados por print-tip group):
antes
después
Las medias y las medianas
quedaron cerca de cero en todos
los arrays del experimento!!!!
60
Normalización dentro de cada microarray: veámoslo en gráficos
par(mfrow=c(2,1))
Generamos un espacio gráfico de 2 gráficos independientes.
boxplot(crudosexperimentospumilio[,1], xvar = "maPrintTip", yvar =
"maM",main="crudosexperimentospumilio")
Generamos un gráfico boxplot MA ( yvar = "maM" ) de los datos sin normalizar
de crudosexperimentospumilio, divididos por grupo de print-tip ( xvar =
"maPrintTip“ ) del primer microarray ( [ ,1 ] ).
boxplot(experimentospumilionorm[,1], xvar = "maPrintTip", yvar =
"maM",main="experimentospumilionorm")
Generamos un gráfico boxplot MA ( yvar = "maM" ) de los datos normalizados
por print-tip de experimentospumilionorm, divididos por grupo de print-tip ( xvar
= "maPrintTip“ ) del primer microarray ( [ ,1 ] ).
61
Normalización dentro de cada microarray: veámoslo en gráficos
Antes de normalizar
por Print-tip
Hay variaciones entre las
medias o medianas, en los
spots spoteados por
diferentes print-tips.
La mayoría de los valores
está debajo de cero
Luego de normalizar
por Print-tip
Transformamos los
datos de manera que
las medias o medianas
de cada print tip son
iguales ( 0) !!!!
62
Normalización dentro de cada microarray: veámoslo en gráficos
par(mfrow=c(2,1))
Generamos un espacio gráfico de 2 gráficos independientes.
boxplot(crudosexperimentospumilio, yvar = "maM", main = "pre
normalización")
Generamos un gráfico boxplot MA ( yvar = "maM" ) de los datos sin
normalizar de crudosexperimentospumilio, dividos por microarray.
boxplot(experimentospumilionorm, yvar = "maM", main = "post
normalización")
Generamos un gráfico boxplot MA ( yvar = "maM" ) de los datos
normalizados de experimentospumilionorm, divididos por microarray.
63
Normalización dentro de cada microarray: veámoslo en gráficos
Antes de normalizar
Hay variaciones entre las
medias o medianas en los
datos de los diferentes
microarrays.
La mayoría de los valores
está debajo de cero
Luego de normalizar
Transformamos los datos
de manera que las medias
o medianas de todos los
microarrays son similares
( 0) !!!!
64
Normalización dentro de cada microarray:
veámoslo en gráficos ( MA plot )
par(mfrow=c(2,1))
Generamos un espacio de 2 gráficos.
maPlot(crudosexperimentospumilio, legend.func = "NULL")Generamos un
gráfico Generamos un gráfico MA plot de los datos sin normalizar de
crudosexperimentospumilio; las líneas de colores representan las curvas
loess para cada print tip group.
maPlot(experimentospumilionorm, legend.func = "NULL")
Generamos un gráfico MA plot de los datos normalizados de
experimentospumilionorm; las líneas de colores representan las curvas
loess para cada print tip group.
65
Normalización dentro de cada microarray: veámoslo en gráficos ( MA plot )
Antes de
normalizar
Luego de
normalizar
66
Normalización entre microarrays ( por escala ) :
experimentospumilionormescala <- maNormScale(experimentospumilionorm)
Generamos un objeto que contiene los datos normalizados por escala (
función maNormScale ) a partir de objeto de datos prenormalizados por printtip group experimentospumilionorm.
par(mfrow=c(2,1))
Generamos un espacio gráfico de 2 gráficos.
boxplot(experimentospumilionorm, yvar = "maM", main = "pre
normalización")
Generamos un gráfico boxplot de M ( yvar = "maM" ) de los datos sin
normalizar por escala de experimentospumilionorm, dividos por microarray.
boxplot(experimentospumilionormescala, yvar = "maM", main = "post
normalización")
Generamos un gráfico boxplot de M ( yvar = "maM" ) de los datos
normalizados por escala de experimentospumilionormescala, divididos por
67
microarray.
Normalización entre microarrays ( por escala ): veámoslo en gráficos
Antes de normalizar
por escala
Hay variaciones entre
las distribuciones de los
datos de los diferentes
microarrays
Luego de normalizar
por escala
Luego de normalizar
por escala, las cajas
de los box-plots
tienen anchos casi
idénticos!!!!
68
Normalización entre microarrays ( por cuantiles ) :
library(convert)
Necesitamos cargar el paquete convert para transformar los objetos que obtuvimos
en el paquete marray ( crudosexperimentospumilio, experimentospumilionorm ) a
objetos legibles en el paquete limma, sobre el cuál realizaremos las normalizaciones
por cuantiles.
crudosexperimentospumilio.limma <- as(crudosexperimentospumilio,"RGList")
Creamos un objeto crudosexperimentospumilio.limma que contiene los datos del
objeto crudosexperimentospumilio de marray.
experimentospumilionorm.limma <- as(experimentospumilionorm,"MAList")
Creamos un objeto experimentospumilionorm.limma que contiene los datos del objeto
experimentospumilionorm de marray.
experimentospumilionormescala.limma <as(experimentospumilionormescala,"MAList")
Creamos un objeto experimentospumilionormescala.limma que contiene los datos del
objeto experimentospumilionormescala de marray.
experimentospumilionormporquantiles.limma <normalizeBetweenArrays(experimentospumilionorm.limma,method="quantile")
Creamos un objeto experimentospumilionormporquantiles.limma que contiene los
datos de la normalización por cuantiles de los datos ya normalizados por print-tip
group.
69
Normalización entre microarrays ( por cuantiles ): veámoslo en gráficos
par(mfrow=c(2,1))
Generamos un espacio gráfico de 2 gráficos independientes.
boxplot(experimentospumilionormescala.limma$M~col(experimentospumilionorm
escala.limma$M), names=colnames(experimentospumilionormescala.limma$M),
main="norm. por escala", col="red")
Generamos un gráfico boxplot de M de los datos normalizados por escala
experimentospumilionormescala.limma, dividos por microarray.
boxplot(experimentospumilionormporquantiles.limma$M~col(experimentospumili
onormporquantiles.limma$M),names=colnames(experimentospumilionormporqua
ntiles.limma$M), main="norm. por quantiles", col="red")
Generamos un gráfico boxplot de M de los datos normalizados por escala
experimentospumilionormporquantiles.limma, dividos por microarray.
70
Normalización entre microarrays ( por cuantiles ): veámoslo en gráficos
Antes de normalizar
por cuantiles
Hay variaciones entre las
distribuciones
de
los
datos de intensidad de
los diferentes microarrays
Luego de normalizar
por cuantiles
Las diferencias entre
las distribuciones se
reducen!!!!
71
Podemos también normalizar enteramente en el paquete limma:
Para ello generaremos un archivo .txt con información acerca de los
experimentos que estamos analizando, que tenga la siguiente
estructura de datos:
SlideNumber
075-dm0603B2
079-dm0603B2
085-dm0603B2
097-dm0603B2
FileName
52987.gpr
53530.gpr
54253.gpr
54596.gpr
Cy3
RNA total
RNA total
RNA total
RNA total
Cy5
RNA IP
RNA IP
RNA IP
RNA IP
Date
2004/06/11
2004/07/06
2004/07/27
2004/08/09
Llamaremos a este archivo PumilioSample.txt , más tarde nos servirá
para poder trabajar con los datos.
72
Normalización entre microarrays ( en limma ):
library(limma)
Cargamos el paquete limma para trabajar sobre el mismo.
targetspumilio<-readTargets("PumilioSample.txt")
Generamos un objeto targetspumilio que es un dataframe, generado por
la función readTargets, con los datos de PumilioSample.txt.
f <- function(x) as.numeric(x$Flags > -99)
Generamos una función f que se utilizará para fijar un filtro de manera
que cualquier spot marcado con ( flag de) -99 o menos tenga peso cero.
RGpumilio <- read.maimages(targetspumilio$FileName,
source="genepix", wt.fun=f)
Los nombres de los archivos guardados en targetspumilio permiten
identificar la información de intensidades contenida en los archivos de
salida del procesamiento de imágenes, en este caso son archivos con
extensión .gpr
73
experimentospumiliocorregidosporRMA <- backgroundCorrect(RGpumilio,
method="rma")
Generamos un objeto experimentospumiliocorregidosporRMA que contiene
los datos de la corrección por “rma”, que es la recomendada por la guía de
usuarios de limma para los datos de GenePix.
experimentospumilionormalizadosporprinttip <normalizeWithinArrays(experimentospumiliocorregidosporRMA) # por defecto
es print tip loess
Generamos un objeto experimentospumilionormalizadosporprinttip que
contiene los datos de la normalización por grupo de print-tip de los datos
previamente corregidos por rma (experimentospumiliocorregidosporRMA ).
experimentospumilionormalizadosporprinttipycuantiles<normalizeBetweenArrays(experimentospumilionormalizadosporprinttip,RG$pri
nter,method="quantile")
Generamos un objeto experimentospumilionormalizadosporprinttipycuantiles
que contiene los datos de la normalización por cuantiles de los datos
experimentospumilionormalizadosporprinttip.
experimentospumilionormalizadosporvsn<normalizeBetweenArrays(experimentospumiliocorregidosporRMA,RG$printer,
method="vsn")
Generamos un objeto experimentospumilionormalizadosporprinttipyvsn que
contiene los datos de la normalización por cuantiles de los datos
74
experimentospumilionormalizadosporprinttip.
Normalización entre microarrays ( en limma ): veámoslo en gráficos
par(mfrow=c(2,2))
Generamos espacio de 2 x 2 gráficos.
plotDensities(experimentospumiliocorregidosporRMA)
text(15,0.8,"RMA")
Generamos un gráfico de las densidades en función de la intensidad de
experimentospumiliocorregidosporRMA.
plotDensities(experimentospumilionormalizadosporprinttip)
text(15,0.45,"print tip")
Generamos un gráfico de las densidades en función de la intensidad de
experimentospumilionormalizadosporprinttip.
plotDensities(experimentospumilionormalizadosporprinttipycuantiles)
text(15,0.25,"cuantiles")
Generamos un gráfico de las densidades en función de la intensidad de
experimentospumilionormalizadosporprinttipycuantiles.
plotDensities(experimentospumilionormalizadosporvsn)
text(15,0.12,"vsn")
Generamos un gráfico de las densidades en función de la intensidad de
experimentospumilionormalizadosporprinttipyvsn.
75
RG densities
0.6
5
6
0.4
0.2
10
15
Intensity
RG densities
RG densities
cuantiles
4
5
Intensity
0.00 0.10 0.20 0.30
Print-tip + cuantiles
15
8
10
Intensity
14
0.00 0.05 0.10 0.15
Nos quedamos con esta
normalización:
10
Density
0
Density
print tip
0.0
Density
0.8
0.4
RMA
0.0
Density
1.2
RG densities
vsn
0
5
10
Intensity
15
20
76
Identificación de genes
diferencialmente expresados
Gráficos y tablas de genes
diferencialmente expresados
77
OK, ahora que ya tenemos los datos normalizados,
necesitamos saber en que canales ( Cy3 – Cy5 ) está hibridada
cuál de las dos muestras ( RNA total – RNA inmunoprecipitado
con pumilio ) en cada uno de los microarrays que analizamos,
clickeamos el botón view de la página de Stanford de la cuál
bajamos nuestros datos:
78
RNA total Cy3 = canal 1
RNA IP Cy5 = canal 2
79
Exp ID
Cy3
Cy5
52987
RNA total
RNA IP
53530
RNA total
RNA IP
54253
RNA total
RNA IP
54596
RNA total
RNA IP
Aquí observamos un gran error que hicieron nuestros amigos
de Stanford en cuanto al diseño experimental: NO hicieron
dye swapping en este experimento,
RECUERDEN REALIZAR DYE SWAPPING EN SUS
EXPERIMENTOS!!!!!
80
Podemos adoptar el siguiente modelo lineal para nuestro
experimento:


RNA _ total 
Cy5 
  E  log 2
  E M   
E  log 2
RNA _ IP 
Cy3 


o, equivalentemente:

Cy3 
  E  M    
E  log 2
Cy5 

La matriz de diseño
para un experimento
con dye-swap sería:
para un dye-swap:
Cy3
Cy5
matriz de
diseño
RNA IP
RNA total
1
RNA total
RNA IP
-1
81
Si asignamos β = E(log2 RNA totalCy5/RNA IPCy3) nos queda
una grilla de este tipo:
-M
M
Exp ID
Cy3/Cy5
Cy5/Cy3
52987
-1
1
53530
-1
1
54253
-1
1
54596
-1
1
Con estos datos creamos una matriz de diseño del experimento que será
la siguiente : design = c(-1,-1,-1,-1) para estimar la relación
RNA_total/RNA_IP (el inverso de M, ya que aquí marcamos con Cy5 a
RNA IP).
82
pumiliofit1 <-lmFit(experimentospumilionormalizadosporprinttipycuantiles,
design=c(-1,-1,-1,-1))
Generamos un objeto pumiliofit1 que contiene los datos de M ajustados
de cada gen (previamente normalizados) con la función lmFit. Ésta, por
defecto realiza un ajuste por cuadrados mínimos ordinario. Tiene opciones
para incluir una estructura de covarianzas de los errores y un ajuste
robusto. La matriz design indica cuáles arreglos son dye-swaps.
pumiliofit <- eBayes(pumiliofit1)
La función eBayes calcula el estadístico t moderado (Lönnstedt and Speed
2001) y un p-valor, para cada gen.
summary(pumiliofit)
Nos muestra las características de pumiliofit.
83
Tablas de genes diferencialmente expresados:
tablatop10Bpumilio <topTable(pumiliofit,coef=1,number=10,genelist=pumiliofit$genes,
adjust.method="BH",sort.by="B",resort.by=NULL)
Generamos un objeto que contiene una tabla con los diez primeros genes
con score mayor para el valor estadístico B
tablatop10Bpumilio
Visualizamos tablatop10Bpumilio.
tablatop10tpumilio <topTable(pumiliofit,coef=1,number=10,genelist=pumiliofit$genes,
adjust.method="BH",sort.by="t",resort.by=NULL)
Generamos un objeto que contiene una tabla con los diez primeros genes
con score mayor para el valor estadístico t (t de student).
tablatop10Mpumilio <topTable(pumiliofit,coef=1,number=10,genelist=pumiliofit$genes,
adjust.method="BH",sort.by="M",resort.by=NULL)
Generamos un objeto que contiene una tabla con los diez primeros genes
84
con score mayor para el valor de M.
85
tablatop17000Bpumilio <topTable(pumiliofit,coef=1,number=17000,genelist=pumiliofit$genes,
adjust.method="BH",sort.by="B",resort.by=NULL)
Generamos un objeto que contiene una tabla con los 17000 primeros
genes con score mayor para el valor estadístico B (Bonferroni).
pumiliomenoresa00000001<subset(tablatop17000Bpumilio,P.Value<0.0000001,c(ID,
Name,M,A,t,P.Value,B))
Generamos un objeto que contiene una tabla con todos los genes
con un estadístico p menor a 0.0000001.
dim(pumiliomenoresap00000001)
Nos dice las dimensiones de ese objeto ( cuantos genes contiene ).
Pumiliomenoresap00000001
Visualizamos pumiliomenoresap00000001.
86
87
Puedo ver en un gráfico MA cuáles son los genes top score de expresión
diferencial, con respecto a diferentes criterios de selección ( M, t o B):
M<-pumiliofit$coefficients
Generamos un objeto M que contenga los datos del espacio coefficients ( valores
M ) de pumiliofit.
A<-pumiliofit$Amean
Generamos un objeto A que contenga los datos del espacio Amean ( valores A ) de
pumiliofit.
plot(A,M,pch=".",col="lightblue",cex=3,main="Valores de M ")
Generamos un gráfico de M vs. A .
points(tablatop10Bpumilio$A[1],tablatop10Bpumilio$M[1],pch=“?",col="blue",cex=2)
;text(12.27,4.738,"CG2229")
Generamos un punto en el gráfico anterior, que representa al gen con mejor score
de tablatop10Bpumilio, ubicándolo espacialmente por su A ( 12.27 ) y M ( 4.738 ).
points(tablatop10Mpumilio$A[1],tablatop10Mpumilio$M[1],pch=“?",col="blue",cex=2
);text(9.206855,-6.772648,"CG7979")
Generamos un punto en el gráfico anterior, que representa al gen con mejor score
de tablatop10Mpumilio, ubicándolo espacialmente por su A (9.206855 ) y M (- 88
6.772648 ).
89
También puedo buscar que datos estadísticos de expresión diferencial
presenta un gen x del cuál conozco su nombre o ID:
genquemeinteresa<subset(tablatop17000Bpumilio,ID=="CT4618.a",c(ID,Name,M,A,t,P.Value,B))
Genero un objeto genquemeinteresa que contenga los datos y estadísticos
siguientes: ID,Name,M,A,t,P.Value,B del gen con ID CT4618.a .
genquemeinteresa
Visualizo los datos del objeto genquemeinteresa.
90
También podemos intersectar los 100 genes con mayores valores de M,
B y t, los cuáles van a ser los mayores candidatos a estar
diferencialmente expresados:
ordinary.t <- pumiliofit$coef / pumiliofit$stdev.unscaled / pumiliofit$sigma
Generamos un objeto que contiene los datos de coef, stdev.unscaled y sigma de
pumiliofit.
ordM <- order(abs(pumiliofit$coefficients),decreasing=TRUE)
Generamos un objeto que contiene los datos ordenados de M (coefficients) en orden
decreciente (valores absolutos).
ordB <- order(pumiliofit$lods,decreasing=TRUE)
Generamos un objeto que contiene los datos ordenados de B (lods) en orden
decreciente.
ordt <- order(abs(ordinary.t),decreasing=TRUE)
Generamos un objeto que contiene los datos ordenados de t (ordinary.t) en orden
decreciente (valores absolutos).
vector <- rep(0,17664)
Generamos un vector de tantos valores como spots tengael array ( 0 a 17664 en éste
caso).
top100M <- ordM[1:100]
Generamos un objeto que contenga los 100 primeros valores de ordM.
top100B <- ordB[1:100]
Generamos un objeto que contenga los 100 primeros valores de ordB.
91
top100t <- ordt[1:100]
Generamos un objeto que contenga los 100 primeros valores de ordt.
M <- vector
Asigno a M el objeto vector.
B <- vector
Asigno a B el objeto vector.
t <- vector
Asigno a t el objeto vector.
B[top100B]<-1
A los elementos de top100B que estén en el vector B, les asigno un valor
de 1.
t[top100t]<-1
A los elementos de top100t que estén en el vector t, les asigno un valor de
1.
M[top100M]<-1
A los elementos de top100M que estén en el vector M, les asigno un valor
de 1.
interseccion <- (B*t*M)/(B*t*M)
Genero un objeto interseccion que contiene los datos resultantes de dicha
operación entre vectores.
interseccion
Visualizamos intersección.
92
gen 10020
gen 10328
93
¿Cómo se qué gen es el 10020?
Names<-pumiliofit$genes$Name
Generamos un objeto Names que contiene los datos de nombres de genes de
pumiliofit.
diferencial10020<-Names[10020]
Generamos un objeto que contenga el nombre del gen 10020.
diferencial10020
Visualizamos a diferencial10020.
94
¿Cuáles son los genes de la intersección entre t, B y M?:
pumiliointerseccion2 <- data.frame (interseccion, Names)
Generamos un objeto pumiliointerseccion2 que contiene dos columnas de datos,
en la primera interseccion y en la segunda names.
subset (pumiliointerseccion2 , interseccion = =1)
Visualizamos los elementos de pumiliointerseccion2 que tienen un valor de la
columna correspondiente a interseccion = 1 ( recordar que los elementos de
intersección entre t, B y M tienen valor =1).
95
subset (pumiliointerseccion2, Names==“Vha16”)
Exploramos los elementos de pumiliointerseccion2 que contienen datos
acerca del gen de nombre “Vha16”.
Existen varios spots del gen Vha16, uno de los cuáles está en la
intersección de los estadísticos t, B y M.
96
Graficamos los genes de la intersección sobre MA y tA:
M1 <- pumiliofit$coefficients
Generamos un objeto M1 que contenga los datos de M (coefficients) de
pumiliofit.
A <- pumiliofit$Amean
Generamos un objeto A que contenga los datos de A (Amean) de pumiliofit.
par(mfrow=c(1,2))
Generamos un espacio de 2 gráficos.
plot(A,M1,pch='.', col="green")
Generamos un gráfico de MA con los objetos M1 y A.
points(A*interseccion,M1*interseccion, pch="*")
Generamos puntos * en los datos correspondientes a la intersección.
points(tablatop10Mpumilio$A[1],tablatop10Mpumilio$M[1],pch="*",col="blue",ce
x=2);text(9.206855,-6.772648,"CG7979“)
Generamos un punto * azul en el dato correspondiente al primer gen de
tablatop10Mpumilio.
plot(A, ordinary.t, ylab="t",pch=".",col="green")
Generamos un gráfico de tA.
points(A*interseccion,ordinary.t*interseccion , pch="*")
Generamos puntos * en los datos correspondientes a la intersección.
97
points(tablatop10tpumilio$A[1],tablatop10tpumilio$ordinary.t[1],pch="*",col="blue",c
ex=2);text(12.267057,4.738344,"CG2229")
Generamos un punto * azul en el dato correspondiente al primer gen de
tablatop10tpumilio.
98
Gracias por venir!!!!!!!!!!!!!!!!!!
Está claro que ahora no van a tener más
dudas……jajaja……….pero en caso de
que las tengan…………….
http://www.dm.uba.ar/materias/analisis_expl_y_conf_de_datos_de_exp_de
_marrays_Mae/2006/1/
Página de la materia Análisis Exploratorio y Confirmatorio de Datos de
Experimentos de Microarrays dictada por la Dra. Diana Kelmansky.
99