---
title: "Datos de RNA-Seq"
bibliography: "/home/gag/Nextcloud/BIBLIOGRAFIA/Bibliografia.bib"  
author: "Guillermo Ayala Gallego"  
date: "`r Sys.Date()`"  
address: Departamento de Estadística e I.O. Universidad de Valencia  
format:  
    revealjs:  
      scrollable: true
      echo: true
      embed-resources: true  
    pdf:
       toc: true
       number-sections: false
       colorlinks: true
---

# Introducción

## Flujo de trabajo


1. Diseño experimental.
2.  Protocolos de extracción del RNA. 
3. Preparación de las librerías. Se convierte el RNA en cDNA y se 
añaden los adaptadores para la secuenciación. 
4. Se secuencian las lecturas cDNA utilizando una plataforma de 
secuenciación.
5. **Alineamiento de las lecturas secuenciadas a un genoma de referencia.**
6. **Resumen del número de lecturas alineadas a una región.**
7. **Normalización de las muestras para eliminar diferencias técnicas en 
la preparación.** 
8. **Estudio estadístico de la expresión diferencial incluyendo
en lo posible un modelo.** 
9. Interpretación de los resultados desde el punto de vista biológico. 


## Objetivos

1. ¿Dónde podemos obtener datos de secuenciación? De otro modo:
	¿qué repositorios podemos utilizar?
2. ¿Cómo podemos realizar búsquedas en los metadatos con objeto de 
		obtener experimentos que tenga que ver con lo que me interesa?
3. ¿Cómo bajarlos?
4. ¿Cómo contar las lecturas alineadas sobre regiones genómicas (genes, exones
o $\ldots$)?

## Repositorios  

1. [NCBI SRA](http://www.ncbi.nlm.nih.gov/sra)
2. [EBI ENA](http://www.ebi.ac.uk/ena)
3. [DDBJ](http://www.ddbj.nig.ac.jp)

## Formato FASTA

- El formato **fasta** está basado en texto.
- Se utiliza representar secuencias bien de nucleótidos bien de aminoácidos.
- Tanto unos como otros son representados por una sola letra. 
- También tiene símbolos para representar un hueco (gap) o parada en la 
  traducción o bien que no se sabe el nucleótido o aminoácido. 
- Tiene una línea que comienza con el símbolo $>$ al que sigue una 
descripción de la secuencia.
- En la siguiente línea empieza la secuencia de bases o aminoácidos. Se 
recomienda que no tener más de 80 columnas y se pueden tener todas las filas
que se precisen.
- Un ejemplo 

    `>MCHU - Calmodulin - Human, rabbit, bovine, rat, and chicken  
    MADQLTEEQIAEFKEAFSLFDKDGDGTITTKELGTVMRSLGQNPTEAELQDMINEVDADGNGTIDF
    PEFLTMMARKMKDTDSEEEIREAFRVFDKDGNGYISAAELRHVMTNLGEKLTDEEVDEMIREADI*`


## Formato FASTQ
- El formato **fastq** es el más popular para datos de secuencias.
- Consiste de cuatro líneas por lectura.
- La primera que comienza con el carácter `@`
  y contiene el **nombre de la secuencia** con alguna descripción opcional
  de la misma.
- La segunda línea contiene la secuencia con las letras
  que correspondan dependiendo del tipo (nucleótidos, aminoácidos).
-  La tercera línea que comienza con $+$ contiene información opcional
  sobre la secuencia.
- La cuarta línea cuantifica la confianza o calidad en la determinación
  de cada base recogida en la segunda línea (**Phred**).
 
 - Un ejemplo

    `@SRR1293399.1 ILLUMINA-545855_0026_FC629BG:6:1:1022:5049 length=50  
    ACAGGGACGCCATCGAATCCGGATCNTNNNNNNNNNNNNANNNNNNNNNN`  
	
    `+SRR1293399.1 ILLUMINA-545855_0026_FC629BG:6:1:1022:5049 length=50  
    dee\edYcdc`bbY`S]bb_]]Ua^BBBBBBBBBBBBBBBBBBBBBBBBB`  

##  [Phred](https://en.wikipedia.org/wiki/Phred_base_calling)
- ¿Cómo se cuantifica la confianza o precisión?

	1. Phred asigna los picos de fluorescencia a una de las cuatro bases: 
	**base call**.
	2. $P$: probabilidad para una base dada de ser mal asignada o clasificada
	3. Se muestra:
	$$
	Q = - 10 \log_{10} P.
	$$
	4. Una probabilidad $P$ muy pequeña de clasificación incorrecta
	se traduce en un valor grande de $Q$.
	5. Se muestra el caracter ASCII que ocupa la posición $33+Q$.
     
	 
## Phred	 
	 
- Simulamos las probabilidades de clasificación incorrecta y generamos al 
azar las bases.
```{r}
x = sample(Biostrings::DNA_ALPHABET[1:4],10,replace=TRUE)
y = runif(10,min=0,max=.001)
names(y) = x
y
```  
- Calculamos los `Q` valores.  
```{r}
(Q = round(-10 * log10(y)))
```
- Codificación Sanger   
```{r}
intToUtf8(Q+33)
```

## Ejemplos  

- `tamidata::PRJNA266927`  
- `tamidata::PRJNA297664`  
- `tamidata::PRJNA297798`  

## `SummarizedExperiment`  

- Cargamos el paquete.
```{r}
#| message: false
library(SummarizedExperiment)
```
- Leemos unos datos de `tamidata`  
```{r}
data(PRJNA297664,package = "tamidata")  
```
- ¿Qué clase es?  
```{r}
class(PRJNA297664)
```
-  ¿Cuántos genes y muestras tenemos?  
```{r}
dim(PRJNA297664)
```
- La matriz de conteos es   
```{r}
#| results: false
head(assay(PRJNA297664))
```
- ¿Qué devuelve `assay`?
```{r}
class(assay(PRJNA297664))
head(which(apply(assay(PRJNA297664),1,sum) > 10))
```
- Variables fenotípicas
```{r}
colData(PRJNA297664)
```
- ¿Y es?
```{r}
class(colData(PRJNA297664))
```
- Nombres de las variables fenotípicas.
```{r}
names(colData(PRJNA297664))
```
- Y accedemos a los valores de la variable `treatment` con
```{r}
colData(PRJNA297664)[,"treatment"]
```
- Datos relativos a las filas o características o genes. 
```{r}
head(rownames(PRJNA297664))
```
- ¿Hay más? Pues sí.  
```{r}
head(rowData(PRJNA297664))  
```
- ¿Algo más? De hecho, mucho más. 

```{r}
rowRanges(PRJNA297664)[3]
```
