---
title: "Limma"
author: "Guillermo Ayala Gallego"  
date: "`r Sys.Date()`"  
address: Departamento de Estadística e I.O. Universidad de Valencia  
bibliography: /home/gag/ownCloud/BIBLIOGRAFIA2/Bibliografia.bib
format: 
   revealjs:  
      scrollable: true
      echo: true
      embed-resources: true  
   pdf:
       toc: true
       number-sections: false
       colorlinks: true
---


# El método limma  

## Lo previo  
- Propuesto en @GKSmyth04  
  
- La respuesta es de caracter continuo (DNA microarrays,
Methylation array, ...) 
- Respuesta aleatoria: 
  $$
  \mathbf Y_i = (Y_{i1},\ldots,Y_{in})^T. 
  $$ 
- Valores observados:
   $$
   \mathbf y_i = (y_{i1},\ldots,y_{in})^T. 
   $$  
- Para cada muestra nuestros predictores son las variables fenotípicas
observadas en dicha muestra:  $({\mathbf x}_1,\ldots,{\mathbf x}_n)$,
	con ${\mathbf x}_j \in \mathbb R^p$.
		- Denotamos: $E\mathbf Y_i = (E Y_{i1},\ldots,E Y_{in})^T$.
			- Asumimos:  
$$
E Y_{ij} = \beta_{i1} x_{j1}+ \ldots  + \beta_{ip}  x_{jp} = \mathbf x_j^T 
\mathbf \beta_i,
$$
siendo:
- ${\mathbf x}_j = (x_{j1},\ldots,x_{jp})^T$ es el vector de covariables para
		la $j$-ésima muestra,  
- ${\mathbf \beta}_i = (\beta_{i1},\ldots,\beta_{ip})^T$ el vector de 
coeficientes.  
- Sobre la matriz de varianzas de $\mathbf Y_i$
$$
 var(\mathbf Y_i) = \sigma^2_i \mathbf V_i,
$$
donde ${\mathbf V}_i$ es una matriz definida positiva conocida.  

## Contrastes  

- Estamos interesados en ciertos contrastes sobre el vector de coeficientes 
${\mathbf \beta}_i$ que corresponden al diseño utilizado en la experiencia.
- Un contraste tiene la forma general
$$
\mathbf a_i = \mathbf C^T \mathbf \beta_i.
$$
- Se quiere contrastar
$$
H_0:  a_{ij} = 0,
$$
$$
H_1:  a_{ij} \neq 0.
$$  

##  

- **Ajustamos un modelo lineal para cada gen** y obtenemos:
- Los estimadores $\hat{\mathbf \beta}_i$,
- el estimador $s^2_i$ de $\sigma^2_i$ y
- la matriz de covarianzas estimada
$$
  \widehat{var(\hat{\mathbf \beta}_i)} =  s^2_i \mathbf V_i
$$
siendo $\mathbf V_i$ una matriz definida positiva que no depende de $s^2_i$.

- Los estimadores de los contrastes son
$$
\hat{\mathbf a}_i = \mathbf C^T \hat{\mathbf \beta}_i,
$$
y de su matriz de covarianzas es 
$$
var(\hat{\mathbf a}_i) = \sigma^2_i \mathbf C' \mathbf V_i \mathbf C.
$$
La correspondiente matriz de covarianzas estimada es 
la obtenemos sustituyendo $\sigma_i^2$ por $s_i^2$.

## Hipótesis sobre los estimadores  

- No se asume necesariamente normalidad para ${\mathbf Y}_i$.
- Tampoco asumimos que el ajuste de los modelos lineales se hace
  mediante el procedimiento de los mínimos cuadrados. 
- **Sí** se asume que:
    + $\hat{\mathbf a}_i$ tiene una distribución aproximadamente normal
	con vector de medias $\mathbf a_i$ y matriz de covarianzas
	$\sigma^2_i \mathbf C^T \mathbf V_i \mathbf C$.
- También asumimos que $s^2_i$ siguen aproximadamente una distribución 
  ji-cuadrado escalada.

## Consecuencias  

- Si denotamos por $v_{ij}$ el $j$-ésimo elemento de la diagonal de
  $\mathbf C' \mathbf V_i \mathbf C$ entonces todas las hipótesis 
  indicadas previamente se traducen en los siguientes resultados.
