trabajando con Gene Ontologies (GOs) he empezado a hacer unas pruebas con una biblioteca de Bioconductor llamada topGO. Esta sirve para comprobar si los términos de GOs de un conjunto de genes son diferentes de los que encontramos en otro conjunto.
Por ejemplo, si tenemos los genes de un organismo y nos interesan los que en nuestro tratamiento aumentan su expresión, podemos comprobar si los términos GO asociados a estos últimos son diferentes de los que tenemos para el conjunto de todos los genes del organismo.
Pero uno de los pasos previos, antes de hacer tests estadísticos con topGO, es asociar nuestros genes a términos GO. En Bioconductor y otras bases de datos hay listas ya precalculadas que podemos utilizar. Sin embargo, por distintos motivos (por ejemplo, al trabajar con un organismo para el cual no está tan estudiada la anotación funcional), podemos querer cargar nuestras propias asociaciones para que topGO trabaje sobre ellas.
Ya que un gen puede tener asociado más de un término GO, los desarrolladores de topGO han decidido que la estructura de datos a utilizar sea una lista de vectores de caracteres (técnicamente una "named list of character vectors"). Dicho de otra forma, una lista de genes, con cada una de las entradas de la lista asociada a un vector de nombres de términos GO. Una opción, si queremos trabajar con nuestra propia anotación, es utilizar la función 'readMappings' de topGO, que crea esta estructura si le damos un formato de fichero adecuado:
"It consists of one line for each gene with the following syntax:
gene_ID/TAB/GO_ID1, GO_ID2, GO_ID3, ...."
La otra opción es que creemos nosotros mismos la lista de vectores a partir de los datos que queramos. A veces, cuesta menos transformar un fichero de datos a otro formato (y usar la función 'readMappings', por ejemplo), y otras es más conveniente mantener nuestro formato de fichero y manipularlo como estructura de datos, sin crear otro fichero. En este caso, es además una buena oportunidad para probar algunas estructuras de datos y sentencias de R, que es lo que haremos nosotros.
¿Cómo es el formato de datos que tenemos en el fichero inicial? Obviamente no es el formato requerido por topGO (pues usaríamos la función 'readMappings' sin más). Nosotros vamos a partir de una tabla con 2 columnas (separadas por tabulador): una para el identificador del gen, la otra para un identificador de término GO. De esta forma, los distintos términos GO de cada gen quedan realmente en varias líneas:
gene_1 GO:000001
gene_1 GO:000004
gene_2 GO:000003
...
Esto obviamente puede cambiar según el origen de los datos, con qué programa se hayan obtenido las anotaciones, etc.
No vamos a entrar en detalles aquí sobre cómo leer el fichero con R o en dar ejemplos concretos de los datos. Es muy fácil generarlos y no es el objetivo de esta entrada. Lo importante para entender el código es que vamos a trabajar con la salida de read.table, en la variable 'go_tab', que será de tipo data.frame. En 'go_tab', los campos de las dos columnas se llaman 'identifier' (la columna con genes) y 'GOterm' (la columna con términos GO).
Vamos a ver varias alternativas para transformar éste 'data.frame' en una lista de vectores. En primer lugar, probaremos creando directamente una lista que iremos poblando con un bucle 'for'. Éste código es muy legible e intuitivo y programadores de muchos lenguajes diferentes pueden entenderlo fácilmente con solo echarle un vistazo.
El código completo de éste fragmento:
gene2GO = list()
for (i in seq(1:nrow(go_tab))) {
identifier = go_tab$identifier[i]
go_term = go_tab$GOterm[i]
if (identifier %in% names(gene2GO)){
gene2GO[[identifier]] = c(gene2GO[[identifier]], as.character(go_term))
} else {
gene2GO[[identifier]] = as.character(go_term)
}
}
En éste ejemplo, por tanto, vamos a recorrer la tabla de datos (go_tab) importada del fichero, e iremos construyendo nuestra lista final (gene2GO) concatenando nuevos términos GO para cada gen:
gene2GO[[identifier]] = c(gene2GO[[identifier]], as.character(go_term))
Vemos aquí una de las peculiaridades de R en el manejo de listas, que es el uso de corchetes dobles '[['. Si hay dudas de por qué se usa aquí '[[' en lugar de '[', la explicación rápida es que en R para acceder a datos en una lista hay que usar '[['. Para ampliar esto, está detallado en http://adv-r.had.co.nz/Subsetting.html
Sin embargo, éste método es bastante lento en mi equipo (Dell Optiplex 9010) y con mis datos (unas 90.000 filas en el fichero de entrada). Pongamos, por referencia y para comparar con métodos que veremos a continuación, que éste ejemplo en el que usamos una lista y un bucle 'for' tarda 1' (54'' en promedio, realmente).
Parte de por qué es así puede tener que ver con el bucle 'for'. Más adelante veremos alguna alternativa al bucle. Sin embargo, también tiene que ver el rendimiento de la búsqueda de los elementos en la lista.
Una alternativa al operador '%in%' podría ser:
if (exists(identifier, gene2GO)){
Pero en principio esta opción es más lenta, en mi caso unos 2' (1'57'').
En ambos casos, la comprobación la hacemos para acceder después al gen en cuestión. Por tanto, se está buscando 2 veces el mismo elemento en la lista: una en la comprobación, otra al obtener el elemento. Una forma de evitar esto sería cogiendo directamente el objeto (1 sola búsqueda) y comprobando si existe o no (NULL):
current = gene2GO[[identifier]]
if (is.null(current)){
(Nota: ojo que ahora el if comprueba si NO existe, no si existe como en los ejemplos anteriores).
Vemos que con esta otra forma tenemos un tiempo promedio de unos 26''. Aproximadamente la mitad que con el uso de '%in%', como esperaríamos.
Sin embargo, ya hemos comentado antes que el uso de un bucle 'for' en R también puede tener algo que ver con la velocidad de nuestro código (por ejemplo ver: http://faculty.nps.edu/sebuttre/home/R/apply.html). En R se suele preferir utilizar funciones que van recorriendo las estructuras de datos y aplicando una función sobre cada elemento de la estructura. Es el caso de la familia de funciones 'apply' (ver por ejemplo https://nsaunders.wordpress.com/2010/08/20/a-brief-introduction-to-apply-in-r/). Veamos el equivalente al código anterior usando la función 'apply':
gene2GO = list()
add_to_list <- function(go_tab){
identifier = go_tab[[1]]
go_term = go_tab[[2]]
current = gene2GO[[identifier]]
if (is.null(current)){
gene2GO[[identifier]] <<- as.character(go_term)
} else {
gene2GO[[identifier]] <<- c(gene2GO[[identifier]], as.character(go_term))
}
}
invisible(apply(go_tab, 1, add_to_list))
Ahora está tardando unos 21'' en mi máquina. Una ligera mejora respecto al uso del bucle 'for', que podría tener impacto en conjuntos de datos muy grandes o cuando queremos correr muchos análisis.
Vemos varias peculiaridades en el código anterior. La función 'apply' va a recorrer nuestra lista y llamará a la función 'add_to_list' para cada elemento. Sin embargo, en nuestro caso no queremos simplemente transformar los elementos de la lista (por ejemplo, formateando el nombre del gen), si no que queremos poblar una estructura de datos diferente a partir de los valores de la lista recorrida. Por eso se utiliza el operador <<-:
"The operators are normally only used in functions, and cause a search to made through parent environments for an existing definition of the variable being assigned."
Otra posible opción sería utilizar un parámetro opcional en la función, pero no me queda claro que vaya a funcionar para escribir en él, ya que el parámetro (la lista de vectores que estamos creando) necesitaría pasarse por referencia, y me suena que en R no he pasado parámetros así anteriormente ¿Os suena? ¿Algún consejo sobre esto? ¿Alguna otra opción?
La otra peculiaridad es la función 'invisible'. Ya que las funciones en R devuelven un valor y 'apply' lo propaga para obtener la posible transformación de 'go_tab' en otra cosa, si no asignamos el resultado de 'apply' a una variable tendremos el resultado en nuestra salida. Para evitar esto, podemos asignar el resultado a una variable "dummy" o utilizar el método 'invisible', que a mi me ha parecido más adecuado y legible.
Otra diferencia, más básica, es cómo accedemos a nuestros datos originales. En el primer ejemplo, lo que hacíamos es acceder a los campos del 'data.frame' y, en un campo dado, obtener el dato correspondiente según el bucle for:
identifier = go_tab$identifier[i]
Con 'apply' en cambio estamos recibiendo cada una de las líneas del 'data.frame' como un 'vector', por lo que accedemos directamente a los datos en éste vector:
identifier = go_tab[[1]]
Por ahora hemos mejorado nuestro algoritmo cambiando la forma de comprobar la existencia de un elemento en la lista para luego obtenerlo, y pasando de un bucle 'for' a usar la función 'apply'. Pero aún podemos ir más allá. ¿Por qué utilizar una lista? Quizás podamos mejorar el rendimiento utilizando un entorno con un mapeado "hash" (https://stat.ethz.ch/R-manual/R-devel/library/base/html/environment.html).
gene2GO = new.env()
add_to_list <- function(go_tab){
identifier = go_tab[[1]]
go_term = go_tab[[2]]
current = gene2GO[[identifier]]
if (is.null(current)){
gene2GO[[identifier]] <<- as.character(go_term)
} else {
gene2GO[[identifier]] <<- c(gene2GO[[identifier]], as.character(go_term))
}
}
invisible(apply(go_tab, 1, add_to_list))
gene2GO = as.list(gene2GO)
Aquí, en lugar de una lista creamos un entorno ('environment'), que por defecto utiliza una función "hash" para mapear su contenido. Por lo demás, como vemos, lo estamos manejando aquí de forma muy similar a la lista, aunque hemos añadido al final una conversión a formato de lista (en una sola línea pasamos de nuestro "hash" a una lista de vectores).
Éste código tarda unos 7'', incluída la conversión. Hay que tener en cuenta que el uso de un "hash" puede ser más o menos aconsejable según el volumen de datos con el que se trabaje. Desde luego, con nuestros datos, hemos mejorado bastante el rendimiento prestando atención a distintas alternativas que ofrece R y cuál sería la mejor combinación de las que hemos probado.
Y ya tendríamos nuestra lista creada. Nuestra intención era utilizar esta lista para informar a topGO de qué términos de GOs se relacionan con nuestros genes ¿recordáis? Le pasaremos esta lista a topGO como parámetro 'gene2GO' en la función en la creamos el primer objeto necesario con éste paquete, e indicaremos que la función de anotación es del tipo "annFun.gene2GO". Para eso, y terminar así la preparación de los datos, tendríamos que importar las listas de genes a comparar, o una sola lista con todos pero con algún criterio que permita diferenciar un subconjunto de otro. Crearíamos entonces el objeto con:
GOdata <- new("topGOdata", description="whatever", ontology = "BP", allGenes = background, geneSel = test,
annot = annFUN.gene2GO, gene2GO = gene2GO, nodeSize = 10)
Donde:
'background' es nuestro conjunto de genes de referencia.
'test' es el grupo de genes que estamos interesados en comparar con la referencia.
'gene2GO' la lista que hemos generado.
Seguiríamos entonces con las siguientes fases de topGO: los test de enriquecimiento y el análisis de los resultados. Si te ha entrado el gusanillo de probar, topGO está perfectamente documentado:
http://www.bioconductor.org/packages/release/bioc/vignettes/topGO/inst/doc/topGO.pdf
Un saludo!
http://www.bioconductor.org/packages/release/bioc/vignettes/topGO/inst/doc/topGO.pdf
Un saludo!
