index.Rmd

--- 
title: "Proyecto final RNA-seq"
author: "Salvador González Juárez"
date: "`r Sys.Date()`"
site: bookdown::bookdown_site
---

# <span style="color:darkblue">Proyecto</span>

Universidad Nacional Autónoma de México - Licenciatura en ciencias genómicas

Bioinformática y Estadística 2 2021 - Análisis de datos de secuenciación masiva (Dr. Leonardo Collado)

**Salvador González Juárez**

---

# <span style="color:orange">Especificaciones:</span>

- Con datos de algún estudio disponible vía ```recount3```, hagan un análisis de expresión diferencial.
- Incluyan al menos 3 gráficas en su reporte.
- Su reporte debe ser público y estar listado en el Google Sheet del curso.

---

# <span style="color:orange">Proyecto escogido:</span>

**Firmas transcripcionales específicas del sexo en la depresión humana - SRP115956**

- **Marco teórico:**

El **trastorno depresivo mayor** (**MDD** por sus siglas en inglés) es una enfermedad debilitante crónica que afecta cada año a 350 millones de personas en todo el mundo, lo que genera importantes cargas económicas y médicas en las sociedades. Si bien afecta tanto a hombres como a mujeres, el MDD se caracteriza por un fuerte dimorfismo sexual: las mujeres tienen 2-3 veces más probabilidades de desarrollar MDD, y presentan síntomas de mayor gravedad, mayor deterioro funcional, síntomas depresivos más atípicos y tasas más altas de ansiedad comórbida. Además, los hombres y las mujeres con MDD responden de manera diferente al tratamiento con antidepresivos. Sin embargo, los mecanismos moleculares que subyacen a este dimorfismo sexual siguen siendo en gran parte desconocidos.

- **Métodos:**

Combinando la expresión diferencial y los análisis de redes de coexpresión de genes, los investigadores proporcionaron una caracterización completa de los perfiles transcripcionales masculinos y femeninos asociados con el MDD en 6 regiones del cerebro. 

- **Resultados:**

Los resultados que presentaron los autores del estudio muestran una reordenación importante de los patrones de transcripción en el MDD, con una superposición limitada entre hombres y mujeres, el cual es un efecto observado en humanos deprimidos. Identificaron reguladores clave de las redes de genes específicos del sexo que subyacen al MDD y confirmaron su impacto específico del sexo como mediadores de la susceptibilidad al estrés. Por ejemplo, la regulación a la baja del gen *hub* específico de la mujer *DUSP6* en la corteza prefrontal (CPF) imita la susceptibilidad al estrés en las mujeres solo al aumentar la señalización ERK y la excitabilidad de las neuronas piramidales. Tal regulación a la baja de DUSP6 también recapitula la remodelación transcripcional que ocurre en la CPF de mujeres deprimidas. 

- **Conclusiones:**

Estos hallazgos revelan un dimorfismo sexual dramático a nivel transcripcional en el MDD y resaltan la importancia de estudiar tratamientos específicos por sexo para este trastorno.

- **Diseño general:**

Datos de secuenciación de ARN de 6 regiones cerebrales post mortem humanas en hombres y mujeres con y sin depresión mayor.

---

# <span1 style="color:orange">Datos de SRP115956</span>

Voy a usar datos de [10.1038/nm.4386](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5734943/) procesados con **recount3**. Por lo tanto, descargo los datos con los siguientes comandos:

```{r}
# Cargar el paquete de recount3
library("recount3")

# Encontrar el proyecto de interes
proj_info <- subset(available_projects(),
    project == "SRP115956" & project_type == "data_sources")

# Crear un objetio de tipo RangedSummarizedExperiment (RSE) con la informacion a nivel de genes
rse_gene <- create_rse(proj_info)

# Explorar el objeto RSE
rse_gene

# Guardar el objeto RSE original
rse_gene_original <- rse_gene

# Usar para restablecer el objeto RSE a su forma original
#rse_gene <- rse_gene_original
```

