1 Introducción y Objetivos

Entre las infecciones nosocomiales más frecuentes, se encuentran aquellas por Staphylococcus aureus. Existen cepas resistentes a meticilina (MRSA) cuya patogénesis parece estar relacionada con la producción de determinadas toxinas.

Se ha observado que la interrupción de la síntesis de toxina por inhibidores de síntesis de proteínas, podría influir en la mejoría clínica de los pacientes. Linezolid es un potente inhibidor de la síntesis proteica por lo que sería interesante estudiar su efecto en infecciones mediadas por toxinas.

Sin embargo, la vancomicina es un agente activo de la pared celular, los cuales suelen inducir la producción de toxinas bacterianas.

Estos dos tipos de antibióticos pueden tener efectos muy distintos en la expresión de genes de toxinas bacterianos. Los objetivos consistirán en realizar las comparaciones de los perfiles de expresión génica entre:

  • Ratones infectados y no infectados sin tratamiento.

  • Ratones infectados y no infectados tratados con linezolid.

  • Ratones infectados y no infectados tratados con vancomicina.

Aunque este análisis se ha realizado en el artículo tanto en el huésped (para lo cual se ha empleado un modelo murino) como en el patógeno, en este informe nos centraremos en el estudio de expresión génica en el huésped (ratón).

2 Métodos

Para el análisis de los perfiles de expresión génica, se han tomado muestras de ARN de sangre de ratones, resultando en un total de 35 muestras con calidad suficiente para ser usadas en análisis de microarray. Se han utilizado arrays de Affymetrix GeneChip® Mouse Genome 430 2.0, los cuales constan de 45,101 sondas.

Tras realizar un análisis exploratorio y de calidad de los datos, se utilizará limma para seleccionar los genes diferencialmente expresados (GDE) para cada una de las comparaciones. A continuación, se realizará la anotación de estos genes y se procederá con un análisis de enriquecimiento (ORA), para identificar las categorías enriquecidas.

3 Datos para el análisis

En primer lugar, descargaremos los archivos CEL con los datos crudos (GSE38531_RAW.tar) a partir del siguiente enlace: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE38531

De las 35 muestras, seleccionaremos un subconjunto de 24 muestras de forma aleatoria. Estos datos crudos serán preprocesados usando robust multiarray average analysis (RMA), que realiza una corrección del ruido de fondo, una compresión (summarization) y una transformación logarítmica. Además, nos proporciona un objeto ‘ExpressionSet’ con el que trabajaremos.

4 Análisis exploratorio y control de calidad

Podemos observar los datos crudos de expresión y los datos tras el pre-procesado con RMA sin normalizar.

Biobase::exprs(gse_affy)[1:5, 1:5] # Datos crudos
##   GSM944864 GSM944836 GSM944843 GSM944857 GSM944840
## 1       204       196       212       365       233
## 2     17994     14487     16280     13342     10731
## 3       572       426       452       567       386
## 4     17551     14813     16042     13277     11168
## 5       231       326       260       586       355
Biobase::exprs(rmaData)[1:5, 1:5] # Datos pre-procesados
##              GSM944864 GSM944836 GSM944843 GSM944857 GSM944840
## 1415670_at    7.107638  6.812403  7.152999  7.107432  7.126219
## 1415671_at   10.469632 10.433229 10.700629 10.481766 10.088008
## 1415672_at   11.038996 11.153967 11.385511 11.111767 11.095912
## 1415673_at    6.708517  6.518159  6.594021  6.105470  6.714843
## 1415674_a_at  8.615501  8.440357  8.525299  8.331207  8.001695

Primero, nos interesa visualizar mediante un gráfico PCA si las muestras se encuentran agrupadas adecuadamente. En principio, no podríamos analizar un posible efecto batch, ya que no contamos con réplicas técnicas y desconocemos si existen lotes o qué muestras pertenecen a ellos.

Los gráficos PCA de los datos preprocesados indican que el primer componente principal diferencia de acuerdo a la infección, lo que significa que juega un papel importante en las diferencias de expresión génica.







También podemos representar la intensidad de las sondas con un gráfico boxplot. Al observar el boxplot de los datos preprocesados sin normalizar, podemos ver que aunque las distribuciones de las muestras difieren ligeramente entre ellas, podrían beneficiarse de una normalización.


## Background correcting
## Normalizing
## Calculating Expression




Podemos realizar otro análisis de calidad, llamado Relative Log Expression (RLE), que nos da una idea de cómo la expresión de los transcritos de una muestra varía con respecto a la de otras muestras. En el boxplot, no observamos que ninguno de los arrays difiera significativamente del resto, por lo que no habría outliers.

5 Filtrado de los datos

Filtramos las sondas y seleccionamos aquellas con mayor variabilidad en la expresión entre las muestras. De 45101 sondas que contenía nuestro ExpresionSet inicial, pasamos a 2048 sondas.

# Filtramos features cuyo IQR esté por debajo del percentil 90
(normData_filtered <- nsFilter(normData, var.func = IQR,
                          var.cutoff = 0.9, var.filter = TRUE,
                          require.entrez = TRUE,
                          filterByQuantile = TRUE))
