Descripción del equipo

Contexto

Este entregable continúa a partir del Entregable 3, donde se realizó el control de calidad con FastQC/MultiQC, la eliminación de adaptadores con Trimmomatic, el pseudoalineamiento con Kallisto y el alineamiento con STAR, seguidos de la importación de datos en R mediante tximport. Los objetos dds (Kallisto) y dds_star (STAR) fueron construidos con DESeq2, filtrados (≥10 cuentas en ≥3 muestras) y con diseño ~ genotype (referencia: WT).

En este entregable realizamos:

  1. Normalización de cuentas (VST, z-score, melted counts)
  2. Análisis exploratorio de componentes principales (PCA) y corrección de batch effect
  3. Análisis de expresión diferencial con DESeq2 y shrinkage de log2FC
  4. Visualización: volcano, MA plot, heatmap
  5. Comparación Kallisto vs STAR
  6. Análisis funcional de términos GO con gprofiler2

0. Cargar paquetes

library(DESeq2)        # Normalización y análisis de expresión diferencial
library(tximport)      # Importar datos de Kallisto
library(rhdf5)         # Leer archivos .h5 de Kallisto
library(ggplot2)       # Visualización
library(dplyr)         # Manipulación de datos
library(reshape2)      # Conversión wide → long (melt)
library(pheatmap)      # Heatmaps
library(RColorBrewer)  # Paletas de colores
library(knitr)         # Tablas en el documento
library(kableExtra)    # Formato de tablas
library(ggrepel)       # Etiquetas sin superposición
library(gprofiler2)    # Análisis de enriquecimiento funcional
library(VennDiagram)   # Diagramas de Venn
library(grid)          # Para dibujar Venn
library(limma)         # Corrección de batch effect
library(apeglm)        # Prior para lfcShrink

1. Cargar datos y reconstruir objetos DESeq2

# Rutas
ruta_kallisto <- "/Users/dubsirubs/Desktop/RNAseq_Equipo4/results_final/kallisto"
star_dir      <- "/Users/dubsirubs/Desktop/RNAseq_Equipo4/results_final/STAR"
outdir        <- "/Users/dubsirubs/Desktop/RNAseq_Equipo4/results_final/"
figdir        <- "/Users/dubsirubs/Desktop/RNAseq_Equipo4/results_final/figures/"

samples <- c("SRR33445252","SRR33445253","SRR33445254",
             "SRR33445260","SRR33445261","SRR33445262")

# Metadata
metadata <- read.csv("/Users/dubsirubs/Desktop/RNAseq_Equipo4/metadata.csv")
metadata$genotype <- factor(metadata$genotype, levels = c("WT", "Dp16"))
rownames(metadata) <- metadata$sample

# Kallisto → tximport → DESeq2
archivos_h5 <- file.path(ruta_kallisto, samples, "abundance.h5")
names(archivos_h5) <- samples

tx2gene <- read.csv("/Users/dubsirubs/Desktop/RNAseq_Equipo4/results_final/tx2gene_kallisto.csv")

txi.kallisto <- tximport(archivos_h5,
                         type            = "kallisto",
                         tx2gene         = tx2gene,
                         txOut           = FALSE,
                         ignoreTxVersion = FALSE,
                         ignoreAfterBar  = TRUE)

dds <- DESeqDataSetFromTximport(txi.kallisto,
                                colData = metadata,
                                design  = ~ genotype)
# Filtro: genes con ≥10 cuentas en al menos 3 muestras
keep <- rowSums(counts(dds) >= 10) >= 3
dds  <- dds[keep, ]

# STAR → DESeq2
count_list  <- lapply(samples, function(s) {
  f <- file.path(star_dir, s, paste0(s, ".ReadsPerGene.out.tab"))
  read.table(f, skip = 4, header = FALSE)[, c(1, 2)]
})
counts_star <- do.call(cbind, lapply(count_list, function(x) x[, 2]))
rownames(counts_star) <- count_list[[1]][, 1]
colnames(counts_star) <- samples

dds_star <- DESeqDataSetFromMatrix(countData = counts_star,
                                   colData   = metadata,
                                   design    = ~ genotype)
keep_star <- rowSums(counts(dds_star) >= 10) >= 3
dds_star  <- dds_star[keep_star, ]

cat("Genes retenidos — Kallisto:", nrow(dds), "| STAR:", nrow(dds_star), "\n")
## Genes retenidos — Kallisto: 14715 | STAR: 14091

2. Normalización de las cuentas

Las cuentas crudas de RNA-seq no son comparables entre muestras porque dependen del tamaño total de la librería y, en el caso de Kallisto, también de la longitud de los transcritos. La VST (Variance Stabilizing Transformation) estabiliza la varianza a lo largo de todo el rango de expresión: los genes de baja expresión dejan de tener una varianza artificialmente alta, y los datos quedan en una escala apta para PCA, heatmaps y clustering. Se usa blind = FALSE porque conocemos el diseño experimental — esto permite que la transformación tome en cuenta el genotipo, lo cual es correcto cuando el objetivo es visualización comparativa.

2.1 VST y Kallisto

vsd_kallisto <- vst(dds, blind = FALSE)
cat("Dimensiones VST Kallisto (genes × muestras):", dim(vsd_kallisto), "\n")
## Dimensiones VST Kallisto (genes × muestras): 14715 6

2.2 VST y STAR

vsd_star <- vst(dds_star, blind = FALSE)
cat("Dimensiones VST STAR (genes × muestras):", dim(vsd_star), "\n")
## Dimensiones VST STAR (genes × muestras): 14091 6

2.3 Distribución de cuentas normalizadas

El boxplot muestra la distribución de los valores VST por muestra. Si la normalización fue exitosa, las medianas de todas las cajas deben estar al mismo nivel, indicando que las muestras son comparables entre sí independientemente del tamaño de su librería.

melted_norm_counts <- data.frame(melt(assay(vsd_kallisto)))
colnames(melted_norm_counts) <- c("gene", "sample", "normalized_counts")

melted_norm_counts <- left_join(
  melted_norm_counts,
  data.frame(sample = rownames(metadata), genotype = metadata$genotype),
  by = "sample"
)

ggplot(melted_norm_counts,
       aes(x = sample, y = normalized_counts, fill = genotype)) +
  geom_boxplot(outlier.size = 0.3, alpha = 0.8) +
  scale_fill_manual(values = c("WT" = "steelblue", "Dp16" = "firebrick3")) +
  labs(title = "Distribución de cuentas normalizadas (VST) por muestra",
       x = "Muestra", y = "Cuentas normalizadas (VST)", fill = "Genotipo") +
  theme_bw(base_size = 12) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        plot.background = element_rect(fill = "white"))
Figura 1. Distribución de cuentas normalizadas (VST, Kallisto) por muestra. Las medianas alineadas indican que la normalización fue exitosa y las muestras son comparables.

Figura 1. Distribución de cuentas normalizadas (VST, Kallisto) por muestra. Las medianas alineadas indican que la normalización fue exitosa y las muestras son comparables.

Interpretación: La Figura 1 muestra que las medianas de las seis muestras se encuentran alineadas en torno a 9, lo que confirma que la normalización VST corrigió exitosamente las diferencias de tamaño de librería. Los rangos intercuartílicos son homogéneos entre muestras, sin outliers sistemáticos que indiquen problemas técnicos. Los puntos extremos superiores (~17–21) corresponden a genes altamente expresados presentes en todas las muestras por igual, lo que refleja una característica biológica esperada. La ausencia de separación sistemática entre genotipos (WT vs Dp16) en esta distribución global es correcta: las diferencias específicas por gen se detectarán en el análisis diferencial.

2.4 Z-score

El z-score estandariza cada gen a media = 0 y desviación estándar = 1 a través de todas las muestras. Esto elimina las diferencias de escala absoluta entre genes (algunos tienen expresión basalmente alta y otros baja) y resalta los patrones de expresión relativa, lo que facilita el clustering y la visualización en heatmaps. Se aplica t(scale(t(x))) porque scale() opera por columnas — al transponer dos veces, escalamos por filas (por gen).

counts_raw <- counts(dds)

# Eliminar genes con varianza = 0 (misma expresión en todas las muestras)
var_non_zero    <- apply(counts_raw, 1, var) != 0
filtered_counts <- counts_raw[var_non_zero, ]

# t(scale(t(x))): transponer → escalar por gen → transponer de vuelta
zscores    <- t(scale(t(filtered_counts)))
zscore_mat <- as.matrix(zscores)

cat("Genes con varianza > 0:", nrow(zscore_mat), "\n")
## Genes con varianza > 0: 14715
cat("Muestras:", ncol(zscore_mat), "\n")
## Muestras: 6

De los ~20,000+ genes del genoma de ratón, 14,715 muestran al menos algo de variación entre las 6 muestras. Esto es normal ya que muchos genes housekeeping se expresan de forma constitutiva e idéntica en todas las muestras y quedan fuera. Estos 14,715 son los informativos para clustering y heatmaps.

3. Análisis exploratorio: PCA

El PCA reduce la dimensionalidad de los datos (miles de genes) a componentes que capturan la mayor variación. Su propósito es verificar que las muestras agrupan por la variable biológica de interés (genotipo) y no por variables técnicas (fecha, técnico, equipo de secuenciación), lo que indicaría un batch effect.

3.1 PCA y Kallisto

