Download Análisis de datos de microarrays

Análisis de datos de
microarrays
Disenõ de Experimentos
Alex Sánchez-Pla y M. Carme Ruı́z de Villa
Departament d’Estadı́stica. Universitat de Barcelona.
Facultat de Biologia. Avda. Diagonal 643. 08028 Barcelona. Spain.
asanchez@ub.edu;mruiz_de_villa@ub.edu
xx
PID 00191030
Módulo 5
©
FUOC • PID 00191030 • Módulo XXX
Análisis de datos de microarrays
Índice
I
Preliminares
3
1
Estudio de un caso: Análisis básico usando R . . . . . . . . . . . . . . . . . .
4
1.1
Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
1.1.1
Carga de los datos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
1.1.2
1.2
Exploración de los datos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
Buscando los genes diferencialmente expresados . . . . . . . . . . . . . . . . .
15
1.2.1
Seleccionando los genes diferencialmente expresados . . . . .
17
1.2.2
Anotación de los resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
17
Part I
Preliminares
3
© FUOC • PID 00191030 • Módulo XXX
4
1. Estudio de un caso: Análisis básico usando R
.
1.1
Introducción
En este capı́tulo se llevará a cabo un ejemplo de cómo realizar un análisis de
microarrays usandoR. A la vez que se introducen los conceptos básicos del análisi
exploratorio (preliminar) en los microarrays se verá como puede usarse R para
hacer dicha exploración.
Los ejemplos que se presentaran más adelante en el curso se harán usando las
librerı́as para análisis de microarrays del proyecto Bioconductor
(www.bioconductor.org).
El ejercicio se llevará a cabo con datos de un estudio depositado en la base
de datos pública GEO (“Gene Expresion Omnibus”). El código del estudio es
GDS2744 y puede accederse a él a través del enlace
http://www.ncbi.nlm.nih.gov/sites/GDSbrowser?acc=GDS2744.
Los datos pueden descargarse de la dirección:
ftp://ftp.ncbi.nih.gov/pub/geo/DATA/SeriesMatrix/GSE7765/.
Para facilitar el trabajo los datos se han extraido y depositado en la carpeta “data”
que acopaña el material del texto. Consisten en un archivo de texto denominado
"GSE7765-GPL96-series-matrix.txt".
ˆ
GPL96 indica el tipo de array utilizado ( hgu133-A).La base de datos GEO
contiene un registro de “plataforma” para cada tipo de microarray, y de él
puede descargarse la tabla de las anotaciones en donde se indica, para cada
sonda el código “Entrez”, el sı́mbolo y otros identificadores habituales del gen
asociado. La tabla de anotaciones se encuentra también en la carpeta “data”
o puede descargarse de GEO:
ftp://ftp.ncbi.nih.gov/pub/geo/DATA/SeriesMatrix/GSE7765/GSE7765-GPL96 series matrix.txt.gz
ˆ
GSE7765 indica la “serie” de muestras que se estudian. Ésta contiene 12 muestras de 12 arrays, 6 hibridados con la plataforma “GPL96” –modelo hgu133a–,
los que se estudian en este ejemplo, y otros 6 con la plataforma “GPL97” –
modelo hgu133b-. En cada caso tres de las muestras han sido tratadas con
DMSO y otras tres con Dioxina-
Análisis de datos de microarrays
© FUOC • PID 00191030 • Módulo XXX
Análisis de datos de microarrays
5
Nota importante acerca del diseño experimental
Observando la información contenida en GEO ası́ como en el artı́culo original
(enlazado también desde la misma dirección) se deduce que aunque se trata de 6
muestras se han realizado tan sólo tres experimentos independientes, por lo que
en realidad se dispone de tres pares “control–tratamiento”.
Esta caracterı́stica, que se comprobará en e análisis exploratorio posterior, es
importante porque determina también que, para decidir si hay genes diferencialmente expresados será mejor hacerlo utilizando un test de datos apareados que
suponiendo que los controles y los tratamientos son independientes.
En lo que queda de ejercio se supondrá que las matrices de datos y las anotaciones
se encuentra en un subdirectorio “data” que pende del directorio de trabajo.
El archivo de anotaciones proporcionado se ha obtenido a partir del existente en
GEO pero se ha manipulado un poco para eliminar artefactos que aparecen al
importarlos.
1.1.1
Carga de los datos
Los datos tı́picos de microarrays consisten en una matriz de expresión génica que
tiene los genes en las filas y las muestras en las columnas. El archivo de datos que
se obtiene de GEO no es una tabla simple, sino que contiene primero algunas filas
de información general seguidas a partir de la lı́nea 67 de la matriz de expresión.
La lectura de los datos se lleva a cabo como de costumbre. Se crean dos objetos:
uno llamado “info”, con la información de la base de datos, y otro con la matriz
de expresión.
> # setwd(" ") # Escribir el directorio en el que se trabaja
> datadir <- "./data"
> info <-readLines(file.path(datadir,"GSE7765-GPL96_series_matrix.txt"), n=70)
> x <-read.table(file.path(datadir,"GSE7765-GPL96_series_matrix.txt"),
+ skip=68, header=TRUE, sep="\t",row.names=1, nrows = 22283)
> head(x)
GSM188013 GSM188014 GSM188016 GSM188018 GSM188020 GSM188022
1007_s_at 15630.200 17048.800 13667.500 15138.800 10766.600 15680.800
1053_at
3614.400
3563.220
2604.650
1945.710
3371.290
3406.660
117_at
1032.670
1164.150
510.692
5061.200
452.166
400.477
121_at
5917.800
6826.670
4562.440
5870.130
3869.480
3680.440
1255_g_at
224.525
395.025
207.087
164.835
111.609
130.123
1294_at
799.786
839.787
592.434
593.632
431.526
332.962
© FUOC • PID 00191030 • Módulo XXX
6
Si se inspecciona la información contenida en la cabecera del archivo o bien la información disponible en GEO puede verse que las muestras impares corresponden
a controles (DMSO) y las pares a tratamientos (Dioxina) a controles y tratamientos por lo que podemos etiquetar las columnas de x para trabajar cómodamente.
> colnames(x) <- paste(rep(c("DMSO","Diox"),3),substr(colnames(x),7,9),
+ sep="")
> colnames(x)
[1] "DMSO013" "Diox014" "DMSO016" "Diox018" "DMSO020" "Diox022"
Anotaciones
En general se asume que los nombres de las filas de una matriz de expresión
identifican las sondas, en este caso “probesets” o grupos de sondas de Affymetrix,
con las que se interrogan los genes. Los identificadores de las sondas deben de
asociarse de alguna forma con los de los genes, proceso que en general se denomina
“anotación”.
Para relacionar las sondas con los genes los desarrolladores de los chips –por
ejemplo la empres Affymetrix– suelen proporcionar archivos de anotaciones, que
no son más que tablas que asocian el nombre de cada sonda con algunos identificadores “oficiales”, donde por oficiales se entiende que corresponden a bases de
datos utilizadas por la comunidad.
La base de datos GEO contiene, ademas de datos de expresión, más de 9000
tablas de anotaciones, que denomina “platforms”. Si se conoce el nombre de la
plataforma puede descargarse y utilizarse para inspeccionar los genes interesantes
(de hecho cualquier gen) que aparecen en un estudio.
El archivo ,”GPL96-annotations” contiene las anotaciones para los arrays de tipo
HGU133-A utilizados en este estudio. Las primeras 26 lineas contienen una descripción de las columnas y las siguientes la anotación de cada grupo de sondas
en una fila distinta. Para simplificar el proceso hemos extraı́do las anotaciones
y las hemos grabado en el archivo “GPL96.annot.csv2” separado por puntos y
comas.
Asumiendo que los archivos se han descargado, descomprimido, y almacenado en
el directorio “data” podemos leerlas ambas con las instrucciones siguientes:
> infoAnots <-readLines(file.path(datadir,"GPL96.annot"), n=26)
> head(infoAnots, n=25)
Análisis de datos de microarrays
© FUOC • PID 00191030 • Módulo XXX
Análisis de datos de microarrays
7
[1] "^Annotation"
[2] "!Annotation_date = Jul 29 2010"
[3] "!Annotation_platform = GPL96"
[4] "!Annotation_platform_title = [HG-U133A] Affymetrix Human Genome U133A Array"
[5] "!Annotation_platform_organism = Homo sapiens"
[6] "#ID = ID from Platform data table"
[7] "#Gene title = Entrez Gene name"
[8] "#Gene symbol = Entrez Gene symbol"
[9] "#Gene ID = Entrez Gene identifier"
[10] "#UniGene title = Entrez UniGene name"
[11] "#UniGene symbol = Entrez UniGene symbol"
[12] "#UniGene ID = Entrez UniGene identifier"
[13] "#Nucleotide Title = Entrez Nucleotide title"
[14] "#GI = GenBank identifier"
[15] "#GenBank Accession = GenBank accession"
[16] "#Platform_CLONEID = CLONE_ID from Platform data table"
[17] "#Platform_ORF = ORF from Platform data table"
[18] "#Platform_SPOTID = SPOT_ID from Platform data table"
[19] "#Chromosome location = Entrez gene chromosome and location"
[20] "#Chromosome annotation = Entrez gene chromosome annotation"
[21] "#GO:Function = Gene Ontology Function term"
[22] "#GO:Process = Gene Ontology Process term"
[23] "#GO:Component = Gene Ontology Component term"
[24] "#GO:Function ID = Gene Ontology Function identifier"
[25] "#GO:Process ID = Gene Ontology Process identifier"
> annot <-read.csv2(file.path(datadir,"GPL96.annot.csv2"), head=T)
> dim(annot)
[1] 22283
21
> head(annot[1:5, c(1:4)])
ID
Gene.title Gene.symbol Gene.ID
1 1007_s_at discoidin domain receptor tyrosine kinase 1
DDR1
780
2
1053_at replication factor C (activator 1) 2, 40kDa
RFC2
5982
heat shock 70kDa protein 6 (HSP70B')
HSPA6
3310
3
117_at
4
121_at
paired box 8
PAX8
7849
5 1255_g_at
guanylate cyclase activator 1A (retina)
GUCA1A
2978
> rownames(annot)<-annot$ID
Con esta tabla es posible obtener información de cualquier gen seleccionado en
el estudio para el que se disponga de anotaciones.
© FUOC • PID 00191030 • Módulo XXX
1.1.2
Análisis de datos de microarrays
8
Exploración de los datos
En general la matriz de datos que se obtiene de GEO habrá sido preprocesada
o normalizada de alguna forma que estará descrita ela ficha del estudio. Para
explorar la matriz de expresión puede hacerse como con cualquier objeto de R.
> dim(x)
# Dimensions
[1] 22283
6
> names(x)
# Vector of strings to name the columns
[1] "DMSO013" "Diox014" "DMSO016" "Diox018" "DMSO020" "Diox022"
> round(apply(x,2, summary))
# Column-wise summary statistics,3)
DMSO013 Diox014 DMSO016 Diox018 DMSO020 Diox022
Min.
1st Qu.
1
1
2
2
1
1
327
317
192
233
172
169
Median
1146
1102
774
844
692
686
Mean
3459
3451
3364
3503
3438
3536
3rd Qu.
3052
3055
2968
2880
2853
2860
115400
112600
95910
115300
104000
119300
Max.
Para simplificar puede hacerse un attach de la matriz de datos de forma que no
haga falta usar “x$” para acceder a las columnas de x.
> attach(x)
The following object(s) are masked from 'x (position 4)':
Diox014, Diox018, Diox022, DMSO013, DMSO016, DMSO020
Los datos de expresión suelen trabajarse en escala logarı́tmica en base 2.
> logX <- log2(x)
> round(apply(x,2, summary),3)
# Column-wise summary statistics,3)
© FUOC • PID 00191030 • Módulo XXX
DMSO013
Min.
Diox014
Análisis de datos de microarrays
9
DMSO016
Diox018
DMSO020
Diox022
7.350e-01 8.870e-01
2.009 2.452e+00 5.470e-01 8.060e-01
1st Qu. 3.266e+02 3.173e+02
191.800 2.328e+02 1.723e+02 1.694e+02
Median
1.146e+03 1.102e+03
773.900 8.445e+02 6.917e+02 6.861e+02
Mean
3.459e+03 3.451e+03
3364.000 3.503e+03 3.438e+03 3.536e+03
3rd Qu. 3.052e+03 3.055e+03
2968.000 2.880e+03 2.853e+03 2.860e+03
Max.
1.154e+05 1.126e+05 95910.000 1.153e+05 1.040e+05 1.193e+05
> round(apply(logX,2, summary),3)
# Column-wise summary statistics,3)
DMSO013 Diox014 DMSO016 Diox018 DMSO020 Diox022
Min.
-0.443
-0.173
1.007
1.294
-0.871
-0.311
1st Qu.
8.352
8.310
7.584
7.863
7.429
7.405
Median
10.160
10.110
9.596
9.722
9.434
9.422
Mean
9.988
9.955
9.521
9.678
9.416
9.400
3rd Qu.
11.580
11.580
11.540
11.490
11.480
11.480
Max.
16.820
16.780
16.550
16.810
16.670
16.860
Una mirada a la distribución de los datos muestra como tomar logaritmos actúa
simetrizando la curva:
> opt <- par(mfrow=c(1,2))
> hist(x[,1], main="Escala original",
xlab="Expression",)
> hist(M.1 <- logX[,1], main="Escala logaritmica", xlab="Expresion")
> par(opt)
Escala logaritmica
1000
2000
Frequency
10000
0
0
5000
Frequency
15000
3000
20000
4000
Escala original
0
40000
100000
Expression
> opt<- par(mfrow=c(3,2))
> for(i in 1:6){
0
5
10
Expresion
15
© FUOC • PID 00191030 • Módulo XXX
+
+
10
hist(logX[,i], main=paste("Escala logaritmica. Array",i), xlab="Expresion")
}
> par(opt)
0
5
10
4000
15
0
5
10
15
Escala logaritmica. Array 3
Escala logaritmica. Array 4
Frequency
0
1500
1500
Expresion
3000
Expresion
0
5
10
15
5
10
15
Escala logaritmica. Array 5
Escala logaritmica. Array 6
Frequency
1500
0
0
5
10
2500
Expresion
3000
Expresion
0 1000
Frequency
2000
0
2000
Frequency
4000
Escala logaritmica. Array 2
0
Frequency
Escala logaritmica. Array 1
Frequency
Análisis de datos de microarrays
15
0
Expresion
5
10
15
Expresion
Una representación muy habitual cuando se usan arrays de dos colores es comparar las expresiones de los dos canales –llamémosles “R” y “G” mediante un diagrama de dispersión (R G) o mediante un diagrama “MA-plot” que representa la
log(R)+log(G)
R
∼
expresión relativa frente a la intensidad media de cada gen log G
2
(este gráfico se discute con detalle más adelante).
Con los arrays de un color pueden compararse dos arrays entre si pero lo que
suele hacerse es comparar cada array con un “array promedio” obtenido tomando
la mediana de las expresiones de cada gen en todos los arrays.
En primer lugar se crea este array promedio
> y<- medArray <- apply (x,1,median)
> logY<- log2(y)
A continuación tenéis un ejemplo de un scatterplot y un “MA-plot” del primer
array frente al array medio.
> opt<-par(mfcol=c(2,1))
> plot(X1<-x[,1], Y1<-medArray,
+
xlab="Intensidad (array 1)",
+
ylab="Intensidad (mediana)")
> abline(0,1, col="yellow")
> M.1<-log2(X1)-log2(Y1)
© FUOC • PID 00191030 • Módulo XXX
11
Análisis de datos de microarrays
> A.1 <- 0.5*(log2(X1)+log2(Y1))
> plot(A.1, M.1,
+
xlab="Intensidad media (A)",
+
ylab="Log Ratio (M)"
)
> abline (h=0,col="yellow")
8e+04
0e+00
Intensidad (array 2)
> par(opt)
0e+00
2e+04
4e+04
6e+04
8e+04
1e+05
5
0
−5
Log Ratio (M)
Intensidad (array 1)
5
10
15
Intensidad media (A)
Lo habitual es revisar totdos los arrays con el fin de detectar si alguno presenta
un patrón distinto de los demás.
> opt<- par(mfrow=c(3,2))
> for (i in 1:6){
+
M <-logX[,i]-logY
+
A <- 0.5*(logX[,i]+logY)
+
plot(A, M,
+
xlab="Intensidad media (A)",
+
ylab="Log Ratio (M)", main= colnames(x)[i])
+
abline (h=0,col="yellow")
+ }
>
par(opt)
Otros gráficos muy utilizados són los diagramas de caja –o “boxplots”– ideales
formarse una idea de si los datos necesitan un mayor preprocesado. La información del estudio indica que los arrays 1, 3 y 5 corresponden a un grupo y los
arrays 2, 4 y 6 a otro, por lo que se colorean de forma distinta.
> opt <- par(mfcol=c(1,2))
> boxplot(x, col=c(2,3,2,3,3,3),main="Valores de expresion\n 3 controles y 3 tratados",
© FUOC • PID 00191030 • Módulo XXX
12
Análisis de datos de microarrays
Figura 1. MAPlot.
Diox014
0
Log Ratio (M)
−6
−6
−4
−2
0
−2
−4
Log Ratio (M)
2
2
4
4
DMSO013
5
10
15
5
10
15
Intensidad media (A)
DMSO016
Diox018
−2
0
Log Ratio (M)
0
−2
−6
−4
−4
Log Ratio (M)
2
2
4
4
Intensidad media (A)
5
10
15
5
10
15
Intensidad media (A)
DMSO020
Diox022
0
Log Ratio (M)
−2
−4
−10 −8 −6 −4 −2
Log Ratio (M)
0
2
2
4
4
Intensidad media (A)
5
10
15
5
10
Intensidad media (A)
15
Intensidad media (A)
+ ylab="Expresion", las=2, cex.axis=0.7, cex.main=0.7)
> boxplot(logX, col=c(2,3,2,3,2,3), main="log(Valores de expresion) para\n
+ 3 controles y 3 tratados", ylab="Expresion", las=2, cex.axis=0.7, cex.main=0.7)
> abline(0,0, col="black")
> par(opt)
log(Valores de expresion) para
Valores de expresion
3 controles y 3 tratados
3 controles y 3 tratados
120000
15
100000
80000
Expresion
Expresion
10
60000
40000
5
20000
0
Diox022
Diox018
DMSO020
Diox014
DMSO016
DMSO013
Diox022
Diox018
DMSO020
Diox014
DMSO016
DMSO013
0
Transformaciones de los datos
Todos los gráficos realizados muestran una simetrı́a aceptable lo que sugiere que
los datos se pueden analizar relativamente bien. En capı́tulos siguientes se dicutirá
como transformar los datos adecuadamente para facilitar su análisis.
En general “transformar los datos” es algo que tiene ventajas e inconvenientes. Las
ventajas suelen consistir que los datos se colocan en una escala en la que puede ser
posible apreciar caracterı́sticas no visibles anteriormente. El inconveniente es que
© FUOC • PID 00191030 • Módulo XXX
13
debemos de asegurarnos que dichas caracterı́sticas son propias de los datos y no
de la transformación. EN consecuencia el consejo general suele ser “transformar
si hace falta, pero transformar lo mı́nimo posible”.
En este ejemplo los boxplots muestran que los valores medios de cada array
son ligeramente diferentes. Esto no se debe a la biologı́a del tema sino que es
probablemente un artefecto técnico (unos arrays se iluminan más que otros, por
ejemplo). EN consecuencia para asegurarnos de que este artefacto no afectará a
una eventual comparación entre arrays decidimos centrar cada columna, lo que
hacemis restándole su media. Esta transformación no afectará a la distribución
de las expresiones pero facilitará que se puedan encontrar diferencias entre los
grupos.
> logX <- scale(logX, scale=F)
Exploración multivariante
En general se espera que microarrays pertenecientes a grupos homogéneos esten
muy correlacionados entre si, aunque artefactos técnicos – derivados del proceso
de producción – pueden introducir un sesgo que afecte esta homogeneidad.
Una forma de comprobarlo es aplicar alguna técnica de visualización como un
análisis de componentes principales o un dendrograma basado en una matriz de
distancias entre los arrays.
El análisis de componentes principales
El análisis de componentes prin-
cipales o PCA permite detectar si las muestras se agrupan por su similaridad –
relacionada con el grupo al que pertenecen – o bien se producen agrupaciones
inesperadas (o falta de ellas), lo que puede ser indicativo – aunque no necesariamente – de posibles efectos técnicos.
> pcX<-prcomp(t(logX), scale=FALSE) # Ya se han escalado los datos
> loads<- round(pcX$sdev^2/sum(pcX$sdev^2)*100,1)
Una vez calculadas las componentes principales se suele realizar la representación
de los datos en las dos o tres primeras componentes. En la representación gráfica
siempre es conveniente mostrar el porcentaje de variabilidad explicada por cada
componente. A menor porcentaje menor relavancia de la clasificación.
> xlab<-c(paste("PC1",loads[1],"%"))
> ylab<-c(paste("PC2",loads[2],"%"))
> plot(pcX$x[,1:2],xlab=xlab,ylab=ylab, xlim=c(-150, 150))
Análisis de datos de microarrays
© FUOC • PID 00191030 • Módulo XXX
14
> title("Componentes principales (PCA)")
> text(pcX$x[,1],pcX$x[,2],colnames(x), pos=4)
Componentes principales (PCA)
20
40
Diox022
DMSO020
DMSO013
0
Diox014
−20
−40
−80
−60
PC2 18.4 %
DMSO016
Diox018
−150
−100
−50
0
50
100
150
PC1 48.7 %
Sin embargo, los dendrogramas de un análisis de cluster también resultan informativos.
> clust.euclid.average <- hclust(dist(t(x)),method="average")
> plot(clust.euclid.average, hang=-1)
> par(opt)
DMSO020
DMSO016
Diox022
Diox014
DMSO013
Diox018
150000
Height
250000
350000
Cluster Dendrogram
dist(t(x))
hclust (*, "average")
Los dos gráficos anteriores sugieren que las muestras no se agrupan por tratamiento
sinó que se agrupan por experimento lo que confirma la información de que se
trata de tres pares de muestras obtenidas en tres experimentos independientes.
Análisis de datos de microarrays
© FUOC • PID 00191030 • Módulo XXX
1.2
15
Buscando los genes diferencialmente expresados
En esta sección se asume que los datos ya han sido normalizados, de modo que
podemos proceder a compararlos a ver si hay diferencias en las expresiones genicas
entre dos condiciones.
Existen paquetes que permiten la aplicación de las diferentes opciones de tests
de expresión diferencial, ası́ como las correcciones para el problema del multiple
testing.
La matriz contiene los datos de dos tipos de lı́neas celulares: tratamiento y control.
Esto datos pueden ser almacenados en un vector de etiquetas: logX.cl que indica
cúal de los chips pertenece al grupo de tratamiento y cúal al grupo control.
> logX.cl<-c(0,1,0,1,0,1)
El siguiente bucle realiza una comparación de medias entre los dos grupos. Esta
no es la mejor prueba que se podrı́a hacer aquı́. No obstante, mostramos este test
con fines ilustrativos.
> get.fc=function(x){sum(x[logX.cl==1]-x[logX.cl==0])/3}
> fcs=apply(logX,1,get.fc)
Un enfoque más razonable para comparar estos datos serı́a utilizar un test t–
Student, usando por ejemplo la función t.test. Podemos adaptarlo a nuestras
necesidades, teniendo en cuenta que los datos pueden considerarse apareados.
> ttest=function(x){
+ tt=t.test(x[logX.cl==0],x[logX.cl==1], paired=TRUE)
+ return(c(tt$statistic,tt$p.value))}
Una vez preparada la función se aplica a cada gen usando apply:
> ans=apply(logX,1,ttest)
> ts<- ans[1,]
> pvals<-ans[2,]
Podemos estudiar la distribución de los valores de t obtenidos utilizando un histograma o un gráfico “qq-plot”. En todo caso no debemos olvidar que estamos
observando la distribución de los estadı́sticos de test, no de los datos. El hecho
por ejemplo de que en en el “qq-plot” aparezcan valores por encima o debajo de
Análisis de datos de microarrays
© FUOC • PID 00191030 • Módulo XXX
16
la lı́nea aparte de sugerir desviaciones de la distribución “t” sugiere que podrı́an
aparecer muchos genes cuya expresión sea distinta entre ambas condiciones (es
decir con valores muy altos o muy bajos de t).
> hist(ts, breaks=40)
6000
4000
0
2000
Frequency
8000
10000
Histogram of ts
−100
−50
0
50
100
2
4
ts
> qqnorm(ts)
> qqline(ts)
0
−50
−100
Sample Quantiles
50
100
Normal Q−Q Plot
−4
−2
0
Theoretical Quantiles
Análisis de datos de microarrays
© FUOC • PID 00191030 • Módulo XXX
1.2.1
Análisis de datos de microarrays
17
Seleccionando los genes diferencialmente expresados
Un gen se denomina diferencialmente expresado si su expresión es significativamente diferente entre dos condiciones. Ası́, podemos decidir que hay “expresión
diferencial” si el p–valor del test es menor que un valor fijado. Normalmente este
valor será menor que el habitual 0.05 o 0.01 debido al problema de multiplicidad
de tests que también se discutirá más adelante.
> for (i in c(0.01, 0.001, 0.0001, 0.00001, 0.000001, 0.0000001))
+
print(paste("genes con p-valores menores que", i, length(which(pvals < i))))
[1] "genes con p-valores menores que 0.01 248"
[1] "genes con p-valores menores que 0.001 20"
[1] "genes con p-valores menores que 1e-04 2"
[1] "genes con p-valores menores que 1e-05 0"
[1] "genes con p-valores menores que 1e-06 0"
[1] "genes con p-valores menores que 1e-07 0"
En vez de basarnos tan sólo en el p valor podemos seleccionar los genes que tengan
un p-valor no ajustado pequeño y una diferencia de medias (“fold–change”) de
como mı́nimo 1.
> selGenes<- which(pvals<0.01 & fcs > 1)
> print(selGenes)
201195_s_at 204939_s_at
205623_at
206285_at
4466
5150
5811
6281
6513
208409_at 211531_x_at
212665_at
213839_at
213966_at
214235_at
13219
13345
13614
217686_at 218691_s_at
220505_at
723
7910
10937
12050
214852_x_at
214967_at
216588_at
14227
14342
15958
17051
220880_at 221880_s_at
51158_at
220812_s_at
20176
1.2.2
20244
21240
206755_at 206987_x_at
18055
19869
22068
Anotación de los resultados
Obviamente los nombres de los grupos de sondas no son muy informativos. Ahora
podemos recurrir a la tabla de anotaciones para saber más acerca de estos genes.
La tabla de anotaciones contiene en su primera columna los nombres de las sondas. Si los utilizamos para localizar los genes seleccionados podremos disponer
de información sobre éstos.
© FUOC • PID 00191030 • Módulo XXX
18
Análisis de datos de microarrays
> selAnots<- annot[names(selGenes),]
Para concluir podemos preparar una tabla resumen de los resultados en donde
combinemos los identificadores de los genes, los estadı́sticos de test y sus anotaciones.
> selFC <- fcs[selGenes]
> seltt <- ts[selGenes]
> selp
<- pvals[selGenes]
> selGeneTitle<- substr(selAnots[,"Gene.title"],1,25)
> selGeneSymbol <- substr(selAnots[,"Gene.symbol"],1,10)
> selGeneID <- substr(selAnots[,"Gene.ID"],1,10)
> result<- data.frame(foldChg=selFC, ttest= seltt, pvalue=selp,
+
Gene.Name=selGeneTitle, Gene.Symbol=selGeneSymbol, Gene.ID=selGeneID)
> result<- result[order(result$pvalue),]
> if(!(require(xtable))) install.packages("xtable")
> require(xtable)
> print(xtable(result))
© FUOC • PID 00191030 • Módulo XXX
Análisis de datos de microarrays
19
foldChg
ttest
pvalue
Gene.Name
Gene.Symbol
Gene.ID
216588 at
1.47
-33.95
0.00
ribosomal protein L7 pseu
RPL7P12///
100131085/
204939 s at
2.48
-31.80
0.00
phospholamban
PLN
5350
201195 s at
2.43
-25.64
0.00
solute carrier family 7 (
SLC7A5
8140
208409 at
1.29
-24.33
0.00
solute carrier family 14
SLC14A2
8170
214967 at
2.16
-22.90
0.00
218691 s at
1.87
-22.80
0.00
PDZ and LIM domain 4
PDLIM4
8572
211531 x at
1.02
-21.57
0.00
proline-rich protein BstN
PRB2///PRB
653247///5
206987 x at
1.43
-21.12
0.00
fibroblast growth factor
FGF18
8817
212665 at
1.55
-18.94
0.00
TCDD-inducible poly(ADP-r
TIPARP
25976
221880 s at
1.14
-18.83
0.00
family with sequence simi
FAM174B
400451
220880 at
1.21
-16.41
0.00
51158 at
1.30
-15.76
0.00
family with sequence simi
FAM174B
400451
213966 at
1.97
-15.60
0.00
high-mobility group 20B
HMG20B
10362
220812 s at
1.53
-15.04
0.00
HERV-H LTR-associating 2
HHLA2
11148
214852 x at
1.16
-14.40
0.00
vacuolar protein sorting
VPS13A
23230
214235 at
1.58
-14.37
0.00
cytochrome P450, family 3
CYP3A5
1577
213839 at
1.02
-13.89
0.01
calmin (calponin-like, tr
CLMN
79789
206285 at
1.15
-12.73
0.01
nephronophthisis 1 (juven
NPHP1
4867
205623 at
3.19
-12.16
0.01
aldehyde dehydrogenase 3
ALDH3A1
218
206755 at
1.25
-10.54
0.01
cytochrome P450, family 2
CYP2B6
1555
217686 at
1.61
-10.35
0.01
protein tyrosine phosphat
PTPN1
5770
220505 at
1.05
-10.12
0.01
chromosome 9 open reading
C9orf53
51198