## $eset
## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 2048 features, 24 samples 
##   element names: exprs 
## protocolData
##   sampleNames: GSM944864 GSM944836 ... GSM944841 (24 total)
##   varLabels: ScanDate
##   varMetadata: labelDescription
## phenoData
##   rowNames: GSM944864 GSM944836 ... GSM944841 (24 total)
##   varLabels: sample infection time agent
##   varMetadata: labelDescription
## featureData: none
## experimentData: use 'experimentData(object)'
## Annotation: mouse4302 
## 
## $filter.log
## $filter.log$numDupsRemoved
## [1] 16955
## 
## $filter.log$numLowVar
## [1] 18428
## 
## $filter.log$numRemoved.ENTREZID
## [1] 7657
## 
## $filter.log$feature.exclude
## [1] 13

6 Construcción de la matriz de diseño y de contrastes

Construimos un modelo lineal, con una matriz de diseño y una de constrastes, para el análisis. Ya que nuestro objetivo es comparar infectados y no infectados para los tres tratamientos, crearemos un factor combinado con la infección y el agente. Usando el paquete limma ajustaremos un modelo lineal a los datos normalizados utilizando la matriz de diseño. En segundo lugar, aplicaremos los contrastes de hipótesis para realizar las comparaciones y calcularemos la significancia para cada uno de los genes. El resultado del análisis se guardará en un objeto MArrayLM, a partir del cual construiremos las top tables, que nos mostrarán los GDE para cada contraste.

6.0.0.1 Matriz de diseño

##           Linezolid_Inf Linezolid_Uninf Untreated_Inf Untreated_Uninf
## GSM944864             1               0             0               0
## GSM944836             1               0             0               0
## GSM944843             1               0             0               0
## GSM944857             1               0             0               0
## GSM944840             0               1             0               0
## GSM944847             0               1             0               0
##           Vancomycin_Inf Vancomycin_Uninf
## GSM944864              0                0
## GSM944836              0                0
## GSM944843              0                0
## GSM944857              0                0
## GSM944840              0                0
## GSM944847              0                0

6.0.0.2 Matriz de contrastes

##                   Contrasts
## Levels             Untreated_Inf_vs_Uninf Linezolid_Inf_vs_Uninf
##   Linezolid_Inf                         0                      1
##   Linezolid_Uninf                       0                     -1
##   Untreated_Inf                         1                      0
##   Untreated_Uninf                      -1                      0
##   Vancomycin_Inf                        0                      0
##   Vancomycin_Uninf                      0                      0
##                   Contrasts
## Levels             Vancomycin_Inf_vs_Uninf
##   Linezolid_Inf                          0
##   Linezolid_Uninf                        0
##   Untreated_Inf                          0
##   Untreated_Uninf                        0
##   Vancomycin_Inf                         1
##   Vancomycin_Uninf                      -1

7 Resultados

7.1 Anotación de los genes

Procedemos a anotar los GDE de las top tables. Para ello, recuperamos los identificadores de HUGO y ENTREZ correspondientes a las sondas de cada top table.

7.1.1 Top table anotada Sin tratamiento

head(topTabAnot_Untreated, 10)[1:7]
##         PROBEID ENTREZID SYMBOL    logFC   AveExpr        t      P.Value
## 1    1421262_at    16891   Lipg 7.827995  7.306241 21.93648 4.502151e-16
## 2    1422953_at    14289   Fpr2 4.138060 11.347255 20.34523 2.078065e-15
## 3  1427747_a_at    16819   Lcn2 6.444372 10.704351 19.17606 6.862275e-15
## 4    1440865_at   213002 Ifitm6 4.495230 11.164688 18.32216 1.712446e-14
## 5    1448213_at    16952  Anxa1 3.683435 12.082075 16.83065 9.296367e-14
## 6    1449366_at    17394   Mmp8 4.913963 10.242745 16.48755 1.397413e-13
## 7    1418722_at    18054    Ngp 6.004203 11.044956 16.23213 1.901841e-13
## 8    1421366_at    23845 Clec5a 3.625643  6.373060 15.22615 6.675293e-13
## 9    1417290_at    76905   Lrg1 5.114436 10.157652 15.15640 7.301231e-13
## 10   1437060_at   380924  Olfm4 5.659287  7.705889 14.28373 2.309337e-12

7.1.2 Top table anotada Tratamiento con Linezolid

