Analisis_Compacto.Rmd

---
title: "Analisis de Datos Omicos en COVID"
author: "Alex Sanchez"
date: "Enero 2022"
output:
   html_document:
    code_folding: show
    toc: true
    toc_depth: 3
    theme: darkly
    highlight: textmate
    number_sections: true
bibliography: "ADOreferences.bib"
editor_options: 
  chunk_output_type: console
# bibliography: references.bib
# link-citations: yes
# theme args should be one of: "default", "cerulean", "journal", "flatly", "darkly", "readable", "spacelab", "united", "cosmo", "lumen", "paper", "sandstone", "simplex", "yeti"
# highlight arg should be one of: "default", "tango", "pygments", "kate", "monochrome", "espresso", "zenburn", "haddock", "breezedark", "textmate"
---


```{r class.source = 'fold-hide', setup, include=FALSE}
library(knitr)
library(rmdformats)
source('Rcode/functions.R')

## Global options
options(max.print="75")
opts_chunk$set(echo=TRUE,
	             cache=TRUE,
               prompt=FALSE,
               tidy=TRUE,
               comment=NA,
               message=FALSE,
               warning=FALSE)
opts_knit$set(width=75)
```


# Introducción

Este estudio presenta un análisis básico de RNA-Seq con énfasis en las posibilidades de visualización de datos y resultados^.

En su versión actual los gráficos son correctos pero poco elaborados. Para mejorarlos se utilizaran gráficos con `ggplot2` utilizando como modelo algunas de las propuestas del documento: https://github.com/UofABioinformaticsHub/DataVisualisaton_BIS2016.git


## Infraestructura informática para el análisis

Este estudio se llevará a cabo usando R/Bioconductor por lo que es preciso tener instaladas un conjunto de librerías. Esto puede hacerse siguiendo el procedimiento descrito a continuación.

**El código que se presenta debería ejecutarse tan sólo una vez!**

```{r instalaPaquetes, librerias, echo=TRUE, eval=FALSE}
if(!require(BiocManager)){
      install.packages("BiocManager", dep=TRUE)
}

installifnot <- function (pckgName, BioC=TRUE){
  if(BioC){
    if(!require(pckgName, character.only=TRUE)){
      BiocManager::install(pckgName)
    }
  }else{
    if(!require(pckgName, character.only=TRUE)){
      install.packages(pckgName, dep=TRUE)
    }
  }
}
installifnot("limma")
installifnot("edgeR")
installifnot("org.Hs.eg.db")
installifnot("clusterProfiler")
installifnot("dplyr", BioC=FALSE)
installifnot("gplots", BioC=FALSE)
installifnot("ggvenn", BioC=FALSE)
installifnot("pheatmap")
installifnot("prettydoc", BioC=FALSE)
installifnot("ggnewscale", BioC=FALSE)
# Pels grafics que empren ggplot2
installifnot("locfit", BioC=FALSE)
installifnot("magrittr", BioC=FALSE)
installifnot("statmod", BioC=FALSE)
# Pel HeatMap Interactiu 
installifnot("heatmaply")
# Pel VolcanoPlot interactiu
installifnot("here")
installifnot("janitor") # Cleaning column names  
installifnot("scales") # Transform axis scales   
installifnot("ggrepel")
```

Puede resultar útil, aunque no es para nada imprescindible, disponer de, al menos dos, directorios específicos:
- `datos` (o un nombre similar) donde guardar y de donde cargar  los datos.
- `results`o (o un nombre similar) donde escribir los resultados.

```{r directorios}
if(!dir.exists("datos")) dir.create("datos")
if(!dir.exists("results")) dir.create("results")
if(!dir.exists("figures")) dir.create("figures")
```

# Los datos para el análisis

Este análisis se basa en un estudio publicado recientemente [@Arunachalam2020] que investigaba la respuesta inmune a la infección con SARS-COV-2 desde un a perspectiva de biología de sistemas utilizando tecnología de secuenciación, RNA-Seq.

Los datos se han depositado en el repositorio "[Gene Expression Omnibus (GEO)](https://www.ncbi.nlm.nih.gov/geo/)" con el identificador [GSE152418](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE152418).

## Lectura de la matriz de contajes

Los datos de contajes se han descargado desde l repositorio "[GEO]((https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE152418))" al archivo `Rawcounts.csv` que se utilizará como punto de partida para el análisis.

