---
title: "Expresión diferencial con datos RNA-seq"
author: "Guillermo Ayala Gallego"  
date: "`r Sys.Date()`"  
address: Departamento de Estadística e I.O. Universidad de Valencia  
bibliography: /home/gag/Nextcloud/BIBLIOGRAFIA/Bibliografia.bib
format: 
   revealjs:  
      smaller: true
      scrollable: true
      echo: true
      embed-resources: true  
   pdf:
       toc: true
       number-sections: false
       colorlinks: true
---

# Introducción  

## Paquetes


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

```{r  }
#| label: rna2
pacman::p_load(SummarizedExperiment,edgeR,ggplot2)
```

## Data 

```{r}
#| label: rna3
data(PRJNA297664,package="tamidata")
```

```{r}
#| label: rna4
colData(PRJNA297664)[,"treatment"]
```

# Una muestra por condición

## Un test para comparar proporciones 

- Queremos comparar dos condiciones y disponemos de una muestra en cada una
de las condiciones.
- En este caso para cada gen tendríamos una tabla $2\times2$ como


    Muestra    |       Gen            |    Resto                  |     Total  
-------------  |     --------------   |  -----------------------  |  ------------------   
       1       |        $y_1$         |     $m_1-y_1$             |     $m_1$              
       2       |        $y_2$         |     $m_2-y_2$             |     $m_2$              
     Total     |        $y_1 + y_2$   |     $m_1+m_2 -y_1 - y_2$  |     $m_1 + m_2$  



## edgeR::binomTest

- Para el gen i-ésimo tendríamos los conteos $x_{i1}$ e $x_{i2}$ en las dos
muestras siendo $m_1$ y $m_2$ los tamaños de las librerías.
- Suponemos fijo el número total de muestras para el gen i 
($n_i = x_{i1} + x_{i2}$)
y tamaño total de las dos librerías $m_1$ y $m_2$ 
- Bajo estas dos hipótesis previas vamos a contrastar que 
$H_i: p_1 = m_1/(m_1+m_2)$ frente a $H_i: p_1 \neq m_1/(m_1+m_2)$

- Como ilustración nos fijamos en las muestras 1 y 3. 


```{r}
#| label: rna5
pvalor13 = edgeR::binomTest(assay(PRJNA297664)[,1], assay(PRJNA297664)[,3])
```

```{r}
#| label: rna6
pacman::p_load(ggplot2)
df = data.frame(pvalor13)
ggplot2::ggplot(df,aes(x=pvalor13)) + geom_density()
```


```{r}
#| label: rna7
pares = rbind(c(1,3),c(1,5),c(1,6),c(2,3),c(2,5),c(2,6),c(4,3),
    c(4,5),c(4,6))
pvalores = apply(pares,1,
                 function(x) binomTest(assay(PRJNA297664)[,x[1]],
                                       assay(PRJNA297664)[,x[2]]))
```

```{r}
#| label: rna8
y = apply(pvalores,1,function(x) c(mean(x),min(x),max(x)))
df = data.frame(t(y))
names(df) = c("mean","min","max")
df=df[sort(df[,"mean"],index.return=TRUE)$ix,] 
df1 = reshape2::melt(df,id="mean") 
pp = ggplot(df1,aes(x=mean,y=value,color=variable))
pp + geom_smooth()  
```

# edgeR clásico  

## Introducción  

- Para dos grupos.  
- Es un test exacto.  
- Propuesto en @RobinsonSmyth08 y @RobinsonSmyth07 (en este orden).

## Estimación de una dispersión común

- Consideramos una característica (gen por ejemplo) en $n$ muestras (o librerías) de
  tamaños distintos.
- $m_i$,  tamaño de la $i$-ésima librería (total de lecturas).
- $\lambda$, la proporción que hay en una librería cualquiera de la 
  característica en que estamos interesados.