head(topTabAnot_Linezolid, 10)[1:7]
##         PROBEID ENTREZID SYMBOL    logFC   AveExpr        t      P.Value
## 1    1421262_at    16891   Lipg 6.323699  7.306241 17.72097 3.335812e-14
## 2    1440865_at   213002 Ifitm6 3.528667 11.164688 14.38254 2.021237e-12
## 3  1427747_a_at    16819   Lcn2 4.796674 10.704351 14.27312 2.342720e-12
## 4    1448562_at    22271   Upp1 4.460398  7.874669 11.55319 1.274572e-10
## 5  1419681_a_at    50501  Prok2 4.965052  7.537815 11.48719 1.416495e-10
## 6    1449366_at    17394   Mmp8 3.419552 10.242745 11.47344 1.448097e-10
## 7    1418722_at    18054    Ngp 4.232492 11.044956 11.44238 1.522168e-10
## 8    1417290_at    76905   Lrg1 3.805457 10.157652 11.27730 1.987671e-10
## 9    1427102_at    20558  Slfn4 3.728931 10.744466 11.05364 2.866383e-10
## 10   1424254_at    68713 Ifitm1 3.754836 11.367332 10.42896 8.201048e-10



7.1.3 Top table anotada Tratamiento con Vancomicina

head(topTabAnot_Vancomycin, 10)[1:7]
##         PROBEID ENTREZID SYMBOL    logFC   AveExpr        t      P.Value
## 1    1421262_at    16891   Lipg 6.948247  7.306241 19.47115 5.045338e-15
## 2  1427747_a_at    16819   Lcn2 5.442648 10.704351 16.19530 1.988939e-13
## 3    1440865_at   213002 Ifitm6 3.891800 11.164688 15.86263 2.992505e-13
## 4    1422953_at    14289   Fpr2 3.004137 11.347255 14.77017 1.206972e-12
## 5    1418722_at    18054    Ngp 5.459352 11.044956 14.75915 1.224613e-12
## 6    1449366_at    17394   Mmp8 4.041010 10.242745 13.55858 6.287812e-12
## 7    1424509_at    68891  Cd177 4.756664  9.001514 12.84116 1.768674e-11
## 8  1419681_a_at    50501  Prok2 5.493830  7.537815 12.71058 2.145437e-11
## 9    1437060_at   380924  Olfm4 4.969610  7.705889 12.54303 2.754989e-11
## 10   1434848_at    14761  Gpr27 3.999830  6.007065 12.21805 4.508101e-11

7.2 Visualización de resultados

Para visualizar los resultados de cada comparación, usaremos un volcano plot, donde podremos ver los GDE en la parte superior izquierda (genes sobreexpresados).





Con la función de limma decideTests podemos observar el número de GDE para cada comparación.

##        Untreated_Inf_vs_Uninf Linezolid_Inf_vs_Uninf Vancomycin_Inf_vs_Uninf
## Down                       12                      0                       3
## NotSig                   1961                   1983                    2016
## Up                         75                     65                      29

El diagrama de Venn nos muestra cuántos genes, de los diferencialmente expresados, son comunes entre comparaciones. En este caso, habiendo elegido un log-Fold-Change de 3 y un p-valor de 0.05, obtenemos 87 GDE para los individuos sin tratamiento, 65 para los tratados con linezolid y 32 para los tratados con vancomicina.

7.3 Análisis de la significación biológica

Obtenemos los identificadores de ENTREZ correspondientes a todas las sondas y procedemos a realizar un análisis de enriquecimiento (ORA) a partir de la top table anotada para cada contraste, basado en la Gene Ontology. Nos centraremos únicamente en los genes sobreexpresados ya que hemos seleccionado aquellos con un log-Fold-Change mayor a 3.

Además, con un dotplot de los 10 términos más enriquecidos, podemos comparar visualmente cómo de enriquecidas están las categorías y cuál es su p-valor en el test de enriquecimiento.

7.3.1 Contraste 1: Infectados vs No Infectados Sin tratamiento

Las categorías enriquecidas para este contraste están relacionadas con respuesta a bacterias, a estímulos bióticos o estímulos externos, con interacciones interespecíficas, respuesta inmunitaria, etc.

##       GOBPID       Pvalue OddsRatio  ExpCount Count Size
## 1 GO:0009617 3.510389e-11  5.672410  8.777275    31  203
## 2 GO:0043207 5.987801e-10  4.191347 17.078934    42  395
## 3 GO:0051707 5.987801e-10  4.191347 17.078934    42  395
## 4 GO:0009607 9.933517e-10  4.105343 17.338361    42  401
## 5 GO:0044419 1.146636e-09  4.047745 18.159879    43  420
## 6 GO:0009605 8.062462e-09  3.653683 22.440422    47  519
##                                                                        Term
## 1                                                     response to bacterium
## 2                                      response to external biotic stimulus
## 3                                                response to other organism
## 4                                               response to biotic stimulus
## 5 biological process involved in interspecies interaction between organisms
## 6                                             response to external stimulus
dotplot(ego_unt, showCategory=9)

7.3.2 Contraste 2: Infectados vs No Infectados con Linezolid

Las categorías enriquecidas para este contraste están relacionadas con la regulación negativa de diversas funciones moleculares, como actividad hidrolasa, catalítica, proteolisis, peptidasa, etc.