```{r}
counts <- read.csv("datos/RawCounts.csv", row.names = 1)
selectedCounts <- as.matrix(counts)
```

En la misma web se dispone de información los grupos y otras covariables o variables auxiliares. Con éstos se crea un objeto (habitualmente denominado "targets") que debe estar _sincronizado_ con el anterior, es decir, cuyas filas deben corresponderse con las columnas de la matriz de datos.

```{r creaTargets}
sampleNames <- colnames(selectedCounts)
grupos <- c(rep("COVID", 17), rep("SANO", 17))
colores=c(rep("red", 17), rep("blue", 17))
selectedTargets <- data.frame(samples=sampleNames, group=grupos, cols=colores)
rownames(selectedTargets) <- selectedTargets[,1]
```

# Preprocesado de los datos

## Estandarización de los contajes 

Los datos de secuenciación pueden estar "desbalanceados" en el sentido que distintas muestras pueden contener un número distinto de secuencias, lo que puede inducir erróneamente a pensar que un gen se expresa más en una muestra que en otra, cuando esto se deba a esta diferencia global.

Esto puede evitarse expresando los contajes como  "CPMs" es decir "counts per million", lo que no modificará los resultados de comparaciones posteriores, pero hará que las muestras sean comparables en número , lo que es útil y necesario para los análisis posteriores.

```{r getCPM}
library(edgeR)
selectedCounts[1:5,1:6]
apply(selectedCounts,2, sum)
counts.CPM <- cpm(selectedCounts)
counts.CPM[1:5,1:6]
apply(counts.CPM,2, sum)
```

Una vez los datos estan como CPMs, se procede a filtrarlos, 

## Filtraje de genes poco expresados

Los genes con recuentos muy bajos en todas las librerías (es decir en todas las muestras) proporcionan poca evidencia de expresión diferencial, por lo que es habitual eliminar aquellos genes que, o bien varían muy poco entre grupos, o bien presentan poca o nula expresión en la mayoría de las muestras.

En este caso, siguiendo un criterio habitual, se opta por _conservar únicamente aquellos genes que presentan algún valor en, al menos, tres muestras **de cada grupo**_.

```{r subsetThreshold}
thresh <- counts.CPM > 0
keep <- (rowSums(thresh[,1:10]) >= 3) &
        (rowSums(thresh[,11:20]) >= 3)
counts.keep <- counts.CPM[keep,]
dim(counts.CPM)
dim(counts.keep)
```

Aunque no sea más que un ejemplo basta con ver los dos primeros genes para comprobar como el primero no cumple la condición, en el grupo "SANO", mientras que los siguientes sí que la cumplen. Por lo tanto, al filtrar desaparece el primer gen de la matriz filtrada, pero no los dos siguientes.

```{r comparaMatrices}
round(head(counts.CPM), 3)
round(head(counts.keep), 3)
```

## Uso de clases específicas para manejar los datos

Cuando se trabaja con distintos objetos referidos a unos mismos datos, como la matriz de contajes y el objeto "targets",  es útil disponer de _contenedores_ que permitan trabajar con todos ellos a la vez, lo que no sólo facilita el trabajo sino que ayuda a evitar "desincronizaciones". 

Éste es el caso de la clase `ExpressionSet` habitualmente utilizada con microarrays o de la clase que la generaliza, llamada  `SummarizedExperiment`. 

Para datos de contaje es habitual usar una clase similar a `ExpressionSet` llamada `DGEList`" pensadas para manejar datos de contajes , definida en el paquete `edgeR`. Esta clase, más simple que las anteriores, utiliza listas para almacenar recuentos de "reads" e información asociada de tecnologías de secuenciación. Puede encontrarse información al respecto en la ayuda del paquete `edgeR`.

```{r makeDGEObj0}
dgeObj <- DGEList(counts = counts.keep, 
                  lib.size = colSums(counts.keep),
                  norm.factors = rep(1,ncol(counts.keep)), 
                  samples = selectedTargets,
                  group = selectedTargets$group, 
                  genes = rownames(counts.keep), 
                  remove.zeros = FALSE)
show(dgeObj)
```

Uno de los aspectos interesantes de estas clases es la posibilidad de extraer partes de todos los objetos a la vez con el operador de "subsetting". 

```{r makeDGEObj}
dim(dgeObj)
dgeObjShort<-dgeObj[,c(1:5, 11:15)]
# Library size information is stored in the samples slot
dgeObjShort$samples
colnames(dgeObjShort$counts)
```