En resumen, los datos transcriptómicos estan representados por **63,856 genes** de **263 muestras**. Ahora, es necesario obtener el número de lecturas para cada una de las muestras, en lugar de los datos a nivel de nucleótidos. Además, con el siguiente código podemos explorar la información contenida en las categorias, que representan variables del experimento.

```{r}
# Obtener los numeros de lecturas
assay(rse_gene, "counts") <- compute_read_counts(rse_gene)

# Facilitar el uso de la informacion del experimento
rse_gene <- expand_sra_attributes(rse_gene)

# Explorar los parametros del experimento y algunas de sus variables
colData(rse_gene)[,
    grepl("^sra_attribute", colnames(colData(rse_gene)))]
```

De estas categorias, elegí como variable categórica el fenotipo, el cual se refiere a si el individuo presentaba el MDD o no. Además, me pareció interesante analizar el género de los individuos, ya que precisamente el estudio encontró, en esta categoría, diferencias significativas. Por último, la variable que tambien me llamó la atención es la causa de muerte, ya que algunos de los individuos se suicidaron, y esto podría tener una fuerte relación con el transtorno.

Los datos de estas variables, junto con el de las variables numéricas de edad del individuo y el número de integridad del RNA (RIN por sus siglas en inglés), deben ser corregidas como se muestra en el siguiente código. A continuación, presento el resumen de algunas de las variables, incluyendo las de mi interés.

```{r}
# Pasar de character a nuemric o factor
rse_gene$sra_attribute.phenotype <- factor(rse_gene$sra_attribute.phenotype)
rse_gene$sra_attribute.gender <- factor(rse_gene$sra_attribute.gender)
rse_gene$sra_attribute.age <- as.numeric(rse_gene$sra_attribute.age)
rse_gene$sra_attribute.rin <- as.numeric(rse_gene$sra_attribute.rin)

# Resumen de las variables de interes
summary(as.data.frame(colData(rse_gene)[,
    grepl("^sra_attribute.[phenotype|gender|age|rin]", colnames(colData(rse_gene)))]))
```

La categoría de causa de muerte no es categórica, ya que hay múltiples opciones en ella. Por lo tanto, las voy a asignar como suicidio o no suicidio. Esta clasificación si es apropiada, ya que solo hay dos opciones. 

```{r}
# Encontraremos diferencias entre muestra prenatalas vs postnatales
rse_gene$death <- factor(ifelse(rse_gene$sra_attribute.Cause_of_death == "Suicide", "suicide", "no suicide"))
table(rse_gene$death)
```

Con base en las variables de control de calidad se divide el número de lecturas asignadas a los genes entre el número de lecturas totales, lo cual resulta en la proporción de lecturas asignadas a los genes y facilita la identificación de muestras malas.

```{r}
# Ver el resumen de los niveles de expresion    
rse_gene$assigned_gene_prop <- rse_gene$recount_qc.gene_fc_count_all.assigned / rse_gene$recount_qc.gene_fc_count_all.total
summary(rse_gene$assigned_gene_prop)
```

Una forma útil para interpretarlo es mediante una gráfica RIN como la siguiente. En ella se puede observar que algunas muestras no son lo suficientemente buenas, por lo que se deberá considerar para la limpieza de datos.

```{r}
# Graficar los niveles de expresion RIN
with(colData(rse_gene), plot(assigned_gene_prop, sra_attribute.rin))
abline(v=0.37,col = "red")
```

Por último, analicé como varían los resultados entre los dos fenotipos usando el siguiente código. Al parecer no hay mucha diferencia entre ambas condiciones.       

```{r}
# Checar si hay una diferencia entre los grupos
with(colData(rse_gene), tapply(assigned_gene_prop, sra_attribute.phenotype, summary))
```

---

# <span style="color:orange">Limpieza de datos</span>

A continuación, voy a hacer la limpieza de las muestras poco informativas. Observa en la siguiente gráfica la distribución de las muestras.

```{r}
# Guardar nuestro objeto RSE por si luego cambio de opinión
rse_gene_unfiltered <- rse_gene

# Restablecer el objeto RSE a una instancia antes del filtrado
#rse_gene <- rse_gene_unfiltered

# Graficar la distribucion de las muestras
hist(rse_gene$assigned_gene_prop)
abline(v=0.37,col = "red")
```