- $\hat{a}_{ij} | a_{ij}, \sigma^2_i
\sim N(a_{ij},v_{ij} \sigma^2_i)$.
- $s^2_i | \sigma^2_i \sim \frac{\sigma^2_i}{d_i} \chi^2_{d_i}$
siendo $d_i$ los grados de libertad residuales del modelo ajustado para el 
$i$-ésimo gen. 
- Asumiendo estas hipótesis se tiene que el siguiente estadístico
$$
t_{ij} = \frac{\hat{\beta}_{ij}}{s_i \sqrt{v_{ij}}} 
$$ 
aproximadamente tiene una distribución t de Student con $d_i$ grados de 
libertad. 
- En lo que sigue se asumirá (y no tienen porqué ser cierto) que 
$\hat{\beta}_i$ y $s^2_i$ son independientes para los distintos
genes.
	
## Modelo jerárquico
- Tenemos muchos test simultáneos.
- Se trata de modelizar el comportamiento entre genes.
- Se opta por modelizar mediante distribuciones de 
 probabilidad sobre los parámetros.
- Es un modelo bayesiano.

## Distribuciones a priori sobre los parámetros
- Sobre $\sigma^2_i$ la siguiente distribución,
$$
\frac{1}{\sigma^2_i} \sim \frac{1}{d_0 s^2_0} \chi^2_{d_0}.
$$
- Para un $j$ dado, suponemos que $\beta_{ij}$ es no nula con una probabilidad
conocida
$$
P(\beta_{ij} \neq 0) = p_j.
$$
- Notemos que $p_j$ tiene el sentido de la proporción esperada de genes que se 
expresan diferencialmente.

- Para los coeficientes no nulos se asume
$$
\beta_{ij} \bigg | \sigma^2_i,\beta_{ij} \neq 0 \sim N(0, v_{0j}\sigma^2_i).
$$

## Distribuciones a posteriori 

Asumiendo el modelo jerárquico que acabamos de especificar se tiene que
la media de la distribución a posteriori de $1/\sigma^2_i$ condicionada
a  $s^2_i$ viene dada por 
$$
E \bigg [\frac{1}{\sigma^2_i} \bigg | s^2_i \bigg ] = \frac{1}{\tilde{s}^2_i},
$$
con 
$$
\tilde{s}^2_i = \frac{d_0 s^2_0 + d_i s^2_i}{d_0 + d_i}. 
$$

- Podemos definir el estadístico t moderado como
$$
\tilde{t}_{ij} = \frac{\hat{\beta}_{ij}}{\tilde{s}^2_i \sqrt{v_{ij}}}.
$$
- Se demuestra que los t-estadísticos moderados $\tilde{t}_{ij}$ y las
varianzas muestrales residuales $s^2_i$ se distribuyen independientemente.
- Bajo la hipótesis nula $H_0: \beta_{ij}=0$, el 
  t-estadístico moderado $\tilde{t}_{ij}$ sigue una distribución t de Student
con $d_i + d_0$ grados de libertad.
- Los grados de libertad que estamos
añadiendo $d_0$ expresan la ganancia de información que obtenemos de 
utilizar todos los genes siempre **asumiendo el modelo jerárquico**.
- Los valores $d_0$ y $s_0$ que se suponen conocidos en lo previo serán 
estimados a partir de los datos.
- Tenemos un método empírico bayesiano.

## Odds a posteriori

- Se consideran los **odds a posteriori** dados por 
$$
  O_{ij} = \frac{P(\beta_{ij} \neq 0 | 
  \tilde{t}_{ij,s^2_{i}})}{P(\beta_{ij} = 0 | \tilde{t}_{ij,s^2_{i}})}.
$$
- Estos odds prueban que tienen la siguiente expresión.
$$
O_{ij} = \frac{p_j}{1-p_j} \bigg ( \frac{v_{ij}}{v_{ij}+v_{0j}}\bigg )^{1/2} 
\bigg (
\frac{\tilde{t}^2_{ij} + d_0 + d_i}{\tilde{t}^2_{ij} 
\frac{v_{ij}}{v_{ij}+v_{oj}} + d_o + d_i}
\bigg )^{(1-d_0+d_i)/2}
$$
- Se propone trabajar con los **log-odds** dados por
$$
  B_{ij} = \ln O_{ij}.
$$

## Relación con método SAM

