Análisis de expresión diferencial - Levadura

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install(version = "3.21")

BiocManager::install("DESeq2")
BiocManager::install("biomaRt")
BiocManager::install("rtracklayer")

library(DESeq2)
library(biomaRt)
library(tidyverse)
library(ggplot2)
library(reshape2)
library(pheatmap)
library(ggrepel)
library(rtracklayer)
library(plotly)
library(dplyr)
library(nord)
library(httr)
library(jsonlite)
library(readr)

Preparando datos de alineamiento con STAR

En primer lugar se toman los archivos de ReadsPerGene que genera STAR y se convierten a un archivo de conteos más simple para luego cargarlos a R cuando se lo requiera.

input_dir <- "../results/mapping/alignments_STAR"
output_dir <- "../results/differential_expression/counts_STAR"

if (!dir.exists(output_dir)) {
  dir.create(output_dir, recursive=TRUE)
}

files <- list.files(input_dir, pattern = "\\ReadsPerGene.out.tab$", full.names = TRUE, recursive = TRUE)

for (archivo in files) {
  # Leemos saltándonos las 4 líneas iniciales y nombrando las 4 columnas
  data <- read_tsv(archivo,
                   skip = 4,
                   col_names = c("Gene_ID", "Unstranded", "Forward", "Reverse"))
  # Construimos el nombre de muestra a partir del nombre de fichero
  sample <- str_split(basename(archivo), "Reads")[[1]][1]
  # Exportamos solo Gene_ID + Unstranded
  out_file <- file.path(output_dir, paste0(sample, "_counts.txt"))
  write_tsv(data[, c("Gene_ID", "Unstranded")],
            out_file,
            col_names = FALSE)
    message("Escrito: ", out_file)
}

Preparando tabla de muestras

Se prepara la tabla de muestras de DESeq2, con tres columnas:

• Nombre de muestras (SampleName): El nombre de cada muestra (batch1, 2 y 3, etc.)

• Nombre de Archivo (FileName): El archivo de counts que se guardó en el paso anterior para cada muestra.

• Condition: Es la condición metabólica de cada muestra. Es importante definirla para el análisis de expresión diferencial posterior.

sampletable <- data.frame(
  SampleName = c("batch1", "batch2", "batch3", 
                 "quimiostato1", "quimiostato2", "quimiostato3"),
  FileName = c("batch1_counts.txt", "batch2_counts.txt", "batch3_counts.txt",
               "quimiostato1_counts.txt", "quimiostato2_counts.txt", "quimiostato3_counts.txt"),
  Condition = c(rep("Batch", 3), rep("Quimiostato", 3))
)

rownames(sampletable) <- gsub("_counts.txt", "", sampletable$FileName)

Importación recuentos STAR

Ahora sí, se importan los archivos de recuentos que se crearon antes y se construye la tabla de DESeq2, con el método DESeqDataSetFromHTSeqCount. Se establece que el diseño es por la condición que definimos en la tabla anterior.

setwd("../results/differential_expression")

se_star <- DESeqDataSetFromHTSeqCount(
  sampleTable = sampletable,
  directory = "counts_STAR",
  design = ~ Condition
)

Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
design formula are characters, converting to factors

Filtrado de genes

En este caso se decidió filtrar por genes que tengan al menos 10 reads.

genes_before <- nrow(se_star)

se_star_before <- se_star

se_star <- se_star[rowSums(counts(se_star)) > 10, ]

genes_after <- nrow(se_star)

df <- data.frame(
  Etapa = c("Antes del filtrado", "Después del filtrado"),
  Genes = c(genes_before, genes_after)
)

ggplot(df, aes(x=Etapa, y=Genes, fill=Etapa)) +
  geom_bar(stat="identity", width=0.5) +
  geom_text(aes(label=Genes), vjust=-0.5, size=5) +
  theme_minimal() +
  labs(title="Número de genes antes y después del filtrado", y="Cantidad de genes") +
  theme(legend.position = "none") +
  ylim(0, max(df$Genes) * 1.1)

Figure 1: Número de genes presentes antes y despúes del filtrado.

Como se puede ver en Figure 1 se pasa de 127 genes a 100.

