
---
title: "Informe de Resultados de qPCR"
author: "Elena Garc�a Bravo"
date: '`r format(Sys.time(), "%d de %B de %Y")`'

#output:
#  html_document:
#     toc: true
#     toc_float: true
#     number_sections: yes
#     highlight: default
    
output:
 pdf_document:
   toc: true
   number_sections: yes

header-includes:
- \usepackage[utf8]{inputenc}
- \usepackage[spanish, es-tabla]{babel}
- \usepackage{amsmath,amssymb}

bibliography: bibliografia.bib
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, comment = NULL, error=FALSE, message = FALSE, warning = FALSE)
options(width = 90, digits = 4)
library(knitr)
```


# Pre�mbulo

Este es un documento din�mico que contiene un flujo de trabajo de an�lisis de datos de qPCR. Ha sido dise�ado para leer datos de Cq de placas predise�adas de qPCR miRCURY LNAT (Exiqon), medidas en un termociclador ABI PRISM 7000 Sequence Detection System. 

El formato de los datos deber� cumplir las siguientes caracter�sticas:

* Cada placa de qPCR tendr� la misma distribuci�n: cada gen medido se encuentra en el mismo pocillo en todas las placas. 

* A cada placa le corresponde un �nico fichero de texto con 1 columna y tantas filas como pocillos, con los valores Cq.

* Los ficheros han de ser de texto plano (csv, txt).

* Se necesita un fichero adicional `Genesinfo.txt` con informaci�n sobre el nombre del gen que se eval�a en cada pocillo.

* Se necesita otro fichero adicional `samplesinfo.txt` que relaciona el nombre del fichero con las condiciones experimentales correspondientes.

El documento viene acompa�ado de un set de datos de muestra, en el que se midi� la expresi�n g�nica de miRNAs purificados de exosomas de orina en humanos. Algunos comentarios del informe no son din�micos, sino espec�ficos de este set de datos. En dichos casos se indicar�.


Tras la compilaci�n del documento, se generar� un informe con el an�lisis y los resultados de los datos de qPCR introducidos por el usuario. Adem�s, los gr�ficos y tablas generados se guardar�n autom�ticamente en archivos independientes para que el usuario pueda disponer de ellos.

En este documento se compagina la ejecuci�n de c�digo necesario para el an�lisis (en lenguaje R) con explicaciones de cada etapa. El c�digo instalar� el/los paquetes de R y bioconductor necesarios para ejecutar el pipeline, en el caso de que no estuviesen instalados.



```{r instalar_paquetes, echo=FALSE, include=FALSE}
# listas de paquetes
paquetes_bioc <- c("HTqPCR")
paquetes_nuevos <- c("pander")


# cada paquete ha de descargarse de un repositorio
# respositorio bioconductor
nuevos_bioc <- paquetes_bioc[!(paquetes_bioc %in% installed.packages()[,"Package"])]
if(length(nuevos_bioc)){
  source("https://bioconductor.org/biocLite.R")
  biocLite(nuevos_bioc)
}

# respositorio CRAN
nuevos <- paquetes_nuevos[!(paquetes_nuevos %in% installed.packages()[,"Package"])]
if(length(nuevos)){
  install.packages(nuevos)
}
```


```{r cargar_paquetes, echo=FALSE, results='hide'}
# Una vez instalados los paquetes necesarios, vamos a cargarlos
library(pander)
library(HTqPCR)
```

# Lectura de los datos

## Par�metros a modificar por el usuario

En primer lugar, para la lectura de los datos el usuario deber� editar algunas peque�as partes del c�digo en el documento `pipeline_qPCR.Rmd` para adaptarlo a sus datos. El usuario deber� establecer:

* Directorio de trabajo: el directorio donde se encuentran sus datos. Ser� necesario modificar, en el c�digo, el texto `"./dades"` y cambiarlo por la ruta de carpetas de los datos a analizar.

```{r directorio}
# MODIFICABLE POR EL USUARIO
# ruta de carpetas donde se encuentran los datos a analizar:
ruta <- "./dades"
```
* N�mero de pocillos en cada placa. Podr�n ser 96 � 384.

```{r pocillos}
# MODIFICABLE POR EL USUARIO
# Escoge el n�mero de pocillos de la placa (96/384):
pocillos <- 96
```


```{r setwd, echo=FALSE, include=FALSE}
# Para que R trabaje en esta carpeta:
setwd(ruta) 
```


Una vez puesto a punto el c�digo, los datos se estar�n leyendo desde el directorio `r ruta`.


## Descripci�n

Para la lectura de los datos se utiliza el paquete `HTqPCR` (@dvinge2009htqpcr).

```{r leerCq, echo=FALSE, include=FALSE}
# leemos el fichero con la informacion de las muestras:
targets <- read.table(file.path(ruta, "samplesinfo.txt"),sep = ",", header=TRUE)