mat_k <- assay(vsd_kallisto)
pc_k  <- prcomp(t(mat_k))   # t() porque prcomp espera muestras en filas

var_exp_k <- pc_k$sdev^2 / sum(pc_k$sdev^2) * 100

pca_df_k <- data.frame(
  PC1      = pc_k$x[, 1],
  PC2      = pc_k$x[, 2],
  genotype = metadata$genotype,
  sample   = rownames(metadata)
)

ggplot(pca_df_k, aes(x = PC1, y = PC2, color = genotype, label = sample)) +
  geom_point(size = 4) +
  geom_text_repel(size = 3, show.legend = FALSE) +
  stat_ellipse(aes(group = genotype), linetype = "dashed", level = 0.8) +
  scale_color_manual(values = c("WT" = "steelblue", "Dp16" = "firebrick3")) +
  labs(title = "PCA — Kallisto (VST)",
       x     = sprintf("PC1 (%.1f%% varianza)", var_exp_k[1]),
       y     = sprintf("PC2 (%.1f%% varianza)", var_exp_k[2]),
       color = "Genotipo") +
  theme_minimal(base_size = 13) +
  theme(plot.background = element_rect(fill = "white"))
Figura 2. PCA de muestras normalizadas con VST (Kallisto). Cada punto representa una muestra; el color indica el genotipo (azul = WT, rojo = Dp16). Las elipses muestran la dispersión de cada grupo.

Figura 2. PCA de muestras normalizadas con VST (Kallisto). Cada punto representa una muestra; el color indica el genotipo (azul = WT, rojo = Dp16). Las elipses muestran la dispersión de cada grupo.

Interpretación: La Figura 2 muestra una separación clara entre genotipos a lo largo de PC1(43.9% de la varianza total): las tres réplicas WT se agrupan a valores negativos y las tres réplicas Dp16 a valores positivos. Esto indica que el genotipo es la principal fuente de variación en los datos, lo cual es la señal biológica esperada en un modelo de trisomía. PC2 (24.9%) captura variabilidad dentro de los grupos, consistente con la variación biológica natural entre réplicas individuales con n = 3. No se observa ninguna muestra fuera de su grupo genotípico, lo que descarta batch effects evidentes en este primer análisis exploratorio.

3.2 Scree Plot

El scree plot muestra qué porcentaje de la varianza total captura cada componente principal. Es útil para saber cuántos PCs son informativos: un “codo” pronunciado indica que los primeros PCs concentran la mayor parte de la información.

df_scree <- data.frame(
  Componente = seq_along(var_exp_k),
  Varianza   = var_exp_k,
  Etiqueta   = paste0("PC", seq_along(var_exp_k), ": ",
                      round(var_exp_k, 1), "%")
)

ggplot(df_scree, aes(x = Componente, y = Varianza)) +
  geom_point(color = "steelblue", size = 3) +
  geom_line(color = "steelblue") +
  geom_text_repel(data = df_scree[1:4, ],
                  aes(label = Etiqueta),
                  size = 4, color = "darkred", fontface = "bold",
                  box.padding = 0.4) +
  labs(title = "Scree plot — varianza explicada por PCA (Kallisto)",
       x = "Componente principal",
       y = "% de varianza explicada") +
  theme_minimal(base_size = 12) +
  theme(plot.background = element_rect(fill = "white"))
Figura 3. Scree plot: porcentaje de varianza explicada por cada componente principal (Kallisto). Las etiquetas muestran los primeros cuatro componentes.

Figura 3. Scree plot: porcentaje de varianza explicada por cada componente principal (Kallisto). Las etiquetas muestran los primeros cuatro componentes.

Interpretación: PC1 y PC2 capturan conjuntamente el 68.8% de la varianza total, lo que indica que la representación bidimensional del PCA es una descripción  fiel de la estructura global de los datos. La caída pronunciada entre  PC1 (43.9%) y los componentes siguientes confirma que la mayor parte de  la variación está concentrada en el primer eje, consistente con una diferencia transcriptómica dominante entre genotipos.

3.3 PCA y STAR

mat_s <- assay(vsd_star)
pc_s  <- prcomp(t(mat_s))

var_exp_s <- pc_s$sdev^2 / sum(pc_s$sdev^2) * 100

pca_df_s <- data.frame(
  PC1      = pc_s$x[, 1],
  PC2      = pc_s$x[, 2],
  genotype = metadata$genotype,
  sample   = rownames(metadata)
)

ggplot(pca_df_s, aes(x = PC1, y = PC2, color = genotype, label = sample)) +
  geom_point(size = 4) +
  geom_text_repel(size = 3, show.legend = FALSE) +
  stat_ellipse(aes(group = genotype), linetype = "dashed", level = 0.8) +
  scale_color_manual(values = c("WT" = "steelblue", "Dp16" = "firebrick3")) +
  labs(title = "PCA — STAR (VST)",
       x     = sprintf("PC1 (%.1f%% varianza)", var_exp_s[1]),
       y     = sprintf("PC2 (%.1f%% varianza)", var_exp_s[2]),
       color = "Genotipo") +
  theme_minimal(base_size = 13) +
  theme(plot.background = element_rect(fill = "white"))
Figura 4. PCA de muestras normalizadas con VST (STAR). Se esperan los mismos patrones de agrupamiento que en Kallisto si ambos alineadores son concordantes.

Figura 4. PCA de muestras normalizadas con VST (STAR). Se esperan los mismos patrones de agrupamiento que en Kallisto si ambos alineadores son concordantes.

Interpretación: El PCA de STAR reproduce el mismo patrón que Kallisto: las tres réplicas WT se agrupan a valores negativos de PC1 (44.6%) y las tres réplicas Dp16 a valores positivos. La concordancia entre ambos pipelines confirma que la separación entre genotipos no es un artefacto del método de cuantificación sino una señal biológica reproducible, independientemente de si se usa pseudoalineamiento (Kallisto) o alineamiento directo al genoma (STAR).

3.4 Corrección de batch effect

Aunque no tenemos un batch effect conocido en este experimento, aplicamos limma::removeBatchEffect como buena práctica para visualización. El argumento design protege la señal biológica del genotipo durante la corrección: la función solo elimina la variación que no está explicada por el genotipo.

Nota importante: Esta corrección es únicamente para visualización (PCA y heatmaps). El análisis diferencial con DESeq2 siempre se corre sobre los datos originales sin corregir; si hubiera un batch real, se incluiría como covariable en la fórmula de diseño (~ batch + genotype).

vsd_bc    <- vsd_kallisto   # copia para no modificar el original
mm_design <- model.matrix(~ genotype, colData(vsd_bc))
mat_bc    <- assay(vsd_bc)
mat_bc    <- removeBatchEffect(mat_bc, design = mm_design)
assay(vsd_bc) <- mat_bc

pcaData_bc    <- plotPCA(vsd_bc, intgroup = "genotype", returnData = TRUE)
percentVar_bc <- round(100 * attr(pcaData_bc, "percentVar"))

ggplot(pcaData_bc, aes(PC1, PC2, color = genotype, label = name)) +
  geom_point(size = 4) +
  geom_text_repel(size = 3, show.legend = FALSE) +
  scale_color_manual(values = c("WT" = "steelblue", "Dp16" = "firebrick3")) +
  labs(title = "PCA post-corrección de batch (Kallisto VST)",
       x = paste0("PC1: ", percentVar_bc[1], "% varianza"),
       y = paste0("PC2: ", percentVar_bc[2], "% varianza"),
       color = "Genotipo") +
  theme_minimal(base_size = 13) +
  theme(plot.background = element_rect(fill = "white"))
Figura 5. PCA post-corrección de batch effect (Kallisto VST). Sin batch conocido, la corrección elimina solo la variación residual no explicada por el genotipo, lo que generalmente refuerza la separación entre grupos.

Figura 5. PCA post-corrección de batch effect (Kallisto VST). Sin batch conocido, la corrección elimina solo la variación residual no explicada por el genotipo, lo que generalmente refuerza la separación entre grupos.

Interpretación: Tras la corrección, PC1 aumenta de 43.9% a 67%, lo que indica que la variación asociada al genotipo se concentra aún más en el primer componente al eliminar el ruido técnico residual. La separación entre WT y Dp16 se mantiene clara, confirmando que la señal biológica es sólida y no dependía de efectos técnicos. SRR33445254 continúa ligeramente elevado en PC2, sugiriendo variabilidad biológica propia de esa réplica que no corresponde a un artefacto corregible.

4. Análisis de Expresión Diferencial (DESeq2)

DESeq2 modela las cuentas con una distribución binomial negativa que captura la sobredispersión típica de RNA-seq. El análisis realiza tres pasos: (1) estimación de factores de tamaño para corregir diferencias en el total de lecturas, (2) estimación de la dispersión gen a gen, y (3) ajuste del modelo lineal generalizado y prueba de Wald para obtener log2FC y p-values.

4.1 DESeq2 y Kallisto

# relevel: WT como referencia para que todos los FC se interpreten como "Dp16 vs WT"
dds$genotype <- relevel(dds$genotype, ref = "WT")
dds <- DESeq(dds)
resultsNames(dds)
## [1] "Intercept"           "genotype_Dp16_vs_WT"
res_kallisto <- results(dds,
                        contrast = c("genotype", "Dp16", "WT"),
                        alpha    = 0.05)