##       GOBPID       Pvalue OddsRatio ExpCount Count Size
## 1 GO:0051346 2.342173e-13 17.118056 1.830065    17   56
## 2 GO:0010466 6.245763e-10 14.325758 1.503268    13   46
## 3 GO:0043086 6.846544e-10  9.243251 2.875817    17   88
## 4 GO:0045861 2.053610e-09 11.207161 1.960784    14   60
## 5 GO:0044092 2.015733e-08  6.321665 4.477124    19  137
## 6 GO:0051336 1.828139e-07  5.955247 4.084967    17  125
##                                        Term
## 1 negative regulation of hydrolase activity
## 2 negative regulation of peptidase activity
## 3 negative regulation of catalytic activity
## 4        negative regulation of proteolysis
## 5 negative regulation of molecular function
## 6          regulation of hydrolase activity
dotplot(ego_lin, showCategory=9)

7.3.3 Contraste 3: Infectados vs No Infectados con Vancomicina

Las categorías enriquecidas para este contraste son similares a las del primer contraste e igualmente incluyen categorías relacionadas con respuesta a bacterias, a estímulos bióticos o estímulos externos, con interacciones interespecíficas, etc. Igualmente, este resultado se visualiza en el diagrama de Venn, ya que todos los GDE para esta comparación coinciden con los de la comparación sin tratamiento.

##       GOBPID       Pvalue OddsRatio ExpCount Count Size
## 1 GO:0009617 5.388317e-06  6.363158 3.265963    13  203
## 2 GO:0044419 1.220485e-05  4.973348 6.757164    18  420
## 3 GO:0043207 2.499982e-05  4.734215 6.354952    17  395
## 4 GO:0051707 2.499982e-05  4.734215 6.354952    17  395
## 5 GO:0009607 3.082250e-05  4.642535 6.451483    17  401
## 6 GO:0009605 6.227269e-05  4.258923 8.349925    19  519
##                                                                        Term
## 1                                                     response to bacterium
## 2 biological process involved in interspecies interaction between organisms
## 3                                      response to external biotic stimulus
## 4                                                response to other organism
## 5                                               response to biotic stimulus
## 6                                             response to external stimulus
dotplot(ego_van, showCategory=9)

8 Discusión

De acuerdo al análisis de expresión diferencial, las categorías enriquecidas para el contraste entre individuos no tratados y aquellos tratados con vancomicina son similares. Esto es coherente con la expectativa de que, al comparar individuos infectados con no infectados, se sobreexpresen genes asociados con la respuesta a bacterias o a estímulos externos.

Sin embargo, destaca el contraste correspondiente a los individuos tratados con linezolid, donde se observan categorías enriquecidas distintas entre infectados y no infectados.

Estas categorías son consistentes con el hecho de que linezolid actúa como un potente inhibidor de la síntesis proteica. Por tanto, en individuos infectados se sobreexpresan genes relacionados con la regulación negativa de hidrolasas, peptidasas, actividad catalítica, etc. Esto sugiere que el efecto de linezolid se manifiesta en presencia de una infección, influenciando así la inmunomodulación de la respuesta.

Al analizar las top tables anotadas, se identifican algunos genes que coinciden, como Lipg, Lcn2, Ifitm6, Mmp8, entre otros.

Lipg codifica una proteína con actividad fosfolipasa, implicada en el metabolismo de lipoproteínas. Sin embargo, este resultado no concuerda con lo reportado por Sharma-Kuinkel et al. (2013), quienes señalaron que el pathway del metabolismo de glicerolípidos, en el que Lipg está involucrado, fue observado con significancia exclusivamente en los tratamientos con linezolid.

Por otro lado, Lcn2 codifica una proteína perteneciente a la familia de las lipocalinas, que suele intervenir en procesos de inmunidad innata limitando el crecimiento bacteriano, por lo que también resultaría coherente que esté más sobreexpresada en infectados que no infectados.

9 Referencias

10 Apéndice

10.1 Preparación de los datos

## Leemos el archivo con todas las muestras
allTargets <- read.csv("allTargets.txt", sep = " ")

# Seleccionamos un subconjunto de las muestras
targets <- filter_microarray(allTargets, seed = 79037158)
rownames(targets) <- targets$sample

# Listamos los archivos .CEL
CELfiles <- list.files(path_datos, pattern = "CEL")

# Guardamos los nombres de los archivos .CEL del subconjunto
lista_CELfiles_filtrados <- character(0)

for (target in targets$sample){
  
  CELfiles_filtrados <- CELfiles[grepl(target, CELfiles)]
  
  lista_CELfiles_filtrados <- c(lista_CELfiles_filtrados, CELfiles_filtrados)
}

lista_CELfiles_filtrados <- file.path(path_datos, lista_CELfiles_filtrados)

# Construimos un objeto 'AffyBatch' 
gse_affy <- ReadAffy(filenames = lista_CELfiles_filtrados)

# Pre-procesamos los datos crudos con RMA y convertimos el objeto en un 'ExpressionSet'
rmaData <- affy::rma(gse_affy, target = "core", normalize = FALSE)

# Acortamos los nombres de las muestras
nombres_muestras <- sampleNames(gse_affy)

nombres_nuevos <- sub("_.*", "", nombres_muestras)