- Si $Y \sim NB(\mu,\phi)$ entonces la función de probabilidad viene dada por
$$
f(y|\mu,\phi) = P(Y=y|\mu,\phi) = \frac{\Gamma(y+\phi^{-1})}{\Gamma(\phi^{-1})
\Gamma(y+1)} \bigg( \frac{1}{1 + \mu \phi} \bigg)^{\phi^{-1}}
\bigg( \frac{\mu}{\phi^{-1} + \mu} \bigg)^{y}.
$$
- Se tiene:
    +  $E(Y) = \mu$,  
    +  $var(Y) = \mu + \phi \mu^2$.  
- **Asumimos**:  $Y_i \sim NB(m_i \lambda, \phi)$.

## Si las librerías tienen el mismo tamaño  

- **Si** todas las librerías tienen el mismo tamaño: $m_i = m$ para $i=1,\ldots,n$
entonces 
$$Z = \sum_{i=1}^n Y_i \sim NB(nm\lambda,\phi/n).$$  
- La **logverosimilitud condicionada** al valor $Z=z$ no depende
  de $\lambda$ y podemos estimar el valor de $\phi$.   
- Se estimaría $\phi$ maximizando la logverosimilitud condicionada
  (a la suma de todos los conteos del gen) dada por 
  $$\ell_{\mathbf{y} | z} = \bigg [ \sum_{i=1}^n \log \Gamma(y_i + 1/\phi) \bigg ] + \log \Gamma(n/\phi) - \log \Gamma(z + n/\phi) - n \log \Gamma(1/\phi),$$
  donde ${\mathbf y} =(y_1,\ldots,y_n)'$ son los conteos observados.  
- En el caso en que los tamaños de las librerías son distintos la verosimilitud 
  no tiene una expresión simple.  
- Modificamos los  $y_i$ observados generando unos **pseudodatos** con el mismo 
  tamaño de librería dado por 
  $$m^* = \bigg (\prod_{i=1}^n m_i \bigg )^{\frac{1}{n}}.$$  
- ¿Y cómo transformamos los conteos $y_i$?  

## qCML: Quantile-adjusted condicional maximum likelihood

1. Inicializamos $\phi$ (por ejemplo, con el estimador maxímo verosímil 
condicionado sin realizar ningún ajuste).

2. Dado el valor estimado de  $\phi$, estimamos $\lambda$ maximizando
la verosimilitud para el valor dado de la dispersión.

3. Suponemos que cada conteo $y_i$ es un valor observado de una 
   distribución binomial negativa con media $m_i \lambda$ y parámetro de
   dispersión $\phi$. Calculamos los percentiles
   $$p_i = P(Y < y_i|m_i \lambda,\phi) + \frac{1}{2} P(Y = y_i|m_i \lambda,\phi),$$
   para $i=1,\ldots,n$. 

4. Suponemos ahora una distribución binomial negativa con media $m^* \lambda$ 
   y dispersion $\phi$.
    1. Determinamos qué valor sería el percentil de orden $p_i$ en la nueva 
    distribución. 
    Ya no tiene porqué ser un dato entero e incluso puede ser negativo.
    2. Los pseudodatos para un mismo gen tienen aproximadamente la misma 
    distribución.  

5. Estimamos la dispersión $\phi$ con los pseudodatos utilizando la veromilitud
	condicionada a la (pseudo) suma total de un gen.  

6. Repetimos desde 2 hasta 5 hasta que converja $\phi$.  

##  Un test exacto para dos grupos  
-  $y_{ijk}$ el conteo para la **muestra** $k$ en la **condición** $j$ del gen $i$
      donde $j=1,2$ y $k=1,\ldots,n_j$.  