summary(res_kallisto)
## 
## out of 14715 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 171, 1.2%
## LFC < 0 (down)     : 229, 1.6%
## outliers [1]       : 17, 0.12%
## low counts [2]     : 0, 0%
## (mean count < 5)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

4.2 Shrinkage de log2FC y Kallisto

Genes con baja expresión o pocas réplicas pueden presentar log2FC muy grandes e inestables, no porque realmente cambien mucho, sino porque hay poca información para estimarlos. El shrinkage con apeglm aplica un prior adaptativo que “contrae” esos valores extremos hacia 0, sin modificar los p-values ni el número de genes significativos. El resultado son estimaciones de FC más confiables biológicamente y visualizaciones más limpias.

# coef debe coincidir exactamente con un nombre de resultsNames(dds)
res_kallisto_shrink <- lfcShrink(dds,
                                 coef = "genotype_Dp16_vs_WT",
                                 type = "apeglm")
summary(res_kallisto_shrink)
## 
## out of 14715 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 240, 1.6%
## LFC < 0 (down)     : 297, 2%
## outliers [1]       : 17, 0.12%
## low counts [2]     : 0, 0%
## (mean count < 5)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

4.3 DESeq2 y STAR

dds_star$genotype <- relevel(dds_star$genotype, ref = "WT")
dds_star <- DESeq(dds_star)
resultsNames(dds_star)
## [1] "Intercept"           "genotype_Dp16_vs_WT"
res_star <- results(dds_star,
                    contrast = c("genotype", "Dp16", "WT"),
                    alpha    = 0.05)
summary(res_star)
## 
## out of 14091 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 152, 1.1%
## LFC < 0 (down)     : 211, 1.5%
## outliers [1]       : 10, 0.071%
## low counts [2]     : 0, 0%
## (mean count < 5)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
res_star_shrink <- lfcShrink(dds_star,
                              coef = "genotype_Dp16_vs_WT",
                              type = "apeglm")
summary(res_star_shrink)
## 
## out of 14091 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 237, 1.7%
## LFC < 0 (down)     : 289, 2.1%
## outliers [1]       : 10, 0.071%
## low counts [2]     : 0, 0%
## (mean count < 5)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

4.4 Tabla comparativa de DEGs

sig_kallisto <- as.data.frame(res_kallisto_shrink) %>%
  filter(!is.na(padj) & padj < 0.05)

sig_star <- as.data.frame(res_star_shrink) %>%
  filter(!is.na(padj) & padj < 0.05)

resumen <- data.frame(
  Alineador      = c("Kallisto", "STAR"),
  Genes_totales  = c(nrow(dds), nrow(dds_star)),
  DEGs_total     = c(nrow(sig_kallisto), nrow(sig_star)),
  Up_regulated   = c(sum(sig_kallisto$log2FoldChange > 0),
                     sum(sig_star$log2FoldChange > 0)),
  Down_regulated = c(sum(sig_kallisto$log2FoldChange < 0),
                     sum(sig_star$log2FoldChange < 0))
)

kable(resumen,
      col.names = c("Alineador", "Genes totales", "DEGs (padj<0.05)",
                    "Up-regulated", "Down-regulated"),
      caption = "Tabla 1. Resumen de genes diferencialmente expresados por alineador (Dp16 vs WT, padj < 0.05).") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = FALSE)
Tabla 1. Resumen de genes diferencialmente expresados por alineador (Dp16 vs WT, padj < 0.05).
Alineador Genes totales DEGs (padj<0.05) Up-regulated Down-regulated
Kallisto 14715 400 171 229
STAR 14091 363 152 211

Interpretación: El análisis diferencial identificó 400 DEGs en Kallisto y 363 en STAR (padj < 0.05), con un balance consistente entre pipelines: aproximadamente 57% de los genes están reprimidos en Dp16 respecto a WT y 43% activados. Este patrón sugiere que la trisomía del cromosoma 16 no simplemente incrementa la expresión génica global, sino que genera una reorganización transcriptómica donde la represión domina ligeramente sobre la activación en hígado adulto.

La concordancia entre Kallisto y STAR (~90% en el total de DEGs y balance up/down idéntico) indica que este patrón es robusto e independiente del método de cuantificación. La diferencia en magnitud respecto a los ~1,500 DEGs reportados por Dunn et al. (2026) es atribuible al menor poder estadístico de este análisis (n=3 por genotipo vs n=6 en el paper original); los genes detectados aquí representan los cambios de expresión más robustos y de mayor magnitud del dataset.

5. Top genes más significativos

Top 20 DEGs y Kallisto

top20_k <- sig_kallisto %>%
  arrange(padj) %>%
  head(20) %>%
  select(baseMean, log2FoldChange, lfcSE, pvalue, padj) %>%
  mutate(across(where(is.numeric), ~ round(.x, 4)))

kable(top20_k,
      col.names = c("Expr. media", "log2FC", "SE", "p-value", "p-adj"),
      caption = "Tabla 2. Top 20 genes más significativos — Kallisto (Dp16 vs WT, padj < 0.05). Los genes están ordenados por p-value ajustado de menor a mayor.") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = FALSE)
Tabla 2. Top 20 genes más significativos — Kallisto (Dp16 vs WT, padj < 0.05). Los genes están ordenados por p-value ajustado de menor a mayor.
Expr. media log2FC SE p-value p-adj
ENSMUSG00000097148.2 221.5899 -3.8684 0.4515 0 0
ENSMUSG00000001665.11 2301.2742 1.6308 0.1995 0 0
ENSMUSG00000117780.1 615.5616 -3.1112 0.3911 0 0
ENSMUSG00000032080.6 13254.4287 -1.2916 0.1689 0 0
ENSMUSG00000029417.9 817.9000 -1.8694 0.2546 0 0
ENSMUSG00000006494.11 3493.7405 1.4713 0.2022 0 0
ENSMUSG00000069456.4 3868.0862 1.3199 0.1866 0 0
ENSMUSG00000094708.1 105.7300 10.9872 3.0280 0 0
ENSMUSG00000095937.1 105.7300 10.9872 3.0280 0 0
ENSMUSG00000029822.15 127.7477 -2.3736 0.3493 0 0
ENSMUSG00000057074.6 6673.3986 1.0001 0.1484 0 0
ENSMUSG00000078606.8 488.7620 -1.4457 0.2314 0 0
ENSMUSG00000027596.10 39.5124 -9.6737 2.9575 0 0
ENSMUSG00000037827.5 45.5567 9.5875 2.8760 0 0
ENSMUSG00000095597.2 40.6972 9.4140 2.8425 0 0
ENSMUSG00000024899.7 9075.8872 0.8624 0.1441 0 0
ENSMUSG00000078922.9 427.4633 -2.1887 0.3718 0 0
ENSMUSG00000074240.5 144.8147 -4.3150 0.7433 0 0
ENSMUSG00000090381.2 41.3654 -5.4125 0.9232 0 0
ENSMUSG00000063590.3 94.0080 -4.3512 0.7631 0 0

Top 20 DEGs y STAR

top20_s <- sig_star %>%
  arrange(padj) %>%
  head(20) %>%
  select(baseMean, log2FoldChange, lfcSE, pvalue, padj) %>%
  mutate(across(where(is.numeric), ~ round(.x, 4)))

kable(top20_s,
      col.names = c("Expr. media", "log2FC", "SE", "p-value", "p-adj"),
      caption = "Tabla 3. Top 20 genes más significativos — STAR (Dp16 vs WT, padj < 0.05). Los genes están ordenados por p-value ajustado de menor a mayor.") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = FALSE)
Tabla 3. Top 20 genes más significativos — STAR (Dp16 vs WT, padj < 0.05). Los genes están ordenados por p-value ajustado de menor a mayor.
Expr. media log2FC SE p-value p-adj
ENSMUSG00000117780.1 624.4800 -3.3229 0.3469 0 0
ENSMUSG00000029417.9 772.7797 -1.8147 0.1940 0 0
ENSMUSG00000097148.2 186.4670 -3.8658 0.4495 0 0
ENSMUSG00000001665.11 1001.2573 1.5612 0.1837 0 0
ENSMUSG00000069456.4 3676.6815 1.3136 0.1731 0 0
ENSMUSG00000057074.6 6515.7780 1.0052 0.1393 0 0
ENSMUSG00000032080.6 2017.0016 -1.2743 0.1793 0 0
ENSMUSG00000024411.11 72.6810 -2.6860 0.4011 0 0
ENSMUSG00000021416.11 57.1663 5.8161 0.8390 0 0
ENSMUSG00000015090.13 61.9050 3.7672 0.5926 0 0
ENSMUSG00000040938.16 392.8636 1.3553 0.2179 0 0
ENSMUSG00000056131.13 325.1514 -1.3199 0.2152 0 0
ENSMUSG00000023057.5 858.5779 -0.9842 0.1614 0 0
ENSMUSG00000020926.16 202.3003 2.4592 0.4055 0 0
ENSMUSG00000028359.4 228.2812 2.3011 0.3817 0 0
ENSMUSG00000034220.7 244.6405 -2.5862 0.4312 0 0
ENSMUSG00000074240.5 149.4177 -4.2809 0.7191 0 0
ENSMUSG00000024899.7 9192.1437 0.8076 0.1360 0 0
ENSMUSG00000078922.9 151.4595 -1.9564 0.3323 0 0
ENSMUSG00000027596.10 28.4416 -9.1148 2.8936 0 0

