---
title: "DNA Microarrays"
author: "Guillermo Ayala Gallego"  
date: "`r Sys.Date()`"  
address: Departamento de Estadística e I.O. Universidad de Valencia  
format: 
    revealjs:  
      scrollable: true
      echo: true
      embed-resources: true	  
    pdf:
       toc: true
       number-sections: false
       colorlinks: true	  
---


# Microarrays 


## DNA Microarrays

![](../figures/affymetrix_microarray.jpg)

## Un ejemplo

- [gse21779](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE21779)

## Affymetrix Genechip

```{r}
#| label: dnamicroarray1
#| echo: false
dirTamiData = "/home/gag/ownCloud/alltami/tami-data/"
```


```{r}
#| label: dnamicroarray2
#| echo: false
load(paste0(dirTamiData,"gse21779raw.rda"))
```
- ¿Qué tipo de objeto tenemos? ¿Cuál es su clase?

```{r  }
#| label: dnamicroarray3
class(gse21779raw)
```

- ¿Qué anotación tiene?

```{r }
#| label: dnamicroarray4
pacman::p_load(Biobase,affy)
annotation(gse21779raw)
```

- ¿Cuántas sondas y muestras?

```{r }
#| label: dnamicroarray5
dim(exprs(gse21779raw))
```

## La primera muestra

![gse21779_1](../figures/imagegse21779_1.png)


## AffyID

```{r}
#| label: dnamicroarray6
#| warning: false
head(probeNames(gse21779raw))
(ID = probeNames(gse21779raw)[400])
sum(probeNames(gse21779raw) == ID)
```

- ¿Qué tamaños tienen los grupos de sondas?

```{r }
#| label: dnamicroarray7
counts = table(probeNames(gse21779raw))
table(counts)
```

# Variabilidad intra y entre arrays

## PM y MM: código  

```{r }
#| label: dnamicroarray8
indexProbes(gse21779raw,"both",ID)[[1]]
(posiciones = indexProbes(gse21779raw,"pm",ID)[[1]])
indexProbes(gse21779raw,"mm",ID)[[1]]
```  

## PM y MM: sondas   

![](../figures/imagegse21779_1_con_sondas.png)

## Intra  

```{r }
#| label: dnamicroarray9
library(ggplot2)
pmID = probes(gse21779raw[,1],"pm",ID)
df1 = data.frame(sondas = 1:11,intensidad = pmID[,"GSM542488.CEL.gz"])
ggplot(df1,aes(x=sondas,y=intensidad))+geom_line()
```

## Entre

```{r }
#| label: dnamicroarray10
pm0 = probes(gse21779raw,"pm",ID)
pm0 = data.frame(pm0,sondas=1:11)
pm1 = reshape::melt(pm0,id="sondas")
ggplot2::ggplot(data=pm1,aes(x=sondas,y=value,colour=variable))+
    geom_line()+ylab("Intensidad")

```

## Comparando PM y MM


```{r }
#| label: dnamicroarray11
library(reshape)
pmm = data.frame(sondas=1:11,pm = probes(gse21779raw[,1],"pm",ID),
    mm = probes(gse21779raw[,1],"mm",ID)) 
df2 = reshape::melt(pmm,id="sondas")
levels(df2[,"variable"]) = c("pm","mm")
ggplot(df2,aes(x=sondas,y=value,colour=variable,linetype=variable))+geom_line()
```

# Control de calidad  

## MA plot

![](../figures/gse21779f7.png)

![](../figures/gse21779f8.png)


## Estimadores de densidad

![](../figures/gse21779global_a.png)


![](../figures/gse21779global_b.png)


## Diagramas de cajas  

![](../figures/gse21779global_c.png)



# MAS5


## ¿Qué tenemos?  

- En cada localización tendremos sondas que son oligonucleótidos específicos 
  de un gen.
- De hecho son pares de sonda: una con la secuencia deseada
  y otra en la que se modifica la base en la posicion 13.