- Los tamaños de la librería k de la condición j será $m_{jk}$.  
- Los tamaños de las librerías **no** son iguales dentro de cada clase.  
- Aplicamos el procedimiento qCML dentro  de cada clase.  
- Tenemos un tamaño común $m_j$ en la condición $j$.  
- Estimamos $\phi$ maximizando 
  $$
  l(\phi) = \sum_{i=1}^N \ell_i (\phi).
  $$
  con 
  $$\ell_i (\phi) = \sum_{j=1}^2 \bigg ( \sum_{k=1}^{n_j}
  \log \Gamma(y_{ijk} + \phi^{-1}) +  \log \Gamma(n_j \phi^{-1}) - 
  \log \Gamma(z_{ij} + n_j \phi^{-1}) - n_j
  \log \Gamma(\phi^{-1})\bigg ).$$
  siendo $z_{ij}= y_{ij\cdot} = \sum_{k=1}^{n_j} y_{ijk}$.
  
- El estimador de $\phi$ será $\hat{\phi}_C$.  

## Contraste de hipótesis  
 
- Asumimos: $EY_{ijk} = m_{jk} \lambda_{ij}$.  
- La hipótesis nula, **no hay diferencias entre las medias de los conteos
en las dos condiciones para el $i$-ésimo gen** se formularía como

$H_i: \lambda_{i1}  = \lambda_{i2},$ vs $K_i: \lambda_{i1}  \neq \lambda_{i2}$.  
- Bajo $H_i$,  $\lambda_{i1} = \lambda_{i2} = \lambda_i$.  
- Aplicamos  **qCML** a **todas** las muestras: $y_{ijk}$ son los pseudodatos.  
- Utilizando el estimador $\hat{\phi}_C$ y los conteos $y_{ijk}$
  podemos estimar $\lambda_i$.  
- Bajo la hipótesis nula de no diferencia entre grupos tendríamos que
$Y_{ij\cdot} = \sum_{k=1}^{n_j} Y_{ijk} \sim NB(n_j m^* \hat{\lambda}_i,
\hat{\phi}_C/n_j)$.

- $Y_{i1\cdot}$ e $Y_{i2\cdot}$ son independientes.

- La suma $Y_{i1\cdot} + Y_{i2\cdot}$ también tiene una distribución binomial negativa:
 $Y_{i1\cdot} + Y_{i2\cdot} = \sum_{i=1}^2 \sum_{k=1}^{n_j} Y_{ijk} 
\sim NB((n_1+n_2) m^* \hat{\lambda}_i,\hat{\phi}_C/(n_1+n_2))$

- Podemos considerar la distribución condicionada del vector aleatorio
$(Y_{i1\cdot},Y_{i2\cdot})$ a la suma $Y_{i1\cdot} + Y_{i2\cdot}$ y considerar las 
probabilidades de los conteos conjuntos que **son menos probables que el 
observado**.
- La suma de estas probabilidades nos daría el p-valor del test.


## Dispersiones posiblemente distintas  
- Una misma dispersión sobre una cantidad grande de genes no es muy razonable.
- Buena desde el punto de vista del modelo y su estimación pero nos aleja de 
los datos.
- Supongamos dispersiones posiblemente distintas para los distintos genes,
$\phi_i$.
- Proponen
$$WL(\phi_i) = \ell_i(\phi_i) + \alpha \ell(\phi_i)$$
siendo $\alpha$ el peso que se da a la verosimilitud global.
- Es una función que propone un compromiso entre considerar $l$ como función a
maximizar considerando la misma contribución a todos los genes y $l_i$ en
donde solamente consideramos los conteos del propio gen.
- En @RobinsonSmyth08 proponen un método de estimación de $\alpha$.


# Un análisis con edgeR  

## DGEList  

- Tenemos un  `SummarizedExperiment`.
- Construimos un `DGEList`.

```{r }
#| label: rna10
x = DGEList(counts = assay(PRJNA297664),
            group = colData(PRJNA297664)[,"treatment"])
```
- La clase del nuevo objeto es 
```{r }
#| label: rna11
class(x)
```
Y sus atributos son 
```{r}
#| label: rna12
attributes(x)
```
La matriz de conteos  
```{r}
#| eval: false
#| label: rna13
x$counts
```
Tenemos la componente `samples`. 

