Trabajo Práctico 5 - Herramientas para el análisis de Ontologías de Genes y Análisis de redes de Interacción de Proteínas

Objetivos:

El objetivo de este trabajo práctico es aprender a extraer información biológica a partir de los resultados obtenidos de experimentos de RNA-Seq, mediante el uso de análisis de enriquecimiento de pathways o grupos de genes (gene sets).

Introducción:

Las estrategias de experimentación tradicionales en la investigación biológica típicamente estudian uno o unos pocos genes a la vez. Contrariamente, las técnicas de high throughput como la genómica, transcriptómica y la proteómica permiten medir simultáneamente, en una determinada condición, los cambios y la regulación de genes a nivel del genoma completo. Estas tecnologías dan como resultado largas listas de genes “interesantes”, que muestran una diferencia cuali-cuantitativa entre dos estados biológicos a comparar. Esa información debería indicar cuáles procesos celulares están afectados en las condiciones biológicas analizadas. Sin embargo, interpretar biológicamente dichas listas de cientos o miles de genes no es una tarea fácil. Para ayudarnos en esto, han aparecido nuevas herramientas bioinformáticas que disecan sistemáticamente largas listas de genes y nos cuentan cuales son los procesos biológicos enriquecidos en dichas listas, es decir, aquellos en donde los genes “interesantes” aparecen con una frecuencia mayor a la esperada por azar. En otras palabras, estos softwares mapean sistemáticamente una larga lista de genes “interesantes” a su anotación biológica (i.e. “términos” ontológicos) y luego destacan aquellos términos o anotaciones sobrerrepresentados en forma estadísticamente significativa. De esta manera es posible aumentar la probabilidad de identificar los procesos biológicos más pertinentes a los fenómenos bajo estudio. Además de explorar la anotación de los genes candidatos para encontrar términos enriquecidos, también es posible analizar el enriquecimiento de set de genes con características en común. Estos pueden ser set de genes asociados a un determinado pathway biológico, a una enfermedad, a una firma genética, etc. Existen diversas bases de datos de gene sets que pueden utilizarse para este tipo de análisis.

Gene Ontology

Los procesos biológicos considerados como “términos ontológicos” son aquellos definidos por el Consorcio Gene Ontology (GO). GO clasificó sistemáticamente las funciones de las proteínas de una célula en cuanto a 3 grandes criterios o categorías:

  • PROCESOS BIOLÓGICOS (BP)
  • FUNCIÓN MOLECULAR (MF)
  • COMPONENTE CELULAR (CC)

En cada una de estas categorías se definieron procesos, funciones y componentes celulares únicos de una manera sistemática, consistente y comparable de análisis a análisis, y mediante el uso de un vocabulario controlado.

Estos términos pueden ser utilizados a la hora de analizar el enriquecimiento de determinados procesos o funciones celulares en una lista de genes de interés. Esta lista podría ser, por ejemplo, los genes que mostraron diferencias significativas en su expresión entre dos condiciones de estudio (Ej.: Enfermos vs sanos, o muestras tratadas con una droga vs sin tratamiento).

KEGG y Reactome

Existen otras bases de datos curadas que contienen set de genes agrupados a partir de su participación en distintos pathways biológicos o vías metabólicas. Este es el caso de KEGG y Reactome. Examinando estas bases de datos es posible encontrar los pathways de los que forma parte un determinado gen o una lista de ellos. También cuentan con información de grupos de genes que participan de pathways que se encuentran desregulados en determinadas condiciones, como por ejemplo en ciertas enfermedades. De esta manera es posible utilizar estos set de genes para analizar el enriquecimiento de ciertos pathways en una determinada condición con respecto a un control.

Otras bases de datos de gene sets

Existen también bases de datos que recopilan todo tipo de gene sets, como el caso de MSigDB. Allí podemos encontrar una recopilación de los set de genes mencionados anteriormente (GO, KEGG, Reactome) y muchos otros más. Esta base de datos pretende centralizar colecciones de grupos de genes de distintas fuentes para su utilización en análisis de enriquecimiento.

