# Modelos lineales con limma

Francisco García García (2015-01-27)

### 1. Objetivo del estudio y diseño del experimento

### 2. Lectura y analisis exploratorio de los datos

### 3. Organización de los datos en R

### 4. Deteccion de diferencias significativas: modelos lineales con el paquete limma de Bioconductor

### 5. Bibliografia y enlaces de interes

## 1. Objetivo del estudio y diseño del experimento

El objetivo del estudio es determinar el número de ciclos necesario para detectar diferencias en la intensidad de un grupo de metabolitos

Diseño del experimento. Se realizaron un conjunto de experimentos con la misma estructura: 5 puntos de tiempo en términos de ciclos: 0 (control), 1, 2, 3, 4, 5 donde se cuantificaron la intensidad de 48 metabolitos en 2 voluntarios. Para cada punto de tiempo / metabolito / voluntario se dispone de varias réplicas (3 a 5)

## 2. Lectura y análisis exploratorio de los datos en R

Indicamos el directorio donde esta el input y leemos los datos:

    datos <- read.csv("datos/cycles.csv",sep = "\t", header = T)

    #> Warning in file(file, "rt"): cannot open file 'datos/cycles.csv': No such
    #> file or directory

    #> Error in file(file, "rt"): cannot open the connection

Estructura de los datos:

    str(datos)

    #> Error in str(datos): object 'datos' not found

Dimensión de la matriz de datos (número de filas x número de columnas):

    dim(datos)

    #> Error in eval(expr, envir, enclos): object 'datos' not found

Exploramos la muestra:

    attach(datos)

    #> Error in attach(datos): object 'datos' not found

    summary(datos)

    #> Error in summary(datos): object 'datos' not found

    boxplot(datos[,c(6:53)], col = "green", main = "Box-plots metabolitos")

    #> Error in boxplot(datos[, c(6:53)], col = "green", main = "Box-plots metabolitos"): object 'datos' not found

    table(ciclos)

    #> Error in table(ciclos): object 'ciclos' not found

    table(replicate)

    #> Error in unique.default(x, nmax = nmax): unique() applies only to vectors

    table(VOLUNTARIO)

    #> Error in table(VOLUNTARIO): object 'VOLUNTARIO' not found

Representamos los niveles de los metabolitos por voluntario:

    par(mfrow=c(2,1))
    boxplot(datos[VOLUNTARIO == 1,c(6:53)], col = "green", main = "Box-plots metabolitos, VOLUNTARIO 1")

    #> Error in boxplot(datos[VOLUNTARIO == 1, c(6:53)], col = "green", main = "Box-plots metabolitos, VOLUNTARIO 1"): object 'datos' not found

    boxplot(datos[VOLUNTARIO == 2,c(6:53)], col = "blue", main = "Box-plots metabolitos, VOLUNTARIO 2")

    #> Error in boxplot(datos[VOLUNTARIO == 2, c(6:53)], col = "blue", main = "Box-plots metabolitos, VOLUNTARIO 2"): object 'datos' not found

    par(mfrow=c(1,1))

Representamos los niveles de los metabolitos segun el numero de ciclos:

    par(mfrow=c(3,2))
    boxplot(datos[ciclos == 0,c(6:53)], col = "green", main = "Box-plots metabolitos, 0 ciclos")

    #> Error in boxplot(datos[ciclos == 0, c(6:53)], col = "green", main = "Box-plots metabolitos, 0 ciclos"): object 'datos' not found

    boxplot(datos[ciclos == 1,c(6:53)], col = "blue", main = "Box-plots metabolitos, 1 ciclos")

    #> Error in boxplot(datos[ciclos == 1, c(6:53)], col = "blue", main = "Box-plots metabolitos, 1 ciclos"): object 'datos' not found

    boxplot(datos[ciclos == 2,c(6:53)], col = "red", main = "Box-plots metabolitos, 2 ciclos")

    #> Error in boxplot(datos[ciclos == 2, c(6:53)], col = "red", main = "Box-plots metabolitos, 2 ciclos"): object 'datos' not found

    boxplot(datos[ciclos == 3,c(6:53)], col = "pink", main = "Box-plots metabolitos, 3 ciclos")

    #> Error in boxplot(datos[ciclos == 3, c(6:53)], col = "pink", main = "Box-plots metabolitos, 3 ciclos"): object 'datos' not found

    boxplot(datos[ciclos == 4,c(6:53)], col = "yellow", main = "Box-plots metabolitos, 4 ciclos")

    #> Error in boxplot(datos[ciclos == 4, c(6:53)], col = "yellow", main = "Box-plots metabolitos, 4 ciclos"): object 'datos' not found

    boxplot(datos[ciclos == 5,c(6:53)], col = "grey", main = "Box-plots metabolitos, 5 ciclos")

    #> Error in boxplot(datos[ciclos == 5, c(6:53)], col = "grey", main = "Box-plots metabolitos, 5 ciclos"): object 'datos' not found

    par(mfrow=c(1,1))