sampleNames(rmaData) <- nombres_nuevos
sampleNames(gse_affy) <- nombres_nuevos

# Creamos el AnnotatedDataFrame con la información sobre las muestras 
pheno_data <- AnnotatedDataFrame(data = targets)

# Lo añadimos al 'ExpressionSet'
phenoData(rmaData) <- pheno_data

10.2 Análisis exploratorio y control de calidad

10.2.1 PCA plots

# Infectados vs no infectados sin tratamiento
# Seleccionamos individuos sin tratamiento
untreated <- pData(rmaData)$agent == "untreated"

rmaData_untreated <- rmaData[, untreated]

# Realizamos Análisis Componentes Principales
PCA_untreated <- prcomp(t(exprs(rmaData_untreated)), scale. = FALSE)

# Calculamos porcentaje explicado por cada componente
porcVar <- round(100*PCA_untreated$sdev^2/sum(PCA_untreated$sdev^2),1)

comp_untreated <- data.frame(PC1 = PCA_untreated$x[,1], PC2 = PCA_untreated$x[,2],
                     Infection = pData(rmaData)$infection)

# Representamos el PCA
options(repr.plot.width = 5, repr.plot.height =2)

ggplot(comp_untreated, aes(PC1, PC2)) +
      geom_point(aes(colour = Infection)) +
  ggtitle("Inf vs No inf sin tratamiento") +
  xlab(paste0("PC1, VarExp: ", porcVar[1], "%")) +
  ylab(paste0("PC2, VarExp: ", porcVar[2], "%")) +
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_color_manual(values = c("darkorange", "dodgerblue2"))

# Infectados vs no infectados tratados con LINEZOLID
# Seleccionamos individuos tratados con Linezolid
linezolid <- pData(rmaData)$agent == "linezolid"

rmaData_linezolid <- rmaData[, linezolid]

# Realizamos Análisis Componentes Principales
PCA_linezolid <- prcomp(t(exprs(rmaData_linezolid)), scale. = FALSE)

# Calculamos porcentaje explicado por cada componente
porcVar <- round(100*PCA_linezolid$sdev^2/sum(PCA_linezolid$sdev^2),1)

comp_linezolid <- data.frame(PC1 = PCA_linezolid$x[,1], PC2 = PCA_linezolid$x[,2],
                     Infection = pData(rmaData)$infection)

# Representamos el PCA
ggplot(comp_linezolid, aes(PC1, PC2)) +
      geom_point(aes(colour = Infection)) +
  ggtitle("Inf vs No inf tratados con LINEZOLID") +
  xlab(paste0("PC1, VarExp: ", porcVar[1], "%")) +
  ylab(paste0("PC2, VarExp: ", porcVar[2], "%")) +
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_color_manual(values = c("darkorange", "dodgerblue2"))

# Infectados vs no infectados tratados con VANCOMICINA
# Seleccionamos individuos tratados con Vancomicina
vancomicina <- pData(rmaData)$agent == "vancomycin"

rmaData_vancomicina <- rmaData[, vancomicina]

# Realizamos Análisis Componentes Principales
PCA_vancomicina <- prcomp(t(exprs(rmaData_vancomicina)), scale. = FALSE)

# Calculamos porcentaje explicado por cada componente
porcVar <- round(100*PCA_vancomicina$sdev^2/sum(PCA_vancomicina$sdev^2),1)

comp_vancomicina <- data.frame(PC1 = PCA_vancomicina$x[,1], PC2 = PCA_vancomicina$x[,2],
                     Infection = pData(rmaData)$infection)

# Representamos el PCA
ggplot(comp_vancomicina, aes(PC1, PC2)) +
      geom_point(aes(colour = Infection)) +
  ggtitle("Inf vs No inf tratados con VANCOMICINA") +
  xlab(paste0("PC1, VarExp: ", porcVar[1], "%")) +
  ylab(paste0("PC2, VarExp: ", porcVar[2], "%")) +
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_color_manual(values = c("darkorange", "dodgerblue2"))

10.2.2 Boxplots

# Antes de normalizar
colores_pastel <- rainbow_hcl(ncol(exprs(rmaData)), start = 90, end = 360, l = 85, c = 45)

boxplot(rmaData, las=2, cex.axis=0.6, col = colores_pastel, 
        main = "Boxplot de los datos no normalizados")

# Normalizamos con RMA
normData <- affy::rma(gse_affy, target = "core", normalize = TRUE)
sampleNames(normData) <- nombres_nuevos

# Añadimos phenoData 
phenoData(normData) <- pheno_data

# Después de normalizar
boxplot(normData, las=2, cex.axis=0.6, col = colores_pastel, 
        main = "Boxplot de los datos normalizados")

10.2.3 RLE

# Calculamos y representamos los RLE
row_medians_assayData <- 
  Biobase::rowMedians(as.matrix(Biobase::exprs(rmaData)))

RLE_data <- as.data.frame(sweep(Biobase::exprs(rmaData), 1, row_medians_assayData))

RLE_data_gathered <- 
  tidyr::gather(RLE_data, targets, log2_expression_deviation)