# leemos los ficheros con los valores de Cq:
raw <- readCtData(files = targets$Sample, path = ruta,
                  format = "plain", n.features = pocillos, 
                  column.info = list(position=1, Ct=3), 
                  n.data = 1, 
                  sep = ",", header = TRUE, skip = 2)


```

Durante la lectura de los datos, puede suceder que falten datos de Cq para algunas medidas.  Cuando se da el caso de que faltan algunos valores, para poder continuar con el an�lisis se recomienda (@dvinge2009htqpcr) sustituir estos valores ausentes por un valor num�rico elevado, como por ejemplo 40. Si el umbral de detecci�n est� en el ciclo n�mero 40 se puede considerar que no hay amplificaci�n de ARN que detectar, y la expresi�n ser� por tanto cercana a 0.

En este caso, en el set de datos a analizar faltan datos de Cq de: `r length(exprs(raw)[is.na(exprs(raw))])` pocillos.

```{r}
# Para que no aparezca el aviso de los valores NA (ausentes), 
# se suele sustituir por un valor num�rico elevado como 40
# MODIFICABLE POR EL USUARIO
limitcq <- 40
exprs(raw)[is.na(exprs(raw))] <- limitcq
```


```{r genes, echo=FALSE, include=FALSE}
# leemos y almacenamos los nombres de los genes:
genes <- file.path(ruta, "Genesinfo.csv")
genesinfo <- read.csv(genes, header = TRUE)
genelist <- genesinfo[,1]

# asignamos los nombres de los genes a nuestros datos de qPCR:
featureNames(raw) <- as.vector(genelist)
```

El set de datos actual consiste en un conjunto de `r nrow(targets)` placas de qPCR o muestras. Estas muestras se dividen en `r length(unique(targets$Treatment))` tratamientos experimentales denominados: `r unique(targets$Treatment)`. En cada tratamiento se ha medido `r length(unique(targets$Response))` respuestas, denominadas: `r unique(targets$Response)`. Tenemos, en total, `r length(unique(targets$Treatment))*length(unique(targets$Response))` grupos experimentales distintos.

A continuaci�n se presenta una lista de los miRNA y controles medidos en cada placa. Cada n�mero corresponde a la posici�n del pocillo correspondiente.

```{r genelist, echo=FALSE, include=FALSE}
# Creamos una tabla con la lista de genes:
genes <- paste(c(1:length(genelist)), genelist)
pander(t(matrix(genes, nrow = 8, ncol = 12)), justify="left")
```

## Resumen num�rico de los datos

En la siguiente tabla se muestra un resumen estad�stico de los datos de expresi�n global (Cq o ciclo de detecci�n) de cada muestra. Esta tabla se guardar� autom�ticamente en un archivo de texto denominado `Results_raw_summary.csv`.

```{r resumen, echo=FALSE}
# Creamos una tabla-resumen:
pander(t(as.matrix(summary(raw))))
```

```{r resumencsv, echo=FALSE, include=FALSE}
# salvamos el resumen en un archivo de texto:
write.csv(t(as.matrix(summary(raw))), file="Results_raw_summary.csv") 
```


```{r grupos, echo=FALSE, include=FALSE}
# ESTE TROZO DE C�DIGO ES ESPEC�FICO DE LOS DATOS DEL PROYECTO
# En el dataset 'targets', a�adimos una nueva variable 'Group' 
# que indique grupo y respuesta:
targets$Group <- NA

# recodificamos la variable Group:
targets$Group[targets$Treatment == 1 & targets$Response == "A"] <- "1A"
targets$Group[targets$Treatment == 1 & targets$Response == "B"] <- "1B"
targets$Group[targets$Treatment == 2 & targets$Response == "A"] <- "2A"
targets$Group[targets$Treatment == 2 & targets$Response == "B"] <- "2B"
```



# Control de calidad 
## Control de las muestras

```{r colores, echo=FALSE, include=FALSE}
# preparaci�n de elementos para plots de visualizacion preliminar

# paletas de color
colores2 <- c('lightcoral', 'lightblue4') 
colores4 <- c(
  rep(c('lightcoral', 'lightblue4'), 7), 
  rep(c('peachpuff', 'lightblue'), 7)
  )

# leyenda
leyenda_grupos <- c("bottomright", cex = 0.8,
                   legend = c("1A", "2A", "1B", "2B"),
                   fill = c('lightcoral', 'peachpuff', 'lightblue4', 'lightblue'))