Interpretación de las tablas: baseMean es la expresión promedio del gen en todas las muestras (a mayor baseMean, más abundante es el transcrito). log2FoldChange indica la dirección y magnitud del cambio en Dp16 respecto a WT: un valor positivo significa que el gen está más expresado en Dp16, negativo que está menos expresado. lfcSE es el error estándar del FC después del shrinkage. El padj es el p-value corregido por el método de Benjamini-Hochberg para controlar la tasa de falsos positivos.

6. Visualización de resultados

6.1 Volcano Plot y Kallisto

El volcano plot enfrenta la magnitud del cambio (log2FC, eje X) contra la significancia estadística (-log10(padj), eje Y). Los genes biológicamente más relevantes se encuentran en las esquinas superiores: alto FC y muy significativos. Las líneas discontinuas marcan los umbrales de corte: |log2FC| > 1 (cambio de al menos 2 veces) y padj < 0.05.

df_k <- as.data.frame(res_kallisto_shrink) %>%
  mutate(
    gene       = rownames(.),
    Expression = case_when(
      log2FoldChange >= 1  & padj < 0.05 ~ "Up-regulated",
      log2FoldChange <= -1 & padj < 0.05 ~ "Down-regulated",
      TRUE ~ "Unchanged"
    )
  )

top_labels_k <- df_k %>%
  filter(Expression != "Unchanged") %>%
  arrange(padj) %>%
  head(10)

ggplot(df_k, aes(x = log2FoldChange, y = -log10(padj))) +
  geom_point(aes(color = Expression), size = 0.8, alpha = 0.7) +
  geom_text_repel(data = top_labels_k, aes(label = gene),
                  size = 2.8, max.overlaps = 15, show.legend = FALSE) +
  scale_color_manual(values = c("Up-regulated"   = "firebrick3",
                                "Down-regulated" = "dodgerblue3",
                                "Unchanged"      = "gray70")) +
  geom_vline(xintercept = c(-1, 1), linetype = "dashed", alpha = 0.5) +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", alpha = 0.5) +
  labs(title = "Volcano Plot — Kallisto (Dp16 vs WT)",
       x     = expression("log"[2]*"FC (shrinkage apeglm)"),
       y     = expression("-log"[10]*"padj"),
       color = "Expresión") +
  theme_bw(base_size = 13) +
  theme(plot.background = element_rect(fill = "white"))
Figura 6. Volcano plot — Kallisto (Dp16 vs WT). Rojo: genes up-regulated en Dp16 (log2FC > 1, padj < 0.05). Azul: genes down-regulated (log2FC < -1, padj < 0.05). Gris: sin cambio significativo. Las etiquetas corresponden a los 10 genes con menor padj.

Figura 6. Volcano plot — Kallisto (Dp16 vs WT). Rojo: genes up-regulated en Dp16 (log2FC > 1, padj < 0.05). Azul: genes down-regulated (log2FC < -1, padj < 0.05). Gris: sin cambio significativo. Las etiquetas corresponden a los 10 genes con menor padj.

Interpretación: Los genes más significativos y con mayor cambio se concentran en ambas esquinas superiores, confirmando una respuesta transcriptómica clara en Dp16. Los genes etiquetados con mayor significancia en down-regulated (ENSMUSG00000097148.2, ENSMUSG00000117780.1) tienen log2FC entre -3 y -4, indicando una reducción de expresión de 8 a 16 veces respecto a WT. En up-regulated, ENSMUSG00000001665.11 destaca con alta significancia y log2FC ~1.6.El shrinkage es evidente: la nube gris está comprimida cerca de log2FC = 0, y los genes de baja expresión no dominan la gráfica.

6.2 Volcano Plot y STAR

df_s <- as.data.frame(res_star_shrink) %>%
  mutate(
    gene       = rownames(.),
    Expression = case_when(
      log2FoldChange >= 1  & padj < 0.05 ~ "Up-regulated",
      log2FoldChange <= -1 & padj < 0.05 ~ "Down-regulated",
      TRUE ~ "Unchanged"
    )
  )

top_labels_s <- df_s %>%
  filter(Expression != "Unchanged") %>%
  arrange(padj) %>%
  head(10)

ggplot(df_s, aes(x = log2FoldChange, y = -log10(padj))) +
  geom_point(aes(color = Expression), size = 0.8, alpha = 0.7) +
  geom_text_repel(data = top_labels_s, aes(label = gene),
                  size = 2.8, max.overlaps = 15, show.legend = FALSE) +
  scale_color_manual(values = c("Up-regulated"   = "firebrick3",
                                "Down-regulated" = "dodgerblue3",
                                "Unchanged"      = "gray70")) +
  geom_vline(xintercept = c(-1, 1), linetype = "dashed", alpha = 0.5) +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", alpha = 0.5) +
  labs(title = "Volcano Plot — STAR (Dp16 vs WT)",
       x     = expression("log"[2]*"FC (shrinkage apeglm)"),
       y     = expression("-log"[10]*"padj"),
       color = "Expresión") +
  theme_bw(base_size = 13) +
  theme(plot.background = element_rect(fill = "white"))
Figura 7. Volcano plot — STAR (Dp16 vs WT). Mismos criterios de corte que Kallisto (|log2FC| > 1, padj < 0.05). La comparación visual entre ambos volcano plots permite evaluar la concordancia entre alineadores.

Figura 7. Volcano plot — STAR (Dp16 vs WT). Mismos criterios de corte que Kallisto (|log2FC| > 1, padj < 0.05). La comparación visual entre ambos volcano plots permite evaluar la concordancia entre alineadores.

Interpretación: El patrón es concordante con Kallisto, los mismos genes (ENSMUSG00000117780.1, ENSMUSG00000029417.9, ENSMUSG00000097148.2) aparecen entre los más significativos en ambos alineadores, lo que confirma que son señales biológicas reales y no artefactos de un pipeline específico.

6.3 MA Plot y Kallisto

El MA plot muestra la relación entre la expresión media (baseMean, eje X en log2) y el log2FC (eje Y). Es una herramienta diagnóstica: en datos bien normalizados, los puntos grises (no significativos) deben estar centrados en log2FC = 0 a lo largo de todo el rango de expresión. El efecto del shrinkage es claramente visible: los genes con baja expresión media (izquierda del eje X) tienen FC contraídos hacia 0, mientras que los genes con alta expresión mantienen sus FC originales.

ma_k <- as.data.frame(res_kallisto_shrink) %>%
  mutate(
    significant  = !is.na(padj) & padj < 0.05,
    baseMean_log = log2(baseMean + 1)
  )

ggplot(ma_k, aes(x = baseMean_log, y = log2FoldChange, color = significant)) +
  geom_point(size = 0.6, alpha = 0.6) +
  scale_color_manual(values = c("FALSE" = "gray70", "TRUE" = "steelblue"),
                     labels = c("No significativo", "padj < 0.05")) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  labs(title = "MA Plot — Kallisto (shrinkage apeglm)",
       x     = expression("log"[2]*"(expresión media + 1)"),
       y     = expression("log"[2]*"FC"),
       color = "") +
  theme_bw(base_size = 13) +
  theme(plot.background = element_rect(fill = "white"))
Figura 8. MA plot — Kallisto (shrinkage apeglm). Azul: genes significativos (padj < 0.05). Gris: no significativos. La línea horizontal en log2FC = 0 sirve de referencia: puntos equilibrados sobre y bajo la línea indican ausencia de sesgo sistemático.

Figura 8. MA plot — Kallisto (shrinkage apeglm). Azul: genes significativos (padj < 0.05). Gris: no significativos. La línea horizontal en log2FC = 0 sirve de referencia: puntos equilibrados sobre y bajo la línea indican ausencia de sesgo sistemático.

Interpretación:Los puntos grises están centrados en log2FC = 0 a lo largo de todo el rango de expresión, confirmando que la normalización no introdujo sesgos sistemáticos. Los DEGs significativos (azul) se distribuyen tanto por encima como por debajo de la línea de referencia, consistente con el balance observado de ~57% down y ~43% up. La compresión de puntos hacia 0 en el extremo izquierdo (genes poco expresados) es el efecto esperado del shrinkage apeglm.

6.4 MA Plot y STAR

ma_s <- as.data.frame(res_star_shrink) %>%
  mutate(
    significant  = !is.na(padj) & padj < 0.05,
    baseMean_log = log2(baseMean + 1)
  )

ggplot(ma_s, aes(x = baseMean_log, y = log2FoldChange, color = significant)) +
  geom_point(size = 0.6, alpha = 0.6) +
  scale_color_manual(values = c("FALSE" = "gray70", "TRUE" = "steelblue"),
                     labels = c("No significativo", "padj < 0.05")) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  labs(title = "MA Plot — STAR (shrinkage apeglm)",
       x     = expression("log"[2]*"(expresión media + 1)"),
       y     = expression("log"[2]*"FC"),
       color = "") +
  theme_bw(base_size = 13) +
  theme(plot.background = element_rect(fill = "white"))
Figura 9. MA plot — STAR (shrinkage apeglm). Se esperan patrones similares a Kallisto. Diferencias en la densidad de puntos azules reflejan variaciones en el número de DEGs detectados por cada alineador.

Figura 9. MA plot — STAR (shrinkage apeglm). Se esperan patrones similares a Kallisto. Diferencias en la densidad de puntos azules reflejan variaciones en el número de DEGs detectados por cada alineador.

