Trabajo Práctico 4 - Análisis transcriptómico II

Objetivos:

El objetivo de este trabajo práctico es aprender a realizar un análisis de expresión diferencial entre grupos de muestras partiendo de una tabla de expresión en counts y utilizando paquetes bioinformáticos para el lenguaje R.

Introducción:

Uno de los primeros análisis que se puede realizar a partir de un experimento de transcriptómica es el análisis de expresión diferencial. Estos análisis nos permiten estudiar como se diferencia la expresión génica entre dos o más grupos de muestras. Si bien es posible examinar la expresión de uno o algunos pocos genes y compararla su expresión en distintos grupos mediante pruebas estádisticos clásicas, el objetivo de este tipo de estudios suele ser el análisis a gran escala de todos los genes presentes en el genoma para encontrar aquellos que están aumentados o disminuidos en una determinada condición.

Métodos y paquetes estadísticos

Un análisis de este tipo requiere de métodos estadísticos apropiados para este tipo de datos. Por un lado, se utilizan métodos que permiten modelar la distribución estadística de las counts. El modelado varía según el software a utilizar, y cada uno tiene sus ventajas y desventajas. Dentro de los programas más usados para estos se encuentran los paquetes para R: DeSeq2, EdgeR y el combo de los paquetes Voom y Limma. La principal caracteristica que diferencia a estos paquetes es la manera en la cuál estiman la varianza en la expresión de los genes. Para un mayor detalle del funcionamiento de cada uno encontrarán los papers en la bibliografía al final de esta guía.

Otra cosa que deben tener en cuenta este tipo de análisis, y que es común a todos los análisis realizados para datos ómicos, es el control del error de tipo I (falsos positivos). Esto se debe a que estamos realizando múltiples test al mismo tiempo sobre los datos (un test por cada gen, es decir que para un transcriptoma humano estamos realizando más de 20 mil tests a la vez). Si utilizamos un alfa de 0.05 para definir la significancia estadística, podemos esperar que un 5% de los tests nos arrojen resultados significativos cuando en realidad no hay diferencias. Al analizar 20,000 genes, vamos a esperar obtener 1,000 genes que muestren un p-valor menor a 0.05 por azar, pero que en realidad no tienen diferencias significativas entre los grupos. Para controlar la cantidad de falsos positivos se utilizan diversas metodologías. Una de las más utilizadas en los experimentos de genómica, transcriptómica y proteómcia es el ajuste por False Discovery Rate (FDR). A la hora de utilizar un cutoff para ver que genes muestran diferencias significativas se deben tener en cuenta los p-valores ajustados por FDR. El valor de cutoff utilizado se relaciona con el porcentaje de falsos positivos que estemos dispuestos a aceptar. Para el desarrollo de este TP utilizaremos un FDR < 0.05 para definir la significancia estadística.

Fold change

Los resultados del análisis de expresión diferencial se suelen mostrar como Fold changes (FC) o Log2(FC). Esto se puede interpretar como la magnitud en la que dififiere la expresión de un gen entre dos grupos. Por ejemplo, si la media del gen A en una muestra es 10 y en el control es 5, el FC será 10/5 = 2. Para una mejor interpretación y visualización se suele aplicar el logaritmo en base 2 al FC. De esta manera, en el ejemplo anterior el Log2(FC) quedaría en 1, indicando un aumento del doble en la expresión. Si fuera al revés, 5/10, el Log2(FC) sería -1, lo que indica que el primer grupo tiene la mitad de la expresión que el control. Un Log2 (FC) = 0 indica que no hay cambio.

Datos públicos y tipos de normalización