```

En el siguiente gr�fico se muestra una vista preliminar de los valores de expresi�n Cq en crudo de cada uno de los genes/controles medidos. Este gr�fico se guardar� autom�ticamente en un archivo denominado `Results_raw_plotCtOverview.png`.

```{r plotCt_grupos, echo=FALSE, message=FALSE, results='hide'}

### PLOT EXPRESION Cq RAW

# muestras agrupadas por grupo y respuesta

par(mfrow=c(2,2))


plotCtOverview(raw[1:24], 
               genes = featureNames(raw), groups = (targets$Group),
               ylim = c(0, 45),
               col=colores4,
               legend = FALSE
)

plotCtOverview(raw[25:48], 
               genes = featureNames(raw), groups = (targets$Group),
               ylim = c(0, 45),
               col=colores4,
               legend = FALSE
)

plotCtOverview(raw[49:72], 
               genes = featureNames(raw), groups = (targets$Group),
               ylim = c(0, 45),
               col=colores4,
               legend = FALSE
)

plotCtOverview(raw[72:96], 
               genes = featureNames(raw), groups = (targets$Group),
               ylim = c(0, 45),
               col=colores4,
               legend = FALSE
)

legend("bottomright", cex = 0.5,
                   legend = c("1A", "2A", "1B", "2B"),
                   fill = c('lightcoral', 'peachpuff', 'lightblue4', 'lightblue'))

dev.copy(png,'Results_raw_plotCtOverview.png')
dev.off()
```

En el siguiente gr�fico se muestra una vista preliminar de los valores y la distribuci�n de expresi�n de cada muestra. Este gr�fico se guardar� autom�ticamente en el archivo `Results_raw_plotCtBoxes.png`.


```{r raw_boxplot, echo=FALSE, message=FALSE, results='hide'}
### BOXPLOT de cada paciente
plotCtBoxes(raw, col = colores4, las=2, names = paste("muestra", c(1:28)))

legend("bottomright", cex = 0.8,
        legend = c("1A", "2A", "1B", "2B"),
        fill = c('lightcoral', 'peachpuff', 'lightblue4', 'lightblue')
)

dev.copy(png,'Results_raw_plotCtBoxes.png')
dev.off()
```



En el siguiente gr�fico de jerarqu�a se observa c�mo se agrupan las muestras. Es com�n que las muestras se agrupen por grupos experimentales. Si esto sucede, y hay alguna de ellas muy alejada o agrupada con un grupo al que no corresponder�a, ser� apropiado descartarla para continuar el an�lisis.

```{r, echo=FALSE, results='hide'}
clusterCt(raw, type = "samples")
```

El siguiente gr�fico muestra un an�lisis de componentes principales. Si las muestras se agrupan por grupos experimentales ser� indicativo de que pueden existir diferencias significativas entre los tratamientos del experimento. Este gr�fico se generar� con el nombre de `Results_raw_PCA.png`.

```{r PCA_raw, echo=FALSE, results='hide'}
plotCtPCA(raw, features=FALSE, col = colores4)
legend("bottomright", cex = 0.8,
       legend = c("1A", "2A", "1B", "2B"),
       fill = c('lightcoral', 'peachpuff', 'lightblue4', 'lightblue'))
dev.copy(png,'Results_raw_PCA.png')
dev.off()
```

En este punto podemos decidir si deber�amos descartar alguna muestra que no se agrupe con el resto de su mismo grupo experimental. En nuestro set de datos de ejemplo, podr�a ser una buena idea eliminar la muestra 4 (perteneciente al grupo `r targets$Group[targets$Sample == "sample04.csv"]`).

## Controles adicionales

El protocolo para qPCR miRCURY LNAT (Exiqon) incluye una serie de ARNs control que se adicionan a la muestra en diferentes etapas para monitorizar los procesos de aislamiento del ARN, reversotranscripci�n a cDNA, y rendimiento de la reacci�n de amplificaci�n. Estos controles pueden servir para identificar y excluir muestras at�picas o 'outliers'.

Las placas predise�adas qPCR miRCURY LNAT (Exiqon) contienen secuencias para medir esos controles. La posici�n de estos controles en la placa var�a seg�n el modelo. Por eso, seg�n el modelo de placa utilizado, en esta etapa el usuario habr� de introducir la informaci�n sobre qu� pocillos de la placa contienen estos controles, siguiendo las plantillas de la [p�gina de Exiqon.](http://www.exiqon.com/plate-layout-files)


```{r controles Uni}
# MODIFICABLE POR EL USUARIO SI ES NECESARIO
# Asignar posiciones de los controles en la placa:
UniSp3 <- c(12, 36, 48)
UniSp6 <- 24
CelmiR39 <- 92
UniSp2 <- 93
UniSp4 <- 94
UniSp5 <- 95
Blank <- 96
```


```{r, echo=FALSE, include=FALSE}
# creamos listas con las posiciones y los nombres de los controles
controles_pocillos <- c(UniSp3, UniSp6, CelmiR39, UniSp2, UniSp4, UniSp5, Blank)
controles_nombre <- rep("Control", length(controles_pocillos))