Se descartan las uestras que estan por debajo de un umbral de 0.37, el cual representa el primer cuartil. La distribución resultante es la siguiente.

```{r}
# Eliminar las muestras de menor calidad
rse_gene <- rse_gene[, rse_gene$assigned_gene_prop > 0.37]
hist(rse_gene$assigned_gene_prop)
```

Es momento de hacer limpieza de genes poco informativos. Las estadísticas de todos los genes son las siguientes.

```{r}
# Obtener estadísticas de la expresión de genes
gene_means <- rowMeans(assay(rse_gene, "counts"))
summary(gene_means)
```

Se eliminará de nuevo el primer cuartil. Por último, se índica el porcentaje de muestras conservadas depués del filtrado. En este caso, me quede con aproximadamente el 73.4% de los datos originales. De nuevo muestro la gráfica RIN, pero esta vez con los datos filtrados. Al parecer, la calidad de los datos mejoró considerablemente después de la limpieza.

```{r}
# Filtrar genes
rse_gene <- rse_gene[gene_means > 0.2, ]
round(nrow(rse_gene) / nrow(rse_gene_unfiltered) * 100, 2)
```

```{r}
# Graficar los niveles de expresion RIN
with(colData(rse_gene), plot(assigned_gene_prop, sra_attribute.rin))
```

---

# <span style="color:orange">Normalización de datos</span>

Es necesario normalizar los datos asumiendo que la mayoria de los genes no se estan expresando diferencialmente. El objetivo es redudir la incidencia de falsos positivos.

```{r}
# Construyo un objeto con el cual se podran normalizar los datos
library("edgeR")
dge <- DGEList(
    counts = assay(rse_gene, "counts"),
    genes = rowData(rse_gene))
dge <- calcNormFactors(dge)
```

---

# <span style="color:orange">Expresión diferencial</span>

A continuación, se realizó un análisis mediante gráficas de tipo *boxplot*, el cual es útil para visualizar la diferencia entre la expresión de las muestras bajo distintas condiciones. Analicé la diferencia de expresión entre las muestras que tenían MDD y las que no; las muestras de mujeres y las de hombres; y entre individuos que se suicidaron y los que murieron de otras formas. Existe poca diferencia entre la expresión diferencial de cada uno de los casos, por lo que un análisis más exhaustivo es necesario.

```{r}
# Importa una libreria para generar boxplot
library("ggplot2")

# Grafica la expresion diferencial entre las condiciones MDD y no MDD
ggplot(as.data.frame(colData(rse_gene)), aes(y = assigned_gene_prop, x = sra_attribute.phenotype)) +
    geom_boxplot() +
    theme_bw(base_size = 20) +
    ylab("Assigned Gene Prop") +
    xlab("Phenotype")
```

```{r}
# Grafica la expresion diferencial entre las condiciones hombre y mujer
ggplot(as.data.frame(colData(rse_gene)), aes(y = assigned_gene_prop, x = sra_attribute.gender)) +
    geom_boxplot() +
    theme_bw(base_size = 20) +
    ylab("Assigned Gene Prop") +
    xlab("Gender")
```

```{r}
# Grafica la expresion diferencial entre las condiciones suicida y no suicida
ggplot(as.data.frame(colData(rse_gene)), aes(y = assigned_gene_prop, x = death)) +
    geom_boxplot() +
    theme_bw(base_size = 20) +
    ylab("Assigned Gene Prop") +
    xlab("Cause of death")
```

Generé el modelo estadistico de acuerdo a las variables previamente mencionadas. Además, señalo a cuál instancia de las variables apuntan cada análisis. Estas instancias son MDD, hombre y suicida, respectivamente.

```{r}
# Generar el modelo linear estadistico
mod <- model.matrix(~ sra_attribute.phenotype + sra_attribute.gender + death + assigned_gene_prop,
    data = colData(rse_gene))

# Observar las variables que componen el modelo
colnames(mod)
```