Interpretación:El patrón es prácticamente idéntico al de Kallisto, con DEGs distribuidos simétricamente alrededor de log2FC = 0 y la misma compresión en genes de baja expresión. La densidad ligeramente menor de puntos azules refleja los 37 DEGs menos detectados respecto a Kallisto, atribuible a diferencias en cómo cada alineador maneja lecturas multi-mapeo.

6.5 Heatmap y Top 20 DEGs (Kallisto)

El heatmap visualiza los patrones de expresión de los 20 genes con menor padj. Cada fila es un gen, cada columna es una muestra. Los colores representan la expresión estandarizada por fila (z-score dentro del heatmap): rojo = expresión alta relativa, azul = expresión baja relativa. El clustering jerárquico de filas y columnas se realiza de forma no supervisada — si las muestras del mismo genotipo se agrupan juntas, significa que los genes seleccionados realmente distinguen WT de Dp16.

top_genes_k <- head(order(res_kallisto$padj, na.last = TRUE), 20)
mat_heat_k  <- assay(vsd_kallisto)[top_genes_k, ]

ann_col    <- data.frame(Genotipo = metadata$genotype)
rownames(ann_col) <- rownames(metadata)
ann_colors <- list(Genotipo = c("WT" = "steelblue", "Dp16" = "firebrick3"))

pheatmap(mat_heat_k,
         cluster_rows      = TRUE,
         cluster_cols      = TRUE,
         show_rownames     = TRUE,
         show_colnames     = TRUE,
         annotation_col    = ann_col,
         annotation_colors = ann_colors,
         scale             = "row",
         color             = colorRampPalette(rev(brewer.pal(9, "RdBu")))(100),
         fontsize_row      = 9,
         main              = "Top 20 DEGs — Kallisto (escala por fila)")
Figura 10. Heatmap de los 20 genes con menor padj — Kallisto (VST, escalado por fila/z-score). La barra de color en la parte superior indica el genotipo de cada muestra. El clustering jerárquico es no supervisado.

Figura 10. Heatmap de los 20 genes con menor padj — Kallisto (VST, escalado por fila/z-score). La barra de color en la parte superior indica el genotipo de cada muestra. El clustering jerárquico es no supervisado.

Interpretación: El clustering jerárquico no supervisado agrupa correctamente las tres réplicas Dp16 (izquierda) separadas de las tres réplicas WT (derecha), sin que se le indique el genotipo. Esto confirma que los 20 genes más significativos son suficientes para distinguir ambos genotipos. Se observan dos bloques principales de genes: uno con expresión alta en Dp16 (rojo) que incluye ENSMUSG00000001665.11 y ENSMUSG00000006494.11, y otro con expresión baja en Dp16 (azul) que incluye ENSMUSG00000097148.2 y ENSMUSG00000117780.1.

6.6 Heatmap y Top 20 DEGs (STAR)

top_genes_s <- head(order(res_star$padj, na.last = TRUE), 20)
mat_heat_s  <- assay(vsd_star)[top_genes_s, ]

pheatmap(mat_heat_s,
         cluster_rows      = TRUE,
         cluster_cols      = TRUE,
         show_rownames     = TRUE,
         show_colnames     = TRUE,
         annotation_col    = ann_col,
         annotation_colors = ann_colors,
         scale             = "row",
         color             = colorRampPalette(rev(brewer.pal(9, "RdBu")))(100),
         fontsize_row      = 9,
         main              = "Top 20 DEGs — STAR (escala por fila)")
Figura 11. Heatmap de los 20 genes con menor padj — STAR (VST, escalado por fila). La comparación con el heatmap de Kallisto (Figura 10) permite identificar genes consistentemente significativos en ambos alineadores.

Figura 11. Heatmap de los 20 genes con menor padj — STAR (VST, escalado por fila). La comparación con el heatmap de Kallisto (Figura 10) permite identificar genes consistentemente significativos en ambos alineadores.

Interpretación: El patrón de clustering es idéntico — Dp16 a la izquierda, WT a la derecha — y varios genes se repiten entre ambos top 20 (ENSMUSG00000001665.11, ENSMUSG00000069456.4, ENSMUSG00000029417.9, ENSMUSG00000074240.5), lo que los identifica como los candidatos más robustos del análisis.

7. Comparación Kallisto vs STAR

7.1 Genes compartidos y exclusivos

genes_kallisto <- rownames(sig_kallisto)
genes_star     <- rownames(sig_star)

genes_comunes <- intersect(genes_kallisto, genes_star)
genes_solo_k  <- setdiff(genes_kallisto, genes_star)
genes_solo_s  <- setdiff(genes_star, genes_kallisto)

venn.plot <- draw.pairwise.venn(
  area1      = length(genes_kallisto),
  area2      = length(genes_star),
  cross.area = length(genes_comunes),
  category   = c("Kallisto", "STAR"),
  fill       = c("steelblue", "firebrick3"),
  alpha      = 0.4,
  col        = "white",
  cex        = 1.5,
  cat.cex    = 1.3,
  cat.pos    = c(-20, 20),
  cat.dist   = 0.07
)
grid.draw(venn.plot)
Figura 12. Diagrama de Venn: DEGs compartidos entre Kallisto y STAR (padj < 0.05). El área de intersección representa genes detectados como diferencialmente expresados por ambos alineadores.

Figura 12. Diagrama de Venn: DEGs compartidos entre Kallisto y STAR (padj < 0.05). El área de intersección representa genes detectados como diferencialmente expresados por ambos alineadores. 

venn_resumen <- data.frame(
  Categoría = c("Solo Kallisto", "Compartidos", "Solo STAR"),
  N_genes   = c(length(genes_solo_k), length(genes_comunes), length(genes_solo_s))
)

kable(venn_resumen,
      caption = "Tabla 4. Distribución de DEGs entre alineadores (padj < 0.05).") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = FALSE)
Tabla 4. Distribución de DEGs entre alineadores (padj < 0.05).
Categoría N_genes
Solo Kallisto 136
Compartidos 264
Solo STAR 99

Interpretación: De los 400 DEGs de Kallisto y 363 de STAR, 264 son compartidos (~66% del total). Esto indica una concordancia sólida entre pipelines. Los 136 exclusivos de Kallisto y 99 de STAR probablemente corresponden a genes con múltiples isoformas o regiones de mapeo ambiguo donde pseudoalineamiento y alineamiento directo difieren en la asignación de lecturas.

8. Análisis Funcional y GO terms con gprofiler2

El análisis funcional convierte una lista de genes en significado biológico. gprofiler2 implementa una prueba hipergeométrica que evalúa si la proporción de genes de una categoría GO/KEGG/REAC en nuestra lista de DEGs es mayor de lo esperado por azar dado el universo de todos los genes analizados. Usamos los genes compartidos entre Kallisto y STAR para maximizar la confiabilidad biológica.

8.1 Preparar listas de genes

# Genes compartidos entre alineadores, clasificados por dirección del cambio
up_genes   <- genes_comunes[sig_kallisto[genes_comunes, "log2FoldChange"] > 0]
down_genes <- genes_comunes[sig_kallisto[genes_comunes, "log2FoldChange"] < 0]

up_genes   <- up_genes[!is.na(up_genes)]
down_genes <- down_genes[!is.na(down_genes)]

# Remover número de versión de IDs Ensembl (ENSMUSG00000001234.5 → ENSMUSG00000001234)
up_genes   <- gsub("\\..*", "", up_genes)
down_genes <- gsub("\\..*", "", down_genes)

cat("Genes up-regulated (compartidos):", length(up_genes), "\n")
## Genes up-regulated (compartidos): 113
cat("Genes down-regulated (compartidos):", length(down_genes), "\n")
## Genes down-regulated (compartidos): 151
head(up_genes)
## [1] "ENSMUSG00000001665" "ENSMUSG00000002588" "ENSMUSG00000002847"
## [4] "ENSMUSG00000003477" "ENSMUSG00000005373" "ENSMUSG00000005677"
head(down_genes)
## [1] "ENSMUSG00000000159" "ENSMUSG00000000303" "ENSMUSG00000000915"
## [4] "ENSMUSG00000001506" "ENSMUSG00000001763" "ENSMUSG00000002033"

Interpretación: De los 264 genes compartidos entre alineadores, 113 están sobreexpresados en Dp16 respecto a WT y 151 están reprimidos, manteniendo el balance observado en el análisis completo (~57% down, ~43% up). Los identificadores Ensembl fueron convertidos a formato sin versión (eliminando el sufijo numérico después del punto) para garantizar la compatibilidad con la base de datos de gprofiler2 para Mus musculus. Estas dos listas (genes up-regulated y down-regulated) se analizaron de forma independiente para identificar qué procesos biológicos se activan y cuáles se reprimen específicamente en el hígado del modelo Dp16.

8.2 Enriquecimiento con gprofiler2

sources_db <- c("GO:BP", "GO:MF", "GO:CC", "KEGG", "REAC")

# multi_query = FALSE: devuelve una fila por término por query,
# con columnas normales: query, p_value, intersection_size, source, term_name.
# Con multi_query = TRUE esas columnas son listas y causan errores al filtrar.
multi_gp <- gost(
  query = list(
    "Up-regulated"   = up_genes,
    "Down-regulated" = down_genes
  ),
  organism          = "mmusculus",
  correction_method = "fdr",
  user_threshold    = 0.05,
  multi_query       = FALSE,
  ordered_query     = FALSE,
  sources           = sources_db,
  evcodes           = TRUE
)