# asignamos los nombres de los genes a nuestros datos de qPCR:
featureType(raw)[controles_pocillos] <- controles_nombre
```


### Control del proceso de extracci�n de ARN

Los controles UniSp2, UniSp4 and UniSp5 se adicionan a la muestra previamente a la extracci�n de ARN, y sirven para monitorizar dicho proceso. Cada uno de los controles tiene una concentraci�n final 100 veces menor que el anterior, siendo el m�s concentrado UniSp2 y UniSp5 el m�s diluido. Tras la reacci�n de PCR se espera que las diferencias en Cq entre dichos controles sea de 5-7 Cq. En la gr�fica siguiente se muestra el resultado de amplificaci�n de estos controles. El gr�fico se guardar� como `Results_plotCt_norm_unisp245.png`.


```{r extraccion ARN, echo=FALSE, results='hide'}

# CONTROL DEL PROCESO DE EXTRACCI�N de ARN
# UniSp2, UniSp4 and UniSp5
# stepwise diferences btwn thm is expected to be 5-7 Cq

plotCtOverview(raw[c(UniSp2, UniSp4, UniSp5)], 
               genes = featureNames(raw),
               ylim = c(0, 45),
               las=1,
               col=colores4,
               legend = FALSE
)

legend("bottomright", cex = 0.8,
       legend = c("1A", "2A", "1B", "2B"),
       fill = c('lightcoral', 'peachpuff', 'lightblue4', 'lightblue'))

dev.copy(png,'Results_plotCt_unisp245.png')
dev.off()
```


### Control de s�ntesis de cDNA

Un segundo control consiste en la adici�n de UniSp6 y cel-39-3p tras el aislamiento de ARN y previamente a la s�ntesis de cDNA. Es una manera de diagnosticar problemas en la reacci�n de retrotranscripci�n. El gr�fico generado a continuaci�n permite visualizar estos controles, y ser� almacenado como `Results_plotCt_unisp6CelmiR39.png`.
 

```{r cDNA, echo=FALSE, results='hide'}

plotCtOverview(raw[c(UniSp6, CelmiR39)], 
               genes = featureNames(raw),
               ylim = c(0, 45),
               las=1,
               col=colores4,
               legend = FALSE
)

legend("bottomright", cex = 0.8,
       legend = c("1A", "2A", "1B", "2B"),
       fill = c('lightcoral', 'peachpuff', 'lightblue4', 'lightblue'))

dev.copy(png,'Results_plotCt_unisp6CelmiR39.png')
dev.off()
```

### Control de eficacia de la qPCR

Este control tiene varias funciones:

* Permite distinguir resultados negativos verdaderos de falsos resultados negativos causados por un mal funcionamiento de la reacci�n de PCR.

* Permite calibrar y comparar entre placas. Si todo ha ido bien, los valores Cq deber�an ser muy similares dentro de la misma placa y especialmente entre placas.

A continuaci�n se muestran los niveles de amplificaci�n de estos controles, que se guardar�n en el archivo `Results_plotCt_Uni3.png`.

```{r interplate, echo=FALSE, results='hide'}
# interplate calibrator!!
# UniSp3 IPC
plotCtOverview(raw[UniSp3], 
               genes = featureNames(raw),
               ylim = c(0, 20),
               las=1,
               legend = FALSE,
               col = colores4
)

legend("bottomright", cex = 0.8,
       legend = c("1A", "2A", "1B", "2B"),
       fill = c('lightcoral', 'peachpuff', 'lightblue4', 'lightblue'))