```{r}
#| label: rna14
class(x$samples)
```
Las primeras filas son 

```{r}
#| label: rna15
head(x$samples)
```
Tenemos un `factor`. 
```{r}
#| label: rna16
x$samples[,"group"]
```
Los tamaños de las librerías son 
```{r}
#| label: rna17
x$samples[,"lib.size"]
```
Los factores de normalización son 
```{r}
#| label: rna18
x$samples[,"norm.factors"]
```
Veamos una descripción de los tamaños de las librerías.
```{r}
#| label: rna19
summary(x$samples[,"lib.size"])
```
```{r}
#| label: rna20
ggplot(x$samples,aes(x= lib.size)) + geom_dotplot()
```

## Eliminando genes con conteos bajos  

Utilizamos los conteos por millón CPM en lugar de los conteos originales. 
Obligamos a que se recalculen los tamaños de las librerías.  
```{r }
#| label: rna21
 keep = rowSums(cpm(x)>1) >= 2
 x <- x[keep, , keep.lib.sizes=FALSE]
```
```{r}
#| label: rna23
dge = DGEList(counts=assay(PRJNA297664),
              group=colData(PRJNA297664)[,"treatment"])
dge.c = estimateCommonDisp(dge) ##Estimamos dispersión común
dge.c$common.dispersion
dge.t =  estimateTagwiseDisp(dge.c) ##Dispersiones por gen
et.c = exactTest(dge.c)
et.t = exactTest(dge.t)
```
Y vemos los resultados.  
```{r}
#| label: rna24
topTags(et.c)
topTags(et.t)
```

# edgeR utilizando modelo lineal generalizado  

## Modelo  

- $Y_{ij}$ el conteo aleatorio (número de lecturas alineadas) para
el gen $i$ en la muestra $j$. 
- Denotamos por $m_j = \sum_{i=1}^N y_{\cdot j}$
la profundidad de secuenciación o total de lecturas de la muestra $j$.
- Utilizamos como función de enlace el logarítmo natural.
- Consideramos la profundidad de secuenciación como offset (un modelo de tasas
sobre la profundidad de secuenciación).
- El modelo para la media es 
$$
\ln \mu_{ij} = \mathbf x_j^T \mathbf \beta_i + \ln m_j.
$$
- En el modelo las variables predictoras son comunes a todos los genes. 
- Asumimos que la componente aleatoria sigue una distribución binomial 
negativa (con el parámetro de dispersión conocido).
- Entonces 
$$
var(Y_{ij}) = \mu_{ij} + \phi_i \mu_{ij}^2,
$$
siendo $\phi_i$ el parámetro de dispersión que hemos de asumir conocido
o, de otro modo, tenemos que estimarlo previamente. 

##  

- En @McCarthyChenSmyth2012 muestran cómo estimar por máxima
verosimitud el vector de coeficientes $\mathbf \beta_i$.
- Utilizan una modificación de los mínimos cuadrados iterativamente reponderados
(IRWLS). 
- El parámetro de dispersión se estima maximizando la
logverosimilitud penalizada definida como
$$
APL_i(\phi_i) = \ell(\phi_i;\mathbf y_i,\hat{\mathbf \beta}_i) -
\frac12 \ln |\mathbb I_i|
$$
siendo: 
    - $\mathbf y_i$ los conteos para el gen $i$, 
    - $\hat{\mathbf \beta}_i$ el vector de coeficientes, 
    - $\ell()$ es la función de logverosimilitud
    -  $|\mathbb I_i|$ el determinante de la matriz de información de Fisher
	para el $i$-ésimo gen. 

## TCGA-COAD  

