---
title: "Anotación"  
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
---

## Tipos de paquetes  

- **ChipDb**:  Relativos a una plataforma, a un chip concreto de microarray. 
- **OrgDb**: Se refieren a organismos.
- **TxDb**: Relativos a transcriptomas de un organismo.
- **BSgenome**: Genomas.


## AnnotationDbi  

- Las bases de datos de tipo **ChipDb**, **OrgDb** y **TxDb**
heredan todos los métodos de la clase **AnnotationDbi::AnnotationDb**.
- Podemos aplicar los métodos: **columns**, **keytypes**,**keys** y **select**.

- Cargamos **AnnotationDbi** y un paquete de anotación, por ejemplo,
**hgu133a.db**.  

```{r }
pacman::p_load("AnnotationDbi","hgu133a.db")
```

- La información contenida la podemos ver con

```{r }
ls("package:hgu133a.db")
```

Con el nombre del paquete también tenemos información.

```{r }
hgu133a.db
```

- En los paquetes de anotación tenemos **columns**.
Algunas de estas columnas pueden ser **keys**.
Podemos realizar consultas en la base de datos utilizando
una **key** y pedir que nos devuelva una o más de una
**columns**.

- ¿Qué información podemos recuperar utilizando **select**?

```{r }
columns(hgu133a.db)
```
- Pero: ¿qué información es? 

```{r ,results='hide'}
help("ENTREZID")
```

- No todas las variables que hemos obtenido con **columns** son utilizables
para realizar consultas.

- Aquellas utilizables para las consultas las podemos conocer con
**keytypes**.

- A estas variables las llamamos llaves **keys**.

```{r }
keytypes(hgu133a.db)
```
- ¿Cómo conseguir todas los valores de una llave determinada?

```{r }
head(keys(hgu133a.db,keytype="ENTREZID"))
head(keys(hgu133a.db,keytype="ENSEMBL"))
```

- Supongamos que tenemos algunos identificadores Affymetrix (**PROBEID**) y 
pretendemos conocer sus identificadores **ENTREZID**.

```{r }
(ids = sample(keys(hgu133a.db,keytype="PROBEID"),5))
AnnotationDbi::select(hgu133a.db,keys=ids,columns=c("ENTREZID","ENSEMBL",
                                              "SYMBOL"),keytype="PROBEID")
```


## ChipDb
- Si trabajamos con microarrays nuestras características serán las sondas.
- Hemos de poder hacer corresponder este identificador
con el gen al que corresponde.

```{r }
pacman::p_load(hgu133a.db)
```

¿Qué sondas tienen correspondencia en Entrez?

```{r }
mappedProbes = mappedkeys(hgu133aENTREZID)
```

- Lo guardamos en forma de lista.  

```{r }
mappedProbesList = as.list(hgu133aENTREZID[mappedProbes])
```
Por ejemplo, la primera posición de la lista nos da el identificador de 
Affymetrix y su identificador Entrez.

```{r }
mappedProbesList[1]
```

- Podemos ver todas las correspondencias que nos ofrece el paquete. 

```{r }
ls("package:hgu133a.db")
```

- Otro ejemplo de determinar la correspondencia entre PROBEID e identificadores.

```{r }
ids = c("39730_at", "1635_at", "1674_at", "40504_at", "40202_at")
pacman::p_load("hgu95av2.db")
columns(hgu95av2.db)
keytypes(hgu95av2.db)
columns = c("PFAM","SYMBOL")
AnnotationDbi::select(hgu95av2.db, keys=ids, columns, keytype="PROBEID")
```

## OrgDb

```{r }
pacman::p_load(org.Hs.eg.db)
```

- ¿Qué tipo de cosas o qué claves podemos manejar?

```{r }
columns(org.Hs.eg.db)
keytypes(org.Hs.eg.db)
```

- Los primeros identificadores de genes utilizando el código ENTREZID. 

```{r }
head(keys(org.Hs.eg.db, keytype="ENTREZID"))
```

- Los mismos pero en ENSEMBL.

