Trabajo de Fin de Curso R Vicente Fernández Rodero

Script para el proyecto de fin del Curso de Programación y análisis de datos con R de la USAL (Abril de 2025) - Vicente Fernández Rodero.

Este Script pretende analizar datos de RNAseq de Neuroblastoma (GSE224037) https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE224037

Los datos que utilizo son datos de RNAseq de acceso público a través de NCBI. Los datos utilizados tienen datos de secuenciación de dos líenas celulares de Neuroblastoma (BE2C y Kelly) en las que se ha silenciado AhR (Aryl Hydrocarbon Receptor) a través del uso de dos shRNAs. A efectos de este script utilizo solo para el análisis los datos de la línea celular BE2C. En un caso real sería conveniente analizar ambas líneas y sacar conclusiones con los datos de ambas (Como checkeo extra incluyo también una copia del script para ver si el mismo análisis es reproducible para la línea celular Kelly ( en carpeta “/Kelly cells reproducible”).

**NOTA: En este archivo alguna de las fuciones están modificadas para que en lugar de guardar las gráficas en PDF haga print de ellas en el documento. Además para los análsis GO y GSEA también se ha quitado la parte de las funciones que guardan los datos en un archivo excel.

1. Directorio de trabajo y carga de paquetes a utilizar

Lo primero que hago es definir el directorio en el que voy a trabajar:

setwd(dir = "~/Desktop/R/Trabajo Fin Curso R/")

Luego cargo los paquetes que utilizaré para analizar los datos seleccionados:

# Paquetes generales:
library(tidyverse)
library(furrr)
library(openxlsx)
library(ggrepel)
# Para importar datos de NCBI:
library(GEOquery)
# Para expresión diferencial:
library(DESeq2)
# Dispersión muestras
library(FactoMineR)
library(factoextra)
#Venn Diagrams
library(ggvenn)
# GO & GSEA 
library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)
library(AnnotationDbi)
library(ggridges)
# Heatmap
library(ComplexHeatmap)
library(colorRamp2)

2. Importar datos de NCBI de nuestro problema (“GSE224037”)

Como los datos a utilizar se pueden acceder desde NCBI, hago uso del paquete GEOquery que permite acceder a los datos directamente desde R. Además, modifico y selecciono los datos de interés para que sea más sencillo utilizarlos más adelante.

GEOquery::getGEOSuppFiles(GEO ="GSE224037", fetch_files = TRUE)
  # Obtener metadata y seleccionar los datos que nos interesan
metadata <- GEOquery::getGEO(GEO ="GSE224037")
  # Trato los datos para quedarme la info  y formato/nombre de interés 
sampledata <- metadata$GSE224037_series_matrix.txt.gz@phenoData@data %>% 
  dplyr::select(title, geo_accession, `cell line:ch1`, `genotype:ch1`, `treatment:ch1`) %>% 
  dplyr::rename(Cell =`cell line:ch1`, Genotype =`genotype:ch1`, Treatment=`treatment:ch1` ) %>% 
  mutate(title = gsub("-", ".", title)) %>% # *NOTA
  mutate(Treatment = gsub("AhR-", "", Treatment)) 
rownames(sampledata) <- NULL
    # *NOTA: en los counts hay separación por puntos y en la info por guión (sustituyo para que el nombre sea igual)

Separo los datos de ambos shRNAs y los modifico para que tengan el formato recomendado para su análisis con DESeq2

  # Datos de la línea BE2C celular separando ambos shRNAs 
BE2C.gz <- gzfile("GSE224037/GSE224037_BE2C_RawCounts.csv.gz")
BE2C.csv <- as.data.frame(read.csv(BE2C.gz) %>% dplyr::rename(Symbol=X) %>% column_to_rownames(var="Symbol"))
      #Para el shRNA 1