Importación de anotación

En este caso, se trabajo con un archivo GTF (../data/reference/sacCer_genes.gtf) brindado por la catedra, específico para esta cepa.

# ---- Preparar anotación con archivo GTF ----

# Archivo GTF utilizado en el alineamiento 
gtf_file <- "../data/reference/sacCer_genes.gtf"

gtf <- import(gtf_file)

# Filtro por exon, no por CDS por si tengo algún ARN que no codifica para una proteína
genes_gtf <- gtf[gtf$type == "exon"]

# Convertir a data frame para manejarlo
annot <- as.data.frame(mcols(genes_gtf))

# Seleccionar columnas de interés
annot_df <- data.frame(
  gene_id = annot$gene_id,
  external_gene_name = annot$gene_name,
  chromosome_name = as.character(seqnames(genes_gtf)),
  start_position = start(genes_gtf),
  end_position = end(genes_gtf),
  description = if ("description" %in% colnames(annot)) annot$description else NA
)

dim(annot_df)

[1] 7557    6

knitr::kable(head(mcols(genes_gtf)), caption = "Metadatos disponibles (primeras filas)")

Metadatos disponibles (primeras filas)

source	type	score	phase	gene_id	transcript_id	exon_number	gene_name	p_id	protein_id	transcript_name	tss_id
protein_coding	exon	NA	NA	R0010W	R0010W	1	FLP1	P6650	NA	FLP1	TSS1576
protein_coding	exon	NA	NA	R0020C	R0020C	1	REP1	P2013	NA	REP1	TSS6001
protein_coding	exon	NA	NA	R0030W	R0030W	1	RAF1	P3033	NA	RAF1	TSS6681
protein_coding	exon	NA	NA	R0040C	R0040C	1	REP2	P1202	NA	REP2	TSS1452
protein_coding	exon	NA	NA	YAL069W	YAL069W	1	YAL069W	P1259	NA	YAL069W	TSS1092
protein_coding	exon	NA	NA	YAL068W-A	YAL068W-A	1	YAL068W-A	P3321	NA	YAL068W-A	TSS5453

Ajuste modelo estadístico DESeq2

Para el análisis de la expresión diferencial se va a usar DESeq2, por ello se procede a ajustar el modelo, con los datos ya transformados y filtrados previamente.

Se calcula el conteo normalizado de los datos, con una transformación log2 + 1 y finalmente se agrega la anotación que se cargo desde el GTF.

Estos datos son en suma almacenados en normalized_counts_log2_star.txt

#---- Ajustar modelo estadístico DESeq2----

se_star2 <- DESeq(se_star)

estimating size factors

estimating dispersions

gene-wise dispersion estimates

mean-dispersion relationship

final dispersion estimates

fitting model and testing

#Calcular conteo normalizado (transformación log2); + 1
norm_counts <- log2(counts(se_star2,
                           normalized = TRUE)+1)
#Agregar anotación
norm_counts_symbols <- merge(data.frame(ID=rownames(norm_counts),
                                        norm_counts,
                                        check.names=FALSE),
                             annot,
                             by.x="ID",
                             by.y="gene_id",
                             all=F)

#Escribir recuentos normalizados en un archivo
write.table(norm_counts_symbols, "normalized_counts_log2_star.txt", quote=F, col.names=T, row.names=F, sep="\t")

Análisis de expresión diferencial

Se extraen los resultados de la expresión diferencial con el contraste Batch vs Condition.

#----Análisis de expresión diferencial----

#Chequeo de nombre de resultados: depende de qué se haya modelado. Aquí se modeló la "Condition"
resultsNames(se_star2)

[1] "Intercept"                      "Condition_Quimiostato_vs_Batch"

#Extraer resultados para Batch vs Quimiostato
de <- results(se_star2,
              contrast = c("Condition", "Batch", "Quimiostato"))


#Chequeo de filas
head(de)

log2 fold change (MLE): Condition Batch vs Quimiostato 
Wald test p-value: Condition Batch vs Quimiostato 
DataFrame with 6 rows and 6 columns
           baseMean log2FoldChange     lfcSE       stat       pvalue
          <numeric>      <numeric> <numeric>  <numeric>    <numeric>