ggplot2::ggplot(RLE_data_gathered, aes(targets,
                                       log2_expression_deviation)) +
  geom_boxplot(outlier.shape = NA) + 
  ylim(c(-2, 2)) + theme(axis.text.x = element_text(angle = 60, size = 6.5, hjust = 1 , 
                                                    face = "bold"))

10.2.4 ArrayQualityMetrics

output_dir <- paste0(path_PEC2,"arrayQualityMetrics") 

arrayQualityMetrics(expressionset = rmaData,
                    outdir = output_dir,
    force = TRUE, do.logtransform = TRUE,
    intgroup = c("infection", "time", "agent"))

10.3 Matriz de diseño y contrastes

agent <- pData(normData_filtered$eset)$agent
infection <- pData(normData_filtered$eset)$infection
               
# Creamos un factor combinado
grupo <- factor(paste(agent, infection, sep = "_"))  

# Creamos la matriz de diseño
design <- model.matrix(~ 0 + grupo)

colnames(design) <- c("Linezolid_Inf", "Linezolid_Uninf", "Untreated_Inf", "Untreated_Uninf", 
                      "Vancomycin_Inf", "Vancomycin_Uninf")
rownames(design) <- sampleNames(normData)

head(design)

# Creamos la matriz de constrastes
(cont.matrix <- makeContrasts(
  Untreated_Inf_vs_Uninf = Untreated_Inf - Untreated_Uninf,
  Linezolid_Inf_vs_Uninf = Linezolid_Inf - Linezolid_Uninf,
  Vancomycin_Inf_vs_Uninf = Vancomycin_Inf - Vancomycin_Uninf,
  levels = design
))

# Estimamos el modelo 
fit <- lmFit(normData_filtered$eset, design)
fit.main <- contrasts.fit(fit, cont.matrix)
fit.main <- eBayes(fit.main)

10.4 Top tables

topTab_Untreated <- topTable(fit.main, number=nrow(fit.main), coef="Untreated_Inf_vs_Uninf", 
                              adjust="fdr", p.value=0.05, lfc=3)

dim(topTab_Untreated)

topTab_Linezolid <- topTable(fit.main, number=nrow(fit.main), coef="Linezolid_Inf_vs_Uninf", 
                               adjust="fdr", p.value=0.05, lfc=3)

dim(topTab_Linezolid)

topTab_Vancomycin <- topTable(fit.main, number=nrow(fit.main), coef="Vancomycin_Inf_vs_Uninf", 
                               adjust="fdr", p.value=0.05, lfc=3)

dim(topTab_Vancomycin)

10.5 Anotación de los genes

# Anotamos las sondas 
anot_Untreated <- AnnotationDbi::select(mouse4302.db, keys = rownames(topTab_Untreated), 
                                        columns = c("ENTREZID", "SYMBOL"))

# Añadimos las anotaciones a la top table
topTabAnot_Untreated <- topTab_Untreated %>%
    mutate(PROBEID = rownames(topTab_Untreated)) %>%
    left_join(anot_Untreated) %>%
    arrange(P.Value) %>%
    select(7, 8, 9, 1:5)

anot_Linezolid <- AnnotationDbi::select(mouse4302.db, keys = rownames(topTab_Linezolid), 
                                        columns = c("ENTREZID", "SYMBOL"))

topTabAnot_Linezolid <- topTab_Linezolid %>%
    mutate(PROBEID = rownames(topTab_Linezolid)) %>%
    left_join(anot_Linezolid) %>%
    arrange(P.Value) %>%
    select(7, 8, 9, 1:5)

anot_Vancomycin <- AnnotationDbi::select(mouse4302.db, keys = rownames(topTab_Vancomycin), 
                                         columns = c("ENTREZID", "SYMBOL"))

topTabAnot_Vancomycin <- topTab_Vancomycin %>%
    mutate(PROBEID = rownames(topTab_Vancomycin)) %>%
    left_join(anot_Vancomycin) %>%
    arrange(P.Value) %>%
    select(7, 8, 9, 1:5)

head(topTabAnot_Untreated)

head(topTabAnot_Linezolid)

head(topTabAnot_Vancomycin)

10.6 Visualización de los resultados

genenames <- AnnotationDbi::select(mouse4302.db, rownames(fit.main),
    c("SYMBOL"))$SYMBOL

# Función para generar volcanoplots
generar_volcanoplot <- function(index){
  
  volcanoplot(fit.main, coef = index, highlight = 6, names = genenames, 
              main = colnames(cont.matrix)[index])

  abline(v=c(-3,3))
}

# Representamos volcanoplots para cada contraste
generar_volcanoplot(1)
generar_volcanoplot(2)
generar_volcanoplot(3)

# Resumen genes diferencilamente expresados en cada contraste
summa.fit <- decideTests(fit.main, p.value = 0.05, lfc = 3)
summary(summa.fit)

# Diagrama de Venn
vc<- vennCounts(summa.fit)
vennDiagram(vc, include=c("up", "down"),
    counts.col=c("red", "blue"),
    circle.col = c("red", "blue", "green3"), cex=c(1,1,1))