```{r FRC0}
# https://hbctraining.github.io/DGE_workshop_salmon/lessons/AnnotationDbi_lesson.html
annotations_orgDb <- AnnotationDbi::select(org.Hs.eg.db, # database
                                           keys = dgeObj$genes[,1],  # data to use for retrieval
                                           columns = c("SYMBOL", "ENTREZID","GENENAME"), # information to retreive for given data
                                           keytype = "ENSEMBL")

length(which(is.na(annotations_orgDb$SYMBOL) == FALSE))
non_duplicates_idx <- which(duplicated(annotations_orgDb$SYMBOL == FALSE))
annotations_orgDb<- annotations_orgDb[non_duplicates_idx, ]
```


```{r FRC00}
#https://hbctraining.github.io/Intro-to-R-flipped/lessons/09_reordering-to-match-datasets.html
cbind(dgeObj$genes[,1],annotations_orgDb$ENSEMBL)[1:100,]
reorder_idx <- match(dgeObj$genes[,1],annotations_orgDb$ENSEMBL)
annotations_orgDb_reordered <- annotations_orgDb[reorder_idx,]  
```

```{r FRC000}
cbind(dgeObj$genes[,1],annotations_orgDb_reordered$ENSEMBL)[1:100,]
```


```{r FRC0000}
dgeObj$genes<-annotations_orgDb_reordered
```


```{r FRC00000}
dgeObj$genes[1:100,]
```

Aunque podríamos haber creado el objeto a partir de todas las muestras, y haber realizado la extracción de genes y muestras posteriormente, hemos optado por no hacerlo para facilitar el seguimiento del proceso.

## Normalización

Además de estandarizar los contajes, es importante eliminar otros sesgos de composición entre librerías. Esto puede hacerse aplicando la normalización por el método TMM que genera un conjunto de factores de normalización, tal que producto de estos factores y los tamaños de librería   (el número de secuencias de cada muestra) definen el _tamaño efectivo_ de dichas muestras, es decir el peso real que se les asignará en las comparaciones posteriores.

Aunque esto puede parecer artificial, no lo es porque la normalización tiene en cuenta otros factores, como el sesgo de composición entre librerías, que podrían hacer que los mismos valores en distintas muestras no reflejaran su importancia relativa.

La función `calcNormFactors`, de la librería `edgeR`, calcula los factores de normalización mencionados.

```{r calcNormFactors}
library(edgeR)
dgeObj_norm <- calcNormFactors(dgeObj)
head(dgeObj_norm$samples, 10)
```

Esto no modificará la matriz de contajes, pero actualizará los factores de normalización en el objeto DGEList  (sus valores predeterminados son 1).

```{r}
head(dgeObj$samples)
head(dgeObj_norm$samples)
```


Es decir, _aunque no se observen cambios en la matriz de contajes_, cuando se utilizan estos factores de normalización en algún cálculo la importancia de las distintas columnas se tendrá en cuenta.

## Transformación logarítmica

Las transformaciones anteriores buscan compensar el tamaño distinto de las librerías o la distinta composición de éstas, pero las distribuciones de los contajes en cada muestra son asimétricas.

```{r}
boxplot(dgeObj_norm$counts, 
        col=dgeObj_norm$samples$cols,
        las=2, cex.axis=0.7, main="Contajes normalizados", ylim=c(0,10000))
```



Para finalizar el preprocesado se toman logaritmo base dos de los contajes.

```{r log2count_norm}
log2count_norm <- cpm(dgeObj_norm, log=TRUE)
```

```{r eval=FALSE, echo=FALSE}
boxplot(log2count_norm, col=dgeObj_norm$samples$cols,
        las=2, cex.axis=0.7, 
        main="Contajes normalizados (log2))")
```

```{r niceBoxPlot, warning=FALSE, out.width="100%"}
source("Rcode/niceBoxPlot.R")    
gg1 <- niceBoxPlot(log2count_norm) 
ggsave(filename = 'figures/boxplot.png',gg1,bg = "transparent") # ,dpi = 300
gg1
```

También se puede hacer usando plotly.

```{r niceBoxPlotly, warning=FALSE, out.width="100%"}
niceBoxPlotly(log2count_norm)  
```

**Esta será nuestra matriz de partida para los análisis siguientes**

# Exploración de los datos