YAL064W-B    2.1466      -0.419267  1.231087  -0.340566  7.33430e-01
YAL063C     16.4781      -0.707654  0.464899  -1.522168  1.27967e-01
YAL062W    677.8951      -4.079542  0.131949 -30.917537 6.94246e-210
YAL061W    599.0764      -3.001234  0.135487 -22.151453 1.01008e-108
YAL060W   1644.8971       0.160265  0.103134   1.553956  1.20195e-01
YAL059W    115.3256       1.015146  0.192330   5.278144  1.30499e-07
                  padj
             <numeric>
YAL064W-B           NA
YAL063C    1.70623e-01
YAL062W   3.33238e-208
YAL061W   1.93935e-107
YAL060W    1.62517e-01
YAL059W    3.68468e-07

#Agregar la anotación
de_symbols <- merge(data.frame(ID=rownames(de),
                               de,
                               check.names=FALSE),
                    annot,
                    by.x="ID",
                    by.y="gene_id",
                    all=F)
#Escribir los resultados de expresión diferencial en el archivo
write.table(de_symbols,
            "deseq2_results.txt",
            quote=F,
            col.names=T,
            row.names=F,
            sep="\t")

Se seleccionan los genes diferencialmente expresados con el criterio:

• \(p_adj < 0.05\)

• \(|log_2(FoldChange)| > 0.5\)

#----Selección de genes----
de_select <- de_symbols[de_symbols$padj < 0.05 & !is.na(de_symbols$padj) & abs(de_symbols$log2FoldChange) > 0.5,]

write.table(de_select,
            "deseq2_selection.txt",
            quote=F,
            col.names=T,
            row.names=F,
            sep="\t")

Visualización

#---- Visualización----
# Transformar los recuentos sin procesar para poder visualizar los datos
se_rlog <- rlog(se_star2)

Correlación de muestras

#---- Visualización - Correlación de muestras----

#Calcular la matriz de distancia entre muestras
sampleDistMatrix <- as.matrix(dist(t(assay(se_rlog))))
#Preparar metadata
metadata <- data.frame(Condition = sampletable$Condition)
rownames(metadata) <- sampletable$SampleName

pheatmap(sampleDistMatrix,
         annotation_col = metadata,
         main = "Matriz de distancia entre muestras (transformación rlog)",
         fontsize = 10,
         fontsize_row = 10,
         fontsize_col = 10,
         border_color = NA,
         color = colorRampPalette(rev(RColorBrewer::brewer.pal(n = 9, name = "RdBu")))(100),
         angle_col = "45")

Figure 2: Matriz de distancia para las 6 muestras.

Como se puede ver en Figure 2, las muestras se agrupan correctamente según su condición, en la jerarquía también se pueden ver correctamente diferenciadas.

PCA

# ---- Visualización - PCA personalizado sin líneas en etiquetas ----

# Obtener datos de PCA con DESeq2
pcaData <- plotPCA(se_rlog, intgroup = "Condition", returnData = TRUE)

using ntop=500 top features by variance

# Calcular porcentaje de varianza explicada
percentVar <- round(100 * attr(pcaData, "percentVar"))

# Graficar
ggplot(pcaData, aes(x = PC1, y = PC2, color = Condition, label = name)) +
  geom_point(size = 4, alpha = 0.85) +
  geom_text_repel(size = 3.5, segment.color = NA) +
  scale_color_brewer(palette = "Set1") +
  labs(
    title = "Análisis de componentes principales (PCA)",
    x = paste0("PC1: ", percentVar[1], "% varianza"),
    y = paste0("PC2: ", percentVar[2], "% varianza")
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    panel.grid = element_blank(),
    panel.border = element_rect(color = "black", fill = NA),
    legend.position = "bottom"
  )

Figure 3: PCA para los datos de DESeq2 diferenciando por condición

Como se puede observar en la Figure 3, los grupos se encuentran bien diferenciados en las dos condiciones en el eje de la primer componente principal. Por otro lado, en el eje de PC2 la varianza es muy baja.

Volcano Plot