cat("Columnas del resultado:", paste(colnames(multi_gp$result), collapse = ", "), "\n")
## Columnas del resultado: query, significant, p_value, term_size, query_size, intersection_size, precision, recall, term_id, source, term_name, effective_domain_size, source_order, parents, evidence_codes, intersection
cat("Queries:", paste(unique(multi_gp$result$query), collapse = ", "), "\n")
## Queries: Down-regulated, Up-regulated
cat("Términos enriquecidos totales:", nrow(multi_gp$result), "\n")
## Términos enriquecidos totales: 720

Interpretación: El análisis con gprofiler2 identificó 720 términos enriquecidos en total entre genes up-regulated y down-regulated (FDR < 0.05), distribuidos en cinco bases de datos: GO:BP, GO:MF, GO:CC, KEGG y REAC. El alto número de términos significativos refleja que los 264 genes compartidos, aunque representan una fracción pequeña del transcriptoma analizado, convergen en procesos biológicos bien definidos y estadísticamente sobre-representados en ambas direcciones de cambio.

8.3 Manhattan plot

El manhattan plot de gprofiler2 muestra todos los términos enriquecidos organizados por base de datos (color) y ordenados en el eje X. El eje Y es el -log10(FDR): puntos más altos = más significativos. Es útil para tener una visión panorámica de qué bases de datos están más enriquecidas y dónde se concentran los términos significativos.

# gostplot() es la función nativa de gprofiler2 para este tipo de gráfico.
# Es la más robusta porque maneja internamente la estructura del objeto multi_gp.
gostplot(multi_gp, interactive = FALSE)
Figura 14. Manhattan plot de términos GO (BP, MF, CC), KEGG y REAC enriquecidos en los genes DEG compartidos (FDR < 0.05). Cada color representa una base de datos diferente.

Figura 14. Manhattan plot de términos GO (BP, MF, CC), KEGG y REAC enriquecidos en los genes DEG compartidos (FDR < 0.05). Cada color representa una base de datos diferente.

Interpretación: En los genes up-regulated, GO:MF presenta el término individual más significativo de todo el análisis (-log10 ~15), mientras que GO:BP concentra el mayor número de términos enriquecidos, varios de ellos superando el límite de >16. KEGG y REAC también muestran señal notable, indicando que las rutas metabólicas activadas están anotadas en múltiples bases de datos. En los genes down-regulated, GO:CC es la categoría más significativa con un término que alcanza -log10 ~10, mientras que GO:BP tiene muchos términos pero de significancia moderada (~2-4). KEGG y REAC tienen señal mínima en down-regulated, lo que indica que la represión inmune no converge en rutas metabólicas conocidas tan claramente como la activación.

8.4 Top términos enriquecidos y Up-regulated

if (!is.null(multi_gp$result)) {
  gost_df <- as.data.frame(multi_gp$result)

  bar_up <- gost_df %>%
    filter(query == "Up-regulated") %>%
    arrange(p_value) %>%
    head(20) %>%
    select(term_name, source, term_size, intersection_size, p_value) %>%
    mutate(p_value = formatC(p_value, format = "e", digits = 2))

  kable(bar_up,
        col.names = c("Término GO", "Base de datos", "Tamaño del término",
                      "Genes en overlap", "p-value (FDR)"),
        caption = "Tabla 5. Top 20 términos enriquecidos en genes up-regulated en Dp16 vs WT. 'Tamaño del término' indica cuántos genes del genoma pertenecen a esa categoría; 'Genes en overlap' es cuántos de nuestros DEGs caen en ella.") %>%
    kable_styling(bootstrap_options = c("striped","hover","condensed"),
                  full_width = FALSE)
}
Tabla 5. Top 20 términos enriquecidos en genes up-regulated en Dp16 vs WT. ‘Tamaño del término’ indica cuántos genes del genoma pertenecen a esa categoría; ‘Genes en overlap’ es cuántos de nuestros DEGs caen en ella.
Término GO Base de datos Tamaño del término Genes en overlap p-value (FDR)
small molecule metabolic process GO:BP 1549 44 1.54e-23
lipid metabolic process GO:BP 1266 40 7.01e-23
oxoacid metabolic process GO:BP 783 32 5.42e-21
organic acid metabolic process GO:BP 789 32 5.42e-21
carboxylic acid metabolic process GO:BP 769 30 4.36e-19
monocarboxylic acid metabolic process GO:BP 549 26 2.46e-18
catalytic activity GO:MF 5584 68 1.16e-15
fatty acid metabolic process GO:BP 357 18 1.33e-12
small molecule biosynthetic process GO:BP 500 20 2.43e-12
Metabolic pathways KEGG 1646 40 1.00e-11
organic acid biosynthetic process GO:BP 251 15 2.01e-11
steroid metabolic process GO:BP 277 15 7.61e-11
alcohol metabolic process GO:BP 287 15 1.16e-10
carboxylic acid biosynthetic process GO:BP 247 14 2.21e-10
Drug ADME REAC 105 13 2.79e-10
oxidoreductase activity GO:MF 742 22 5.97e-10
lipid biosynthetic process GO:BP 610 19 6.20e-10
sulfur compound metabolic process GO:BP 224 13 9.03e-10
monocarboxylic acid biosynthetic process GO:BP 176 12 9.03e-10
Metabolism REAC 1740 39 1.61e-09

Interpretación: Los genes activados en Dp16 se enriquecen de forma masiva en procesos metabólicos hepáticos. Los dos términos más significativos “small molecule metabolic process” y “lipid metabolic process” tienen FDR de 1.54e-23 y 7.01e-23 respectivamente, valores extremadamente por encima del esperado por azar. El patrón es coherente a lo largo de los 20 términos: metabolismo y síntesis de lípidos, ácidos grasos, ácidos orgánicos y esteroides forman un bloque funcional activado en conjunto. El término “catalytic activity” de GO:MF (FDR 1.16e-15) con 68 genes en overlap indica que gran parte de los genes activados codifican enzimas metabólicas. Esto es directamente concordante con lo reportado por Dunn et al. (2026), que identifica disfunción del metabolismo lipídico y de ácidos biliares como el fenotipo hepático principal del modelo Dp16.

8.5 Top términos enriquecidos y Down-regulated

if (!is.null(multi_gp$result)) {
  gost_df <- as.data.frame(multi_gp$result)

  bar_down <- gost_df %>%
    filter(query == "Down-regulated") %>%
    arrange(p_value) %>%
    head(20) %>%
    select(term_name, source, term_size, intersection_size, p_value) %>%
    mutate(p_value = formatC(p_value, format = "e", digits = 2))

  kable(bar_down,
        col.names = c("Término GO", "Base de datos", "Tamaño del término",
                      "Genes en overlap", "p-value (FDR)"),
        caption = "Tabla 6. Top 20 términos enriquecidos en genes down-regulated en Dp16 vs WT. Procesos reprimidos en el modelo de trisomía 16 respecto al ratón silvestre.") %>%
    kable_styling(bootstrap_options = c("striped","hover","condensed"),
                  full_width = FALSE)
}
Tabla 6. Top 20 términos enriquecidos en genes down-regulated en Dp16 vs WT. Procesos reprimidos en el modelo de trisomía 16 respecto al ratón silvestre.
Término GO Base de datos Tamaño del término Genes en overlap p-value (FDR)
extracellular region GO:CC 2340 44 5.80e-11
extracellular space GO:CC 1136 24 2.79e-06
membrane attack complex GO:CC 7 4 3.55e-06
pore complex GO:CC 26 5 2.49e-05
Terminal pathway of complement REAC 8 4 1.14e-04
immune system process GO:BP 2557 33 4.06e-04
defense response GO:BP 1633 25 4.06e-04
response to stress GO:BP 3596 41 4.06e-04
immune response GO:BP 1726 26 4.06e-04
PPAR signaling pathway KEGG 89 7 7.17e-04
Complement and coagulation cascades KEGG 94 7 7.17e-04
endoplasmic reticulum GO:CC 1758 25 8.25e-04
cell surface GO:CC 753 15 1.13e-03
humoral immune response mediated by circulating immunoglobulin GO:BP 43 5 1.32e-03
endomembrane system GO:CC 4005 41 2.22e-03
response to external stimulus GO:BP 2346 29 2.27e-03
response to stimulus GO:BP 9759 74 2.49e-03
positive regulation of response to stimulus GO:BP 2034 26 2.49e-03
response to other organism GO:BP 1566 22 2.49e-03
response to external biotic stimulus GO:BP 1570 22 2.49e-03

Interpretación: Los genes reprimidos en Dp16 muestran dos patrones distintos en la tabla. El término más significativo es “extracellular region” de GO:CC (FDR 5.80e-11, 44 genes en overlap), seguido de “extracellular space” y “membrane attack complex” lo que indica que los genes reprimidos codifican principalmente proteínas secretadas o de la matriz extracelular. El segundo patrón es inmune: “immune system process”, “defense response”, “immune response” y “response to stress” comparten el mismo FDR (4.06e-04), y la presencia de “Terminal pathway of complement” en REAC y “Complement and coagulation cascades” en KEGG apunta específicamente a la represión de la cascada del complemento en el hígado Dp16.

Comparado con los términos up-regulated, los FDR son considerablemente menos extremos (máximo 5.80e-11 vs 1.54e-23), lo que indica que la señal metabólica activada es estadísticamente más robusta que la señal de represión. Esto es consistente con lo reportado por Dunn et al. (2026), que describe disfunción metabólica e inflamación como los fenotipos dominantes del hígado Dp16.