Una vez descartados los genes poco expresados y con los recuentos almacenados en un objeto `DGEList`,  podemos`proceder a realizar algunos gráficos exploratorios para determinar si los datos aparentan buena calidad y/o si presentan algun problema.

## Distribución de los contajes

Un diagrama de cajas con los datos, normalizados o no, muestra que la distribución de los contajes es muy asimétrica, lo que justifica la decisión de trabajar con los logaritmos de los datos.

La transformación logarítmica puede hacerse directamente pero es mejor usar la función `cpm`, como se ha hecho, que agrega una pequeña cantidad para evitar tomar logaritmos de cero.

```{r distriCounts1, fig.height=6, out.width="100%"}
par(mfrow=c(2,1))
rawCounts <- dgeObj_norm$counts
boxplot(rawCounts, ylab="CPM",las=2, xlab="", col = dgeObj$samples$cols, cex.axis=0.6, main="Distribución de contajes")
boxplot(log2count_norm, ylab="Log2-CPM",las=2, xlab="", col=dgeObj$samples$cols, cex.axis=0.6, main="Distribución de log(contajes)")
abline(h=median(log2count_norm), col="blue")
```

## Análisis de similaridad entre las muestras

### Distancia entre muestras

La función `dist` permite calcular una _matriz de distancias_ que contiene las comparaciones dos a dos entre todas las muestras. Por defecto se utiliza una distancia euclídea.

```{r}
sampleDists <- dist(t(log2count_norm))
round(sampleDists,1)
```

Las matrices de distancias se pueden visualizar directamente con un heatmap.

```{r heatmap1, out.width="100%"}
# par(mfrow=c(1,1))
library(factoextra)
fviz_dist(sampleDists)
```

Como puede verse _las muestras tienden a agruparse por el factor SANO/COVID, aunque una de las muestras. COV155 se separa del resto de las del grupo COVID.

### Agrupamiento jerárquico

Un agrupamiento jerárquico proporciona una representación alternativa, también basada en la matriz de distancias.

```{r uglydedndogram, out.width="100%"}
hc <-hclust(sampleDists)
plot(hc,labels = colnames(log2count_norm),main = "Agrpamiento jerárquico de las muestras", cex=0.8)
```


El siguiente es un gráfico mejor.

```{r  niceDendrogram, out.width="100%"}
source("Rcode/niceDendrogram.R")
library("ggdendro", "dendextend")
hc <- hclust(sampleDists,method = 'ward.D2')
gg2 <- niceDendrogram (hc)
ggsave('figures/dendograma.png',gg2, bg = "transparent") # , dpi = 300
gg2
```

Una de las muestras COVID parece más similar a las saludables que a las otras COVID.

### Visualización en dimensión reducida


```{r uglyMDS, out.width="100%"}
col.status <- dgeObj_norm$samples$cols
plotMDS(log2count_norm,col=col.status, main="Status", cex=0.7)
```

Con el paquete [emojifont](https://cran.r-project.org/web/packages/emojifont/vignettes/emojifont.html) se pueden usar emoticonos para distinguir los pacientes sanos de los que tienen COVID.

```{r niceMDS1, out.width="100%"}
source("Rcode/niceMDS.R")   
gg3 <- niceMDS(dgeObj_norm)  
ggplot2::ggsave('figures/MDS_emotis.png',gg3) # ,dpi=300, bg = "transparent"
gg3
```

Como puede verse, el gráfico muestra la misma agrupación "natural" y el mismo comportamiento atípico de una muestra.

# Análisis de expresión diferencial

El objetivo del análisis de expresión diferencial es seleccionar genes cuya expresión difiere entre grupos. 

## Selección de genes usando limma-Voom

La ventaja principal de esta aproximación es que permite trabajar con toda la flexibilidad de los modelos lineales para representar diseños experimentales, y, en muchos casos , aprovechar la experiencia previa del usuario en el manejo de limma.

### Matriz de diseño y de contrastes

Utilizando la variable `group` podemos definir una _matriz de diseño_ y, sobre ésta, los _contrastes_ que nos interesan.

```{r matrizDisenyo}
group = as.factor(dgeObj_norm$samples$group)
design = model.matrix(~ 0 + group)
colnames(design) = gsub("group", "", colnames(design))
row.names(design) = sampleNames
design
```

Dado que estamos interesados en las diferencias entre los grupos, necesitamos especificar qué comparaciones queremos llevar a cabo. Las comparaciones de interés se puede especificar utilizando la función `makeContrasts`. La matriz de contraste indica  qué columnas de la matriz `design` vamos a comparar. En este caso tan sólo se llevará a cabo una comparación.

```{r matrizContrastes}
cont.matrix = makeContrasts(CONTROLvsCOVID = COVID - SANO,
levels=colnames(design))
cont.matrix
```

### Transformación de los contajes

```{r voom}
voomObj <- voom(dgeObj_norm, design)
voomObj
```

### Selección de genes diferencialmente expresados

Como en el caso de los microarrays el objeto `voomObj` y las matrices de diseño y contrastes se utilizaran para ajustar un modelo y, a continuación realizar las comparaciones especificadas sobre el modelo ajustado. El proceso finaliza con la regularización del estimador del error usando la función ` eBayes`.

```{r ajusteLM}
fit <- lmFit(voomObj)
fit.cont <- contrasts.fit(fit, cont.matrix)
fit.cont <- eBayes(fit.cont)
```

### Top tables

Los resultados de un análisis de expresión diferencial se pueden extraer con la  función `topTable`. Esta función genera una tabla de resultados cuyas columnas contienen información acerca de los genes y la diferencia entre los grupos comparados. Concretamente:

```{r topTable}
toptab <- topTable(fit.cont,coef=1,sort.by="p", number=nrow(fit.cont))
head(toptab)
```
### Visualización de los resultados

#### Significación biológica vs. Significación Estadística

Para visualizar los resultados podemos usar un `volcanoPlot`:

```{r uglyvolcano, out.width="100%"}
volcanoplot(fit.cont,coef=1,highlight=100, main="COVID vs SANO")
```

Es posible también crear Volcanoplots interactivos

```{r niceInteractiveVolcano}
source("Rcode/niceInteractiveVolcano.R") 
niceInteractiveVolcano(fit.cont)