# Buscar el mínimo pvalue positivo (no cero)
min_nonzero_p <- min(de_symbols$pvalue[de_symbols$pvalue > 0], na.rm = TRUE)

# Reemplazar pvalues 0 por un valor pequeño (por ejemplo min_nonzero_p / 10)
de_symbols <- de_symbols %>%
  mutate(pvalue = ifelse(pvalue == 0, min_nonzero_p / 10, pvalue))



de_clean <- de_symbols %>%
  filter(is.finite(log2FoldChange), is.finite(pvalue)) %>%
  mutate(diffexpressed = case_when(
    log2FoldChange > 0.5 & pvalue < 0.05 ~ "UP",
    log2FoldChange < -0.5 & pvalue < 0.05 ~ "DOWN",
    TRUE ~ "NO"
  )) %>%
  mutate(diffexpressed = factor(diffexpressed, levels = c("DOWN", "NO", "UP")))

# Etiquetas para genes DE
de_clean$delabel <- NA
idx <- which(de_clean$diffexpressed != "NO" & !is.na(de_clean$ID))
de_clean$delabel[idx] <- de_clean$ID[idx]

plotly_volcano <- de_clean %>%
  plot_ly(
    x = ~log2FoldChange,
    y = ~-log10(pvalue),
    color = ~diffexpressed,
    colors = c("blue", "black", "red"),
    text = ~paste("Gene ID:", ID,
                  "<br>log2FC:", round(log2FoldChange, 2),
                  "<br>p-value:", signif(pvalue, 3)),
    hoverinfo = "text",
    type = "scatter",
    mode = "markers"
  ) %>%
  layout(
    title = "Volcano plot interactivo",
    xaxis = list(title = "log2 Fold Change"),
    yaxis = list(title = "-log10(p-value)"),
    legend = list(title = list(text = "Expresión"))
  ) %>%
  add_segments(x = -0.5, xend = -0.5,
               y = 0, yend = max(-log10(de_clean$pvalue), na.rm = TRUE),
               line = list(dash = "dash", color = "red"), showlegend = FALSE) %>%
  add_segments(x = 0.5, xend = 0.5,
               y = 0, yend = max(-log10(de_clean$pvalue), na.rm = TRUE),
               line = list(dash = "dash", color = "red"), showlegend = FALSE) %>%
  add_segments(x = min(de_clean$log2FoldChange), xend = max(de_clean$log2FoldChange),
               y = -log10(0.05), yend = -log10(0.05),
               line = list(dash = "dash", color = "red"), showlegend = FALSE)

plotly_volcano

Figure 4: Volcano plot para los genes diferencialmente expresados con el tratamiento Batch vs. Quimiostato

gg_volcano <- ggplot(de_clean, aes(x = log2FoldChange, y = -log10(pvalue), color = diffexpressed, label = delabel)) +
  geom_point(alpha = 0.6, size = 2) +
  geom_text_repel(segment.color = "grey50", max.overlaps = 20, size = 3.5) +
  scale_color_manual(values = c("DOWN" = "blue", "NO" = "black", "UP" = "red")) +
  geom_vline(xintercept = c(-0.5, 0.5), color = "red", linetype = "dashed") +
  geom_hline(yintercept = -log10(0.05), color = "red", linetype = "dashed") +
  labs(
    title = "Volcano Plot",
    subtitle = "log2 Fold Change vs -log10(p-value)",
    x = "log2FoldChange",
    y = "-log10(p-value)",
    color = "Expresión"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    legend.position = "bottom",
    # plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5),
    panel.grid = element_blank(),
    panel.border = element_rect(color = "black", fill = NA, size = 1)
  )

Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
ℹ Please use the `linewidth` argument instead.

print(gg_volcano)

Warning: Removed 49 rows containing missing values or values outside the scale range
(`geom_text_repel()`).

Warning: ggrepel: 37 unlabeled data points (too many overlaps). Consider
increasing max.overlaps

Figure 5: Volcano plot para los genes diferencialmente expresados con el tratamiento Batch vs. Quimiostato