Tal como se mencionó en el TP anterior, es posible utilizar bases de datos públicas con datos transcriptómicos de interés. En la mayor parte de los casos se podrán encontrar las tablas de counts, es decir que nos evitamos hacer el pipeline del TP anterior. Esto tiene ventajas y desventajas, ya que si bien estamos limitados por el análisis que hizo alguien más, podemos acceder a una fuente de información inmensa a un mínimo costo económico y computacional. Es importante destacar que los datos preprocesados pueden venir o no normalizados. Si lo que tenemos es la tabla de counts (lo que obtuvimos en el TP3) decimos que tenemos la tabla crudao sin normalizar. Esto es lo ideal ya que los programas de análisis de expresión diferencial parten de las counts porque utilizan sus propios métodos de normalización. Además en el caso de necesitar algún tipo de normalización particular podemos hacerlo fácilmente desde estos datos. Otra posibilidad es que los datos vengan en FPKM o TPM, en este caso decimos que tenemos una tabla normalizada. El problema de esto es que estos datos no se pueden utilizar con la mayoría de los paquetes de análisis de expresión diferencial. Solo son útiles para algunos análisis particulares y para visualización. Si no podemos acceder a las counts y solo tenemos estos datos normalizados se deberá recurrir un pipeline alternativo para el análisis de expresión diferencial, aunque los resultados no serán del todo confiables.

Bioconductor y R