- La intensidad (o expresión) observada con el oligonucleótido correcto le 
llamamos PM (perfect matching).
- La expresión observado sobre el modificado le llamamos MM (miss matching).
- Una cierta proporción de localizaciones muestran un valor MM mayor que 
el valor PM.
- <http://media.affymetrix.com/support/technical/whitepapers/sadd_whitepaper.pdf>


## Corrección de fondo  

- Se divide el array en K zonas rectangulares: $Z_k$ con $k=1,\ldots,K$.
- En cada $Z_k$ se ordenan las intensidades observadas en las distintas 
    celdas y se considera el 2\% con intensidades menores. 
- Calculamos su media  y desviación estándar muestrales: $b(Z_k)$ y $s(Z_k)$.
- Sea la celda localizada en $(x,y)$. 
- $d_k(x,y)$=  la distancia desde $(x,y)$ al centro de la región $Z_k$.
-  Consideramos la siguiente cantidad
  $$
  w_k(x,y) = \frac{1}{d^2_k(x,y) + s_0}
  $$
  donde $s_0$ nos evita que el cociente se pueda anular y por defecto 
  $s_0 = 100$.
-  Para la celda localizada en $(x,y)$ consideramos
  $$
  b(x,y) = \frac{1}{\sum_{k=1}^K w_k(x,y)} \sum_{k=1}^K w_k(x,y) b(Z_k)
  $$

## Ruido local

- Vamos a restar a las intensidades una medida de ruido
  de fondo local. 
- El problema  es que nos encontremos con valores negativos. 
- Este valor de ruido en la celda localizada en $(x,y)$ se calcula como
  $$
  n(x,y) = \frac{1}{\sum_{k=1}^K w_k(x,y)} \sum_{k=1}^K w_k(x,y) s(Z_k).
  $$
- La idea ahora es quitar a las intensidades originales el fondo
  ($b(x,y)$) que hemos calculado pero asegurándonos que no nos surja
  valores negativos. 
- La intensidad original se denota como $I(x,y)$ en la localización 
  $(x,y)$. El valor ajustado de la intensidad viene dado por
  $$
  A(x,y) = \max\{I(x,y) - b(x,y), f_0 n(x,y) \}
  $$
- El parámetro $f_0$ se suele fijar en $f_0 =0.5$. 

## Cálculo del valor de expresión

- Por el procemiento descrito anteriormente podemos ajustar los valores 
de PM y de MM. 
- En lo que sigue se supone que hemos ajustado estos valores.
- Lo esperable y correcto es que en una celda tengamos un valor de PM mayor 
que el correspondiente valor MM.
- Esto no es así en muchas ocasiones.
- Se define un valor IM (Ideal Mismatch).

## Ideal Mismatch

- Denotamos: $(PM_{ij},MM_{ij})$ el valor de PM y de MM de la j-ésima 
sonda perteneciente al i-ésimo conjunto de sondas.
- Se calcula un valor que llaman background específico para cada conjunto
de sondas. 
$$
B_i = T_{bi} \{\log_2(PM_{ij}) -  \log_2(MM_{ij}): j= 1, \ldots,n_i \},
$$
donde $T_{bi}$ denota el algoritmo bipesado en un solo paso. 