En Figure 5 se pueden apreciar los genes diferencialmente expresados para las muestras, con el tratamiento Batch vs. Condition. A la derecha podemos ver los genes que están expresados positivamente en Batch, a la izquierda el caso contrario. En el eje y (\(-log10(pvalue)\)) podemos ver también la confianza que se tiene de esa expresión diferencial. Además tenemos ciertos valores umbrales (las líneas interlineadas rojas) que se tienen en cuenta para ser considerados como diferencialmente expresados.

Volcano plot para genes que no codifican proteínas

de_clean_no_p <- de_clean %>% 
  filter(is.na(p_id))

gg_volcano <- ggplot(de_clean_no_p, aes(x = log2FoldChange, y = -log10(pvalue), color = diffexpressed, label = delabel)) +
  geom_point(alpha = 0.6, size = 2) +
  geom_text_repel(segment.color = "grey50", max.overlaps = 20, size = 3.5) +
  scale_color_manual(values = c("DOWN" = "blue", "NO" = "black", "UP" = "red")) +
  geom_vline(xintercept = c(-0.5, 0.5), color = "red", linetype = "dashed") +
  geom_hline(yintercept = -log10(0.05), color = "red", linetype = "dashed") +
  labs(
    title = "Volcano Plot",
    subtitle = "log2 Fold Change vs -log10(p-value)",
    x = "log2FoldChange",
    y = "-log10(p-value)",
    color = "Expresión"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    legend.position = "bottom",
    # plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5),
    panel.grid = element_blank(),
    panel.border = element_rect(color = "black", fill = NA, size = 1)
  )

