Primero debemos instalar el paquete “biomaRt” y el paquete “topGO” con el instalador “BiocManager”.
if (!require("BiocManager", quietly = TRUE)) # Instalamos SI FUERA NECESARIO
install.packages("BiocManager") # el instalador de bioconductor
BiocManager::install("biomaRt") # Instalamos los paquetes usando
BiocManager::install("topGO") # el instalador de bioconductor
BiocManager::install("GOstats", force = TRUE)Como veremos mas adelante, una alternativa de usar la base de datos en linea que es “biomaRt”, podemos usar una familia de bases de datos descargables, mucho mas pequeñas pero suficientemente útiles para nuestros análisis como son los paquetes “org.XX.eg/ipi.db”. Este nombre define de que especie queremos descargar la información (“XX”) y si queremos descargar la base de datos de genes (“eg”) o proteínas (“ipi”)
BiocManager::install("org.Hs.eg.db")library(biomaRt)
library(topGO)
library(GOstats)
library(org.Hs.eg.db)Este paquete nos ofrece la oportunidad de enlazarnos y descargar la
bases de datos de los genomas de varias especies, distribuidas en tablas
de dos dimensiones donde las filas corresponden a los genes y las
columnas a diferentes caracteristicas de los genes, como longitud,
secuencia, terminos GO asociados, sitio de inicio, sitio de termino,
entre varios otros. Para saber a que bases de datos que podemos acceder
con este paquete son cuatro y los podemos visualizar con la función
biomaRt::listMarts()
biomaRt::listMarts()## biomart version
## 1 ENSEMBL_MART_ENSEMBL Ensembl Genes 109
## 2 ENSEMBL_MART_MOUSE Mouse strains 109
## 3 ENSEMBL_MART_SNP Ensembl Variation 109
## 4 ENSEMBL_MART_FUNCGEN Ensembl Regulation 109Podemos acceder a estas mismas bases de datos en un navagafor
directamente a la base de datos relacional e!Ensembl.
Para esta clase vamos a trabajar con la primera base de datos llamada
“ENSEMBL_MART_ENSEMBL” (o simplemente
“ensembl”). Para selecionar esta base de datos usaremos
la función biomaRt::useMart(), y como parámetro
mart= el nombre de la base de datos.
mart <- useMart("ensembl") # definimos la base de datosSi bien es cierto, acabomos de definir la base de datos a la que
queremos acceder, esta por si misma hace referencia o contiene multiples
“datasets” como por ejemplo el del ser humano
(Homo sapiens) o zebra fish o delfin entre otros. Para
saber que “datasets” están disponibles podemos usar la
función biomaRt::listDatasets() del objeto
mart en el que fue definida la base de datos para ver
una lista de los datasets.
knitr::kable(head(biomaRt::listDatasets(mart = mart)))| dataset | description | version |
|---|---|---|
| abrachyrhynchus_gene_ensembl | Pink-footed goose genes (ASM259213v1) | ASM259213v1 |
| acalliptera_gene_ensembl | Eastern happy genes (fAstCal1.2) | fAstCal1.2 |
| acarolinensis_gene_ensembl | Green anole genes (AnoCar2.0v2) | AnoCar2.0v2 |
| acchrysaetos_gene_ensembl | Golden eagle genes (bAquChr1.2) | bAquChr1.2 |
| acitrinellus_gene_ensembl | Midas cichlid genes (Midas_v5) | Midas_v5 |
| amelanoleuca_gene_ensembl | Giant panda genes (ASM200744v2) | ASM200744v2 |
Ahora que ya sabemos el contenido de datasets en mi base de datos, ahora puedo seleccionar el que deseamos, en nuestro caso el de Homo sapiens y definirlo en nuestro objeto mart
mart <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")Para ello, debemos seguir los siguientes pasos
biomaRt::useMart()biomaRt::getBM()# Paso 1:
mart <- useMart("ensembl", dataset = "hsapiens_gene_ensembl") # seleccionamos la de Hs
# Paso 2:
immunome<-read.csv("/Users/alfredocardenasrivera/Project/bioinformatics_2022/human_gene_set.txt")
# Paso 3:
goannot<-biomaRt::getBM(attributes = c("entrezgene_id", "go_id", "name_1006"),
filters="entrezgene_id", values=immunome, mart=mart)El output o el resultado de la función biomaRt::getBM()
es un sencillo data.frame, en el cual las columnas son las variables que
nosotros definimos en el parametro “attributes” de la
funci´øn biomaRt::getBM(). En este caso deberían ser 3:
Podemos explorar este data.frame con funciones basicas como
str(goannot)## 'data.frame': 27144 obs. of 3 variables:
## $ entrezgene_id: int 100 100 100 100 100 100 100 100 100 100 ...
## $ go_id : chr "" "GO:0019239" "GO:0009168" "GO:0005886" ...
## $ name_1006 : chr "" "deaminase activity" "purine ribonucleoside monophosphate biosynthetic process" "plasma membrane" ...summary(goannot)## entrezgene_id go_id name_1006
## Min. : 2 Length:27144 Length:27144
## 1st Qu.: 2534 Class :character Class :character
## Median : 4598 Mode :character Mode :character
## Mean : 180143
## 3rd Qu.: 7520
## Max. :100133941knitr::kable(head(goannot,20))| entrezgene_id | go_id | name_1006 |
|---|---|---|
| 100 | ||
| 100 | GO:0019239 | deaminase activity |
| 100 | GO:0009168 | purine ribonucleoside monophosphate biosynthetic process |
| 100 | GO:0005886 | plasma membrane |
| 100 | GO:0005829 | cytosol |
| 100 | GO:0070161 | anchoring junction |
| 100 | GO:0004000 | adenosine deaminase activity |
| 100 | GO:0046936 | 2’-deoxyadenosine deaminase activity |
| 100 | GO:0016814 | hydrolase activity, acting on carbon-nitrogen (but not peptide) bonds, in cyclic amidines |
| 100 | GO:0060205 | cytoplasmic vesicle lumen |
| 100 | GO:0005737 | cytoplasm |
| 100 | GO:0016020 | membrane |
| 100 | GO:0031410 | cytoplasmic vesicle |
| 100 | GO:0046872 | metal ion binding |
| 100 | GO:0008270 | zinc ion binding |
| 100 | GO:0016787 | hydrolase activity |
| 100 | GO:0007155 | cell adhesion |
| 100 | GO:0005764 | lysosome |
| 100 | GO:0009986 | cell surface |
| 100 | GO:0009117 | nucleotide metabolic process |
# knitr::kable(head(table(goannot$entrezgene_id, goannot$go_id)))Con esta data.frame ya podemos hacer los análisis de enriquecimiento como el de ORA, con el paquete topGO
Cuando la lista de genes es muy grande y estamos trabajamos con genomas este trabajo se pude ver enlentecido (especialemnte cuando el ancho de banda de internet no es bueno). Tenemos alternativas a la base de datos en línea de biomaRt como son la base de datos de “AnnotationDbi”, que esta dividida en dos grupos:
Notese que los nombres de las bases de datos AnnotationDbi tiene una estructura particular (separada por puntos). Tenemos cuatro partes
Para construir este data.frame vamos a utilizar la función
biomaRt::select(). Esta función tiene 4 parámetros:
Para saber que columnas tiene nuestra base de datos, podemos usar la
función biomaRt::columns()) y para saber cuales son los
keytype disponibles podemos usar la función
biomaRt::keytypes()
biomaRt::columns(org.Hs.eg.db)## [1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS"
## [6] "ENTREZID" "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME"
## [11] "GENETYPE" "GO" "GOALL" "IPI" "MAP"
## [16] "OMIM" "ONTOLOGY" "ONTOLOGYALL" "PATH" "PFAM"
## [21] "PMID" "PROSITE" "REFSEQ" "SYMBOL" "UCSCKG"
## [26] "UNIPROT"biomaRt::keytypes(org.Hs.eg.db)## [1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS"
## [6] "ENTREZID" "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME"
## [11] "GENETYPE" "GO" "GOALL" "IPI" "MAP"
## [16] "OMIM" "ONTOLOGY" "ONTOLOGYALL" "PATH" "PFAM"
## [21] "PMID" "PROSITE" "REFSEQ" "SYMBOL" "UCSCKG"
## [26] "UNIPROT"Ahora, con esta información, ya podemos llamar/ejecutar la función
biomaRt::select(), sabiendo el calor para cada uno de los
parámetros
goannot2 <- biomaRt::select(x = org.Hs.eg.db,
keys = as.character(immunome$EntrezGeneID),
columns = "GO",
keytype = "ENTREZID")## 'select()' returned 1:many mapping between keys and columnsLe debe aparecer una ADVERTENCIA, como
'select()' returned 1:many mapping between keys and columns.
Esto se debe a que cuando se hizo la asociación no se hizo uno a uno,
sino que para cada elemento de nuestra lista de interes se le asoció
varios elementos de la base de datos org.Hs.eg.db. Podemos revisar como
está estructurado el output con algunas funciones básicas de R
class(goannot2)## [1] "data.frame"str(goannot2)## 'data.frame': 26757 obs. of 4 variables:
## $ ENTREZID: chr "100" "100" "100" "100" ...
## $ GO : chr "GO:0000255" "GO:0001666" "GO:0001829" "GO:0001889" ...
## $ EVIDENCE: chr "IEA" "IDA" "IEA" "IEA" ...
## $ ONTOLOGY: chr "BP" "BP" "BP" "BP" ...knitr::kable(head(goannot2))| ENTREZID | GO | EVIDENCE | ONTOLOGY |
|---|---|---|---|
| 100 | GO:0000255 | IEA | BP |
| 100 | GO:0001666 | IDA | BP |
| 100 | GO:0001829 | IEA | BP |
| 100 | GO:0001889 | IEA | BP |
| 100 | GO:0001890 | IEA | BP |
| 100 | GO:0002314 | IEA | BP |
El resultado es un data frame de casi 27 filas, con 4 columnas: la primera es el ID que corrresponde al ENTREZGEN database, el segundo el código GO, el tercero es el nivel de evidencia y finalmente el tipo de termino GO que se define para cada línea. Los terminos GO tiene tres categorias:
Hasta ahora solo hemos asociado nuestra lista de genes de interes con los terminos GO. Pero también podemos hacer el análisis de enriquecimiento con la base de datos KEGG. Para ello vamos a usar la columna “PATH” de la base de datos “org.Hs.eg.db”.
pathannot <- biomaRt::select(x = org.Hs.eg.db,
keys = as.character(immunome$EntrezGeneID),
columns = "PATH",
keytype = "ENTREZID")## 'select()' returned 1:many mapping between keys and columnsSi la exploramosn veremos que es un data frame de 2723 filas y solo dos columnas, una para el ID de la base de datos ENTREZGENE y una para el ID de la base de datos KEGG.
class(pathannot)## [1] "data.frame"str(pathannot)## 'data.frame': 2723 obs. of 2 variables:
## $ ENTREZID: chr "100" "100" "100" "1003" ...
## $ PATH : chr "00230" "01100" "05340" "04514" ...knitr::kable(head(pathannot))| ENTREZID | PATH |
|---|---|
| 100 | 00230 |
| 100 | 01100 |
| 100 | 05340 |
| 1003 | 04514 |
| 1003 | 04670 |
| 10068 | NA |
De la misma forma que hemos asociado nuestra lista de genes de interes con los codigos de GO, o los ID para los caminos KEGG, podemos asociarlos con las otras columnas de la base de datos “org.Hs.eg.db”, como es el nombre del gen “GENENAME, codigos UNIPROT ”UNIPROT”, el número de acceso en la base de datos de genbank ACCNUM”, entre otros.
Primero definiremos nuestra lista de genes y nuestros genes blanco/objetivo
immunome <- read.csv("/Users/alfredocardenasrivera/Project/bioinformatics_2022/human_gene_set.txt")
target <- read.csv("/Users/alfredocardenasrivera/Project/bioinformatics_2022/human_gene_selected_set.txt")Ahora necesitamos usar un metodo constructor para crear un objeto de la clase “GOHyperGParams”
paramsGO <- new("GOHyperGParams",
geneIds = target$EntrezGeneID,
universeGeneIds =immunome$EntrezGeneID,
annotation = "org.Hs.eg.db",
ontology = "BP",
pvalueCutoff = 0.01,
conditional = FALSE,
testDirection = "over")Este objeto es el input de la función
GOstats::hyperGTest(). Con esta función podemos hacer el
analisis de enriquecimiento. Para ver el output, de una forma mas
amigable podemos usar la función summary y solo seleccionar
las primeros 10 resultados anidandolo dentro de la función
head.
over.GO <- GOstats::hyperGTest(paramsGO)
knitr::kable(head(summary(over.GO)))| GOBPID | Pvalue | OddsRatio | ExpCount | Count | Size | Term |
|---|---|---|---|---|---|---|
| GO:0045087 | 0 | 9.415862 | 53.10882 | 123 | 268 | innate immune response |
| GO:0098542 | 0 | 6.428837 | 65.79152 | 126 | 332 | defense response to other organism |
| GO:0050776 | 0 | 6.155285 | 54.89233 | 113 | 277 | regulation of immune response |
| GO:0043207 | 0 | 6.070538 | 76.69072 | 134 | 387 | response to external biotic stimulus |
| GO:0051707 | 0 | 6.070538 | 76.69072 | 134 | 387 | response to other organism |
| GO:0009607 | 0 | 5.922578 | 77.48339 | 134 | 391 | response to biotic stimulus |
Tambien podemos hacer el analisis de las rutas metabólicas a partir
del objeto “KEGGHyperGParams”. Para ello primero
usaremos el método constructor new y, luego, haremos el
analisis hipergeométrico.
paramsPATH <- new("KEGGHyperGParams",
geneIds = target$EntrezGeneID,
universeGeneIds =immunome$EntrezGeneID,
annotation = "org.Hs.eg.db",
pvalueCutoff = 0.01,
testDirection = "over")
over.PATH <- GOstats::hyperGTest(paramsPATH)
knitr::kable(head(summary(over.PATH)))| KEGGID | Pvalue | OddsRatio | ExpCount | Count | Size | Term |
|---|---|---|---|---|---|---|
| 04610 | 0.00e+00 | 19.080107 | 8.653465 | 31 | 38 | NA |
| 04620 | 0.00e+00 | 7.847222 | 10.475248 | 30 | 46 | NA |
| 05215 | 2.00e-07 | 24.232000 | 3.415842 | 13 | 15 | NA |
| 04622 | 1.80e-06 | 9.390244 | 4.782178 | 15 | 21 | NA |
| 04010 | 3.10e-06 | 5.496368 | 7.742574 | 20 | 34 | NA |
| 04650 | 1.78e-05 | 3.235736 | 15.257426 | 30 | 67 | NA |
Para utilizar la función runTest(), necesitamos editar
el data.frame resultado de la función select() a un objeto
topGOdata (para lo cual usaremos nuevamente la función
constructura new)
Para crear el objeto topGOdata se necesita que la
información esté en un formato específico. Para ello, debemos crear
basicamente dos objetos:
La condición es que ambos objetos deban coincidir.
Para crear ambos objetos podemos usar un loop for{}
entrezID <- c() # creamos un objeto vacío para guardar los ENTREZGENE
listGO <- c() # creamos un objeto vacío para guardar la lista de terminos GO
for (i in 1:nrow(immunome)) {
entrezID[i] <- immunome$EntrezGeneID[i] # llenamos el vector con los codigos
listGO[i] <- list(goannot2[goannot2$ENTREZID==entrezID[i], "GO"]) # llenamos la lista con los terminos
names(listGO)[i] <- immunome$EntrezGeneID[i]
}Lo que necesitamos crear es un “factor” de “unos y ceros”, en el cual se defina cual de los elementos del objeto lista listGO son los genes de nuestra muestra objetivo (target, los genes que estan diferencialmente expresados, nuestros genes de interes). En nuestro caso, estos genes están en el objeto target
factVect <- as.factor(immunome$EntrezGeneID%in%target$EntrezGeneID) # Creamos el vector y lo definimos como factor
levels(factVect) # verificamos los niveles## [1] "FALSE" "TRUE"levels(factVect) <- c(0,1) # cambiamos el nombre de los niveles
names(factVect) <- immunome$EntrezGeneID # ponemos nombre a cada posiciónAhora ya podemos crear el objeto “topGOdata”
GOdata.MF <- new("topGOdata", # definimos la clase del objeto
ontology = "MF", # definimos la categoria del GO term
description = "Practica de clase Enriquecimiento", # una breve descripcion
allGenes = factVect, # factor con los genes, con nombre en cada posición
annot = annFUN.gene2GO, # función que utilizamos para anotar el GO
gene2GO = listGO) # lista con los terminos GO agrupados para cada gen##
## Building most specific GOs .....## ( 866 GO terms found. )##
## Build GO DAG topology ..........## ( 1153 GO terms and 1444 relations. )##
## Annotating nodes ...............## ( 873 genes annotated to the GO terms. )Ya con el objeto