BE2C.sh1 <- BE2C.csv %>%  dplyr::select(-BE2C.shAhR2.1, -BE2C.shAhR2.2, -BE2C.shAhR2.3)
BE2C.sh1.info <- sampledata %>% subset(Cell == "BE2C" & Treatment != "sh2" )
BE2C.sh1.info$Treatment <- as.factor(BE2C.sh1.info$Treatment) # DESeq recomienda que el diseño sea factor y no character
ncol(BE2C.sh1) == nrow(BE2C.sh1.info) # Checkear mismo nº de filas para crear DESeq2 object

      #Para el shRNA 2
BE2C.sh2 <- BE2C.csv %>%  dplyr::select(-BE2C.shAhR1.1, -BE2C.shAhR1.2, -BE2C.shAhR1.3)
BE2C.sh2.info <- sampledata %>% subset(Cell == "BE2C" & Treatment != "sh1" )
BE2C.sh2.info$Treatment <- as.factor(BE2C.sh2.info$Treatment) # DESeq recomienda que el diseño sea factor y no character
ncol(BE2C.sh2) == nrow(BE2C.sh2.info) # Checkear mismo nº de filas  = nº columnas para crear DESeq2 object

Utilizo DESeq2 para crear para ambos shRNAs un objeto de DESeq2 que permite analizar los datos de RNAseq definiendo los grupos a comparar y cual utilizar como control de la comparación de los datos de entrada (counts).

  # Crear DESeqDataSet de cada una de las líneas celulares 
BE2C_sh1.dds <- DESeq2::DESeqDataSetFromMatrix(countData = BE2C.sh1,
                                           colData = BE2C.sh1.info,
                                           design = ~ Treatment)

BE2C_sh2.dds <- DESeq2::DESeqDataSetFromMatrix(countData = BE2C.sh2,
                                               colData = BE2C.sh2.info,
                                               design = ~ Treatment)

Uno ambos objetos DESeq2 creados en una lista para después modificarlos y procesarlos a la vez aplicando funciones.

  # Unir los datos de ambos sh en una lista para modificaciones con purrr y furrr
BE2C.counts <- list(BE2C_sh1.dds, BE2C_sh2.dds)
names(BE2C.counts) <- list("BE2C_sh1.dds", "BE2C_sh2.dds")

3. Realizar pre-procesamiento

Primero creo una función que me permita filtrar y no tener en cuenta genes que no tengan relevancia biológica para este problema, para ello selecciono como corte (no muy estricto) que en total de las tres replicas del control + las tres replicas del shRNA haya al menos un total de 10 counts. DEfino también que la referencia de las comparaciones sea el control (pLKO: plásmido sin la secuencia del shRNA). Por último aplico a la lista la función creada de manera que se modifiquen las comparaciones de ambos shRNAs de un solo paso.

  # Filtrar filas (genes) con al menos 10 cuentas en total y determinar el 
  #factor/caracteristica de diseño (control) de referencia para comparar
filter.counts <- function(data){
  a <- rowSums(counts(data)) >= 10
  data <- data[a,]
  data$Treatment <- as.factor(data$Treatment)
  data$Treatment <- relevel(data$Treatment, ref = "pLKO")
return(data)
}
BE2C.counts <- purrr::map(.x=BE2C.counts, .f= filter.counts)

4. Realizar expresión diferencial con DESeq2

Una vez pre-procesados los datos contenidos en los objetos DESeq2 realizo el análisis de expresión diferencial. Para ello creo dos tipos de funciones de manera que pueda tener los datos en crudo (de cada uno de los shRNAs) y los datos filtrados de solo los genes que se expresan diferencialmente de forma significativa.

  # Expresión diferencial sin filtrar
DESeq2.res <- function(data, alpha =0.05){
  data <- DESeq(data)
  result <- DESeq2::results(data, alpha = alpha)
  result <- as.data.frame(result) %>% rownames_to_column(var="Symbol") %>% 
    mutate(negLog10Pvalue = -log10(pvalue),
           rank = rank(stat)) %>% 
    dplyr::arrange(desc(rank)) %>% 
  return(result)
}

 # Expresión diferencial filtrando significativas p<0.05 y -0.5 > FC > 0.5
