---
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:  
      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)
```

# 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_{ijk}) = \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
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]
```
- 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"[to_keep]
                       + colData(tcga_coad)$"age_at_diagnosis"[to_keep])
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