```

O volcano plots más sofisticados

```{r FRC1}
library("gridExtra")
library("plotly")
```


```{r FRC2}
diseased_vs_healthy <- data.frame(fit.cont$genes[,1], fit.cont$coefficients[,1],fit.cont$p.value[,1])
colnames(diseased_vs_healthy)<-c("geneid","foldchange","adjpval")
```


```{r FRC3, include=FALSE}
vol_plot <- diseased_vs_healthy %>%
  ggplot(aes(x = foldchange,
             y = -log10(adjpval))) + 
  geom_point() 

vol_plot + 
  geom_hline(yintercept = -log10(0.05),
             linetype = "dashed") + 
  geom_vline(xintercept = c(log2(0.5), log2(2)),
             linetype = "dashed") + xlim(-6.5, 6.5) 


px<-vol_plot+ theme(rect = element_rect(fill = "transparent"))
```

```{r FRC4}
diseased_vs_healthy <- diseased_vs_healthy %>%
  mutate(gene_type = case_when(foldchange >=2 & adjpval <= 0.05 ~ "up",
                               foldchange <= -2 & adjpval <= 0.05 ~ "down",
                               TRUE ~ "ns"))   

# Obtain gene_type counts ------------------------------------------------------           
diseased_vs_healthy %>%
  count(gene_type)
```


```{r FRC5, include=FALSE}
diseased_vs_healthy %>%
  distinct(gene_type) %>%
  pull()  
#> [1] "down" "up"   "ns"    
# Add colour, size and alpha (transparency) to volcano plot --------------------
cols <- c("up" = "#ffad73", "down" = "#26b3ff", "ns" = "grey") 
sizes <- c("up" = 3, "down" = 3, "ns" = 1) 
alphas <- c("up" = 1, "down" = 1, "ns" = 0.5)

diseased_vs_healthy %>%
  ggplot(aes(x = foldchange,
             y = -log10(adjpval),
             fill = gene_type,    
             size = gene_type,
             alpha = gene_type)) + 
  geom_point(shape = 21, # Specify shape and colour as fixed local parameters    
             colour = "black") + 
  geom_hline(yintercept = -log10(0.05),
             linetype = "dashed") + 
  geom_vline(xintercept = c(-2, 2),
             linetype = "dashed") +
  scale_fill_manual(values = cols) + # Modify point colour
  scale_size_manual(values = sizes) + # Modify point size
  scale_alpha_manual(values = alphas) + # Modify point transparency
  scale_x_continuous(breaks = c(seq(-6.5, 6.5, 1)),       
                     limits = c(-7, 7))  