```{r }
head(keys(org.Hs.eg.db, keytype="ENSEMBL"))
```
- Y en Gene Ontology.

```{r }
head(keys(org.Hs.eg.db, keytype="GO"))
```

- Supongamos que elegimos un sistema de identificación, por ejemplo,
ENTREZID.

```{r }
(ids = keys(org.Hs.eg.db, keytype="ENTREZID")[1:5])
```

- A partir de estos identicadores podemos obtener el resto.

```{r }
AnnotationDbi::select(org.Hs.eg.db,keys=ids,column="SYMBOL",keytype='ENTREZID')
```

- Supongamos que nos fijamos en el gen con código **ENSG00000000003**
en Ensembl. 

```{r }
id = "ENSG00000000003"
```
- Buscamos su correspondencia en Gene Ontology.

```{r }
(res =  AnnotationDbi::select(org.Hs.eg.db, keys=id, columns="GO", keytype="ENSEMBL"))
```

- Como vemos no tenemos una correspondencia 1-1.
- Al mismo gen le corresponden distintos términos GO.
- Si solamente tenemos interés en ellos podemos hacer
```{r  }
res[,"GO"]
```

- Puesto que tenemos identificadores **GO** podemos utilizar el 
paquete **GO.db** para obtener los términos **GO** correspondientes.

```{r  }
pacman::p_load("GO.db")
```
Y los términos GO serían
```{r  }
AnnotationDbi::select(GO.db, keys=res[,"GO"], columns="TERM", keytype="GOID")
```

##  TxDb  

- Los paquetes TxDb están centrados en el genoma. 

```{r  }
pacman::p_load("TxDb.Dmelanogaster.UCSC.dm3.ensGene")
txdb = TxDb.Dmelanogaster.UCSC.dm3.ensGene
```
- ¿De qué clase es este objeto?

```{r  }
class(txdb)
```

- Las opciones que tenemos

```{r  }
columns(txdb)
keytypes(txdb)
```
- Para el manejo de este tipo de bases de datos es útil **GenomicFeatures**.

```{r  }
pacman::p_load(GenomicFeatures)
```
- Por ejemplo: ¿Qué cromosomas tenemos?

```{r  }
seqlevels(txdb)
```

- Cuando se carga la base de datos todos los cromosomas están activos
y lo que hagamos nos dará información sobre todos ellos.  
- Queremos, por ejemplo, trabajar solamente con **chr2L**.  

```{r  }
seqlevels(txdb) = "chr2L"
```

- Supongamos que queremos conocer los **GENEID** de los primeros genes.

```{r  }
(keysGENEID = head(keys(txdb, keytype="GENEID"),n=3))
```

```{r  }
columns = c("TXNAME", "TXSTART","TXEND","TXSTRAND")
AnnotationDbi::select(txdb, keysGENEID, columns, keytype="GENEID")
```

- Los objetos **TxDb** nos permiten obtener las anotaciones como 
**GRanges**.

```{r  }
(txdb.tr = transcripts(txdb))
```

- También podemos obtener información sobre los exones.

```{r  }
(txdb.ex = exons(txdb))
```

- Podemos incluir metadatos adicionales como puede ser el identificador del gen.

```{r  }
transcripts(txdb, columns = c("tx_id","tx_name","gene_id"))
```
- Obtenemos las regiones CDS con

```{r  }
(txdb.cds = cds(txdb))
```

- En los tres casos hemos obtenido un objeto **GRanges**.

- A este tipo de objetos podemos aplicar otros métodos que nos dan información
adicional.

```{r  }
length(txdb.cds)
strand(txdb.cds)
```

- A partir de un objeto de clase  **TxDb** podemos obtener un 
**GRangesList** en la cual separamos por alguna característica.
- Por ejemplo, los objetos **GRanges** para los distintos genes.

```{r  }
transcriptsBy(txdb, by="gene")
```
- Podemos agrupar los exones por gen.

```{r  }
exonsBy(txdb, by="gene")
```
- Podemos tener los CDS agrupados por transcrito.