8.6 Barplot GO:BP y Up-regulated

if (!is.null(multi_gp$result)) {
  gost_df <- as.data.frame(multi_gp$result)

  bar_up_plot <- gost_df %>%
    filter(query == "Up-regulated", source == "GO:BP") %>%
    arrange(p_value) %>%
    head(20) %>%
    mutate(log10p = round(-log10(p_value), 2),
           num    = seq_len(n()))

  if (nrow(bar_up_plot) > 0) {
    ggplot(bar_up_plot,
           aes(x = log10p, y = reorder(term_name, log10p))) +
      geom_bar(stat = "identity", fill = "firebrick3", alpha = 0.8) +
      geom_text(aes(label = log10p), hjust = -0.1, size = 3) +
      labs(title = "GO:BP — Up-regulated (Dp16 vs WT)",
           x = "-log10(FDR)", y = NULL) +
      theme_classic(base_size = 12) +
      theme(plot.background = element_rect(fill = "white"))
  }
}
Figura 15. Top 20 procesos biológicos (GO:BP) enriquecidos en genes up-regulated en Dp16 vs WT. El eje X muestra la significancia en escala -log10(FDR): barras más largas = mayor enriquecimiento.

Figura 15. Top 20 procesos biológicos (GO:BP) enriquecidos en genes up-regulated en Dp16 vs WT. El eje X muestra la significancia en escala -log10(FDR): barras más largas = mayor enriquecimiento.

Interpretación: Los 20 procesos GO:BP más enriquecidos en genes up-regulated forman un bloque funcional coherente centrado en el metabolismo hepático. Los seis términos más significativos metabolismo de pequeñas moléculas, lípidos, oxoácidos, ácidos orgánicos, ácidos carboxílicos y ácidos monocarboxílicos todos superan -log10(FDR) de 17, indicando una activación masiva y coordinada de la maquinaria metabólica. La presencia de “fatty acid metabolic process”, “steroid metabolic process” y “acyl-CoA metabolic process” en los últimos puestos confirma que la activación no es de un solo proceso sino de toda la red metabólica hepática. Esto es directamente concordante con lo reportado por Dunn et al. (2026).

8.7 Barplot GO:BP — Down-regulated

if (!is.null(multi_gp$result)) {
  gost_df <- as.data.frame(multi_gp$result)

  bar_down_plot <- gost_df %>%
    filter(query == "Down-regulated", source == "GO:BP") %>%
    arrange(p_value) %>%
    head(20) %>%
    mutate(log10p = round(-log10(p_value), 2),
           num    = seq_len(n()))

  if (nrow(bar_down_plot) > 0) {
    ggplot(bar_down_plot,
           aes(x = log10p, y = reorder(term_name, log10p))) +
      geom_bar(stat = "identity", fill = "dodgerblue3", alpha = 0.8) +
      geom_text(aes(label = log10p), hjust = -0.1, size = 3) +
      labs(title = "GO:BP — Down-regulated (Dp16 vs WT)",
           x = "-log10(FDR)", y = NULL) +
      theme_classic(base_size = 12) +
      theme(plot.background = element_rect(fill = "white"))
  }
}
Figura 16. Top 20 procesos biológicos (GO:BP) enriquecidos en genes down-regulated en Dp16 vs WT. Procesos represados en el modelo de trisomía, ordenados por significancia.

Figura 16. Top 20 procesos biológicos (GO:BP) enriquecidos en genes down-regulated en Dp16 vs WT. Procesos represados en el modelo de trisomía, ordenados por significancia.

Interpretación: Los procesos GO:BP reprimidos en Dp16 se concentran en respuesta inmune y defensa. Los cuatro términos más significativos “response to stress”, “immune system process”, “immune response” y “defense response” comparten exactamente el mismo FDR (-log10 = 3.39), lo que indica que forman un bloque co-reprimido. La presencia de “complement activation” por vía clásica y alternativa, y “humoral immune response mediated by circulating immunoglobulin” apunta específicamente a la represión de la inmunidad humoral en el hígado Dp16. Comparando ambas figuras, la escala es muy diferente: el máximo en down-regulated es 3.39 vs 22.81 en up-regulated, lo que refleja que la activación metabólica es el fenotipo dominante y estadísticamente más robusto en el hígado Dp16.

9. Guardar resultados

if (!dir.exists(figdir)) dir.create(figdir, recursive = TRUE)

# Resultados DEG completos (con shrinkage apeglm)
write.csv(as.data.frame(res_kallisto_shrink),
          file = paste0(outdir, "DESeq2_kallisto_results.csv"))
write.csv(as.data.frame(res_star_shrink),
          file = paste0(outdir, "DESeq2_STAR_results.csv"))

# Solo genes significativos (padj < 0.05)
write.csv(sig_kallisto,
          file = paste0(outdir, "DESeq2_kallisto_significant.csv"))
write.csv(sig_star,
          file = paste0(outdir, "DESeq2_STAR_significant.csv"))

# Resultados GO terms
if (!is.null(multi_gp$result)) {
  gost_export <- multi_gp$result
  # Convertir columnas tipo lista a texto para exportar en CSV
  gost_export$parents        <- sapply(gost_export$parents, paste, collapse = ";")
  gost_export$evidence_codes <- sapply(gost_export$evidence_codes, paste, collapse = ";")
  write.csv(gost_export,
            file      = paste0(outdir, "GOterms_gprofiler2_results.csv"),
            row.names = FALSE)
}

# Objetos R: dds, vsd, resultados, z-score, melted counts
save(dds, dds_star,
     vsd_kallisto, vsd_star,
     res_kallisto_shrink, res_star_shrink,
     zscore_mat, melted_norm_counts,
     file = paste0(outdir, "DEG_objects.RData"))

cat("Todos los archivos guardados en:", outdir, "\n")
## Todos los archivos guardados en: /Users/dubsirubs/Desktop/RNAseq_Equipo4/results_final/

10. Conclusiones

  1. Normalización exitosa: La transformación VST produjo distribuciones de expresión comparables entre todas las muestras (Figura 1). El PCA separó claramente los genotipos WT y Dp16 (Figuras 2 y 4), confirmando que el genotipo es la principal fuente de variación en los datos y que no existe un batch effect evidente. El scree plot (Figura 3) indica qué fracción de la varianza total está capturada en los dos primeros componentes.

  2. Expresión diferencial robusta: Se identificaron 400 DEGs con Kallisto y 363 con STAR (padj < 0.05), con predominio de genes reprimidos (~57% down) en ambos pipelines. El shrinkage con apeglm estabilizó los log2FC de genes con baja expresión, como se observa en los MA plots (Figuras 8 y 9). Los heatmaps (Figuras 10 y 11) muestran que los 20 genes más significativos son suficientes para separar ambos genotipos de forma no supervisada.

  3. Alta concordancia Kallisto vs STAR: El diagrama de Venn (Figura 12) muestra que 264 genes son compartidos entre ambos pipelines (~66% de Kallisto, ~73% de STAR), validando que la firma transcriptómica de Dp16 es reproducible independientemente del método de cuantificación.

  4. Procesos biológicos alterados: El análisis funcional con gprofiler2 identificó 720 términos enriquecidos en los 264 genes compartidos. Los genes up-regulated se enriquecen masivamente en metabolismo de lípidos, ácidos grasos y ácidos orgánicos (FDR hasta 1.54e-23), mientras que los down-regulated se enriquecen en respuesta inmune y activación del complemento (FDR hasta 4.06e-04). Estos resultados son consistentes con la disfunción metabólica e inflamatoria hepática reportada por Dunn et al. (2026) en el modelo Dp16.

11. Discusión

11.1 Comparación con el artículo de referencia (Dunn et al., 2026)

Nuestros resultados son consistentes con los hallazgos centrales reportados por Dunn et al. (2026) en el mismo modelo Dp(16)1Yey y el mismo tejido (hígado de ratón adulto). El artículo de referencia reporta disfunción metabólica hepática como fenotipo dominante en Dp16, caracterizada por alteración del metabolismo de ácidos biliares, fibrosis e inflamación. En nuestro análisis, los genes up-regulated en Dp16 se enriquecen masivamente en términos GO relacionados con metabolismo de lípidos, ácidos grasos y ácidos orgánicos (FDR hasta 1.54e-23), mientras que los down-regulated se concentran en respuesta inmune y cascada del complemento (FDR hasta 5.80e-11). Este patrón de activación metabólica combinada con represión de la inmunidad humoral replica la firma transcriptómica descrita por Dunn et al. (2026).

Una diferencia notable es el número de DEGs detectados: nuestro análisis identificó 400 genes (Kallisto) y 363 (STAR) con padj < 0.05, mientras que Dunn et al. reportan ~1,500 DEGs con el dataset completo. Esta diferencia es atribuible directamente al menor poder estadístico de usar solo 3 réplicas por genotipo en lugar de las 6 disponibles. A pesar de ello, la dirección de los cambios y los procesos biológicos identificados son concordantes, lo que valida la robustez de nuestra aproximación analítica.

11.2 Contexto biológico: metabolismo hepático y trisomía 16