Diff.exp <- function(data, alpha =0.05){
  data <- DESeq(data)
  result <- DESeq2::results(data, alpha = alpha)
  result <- as.data.frame(result) %>% rownames_to_column(var="Symbol") %>% 
    subset(padj<alpha) %>% 
    subset( log2FoldChange < -1 | log2FoldChange > 1 )
  return(result)
}
BE2C.res <- purrr::map(.x=BE2C.counts, .f= Diff.exp, .progress = TRUE)
names(BE2C.res) <- list("BE2C_sh1.res", "BE2C_sh2.res")

list2env(BE2C.res, envir = .GlobalEnv)

5. Explorar genes diferencialmente expresados (DE) en ambos shRNA

Una vez ya tengo los datos de expresión diferencial genero diferentes formas de visualizarlos y sacar diferentes conclusiones de ellos.

Lo primero que hago es ver como se distribuyen mis muestras, es decir como se distribuyen las replicas del control, sh1 y sh2 en función de dos componentes que deberían diferenciarlos.

 # Realizar diferentes representaciones de los datos de Diff exp.
    # Dispersión de las muestras
  log_TPM <- log2(BE2C.csv + 1)  # Avoid log(0) issues
  pca_res <- prcomp(t(log_TPM))  # Transpose to have samples as rows
  print(fviz_eig(pca_res)) # To see variances and choose for presenting

  print(fviz_pca_ind(pca_res,
                   geom = "text",
                   col.ind = "black", 
                   repel = TRUE) )

Como se puede observar al comparar las muestras vemos que para cada uno de los tipos (control, sh1 y sh2) sus replicas se agrupan entre sí pero no entre ellas, lo que nos puede indicar que al sileciar AhR las células se diferencian del control. Sin embargo, en base a estos dos componentes vemos que los sh1 y sh2 no se agrupan lo que también nos puede reflejar que quizás la eficiencia de silenciamiento también varía entre ambos (cosa que se puede apreciar en el paper de origen).

Por otro lado, exploro los genes diferencialmente a través de Volcano Plots y haciendo un ranking de los genes en función del estadístico stat que calcula DESeq2 (Wald). Para ello como utilizo en ambas representaciones los datos sin filtrar creo una función que me represente ambas gráficas con ggplot.

    # Volcano y Gene Score plots de los datos de expresión 
      # Utilizo los resultados sin filtrar de DESeq2 
DESeq2.DE <- purrr::map(.x=BE2C.counts, .f= DESeq2.res, .progress = TRUE)
names(DESeq2.DE) <- list("SH1_DE", "SH2_DE")
      # Creo una función para representar de un solo paso todo (lo llamaré VS.plot)
  VS.plot <- function(data, pthreshold = 1.3 , xlim= 6, ylim =20){
    print(ggplot(data, aes(x = log2FoldChange, y = negLog10Pvalue, color = Symbol)) +
            geom_point(color = "grey75" ) +
            geom_point(data = subset(data, log2FoldChange > 1 &  negLog10Pvalue > pthreshold), color = "red") +
            geom_point(data = subset(data, log2FoldChange < -1 &  negLog10Pvalue > pthreshold), color = "blue") +
            geom_text_repel(data = dplyr::slice_max(data, rank, n = 20),
                            aes(label = Symbol),
                            color = "black", size = 3, box.padding = 0.5,
                            point.padding = 0.3, max.overlaps = Inf)  +
            geom_text_repel(data = dplyr::slice_min(data, rank, n = 20),
                            aes(label = Symbol),
                            color = "black", size = 3, box.padding = 0.5,
                            point.padding = 0.3, max.overlaps = Inf)  +
      theme_classic() +
      labs(title = "Volcano Plot", x = "Log2 Fold Change", y = "-log10(p-value)") +
      geom_hline(yintercept = 1.3, linetype = "dashed", color = "black", size = 0.5) +  
      geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "black", size = 0.5) + 
      xlim(c(-xlim, xlim)) + ylim(c(0, ylim)) ) 
    
    print(ggplot(data, aes(x = rank, y = stat, color = Symbol)) +
            geom_point(color = "grey75" ) +
            geom_point(data = subset(data, log2FoldChange > 1 ), color = "red") +
            geom_point(data = subset(data, log2FoldChange < -1 ), color = "blue") +
            geom_text_repel(data = dplyr::slice_max(data, rank, n = 20),
                      aes(label = Symbol),
                      color = "black", size = 3, box.padding = 0.5,
                      point.padding = 0.3, max.overlaps = Inf)  +
            geom_text_repel(data = dplyr::slice_min(data, rank, n = 20),
                            aes(label = Symbol),
                            color = "black", size = 3, box.padding = 0.5,
                            point.padding = 0.3, max.overlaps = Inf)  +
      theme_classic() +
      labs(title = "Gene Score Distribution", x = "Rank", y = "Stat") )
  }