Tipos de análisis de enriquecimiento

  • Over Representation Analysis (ORA): Es un método muy utilizado que permite determinar si la anotación de un proceso biológico está sobrerepresentada en una lista de genes de interés. Esta lista en general proviene de un análisis previo, como podría ser la lista de genes expresados diferencialmente obtenidas de un experimento de RNA-Seq. El p-valor para el enriquecimiento de genes analizado es calculado a través de una distribución hipergeométrica. Se comparar la presencia de genes de un determinados set en la lista de interés contra lo que se esperaría por una selección de genes al azar. Es decir que se debe tomar en cuenta el background o el total de genes del que se obtuvo esta lista. Un ejemplo de este tipo de análisis es el que realiza la aplicación web de DAVID https://david.ncifcrf.gov.

  • Gene Set Enrichment Analysis (GSEA): El análisis de enriquecimiento por sobre representación de set de genes es útil pero posee algunas limitaciones. Por empezar, al utilizar una lista previamente seleccionada de genes, es muy dependiente del cutoff de significancia elegido luego del análisis de expresión diferencial. Además se pierde información de genes que quizás solo aumentan o disminuyen en levemente pero de manera coordinada con otros genes pertenecientes a los mismos procesos biológicos. Para una mejor aproximación a este tipo de problemas se utilizan los análisis de enriquecimiento como el GSEA, que utilizan la lista completa de genes analizados evitando el bias inducido por la selección manual de un grupo de genes. Los genes son rankeados en base a su significancia estadística y fold change, o en base a la correlación que muestran con el fenotipo de estudio, dependiendo del método utilizado. Más allá de la estrategia de ranking elegida, lo que se termina ingresando como input al programa es una lista de todos los genes rankeados donde en un extremo tengo los genes que más aumentan en una condición del estudio y en el otro extremo los que más disminuyen (es decir, los que aumentan en la otra condición o en el control). Luego se calcula un score de enriquecimiento que representa el grado en el que un gene set está enriquecido en uno de los extremos de la lista. Dependiendo en qué extremo de la lista esté enriquecido podremos decir si este set de genes o pathway está aumentado o disminuido en una determinada condición. Luego se utiliza un test de permutación para asignar un p-valor y se ajusta por múltiples test con FDR. Como resultado entonces tendremos listas de gene sets que representan procesos biológicos aumentados de manera significativa en una de las dos condiciones del estudio.

  • Análisis de pathways basdos en redes topológicas : Este tipo de análisis incorpora al análisis de gene sets la topología de las redes de interacción de proteínas. Esto quiere decir que no solo va a tomar en cuenta los genes de un determinado pathway que se encuentran aumentados o disminuidos en una determinada condición, si no que también va a considerar como es la relación entre estos genes y en que posición de la vía se encuentra. Por ejemplo, podríamos tener un gen que se encuentra bien arriba (upstream) de un pathway y es clave en su regulación está aumentado en una determinada condición. Pero como este es el único gen del pathway que muestra una expresión upregulada un análisis de enriquecimiento como los vistos anteriormente no van a detectarlo como enriquecido. En cambio, un método que tome en cuenta la topología de las redes de interacción de las proteínas en los distintos procesos biológicos va a detectarlo como aumentado, ya que probablemente si ese gen es clave para inducir esa vía todo el proceso se encuetre desregulado. Si bien este tipo de métodos es el que mejor modela los patwhays biológicos, son bastante más complicados de realizar, con más costo computacional, y requieren de una base de datos con las interacciones entre proteínas de cada proceso biológico bien curada. Por este motivo al día de hoy se siguen utilizando en mayor proporción los métodos de sobre representación, y principalmente los de tipo GSEA.

Resumen del trabajo práctico: Análisis de enriquecimiento de gene sets en muestras de pacientes con COVID19