```

```{r FRC5x, out.width="100%"}
p<-  ggplot(diseased_vs_healthy,aes(x = foldchange,
             y = -log10(adjpval),
             fill = gene_type,    
             size = gene_type,
             alpha = gene_type)) + 
  geom_point(shape = 21, # Specify shape and colour as fixed local parameters    
             colour = "black") + 
  geom_hline(yintercept = -log10(0.05),
             linetype = "dashed") + 
  geom_vline(xintercept = c(-2, 2),
             linetype = "dashed") +
  scale_fill_manual(values = cols) + # Modify point colour
  scale_size_manual(values = sizes) + # Modify point size
  scale_alpha_manual(values = alphas) + # Modify point transparency
  scale_x_continuous(breaks = c(seq(-6.5, 6.5, 1)),       
                     limits = c(-7, 7))  
p1<-p+theme(
         panel.background = element_rect(fill='transparent'),
         plot.background = element_rect(fill='transparent', color=NA),
         panel.grid.major = element_blank(),
         panel.grid.minor = element_blank(),
         legend.key = element_rect(colour = NA, fill = NA),
         legend.background = element_rect(fill='transparent'),
         legend.box.background = element_rect(fill='transparent'))+
  annotate("text", x=6.9, y=7.8, label= "IFI27")+
  annotate("text", x=5.8, y=11.7, label= "IGHG1") +
  annotate("text", x=5.7, y=10.5, label= "IGHG1-24") + 
  annotate("text", x=5.5, y=4.3, label= "ALAS2") + 
  annotate("text", x=5.2, y=2.9, label= "HBB") + 
  annotate("text", x=-1.0, y=11.4, label= "CACNA2D3") + 
  annotate("text", x=-1.5, y=7.4, label= "UICLM") + 
  annotate("text", x=-1.45, y=5.2, label= "HLA-DQB2") + 
  annotate("text", x=-1.4, y= 4.0, label= "HLA-DQA1") 


#display boxplot
p1
```

```{r FRC5xx}
ggsave('figures/volcanoplot1.png', p1, dpi = 300, bg='transparent')
```


```{r FRC6}
# add a grouping column; default value is "not significant"

diff_df<-diseased_vs_healthy
colnames(diff_df)<-c("geneid","foldchange","adjpval","genetype")
diff_df["group"] <- "NotSignificant"

# for our plot, we want to highlight 
# FDR < 0.05 (significance level)
# Fold Change > 1.5

# change the grouping for the entries with significance but not a large enough Fold change
diff_df[which(diff_df['adjpval'] < 0.05 & abs(diff_df['foldchange']) < 1.5 ),"group"] <- "Significant"

# change the grouping for the entries a large enough Fold change but not a low enough p value
diff_df[which(diff_df['adjpval'] > 0.05 & abs(diff_df['foldchange']) > 1.5 ),"group"] <- "FoldChange"

# change the grouping for the entries with both significance and large enough fold change
diff_df[which(diff_df['adjpval'] < 0.05 & abs(diff_df['foldchange']) > 1.5 ),"group"] <- "Significant&FoldChange"


# Find and label the top peaks..
top_peaks <- diff_df[with(diff_df, order(foldchange, adjpval)),][1:5,]
top_peaks <- rbind(top_peaks, diff_df[with(diff_df, order(-foldchange, adjpval)),][1:5,])

# Add gene labels for all of the top genes we found
# here we are creating an empty list, and filling it with entries for each row in the dataframe
# each list entry is another list with named items that will be used by Plot.ly
a <- list()
for (i in seq_len(nrow(top_peaks))) {
  m <- top_peaks[i, ]
  a[[i]] <- list(
    x = m[["foldchange"]],
    y = -log10(m[["adjpval"]]),
    text = m[["geneid"]],
    xref = "x",
    yref = "y",
    showarrow = TRUE,
    arrowhead = 0.5,
    ax = 20,
    ay = -40
  )
}

```


```{r FRC7, eval=FALSE, out.width="100%"}
# make the Plot.ly plot
p <- plot_ly(data = diff_df, x=~foldchange, y=~-log10(adjpval), text = ~geneid, mode = "markers", color = ~group, xaxis=~c(-7, 7)) %>%  layout(title ="") %>% layout(annotations = a) %>% layout(xaxis = list(range=c(-7,7) )) 
 