Para el sh1:

VS.plot(DESeq2.DE$SH1_DE)

Para el sh2:

VS.plot(DESeq2.DE$SH2_DE)

Por otro lado, haciendo uso de los datos filtrados, represento en un driagrama de Venn los genes diferencialmente expresados de forma significativa para ambos shRNAs, convergiendo en el centro los que se expresan diferencialmente en ambos shRNAs al compararlos con el control.

    # Venn Diagram para ver el número de genes DE comunes
Venn <-  BE2C.csv %>% rownames_to_column(var = "Symbol") %>% 
  dplyr::mutate(sh1.DE= ifelse(Symbol %in% BE2C_sh1.res$Symbol, yes=TRUE, no= FALSE),
                sh2.DE= ifelse(Symbol %in% BE2C_sh2.res$Symbol, yes=TRUE, no= FALSE)) %>% 
  dplyr::select(Symbol, sh1.DE, sh2.DE)

print(ggvenn(Venn, columns =  c("sh1.DE", "sh2.DE"), show_outside = "none", auto_scale = FALSE) ) # A efectos de visualización pongo autoescalado FALSE (en el PDF adjunto está como TRUE para ver también la diferencia en genes DE en el tamaño)

Como podemos ver en los la representación del Volcano Plot, Gene Score y el Venn Diagram, el sh2 tiene como resultado un mayor número de genes DE, lo que puede explicar también las diferencias en la distribución de las muestras comentadas anteriormente. En total 474 genes en total están DE.

6. Ver GeneOntology terms enriquecidos en nuestro set de genes

Utilizando el set de genes que de forma común se encuentran DE en ambos shRNAs, realizo GeneOntology para ver a que términos están asociados estos genes. Hago el análisis para el total de los 474 genes, de igual manera también hago este análisis para los genes que dentro de estos 474 suben o bajan en ambos shRNAs. Para ello, creo una función que me permita realizar el análisis un un menor número de pasos en las tres aproximaciones.

  # Crear una lista de genes comunmente DE entre ambos shRNA

BE2C.DE.list <- inner_join(x = BE2C_sh1.res, y = BE2C_sh2.res, by= "Symbol") 

colnames(BE2C.DE.list) <- gsub(pattern = ".x", replacement = "_sh1", colnames(BE2C.DE.list), fixed = TRUE)
colnames(BE2C.DE.list) <- gsub(pattern = ".y", replacement = "_sh2", colnames(BE2C.DE.list), fixed = TRUE)
  
  # Seleccionar Lista de genes que comunmente en los sh están up- or down- regulated
BE2C.DE.up <- BE2C.DE.list %>% subset(BE2C.DE.list$log2FoldChange_sh1 > 1 & BE2C.DE.list$log2FoldChange_sh2 > 1 )