## 
- Si el valor de 
$$
IM_{ij}  = \left\{ \begin{array}{cc}
    MM_{ij} & \mbox{si $MM_{ij}< PM_{ij}$} \\
    \frac{PM_{ij}}{2^{ B_i}} &\mbox{si $MM_{ij} \geq PM_{ij}$ y 
      $ B_i > \alpha $} \\
    \frac{PM_{ij}}{2^{\frac{\alpha}{1 + (\alpha -  B_i)/\beta}}} &
\mbox{si $MM_{ij} \geq PM_{ij}$ y $ B_i \leq \alpha $} 
  \end{array} \right.
$$
- Por defecto: $\alpha = 0.03$ y $\beta = 10$.


## Valor resumen

- Calculamos 
$$
V_{ij} = \max\{PM_{ij} - IM_{ij},d\}
$$
siendo por defecto $d = 2^{-20}$.
- Transformamos 
$$
{\mathbb V}_{ij} = \log_2 (V_{ij})
$$
con $j = 1, \ldots, n_i$. 
- Y volvemos a aplicar un estimador biponderado en un solo paso. 
$$
SignalLogValue_i = T_{bi}({\mathbb V}_{i1}, \ldots, {\mathbb V}_{in_i}). 
$$

## 

- Fijamos una señal objetivo, Sc: por defecto, $Sc = 500$.
- Se calcula el siguiente factor de escala para el grupo de sondas i.
$$
sf_i = \frac{Sc}{MediaAjustada\{2^{SignalLogValue_i},0.02,0.98\}}
$$
- **MediaAjustada** indica la media ajustada 
(media de los valores quitando una cierta proporción de los más grandes
y de los más pequeños) en donde se quita el 2\% de los más grandes y 
el 2\% de los más pequeños. 
- El valor final para el conjunto de sondas i es 
$$
ReportedValue(i) = sf_i \times 2^{SignalLogValue_i}
$$
- Se incorpora una posibilidad de normalización que no indicamos.


##  Algoritmo biponderado de Tukey

- Tenemos unos datos $x_1,\ldots, x_n$.
- Calculamos la mediana $m_x$ de los datos $x_1,\ldots, x_n$.
- Determinamos las distancias $d_i = |x_i - m_x|$. 
- Calculamos la mediana $m_d$ de $d_i$ con $i = 1, \ldots, n$
 (MAD, median absolute deviation).
- Consideramos
  $$
  u_i = \frac{x_i - m_x}{c m_d + \epsilon},
  $$
  con $i = 1, \ldots, n$. Los valores por defecto son $c= 5$ y $e =0.0001$
  (con el valor de e evitamos la división por cero).

## 
- Consideremos una función bicuadrada (bisquare) definida como
    $$
w(u)  = \left\{ \begin{array}{cl}
 (1- u^2)^2 &\mbox{si  $|u| \leq 1$} \\
  0 &\mbox{si $|u| >1$}
       \end{array} \right.$$
       
- Se define el estimador en un solo paso como
 $$
  T_{bi}(x_1,\ldots,x_n) = \frac{\sum_{i=1}^n w(u_i) x_i}{\sum_{i=1}^n w(u_i)}. 
  $$
       
- Si sustituimos el valor $m_x$ por $T_b$ podríamos aplicar
iterativamente el procedimiento hasta que se estabilice.



# Robust Multichip Average (RMA)


## RMA 

- En este método trabajamos con todos los arrays simultáneamente.
- A diferencia del método anterior donde se procesaba cada array 
independientemente. 
- El procedimiento solamente utiliza los valores PM y no utiliza para nada
los valores MM. 

## Corrección de fondo 

- El valor observado aleatorio $S$ sería
$$
S = X + Y
$$
con $X \sim Exp(\alpha)$ e $Y \sim N(\mu,\sigma^2)$.

- Se tiene
$$
E(X | S = s) = a + b \frac{\phi(\frac{a}{b})}{\Phi(\frac{a}{b})}
$$
siendo $a = s - \mu - \sigma^2 \alpha$ y $b= \sigma$ mientras que   
$\phi$ y $\Phi$ denotan las funciones de densidad y de distribución
de una normal estándar (con media 0 y varianza 1). 

- Los parámetros del modelo ($\alpha, \mu, \sigma^2$) se estiman mediante
un procedimiento no paramétrico.

## Normalización de cuantiles

- Tomamos n vectores de longitud N y construimos la matriz $X$  
  $N \times n$ que los tiene por vectores columna.
- Ordenamos cada columna de $X$ separadamente y tenemos $X_s$.
- Calculamos la media por filas de la matrix $X_s$ y creamos
  $X'_s$ una matriz con la misma dimensión que $X$, tal que en la fila j
  tenemos la media de la fila de $X_s$ repetida n veces.
- Obtenemos $X_t$ en donde cada columna de $X'_s$ recupera el orden 
  original. 

## Resumen

- En cada array tendremos $N$ sondas.
- Suponemos que tenemos $n$ arrays.
- La matriz de expresión a nivel de sonda será: 
$$ {\mathbf x} = [x_{ij}]_{i=1,\ldots,N; j =1,\ldots,n}$$
-   La sonda $i$-ésima en el array $j$-ésimo tiene una expresión 
o intensidad $x_{ij}$. 
- Denotamos el grupo $k$-ésimo de sondas con $S_k$.
- $S_k \subset \{1,\ldots,N\}$
- $S_k = \{i_1,\ldots,i_{|S_k|}\}$
- Denotamos
$${\mathbf y} = [y_{ij}]_{i=1,\ldots,|S_k|,j=1,\ldots,n} = 
[x_{ij}]_{i \in S_k,j =1,\ldots,n}.$$ 

## Median polish
- Consideramos la matriz para **cada grupo de sondas**. 
$$
 \begin{array}{ccc|c}
  \delta_{1,1} &  \cdots & \delta_{1,n} & a_1 \\
   \vdots  &  \ddots & \vdots           & \vdots \\
  \delta_{|S_k|,1} &   \cdots & \delta_{|S_k|,n} & a_{|S_k|}\\
\hline
 b_1 &   \cdots & b_{n} & m
 \end{array}
$$

- Fijamos $\delta_{ij} = \log_2(y_{ij})$, $a_i = b_i = 0 \forall i,j$.
- Calculamos la mediana para cada fila a lo largo de las columnas
ignorando la última columna (separada por la línea).
- Restamos a cada elemento de la fila la mediana correspondiente.
Esta mediana se suma a la última columna (formada por 
$a_1, \ldots, a_{|S_k|},m$.
- Calculamos la mediana para cada columna de los valores correspondientes
a las distintas filas sin considerar la última fila.
- Restamos la mediana calculada en el paso anterior a cada elemento
de la columna exceptuando la última fila. A la última fila sumamos
las medianas calculadas. 
- El proceso descrito en los cuatro pasos anteriores continua iterando
sucesivamente entre filas y columnas hasta que los cambios que se producen
son nulos o muy pequeños. 

## Estimaciones
Una vez finalizado el proceso iterativo tendremos:  

- $\hat{\mu} = m$
- $\hat{\theta}_j = b_j$
- $\hat{\alpha}_i = a_i$ 
- $\hat{\epsilon}_{ij} = \delta_{ij}$
- ¿Y cómo estimamos la expresión de un grupo de sondas en una array?
 $$
\hat{\mu} + \hat{\theta}_j
$$

- **Y esto lo repetimos para cada grupo de sondas.** 

# affy

## MAS5 y RMA
```{r}
#| label: dnamicroarray12
library(affy)
```

- El método **MAS5** 

```{r ,eval = FALSE}
#| label: dnamicroarray13
#| eval: false
gse21779_mas5 = mas5(gse21779raw)
```

- El método **RMA**

```{r}
#| label: dnamicroarray14
#| eval: false
gse21779_rma = rma(gse21779raw)
```


```{r}
#| label: dnamicroarray15
#| echo: false
load(paste0(dirTamiData,"gse21779_mas5.rda"))
load(paste0(dirTamiData,"gse21779_rma.rda"))
```

## Densidades: MAS5

```{r}
#| label: dnamicroarray16
library(geneplotter)
geneplotter::multidensity(exprs(gse21779_mas5))
```


## Densidades: RMA

```{r }
#| label: dnamicroarray17
geneplotter::multidensity(exprs(gse21779_rma))
```


## Diagramas de cajas: MAS5

```{r }
#| label: dnamicroarray18
graphics::boxplot.matrix(exprs(gse21779_mas5))
```

## Diagramas de cajas: RMA

```{r}
#| label: dnamicroarray19
graphics::boxplot.matrix(exprs(gse21779_rma))
```
 
## `ExpressionSet`

```{r }
#| label: dnamicroarray20
#| message: false  
pacman::p_load(ALL)
```

```{r}
#| label: dnamicroarray21
data(ALL)
```
La clase del objeto la podemos ver con 
```{r}
#| label: dnamicroarray21b
class(ALL)
```
El método `print` nos da un breve resumen de la clase.
```{r}
#| label: dnamicroarray22b
#| eval: false
ALL
```
La matriz de expresión la obtenemos con el accessor  `Biobase::exprs`.
```{r}
#| label: dnamicroarray22
#| echo: false
head(exprs(ALL))
```

Los datos fenotípicos los tendremos con `Biobase::pData`.

```{r }
#| label: dnamicroarray23
#| eval: false
pData(ALL)
```

¿Cuál fue chip que se utilizó para obtener estos datos?

```{r }
#| label: dnamicroarray24
annotation(ALL)
```

Los identificadores de las filas los tenemos con 

```{r}
#| label: dnamicroarray25
#| eval: false  
featureNames(ALL)
```


# Construyendo un ExpressionSet (Neve 2006)

## Lo necesario

- Tenemos los datos de expresión en un fichero.
- En **otro** fichero tenemos información sobre los genes
  que tenemos en las filas de la matriz de expresión.
- Tenemos información sobre las distintas muestras que componen toda la 
experimentación.
- Y finalmente sabemos cosas sobre toda la experimentación realizada:
autores, publicación principal que se derivó, objetivo de la experimentación.
- Y todo por separado, en distintos ficheros.
- ¿Cómo construimos el **ExpressionSet**?


##  Lectura datos 

- Cargamos Biobase.

```{r}
#| label: dnamicroarray26
library(Biobase) 
```
- Leemos datos de expresión.

```{r}
#| label: dnamicroarray27
library(tamidata)
finput = system.file("extdata","n06_expr.txt",package="tamidata")
exprs0 = read.table(file = finput,header = TRUE,sep="\t")
```
- Tiene que ser una matriz, `matrix`.

```{r}
#| label: dnamicroarray28
exprs0 = as.matrix(exprs0[,-1])
```
- Como matriz que es tendrá `rownames` y `colnames`. 
```{r}
#| label: dnamicroarray29
colnames(exprs0)
rownames(exprs0)
```
- Como vemos está vacío el elemento `rownames`.

## ...

- Leemos variables fenotípicas.

```{r }
#| label: dnamicroarray30
finput = system.file("extdata","n06_sample_attributes.txt",package="tamidata")
pd0 = read.table(file = finput,header = TRUE,sep="\t")
```
- Podemos ver el número de muestras (que ha de coincidir
con el número de columnas de la matriz de expresión) y de covariables.

```{r }
#| label: dnamicroarray31
dim(pd0)
dim(exprs0)
```
- Es **necesario** que los nombres de las filas de **pd0** coincidan
con los nombres de las columnas de **exprs0**.
```{r }
#| label: dnamicroarray32
rownames(pd0) = colnames(exprs0)
```

##  
- Etiquetamos los metadatos: la columna `labelDescription` debe de estar
presente.

```{r }
#| label: dnamicroarray33
metadatos = data.frame(labelDescription = c("Identificador",
"Número de muestra", "Tipo","En CGH","ALL","Método CDH1","SFRP1",
"HER2","KIT","Presencia de p53","BRCA","Orden","Adherencia",
"per3"),
row.names=colnames(pd0))
```

- Reunimos etiquetas de metadatos y metadatos en un **AnnotatedDataFrame**.
```{r }
#| label: dnamicroarray34
datosfenotipo = new("AnnotatedDataFrame", data = pd0,
    varMetadata = metadatos)
```

##  

- Una vez tenemos este `AnnotatedDataFrame` podemos ver los nombres de las
muestras.

```{r }
#| label: dnamicroarray35
sampleNames(datosfenotipo)
```
O también ver los valores de las covariables.
```{r }
#| label: dnamicroarray36
head(pData(datosfenotipo))
```

##  

- Los experimentos los obtienen investigadores y hay que incluir esta 
información en los datos. 

```{r }
#| label: dnamicroarray37
datosexperimento = new('MIAME',name='Neve et alt.',lab='Varios',
                       contact ='rmneve@lbl.gov',
                       title = 'A collection of breast cancer cell lines 
                       for the study of
functionally distinct cancer subtypes',
abstract = 'An example ExpressionSet',
url = '10.1016/j.ccr.2006.10.008',
other = list(notes = 'Creado a partir  de ficheros de texto'))
```
El formato MIAME (Minimum information about a microarray experiment)
puede consultarse en  <http://fged.org/projects/miame/>

##  Y montamos las piezas

```{r }
#| label: dnamicroarray38
jcliment  = new("ExpressionSet",exprs=exprs0,phenoData = datosfenotipo,
  experimentData = datosexperimento,annotation = "hgu95av2")
```

Guardamos el `ExpressionSet`.

```{r}
#| label: dnamicroarray39
#| eval: false  
save(jcliment,file="jcliment.rda")
```
```{r}
#| label: dnamicroarray39b
#| echo: false
save(jcliment,file=paste0(dirTamiData,"jcliment.rda"))
```

# ArrayExpress

## Otra vez Neve (2006)

- Hemos visto cómo construir un `ExpressionSet` cuando tenemos los datos de
expresión, los fenotípicos y los de anotación.
- Los autores han subido toda la información del experimento a 
<https://www.ebi.ac.uk/arrayexpress/>.

```{r}
#| label: dnamicroarray40
#| eval: false
library(ArrayExpress)
rawset = ArrayExpress("E-TABM-157")
```
Nos devuelve un objeto `affy::AffyBatch`. 

# GSE21779 (de dos formas)

## gse21779:  Gene Expression Omnibus (GEO)  
- <http://www.ncbi.nlm.nih.gov/geo/>  
- Bajamos los datos GSE21779.  
- Los ficheros CEL los tenemos en  
<http://www.ncbi.nlm.nih.gov/sites/GDSbrowser?acc=GDS4128>

- Podemos usar **GEOquery**

```{r}
#| label: dnamicroarray41
#| message: false
library(GEOquery)
```
- Sabiendo el identificador podemos bajar los datos.

```{r}
#| label: dnamicroarray42
#| eval: false
gcel = getGEOSuppFiles("GSE21779")
```

```{r}
#| label: dnamicroarray43
#| eval: false
setwd("GSE21779/")
system("tar xvf GSE21779_RAW.tar")
library(affy)
gse21779raw = ReadAffy()
setwd("../")
system("rm -fr GSE21779")
gse21779rma = affy::rma(gse21779raw)
pData(gse21779rma)
gse = getGEO("GSE21779")
class(gse)
length(gse)
class(gse[[1]])
pData(gse[[1]])
rownames(pData(gse21779rma)) =
    sapply(rownames(pData(gse21779rma)),function(x)  unlist(strsplit(x,split=".CEL."))[1])
all(rownames(pData(gse21779rma)) == rownames(pData(gse[[1]])))
pData(gse21779rma) = pData(gse[[1]])
gse21779 = gse21779rma
```

## geod21779: ArrayExpress

```{r}
#| label: dnamicroarray44
#| eval: false
library(ArrayExpress)
geod21779raw = ArrayExpress("E-GEOD-21779")
geod21779 = rma(geod21779raw)
```

# Correspondencias múltiples  

## El problema  

- Necesitamos conocer la correspondencia entre sondas y genes.  
-  La opción más simple y primera a probar es utilizar un paquete de anotación. 
- <https://www.bioconductor.org/packages/release/data/annotation/>

```{r}  
#| label: dnamicroarray45
#| echo: false
dirTamiData = "/home/gag/ownCloud/alltami/tami-data/"
load(paste0(dirTamiData,"gse21779raw.rda"))
```  
- En Aula Virtual tenemos los datos. 
```{r}
#| label: dnamicroarray46
#| eval: false
load("gse21779raw.rda")
```
- Normalizamos los datos.  

```{r}  
#| label: dnamicroarray47
#| eval: false
library(Biobase)
dim(gse21779raw)
gse21779rma = affy::rma(gse21779raw)
save(gse21779rma,file=paste0(dirTamiData,"gse21779rma.rda"))
```

```{r}
#| label: dnamicroarray48
#| echo: false
load(paste0(dirTamiData,"gse21779rma.rda"))
```
- ¿Cuántas sondas tenemos después de la normalización?  
```{r}
#| label: dnamicroarray49
class(gse21779rma)
nrow(gse21779rma)
ncol(gse21779rma)
dim(gse21779rma)
annotation(gse21779rma)
```
- Establecemos las correspondencias con `Entrez`
con `AnnotationDbi::select()`.
```{r}
#| label: dnamicroarray50
pacman::p_load(hgu133plus2.db)
a = AnnotationDbi::select(hgu133plus2.db,
                          keys=featureNames(gse21779rma),
                          columns=c("ENTREZID"),
                          keytype="PROBEID")
dim(a)
```
- Vemos que tenemos más filas en `a` que en `gse21779`.
```{r}
#| label: dnamicroarray51
dim(a)
```
- Tenemos las correspondencias en forma de `data.frame`.  
```{r}
#| label: dnamicroarray52
class(a)
```
- ¿Cuántas y qué sondas tienen correspondencia con más de un gen?

```{r}
#| label: dnamicroarray53
b = which(table(a[,1]) > 1)
length(b)
sel = is.element(a[,1],names(b))
table(table(a[sel,1]))
```  
- Una posibilidad puede ser quedarnos con la primera aparición.  

```{r}
#| label: dnamicroarray54
c1 = match(unique(a[,1]),a[,1])
a1 = a[c1,]
fData(gse21779rma) = a1
```  
- Vamos a ver si cada identificador `ENTREZID` tiene una sola
correspondencia `PROBEID`.  
```{r}
#| label: dnamicroarray55
table(table(a1[,2]) >1)
```  
- Por ello hemos de quedarnos con una sola sonda por identificador 
`ENTREZID`.  
- Podemos repetir el código anterior para estos identificadores.  
```{r}
#| label: dnamicroarray56
c2 = match(unique(a1[,2]),a1[,2])
a2 = a1[c2,]
dim(a2)
gse21779 = gse21779rma[c2,]
```  

- Podemos comprobar que ahora sí hay una correspondencia 1-1.  
```{r}
#| label: dnamicroarray57
length(unique(a2[,1])) == length(unique(a2[,2]))  
```  
- Llevamos de nuestro atrevimiento queremos establecer la correspondencia
de los identificadores `PROBEID`  con `ENTREZID` y con `ENSEMBL`.  
```{r}
#| label: dnamicroarray58
a = AnnotationDbi::select(hgu133plus2.db,
                          keys=featureNames(gse21779rma),
                          columns=c("ENTREZID","ENSEMBL"),
                          keytype="PROBEID")
head(a)
c1 = match(unique(a[,1]),a[,1])
a1 = a[c1,]
c2 = match(unique(a1[,2]),a1[,2])
a2 = a1[c2,]
```
- Sin embargo, hemos de tener en cuenta que no tenemos identificadores
únicos `ENSEMBL`. 
```{r}
#| label: dnamicroarray59
table(table(a1[,3]) >1)
```  
- Si queremos una correspondencia 1-1 hemos de decidir entre qué par de 
identificadores lo hacemos.


# Un ejemplo Agilent   

## GSE104645 

```{r}
#| label: dnamicroarray60
#| eval: false
wd = getwd()
setwd(dirTamiData)
GEOquery::getGEOSuppFiles("GSE104645")
setwd("GSE104645")
system("tar xvf GSE104645_RAW.tar")
system("gzip -d *.gz")
GSE104645raw = limma::read.maimages(dir(".", "txt"), "agilent",
                                    green.only = TRUE,
                                    other.columns="gIsWellAboveBG")
save(GSE104645raw,file=paste0(dirTamiData,"GSE104645raw.rda"))
setwd(wd)
```  

```{r}
#| label: dnamicroarray61
#| echo: false  
load(paste0(dirTamiData,"GSE104645raw.rda"))
```  
Eliminamos correspondencias múltiples.
```{r}
#| label: dnamicroarray62
#| eval: false
pacman::p_load(hgug4112a.db)
a = AnnotationDbi::select(hgug4112a.db,
                   keys=GSE104645raw$genes[,"ProbeName"],
                   column=c("ENTREZID","ENSEMBL"),
                   keytype="PROBEID")
## Eliminamos sondas sin ENTREZID
a = a[!is.na(a[,"ENTREZID"]),] 
c1 = match(unique(a[,1]),a[,1])
a1 = a[c1,]
c2 = match(unique(a1[,2]),a1[,2])
a2 = a1[c2,]
a2 = na.omit(a2)
GSE104645raw2 =
    GSE104645raw[match(a2[,1],
                       GSE104645raw$genes[,"ProbeName"]),]
dim(GSE104645raw2)
rownames(GSE104645raw2) = a2[,1]
save(GSE104645raw2,file=paste0(dirTamiData,
                               "GSE104645raw2.rda"))
```  
Hacemos la corrección de fondo.
```{r}
#| label: dnamicroarray63
#| eval: false
GSE104645el = limma::backgroundCorrect(GSE104645raw2,
                                       method="normexp")
```  
Aplicamos una normalización de cuantiles.
```{r}
#| label: dnamicroarray64
#| eval: false
GSE104645el = limma::normalizeBetweenArrays(GSE104645el,
                                            method="quantile")
```  
Guardamos los datos.
```{r}
#| label: dnamicroarray65
#| eval: false
save(GSE104645el,file=paste0(dirTamiData,"GSE104645el.rda"))
```  

```{r}
#| label: dnamicroarray66
#| echo: false
load(paste0(dirTamiData,"GSE104645el.rda"))
```  
Contruimos el `Biobase::ExpressionSet`.
```{r}
#| label: dnamicroarray67
pd0 = data.frame(case = 1:ncol(GSE104645el))
rownames(pd0) = colnames(GSE104645el$E)
metadata = data.frame(labelDescription = c("case"),
                      row.names=colnames(pd0))
phenotypedata = new("AnnotatedDataFrame", data = pd0,
                    varMetadata = metadata)
experimentdata = new('MIAME',name="Okita A, Takahashi S,
 Ouchi K, Inoue M et al.",
 lab="Tohoku University Hospital",
 contact ="akira.okita.d8@tohoku.ac.jp",
 title = "Consensus molecular subtypes classification
 of colorectal cancer as a predictive factor for 
chemotherapeutic efficacy against metastatic colorectal 
cancer",
abstract = "Gene expression profile from 193
 formalin-fixed and paraffin embedded primary tumor 
samples.",
url = "https://www.ncbi.nlm.nih.gov/pubmed/29721154",
other = list(notes = "An example of Agilent microarray"))

GSE104645 = new("ExpressionSet",exprs=GSE104645el$E,
                phenoData = phenotypedata,
                experimentData = experimentdata,
                annotation = "hgug4112a.db")

```  
Añadimos el `ExpressionSet::fData` con los identificadores.
```{r}
#| label: dnamicroarray68
fData(GSE104645) = a2[,1:2]
```  

```{r}
#| label: dnamicroarray69
#| echo: false
save(GSE104645,file=paste0(dirTamiData,"GSE104645.rda"))
```  