En este trabajo práctico utilizaremos los resultados del análisis transcriptómico obtenidos en el TP anterior para hacer análisis de enriquecimiento que nos permitan interpretar biológicamente las diferencias observadas entre las muestras de pulmón de los pacientes infectados con SARS-CoV-2 y los donantes sanos.

Para ello utilizaremos análisis de sobre representación de genes (ORA) y GSEA, utilizando paquetes de Bioconductor para R. Recordemos que para el análisis ORA necesitamos una lista de genes con expresión diferencial. Para ello tomaremos los genes con un p-ajustado (FDR) < 0.05. Para el caso de GSEA se utiliza toda la lista de genes rankeados. Los vamos a rankear en base al p-valor y al signo del logFC, de manera que queden en el extremo superior los genes más significativos aumentados en los pacientes enfermos y los más significativos down-regulados. En el medio de la lista quedan los genes cuya expresión no cambia entre las condiciones.

Para cada uno de estos análisis generaremos tablas y gráficos que nos permitan interpretar de mejor manera los resultados del experimento.

Desarrollo

Para realizar este trabajo práctico se deberá utilizar la máquina virtual

Entrar a la carpeta TP4-5 y hacer doble click en el archivo con extensión .Rproj. Esto abrirá de manera automática el proyecto guardado en el TP anterior.

Se puede continuar escribiendo en el script guardado del TP anterior o crear uno nuevo. Los objetos generados anteriormente se cargan automáticamente al abrir el proyecto, como pueden ver en el cuadro de arriba a la derecha.

En primer lugar vamos a realizar un análisis ORA con el paquete goseq. Para ello cargamos la librería correspondiente y generamos un vector de 0 y 1 a partir de los datos de la tabla top.table. Recuerden que este objeto tenía la tabla de todos los genes obtenidos del RNA-Seq con los valores estadísticos y el logFC entre ambos grupos. El nuevo vector tendrá un 0 si el gen correspondiente a esa posición en la lista no es significativo, o un 1 si lo es (con p-ajustado < 0.05). Luego la función nullp aplica un modelo que permite inferir la relación entre el largo de los genes y la probabilidad de que tenga diferencias significativas en su expresión. Esto es necesario ya que este paquete toma esto en cuenta a la hora de calcular el enriquecimiento de los sets de genes:

library(goseq)

# Fit the Probability Weighting Function (PWF)

isSigGene <- top.table$adj.P.Val < 0.05 & !is.na(top.table$adj.P.Val)
genes <- as.integer(isSigGene)
names(genes) <- top.table$Gene
pwf <- nullp(genes, "hg19", "geneSymbol")

Luego realizamos el análisis de enriquecimiento sobre los genes signficativos y generamos un plot para visualizar los resultados:

goResults <- goseq(pwf, "hg19","geneSymbol")

# Plot the top 10
library(dplyr)
library(ggplot2)

goResults %>% 
    top_n(50, wt=-over_represented_pvalue) %>% 
    mutate(hitsPerc=numDEInCat*100/numInCat) %>% 
    ggplot(aes(x=hitsPerc, 
               y=term, 
               colour=over_represented_pvalue, 
               size=numDEInCat)) +
        geom_point() +
        expand_limits(x=0) +
        labs(x="Hits (%)", y="GO term", colour="p value", size="Count")
  • Qué pueden observar de los resultados?
  • Cuál es la interpretación biológica de lo que están viendo?
  • Qué procesos biológicos están enriquecidos en la lista de genes expresados diferencialemente?
  • Con este análisis, se puede saber si los procesos están aumentados o diminuídos en los pacientes enfermos?

Ahora realizaremos el análisis de enriquecimiento con GSEA utilizando la implementación para R del paquete fgsea. En primer lugar, cargamos las librerías necesarias y generamos un ranking de todos los genes (no solo los significativos esta vez). Este ranking se basa en la significancia estadística (p-valor) y en la interpretación biológica (el signo del logFC). Además el p-valor se transforma con la función -log10 para que cuanto mayor sea la significancia estadística mayor sea el score para el ranking. Luego esto se divide por el signo del logFC (1 o -1) para que quede en un extremo u otro de la lista en base a si está aumentado o diminuido ese gen con respecto al control de pacientes sanos:

library(fgsea)
library(msigdbr)

# Ranked list
rank_list <- top.table[,c(1,3,6)]
rank_list$fcSign=sign(rank_list$logFC)
rank_list$logP=-log10(rank_list$P.Value)
rank_list$metric=rank_list$logP/rank_list$fcSign

ranks <- rank_list$metric
names(ranks) <- rank_list$Gene

En el siguiente paso vamos a descargar los gene sets que utilizaremos en el análisis desde la base de datos de MSigDB. Descargaremos los que corresponden a las bases de datos de Reactome, KEGG, GO, y los hallmarks (gene sets curados con los principales pathways)

# Gene sets
m_df = msigdbr(species = "Homo sapiens", category = "C2", subcategory = "CP:KEGG")
m_df2 = msigdbr(species = "Homo sapiens", category = "C2", subcategory = "CP:REACTOME")
m_df3 = msigdbr(species = "Homo sapiens", category = "C5")
m_df4 = msigdbr(species = "Homo sapiens", category = "h")
m_df_all <- rbind(m_df,m_df2,m_df3,m_df4)

m_list = m_df_all %>% split(x = .$gene_symbol, f = .$gs_name)

Ahora ya podemos correr el análisis de GSEA sobre la lista rankeada utilizando los gene sets descargados en el paso anterior. La función head nos muestra los primeros 10 gene sets ordenados de por el score (NES):

fgseaRes <- fgsea(m_list, ranks, minSize=15, maxSize = 500)
head(fgseaRes[order(padj, -abs(NES)), ], n=10)
  • Qué procesos biológicos aparecen más enriquecidos en este análisis?

Podemos visualizar alguno de los pathways enriquecidos con un enrichment score plot:

plotEnrichment(m_list[["GO_RESPONSE_TO_INTERFERON_GAMMA"]], ranks)

Este gráfico muestra el enriquecimiento del pathway (eje y) a lo largo de toda la lista de genes rankeada (eje x). Se observa claramente que hay un enriquecimiento de genes pertenecientes a este gene set en el extremo izquierdo del gráfico. Esto indica que la mayor parte de estos genes están aumentados significativamente en los pacientes con COVID19.

Podemos generar una tabla con los plots para los pathways más enriquecidos:

topUp <- fgseaRes %>% 
    filter(ES > 0) %>% 
    top_n(n=5, wt=-padj)
topDown <- fgseaRes %>% 
    filter(ES < 0) %>% 
    top_n(5, wt=-padj)
topPathways <- bind_rows(topUp, topDown) %>% 
    arrange(-ES)
plotGseaTable(m_list[topPathways$pathway], 
              ranks, 
              fgseaRes, 
              gseaParam = 0.5)

Agranden la figura con el botón “Zoom” arriba del cuadro donde les aparece el plot, en la esquina inferior derecha de la pantalla.

  • Qué observa? Qué puede decir acerca de los procesos que ve aumentados en los pacientes enfermos?

Por último utilizaremos otro paquete, el clusterProfiler, que permite hacer análisis de tipo ORA conectándose de manera directa con la base de datos KEGG. Este paquete contiene permite obtener buenas visualizaciones de los pathways que ayudan a interpretar los resultados.

Primero cargamos las librerías necesarias. Este paquete utiliza IDs de ENTREZ como identificadores de los genes en lugar de los nombres (gene symbols). Por eso utilizaremos el paquete biomaRt para traducir de uno a otro, solo para los genes significativos (recuerde que este paquete es para análisis de tipo ORA).

library(clusterProfiler)
library(biomaRt)

# To use this package is necessary to translate gene symbols to ENTREZ IDs

mart <- useMart(biomart = "ensembl", dataset = "hsapiens_gene_ensembl") # Database conexion.