BE2C.DE.down <- BE2C.DE.list %>% subset(BE2C.DE.list$log2FoldChange_sh1 < -1 & BE2C.DE.list$log2FoldChange_sh2 < -1 )

  # Realizar GO Analysis de los genes que comunmente en los sh están up- or down- regulated

    # PerformGO: Realiza GO enrichment
    PerformGO <- function(data, universe) {
    data <- enrichGO( 
    gene = data,
    OrgDb = org.Hs.eg.db,
    keyType = "SYMBOL",
    ont = "ALL",
    universe = universe,
    pAdjustMethod = "BH",
    pvalueCutoff = 0.05
  )
  data <- dplyr::arrange(data, desc(RichFactor))
  return(data)
  }

    # GO.visualization: Guarda gráficas como PDF
    GO.visualization <- function(data, PDF_name) {
  enrich_result <- pairwise_termsim(data)

  print(dotplot(data, showCategory = 20, font.size = 7.5, title = "TOP GO terms", x = "GeneRatio"))
  print(barplot(data, showCategory = 20, font.size = 6, title = " TOP GO terms"))
  print(emapplot(enrich_result, showCategory = 20))
    }

    # GO.analysis: Combina ambas funciones
    GO.analysis <- function(data, data_name,  universe) {
  # Realizar GO analysis
  GO.res <- PerformGO(data$Symbol,  universe = rownames(BE2C.sh1))
  # Generar y guardar las gráficas
  GO.visualization(GO.res, PDF_name = data_name)
    }

Para el total de los 474 genes DE:

  GO.analysis(data = BE2C.DE.list, data_name = "All DE genes GO")

Para los genes que comunmente suben su expresión:

  GO.analysis(data = BE2C.DE.up, data_name = "UP DE genes GO")

Para los genes que comunmente bajan su expresión:

  GO.analysis(data = BE2C.DE.down, data_name = "DOWN DE genes GO")

Como podemos ver de los tres análisis, tanto en el análisis del total de 474 genes como en los genes que bajan en expresión vemos términos relacionados con el ciclo y división celular. Lo que tendremos en cuenta para ejemplificar la aproximación de como proceder en la última parte de este documento.

7. Analisis GSEA de nuestro set de genes DE

De forma similar, realizo un análisis GSEA de todos los genes (utilizando de nuevo los datos sin filtrar y teniendo en cuenta el ranking en base a stat). Para ello creo de nuevo una función para aplicar a ambos shRNAs.

  DESeq2.gsea <- purrr::map(.x=BE2C.counts, .f= DESeq2.res, .progress = TRUE)
  names(DESeq2.gsea) <- list("SH1_GSEA", "SH2_GSEA")

  GSEA.analysis <- function(data){
      GSEA.res <- gseGO(
        geneList = setNames(data$rank, data$Symbol), 
        OrgDb = org.Hs.eg.db, 
        keyType = "SYMBOL", 
        ont = "ALL", 
        pAdjustMethod = "BH", 
        pvalueCutoff = 0.05 )
      
      gseaplot <- lapply(1:10, function(i) {
        gseaplot2(GSEA.res, 
                  geneSetID = i, 
                  title = GSEA.res@result$Description[i], 
                  subplots = 1:2)
     })
      
       # Generar y guardar las gráficas
     print(gseaplot)
     print(dotplot(GSEA.res, showCategory = 20, title = "Top GSEA Enriched GO Terms", font.size = 8))
     print(ridgeplot(GSEA.res) + ggtitle("Distribution of Enrichment Scores") +theme(axis.text.y = element_text(size = 8) ))
     print(emapplot(pairwise_termsim(GSEA.res), showCategory = 20, label_format = 20))
       
  }

Para el sh1:

GSEA.analysis(DESeq2.gsea$`SH1_GSEA`)    

Para el sh2:

GSEA.analysis(DESeq2.gsea$`SH2_GSEA`)    

8. Conclusiones y evaluación del problema