- Es similar en espíritu pero no es lo mismo.
- Si los grados de libertad verifican $d_o < +\infty$ y $d_i > 0$
entonces los t-estadísticos moderados tienen la siguiente expresión
$$
\tilde{t}_{ij} = \bigg ( \frac{d_0 + d_i}{d_i} \bigg )^{1/2}
\frac{\hat{\beta}_{ij}}{\sqrt{s^2_{*,i} v_{ij}}},
$$
siendo
$$
s^2_{*,i} = s_i^2 + \frac{d_0}{d_i} s^2_0. 
$$
- Cada varianza muestral es incrementada un valor positivo en principio
distinto para cada gen.
-  La idea de SAM consiste en considerar
$$
t_{*,ij} = \frac{\hat{\beta}_{ij}}{(s_i + a) \sqrt{v_{ij}}}
$$
donde $a$ es un valor positivo determinado por un método ad-hoc. En este
método lo que se incrementa es la desviación estándar y no la varianza. 
- El t-estadístico moderado incrementa la varianza.
- Además el incremento puede ser distinto para cada gen. 

# gse25171  

## Datos  

```{r}
#| label: limma1
pacman::p_load(Biobase)
data(gse25171,package="tamidata2")
dim(gse25171)
pData(gse25171)
```

- ¿Cómo pueden influir el tiempo de observación
de la muestra `time` y la presencia o no de fósfatos en la expresión
del gen `Pi` en la expresión de cada una de las sondas que tenemos
en el microarray?

## `time`  

```{r}
#| label: limma2
design = model.matrix(~ pData(gse25171)[,"time"])
head(design)
```
- Cambiamos los nombres de la columnas por razones estéticas.
```{r}
#| label: limma3
colnames(design) = c("intercept","time")
```
- Ajustamos todos los modelos lineales: mismas predictoras y como
respuesta cada sonda.  

```{r}
#| label: limma4
fit = limma::lmFit(gse25171,design)
```

```{r}
#| label: limma5
#| results: false
fit$coefficients
coef(fit)
coefficients(fit)
```
- Vemos los primeros coeficientes y errores estándar estimados.

```{r}
#| label: limma6
head(coef(fit),n=2)
head(fit$sigma,n=2)
```
- Aplicamos el método limma.  

```{r}
#| label: limma7
fit1 = limma::eBayes(fit)
```
- Tenemos los t-estadísticos moderados con
```{r}
#| label: limma8
head(fit1$t)
```
- Los valores $B$ o logaritmo natural del cociente de odds a posteriori.

```{r}
#| label: limma9
head(fit1$lods)
```
- El valor estimado de $s_0^2$ y de $d_0$.
```{r}
#| label: limma10
fit1$s2.prior
fit1$df.prior
```
- Tenemos los estimadores de los errores estándar a posteriori.  

```{r}
#| label: limma11
head(fit$sigma)
head(fit1$s2.post)
```

- La función `limma::topTable()` nos hace un resumen de los
ajustes.
- Ordena (`sort.by`) las sondas por su p-valor original.

- Tenemos el coeficiente estimado en el ajuste `logFC` cuyo nombre
induce a confusión.

```{r}
#| label: limma12
limma::topTable(fit1,coef=2,adjust ="BH")
```
- ¿Cuántas sondas tienen un coeficiente significativamente
no nulo?

```{r}
#| label: limma13
padj = limma::topTable(fit1,coef=2,adjust ="BH",
                       number=nrow(gse25171))[,"adj.P.Val"]
table(padj < 0.05)
```

## time + Pi  

- Construimos la matriz de modelo.
```{r}
#| label: limma14
design = model.matrix(~  pData(gse25171)[,"time"] +
                          pData(gse25171)[,"Pi"] )
colnames(design) = c("constante","time","Pi")
head(design)
```

- Ajustamos los modelos. 

```{r}
#| label: limma15
fit = limma::lmFit(gse25171,design)
fit1 = limma::eBayes(fit)
```
- Evaluamos si los coeficientes correspondientes a la columna
2 `coef=2` son nulos.
- La columna 2 corresponde con la variable `time`.
```{r}
#| label: limma16
tt2 = limma::topTable(fit1,coef="time",adjust ="BH",
                      sort.by="none",number=nrow(gse25171))
head(tt2,n=20)
```

-  Evaluamos si los coeficientes correspondientes a la columna
3 `coef=3` son nulos. 
- La columna 3 corresponde con la variable `Pi`.