p
```


#### Perfiles de expresión

Con el fin de observar si existen perfiles de expresión diferenciados podemo realizar un mapa de colores con los genes más diferencialmente expresados.

Es decir, fijamos un criterio de selección de genes y retenemos aquellos componentes de la tabla de resultados que lo cumplen. Por ejemplo: Genes con un p-balor ajustado inferior a 0.001 y un `fold-change' superior a 2.

```{r deg1VOOM}
topGenesBas <- rownames(subset(toptab, (abs(logFC)> 2) & (adj.P.Val < 0.01)))
length(topGenesBas)
```

Con la matriz de expresión de los genes que verifican dicha condición se puede construir un heatmap.

```{r mapaDeColores}
library(pheatmap)
mat  <- log2count_norm[topGenesBas, ]
mat  <- mat - rowMeans(mat)
pheatmap(mat, fontsize_col = 14, fontsize_row = 1)
png("figures/pheatmap_transp.png", units="in", width=30, height=22, res=600, bg="transparent")
par(bg=NA)
pheatmap(mat, fontsize = 30, fontsize_col = 30, fontsize_row = 1, treeheight_row = 60, treeheight_col = 180)
dev.off()
```

También es posible crear un Heatmap interactivo:

```{r niceInteractiveHeatMap}
source("Rcode/niceInteractiveHeatMap.R") 
mat  <- log2count_norm[topGenesBas, ]
mat  <- mat - rowMeans(mat)
niceInteractiveHeatMap(mat)
```


Los dos grupos estan diferenciados, sobre todo en un subconjunto de genes en donde las expresiones toman signos (colores) distintos. En otro de los grupos parece que, en el grupo de individuos sanos los genes diferencialmente expresados se encuentran sobre-expresados en el grupo COVID, y apenas expresados en el grupo sanos.


## Análisis de expresión diferencial usando el paquete `edgeR`

El análisis con `edgeR` es similar al anterior (se originan en el mismo equipo de investigación) pero la modelización es distinta.

El análisis utiliza un GLM pero, en una forma que recuerda a lo que se hace con `limma`, realiza un paso adicional,en el que se calcula una estimación mejorada de la dispersión (variabilidad) de las muestras que integra las estimaciones individuales y la global mediante estimación Bayes empírica.


```{r}
y = estimateDisp(dgeObj_norm, design, robust=TRUE)
plotBCV(y)
```

Con este objeto, que añade a los contajes normalizados los estimadores mejorados de dispersión, se ajusta un modelo lineal generalizado con distribución binomial para los errores.

```{r}
fit <- glmQLFit(y, design, robust = TRUE)
```

Una vez ajustado el modelo se procede a construir el contraste y realizar el test.

De hecho podemos usar la matriz de contrastes que construímos para limma voom, en la misma forma que hemos reutilizado la de diseño.

```{r}
res <- glmQLFTest(fit, contrast = cont.matrix)
```

Los resultados se almacenan en un objeto similar a la ' topTable' de `limma`.

```{r}
topTags_edge <- topTags(res, n=dim(log2count_norm)[1]) # todos los genes
head(topTags_edge)
```
 Podemos seleccionar los genes más diferencialmente expresados de la misma forma que hicimos con limma-voom
 
```{r}
topGenes_edge <- rownames(subset(topTags_edge$table, (abs(logFC)> 2) & (FDR < 0.01)))
length(topGenes_edge)
```
 
Obsérvese que la lista de genes seleccionados es muy similar en ambos casos.

```{r comparaVoomEdgeR}
library(ggvenn)
x = list(LimmaVoom = topGenesBas, edgeR = topGenes_edge)
ggvenn(x, fill_color = c("#0073C2FF", "#EFC000FF"), stroke_size = 0.5, set_name_size = 3)
```

# Anotación de resultados y análisis de significación biológica

Para el análisis de significación se utilizan dos listas de transcritos:

- La lista de transcritos diferencialmenete expresados
- La lista de todos los tránscritos o "Universo"


```{r}
topGenes <- union(topGenesBas, topGenes_edge)
length(topGenes)
universe <- rownames(toptab)
length(universe)
```

## Anotación de los identificadores

Esto es posible, y de hecho sencillo de llevar a cabo, usando el paquete `annotate`.

```{r anotaTop}
library(org.Hs.eg.db)
AnnotationDbi::keytypes(org.Hs.eg.db)
topAnots = AnnotationDbi::select(org.Hs.eg.db, topGenes, c("SYMBOL", "ENTREZID", "GENENAME"),
keytype = "ENSEMBL")
head(topAnots)
dim(topAnots)
```