Finalmente, para ejemplificar como finalizaría el análisis de estos datos, analizaría alguna de las vías/terminos de interés y/o que saliera enriquecida en los análisis anteriores. Para ejemplificar utilizo un set de genes relacionados con Ciclo celular y División para ver como se expresan en un heatmap en el que represento log2FC de cada shRNA respecto al control (alternativamente con los counts podría representar TPM o calcular y represental z-Scores de las muestras). De esta forma de ese set de genes de traget (de ciclo/división) represento los que están comunmente DE en ambos shRNAs (reduciendo su número). Veo por tanto que acorde con lo que se veía en el GO análisis vemos que al silenciar AhR hay reguladores clave del ciclo y división celular alterados, lo que podría usarse en este ejemplo como base para una validación funcional de esta posible regulación de AhR.

 # Crear una lista de targets de AhR y anotación de su función principal

Target.list <- data.frame(Targets = c("CDK1", "CDK2", "CDK4", "CDK6", "CCNA2", "CCNB1",
                                                           "CCND1", "CCNE1", "RB1", "TP53", "E2F1", "MYC", 
                                                           "AURKA", "AURKB", "PLK1", "BUB1", "MAD2L1", "CDC20",
                                                           "CDH1", "SKP2", "APC", "PTTG1", "TTK", "CENPE", 
                                                           "KIF11", "KIF2C", "RACGAP1", "CYLD", "CHEK1", 
                                                           "CHEK2", "PCNA", "RFC1", "CDC6", "MCM2", "MCM7",
                                                           "ORC1", "CCNH", "CDK7", "CDK9", "KIF23", "NEK2",
                                                           "ESPL1", "SGO1", "SGO2", "INCENP", "CDK11B", "CCNH",
                                                           "CDK12", "CDK13", "FOXM1", "UBE2S", "CDK5", "CDK10",
                                                           "CCND2", "CCND3", "CCNA1", "CCNB2", "CCNB3", "CCNE2",
                                                           "RBL1", "RBL2")) %>% 
  dplyr::mutate(Function = "Cell cicle")


  # Seleccionar de ambos shRNA los log2FC respecto al control para un HeatMap
BE2C.HM <- BE2C.DE.list %>% dplyr::select(Symbol, log2FoldChange_sh1, log2FoldChange_sh2) 
  # Selecciono de ambos sh los genes Targets que sean significativos en ambos shRNA
Targets.HM <- BE2C.DE.list %>% dplyr::select(Symbol, log2FoldChange_sh1, log2FoldChange_sh2) %>% 
  subset(Symbol %in% Target.list$Targets) 
rownames(Targets.HM) <- NULL
Targets.HM <- Targets.HM %>% column_to_rownames(var="Symbol") %>% 
  dplyr::rename(SH1 = log2FoldChange_sh1, SH2 = log2FoldChange_sh2)
Target.HM.ann <- as.data.frame(Target.list) %>% subset(Targets %in% rownames(Targets.HM)) 


  # Genero un Heatmap con los targets que son significativos para ambos shRNAs
      # Crear matriz y annotation input para Complex Hatmap
  Targets.HM <- as.matrix(Targets.HM)
  Target.HM.ann <- HeatmapAnnotation(df = Target.HM.ann$Function, name = "Target Function", which= "row")
  
  print( ComplexHeatmap::Heatmap(Targets.HM, 
                                  cluster_columns = FALSE, 
                                  cluster_rows = TRUE,
                                  row_title = "AhR target Expression",
                                  row_title_side = "left",
                                  column_names_rot = 45,
                                  column_names_gp = gpar(fontsize = 12),
                                  row_names_gp = gpar(fontsize = 10), 
                                  right_annotation = Target.HM.ann,
                                  name= "Log2FC", 
                                  col =colorRamp2(c(-4, 0, 4), c("blue", "white", "red")),
                                  heatmap_width = unit(2, "npc") , heatmap_height = unit(8, "npc")
                                  ))

**NOTA: Al comprobar si el mismo script podría utilizarse en la línea celular Kelly, no se ve un efecto tan claro en el ciclo por lo que en un caso real, se deberían combinar kos resultados de ambas líneas celulares para buscar targets que validar funcionalmente.

9. Guardar el entorno de trabajo

save.image(file = "Trabajo Fin Curso R Vicente Fernández Rodero RMD.RData")