```{r}
#| label: limma17
tt3 = limma::topTable(fit1,coef="Pi",adjust ="BH",
                      sort.by="none",number=nrow(gse25171))
```
- ¿Cualés y cuántas sondas tienen el coeficiente correspondiente a
`time` significativamente no nulo ajustando por el método
de Benjamini-Hochberg?
```{r}
#| label: limma18
tt2.row = which(tt2[,"adj.P.Val"] < .05)
length(tt2.row)
```
- ¿Cualés y cuántas sondas tienen el coeficiente correspondiente a
`Pi` significativamente no nulo ajustando por el método
de Benjamini-Hochberg?
```{r}
#| label: limma19
tt3.row = which(tt3[,"adj.P.Val"] < .05)
length(tt3.row)
```
- ¿Y cuántas sondas tienen coeficiente significativamente no nulo
tanto para una como para la otra variable?
```{r}
#| label: limma20
intersect(tt2.row,tt3.row)
```
- ¿Qué sondas son?
```{r}
#| label: limma21
tt2[intersect(tt2.row,tt3.row),"PROBEID"]
```
¿A qué genes corresponden utilizando su código `ENTREZID`?
```{r}
#| label: limma22
tt2[intersect(tt2.row,tt3.row),"ENTREZID"]
```

## time * Pi  

- ¿Hay interacción entre las variables predictoras `time` y `Pi`? 
- Observemos que sustituimos `+` por `*`.

- La nueva columna de la matriz de modelo se define como el producto
de las columnas 2 y 3 correspondientes a las variables originales
`time` y `Pi`.  

- El modelo es distinto, los coeficientes estimados para la nueva variable
y para las ya existentes son distintos y no necesariamente tendremos
las mismas sondas.  

```{r}
#| label: limma23
design = model.matrix(~ pData(gse25171)[,"time"] *
                          pData(gse25171)[,"Pi"])
colnames(design) = c("constante","time","Pi","time:Pi")
head(design)
```
- ¿Qué sondas tienen interacciones significativas?  
```{r}
#| label: limma24
fit = limma::lmFit(gse25171,design)
fit1 = limma::eBayes(fit)
tt2 = limma::topTable(fit1,coef=2,adjust ="BH",
                      sort.by="none",number=nrow(gse25171))
tt3 = limma::topTable(fit1,coef=3,adjust ="BH",
                      sort.by="none",number=nrow(gse25171))
tt4 = limma::topTable(fit1,coef=4,adjust ="BH",
                      sort.by="none",number=nrow(gse25171))
```
- Ahora podemos evaluar los que tienen interacciones significativas.  
- Cuando una sonda tiene una interacción positiva se suele asumir que los
efectos principales no se deben eliminar del modelo.   
- Si hay interacción quiere decir que influyen significativamente y 
además de un modo distinto.  
```{r}
#| label: limma25
tt4.row = which(tt4[,"adj.P.Val"] < .05)
tt4[tt4.row,c("PROBEID","ENTREZID")]
```

## time2 * Pi  

- Hemos utilizado la variable `time` como numérica (y lo es).
- Podemos ver que en las variables fenotípicas hay una
variable `time2` categórica con dos niveles.

```{r}
#| label: limma26
pData(gse25171)[,"time2"]
```
- Ajustamos modelo con `time2` y `Pi`.

```{r}
#| label: limma27
design = model.matrix(~ pData(gse25171)[,"time2"] *
                          pData(gse25171)[,"Pi"])
colnames(design) = c("constante","time2","Pi","time2:Pi")
head(design)
```

```{r}
#| label: limma28
fit = limma::lmFit(gse25171,design)
fit1 = limma::eBayes(fit)
tt2 = limma::topTable(fit1,coef=2,adjust ="BH",
                      sort.by="none",number=nrow(gse25171))
tt3 = limma::topTable(fit1,coef=3,adjust ="BH",
                      sort.by="none",number=nrow(gse25171))
tt4 = limma::topTable(fit1,coef=4,adjust ="BH",
                      sort.by="none",number=nrow(gse25171))
tt4.row = which(tt4[,"adj.P.Val"] < .05)
tt4[tt4.row,c("PROBEID","ENTREZID")]
```
- Evaluamos si no hay interacción pero sí efecto de `time2`.

```{r}
#| label: limma29
sel = which(tt4[,"adj.P.Val"] >= .05 & tt2[,"adj.P.Val"] < .05)
length(sel)
head(tt4[sel,c("PROBEID","ENTREZID")])
```

- Podemos evaluar si no hay interacción pero sí efecto de `Pi`.
```{r}
#| label: limma30
sel = which(tt4[,"adj.P.Val"] >= .05 & tt3[,"adj.P.Val"] < .05)
length(sel)
```
- Esta forma de plantear el análisis supone cada uno de los coeficientes se ha
de interpretar como la modificación respecto de un grupo de referencia.
- Sería cuando las variables categóricas toman el primer nivel.