El perfil transcriptómico observado en el hígado Dp16 se puede interpretar en el contexto más amplio de la fisiopatología del Síndrome de Down. Tian et al. (2023) realizaron un análisis metabólico multi-tejido exhaustivo del modelo Dp16, incluyendo hígado, músculo esquelético y tejido adiposo, y demostraron que la trisomía genera resistencia a la insulina, intolerancia a la glucosa y dislipidemia en ambos sexos. Sus datos transcriptómicos multi-tejido muestran activación de vías del metabolismo lipídico y supresión de la función mitocondrial como patrones compartidos. Nuestros resultados en hígado son compatibles con este patrón sistémico: la activación masiva de genes del metabolismo de lípidos y ácidos grasos que observamos en Dp16 refleja esta disfunción metabólica sistémica. La concordancia entre el análisis de un solo tejido (nuestro estudio) y el análisis multi-tejido de Tian et al. (2023) sugiere que el hígado es uno de los principales efectores del fenotipo metabólico de la trisomía 16.

11.3 Métodos de cuantificación: Kallisto vs STAR

La comparación entre pipelines que realizamos tiene respaldo metodológico en el artículo fundacional de DESeq2 (Love et al., 2014), donde se enfatiza que la precisión del análisis diferencial depende críticamente de la calidad de las estimaciones de conteo. Nuestro análisis muestra que 264 de los DEGs son compartidos entre Kallisto y STAR (~66–73% de overlap), lo que indica alta concordancia entre pseudoalineamiento y alineamiento directo al genoma para este dataset. La pequeña fracción de DEGs exclusivos de cada pipeline puede explicarse por las diferencias metodológicas: Kallisto trabaja a nivel de transcriptoma y propaga incertidumbre de asignación mediante bootstrap, mientras que STAR alinea al genoma completo y maneja lecturas multi-mapeo de forma diferente. Para el análisis de expresión diferencial en tejido hepático con genes bien anotados en GRCm39, ambas estrategias son equivalentes en términos prácticos, en línea con lo reportado por Love et al. (2014) sobre la robustez del modelo binomial negativo frente a distintas estrategias de cuantificación.

11.4 Represión de la cascada del complemento en Dp16

Un hallazgo particularmente relevante de nuestro análisis funcional es la represión de genes del complemento y la inmunidad humoral en el hígado Dp16. El hígado es el principal sitio de síntesis de proteínas del complemento, por lo que su represión en el contexto de Dp16 podría contribuir a la vulnerabilidad inmunológica observada clínicamente en personas con Síndrome de Down. Bray et al. (2016) demostraron que el pseudoalineamiento de Kallisto preserva la sensibilidad para detectar cambios en genes con expresión moderada-baja, que es precisamente el rango donde se ubican varios genes del complemento. Esto refuerza la confiabilidad de este hallazgo en nuestro análisis. El patrón de represión del complemento identificado aquí es consistente con lo reportado por Dunn et al. (2026) en sus análisis de enriquecimiento de vías, y sugiere que la disfunción hepática en trisomía 21 tiene consecuencias que van más allá del metabolismo y se extienden a la función inmune sistémica.

11.5 Limitaciones del análisis y perspectivas futuras

El principal factor limitante de este estudio es el tamaño de muestra (n=3 por genotipo), seleccionado por criterios de profundidad de secuenciación y porcentaje de duplicación a partir del dataset completo de 6 réplicas. Aunque el diseño estadístico de DESeq2 maneja adecuadamente réplicas triples cuando la señal biológica es fuerte (Love et al., 2014), el poder para detectar genes con cambios moderados (|log2FC| < 1) es reducido, lo que explica la diferencia en el número de DEGs respecto a Dunn et al. (2026). Adicionalmente, al ser análisis de bulk RNA-seq, no es posible distinguir qué tipos celulares específicos dentro del hígado contribuyen a cada cambio transcriptómico; este punto sería abordable con los datos de scRNA-seq del mismo bioproyecto. Finalmente, el análisis se realizó únicamente en ratones adultos de 19–26 semanas, por lo que los resultados no pueden extrapolarse directamente a otros estadios del desarrollo. Como perspectiva futura, incluir las 6 réplicas completas, integrar los datos de scRNA-seq disponibles en el mismo bioproyecto, y comparar el perfil transcriptómico hepático con el análisis multi-tejido de Tian et al. (2023) permitiría una caracterización más completa del fenotipo Dp16.

12. Outlogs del pipeline

Los archivos de log de cada paso del pipeline están disponibles en el repositorio en scripts/out_logs/. La tabla enlaza directamente a cada archivo:

Paso Script Outlogs
Descarga SRA 01_download_adapters.sh download.out · download.err · download_Dp16.out · download_Dp16.err
Trimmomatic trimmomatic_final.sh trimm_final_1611.out · trimm_final_1611.err
Kallisto (×6 muestras) 03_kallisto.sh _0 · _1 · _2 · _3 · _4 · _5
STAR index 04_STAR_index.sh STAR_index.out · STAR_index.err
STAR align (×6 muestras) 05_STAR_align.sh _0 · _1 · _2 · _3 · _4 · _5

13. Referencias

  1. Dunn, A. et al., 2026. Altered hepatic metabolism in Down syndrome. Cell Reports. DOI: 10.1016/j.celrep.2025.116835

  2. Love, M.I., Huber, W., & Anders, S., 2014. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15, 550. DOI: 10.1186/s13059-014-0550-8

  3. Tian, Y. et al., 2023. Gene dosage imbalance disrupts systemic metabolism in the Dp16 Down syndrome mouse model. bioRxiv / eLife. DOI: 10.64898/2026.01.13.699318

  4. Bray, N.L., Pimentel, H., Melsted, P., & Pachter, L., 2016. Near-optimal probabilistic RNA-seq quantification. Nature Biotechnology, 34, 525–527. DOI: 10.1038/nbt.3519

sessionInfo()
## R version 4.5.2 (2025-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS Tahoe 26.3
## 
## Matrix products: default
## BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/Mexico_City
## tzcode source: internal
## 
## attached base packages:
## [1] grid      stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] apeglm_1.32.0               limma_3.66.0               
##  [3] VennDiagram_1.8.2           futile.logger_1.4.9        
##  [5] gprofiler2_0.2.4            ggrepel_0.9.8              
##  [7] kableExtra_1.4.0            knitr_1.51                 
##  [9] RColorBrewer_1.1-3          pheatmap_1.0.13            
## [11] reshape2_1.4.5              dplyr_1.2.1                
## [13] ggplot2_4.0.3               rhdf5_2.54.1               
## [15] tximport_1.38.2             DESeq2_1.50.2              
## [17] SummarizedExperiment_1.40.0 Biobase_2.70.0             
## [19] MatrixGenerics_1.22.0       matrixStats_1.5.0          
## [21] GenomicRanges_1.62.1        Seqinfo_1.0.0              
## [23] IRanges_2.44.0              S4Vectors_0.48.1           
## [25] BiocGenerics_0.56.0         generics_0.1.4             
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.1     viridisLite_0.4.3    farver_2.1.2        
##  [4] bitops_1.0-9         S7_0.2.2             RCurl_1.98-1.18     
##  [7] fastmap_1.2.0        lazyeval_0.2.3       digest_0.6.39       
## [10] lifecycle_1.0.5      statmod_1.5.1        magrittr_2.0.5      
## [13] compiler_4.5.2       rlang_1.2.0          sass_0.4.10         
## [16] tools_4.5.2          yaml_2.3.12          data.table_1.18.4   
## [19] lambda.r_1.2.4       labeling_0.4.3       S4Arrays_1.10.1     
## [22] htmlwidgets_1.6.4    curl_7.1.0           DelayedArray_0.36.1 
## [25] plyr_1.8.9           xml2_1.5.2           abind_1.4-8         
## [28] BiocParallel_1.44.0  numDeriv_2016.8-1.1  withr_3.0.2         
## [31] purrr_1.2.2          Rhdf5lib_1.32.0      MASS_7.3-65         
## [34] scales_1.4.0         mvtnorm_1.3-7        bbmle_1.0.25.1      
## [37] cli_3.6.6            rmarkdown_2.31       otel_0.2.0          
## [40] rstudioapi_0.18.0    httr_1.4.8           bdsmatrix_1.3-7     
## [43] cachem_1.1.0         stringr_1.6.0        parallel_4.5.2      
## [46] formatR_1.14         XVector_0.50.0       vctrs_0.7.3         
## [49] Matrix_1.7-5         jsonlite_2.0.0       systemfonts_1.3.2   
## [52] locfit_1.5-9.12      plotly_4.12.0        jquerylib_0.1.4     
## [55] tidyr_1.3.2          glue_1.8.1           emdbook_1.3.14      
## [58] codetools_0.2-20     stringi_1.8.7        gtable_0.3.6        
## [61] tibble_3.3.1         pillar_1.11.1        htmltools_0.5.9     
## [64] rhdf5filters_1.22.0  R6_2.6.1             textshaping_1.0.5   
## [67] evaluate_1.0.5       lattice_0.22-9       futile.options_1.0.1
## [70] bslib_0.10.0         Rcpp_1.1.1-1.1       coda_0.19-4.1       
## [73] svglite_2.2.2        SparseArray_1.10.10  xfun_0.57           
## [76] pkgconfig_2.0.3

Documento generado con R Markdown.
Análisis realizado en el cluster ken.lavis.unam.mx
Ruta del proyecto: /mnt/data/bioinfo-estadistica-2/RNAseq_2026/equipos/Equipo4/