Usando un modelo de regresión lineal, es posible obtener estimados de la desviación estandar. La siguiente gráfica sintetiza las desviaciones estándar y valor de expresión logarítmica para cada gen. La forma que toma la curva roja indica que los datos son lo suficientemente robustos como para avanzar en el análisis.

```{r}
# Generar la grafica para visualizar la desviacion estandar
library("limma")
vGene <- voom(dge, mod, plot = TRUE)
```

Ahora, genero el modelo de la regresión lineal y muestro la cantidad de genes que poseen una expresión diferencial significativa. El resultado son 28 de 63,828 genes.

```{r}
# Generar el modelo de regresion lineal
eb_results <- eBayes(lmFit(vGene))

# Indicar el coeficiente del modelo, que en este caso es la variable de interes 
de_results <- topTable(
    eb_results,
    coef = 2,
    number = nrow(rse_gene),
    sort.by = "none"
)
# Mostrar los genes diferencialmente expresados entre control y MDD con FDR < 5%
table(de_results$adj.P.Val < 0.05)
```

Ahora, utilizo una gráfica con la cual puedo visualizar el cambio en los niveles de expresión entre las muestras MDD y no MDD. Los valores positivos indican que la expresión es mas alta en el caso MDD y los valores negativos indican que es más alto en el caso no MDD.

```{r}
# Visualizar los resultados estadísticos
plotMA(eb_results, coef = 2)
```

Mediante el uso de una gráfica de tipo *volcano*, se indican los genes con mayor expresión diferencial y con mejor valor de p-value.

```{r}
volcanoplot(eb_results, coef = 2, highlight = 5, names = de_results$gene_name)
```

A continuación, se muestra la información de los genes con mayor expresión diferencial. Es importante remarcar que hay dos genes llamados TRPM1, pero el significativo solo es uno.

```{r}
# Mostrar la informacion de los genes
de_results[de_results$gene_name %in% c("PAX8-AS1", "TRPM1", "FAM182A", "PAX8", "NKAPP1"), ]
```

En la siguiente gráfica de tipo *heatmap* se observa la expresión de los genes dependiendo de todas las condiciones que se analizaron en el modelo. La información solo se extrajo de los primeros 25 genes mas significativos. 

```{r}
# Extraer valores de los genes de interes
exprs_heatmap <- vGene$E[rank(de_results$adj.P.Val) <= 25, ]

# Crear una tabla con informacion de las muestras y con nombres de columnas mas amigables
df <- as.data.frame(colData(rse_gene)[, 
        c("sra_attribute.gender", "death", "sra_attribute.phenotype")])
colnames(df) <- c("Gender", "Cause of death", "Phenotype")


# Indicar el nombre de cada uno de los genes
rownames(exprs_heatmap) <- rowRanges(rse_gene)$gene_name[match(rownames(exprs_heatmap), rownames(rse_gene))]

# Importar la libreria requerida para graficar el heatmap
library("pheatmap")

# Crear el heatmap
pheatmap(
    exprs_heatmap,
    cluster_rows = TRUE,
    cluster_cols = TRUE,
    show_rownames = TRUE,
    show_colnames = FALSE,
    annotation_col = df,
    fontsize_row = 6,
)
```

Por último, grafique la expresión diferencial dependiendo de la categoria a la que pertenecen las muestras. En el caso entre mujeres y hombres existe una gran diferencia de genes diferencialmente expresados, como lo indicaban en el estudio. Esto se interpreta gracias a la formación de clusters correspondientes a cada instancia de las variables en la siguiente gráfica.

```{r}
# Importar una libreria requerida para asignar colores
library("RColorBrewer")

# Convertir los valores de genero a colores
col.gender <- df$Gender
levels(col.gender) <- brewer.pal(nlevels(col.gender), "Set1")
col.gender <- as.character(col.gender)

## MDS por genero
plotMDS(vGene$E, labels = df$Gender, col = col.gender)
```

Sin embargo, encontre un panorama mucho menos informativo al analizar en el caso entre individuos con y sin MDD, y entre individuos que fueron suicidas y los que no. Se interpreta que no hay muchos genes diferencialmente expresados ya que no se forman clusters correspondientes a cada instancia de la variable. Esto se observa en las siguientes dos gráficas.