```{r  }
cdsBy(txdb, by="tx")
```

- Los intrones agrupados por transcrito.

```{r  }
intronsByTranscript(txdb)
```

Las regiones UTR 5' y 3' agrupadas por transcrito.

```{r  }
fiveUTRsByTranscript(txdb)
threeUTRsByTranscript(txdb)
```

## BSgenome  

```{r  }
pacman::p_load(BSgenome.Dmelanogaster.UCSC.dm3)
tx2seqs = extractTranscriptSeqs(BSgenome.Dmelanogaster.UCSC.dm3,TxDb.Dmelanogaster.UCSC.dm3.ensGene)
```

- La secuencia correspondiente al primer gen sería

```{r  }
tx2seqs[[1]]
```
- Si queremos conocer la secuencia entre las posiciones 1000 y 1020 entonces

```{r  }
tx2seqs[[1]][1000:1020]
```
- También podemos trasladar estas secuencias a proteínas con

```{r  }
suppressWarnings(translate(tx2seqs[[1]]))
```

- Para todos los transcritos lo hacemos con

```{r  }
suppressWarnings(translate(tx2seqs))
```

- No todo lo que se transcribe se traslada. Para obtener obtener las que 
realmente se trasladan se puede hacer

```{r  }
cds2seqs = extractTranscriptSeqs(BSgenome.Dmelanogaster.UCSC.dm3,cdsBy(txdb, by="tx"))
translate(cds2seqs)
```


##  OrganismDb   

- Un paquete de tipo OrganismDb nos permite combinar información
(para el organismo con que estemos trabajando) de **GO.db** con el correspondiente
**TxDb** y **OrgDb**.

```{r  }
pacman::p_load(Homo.sapiens)
```

```{r  }
(keys = head(keys(Homo.sapiens, keytype="ENTREZID"), n=2))
columns = c("SYMBOL","TXNAME")
AnnotationDbi::select(Homo.sapiens, keys, columns, keytype="ENTREZID")
```
- También podemos obtener \verb|GRanges|.
```{r  }
transcripts(Homo.sapiens, columns=c("TXNAME","SYMBOL"))
```

## biomaRt  

- El paquete **biomaRt** es un interfaz para poder acceder a 
una serie de bases de datos que implementan [BioMart](http://www.biomart.org "BioMart")

```{r  }
pacman::p_load(biomaRt)
```

- Podemos ver la lista de **mart**'s que tenemos

```{r ,eval=FALSE }
listMarts(host="www.ensembl.org")
```

- Elegimos utilizar **ensembl**.

```{r  ,eval=FALSE }
(ensembl = useMart("ENSEMBL_MART_ENSEMBL",host="www.ensembl.org"))
```
- Ahora hemos de elegir el conjunto de datos a utilizar. 
- Podemos ver los disponibles con 

```{r   ,eval=FALSE}
head(listDatasets(ensembl))
```
- Elegimos la correspondiente a ser humano.

```{r  ,eval=FALSE }
(ensembl = useMart("ENSEMBL_MART_ENSEMBL",dataset="hsapiens_gene_ensembl",
     host="www.ensembl.org"))
```
- Podemos ver los atributos con

```{r   ,eval=FALSE}
head(listAttributes(ensembl))
```
- De hecho, son 
```{r  ,eval=FALSE }
nrow(listAttributes(ensembl))
```
- Podemos comprobar que hay muchos que identifican el gen con las sondas
de Affymetrix, en concreto, con los AffyID.

```{r   ,eval=FALSE}
affyids=c("202763_at","209310_s_at","207500_at")
```

```{r  ,eval = FALSE}
getBM(attributes=c('affy_hg_u133_plus_2', 'entrezgene'), 
                    filters = 'affy_hg_u133_plus_2', 
                    values = affyids, mart = ensembl)
```
- Podemos obtener todos los identificadores.

```{r  ,eval = FALSE}
head(getBM(attributes='affy_hg_u133_plus_2', mart = ensembl))
```