dev.copy(png,'Results_plotCt_Uni3.png')
dev.off()
```




# Pretratamiento y filtrado

En esta etapa se asignan las categor�as 'Ok', 'Unreliable' o 'Undetermined' a cada medici�n de Cq con la funci�n `setCategory`. El usuario establece un rango de detecci�n de Cq, con unos limites m�nimo y m�ximo. 

El l�mite m�ximo est� establecido en 35 por defecto. Las medidas que superen el n�mero m�ximo de ciclos ser�n catalogadas como 'Undetermined', y ser�n aquellas para las que no se ha detectado expresi�n. Aqu� se incluyen los valores NA que previamente se hab�an establecido en `r limitcq` Cq. 

El l�mite m�nimo est� establecido en 10 por defecto. Las medidas que queden por debajo del n�mero m�ximo de ciclos ser�n catalogadas como 'Unreliable'.

Ambos l�mites pueden ser modificados por el usuario.

```{r}
# MODIFICABLE POR EL USUARIO
Cqmax <- 35
Cqmin <- 10
```

En el siguiente gr�fico podemos ver la proporci�n de muestras de cada una de las categor�as. En este estado tambi�n se puede tomar la decisi�n de descartar una muestra poco fiable. Por ejemplo, en el set de datos de muestra, ser�a buena idea descartar la muestra 4 como hab�amos se�alado anteriormente. El gr�fico se guardar� en `Results_categories.png`.

```{r, echo=FALSE, results='hide'}
raw.cat <- setCategory(raw, groups = targets$Group, Ct.max = 38, Ct.min = 1)
plotCtCategory(raw.cat)
dev.copy(png,'Results_categories.png')
dev.off()
```


# Normalizaci�n
## Elecci�n de los genes normalizadores

El gen normalizador es un gen cuya expresi�n sirve de base para calcular la expresi�n relativa del resto de genes a analizar. Normalmente, habr� sido sido escogido por el usuario durante el dise�o del experimento. En tal caso, la elecci�n debe estar debidamente justificada y dicho gen deber� estar validado experimentalmente como un buen gen de referencia @bustin2009miqe.

Si no se dispone de un gen de referencia, este pipeline calcula los genes que menos variaci�n presentan entre muestras, es decir, cuya desviaci�n est�ndar relativa es menor.

El usuario deber� escoger cu�ntos genes normalizadores quiere encontrar. Por defecto se han escogido 3.

```{r}
# MODIFICABLE POR EL USUARIO
numero_genes_ref <- 3
```


```{r, echo=FALSE}
# c�lculo de expresi�n media 
exp_desvest <- apply(X = exprs(raw), FUN = sd, MARGIN = 1)
exp_mean <- apply(X = exprs(raw), FUN = mean, MARGIN = 1)

# c�lculo de la desviaci�n est�ndar relativa a la expresi�n media
rel_desvest <- sort(exp_desvest/exp_mean)

candidatos <- setdiff(names(rel_desvest), as.character(genelist)[controles_pocillos])

genes_ref <- candidatos[1:numero_genes_ref]

```


```{r, echo=FALSE}
pander(as.matrix(rel_desvest[1:10]), justify="left")
```


En la lista precedente se muestran los 10 primeros genes que menos var�an en su expresi�n entre diferentes muestras o placas. Lo normal es que aparezcan algunos de los controles que se utilizan para evaluar la calidad de la qPCR. �stos deben descartarse como genes normalizadores, pues provienen de ARN que no formaba parte de la muestra original.

Se podr�n usar como genes de referencia para normalizar los `r numero_genes_ref` miRNA cuya expresi�n es m�s constante a lo largo de las muestras: `r genes_ref`.




## M�todo de normalizaci�n 

En este paso del an�lisis el usuario debe escoger un m�todo de normalizaci�n de la expresi�n g�nica, para hacer comparables los valores de expresi�n entre genes y muestras. La normalizaci�n se llevar� a cabo con la funci�n `normalizeCtData` del paquete `HTqPCR`.

 * **M�todo 1: $\Delta$ Cq**
 
Este m�todo de normalizaci�n calcula el $\Delta$ Cq seg�n describieron Livak y Schmittgen @livak2001analysis. El argumento `deltaCt` del paquete `HTqPCR` calcula el nivel de expresi�n promedio (Cq) de los genes de referencia escogidos en el apartado anterior, y calcula la diferencia de este promedio con respecto al resto de valores de expresi�n g�nica.

Si se escoge este m�todo de normalizaci�n se toman como referencia el o los genes: `r genes_ref`, tal y como se indic� en el apartado anterior.


 * **M�todo 2:  Media geom�trica**
 
Otro de los m�todos disponibles es la normalizaci�n por media geom�trica. El argumento `geometric.mean` calcula el valor Cq promedio para cada muestra, y reescala todos los valores Cq de la muestra con respecto a ese promedio. Algunos estudios han demostrado que este m�todo es el id�neo para normalizar expresi�n de miRNA @mestdagh2009novel.


A continuaci�n, el usuario deber� especificar qu� m�todo de normalizaci�n prefiere. Por defecto se aplicar� el m�todo $\Delta$Cq.

```{r}
# MODIFICABLE POR EL USUARIO
tipo_de_normalizacion <- "deltaCt" # opciones: "geometric.mean", "deltaCt"