```{r}
## Convertir los valores de MDD a colores
col.phenotype <- df$Phenotype
levels(col.phenotype) <- brewer.pal(nlevels(col.phenotype), "Dark2")
col.phenotype <- as.character(col.phenotype)

## MDS por MDD
plotMDS(vGene$E, labels = df$Phenotype, col = col.phenotype)
```

```{r}
## Conviertiendo los valores de Sex a colores
col.death <- df$`Cause of death`
levels(col.death) <- brewer.pal(nlevels(col.death), "Set2")

col.death <- as.character(col.death)

## MDS por sexo
plotMDS(vGene$E, labels = df$`Cause of death`, col = col.death)
```

---

# <span style="color:orange">Conclusiones</span>

Al igual que el estudio donde se obtuvieron los datos con los que trabaje, encontre una gran divergencia, entre hombres y mujeres, en la caracterización molecular de las firmas transcripcionales; esto ocurre a pesar de que los hombres y las mujeres con MDD exhiben una sintomatología similar y a la vez diferente. Los hallazgos de los investigadores sugieren que el MDD en hombres y mujeres puede surgir en parte debido a las acciones de módulos genéticos similares, que comparten especificidad celular y biológica, pero que se organizan y expresan de manera diferente a través de las regiones del cerebro en los dos sexos. Las alteraciones en esta estructura transcripcional precisa pueden interferir con la actividad coordinada de varias regiones del cerebro en hombres y mujeres y, en consecuencia, alterar las estrategias normalmente utilizadas para hacer frente al estrés.

Sin embargo, yo encontre genes distintos a los que los investigadores señalaron como relevantes, los cuales fueron *DUSP6* y *EMX1*. En su lugar yo encontre que los genes *PAX8-AS1*, *TRPM1*, *FAM182A*, *PAX8* y *NKAPP1* poseen niveles de expresión diferencial significativos. La razón de esta aparente discordancia, podría deberse al uso de las cuentas en lugar de usar RPKMs o CPMs.

La heterogeneidad de los trastornos depresivos es uno de los principales impedimentos para comprenderlos mejor y desarrollar mejores tratamientos. Los resultados del estudio original destacan posibles nuevas vías para desarrollar estrategias terapéuticas más específicas para el tratamiento del MDD en hombres y mujeres. Además, el análisis de más variables podría ayudar a esclarecer la penumbra que rodea a este transtorno. Entre estas variables yo incluí el suicidio, y aunque no mostró ser significativo, hay muchas más opciones; como por ejemplo, el consumo de drogas, alcohol y tabaco, la medicación y la edad, entr otras. 

---

# <span style="color:orange">Referencias</span>

[Labonté, B., Engmann, O., Purushothaman, I., Menard, C., Wang, J., Tan, C., Scarpa, J. R., Moy, G., Loh, Y.-H. E., Cahill, M., Lorsch, Z. S., Hamilton, P. J., Calipari, E. S., Hodes, G. E., Issler, O., Kronman, H., Pfau, M., Obradovic, A. L. J., Dong, Y., … Nestler, E. J. (2017). Sex-specific transcriptional signatures in human depression. Nature Medicine, 23(9), 1102-1111.](https://doi.org/10.1038/nm.4386)

[Collado-Torres L (2021). Explore and download data from the recount3 project. doi: 10.18129/B9.bioc.recount3](https://doi.org/10.18129/B9.bioc.recount3)

Robinson MD, McCarthy DJ and Smyth GK (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139-140

H. Wickham. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York, 2016.

[Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK (2015). “limma powers differential expression analyses for RNA-sequencing and microarray studies.” Nucleic Acids Research, 43(7), e47. doi: 10.1093/nar/gkv007.](https://doi.org/10.1093/nar/gkv007)

[Raivo Kolde (2019). pheatmap: Pretty Heatmaps. R package version 1.0.12.](https://CRAN.R-project.org/package=pheatmap)

[Erich Neuwirth (2014). RColorBrewer: ColorBrewer Palettes. R package version 1.1-2.](https://CRAN.R-project.org/package=RColorBrewer)