Capítulo 5 Exploración de los datos, control de calidad y preprocesado
5.1 Introducción
Los estudios realizados con microarrays, sea cual sea la tecnología en que se basan, tienen una característica común: generan grandes cantidades de datos a través de una serie de procesos, que hacen que su significado no siempre sea completamente intuitivo.
Como en todo tipo de análisis, antes de empezar a trabajar con los datos, debemos de asegurarnos de que éstos son fiables y completos y de que se encuentran en la escala apropiada para proporcionar la información que pretendemos obtener de ellos.
En el caso de los microarrays solemos distinguir dos fases previas al análisis de los datos:
Exploración y control de calidad Mediante gráficos y/o resúmenes numéricos estudiamos la estructura de los datos con el fin de decidir si parecen correctos o presentan anomalías que deben ser corregidas.
Preprocesado Incluso si son correctos, los datos “crudos” no sirven para el análisis sino que necesitan preprocesados, lo que puede interpretarse como uno o más de los procesos siguientes:
- “Limpiados,” para eliminar la parte de la señal no atribuible a la expresión, llamada “background” o ruído.
- “Normalizados” para hacerlos comparables entre muestras y eliminar posibles sesgos técnicos.
- Resumidos o “sumarizados” de forma que se tenga un sólo valor por gen.
- “Transformados” de forma que la escala sea razonable y facilite el análisis.
A menudo -por un cierto abuso de lenguaje- se denomina “normalización” al preprocesado descrito en las etapas anteriores.
5.1.1 Nivel de análisis y tipo de microarray
Desde la generalización del uso de los microarrays se han desarrollado muchas formas para visualizar los datos y decidir acerca de su calidad. Algunas trabajan sobre los datos obtenidos del escáner, otras lo hacen con los datos normalizados. Algunas sirven tanto para arrays de dos colores como arrays de un color. Otras son específicas de la tecnología. Con el fin de organizar esta multitud de opciones podemos diferenciar:
- Datos de bajo nivel Son los datos proporcionados por el escáner contenidos por ejemplo en archivos .gpr (dos colores) o .cel (un color). Estos últimos son binarios por lo que no es posible ni tan sólo visualizarlos sin programas específicos.
- Datos de alto nivel Son los datos resultantes del (pre)procesado de los datos de bajo nivel. Básicamente se corresponden con los datos de expresión ya resumidos (“sumarizados”), normalizados o no.
La exploración y el control de calidad pueden basarse en datos de bajo o de alto nivel. En cada caso se pueden aplicar ya sean técnicas generales de visualización de datos, como histogramas, diagramas de caja o visualización en dimensiones reducidas (PCA u otras), o bien técnicas “ad-hoc” para cada tipo de datos como el MA-Plot u otras que se discuten a continuación.
5.1.2 Datos de partida
La estructura de los datos de microarrays de un color o de dos colores difiere considerablemente, tanto a nivel físico (los chips y las imágenes que de ellos se obtiene) como informático, es decir en la forma en que se representan.
5.1.2.1 Arrays de dos colores (o de “cDNA”)
Tradicionalmente los arrays de dos colores o de cDNA se realizaban de forma menos automatizada que los de un color o de Affymetrix. Esto implica que, tras obtener la imagen, el escaneado del archivo “.TIFF” resultante pudiera ser llevado a cabo mediante un software independiente como Genepix, un programa que convierte las imágenes en números y genera un archivo de información (con extensión “.gpr”) a partir del cual pueden calcularse las expresiones relativas, así como valores de calidad para cada “spot” o punto escaneado en la imagen.
Para cada imagen (o sea para cada microarray) hay un archivo “.gpr” que contiene una fila por gen y varias columnas con distintos valores, por ejemplo:
- la intensidad para cada canal,
- valores resumen de las intensidades y
- resultados de controles de calidad (“FLAGS”).
La tabla siguiente muestra lo que serían las primeras filas y columnas de un archivo “.gpr” obtenido mediante el programa Genepix.
Gen | Señal | Fondo | Señal | Fondo | “Flag” | Otros |
---|---|---|---|---|---|---|
Cy5 (R) | Cy5 (Rb) | Cy3 (G) | Cy3 (Gb) | |||
gen-1 | 3.7547 | 1.8128 | 5.0672 | 1.8496 | 1 | … |
gen-2 | 0.8331 | 0.9175 | 1.1536 | 0.9995 | 0 | … |
gen-3 | 9.8254 | 0.2781 | 0.6921 | 0.5430 | 1 | … |
gen-4 | 9.1539 | 0.1918 | 3.8290 | 0.0014 | 0 | … |
gen-5 | 4.8603 | 0.2377 | 0.5338 | 0.3335 | 0 | … |
gen-6 | 7.8567 | 1.3941 | 1.3050 | 1.4876 | 1 | … |
gen-7 | 2.7619 | 0.8108 | 8.4916 | 1.6518 | 0 | … |
gen-8 | 3.6618 | 0.1918 | 0.3770 | 1.1842 | 0 | … |
gen-9 | 4.4300 | 0.5888 | 2.8161 | 0.8707 | 1 | … |
… | … | … | … | … | … | … |
Los valores de intensidad se convierten en una única matriz de expresión que contiene una columna por chip con los valores de intensidad relativa obtenidas por ejemplo con una “sumarización” del tipo: \[ \log\frac{R-Rb}{G-Gb}, \] y una fila por gen (mismas filas que archivos “.gpr”).
La tabla siguiente muestra lo que sería una matriz de expresión derivada de cuatro archivos “.gpr” como los de la tabla anterior.
gen | Array-1 | Array-2 | Array-3 | Array-4 |
---|---|---|---|---|
gen-1 | 1.65695 | -0.10820 | 1.69515 | 8.25137 |
gen-2 | -1.82305 | 0.11350 | 1.58807 | 0.95676 |
gen-3 | 0.01561 | 0.47682 | -27.88036 | -1.39078 |
gen-4 | 0.42709 | 4.29319 | -4.31366 | 30.96866 |
gen-5 | 0.04332 | 0.24126 | -10.27675 | 0.26680 |
gen-6 | -0.02825 | 0.68408 | 2.04163 | 0.99554 |
gen-7 | 3.50549 | -8.04635 | 0.18286 | 0.22647 |
gen-8 | -0.23261 | -1.08477 | 0.19582 | 1.11561 |
gen-9 | 0.50643 | 1.52147 | 0.99092 | 0.23441 |
… | … | … | … | … |
5.1.2.2 Arrays de un color (Affymetrix)
El resultado de escanear la imagen de un array de affymetrix es un archivo de extensión “.CEL” que, a diferencia de los arrays de dos colores, está en formato binario es decir que solo puede ser leído con programas específicos para ello.
De forma similar a los arrays de dos colores, existe un archivo “.CEL” por cada microarray.
A partir de las intensidades de los archivos “.CEL” se genera la matriz de expresión, que contiene una columna por chip con los valores de intensidad absoluta, y una fila por grupo de sondas. En el caso de arrays de affymetrix existe una gran variedad de algoritmos de “sumarización” y, según cual se utilice, se obtendrá unos u otros valores de expresión. Ahora bien, éstos serán siempre medidas absolutas, es decir independientes del resto de muestras y en una escala arbitraria.
La tabla siguiente muestra las primeras filas de una matriz de expresión sumarizada correspondiente a los primeros genes de uno de los casos resueltos.
ID_REF | GSM188013 | GSM188014 | GSM188016 | GSM188018 |
---|---|---|---|---|
1007_s_at | 15630.2 | 17048.8 | 13667.5 | 15138.8 |
1053_at | 3614.4 | 3563.22 | 2604.65 | 1945.71 |
117_at | 1032.67 | 1164.15 | 510.692 | 5061.2 |
121_at | 5917.8 | 6826.67 | 4562.44 | 5870.13 |
1255_g_at | 224.525 | 395.025 | 207.087 | 164.835 |
… | … | … | … | … |
5.1.2.3 Datos de ejemplo
Los ejemplos de este capítulo se basarán en dos de los conjuntos de datos descritos en el capítulo 3.
- Por un lado, trabajaremos con los datos del conjunto referidos al efecto del tratamiento con estrógenos en cancer de mama descrito en @ref{estrogen}. La exploración se basará en los datos normalizados pero también en los datos “crudos,” de bajo nivel, contenidos en los archivos “.cel.” Estos datos pueden obtenerse directamente del paquete de Bioconductor . Un aspecto interesante de este conjunto de datos es que, junto con 8 muestras “buenas” se proporciona un array defectuoso, denominado “” que facilita el ver como aparecen ciertos gráficos cuando hay problemas.
- Por otro lado, usaremos los datos del conjunto , relativos al efecto del tratamiento con CCL4 en la expresión de los hepatocitos. Los datos resultan accesibles al instalar el paquete .
if(!require(BiocManager)) install.packages("BiocManager")
if (!(require(CCl4))){
::install("CCl4")
BiocManager
}if (!(require(estrogen))){
::install("estrogen")
BiocManager
}if (!(require(affy))){
::install("affy")
BiocManager
}if (!(require(affyPLM))){
::install("affyPLM")
BiocManager }
Aunque, en la actualidad (año 2021) el paquete más utilizado para arrays de este tipo es el paquete oligo
este ejemplo utiliza el paquete affy
.
library(estrogen)
library(affy)
<- system.file("extdata", package = "estrogen")
affyPath = read.AnnotatedDataFrame("phenoData.txt", sep="", path=affyPath)
adfAffy = pData(adfAffy)
affyTargets $filename = file.path(affyPath, row.names(affyTargets))
affyTargets<- read.affybatch(affyTargets$filename, phenoData=adfAffy)
affyRaw # show(affyRaw)
<- getwd()
actualPath setwd(affyPath)
<- ReadAffy()
allAffyRaw setwd(actualPath)
El resultado de la lectura es un objeto “affyraw” de clase “affyBatch”:
class(affyRaw)
## [1] "AffyBatch"
## attr(,"package")
## [1] "affy"
print(affyRaw)
## AffyBatch object
## size of arrays=640x640 features (22 kb)
## cdf=HG_U95Av2 (12625 affyids)
## number of samples=8
## number of genes=12625
## annotation=hgu95av2
## notes=
El paquete limma
es conocido como paquete para la selección de genes diferencialmente expresados pero también incluye algunas funciones para la lectura y preprocesado de arrays de dos colores.
library("limma")
library("CCl4")
= system.file("extdata", package="CCl4")
dataPath = read.AnnotatedDataFrame("samplesInfo.txt",
adf path=dataPath)
#adf
= pData(adf)
targets $FileName = row.names(targets)
targets<- read.maimages(targets, path=dataPath, source="genepix")
RG attach(RG$targets)
<-paste(substr(Cy3,1,3),substr(Cy5,1,3),substr(FileName,10,12), sep="")
newNames colnames(RG$R)<-colnames(RG$G)<-colnames(RG$Rb)<-colnames(RG$Gb)<-rownames(RG$targets)<- newNames
# show(RG)
El resultado de la lectura es un objeto de classe “RGlist.”
class(RG)
## [1] "RGList"
## attr(,"package")
## [1] "limma"
print(RG)
## An object of class "RGList"
## $R
## DMSCCl319 DMSCCl320 DMSCCl321 DMSCCl329 CClDMS330 CClDMS331 CClDMS332
## [1,] 47 46 40 59 56 40 45
## [2,] 1095 46 39 54 56 40 47
## [3,] 47 46 41 55 52 39 46
## [4,] 46 45 41 56 54 43 47
## [5,] 55 53 45 61 71 47 63
## DMSCCl333 CClDMS379 CClDMS380 CClDMS381 DMSCCl382 DMSCCl384 DMSCCl389
## [1,] 70 46 65 43 59 46 63
## [2,] 69 45 61 43 62 48 68
## [3,] 81 44 57 41 58 50 63
## [4,] 78 44 60 42 59 48 62
## [5,] 91 51 83 56 71 52 69
## DMSCCl390 CClDMS391 CClDMS393 CClDMS394
## [1,] 77 40 91 32
## [2,] 71 41 85 32
## [3,] 70 40 39 32
## [4,] 81 41 51 32
## [5,] 81 57 95 37
## 44285 more rows ...
##
## $G
## DMSCCl319 DMSCCl320 DMSCCl321 DMSCCl329 CClDMS330 CClDMS331 CClDMS332
## [1,] 77 78 45 45 87 65 128
## [2,] 4558 79 44 47 81 66 131
## [3,] 79 75 46 47 84 66 122
## [4,] 78 82 44 71 86 63 135
## [5,] 93 95 92 57 86 66 132
## DMSCCl333 CClDMS379 CClDMS380 CClDMS381 DMSCCl382 DMSCCl384 DMSCCl389
## [1,] 107 57 90 99 84 56 68
## [2,] 110 54 89 96 168 54 72
## [3,] 114 54 84 93 78 55 71
## [4,] 113 118 84 97 83 56 75
## [5,] 133 59 95 99 102 64 84
## DMSCCl390 CClDMS391 CClDMS393 CClDMS394
## [1,] 71 74 83 59
## [2,] 71 75 81 60
## [3,] 73 78 62 59
## [4,] 79 77 70 61
## [5,] 99 83 92 63
## 44285 more rows ...
##
## $Rb
## DMSCCl319 DMSCCl320 DMSCCl321 DMSCCl329 CClDMS330 CClDMS331 CClDMS332
## [1,] 41 39 35 47 44 36 40
## [2,] 43 39 35 46 44 36 39
## [3,] 41 40 35 48 44 35 39
## [4,] 40 40 35 49 43 35 39
## [5,] 41 40 36 48 44 35 39
## DMSCCl333 CClDMS379 CClDMS380 CClDMS381 DMSCCl382 DMSCCl384 DMSCCl389
## [1,] 55 39 49 37 46 42 50
## [2,] 55 37 50 37 46 41 52
## [3,] 59 39 47 37 46 41 52
## [4,] 63 39 48 36 46 41 52
## [5,] 60 39 47 37 46 41 52
## DMSCCl390 CClDMS391 CClDMS393 CClDMS394
## [1,] 60 36 34 30
## [2,] 59 36 35 31
## [3,] 55 35 35 30
## [4,] 58 36 34 30
## [5,] 56 36 35 30
## 44285 more rows ...
##
## $Gb
## DMSCCl319 DMSCCl320 DMSCCl321 DMSCCl329 CClDMS330 CClDMS331 CClDMS332
## [1,] 44 44 32 34 40 39 48
## [2,] 46 43 32 35 41 38 48
## [3,] 43 43 32 35 40 38 49
## [4,] 43 43 32 35 41 38 49
## [5,] 43 42 32 35 41 38 48
## DMSCCl333 CClDMS379 CClDMS380 CClDMS381 DMSCCl382 DMSCCl384 DMSCCl389
## [1,] 47 35 50 50 45 37 38
## [2,] 48 35 50 49 45 37 38
## [3,] 48 35 51 49 45 37 38
## [4,] 49 35 51 49 46 37 37
## [5,] 49 35 50 49 45 37 38
## DMSCCl390 CClDMS391 CClDMS393 CClDMS394
## [1,] 38 40 37 36
## [2,] 38 40 37 36
## [3,] 39 40 37 36
## [4,] 38 40 37 36
## [5,] 39 40 37 36
## 44285 more rows ...
##
## $targets
## Cy3 Cy5 RIN.Cy3 RIN.Cy5 FileName
## DMSCCl319 DMSO CCl4 9.0 9.7 251316214319_auto_479-628.gpr
## DMSCCl320 DMSO CCl4 9.0 5.0 251316214320_auto_478-629.gpr
## DMSCCl321 DMSO CCl4 9.0 2.5 251316214321_auto_410-592.gpr
## DMSCCl329 DMSO CCl4 9.0 2.5 251316214329_auto_429-673.gpr
## CClDMS330 CCl4 DMSO 9.7 9.0 251316214330_auto_457-658.gpr
## 13 more rows ...
##
## $genes
## Block Row Column ID Name
## 1 1 1 1 BrightCorner BrightCorner
## 2 1 1 2 BrightCorner BrightCorner
## 3 1 1 3 (-)3xSLv1 NegativeControl
## 4 1 1 4 A_44_P317301 AW523361
## 5 1 1 5 A_44_P386163 NM_001007719
## 44285 more rows ...
##
## $source
## [1] "genepix"
##
## $printer
## $ngrid.r
## [1] 1
##
## $ngrid.c
## [1] 1
##
## $nspot.r
## [1] 430
##
## $nspot.c
## [1] 103
5.2 Exploración y control de calidad
Los gráficos son útiles para comprobar la calidad de los datos de microarrays, obtener información sobre cómo se deben preprocesar los datos y comprobar, finalmente, que el preprocesado se haya realizado correctamente.
Siguiendo el esquema presentado en la tabla @ref{tab:tablaDiagnosticos} se presentan a continuación los distintos gráficos utilizados con una breve descripción de lo que representa cada uno y como interpretarlos adecuadamente. El código para generarlos se presenta al final del capítulo como un apéndice.
5.2.1 Control de calidad con gráficos estadísticos generales
5.2.1.1 Histogramas y gráficos de densidad
Estos gráficos permite hacerse una idea de si las distribuciones de los distintos arrays son similares en forma y posición. La figura 5.1 muestra los histogramas correspondientes a los 9 arrays del conjunto de datos .
<- rownames(pData(allAffyRaw))
affySampleNames <- c(1,2,2,3,3,4,4,8,8)
affyColores <- c(1,2,2,2,2,3,3,3,3)
affyLineas hist(allAffyRaw, main="Signal distribution", col=affyColores, lty=affyLineas)
legend (x="topright", legend=affySampleNames , col=affyColores, lty=affyLineas, cex=0.7)
Tanto los gráficos de densidad como los diagramas de cajas de la figura siguiente muestran un desplazamiento global de los valores en todas las muestras del grupo de 48 horas lo que sugiere un efecto técnico, puesto que no es biológicamente plausible que todas los genes se vean más afectados por la mayor tiempo de exposición.
5.2.1.2 Diagramas de caja o “boxplots”
Como los histogramas los diagramas de caja –basados en los distintos cuantiles de las valores– dan una idea de la distribución de las intensidades. La figura 5.2 muestra los diagramas de caja correspondientes a los 9 arrays del conjunto de datos .
boxplot(allAffyRaw, main="Distribución de las expresiones", col=affyColores, las=2)
5.2.1.3 Gráficos de componentes principales
El análisis de componentes principales puede servir para detectar si las muestras se agrupan de forma “natural” es decir con otras muestras provenientes del mismo grupo o si no hay correspondencia clara entre grupos experimentales y proximidad en este gráfico. Cuando esto sucede no significa necesariamente que haya un problema pero puede ser indicativo de efectos técnicos -como el conocido efecto “batch”- que podría ser necesario corregir.
<- function ( X, labels=NULL, colors=NULL, dataDesc="", scale=FALSE)
plotPCA
{<-prcomp(t(X), scale=scale) # o prcomp(t(X))
pcX<- round(pcX$sdev^2/sum(pcX$sdev^2)*100,1)
loads<-c(paste("PC1",loads[1],"%"))
xlab<-c(paste("PC2",loads[2],"%"))
ylabif (is.null(colors)) colors=1
plot(pcX$x[,1:2],xlab=xlab,ylab=ylab, col=colors,
xlim=c(min(pcX$x[,1])-10, max(pcX$x[,1])+10),
ylim=c(min(pcX$x[,2])-10, max(pcX$x[,2])+10),
)text(pcX$x[,1],pcX$x[,2], labels, pos=3, cex=0.8)
<- ifelse(dataDesc=="", "Visualización de las dos primeras componentes", dataDesc)
titulo title(titulo, cex=0.8)
}
if (!file.exists("datos/affyNorm")){
<- rma(allAffyRaw)
allAffyNorm<- rma(affyRaw)
affyNorm save(allAffyNorm, affyNorm, file="datos/affyNorm.Rda")
else{
}load(file="datos/affyNorm.Rda")
}
Las figuras 5.3 y 5.4 muestran dos diagramas de componentes principales realizados a partir de los datos normalizados del conjunto de datos . El gráfico de la parte superior que incluye el array defectuoso ilustra que la principal fuente de variabilidad es la diferencia de este array con el resto. Cuando se repite el análisis omitiendo esta muestra puede verse como la principal fuente de variación (eje X) se asocia con el tiempo de exposición (alto a la derecha, bajo (10h) a la izquierda, mientras que la segunda fuente de variación se asocia con la exposición a los estrógenos (alto arriba, bajo abajo).
plotPCA(exprs(allAffyNorm), labels=affySampleNames)
plotPCA(exprs(affyNorm), labels=colnames(exprs(affyNorm)))
5.2.1.4 Imagen del chip
Otra forma de ver si las muestras se agrupan según los grupos experimentales, o mediante otros criterios es usando un cluster jerárquico que realiza una agrupación básica de las muestras por grado de similaridad según la distancia que se utilice.
Como en el caso de las componentes principales si las muestras se agrupan según las condiciones experimentales es una buena señal pero si no es así puede deberse a la presencia de otra fuente de variación o bien al hecho de que se trata de un gráfico basado en todo los datos y las condiciones experimentales pueden haber afectado un pequeño número de genes.
La figura 5.5 muestra como se agrupan los datos del conjunto en base a un cluster jerárquico. Como en el caso de las componentes principales tras eliminar el array defectuoso las muestras se separan, primero por el tiempo de exposición y luego por niveles de estrógeno suministrado.
<- hclust(dist(t(exprs(affyNorm))),method="average")
clust.euclid.average plot(clust.euclid.average, labels=colnames(exprs(affyNorm)), main="Hierarchical clustering of samples", hang=-1, cex=0.7)
5.2.2 Gráficos de diagnóstico para microarrays de dos colores
El diagnóstico de arrays de dos canales se basa principalmente en la imagen y en diferentes tipos de gráficos.
5.2.2.1 Diagramas de dispersión y “MA-plots”
La normalización discutida en este mismo capítulo es un punto clave en el proceso de análisis de microarrays y se ha dedicado un gran esfuerzo a desarrollar y probar diferentes métodos (Quackenbush (2002),Yang et al. (2002)). Una razón para ello es que hay diferentes artefactos técnicos que deben ser corregidos para poder ser utilizados, y no cualquier método puede funcionar con todos ellos.
En general, los métodos de normalización se basan en el siguiente principio: _La mayor parte de los genes en un array se pueden expresar o no expresar ante cualquier condición, pero se espera que sólo una pequeña cantidad de genes muestre cambios de expresión entre condiciones}.
Esto da una idea de como debería ser un gráfico de intensidades. Por ejemplo, si no hubiese artefactos técnicos, en arrays de dos canales, una gráfica de dispersión de intensidad del rojo frente al verde debería dejar la mayor parte de los puntos alrededor de una diagonal. Cualquier desviación de esta situación debería ser atribuible a razones técnicas, no biológicas, y por tanto, debería ser eliminada. Esto ha conducido a un método de normalización muy popular consistente en estimar la transformación a aplicar, como una función de las intensidades utilizando el método lowess en la representación transformada de la gráfica de dispersión conocida como el gráfico MA–plot.
La figura 5.6 muestra, en su parte izquierda un gráfico de dispersión del canal
rojo frente al verde en un array de dos colores. En concreto se trata del primer array del dataset CCl4
presentado en la introducción.
El hecho de que los datos estén centrados alrededor de la diagonal sugiere que puede no ser necesario normalizar los datos.
Una representación muy popular que ayuda a visualizar mejor esta relación es lo que se conoce como MA-plot, que aparece en la parte derecha (b) de la figura 5.6. Geométricamente representa una rotación del gráfico de dispersión, en la que el significado de los nuevos ejes es:
- \(A=\displaystyle \frac{1}{2}(\log_2 (R*G))\): El logaritmo de la intensidad media de los dos canales,
- \(M=\log_2 \displaystyle \frac{R}{G}\): El logaritmo de la expresión relativa entre ambos canales (normalmente conocido como “log–ratio”).
<- RG$R[,"DMSCCl319"]
R <- RG$G[,"DMSCCl319"]
G <- log(R)
logR <- log(G)
logG <- logR-logG
M <- 0.5*(logR+logG)
A <- par(mfrow=c(1,2))
optplot(logR~ logG, main= "Scatterplot")
abline(lm(logR~logG), col="yellow")
plot(M~A, main= "MA-Plot")
abline(h=0, col="yellow")
par(opt)
Si en vez de explorar un array mas moderno como el de la figura anterior lo hacemos con uno de los primeros que se utilizaron, el dataset swirl
vemos que el resultado es menos simetrico ,o que sugiere la necesidad de normalizacion.
Este dataset se encuentra disponible en el paquete swirl
pero también aparece como ejemplo en el paquete marray
. En la viñeta de dicho paquete se realiza una exploración completa de este dataset. Aquí nos limitaremos a reproducir los dos gráficos anteriores para ilustrar un caso en que sí que se pone de manifiesto la necesidad de normalizar los datos.
La clase del objeto “swirl” es distinta a la del objeto “RG” por lo que elk acceso es ligeramente distinto
library(marray)
data(swirl)
<- swirl@maRf[,1]
R <- swirl@maGf[,1]
G <- log(R)
logR <- log(G)
logG <- logR-logG
M <- 0.5*(logR+logG)
A <- par(mfrow=c(1,2))
optplot(logR~ logG, main= "Scatterplot")
abline(lm(logR~logG), col="yellow")
plot(M~A, main= "MA-Plot")
abline(h=0, col="yellow")
par(opt)
Aunque el diagrama de dispersión no lo muestra claramente, en el diagrama “MA” sí que se puede ver como los puntos estan claramente por debajo de lo que sería el eje de simetría horizontal
5.2.2.2 Imagen del array
La imagen del chip (véase la figura ??, izquierda) ofrece una visión rápida de la calidad del array, proporcionando información acerca del balance del color, la uniformidad en la hibridación y en los spots, de si el background es mayor del normal y dela existencia de artefactos como el polvo o pequeñas marcas (rasguños).
5.2.2.3 Histogramas de señales y de la relación señal–ruído
Estos gráficos (véase la figura 5.8, derecha) son útiles para detectar posibles anormalidades o un background excesivamente alto .
5.2.2.4 Boxplots
Un gráfico muy utilizado es el diagrama de cajas o “boxplot” múltiple con una caja por cada chip. Del alineamiento (o falta de él) y la semejanza (o disparidad) entre las cajas, se deduce si hace falta, o no, normalizar entre arrays.
En el caso de arrays de dos colores pueden utilizarse diagrams de cajas “dentro de arrays” (entre distintos sectores del mismo chip) y “entre arrays.”
5.2.3 Gráficos de diagnóstico para microarrays de un color
5.2.3.1 Imagen del chip
Los arrays de affymetrix contienen millones de sondas por lo que no pueden examinarse a simple vista. A pesar de ello hay diversas formas de obtener una imagen que, en caso de presentar irregularidades pueden indicar algún tipo de problemas como burbujas, arañazos, etc. La figura @ref(c06plotAffy} muestra algunas de las imágenes
La imagen del array de Affymetrix sólo es útil para evidenciar grandes problemas como burbujas, arañazos, etc. En este ejemplo los dos arrays de la izquierda se considerarían aceptables y los de la derecha defectuosos.
5.2.3.2 Gráfico “M-A”
En los chips de dos colores el MA–plot se utliza para comparar los dos canales en cada array (rojo y verde). En cambio, en los chips de Affymetrics, en que sólo hay un canal en cada array, la única forma de definir M (el log ratio) es a partir de la comparación entre pares de de valores, ya sea los arrays dos a dos o bien cada array respecto un valor de referencia que puede ser la mediana, punto a punto, de todos los arrays (véase por ejemplo la figura 5.10).
\(M=\log_2(I_1) - \log_2(I_2)\): log ratio \(A=\displaystyle \frac{1}{2}(\log_2 (I_1)+\log_2(I_2))\): log de intensidades Donde \(I_1\) es la intensidad del array de estudio, e \(I_2\) es la intensidad media de arrays. Por lo general, se espera que la distribución en el gráfico se concentre a lo largo del eje M = 0.
5.2.3.3 Modelos de bajo nivel (“Probe-Level-Models” o PLM)
Los modelos de bajo nivel (“Probe-Level-Models” o PLM) ajustan a los valores de intensidad –a nivel de sondas, no de valores totalizados de gen– un modelo explicativo. Los valores estimados por este modelo se comparan con los valores reales y se obtienen los errores o “residuos” del ajuste. El análisis de dichos residuos procede de forma similar a lo que se realiza al analizar un modelo de regresión: Si los errores no presentan ningún patrón especial supondremos que el modelo se ajusta relativamente bien. Si, en cambio, observamos desviaciones de esta presunta aleatoriedad querrá decir que el modelo no explica bien las observaciones, lo cual se atribuirá a la existencia de algún problema con los datos.
Con los valores ajustados del modelo se calculan dos medidas:
- La expresión relativa en escala logarítmica ” Relative Log Expression” (RLE) es una medida estandarizada de la expresión. No es de gran utilidad pero debería presentar una distribución similar en todos los arrays.
- El error no estandarizado y normalizado o “NUSE” es el más informativo ya que representa la distribucián de los residuos a la que hacíamos referencia más arriba. Si un array es problemático la caja correspondiente en el boxplot aparece desplazada hacia arriba o abajo de las demás.
stopifnot(require(affyPLM))
<- fitPLM(allAffyRaw) Pset
La figura 5.11 muestra los gráficos RLE y NUSE para el conjunto de datos estrogen. En ambos gráficos puede verse como el array defectuoso “” queda claramente diferenciado del resto.
<- par(mfrow=c(2,1))
optRLE(Pset, main = "Relative Log Expression (RLE)",
names=rownames(pData(allAffyRaw)), las=2, cex.axis=0.6)
NUSE(Pset, main = "Normalized Unscaled Standard Errors (NUSE)",
names=rownames(pData(allAffyRaw)), las=2, cex.axis=0.6)
par(opt)
5.3 Normalización de arrays de dos colores
La palabra normalización describe las técnicas utilizadas para transformar adecuadamente los datos antes de que sean analizados. El objetivo es corregir diferencias sistemáticas entre muestras, en la misma o entre imágenes, lo que no representa una verdadera variación entre las muestras biológicas.
Estas diferencias sistemáticas pueden deberse, entre otras, a:
- Cambios en la tinción que modifican la intensidad del spot.
- La ubicación en el array que puede afectar el proceso de lectura.
- Un problem en la placa de origen.
- La existencia de diferencias en la calidad de la impresión: pueden presentarse variaciones en los pins y el tiempo de impresión
- Camio en los parámetros de la digitalización (escaneo).
A veces puede ser difícil detectar estos problemas , aunque existen algunas formas de saber si es necesaria realizar una normalización. Aqui destacamos dos posibilidades:
- Realizar una auto-hibridación. Si hibridamos una muestra con ella misma, las intensidades deberían ser las mismas en ambos canales. Cualquier desviación de esta igualdad, significa que hay un sesgo sistemático.
- Detectar artefactos espaciales en la imagen o en la tinción de los gráficos de diagnóstico
5.3.1 Normalización global
Este método esta basado en un ajuste global, es decir en modificar todos los valores una cantidad c, estimada de acuerdo a algún criterio. \[\begin{equation} \log_2 R/G \rightarrow \log_2 R/G-c=\log_2 R/(Gk) \end{equation}\] opciones para \(k\) o \(c= \log_2k\) son
\(c\)= mediana o media de log ratio para un conjunto concreto de genes o genes control o genes housekeeping.
La intensidad total de la normalización
\(k=\sum R_i/\sum G_i\)
5.3.2 Normalización dependiente de la intensidad
En este caso se realiza una modificación específica para cada valor. Esta modificación se obtiene como una función de la intensidad total del gen (\(c=c(A)\)). \[\begin{equation} \log_2 R/G \rightarrow \log_2 R/G-c(A)=\log_2 R/(Gk(A)) \end{equation}\]
Una posible estimación de esta función puede hacerse utilizando la función lowess (LOcally WEighted Scatterplot Smoothing).
5.4 Resumen y normalización de microarrays de Affymetrix
En los arrays de Affymetrix, como en todos los tipos de microarrays, tras escanear la imagen se obtiene una serie de valores de intensidad de cada elemento del chip. En el caso de estos arrays sabemos que cada valor no corresponde a la expresión de un gen:
- Hay múltiples valores (sondas o “probes”) por cada gen, que originan un probeset.
- Cada grupo de sondas consiste en múltiples sondas
Originalmente, en los primeros arrays de Affymetrix se introdujo un sistema para estimar el “ruído” y diferenciarlo de la “señal.” Éste consistía en que cada sonda era en ralidad un par de sondas:
- Una que coincidía exactamente con el fragmento del gen al que corresponde la sonda y que se denominaba “perfect match.”
- Otra, que coincidía con la anterior salvo por el valor central que se sustituía por el nucleótido complementario. Esta segunda sonda se denominaba un “mismatch.”
La idea original consistía en utilizar estos “mismatches’ para tener una medida de hibridación no específica pero en las versiones más recientes se han abandonado, por lo que también se ha abandonado la utilización de los métodos”MAS-4” y “MAS-5” que utilizaban estos valores.
El proceso que convierte las señales individuales en valores de expresión normalizados para cada gen consta de tres etapas:
- Corrección del ruido de fondo o “background”
- Normalización para hacer los valores comparables
- “Sumarización”(Resumen) o concentración de los valores de cada grupo de sondas en un único valor de expresión absoluto normalizado para cada gen.
A menudo los tres pasos se denominan genérica -y erróneamente- “normalización.”
A diferencia de los chips de ADNc, aquí las medidas de expresión son absolutas (no se compara una condición contra otra) dado que cada chip se hibrida con un única muestra.
Hay muchos métodos para estimar la expresión (más de 30 publicados). Cada método contempla de forma explícita o implícita las tres formas de preprocesado: corrección del fondo, normalización y resumen.
Aquí nos centraremos únicamente en el método más popular que continua siendo utilizado trece años después de su introducción: El “Robust Multi-chip Average” o RMA
5.4.1 El método RMA (Robust Multi-Array Average)
Para compensar algunas deficiencias de los primeros métodos de resumen y normalización de arrays de Affymetrix, Irizarry y sus colegas introdujeron en 2003 (~) un método basado en la modelización de las intensidades de las sondas que, en vez de basarse en las distintas sondas de un gen dentro de un mismo array se basa en los distintos valores de la misma sonda entre todos los arrays disponibles,
Esquemáticamente los pasos que realiza este método son:
- Ajusta el ruído de fondo (background) basándose solo en los valores PM y utilizando un modelo estadístico complejo en el que combina la modelización de la señal mediante una distribución exponencial con la del ruído mediante una distribución normal.
- Toma logaritmos base 2 de cada intensidad ajustada por el background.
- Realiza una normalización por cuantiles de los valores del paso 2 consistente en substituir cada valor individual por el que tendría la misma posición en la distribución empírica estimada sobre todas las muestras, es decir los promedios de las distribuciones de los valores ordenados de cada array (véase figura )
- Estima las intensidades de cada gen separadamente para cada conjunto de sondas. Para ello realiza una técnica similar a una regresión robusta denominada “pulido de medianas” (median polish) sobre una matriz de datos que tiene los arrays en filas y los grupos de sondas en columnas.
Como resultado final de todos los pasos anteriores se obtiene la matriz con los datos sumarizados y normalizados. A pesar de no estar exento de críticas como la que afirma que este procedimiento “compacta” los valores reduciendo su variabilidad natural, este método se ha convertido en el estándar “de facto” actualmente por muchos usuarios de Bioconductor.
5.5 Filtraje no específico
El objetivo del filtraje es eliminar aquellos genes o aquellos spots cuyas imágenes o señales sean erróneas por diferentes motivos, con el fin de reducir el ruído de fondo.
Los principales procesos de filtrado son:
- Eliminación de los spots marcados como erróneos mediante flags y que son debidos a problemas en la hibridación o en el escaneo.
- Eliminación de spots o genes con señales muy bajas debido a problemas en el spotting o a que no ha habido hibridación en ese spot (Filtraje por señal).
- Eliminación de genes que no presenten una variación significativa en su señal, entre distintas condiciones experimentales (Filtraje por variabilidad).
Ante la duda, se suele aconsejar ser conservador y reducir la operación de filtraje al mínimo, puesto que cuesta mucho convencerse, y convencer a los investigadores, de que el beneficio del filtraje compensa la posibilidad de eliminar, por error, un gen o una señal interesantes.