```{r}
#| label: rna25
pacman::p_load(edgeR,SummarizedExperiment)
load(paste0(dirTamiData,"tcga_coad.rda"))
```
- Nos centramos en las variables fenotípicas  `age_at_diagnosis` y 
`tissue_or_organ_of_origin`.  
```{r}
#| label: rna26
colData(tcga_coad)
summary(colData(tcga_coad)[,"age_at_diagnosis"])
table(colData(tcga_coad)[,"tissue_or_organ_of_origin"])
```  
- Hemos de eliminar aquellas muestras que tienen las variables predictoras
con datos faltantes ya que las funciones que siguen no los admiten.  

```{r}
#| label: rna27
torm1 = which(is.na(colData(tcga_coad)$"age_at_diagnosis"))
torm2 = which(is.na(colData(tcga_coad)$ "tissue_or_organ_of_origin"))
toremove = union(torm1,torm2)
tcga_coad = tcga_coad[,-toremove]
```
- Construimos el objeto `DGEList` sin indicar ninguna variable
`group` ni ninguna matriz de modelo y eliminamos genes con conteos bajos. 
```{r}
#| label: rna28
dge = DGEList(counts=assay(tcga_coad))
to_keep = rowSums(cpm(dge) > 0.5) > 20
dge =  dge[to_keep,keep.lib.sizes=FALSE]
dim(dge)
```
- Construimos la matriz de modelo con las dos variables predictoras,
una de caracter categórico y la otra numérica. 
- Cambiamos los nombres de las columnas de la matriz de modelo.
```{r}
#| label: rna29
design0 = model.matrix(~  0 +
         colData(tcga_coad)$"tissue_or_organ_of_origin"
       + colData(tcga_coad)$"age_at_diagnosis")
y = levels(colData(tcga_coad)$"tissue_or_organ_of_origin")
y = sapply(y,function(x) gsub(" ","_",x)) ## Eliminamos espacios
y = sapply(y,function(x) gsub(",","_",x)) ## Eliminamos las comas
colnames(design0) = c(y,"age_at_diagnosis")
```
- Estimamos las dispersiones por tres métodos distintos:  
    - Asumiendo una dispersión común,  
    - una por gen y  
     - con una relación media-varianza.  
```{r}
#| label: rna29b
#| eval: false 
dge = estimateDisp(dge,design=design0)
```
Si solo queremos una de las tres opciones podemos usar las funciones
`estimateGLMCommonDisp()`, `estimateGLMTagwiseDisp()` y
`estimateGLMTrendedDisp()`.
```{r}
#| label: rna30
#| echo: false
load(paste0(dirTamiData,"tcga_coad_dge_design0.rda"))
```
```{r}
#| label: rna31
plotBCV(dge)
```
```{r}
#| label: rna32
#| eval: false 
fit = glmFit(dge, design=design0)
```
```{r}
#| label: rna33
#| echo: false
load(paste0(dirTamiData,"tcga_coad_dge_design0_fit.rda"))
```
- Veamos si influye la variable `age_at_diagnosis`. 
- Si observamos la matriz  de modelo `design0` corresponde con
la columna 10 de la matriz de modelo. 
- Se realiza un test del cociente de verosimilitudes. 

```{r}
#| label: rna34
lrt1 = glmLRT(fit,coef="age_at_diagnosis")
lrt1 =  glmLRT(fit,coef=10) ## Equivalente a la línea anterior
topTags(lrt1)
```
- Podemos evaluar toda la variable `tissue_or_organ_of_origin`.
```{r}
#| label: rna35
lrt2 = glmLRT(fit,coef=1:9)
topTags(lrt2)
```
- Y elegir los contraste que queramos. 
- Mostramos una comparación entre dos grupos.

```{r}
#| label: rna36
AD = makeContrasts(contrast1 = Ascending_colon - Descending_colon,
                   levels=design0)
lrt3 = glmLRT(fit,contrast = AD)
topTags(lrt3)
```

## Bibliografía  