- El grupo de referencia corresponde `time2 = Short` y `Pi = Control`.

```{r}
#| label: limma31
levels(pData(gse25171)[,"time2"])
levels(pData(gse25171)[,"Pi"])
```

## time2Pi  

- Vamos a trabajar con contrastes.

- Construimos variable que combina ambos factores.

```{r}
#| label: limma32
time2Pi = vector("list",ncol(gse25171))
for(i in seq_along(time2Pi))
  time2Pi[[i]] = paste0(pData(gse25171)[,"time2"][i],
                          pData(gse25171)[,"Pi"][i])
(time2Pi = factor(unlist(time2Pi)))
```  

- Vamos a ajustar un modelo en donde no vamos a incorporar una
constante sino que las variables predictoras nos van a indicar cada
una de las categorías que tenemos.  

```{r}
#| label: limma33
design = model.matrix(~  0 + time2Pi) ## Quitamos constante (también -1)
colnames(design) = levels(time2Pi)
head(design)
``` 
- Realizamos los ajustes.  

```{r}
#| label: limma34
fit = limma::lmFit(gse25171,design)
```  
- Definimos los contrastes.  

```{r}
#| label: limma35
cont.matrix = limma::makeContrasts(
    dif1 =  (MediumControl + MediumTreatment)-
        (ShortControl + ShortTreatment),
    dif2 =  (MediumControl + ShortControl)-
        (MediumTreatment + ShortTreatment),
    dif3 = (MediumControl - ShortControl),
    dif4 = (MediumTreatment - ShortTreatment),
    dif5 = (MediumControl - ShortControl) -
        (MediumTreatment - ShortTreatment),
    levels = design)
```  
- Hemos construido la siguiente matriz de contrastes.
```{r}
#| label: limma36
cont.matrix
```  
- Aplicamos limma.

```{r}
#| label: limma37
fit2 = limma::contrasts.fit(fit,cont.matrix)
fit2 = limma::eBayes(fit2)
``` 
- Obtenemos resúmenes.
```{r}
#| label: limma38
tt1 = limma::topTable(fit2,coef=1,adjust ="BH",
                      sort.by="none",number=nrow(gse25171))
tt2 = limma::topTable(fit2,coef=2,adjust ="BH",
                      sort.by="none",number=nrow(gse25171))
tt3 = limma::topTable(fit2,coef=3,adjust ="BH",
                      sort.by="none",number=nrow(gse25171))
tt4 = limma::topTable(fit2,coef=4,adjust ="BH",
                      sort.by="none",number=nrow(gse25171))
tt5 = limma::topTable(fit2,coef=5,adjust ="BH",
                      sort.by="none",number=nrow(gse25171))
```  
- Nos fijamos como ilustración en el contraste `dif5` que evalua 
la interacción.

```{r}
#| label: limma39
sel = which(tt5[,"adj.P.Val"] < .05)
length(sel)
head(tt4[sel,c("PROBEID","ENTREZID")])
```  

## ANOVA 

- Podemos plantearnos si las medias en las cuatro poblaciones
definidas por `time2Pi` son iguales.

- Definimos la matriz de modelo.  

```{r}
#| label: limma39b
design = model.matrix(~  time2Pi) ## Con la constante  
colnames(design) = c("intercept",levels(time2Pi)[-1])
head(design)
```
- Ajustamos los modelos. 

```{r}
#| label: limma39c
fit = limma::lmFit(gse25171,design)
fit = limma::eBayes(fit)
tt = limma::topTable(fit,coef=2:4,adjust ="BH",
                      sort.by="none",number=nrow(gse25171))
head(tt,n=2)
```
- ¿Cuántos genes tienen p-valor significativo?

```{r}
#| label: limma39d
table(tt$adj.P.Val < .05)
```

# gse44456  

## Datos  

- Leemos los datos `tamidata::gse44456`.
```{r}
#| label: limma40
library(Biobase)
data(gse44456,package="tamidata")
```
- Vamos a considerar la variable que nos da el tiempo desde la muerte 
de la persona.
```{r}
#| label: limma41
pData(gse44456)[,"postmortenint"]
```
- Cargamos Limma.

