---
title: "Expresión diferencial (marginal)"
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  
---  


## El problema

- ¿Hay diferencias de la expresión de genes entre dos (o más) grupos distintos?
- Dicho de otro modo, la expresión de un gen es distinta bajo distintos tratamientos. 
- Se adopta la aproximación gen-a-gen, esto es, de momento solamente buscamos 
genes que se expresan diferencialmente sin atender al comportamiento conjunto 
que puedan tener. 

# Log fold-change

## El logaritmo del cociente de las medias

- Queremos comparar dos grupos.
- $x_{ij}$ ($y_{ij}$) es la expresión de la i-ésima característica en 
la j-ésima muestra del primer grupo (segundo grupo). 
- Tomamos los logaritmos (en base 2 o `log2`) de los valores originales: 
$$
u_{ij} = \log_2(x_{ij}), \ v_{ij} = \log_2(y_{ij}).
$$ 
- Para cada gen, tendremos las expresiones medias  para la i-ésima característica
$$\bar{x}_{i \cdot} = \sum_{j=1}^{n_1} \frac{x_{ij}}{n_1}$$ 
$\bar{y}_{i \cdot}$, $\bar{u}_{i \cdot}$, $\bar{v}_{i \cdot}$. 
- ¿Qué se  entiende por **fold-change**?
- Dos son las interpretaciones **distintas** de este valor. 
$$
  FC^{(1)}_i = \frac{\bar{x}_{i \cdot}}{\bar{y}_{i \cdot}}.
$$
y 
$$
  FC^{(2)}_i = \bar{u}_{i \cdot} -  \bar{v}_{i \cdot},
$$

 - El log-fold change es el logaritmo en base 2 de los fold-change que acabamos
 de definir. 
- Si el cociente anterior es mayor que `c` entonces diríamos que el gen
se expresa de un modo diferencial en ambos grupos.
- Tenemos una expresión diferencial del gen $i$ si

$$\bigg |\log_2 \big ( FC_i \big ) \bigg| \geq c$$

# Comparando dos grupos y un solo gen

## Dos grupos 

- Utilizamos los datos gse21942.

```{r }
#| label: edm19
data(gse21942,package="tamidata")
```
- Nos fijamos en un gen.

```{r }
#| label: edm20
library(Biobase)
x0 = exprs(gse21942)[2058,]
y0 = pData(gse21942)[,"FactorValue..DISEASE.STATE."]
```

## Perfil de expresión


```{r }
#| label: edm21
pacman::p_load(ggplot2)
df = data.frame(id=1:length(x0),x0,y0)
ggplot(df,aes(x=id,y=x0)) + geom_point(aes(colour=y0))
```

## Un diagrama de cajas

```{r}
#| label: edm22
df = data.frame(id=1:length(x0),x0,y0)
ggplot(df,aes(x=y0,y=x0)) + geom_boxplot()
```

## Comparamos dos condiciones

- Tenemos datos relativos a niveles de expresión bajo dos condiciones
experimentales.
- Denotamos por $X$ el nivel de expresión aleatorio que observamos 
bajo la primera condición y por $Y$ lo mismo pero  con la segunda condición.
- $X \sim N(\mu_X,\sigma^2_X)$
- $Y \sim N(\mu_Y,\sigma^2_Y)$
- Nos planteamos el contraste
$$
H_0: \mu_X = \mu_Y
$$
$$
H_1:  \mu_X \neq \mu_Y.
$$

## t-test con varianzas iguales

- Calculamos 
$$T = \frac{\bar{X} - \bar{Y}}{S_p \sqrt{\frac{1}{n}+\frac{1}{m}}}$$
con
$$S_p = \hat{\sigma}^2 = \frac{(n-1) S^2_X + (m-1) S^2_Y}{n+m-2}$$
siendo
$$S^2_X = \sum_{i=1}^n \frac{(X_i - \bar{X})^2}{n-1}.$$
y $S^2_Y$ lo mismo sustituyendo $X_i$ por  $Y_i$.



- Rechazamos la hipótesis nula cuando 
$$
|T| > t_{\nu,1-\alpha/2}
$$
- El p-valor viene dado por 
$$
p= P(|T|  \geq t_0) 
$$
siendo $t_0$ el valor observado de $T$ en cada caso.

## t.test  

```{r}
#| label: edm22b
t.test(x0 ~ y0,var.equal=TRUE)
```

## ...  

```{r,echo=FALSE}
#| label: edm23
tt = t.test(x0 ~ y0,var.equal=TRUE)
nhstplot::plotttest(t = tt$statistic, df = tt$parameter)
```


## Comparación simultánea de todos los genes

Hacemos un t-test para todos los genes

```{r}
#| label: edm24
library(genefilter)
tt = genefilter::rowttests(gse21942,pData(gse21942)[,"FactorValue..DISEASE.STATE."])
```

```{r}
#| label: edm25
head(tt)
```

## 
¿Cuándo tests mostrarían una expresión diferencial entre ambos grupos?

```{r}
#| label: edm26
table(tt$p.value < 0.01)
```
- Como vemos tenemos demasiados genes que muestran expresión diferencial.
- Y no parece muy razonable esto. 
- El nivel de significación $\alpha$ es una cota superior
del error tipo I cuando realizamos **un solo test**.
- Pero no estamos aplicando un solo test.
- Estamos rechazando muchas hipótesis nulas (no diferencia entre grupos
para cada gen) cuando no deberíamos de hacerlo. 