Repetimos el grafico anterior pero con transformacion logaritmica para visualizar mejor los datos:

    par(mfrow=c(3,2))
    boxplot(log(datos[ciclos == 0,c(6:53)]), col = "green", main = "Box-plots metabolitos, 0 ciclos")

    #> Error in boxplot(log(datos[ciclos == 0, c(6:53)]), col = "green", main = "Box-plots metabolitos, 0 ciclos"): object 'datos' not found

    boxplot(log(datos[ciclos == 1,c(6:53)]), col = "blue", main = "Box-plots metabolitos, 1 ciclos")

    #> Error in boxplot(log(datos[ciclos == 1, c(6:53)]), col = "blue", main = "Box-plots metabolitos, 1 ciclos"): object 'datos' not found

    boxplot(log(datos[ciclos == 2,c(6:53)]), col = "red", main = "Box-plots metabolitos, 2 ciclos")

    #> Error in boxplot(log(datos[ciclos == 2, c(6:53)]), col = "red", main = "Box-plots metabolitos, 2 ciclos"): object 'datos' not found

    boxplot(log(datos[ciclos == 3,c(6:53)]), col = "pink", main = "Box-plots metabolitos, 3 ciclos")

    #> Error in boxplot(log(datos[ciclos == 3, c(6:53)]), col = "pink", main = "Box-plots metabolitos, 3 ciclos"): object 'datos' not found

    boxplot(log(datos[ciclos == 4,c(6:53)]), col = "yellow", main = "Box-plots metabolitos, 4 ciclos")

    #> Error in boxplot(log(datos[ciclos == 4, c(6:53)]), col = "yellow", main = "Box-plots metabolitos, 4 ciclos"): object 'datos' not found

    boxplot(log(datos[ciclos == 5,c(6:53)]), col = "grey", main = "Box-plots metabolitos, 5 ciclos")

    #> Error in boxplot(log(datos[ciclos == 5, c(6:53)]), col = "grey", main = "Box-plots metabolitos, 5 ciclos"): object 'datos' not found

    par(mfrow=c(1,1))

## 3. Organización de los datos en R

Generamos un identificador para cada muestra: ciclo_nreplica_voluntario. Por ejemplo: 0_3_1 indicaria ciclo 0, replica 3, voluntario 1

    datos[,"ID"] <- paste(ciclos, replicate, VOLUNTARIO, sep ="_")

    #> Error in paste(ciclos, replicate, VOLUNTARIO, sep = "_"): object 'ciclos' not found

    head(datos[,c("ciclos", "replicate",  "VOLUNTARIO", "ID")])

    #> Error in head(datos[, c("ciclos", "replicate", "VOLUNTARIO", "ID")]): object 'datos' not found

    ndatos <- datos[, c(6:53)]

    #> Error in eval(expr, envir, enclos): object 'datos' not found

    rownames(ndatos) <- datos[, "ID"]

    #> Error in eval(expr, envir, enclos): object 'datos' not found

Transpongo la matriz de datos, de modo que cada columna correspondera al grupo de todos los metabolitos para un determinado ciclo, replica y voluntario

    tndatos <- t(ndatos)

    #> Error in t(ndatos): object 'ndatos' not found

Esta es la estructura que deben tener los datos para aplicar modelos lineales con el paquete limma:

    head(tndatos[, c(1:5)])

    #> Error in head(tndatos[, c(1:5)]): object 'tndatos' not found

    dim(tndatos)

    #> Error in eval(expr, envir, enclos): object 'tndatos' not found

Volvemos a explorar los datos mediante clustering y PCAplot. Previamente cargamos las siguientes funciones:

    source ("function_arbol_4.r")

    #> Warning in file(filename, "r", encoding = encoding): cannot open file
    #> 'function_arbol_4.r': No such file or directory

    #> Error in file(filename, "r", encoding = encoding): cannot open the connection

    source ("function_pcaGenes_2.r")

    #> Warning in file(filename, "r", encoding = encoding): cannot open file
    #> 'function_pcaGenes_2.r': No such file or directory

    #> Error in file(filename, "r", encoding = encoding): cannot open the connection