10.7 Análisis de significación biológica

all_probes <- rownames(exprs(normData_filtered$eset))

entrezUniverse <- AnnotationDbi::select(mouse4302.db, all_probes,
    "ENTREZID")$ENTREZID

entrezUniverse <- entrezUniverse[!duplicated(entrezUniverse)]

# Infectados vs no infectados sin tratamiento
entrezTop_untreated <- topTabAnot_Untreated$ENTREZID 
entrezTop_untreated <- entrezTop_untreated[!duplicated(entrezTop_untreated)]

GOparams_unt <- new("GOHyperGParams", geneIds = entrezTop_untreated, universeGeneIds = entrezUniverse,
    annotation = "mouse4302.db", ontology = "BP", pvalueCutoff = 0.01)

GOhyper_unt <- hyperGTest(GOparams_unt)

head(summary(GOhyper_unt))

ego_unt <- enrichGO(gene = entrezTop_untreated,
                universe = as.character(entrezUniverse),
                keyType = "ENTREZID",
                OrgDb = org.Mm.eg.db,
                ont = "BP",
                pAdjustMethod = "BH",
                qvalueCutoff = 0.05,
                readable = TRUE)

dotplot(ego_unt, showCategory=9)

# Infectados vs no infectados tratados con linezolid
entrezTop_linezolid <- topTabAnot_Linezolid$ENTREZID 
entrezTop_linezolid <- entrezTop_linezolid[!duplicated(entrezTop_linezolid)]

GOparams_lin <- new("GOHyperGParams", geneIds = entrezTop_linezolid, universeGeneIds = entrezUniverse,
    annotation = "mouse4302.db", ontology = "BP", pvalueCutoff = 0.01)

GOhyper_lin <- hyperGTest(GOparams_lin)

head(summary(GOhyper_lin))

ego_lin <- enrichGO(gene = entrezTop_linezolid,
                universe = as.character(entrezUniverse),
                keyType = "ENTREZID",
                OrgDb = org.Mm.eg.db,
                ont = "BP",
                pAdjustMethod = "BH",
                qvalueCutoff = 0.05,
                readable = TRUE)

dotplot(ego_lin, showCategory=9)

# Infectados vs no infectados tratados con vancomicina
entrezTop_vancomicina <- topTabAnot_Vancomycin$ENTREZID 

entrezTop_vancomicina <- entrezTop_vancomicina[!duplicated(entrezTop_vancomicina)]

GOparams_van <- new("GOHyperGParams", geneIds = entrezTop_vancomicina, universeGeneIds = entrezUniverse,
    annotation = "mouse4302.db", ontology = "BP", pvalueCutoff = 0.01)

GOhyper_van <- hyperGTest(GOparams_van)

head(summary(GOhyper_van))

ego_van <- enrichGO(gene = entrezTop_vancomicina,
                universe = as.character(entrezUniverse),
                keyType = "ENTREZID",
                OrgDb = org.Mm.eg.db,
                ont = "BP",
                pAdjustMethod = "BH",
                qvalueCutoff = 0.05,
                readable = TRUE)