sig_genes_entrez <-  getBM(attributes = c("hgnc_symbol", "entrezgene_id"), filters = "hgnc_symbol", values = sig_genes, bmHeader = T, mart = mart)
sig_genes_entrez <- sig_genes_entrez$`NCBI gene (formerly Entrezgene) ID`

En el siguiente paso corremos el análisis de enriquecimiento sobre los gene sets de la base de datos de pathways KEGG. Observamos los primeros 10 resultados con la función head:

kegg_paths <- enrichKEGG(gene = sig_genes_entrez, organism = 'hsa')
head(kegg_paths, n=10)
  • Qué puede decir acerca de estos resultados? Coincide con lo observado con los métodos anteriores?

Ahora pasaremos a visualizar los resultados obtenidos sobre los mapas de los procesos biológicos que posee KEGG. En primer lugar utilizamos la siguiente función para visualizar alguno de los procesos biológicos que aparecen enriquecidos en la tabla anterior indicándole el ID del mismo. Se abrirá una pestña nueva en el navegador de internet mostrando el pathway solicitado.

library(pathview)

# Choose an enriched pathway ID from the last table:

browseKEGG(kegg_paths, 'hsa04062') # Online at KEGG

Mediante el paquete pathview cargado recién es posible generar un archivo de imagen png que muestre un pathway de KEGG con los genes pintados en distinto color en base a su logFC. Los genes aumentados en los pacientes con COVID19 se verán en rojo y los disminuidos en verde.

En primer lugar tenemos que generar la tabla con los logFC de todos los genes para pasarle a la función que graficará el pathway. Para eso, primero debemos traducir los nombres de los genes IDs de ENTREZ como hicimos anteriormente, solo que esta vez lo hacemos para todos los genes.

mart <- useMart(biomart = "ensembl", dataset = "hsapiens_gene_ensembl")

genes_entrez <-  getBM(attributes = c("hgnc_symbol", "entrezgene_id"), filters = "hgnc_symbol", values = top.table$Gene, bmHeader = T, mart = mart)
top.table_entrez <- merge(top.table,genes_entrez,by.x="Gene",by.y="HGNC symbol", all.x=T)
colnames(top.table_entrez)[8] <- "Entrez"
logFC <- top.table_entrez$logFC
names(logFC) <- top.table_entrez$Entrez

Por último producimos las imágenes que se guardarán en la carpeta del TP. Para esto hay que pasarle a la función pathview la tabla con los logFC y el pathway que queremos graficar. Elegimos uno de los que salieron enriquecidos en el análisis e indicamos su ID en los parámetros de la función:

pathview(gene.data = logFC, 
         pathway.id = "hsa04062", 
         species = "hsa", 
         limit = list(gene=5, cpd=1))

Para ver la imagen puede hacer click directamente en el archivo png generado desde el RStudio en la pestaña Files en la esquina inferior derecha de la pantalla. También se puede navegar por la interfaz gráfica de Ubuntu hasta la carpeta del TP y abrir las imágenes desde ahí. La imagen que contiene los genes coloreados en base a su logFC es la que termina en pathview.png.

  • Qué información puede obtener de esta figura? Le ayuda a interpretar mejor lo que puede estar sucediendo con este proceso biológico?
  • Cómo lo interpreta?
  • Pruebe correr este último comando pero eligiendo otro pathway que también haya salido enriquecido. Qué puede decir de la expresión de los genes de este otro pathway.

Hemos llegado al final de esta serie de TPs en los que analizamos muestras de pacientes infectados con SARS-Cov-2 y las comparamos con tejido de pulmón de pacientes sin infección. Para eso recorrimos todo el pipeline de un típico análisis de RNA-Seq, desde los archivos fastq que obtenemos de la secuenciación hasta la interpretación biológica de los resultados observando los procesos y pathways desregulados en el tejido pulmonar de pacientes con COVID19. Excelente!

Bibliografía

Previous
Next

Last updated on Aug 24, 2023