```


```{r normalizacion, echo=FALSE, include=FALSE}

# normalizaci�n
exp_norm <- normalizeCtData(raw, norm=tipo_de_normalizacion, deltaCt.genes = c(genes_ref))
```


## Resultado de la normalizaci�n

Es buena idea comparar los valores crudos con los valores normalizados, para visualizar c�mo ha ido la correcci�n de los datos. En el siguiente gr�fico se representan los datos de expresi�n cada gen en crudo (eje x) frente a los datos normalizados (eje y).

```{r pre_post_norm, echo=FALSE, results='hide'}
plot(exprs(raw), exprs(exp_norm), pch = 20, main = paste("Normalizaci�n por m�todo", tipo_de_normalizacion), col = rep(colores4))
dev.copy(png,'Results_norm_vs_raw.png')
dev.off()
```

En los siguientes gr�ficos de cajas se muestran los valores crudos y normalizados de cada muestra por separado. Se puede apreciar c�mo cambia la distribuci�n de expresiones en cada una de las muestras antes y despu�s de la normalizaci�n.

```{r bx_raw, echo=FALSE, results='hide'}
### BOXPLOT de cada paciente
plotCtBoxes(raw, col = colores4, 
            las = 2, 
            names = paste("muestra", c(1:28)),
            mar = c(5,5,1,1),
            main = "Muestras sin normalizar",
            )
legend("topright", cex = 0.4,
        legend = c("1A", "2A", "1B", "2B"),
        fill = c('lightcoral', 'peachpuff', 'lightblue4', 'lightblue')
)

dev.off()
```

```{r bx_norm, echo=FALSE, results='hide'}
plotCtBoxes(exp_norm, col = colores4, 
            las = 2, 
            names = paste("muestra", c(1:28)),
            mar = c(5,5,1,1),
            main = "Muestras normalizadas",
            )
legend("topright", cex = 0.4,
        legend = c("1A", "2A", "1B", "2B"),
        fill = c('lightcoral', 'peachpuff', 'lightblue4', 'lightblue')
)

dev.off()
```


# Filtrado de datos no experimentales

En la etapa de laboratorio, los controles se adicionan a la muestra de ARN en distintos pasos, y permiten monitorizar los procesos de extracci�n y amplificaci�n. Una vez comprobado que dichos procesos han ido bien, los controles no aportan ninguna informaci�n sobre la expresi�n g�nica de las muestras, y por tanto se habr�n de excluir.

Tambi�n eliminaremos los datos marcados como 'Unreliable' y/o 'Undetermined'. 

De forma autom�tica, adem�s, se eliminan las medidas de miRNA con m�s de 3 medidas marcadas como "Unreliable". El usuario puede cambiar este l�mite en el siguiente c�digo.

```{r}
#MODIFICABLE POR EL USUARIO
excluir <- 3
```


```{r, echo=FALSE, results='asis'}
# filtro de controles,valores y 'Unreliable' y/o 'Undetermined'
exp_Filt <- filterCtData(exp_norm, remove.type = "Control", remove.category = "Unreliable", n.category = excluir)
```


# Control de calidad tras la normalizaci�n

En este gr�fico 'heatmap' tenemos un resumen visual de c�mo se agrupan las muestras experimentales (horizontal), as� como de los valores de expresi�n relativa de cada una de los genes (vertical) en cada una de las muestras. Se guardar� como `Results_heatmap.png`.

```{r, echo=FALSE, results='hide'}
plotCtHeatmap(exp_Filt, dist = "pearson")
dev.copy(png,'Results_heatmap.png')
dev.off()
```

En el siguiente gr�fico se muestra un an�lisis PCA con las dos principales fuentes de variaci�n entre muestras (x e y). Si los puntos o muestras aparecen formando agrupaciones es indicativo de que existir�n diferencias experimentales entre dichas muestras. Se almacenar� en `Results_norm_PCA.png`.

```{r PCA_exp_Filt, echo=FALSE, results='hide'}
plotCtPCA(exp_Filt, features=FALSE, col = colores4)
legend("bottomright", cex = 0.8,
       legend = c("1A", "2A", "1B", "2B"),
       fill = c('lightcoral', 'peachpuff', 'lightblue4', 'lightblue'))