Para este TP y el siguiente utilizaremos RStudio, un entorno de programación para el lenguaje estadístico R. Este lenguaje es ampliamente utilizado en bioinformática para el análisis de datos biológicos. Su uso extendido se debe en gran parte a la enorme cantidad de paquetes bioinformáticos que se encuentran disponibles para R. La mayor parte de estos paquetes se encuentran dentro de un proyecto denominado Bioconductor ( https://bioconductor.org). En la web pueden encontrar mucha información sobre este proyecto, los paquetes disponibles y como instalarlos en R. Para este trabajo práctico todos los paquetes necesarios ya se encuentran instalados en la máquina virtual.

Resumen del trabajo práctico: Análisis transcriptómico de pacientes con COVID19 II

Para este trabajo continuaremos el pipeline de análisis de RNA-Seq en donde lo dejamos la última vez. En lugar de utilizar la tabla de counts obtenida anteriormente (que solo tenía la expresión de los genes del cromosoma 21) usaremos la tabla real descargada de GEO (GSE147507), que contiene la expresión para todos los genes del genoma, tanto para las dos muestras de pacientes con COVID19 como para los controles.

Desarrollo

Para este TP se debe utilizar la máquina virtual provista por los docentes

Iniciar RStudio desde el ícono ubicado en el dock de Ubuntu, a la izquierda de la pantalla.

Crear un nuevo proyecto haciendo click en File -> New Project -> Existing Directory. Hacer click en Browse y seleccionar la carpeta TP4-5 que se encuentra dentro de la carpeta TPs. Hacr click en Open y luego en Create Project.

Ahora vamos a crear un archivo nuevo en donde escribiremos el código. Para eso hacer click en el símbolo de una hoja en blanco con un “+” arriba a la izquierda, y hacer click en R Script.

Guardar el archivo con el nombre que deseen. El mismo se guardará en la carpeta de trabajo (TP4-5).

Para correr los comandos de la guía deben escribirlos en la ventana del script que crearon recién (esquina superior izquierda), y luego con el cursor parado en la misma línea del comando hacer ctrl + Enter. Esto hará que la línea de código se corra en la terminal y mostrará allí el output. En el caso de que la salida sea un gráfico el mismo se verá en la ventana de la esquina inferior derecha de la pantalla.

R se basa en el uso de objetos. Estos objetos se pueden ver como variables que pueden almacenar todo tipo información, desde números y cadenas de texto hasta tablas u objetos más complejos. Para poder utilizar un objeto se escribe el nombre que se utilizó al crearse. Por ejemplo, si quisiear almacenar el número 6 en un objeto que se llame “nuevo” lo haría de la siguiente manera:

nuevo <- 6

También podemos guardar la salida de una función. Por ejemplo usemos la función mean que calcula la media de una serie de números y guardamos el resultado en una variable “x”:

x <- mean(c(4,6,8,10))

Observen que para pasarle una serie de números (o vector), le pasamos los números separados con comas dentro de c().

Existen muchísimos cursos de R gratis muy buenos en la web. Es muy recomendable que aprendan aunque sea un manjeo básico de este lenguaje ya que brinda la posibilidad de manejar datos de grandes dimensiones con relativa facilidad. Debido a que los experimentos generan cada vez mayor cantidad de datos, su uso resulta indispensable para su manejo y posterior análisis.

Primero cargaremos las librerías necesarias para trabajar (es decir, los paquetes con las funciones que utilizaremos). Luego cargaremos la tabla de counts en R. Recuerde que si hay más de una línea de comandos debe hacer ctrl + Enter en cada línea luego de escribirla en la ventana del script, o puede pintar todas las líneas y hacer ctrl + enter para que se ejecuten automáticamente una después de la otra.

library(edgeR) # As edgeR depends on "limma", both packages are loaded
library(dplyr)
library(ComplexHeatmap)
library(ggpubr)
library(Glimma)
library(viridis)
library(circlize)


data <- read.table("./covid19_human_counts.txt", header = T, sep = "\t", row.names = 1)
clinical_data <- data.frame(row.names = colnames(data), group = c("healthy","healthy","covid19","covid19"))
group <- as.factor(clinical_data$group)

Luego creamos el objeto para ser procesado por los paquetes edgeR, Voom y Limma que utilizaremos.

data_dge <- DGEList(data)
data_dge <- calcNormFactors(data_dge) # Normalization factor estimation
data_dge$genes <- data.frame(rownames(data_dge))

Filtrado de genes con muy baja expresión para todas las muestras.

# filter genes with low expression
cutoff <- 1
drop <- which(apply(cpm(data_dge), 1, max) < cutoff)
data_filtered <- data_dge[-drop,] 
dim(data_filtered) # number of genes left

lcpm <- cpm(data_filtered, log=TRUE) # Log 2 of Counts per Million (CPM)

Plot de MDS. Este gráfico muestra la distancia entre las muestras reduciendo la dimensionalidad a solo dos dimensiones (originalmente tenemos tantas dimensiones como genes hay en el dataset).

plotMDS(lcpm, col = as.numeric(group)) # Multidimensional scaling (MDS) plot
  • Qué observa? Cómo se agrupan las muestras? (Es posible agrandar la imagen haciendo click en “Zoom”)

Ahora utilizaremos una función que grafica un MDS plot pero de manera interactiva. Al ejecutarlo se les abrirá un navegador web.

glMDSPlot(lcpm, label= colnames(lcpm) ,groups = group)

Ahora generaremos el modelo lineal que utilizará voom + limma para el análisis de expresión diferencial:

model_matrix <- model.matrix(~0 + group)
data_voom <- voomWithQualityWeights(data_filtered, model_matrix, plot = T)

Esto generará gráficos asociados con el modelo. Pueden observar la relación entre el nivel de expresión en counts y el desvío estándard en las muestras. Se puede ver que cuanto menor expresión tenga un gen mayor es su desvío estándard. Esto es lo que intentamos corregir con voom. Una vez calculada la corrección podemos utilizar limma para el análisis de expresión diferencial propiamente dicho.

fit <- lmFit(data_voom, model_matrix)

Con esto ajustamos el modelo lineal a nuestros datos.

Ahora procedemos a indicar los contrastes o comparaciones que queremos hacer. El contraste se indica como una diferencia entre dos grupos. Esto es importante a la hora de interpretar los resultados, ya que al poner primero al grupo infectado, cuando veamos un log2FC positivo indicará un aumento en la expresión del gen en este grupo, y cuando sea negativo indicará que el grupo control tiene mayor expresión.

contr <- makeContrasts(groupcovid19 - grouphealthy, levels = colnames(coef(fit)))
contr

Ajustamos el modelo en base a los contrastes indicados en el paso anterior. Los resultados del análisis de expresión diferencial se guardan en el objeto top.table. Con el comando head podemos ver los primeros 20 genes ordenados por p-valor.

tmp <- contrasts.fit(fit, contr)
tmp <- eBayes(tmp)

top.table <- topTable(tmp, sort.by = "P", n = Inf)
head(top.table, 20)
  • Cuáles son los genes que muestran diferencias más significativas entre los grupos? Están aumentados o disminuídos en los pacientes infectados?

Con el siguiente comando podemos ver la cantidad de genes expresados diferencialmente:

dt <-  decideTests(tmp)
summary(dt) # Number of DE genes
  • Cuántos genes expresados diferencialmente encontramos? Cuántos están up-regulados en los pacientes y cuantos en los controles?

Guardamos los resultados en un archivo de texto:

top.table$Gene <- rownames(top.table)
top.table <- top.table[,c("Gene", names(top.table)[1:6])]
write.table(top.table, file = "Covid19vsHealthy.txt", row.names = F, sep = "\t", quote = F)

Podemos generar un plot interactivo que se abrirá en el navegador. Allí podemos visualizar los genes estadísitcamente significativos y buscar cualquier otro gen de interés. También nos muestra un gráfico comparando la expresión de las muestras de los distintos grupos.

glMDPlot(tmp, status=dt, main=colnames(tmp)[1],
         side.main="gene", counts=lcpm, groups=group, launch=T)
  • Pruebe buscar el gen CCL2. Está aumentado o disminuído? Es estadísticamente significativo?

Ahora vamos a generar un heatmap que muestre la expresión de los genes significativos (filas) y las muestras (columnas). Para ello filtramos la tabla de LogCPM con lo genes que salieron diferencialmente expresados. Debemos estandarizar los valores por medio del Z-Score (lo recuerdan de Bioestadística?), esto mejora la visualización de todos los genes juntos.

sig_genes <- filter(top.table, adj.P.Val < 0.05) %>% select(Gene) %>% unlist()
lcpm_sig <- t(scale(t(lcpm[sig_genes,])))


color_hm <- colorRamp2(seq(-2,2,by=0.02),inferno(length(seq(-2,2,by=0.02)))) # Try inferno instead of viridis!

annotation <-HeatmapAnnotation(Group=clinical_data[colnames(lcpm_sig),],
                                           annotation_name_gp = gpar(fontface="bold"),gap = unit(1,"mm"),
                                       border = TRUE,which="column",
                                           col=list(Group=c("healthy"="blue","covid19"="red")))

heatmap<-Heatmap(lcpm_sig,name="Z-score",col=color_hm,
                   border="black",show_column_names = TRUE,show_column_dend = TRUE,show_row_names=TRUE, show_row_dend = TRUE,
                   row_names_gp = gpar(fontsize=3),column_names_gp =gpar(fontsize=5),column_names_rot = 45, 
                 show_parent_dend_line = FALSE,row_gap = unit(0,"mm"),show_heatmap_legend = TRUE,
                   cluster_rows = TRUE,cluster_columns = TRUE,column_gap = unit(1,"mm"),column_title_gp = gpar(fontsize=5),
                   top_annotation = annotation,
                   clustering_method_rows = "ward.D2",clustering_method_columns = "ward.D2")
draw(heatmap, merge_legend = TRUE)
  • Qué observa en el gráfico? Recuerde hacer click en Zoom para verlo mejor.

También podemos generar un plot que muestre la expresión de algún gen en particular. Probemos con el gen que buscamos antes, CCL2:

lcpm_ccl2 <- merge(t(lcpm),clinical_data,by="row.names")
lcpm_ccl2 <- lcpm_ccl2[,c("CCL2","group")]

p <- ggbarplot(lcpm_ccl2, x = "group", y = "CCL2", title = "CCL2",
               color = "group", line.color = "gray", line.size = 1, xlab = F,ylab = "CCL2 mRNA expression (Log2CPM)",
               palette = "aaas",short.panel.labs = T,combine=F,
               add = c("mean_sd","jitter"))+theme(text = element_text(size = 20),plot.title = element_text(hjust = 0.5))  +
  stat_compare_means(comparisons = c("healthy", "covid19"), label = "p.signif", method = "wilcox.test", paired = FALSE,size=8)
ggpar(p, legend="none")
  • Ve lo mismo que lo que vio en el plot interactivo generado anteriormente?

Guardar todo al cerrar RStudio ya que utilizaremos los resultados obtenidos en el próximo TP

Con esto finalizamos la segunda parte del trabajo prácatico de transcriptómica. Partimos de una tabla de counts hasta llegar a una lista de genes expresados diferencialmente entre infectados y donantes sanos. Además generamos un archivo con el LogFC entre los grupos junto con la significancia estadística del análisis (p-valor y FDR) para cada uno de los genes del genoma.

Ahora que tenemos toda esta información necesitamos darle una interpretación biológica. Para eso en el próximo TP realizaremos análisis de enriquecimiento de pathways o vías biológicas y de Gene Ontology!

Bibliografía

Previous
Next

Last updated on Aug 24, 2023