```{r}
#| label: limma42
pacman::p_load(limma)
```
- Determinamos la matriz de diseño.
```{r}
#| label: limma43
design = model.matrix(~ pData(gse44456)[,"postmortenint"])
colnames(design) = c("intercept","postmortenint")
```
- Podemos ver que la matriz `design`.
```{r}
#| label: limma44
head(design)
```
- Ajustamos los modelos lineales.
```{r}
#| label: limma45
fit = lmFit(gse44456,design)
```
Los coeficientes ajustados los tenemos con

```{r}
#| label: limma46
head(fit$coefficients)
```

Ajustamos el modelo empírico bayesiano. 
```{r}
#| label: limma47
fit1 = eBayes(fit)
```

- Si la constante por la que multiplicamos growthrate es nula indicará
que no hay una dependencia del nivel de expresión respecto de la tasa
de crecimiento celular. 
- Veamos los resultados.  
```{r}
#| label: limma48
topTable(fit1,coef=2,adjust ="BH")
```

- Aunque alguno de los p-valores originales pudieran (marginalmente)
considerarse como significativos cuando corregimos por comparaciones múltiples
por el método de Benjamini-Hochberg no lo son.  

## Un diseño cruzado con limma
- Seguimos utilizando tamidata::gse44456.
```{r}
#| label: limma49
pacman::p_load("limma")
data(gse44456,package="tamidata")
```
- Observamos los datos fenotípicos o metadatos.
```{r}
#| label: limma50
head(pData(gse44456))
```
- Analizamos la posible dependencia con el alcoholismo recogida
en la covariable `case`
```{r}
#| label: limma51
alcoholism = pData(gse44456)[,"case"] # Simplificamos el código
design = model.matrix(~ 0 + alcoholism)
colnames(design) = c("control","alcoholic")
```
- Ajustamos el modelo.  

```{r}
#| label: limma52
fit = lmFit(gse44456,design)
```
```{r}
#| label: limma53
(contrast.matrix = makeContrasts(dif = control - alcoholic,levels = design))
```
- Estimamos.
```{r}
#| label: limma54
fit2 = contrasts.fit(fit,contrast.matrix)
fit2 = eBayes(fit2)
```
- Veamos cuáles son significativos. 
```{r}
#| label: limma55
topTable(fit2,coef=1,adjust="BH")
```

- En el segunda análisis nos fijamos en dos covariables, el alcoholismo 
y el sexo de la persona. 
- Es un diseño factorial $2 \times 2$.
- Veamos cuántos datos tenemos en cada celda. 
```{r}
#| label: limma56
table(pData(gse44456)[,"case"],pData(gse44456)[,"gender"])
```

- En un diseño como este las preguntas habituales son las siguientes:
    1. ¿Qué genes muestran un comportamiento diferenciado o relacionado con 
    el alcoholismo? 
    2. ¿Qué genes se comportan de un modo distinto para hombres y mujeres?
    3. ¿Para qué genes los cambios en su expresión según el alcoholismo 
    son distintos en cada uno de los sexos?
- En jerga estadística hablaríamos de los efectos principales de los factores y
de la posible interacción. 

- Una aproximación simple y efectiva es construir un solo factor con todas
las combinaciones de los factores. 
```{r}
#| label: limma57
casegender = vector("list",ncol(gse44456))
for(i in seq_along(casegender))
  casegender[[i]] = paste(pData(gse44456)[,"case"][i],
                          pData(gse44456)[,"gender"][i],sep="")
casegender = factor(unlist(casegender))
```

- Consideremos la siguiente matriz de diseño.
```{r}
#| label: limma58
design = model.matrix(~  0 + casegender)
colnames(design) = levels(casegender)
```
- Podemos ver que cada coeficiente corresponde con la expresión media
para la correspondiente combinación de factores.
```{r}
#| label: limma59
fit = limma::lmFit(gse44456,design)
```
- Construimos los contrastes en que estamos interesados.
```{r}
#| label: limma60
cont.matrix = makeContrasts(
    dif1 = controlmale - alcoholicmale,
    dif2 = controlfemale - alcoholicfemale,
    dif12 = (controlmale - alcoholicmale)- (controlfemale - alcoholicfemale),
    levels = design)
fit2 = contrasts.fit(fit,cont.matrix)
fit2 = eBayes(fit2)
```
- Podemos ejecutar el siguiente código y podemos comprobar que ningún 
contraste es significativo. 
```{r}
#| label: limma61
topTable(fit2,coef=1,adjust="BH")
topTable(fit2,coef=2,adjust="BH")
topTable(fit2,coef=3,adjust="BH")
```

## Bibliografía  