print(gg_volcano)

Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_text_repel()`).

Figure 6: Volcano plot para los genes diferencialmente expresados que no codifican para proteínas con el tratamiento Batch vs. Quimiostato

En este caso se buscó graficar genes que hayan sido contados y no codifiquen para proteínas (se tomo que no tenga p_id como criterio). Sólo se observa snR18 como up regulado, al inspeccionar en el explorador genómico (Figure 7) se ve que se sitúa en un intrón de YAL003W (EFB1) una de las regiones con más cobertura, por lo que puede ser esa la razón del ruido.

Figure 7: Cromosoma uno y los archivos de cobertura para cada muestra, remarcado en rojo EFB1 y dentro snR18

Enriquecimiento funcional con YeastEnrichr

Reajusto los criterios para que sean más estrictos:

• \(padj < 0.05\)
• \(|log_2(FC)| > 1\) (al menos duplicación/división por 2)

Con estos criterior construyo dos listas, los genes upregulados \(log_2(FC) > 1\) y downregulados \(log_2(FC) < -1\)

# ---- Selección precisa de genes DE para análisis funcional ----

# Asegurar que las columnas necesarias estén presentes
de_symbols_clean <- de_symbols %>%
  filter(!is.na(padj), !is.na(log2FoldChange), !is.na(gene_name))

padj_corte <- 0.05
log2FoldChange_corte <- 1

# Aplicar filtros estrictos
de_filtrados_up <- de_symbols_clean %>%
  filter(
    padj < padj_corte,
    log2FoldChange > log2FoldChange_corte
  )

de_filtrados_down <- de_symbols_clean %>% 
  filter(
    padj < padj_corte,
    log2FoldChange < -log2FoldChange_corte
  )

# Extraer solo nombres únicos de genes conocidos
genes_interes_up <- de_filtrados_up %>%
  distinct(gene_name) %>%
  pull(gene_name)

genes_interes_down <- de_filtrados_down %>%
  distinct(gene_name) %>%
  pull(gene_name)

# Verificamos
length(genes_interes_up)

[1] 9

head(genes_interes_up)

[1] "ERP2"      "CYS3"      "PMT2"      "MAK16"     "MYO4"      "YAL037C-A"

length(genes_interes_down)

[1] 14

head(genes_interes_down)

[1] "SSA1"    "FUN14"   "YAL018C" "SNC1"    "FUN19"   "ACS1"   

Se lleva a cabo un enriquecimiento con YeastEnrichr. Esto se realiza a través de su API con el siguiente código. Primero hay que enviar la lista de genes que nos interesa analizar. Una vez que

Upregulated

Envío de lista de genes de interés

# ---- Paso 1: Enviar la lista ----
enviar_a_enrichr <- function(lista_genes, descripcion = "Lista de genes") {
  ENRICHR_ADD_URL <- "http://maayanlab.cloud/YeastEnrichr/addList"
  
  genes_str <- paste(lista_genes, collapse = "\n")
  
  response <- POST(
    url = ENRICHR_ADD_URL,
    body = list(
      list = genes_str,
      description = descripcion
    ),
    encode = "multipart"
  )
  
  if (status_code(response) != 200) {
    stop("Error al analizar la lista de genes")
  }
  
  # Extraer respuesta como texto
  raw_text <- as.character(content(response, as = "text", encoding = "UTF-8"))
  
  # Extraer JSON dentro del texto HTML
  json_text <- sub(".*(\\{.*\\}).*", "\\1", raw_text)
  
  # Parsear JSON
  parsed_result <- fromJSON(json_text)
  
  return(list(
    userListId = parsed_result$userListId,
    shortId = parsed_result$shortId
  ))
}


resultado_up <- enviar_a_enrichr(genes_interes_up, "Lista de genes Upregulated")
resultado_down <- enviar_a_enrichr(genes_interes_down, "Lista de genes Downregulated")

Obtención de resultados

En este caso se corre el enriquecimiento con YeastEnrichr. Se define la función en R para poder correrlo pasando como argumento, el token de la lista de genes que obtuvimos en el paso anterior y el tipo de set de datos contra el que queremos que se contraste. Las bibliotecas disponibles pueden consultarse en YeastEnrichr. En este caso, se enriqueció para:

• GO_Biological_Process_2018
• GO_Cellular_Component_2018
• GO_Molecular_Function_2018
• KEGG_2018

Se seleccionaron estas librerías porque son recursos de Gene Ontology o KEGG con una versión bien establecida, otorgando reproducibilidad.

library(httr)
library(readr)

# Función general para correr enriquecimiento y guardar resultados a un `tsv`
run_enrichr <- function(userListId, gene_set_library, output_dir = "../results/functional_enrichment") {
  ENRICHR_ENRICH_URL <- "http://maayanlab.cloud/YeastEnrichr/enrich"
  
  if (!dir.exists(output_dir)) {
    dir.create(output_dir, recursive=TRUE)
  }
  
  query_url <- paste0(ENRICHR_ENRICH_URL, "?userListId=", userListId,
                      "&backgroundType=", gene_set_library)
  
  res <- GET(query_url)
  if (http_error(res)) stop("❌ Error al obtener los resultados de enriquecimiento para ", gene_set_library)
  
  data <- content(res, as = "parsed")[[gene_set_library]]
  
  if (length(data) == 0) {
    warning("⚠️ No se obtuvieron resultados para ", gene_set_library)
    return(NULL)
  }
  
  data_flat <- lapply(data, function(x) {
    x[[6]] <- paste(unlist(x[[6]]), collapse = ", ")
    unlist(x)
  })
  
  enrichment_df <- as.data.frame(do.call(rbind, data_flat), stringsAsFactors = FALSE)
  colnames(enrichment_df) <- c("Rank", "Term", "P-value", "Z-score", "Combined Score", 
                               "Overlapping Genes", "Adjusted P-value", "Old P", "Old Adj P")
  
  output_file <- file.path(output_dir, paste0("enrichment_", gene_set_library, ".tsv"))
  write_tsv(enrichment_df, output_file)
  
  return(enrichment_df)
}

gene_libraries <- c("GO_Biological_Process_2018", "GO_Cellular_Component_2018", "GO_Molecular_Function_2018", "KEGG_2018")

results_up <- lapply(gene_libraries, function(lib) {
  run_enrichr(resultado_up$userListId, lib, "../results/functional_enrichment/up_regulated")
})

Análisis de resultados

En primer lugar, defino funciones básicas que me permitan generalizar la carga de los datos, transformación y la generación de gráficos para las distintas librerías.

# Función para cargar y preprocesar el archivo de enriquecimiento
load_enrichment <- function(lib_name, base_dir = "../results/functional_enrichment", top_n = 15, wrap_width = 60) {
  filepath <- file.path(base_dir, paste0("enrichment_", lib_name, ".tsv"))
  df <- read_tsv(filepath, show_col_types = FALSE) %>%
    mutate(
      `Adjusted P-value` = as.numeric(`Adjusted P-value`),
      `Combined Score` = as.numeric(`Combined Score`),
      `Overlapping Genes` = str_count(`Overlapping Genes`, ",") + 1,
      Term_wrapped = str_wrap(Term, width = wrap_width)
    ) %>%
    arrange(`Adjusted P-value`) %>%
    slice_head(n = top_n) %>%
    mutate(Term_wrapped = fct_reorder(Term_wrapped, `Combined Score`))
  return(df)
}

# Función para hacer barplot
plot_barplot <- function(enrich_df, lib_name) {
  lib_name_clean <- str_replace_all(lib_name, "_", " ")
  title_text <- paste("Top términos enriquecidos\n", lib_name_clean)
  
  my_colors <- nord::nord_palettes$afternoon_prarie
  my_colors <- nord("aurora", 20)

  ggplot(enrich_df, aes(x = Term_wrapped, y = `Combined Score`, fill= as.factor(`Combined Score`))) +
    geom_col(width = 0.7) +
    coord_flip() +
    labs(
      title = title_text,
      x = "Término GO",
      y = "Combined Score"
    ) +
    scale_fill_manual(values = my_colors) +
    theme_minimal(base_size = 12) +
    theme(
      axis.text.y = element_text(size = 8, margin = margin(r = 10)),
      axis.ticks.y = element_blank(),
      plot.title = element_text(
        size = 14,
        margin = margin(t = 20, b = 10),
        hjust = 0.6  # Ajustá este valor entre 0.5 y 0.7 según convenga
      ),
      legend.position="none"
    ) 
}

# Función para hacer dotplot
plot_dotplot <- function(enrich_df, lib_name) {
  lib_name_clean <- str_replace_all(lib_name, "_", " ")
  title_text <- paste("Top términos enriquecidos\n", lib_name_clean)
  
  enrich_df <- enrich_df %>%
    mutate(`Adjusted P-value` = ifelse(`Adjusted P-value` == 0, 1e-300, `Adjusted P-value`))
  
  ggplot(enrich_df, aes(x = `Combined Score`, y = Term_wrapped)) +
    geom_point(aes(color = `Adjusted P-value`, size = `Overlapping Genes`)) +
    scale_color_gradient(trans = "log10", low = "red", high = "blue", 
                         name = "Adjusted P-value") +
    scale_size_continuous(range = c(2, 8), name = "Cantidad genes") +
    labs(
      title = title_text
      ,
      x = "Combined Score",
      y = "Término GO"
    ) +
    theme_minimal(base_size = 12) +
    theme(
      axis.text.y = element_text(size = 8, margin = margin(r = 10)),
      axis.ticks.y = element_blank()
    )
}

GO Biological Process 2018

library_name <- "GO_Biological_Process_2018"
enrich_df <- load_enrichment(library_name, base_dir = "../results/functional_enrichment/up_regulated")

barplot <- plot_barplot(enrich_df, library_name)
dotplot <- plot_dotplot(enrich_df, library_name)

print(barplot)

Figure 8: Barplot of Biological processes from GO 2018 enriched

print(dotplot)

Figure 9: Dotplot of Biological processes from GO 2018 enriched

GO Cellular Component 2018

library_name <- "GO_Cellular_Component_2018"
enrich_df <- load_enrichment(library_name, base_dir = "../results/functional_enrichment/up_regulated")

barplot <- plot_barplot(enrich_df, library_name)
dotplot <- plot_dotplot(enrich_df, library_name)

print(barplot)

Figure 10: Barplot of Cellular component from GO 2018 enriched

print(dotplot)

Figure 11: Dotplot of Cellular component from GO 2018 enriched

GO Molecular Function 2018

library_name <- "GO_Molecular_Function_2018"
enrich_df <- load_enrichment(library_name, base_dir = "../results/functional_enrichment/up_regulated")

barplot <- plot_barplot(enrich_df, library_name)
dotplot <- plot_dotplot(enrich_df, library_name)

print(barplot)

Figure 12: Barplot of Molecular function from GO 2018 enriched

print(dotplot)

Figure 13: Dotplot of Molecular function from GO 2018 enriched

KEGG pathways

library_name <- "KEGG_2018"
enrich_df <- load_enrichment(library_name, base_dir = "../results/functional_enrichment/up_regulated")

barplot <- plot_barplot(enrich_df, library_name)
dotplot <- plot_dotplot(enrich_df, library_name)

print(barplot)

Figure 14: Barplot de rutas metabólicas en KEGG

print(dotplot)

Figure 15: Dotplot de rutas metabólicas en KEGG

Downregulated

results_down <- lapply(gene_libraries, function(lib) {
  run_enrichr(resultado_down$userListId, lib, "../results/functional_enrichment/down_regulated")
})

GO Biological Process 2018

library_name <- "GO_Biological_Process_2018"
enrich_df <- load_enrichment(library_name, base_dir = "../results/functional_enrichment/down_regulated")

barplot <- plot_barplot(enrich_df, library_name)
dotplot <- plot_dotplot(enrich_df, library_name)

print(barplot)

Figure 16: Barplot of Biological processes from GO 2018 enriched

print(dotplot)

Figure 17: Dotplot of Biological processes from GO 2018 enriched

GO Cellular Component 2018

library_name <- "GO_Cellular_Component_2018"
enrich_df <- load_enrichment(library_name, base_dir = "../results/functional_enrichment/down_regulated")

barplot <- plot_barplot(enrich_df, library_name)
dotplot <- plot_dotplot(enrich_df, library_name)

print(barplot)

Figure 18: Barplot of Cellular component from GO 2018 enriched

print(dotplot)

Figure 19: Dotplot of Cellular component from GO 2018 enriched

GO Molecular Function 2018

library_name <- "GO_Molecular_Function_2018"
enrich_df <- load_enrichment(library_name, base_dir = "../results/functional_enrichment/down_regulated")

barplot <- plot_barplot(enrich_df, library_name)
dotplot <- plot_dotplot(enrich_df, library_name)

print(barplot)

Figure 20: Barplot of Molecular function from GO 2018 enriched

print(dotplot)

Figure 21: Dotplot of Molecular function from GO 2018 enriched

KEGG pathways

library_name <- "KEGG_2018"
enrich_df <- load_enrichment(library_name, base_dir = "../results/functional_enrichment/down_regulated")

barplot <- plot_barplot(enrich_df, library_name)
dotplot <- plot_dotplot(enrich_df, library_name)

print(barplot)

Figure 22: Barplot de rutas metabólicas en KEGG

print(dotplot)

Figure 23: Dotplot de rutas metabólicas en KEGG

Enriquecimiento funcional extendido

Como se aclaró previamente, los datos con los que se ha trabajado hasta el momento son correspondientes a lecturas que alinean con el primer cromosoma. Por ello, se llevó a cabo un enriquecimiento funcional de una lista de genes diferencialmente expresados provista por la cátedra. A diferencia del análisis anterior, aquí no se cuenta con información sobre si están downregulados o upregulados los genes.

Se utilizó la misma metodología de llamados a la API de YeastEnrichr.

genes_extendidos = readLines('../data/genes_diff_exp/Genes_diferencialmente_expresados_signif.tabular')

resultado_extendido <- enviar_a_enrichr(genes_extendidos, "Lista de genes extendida")

results_up <- lapply(gene_libraries, function(lib) {
  run_enrichr(resultado_extendido$userListId, lib, "../results/functional_enrichment/extended_analysis")
})

ruta <- "../results/functional_enrichment/extended_analysis"

# Leer los nombres de los archivos TSV
archivos <- list.files(path = ruta, pattern = "\\.tsv$", full.names = TRUE)

# Leer y graficar
for (librarie in gene_libraries) {
  enrichment_df <- load_enrichment(librarie, base_dir="../results/functional_enrichment/extended_analysis")
  
  enrichment_df <- head(enrichment_df, 5)

  # Hacer el gráfico
  p <- plot_barplot(enrichment_df, librarie)
  
  print(p)
}