dev.copy(png,'Results_norm_PCA.png')
dev.off()
```


# Expressi�n diferencial

En este �ltimo apartado se calcular� la expresi�n g�nica relativa entre los `r length(unique(targets$Treatment))*length(unique(targets$Response))` grupos experimentales. 

Se realizar�n contrastes 2 a 2 de cada uno de los grupos experimentales mediante un test t. El usuario podr� escoger el nivel de significaci�n que deben alcanzar los genes diferencialmente expresados. Por defecto est� establecido en $\alpha$=0,05. Se calcular� un p-valor ajustado por el m�todo de Benjamini y Hochberg (@benjamini1995controlling), debido a las comparaciones m�ltiples que se realizar�n.

```{r}
# MODIFICABLE POR EL USUARIO
sig <- 0.05
```

A continuaci�n, los contrastes realizados son espec�ficos de los datos del proyecto, que contiene los grupos experimentales 1A, 1B, 2A y 2B.

* **Grupo 1A vs Grupo 1B**

```{r, echo=FALSE, include=FALSE}
DE.1A.1B.ttest <- ttestCtData(exp_Filt[,1:14], groups = targets$Group[1:14])
```

A continuaci�n se muestra una lista de los 15 genes m�s diferencialmente expresados entre estos grupos experimentales. Se indica si la expresi�n aumenta o se inhibe en relaci�n al primer grupo del contraste t, en este caso, en relaci�n al grupo 1A. La tabla de resultados se guardar� en `Results_DE_1Avs1B.csv`.

```{r, echo=FALSE}
# creamos nueva variable "REGULATED":
DE.1A.1B.ttest$Regulated <- NA

# recodificamos la variable:
DE.1A.1B.ttest$Regulated[DE.1A.1B.ttest$ddCt>0] <- "up"
DE.1A.1B.ttest$Regulated[DE.1A.1B.ttest$ddCt<0] <- "down"

pander(DE.1A.1B.ttest[1:15,c(1, 5, 6, 7, 12)])
write.csv(DE.1A.1B.ttest, file = "Results_DE_1Avs1B.csv")
```




A continuaci�n se muestran estos 15 genes en un gr�fico de expresi�n g�nica relativa.

```{r, echo=FALSE, results='hide'}
plotCtRQ(DE.1A.1B.ttest, genes = 1:15, cex.lab = 0.1)
dev.copy(png,'Results_1A1B.png')
dev.off()
```

En conclusi�n, presentan expresi�n diferencial significativa un total de `r length(DE.1A.1B.ttest$genes[DE.1A.1B.ttest$adj.p.value < sig])` genes`r if(length(DE.1A.1B.ttest$genes[DE.1A.1B.ttest$adj.p.value < sig]) > 0){paste(" que son:", DE.1A.1B.ttest$genes[DE.1A.1B.ttest$adj.p.value<sig])}`. El nivel de significaci�n est� fijado en $\alpha$ = `r sig`.

* **Grupo 2A vs Grupo 2B**

```{r include=FALSE}
DE.2A.2B.ttest <- ttestCtData(exp_Filt[,15:28], groups = targets$Group[15:28])
```

A continuaci�n se muestra una lista de los 15 genes m�s diferencialmente expresados entre estos grupos experimentales. Se indica si la expresi�n aumenta o se inhibe en relaci�n al primer grupo del contraste t, en este caso, en relaci�n al grupo 2A. La tabla de resultados se guardar� en `Results_DE_2Avs2B.csv`.

```{r echo=FALSE}
# creamos nueva variable "REGULATED":
DE.2A.2B.ttest$Regulated <- NA

# recodificamos la variable:
DE.2A.2B.ttest$Regulated[DE.2A.2B.ttest$ddCt>0] <- "up"
DE.2A.2B.ttest$Regulated[DE.2A.2B.ttest$ddCt<0] <- "down"

pander(DE.2A.2B.ttest[1:15,c(1, 5, 6, 7, 12)])
write.csv(DE.2A.2B.ttest, file = "Results_DE_2Avs2B.csv")
```


A continuaci�n se muestran estos 15 genes en un gr�fico de expresi�n g�nica relativa.

```{r, echo=FALSE, results='hide'}
plotCtRQ(DE.2A.2B.ttest, genes = 1:15, cex.lab = 0.1)
dev.copy(png,'Results_2A2B.png')
dev.off()
```

En conclusi�n, presentan expresi�n diferencial significativa un total de `r length(DE.2A.2B.ttest$genes[DE.2A.2B.ttest$adj.p.value < sig])` genes`r if(length(DE.2A.2B.ttest$genes[DE.2A.2B.ttest$adj.p.value < sig]) > 0){paste(" que son:", DE.2A.2B.ttest$genes[DE.2A.2B.ttest$adj.p.value<sig])}`. El nivel de significaci�n est� fijado en $\alpha$ = `r sig`.

* **Grupo 1A vs Grupo 2A**

```{r include=FALSE}
DE.1A.2A.ttest <- ttestCtData(exp_Filt[,seq(1, 27, 2)], groups = targets$Group[seq(1, 27, 2)])
```

A continuaci�n se muestra una lista de los 15 genes m�s diferencialmente expresados entre estos grupos experimentales. Se indica si la expresi�n aumenta o se inhibe en relaci�n al primer grupo del contraste t, en este caso, en relaci�n al grupo 1A. La tabla se guardar� en `Results_DE_1Avs2A.csv`.

```{r echo=FALSE}
# creamos nueva variable "REGULATED":
DE.1A.2A.ttest$Regulated <- NA