Como puede verse, el número de anotaciones es el mismo que el de identificadores ENSEMBL, lo que podría llevar a pensar que, es posible que, antes de subir los datos a GEO se hayan agrupado los contajes por genes.

Para la anotación del universo se procederá igual.

```{r anots}
univAnots = AnnotationDbi::select(org.Hs.eg.db, universe, c("SYMBOL", "ENTREZID", "GENENAME"), keytype = "ENSEMBL")
head(univAnots)
dim(univAnots)
```

En este caso se observa como hay más anotaciones que transcritos, lo que sugiere que múltiples transcritos han sido mapeados en el mismo gen.

## Análisis de enriquecimiento

El paquete `clusterProfiler` admite identificadores de tipo ENSEMBL y permite gran variedad de análisis complementarios al enriquecimiento por lo que, es una de las mejores opciones para el análisis de significación biológica.

```{r enrichment, out.width="100%"}
library(clusterProfiler)
library(org.Hs.eg.db)
ego = enrichGO(gene = topGenes, 
               universe=universe,
               keyType = "ENSEMBL", 
               OrgDb = org.Hs.eg.db,
               ont="BP",
               pAdjustMethod = "BH",
               pvalueCutoff = 0.05,
               qvalueCutoff = 0.05,
               readable = TRUE)
head(ego[,-8], n=5)
```

Con los resultados del análisis de enriquecimiento se pueden llevar a cabo distintas visualizaciones cuya interpretación exacta puede verse en el manual de clusterProfiler


```{r viewEnrichment1, out.width="100%"}
dotplot(ego, showCategory=7)
```

Un dotplot un poco más trabajado

```{r viewEnrichment2, out.width="100%"}
## sense transparencia
png("figures/dotplot.png", units="in", width=14, height=9, res=600)
dotplot(ego, showCategory=15, font.size = 12)
invisible(dev.off())
## amb transparencia
ego2 = clusterProfiler::simplify(ego, cutoff = 0.5, by = "p.adjust")
p <- dotplot(ego2, showCategory=15, font.size = 20) + 
  scale_colour_gradient2(low = "red", mid = "orange", high = "orange", midpoint = 0.0010) +
  theme(panel.background  = element_rect(fill = "transparent"),             # bg of the panel
        plot.background   = element_rect(fill = "transparent", color = NA), # bg of the plot
        legend.background = element_rect(fill = "transparent"),             # get rid of legend bg
        legend.box.background = element_rect(fill = "transparent"))         # bg legend boxes
ggsave("figures/dotplot_transp.png", width=14, height=9, dpi=300, bg = "transparent")
p
```

```{r viewEnrichment3, out.width="100%"}
library(clusterProfiler)
library(ggplot2)
ego2 = clusterProfiler::simplify(ego, cutoff = 0.5, by = "p.adjust")
png("figures/cnetplot_transp.png", units="in", width=24, height=16, res=600, bg="transparent")
par(bg=NA)
a <- cnetplot(ego2, showCategory=5, cex_category=1, cex_label_category=2.5, 
         cex_gene=1, cex_label_gene=1, circular=FALSE, colorEdge=TRUE)
a
invisible(dev.off())
a
```

```{r viewEnrichment4, out.width="100%"}
library(enrichplot)
goplot(ego2, showCategory=6, cex=0.1)
```

```{r viewEnrichment5, out.width="100%"}
heatplot(ego2)
```


```{r viewEnrichment6, out.width="100%"}
term_similarity_matrix = pairwise_termsim(ego)
emapplot(term_similarity_matrix, showCategory = 15, group_category=TRUE, group_legend=TRUE)
```

# Referencias cortas

- @geoquery
- @oligo
- @pdmogene
- @AQM
- @ggplot2
- @ggrepel
- @pvca
- @rmaIri
- @Hackstadt2009
- @Jeanmougin2010
- @Smyth2005
- @Smyth2004
- @Tusher2001
- @Chrominski2015
- @Tan2011
- @Allison2006
- @GregoryAlvord2007
- @Benjamini1995
- @biobase
- @genefilter
- @xtable
- @Sanchez-Pla2014
- @mogene21
- @anotate
- @orgMm
- @gplots
- @Draghici2012
- @Efron2013
- @Khatri2005
- @Khatri2012
- @Arunachalam2020


# Referencias largas