sessionInfo()
## R version 4.4.2 (2024-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26100)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=Spanish_Spain.utf8  LC_CTYPE=Spanish_Spain.utf8   
## [3] LC_MONETARY=Spanish_Spain.utf8 LC_NUMERIC=C                  
## [5] LC_TIME=Spanish_Spain.utf8    
## 
## time zone: Europe/Madrid
## tzcode source: internal
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] GO.db_3.20.0               mouse4302cdf_2.18.0       
##  [3] dplyr_1.1.4                clusterProfiler_4.14.6    
##  [5] GOstats_2.72.0             graph_1.84.1              
##  [7] Category_2.72.0            Matrix_1.7-1              
##  [9] biomaRt_2.62.1             tibble_3.2.1              
## [11] limma_3.62.2               arrayQualityMetrics_3.62.0
## [13] colorspace_2.1-1           ggplot2_3.5.1             
## [15] affy_1.84.0                mouse4302.db_3.13.0       
## [17] org.Mm.eg.db_3.20.0        AnnotationDbi_1.68.0      
## [19] genefilter_1.88.0          oligo_1.70.0              
## [21] Biostrings_2.74.1          GenomeInfoDb_1.42.3       
## [23] XVector_0.46.0             IRanges_2.40.1            
## [25] S4Vectors_0.44.0           Biobase_2.66.0            
## [27] oligoClasses_1.68.0        BiocGenerics_0.52.0       
## 
## loaded via a namespace (and not attached):
##   [1] splines_4.4.2               ggplotify_0.1.2            
##   [3] bitops_1.0-9                filelock_1.0.3             
##   [5] R.oo_1.27.0                 preprocessCore_1.68.0      
##   [7] affxparser_1.78.0           XML_3.99-0.18              
##   [9] rpart_4.1.23                lifecycle_1.0.4            
##  [11] httr2_1.1.1                 lattice_0.22-6             
##  [13] base64_2.0.2                backports_1.5.0            
##  [15] magrittr_2.0.3              Hmisc_5.2-3                
##  [17] sass_0.4.9                  rmarkdown_2.29             
##  [19] jquerylib_0.1.4             yaml_2.3.10                
##  [21] ggtangle_0.0.6              gcrma_2.78.0               
##  [23] askpass_1.2.1               cowplot_1.1.3              
##  [25] DBI_1.2.3                   RColorBrewer_1.1-3         
##  [27] abind_1.4-8                 zlibbioc_1.52.0            
##  [29] GenomicRanges_1.58.0        purrr_1.0.4                
##  [31] R.utils_2.13.0              RCurl_1.98-1.16            
##  [33] yulab.utils_0.2.0           nnet_7.3-19                
##  [35] rappdirs_0.3.3              GenomeInfoDbData_1.2.13    
##  [37] enrichplot_1.26.6           ggrepel_0.9.6              
##  [39] AnnotationForge_1.48.0      tidytree_0.4.6             
##  [41] setRNG_2024.2-1             annotate_1.84.0            
##  [43] svglite_2.1.3               gridSVG_1.7-5              
##  [45] codetools_0.2-20            DelayedArray_0.32.0        
##  [47] DOSE_4.0.0                  xml2_1.3.6                 
##  [49] tidyselect_1.2.1            aplot_0.2.5                
##  [51] farver_2.1.2                UCSC.utils_1.2.0           
##  [53] matrixStats_1.5.0           BiocFileCache_2.14.0       
##  [55] base64enc_0.1-3             illuminaio_0.48.0          
##  [57] jsonlite_1.9.0              Formula_1.2-5              
##  [59] survival_3.7-0              iterators_1.0.14           
##  [61] systemfonts_1.2.1           foreach_1.5.2              
##  [63] tools_4.4.2                 progress_1.2.3             
##  [65] treeio_1.30.0               Rcpp_1.0.14                
##  [67] glue_1.8.0                  gridExtra_2.3              
##  [69] SparseArray_1.6.2           xfun_0.51                  
##  [71] qvalue_2.38.0               MatrixGenerics_1.18.1      
##  [73] affyPLM_1.82.0              withr_3.0.2                
##  [75] BiocManager_1.30.25         fastmap_1.2.0              
##  [77] latticeExtra_0.6-30         openssl_2.3.2              
##  [79] digest_0.6.37               gridGraphics_0.5-1         
##  [81] BeadDataPackR_1.58.0        R6_2.6.1                   
##  [83] jpeg_0.1-10                 RSQLite_2.3.9              
##  [85] R.methodsS3_1.8.2           tidyr_1.3.1                
##  [87] generics_0.1.3              hexbin_1.28.5              
##  [89] data.table_1.17.0           beadarray_2.56.0           
##  [91] prettyunits_1.2.0           httr_1.4.7                 
##  [93] htmlwidgets_1.6.4           S4Arrays_1.6.0             
##  [95] pkgconfig_2.0.3             gtable_0.3.6               
##  [97] blob_1.2.4                  hwriter_1.3.2.1            
##  [99] htmltools_0.5.8.1           fgsea_1.32.2               
## [101] RBGL_1.82.0                 GSEABase_1.68.0            
## [103] scales_1.3.0                png_0.1-8                  
## [105] ggfun_0.1.8                 knitr_1.49                 
## [107] rstudioapi_0.17.1           reshape2_1.4.4             
## [109] nlme_3.1-166                checkmate_2.3.2            
## [111] curl_6.2.1                  cachem_1.1.0               
## [113] stringr_1.5.1               parallel_4.4.2             
## [115] foreign_0.8-87              vsn_3.74.0                 
## [117] pillar_1.10.1               grid_4.4.2                 
## [119] vctrs_0.6.5                 ff_4.5.2                   
## [121] dbplyr_2.5.0                xtable_1.8-4               
## [123] cluster_2.1.6               htmlTable_2.4.3            
## [125] Rgraphviz_2.50.0            evaluate_1.0.3             
## [127] cli_3.6.4                   compiler_4.4.2             
## [129] rlang_1.1.5                 crayon_1.5.3               
## [131] labeling_0.4.3              interp_1.1-6               
## [133] fs_1.6.5                    plyr_1.8.9                 
## [135] stringi_1.8.4               deldir_2.0-4               
## [137] BiocParallel_1.40.0         munsell_0.5.1              
## [139] lazyeval_0.2.2              GOSemSim_2.32.0            
## [141] patchwork_1.3.0             hms_1.1.3                  
## [143] bit64_4.6.0-1               KEGGREST_1.46.0            
## [145] statmod_1.5.0               SummarizedExperiment_1.36.0
## [147] igraph_2.1.4                memoise_2.0.1              
## [149] affyio_1.76.0               bslib_0.9.0                
## [151] ggtree_3.14.0               fastmatch_1.1-6            
## [153] bit_4.6.0                   gson_0.1.0                 
## [155] ape_5.8-1