## 4. Detección de diferencias significativas: modelos lineales con el paquete limma de Bioconductor

Instalamos el paquete limma en RStudio incluyendo el siguiente codigo en R:

source(“http://bioconductor.org/biocLite.R”)

biocLite(“limma”)

Cargamos la libreria limma

    library(limma)

Cargamos una funcion que ordenara los resultados que obtengamos en limma, preparando varios ficheros de texto de resultados (uno por cada comparacion):

    source ("function_save.limma.results_10.r")

    #> Warning in file(filename, "r", encoding = encoding): cannot open file
    #> 'function_save.limma.results_10.r': No such file or directory

    #> Error in file(filename, "r", encoding = encoding): cannot open the connection

Preparamos las variables factor (son las que proporcionan diferentes grupos en los que queremos detectar posibles diferencias)

    ciclos      <- as.factor(substr(colnames(tndatos), 1,1))

    #> Error in is.data.frame(x): object 'tndatos' not found

    voluntarios <- as.factor(substr(colnames(tndatos), 5,5))

    #> Error in is.data.frame(x): object 'tndatos' not found

    table(ciclos)

    #> Error in table(ciclos): object 'ciclos' not found

    table(voluntarios)

    #> Error in table(voluntarios): object 'voluntarios' not found

Matriz de diseño para el modelo lineal

    des.mat <- model.matrix (~ ciclos +  voluntarios)

    #> Error in eval(expr, envir, enclos): object 'ciclos' not found

    rownames(des.mat) <- colnames(tndatos)

    #> Error in is.data.frame(x): object 'tndatos' not found

    des.mat

    #> Error in eval(expr, envir, enclos): object 'des.mat' not found

¿Que preguntas queremos contestar con el modelo? Utilizaremos los contrastes para representar las comparaciones que queremos valorar:

    cont.mat <- cbind(           ciclos5.vs.ciclos4   =  c(0, 0, 0, 0,-1,1,0),
                                 ciclos5.vs.ciclos3   =  c(0, 0, 0,-1, 0,1,0),
                                 ciclos5.vs.ciclos2   =  c(0, 0,-1, 0, 0,1,0),
                                 ciclos5.vs.ciclos1   =  c(0,-1, 0, 0, 0,1,0),
                                 ciclos4.vs.ciclos3   =  c(0, 0, 0,-1, 1,0,0),
                                 ciclos4.vs.ciclos2   =  c(0, 0,-1, 0, 1,0,0),
                                 ciclos4.vs.ciclos1   =  c(0,-1, 0, 0, 1,0,0),
                                 ciclos3.vs.ciclos2   =  c(0, 0,-1, 1, 0,0,0),
                                 ciclos3.vs.ciclos1   =  c(0,-1, 0, 1, 0,0,0),
                                 ciclos2.vs.ciclos1   =  c(0,-1, 1, 0, 0,0,0),
                                 ciclos5.vs.ciclos0   =  c(0,0,0,0,0,1,0),
                                 ciclos4.vs.ciclos0   =  c(0,0,0,0,1,0,0),
                                 ciclos3.vs.ciclos0   =  c(0,0,0,1,0,0,0),
                                 ciclos2.vs.ciclos0   =  c(0,0,1,0,0,0,0),
                                 ciclos1.vs.ciclos0   =  c(0,1,0,0,0,0,0),
                                 voluntario2.vs.voluntario1 = c(0,0,0,0,0,0,1))
    cont.mat

    #>      ciclos5.vs.ciclos4 ciclos5.vs.ciclos3 ciclos5.vs.ciclos2
    #> [1,]                  0                  0                  0
    #> [2,]                  0                  0                  0
    #> [3,]                  0                  0                 -1
    #> [4,]                  0                 -1                  0
    #> [5,]                 -1                  0                  0
    #> [6,]                  1                  1                  1
    #> [7,]                  0                  0                  0
    #>      ciclos5.vs.ciclos1 ciclos4.vs.ciclos3 ciclos4.vs.ciclos2
    #> [1,]                  0                  0                  0
    #> [2,]                 -1                  0                  0
    #> [3,]                  0                  0                 -1
    #> [4,]                  0                 -1                  0
    #> [5,]                  0                  1                  1
    #> [6,]                  1                  0                  0
    #> [7,]                  0                  0                  0
    #>      ciclos4.vs.ciclos1 ciclos3.vs.ciclos2 ciclos3.vs.ciclos1
    #> [1,]                  0                  0                  0
    #> [2,]                 -1                  0                 -1
    #> [3,]                  0                 -1                  0
    #> [4,]                  0                  1                  1
    #> [5,]                  1                  0                  0
    #> [6,]                  0                  0                  0
    #> [7,]                  0                  0                  0
    #>      ciclos2.vs.ciclos1 ciclos5.vs.ciclos0 ciclos4.vs.ciclos0
    #> [1,]                  0                  0                  0
    #> [2,]                 -1                  0                  0
    #> [3,]                  1                  0                  0
    #> [4,]                  0                  0                  0
    #> [5,]                  0                  0                  1
    #> [6,]                  0                  1                  0
    #> [7,]                  0                  0                  0
    #>      ciclos3.vs.ciclos0 ciclos2.vs.ciclos0 ciclos1.vs.ciclos0
    #> [1,]                  0                  0                  0
    #> [2,]                  0                  0                  1
    #> [3,]                  0                  1                  0
    #> [4,]                  1                  0                  0
    #> [5,]                  0                  0                  0
    #> [6,]                  0                  0                  0
    #> [7,]                  0                  0                  0
    #>      voluntario2.vs.voluntario1
    #> [1,]                          0
    #> [2,]                          0
    #> [3,]                          0
    #> [4,]                          0
    #> [5,]                          0
    #> [6,]                          0
    #> [7,]                          1

Ajustamos el modelo:

    fit <- lmFit(tndatos, design =des.mat)

    #> Error in is(object, "list"): object 'tndatos' not found

    fit <- contrasts.fit (fit, contrasts = cont.mat)

    #> Error in NCOL(fit$coefficients): object 'fit' not found

    fit

    #> Error in eval(expr, envir, enclos): object 'fit' not found

Determinamos los p-valores correspondientes a cada contraste y metabolito

    res <- eBayes (fit)  #getting p-values

    #> Error in ebayes(fit = fit, proportion = proportion, stdev.coef.lim = stdev.coef.lim, : object 'fit' not found

    res

    #> Error in eval(expr, envir, enclos): object 'res' not found

    res$p.adj <- apply (res$p.value, 2, p.adjust, method = "BH")  #determino el p-valor ajustado

    #> Error in apply(res$p.value, 2, p.adjust, method = "BH"): object 'res' not found

    results <- decideTests(res, p.value= 0.05)  #considera  los p-valores ajustados

    #> Error in is(object, "MArrayLM"): object 'res' not found

Resumen de resultados significativos para cada una de las comparaciones indicadas:

- Un valor 0 significa que NO hay diferencias entre los dos grupos comparados.
- Un valor 1 o -1 indica que existen diferencias significativas entre ambos grupos, siendo mayor el nivel del metabolito del grupo 1 que el segundo cuando es 1 y el nivel de metabolito es mayor en el grupo segundo que en el primero cuando el valor es -1
- Un metabolito es significativo (muestra diferencias estadísticamente significativas entre los dos grupos comparados) si presenta un p-valor ajustado < 0.05

    summary(results)

    #> Error in summary(results): object 'results' not found

Con los diagramas de Venn representamos los metabolitos significativos en todas las comparaciones que tienen como segundo grupo la referencia basal: 0 ciclos.

    vennDiagram(results[,c(11:15)], include = "down")

    #> Error in is(object, "VennCounts"): object 'results' not found

    vennDiagram(results[,c(11:15)], include = "up")

    #> Error in is(object, "VennCounts"): object 'results' not found

Más informacion sobre estos diagramas: ?vennDiagram Guardamos los datos correspondientes a los metabolitos:

    cat ("ID\t", file = "tdatos_metabolitos.txt")
    write.table(tndatos, "tdatos_metabolitos.txt", append = TRUE, sep = "\t", quote = F, col.names = T, row.names= T)

    #> Error in is.data.frame(x): object 'tndatos' not found

Guardamos resultados del análisis estadístico para cada una de las comparaciones descritas:

    save.limma.results (res, geneID = "ID", ext = ".txt", decimal = ".", datatofatiscandir = "nosave", includeFoldChange = TRUE, resdir= getwd())

    #> Error in eval(expr, envir, enclos): could not find function "save.limma.results"

## 5. Bibliografía y enlaces de interés:

- Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W and Smyth GK (2015). “limma powers differential expression analyses for RNA-sequencing and microarray studies.” Nucleic Acids Research, 43

- Quick-R Web con recursos para trabajar con R.

- r-tutor An R Introduction to Statistics.