# recodificamos la variable:
DE.1A.2A.ttest$Regulated[DE.1A.2A.ttest$ddCt>0] <- "up"
DE.1A.2A.ttest$Regulated[DE.1A.2A.ttest$ddCt<0] <- "down"

pander(DE.1A.2A.ttest[1:15,c(1, 5, 6, 7, 12)])
write.csv(DE.1A.2A.ttest, file = "Results_DE_1Avs2A.csv")
```


A continuaci�n se muestran estos 15 genes en un gr�fico de expresi�n g�nica relativa.

```{r, echo=FALSE, results='hide'}
plotCtRQ(DE.1A.2A.ttest, genes = 1:15, cex.lab = 0.1)
dev.copy(png,'Results_1A2A.png')
dev.off()
```

En conclusi�n, presentan expresi�n diferencial significativa un total de `r length(DE.1A.2A.ttest$genes[DE.1A.2A.ttest$adj.p.value < sig])` genes`r if(length(DE.1A.2A.ttest$genes[DE.1A.2A.ttest$adj.p.value < sig]) > 0){paste(" que son:", DE.1A.2A.ttest$genes[DE.1A.2A.ttest$adj.p.value<sig])}`. El nivel de significaci�n est� fijado en $\alpha$ = `r sig`.

* **Grupo 1B vs Grupo 2B**

```{r include=FALSE}
DE.1B.2B.ttest <- ttestCtData(exp_Filt[,seq(2, 28, 2)], groups = targets$Group[seq(2, 28, 2)])
```
A continuaci�n se muestra una lista de los 15 genes m�s diferencialmente expresados entre estos grupos experimentales. Se indica si la expresi�n aumenta o se inhibe en relaci�n al primer grupo del contraste t, en este caso, en relaci�n al grupo 1A. La tabla se guardar� en `Results_DE_1Bvs2B.csv`.

```{r echo=FALSE}
# creamos nueva variable "REGULATED":
DE.1B.2B.ttest$Regulated <- NA

# recodificamos la variable:
DE.1B.2B.ttest$Regulated[DE.1B.2B.ttest$ddCt>0] <- "up"
DE.1B.2B.ttest$Regulated[DE.1B.2B.ttest$ddCt<0] <- "down"

pander(DE.1B.2B.ttest[1:15,c(1, 5, 6, 7, 12)])
write.csv(DE.1B.2B.ttest, file = "Results_DE_1Bvs2B.csv")
```


A continuaci�n se muestran estos 15 genes en un gr�fico de expresi�n g�nica relativa.

```{r, echo=FALSE, results='hide'}
plotCtRQ(DE.1B.2B.ttest, genes = 1:15, cex.lab = 0.1)
dev.copy(png,'Results_1B2B.png')
dev.off()
```

En conclusi�n, presentan expresi�n diferencial significativa un total de `r length(DE.1B.2B.ttest$genes[DE.1B.2B.ttest$adj.p.value < sig])` genes`r if(length(DE.1B.2B.ttest$genes[DE.1B.2B.ttest$adj.p.value < sig]) > 0){paste(" que son:", DE.1B.2B.ttest$genes[DE.1B.2B.ttest$adj.p.value<sig])}`. El nivel de significaci�n est� fijado en $\alpha$ = `r sig`.


# Informaci�n de la sesi�n

Como �ltima parte de este documento, se indica la versi�n de R con la que se ha ejecutado el pipeline, as� como los paquetes usados. Se guardar� autom�ticamente un archivo con la informaci�n de salida. �sta informaci�n puede ser de utilidad a la hora de diagnosticar problemas de ejecuci�n.

```{r sesioninfo, echo=FALSE}
sessionInfo()

writeLines(capture.output(sessionInfo()), "Results_sessionInfo.txt")
```

# Referencias