## Cargar paquetes

```{r pqauetes, message=FALSE, warning=FALSE}
library(here)      # Facilitar localización de archivos
library(tidyverse) # Paquetes como dplyr, ggplot, etc
library(tximport)  # Importar datos de Kallisto en R 
library(DESeq2)   # Normalización y análisis de DEG
library(reshape2)
library(ggrepel)
# Manipulacion de GTF
library(GenomicFeatures)
library(txdbmaker) # BiocManager::install("txdbmaker")
```

Establecer mi ruta

```{r}
workdir <- here("Practica_Dia3")
```

## Paso 1. Crear el archivo `tx2gene` a partir del GTF

Descargue el GTF de [Ensembl](https://ftp.ensembl.org/pub/release-115/gtf/homo_sapiens/), ya que el transcriptoma de referencia proviene de esa fuente y nos interesa que los IDs coincidan.

```{r}
# Cargar el archivo GTF.gz
txdb <- txdbmaker::makeTxDbFromGFF(
  paste0(workdir, "/data/Homo_sapiens.GRCh38.115.gtf.gz"),
  format="gtf")
```

Nombre de las columnas del objetivo txdb

```{r}
columns(txdb)
```

Ahora si seleccionamos las columnas de interes de los transcritos (TXNAME) y genes (GENEID)

```{r}
# Extraer la información de transcritos con sus nombres y genes
transcripts <- transcripts(txdb, columns=c("TXNAME", "GENEID"))

# Construir el data.frame con las columnas deseadas
tx2gene <- data.frame(
  TXNAME  = transcripts$TXNAME,
  GENEID  = transcripts$GENEID
) %>% 
  # Seleccionar columnas de interes y renombrar columna GENEID
  dplyr::select(TXNAME, GENEID = GENEID.value)

# Ver los primeros registros
head(tx2gene)
```

## Paso 2. Importar datos de salida de `kallisto` a R

### Opción A. Importar archivo `.h5` de `kallisto` (se cargan a nivel de transcritos)

```{r}
# Definir la carpeta donde están los resultados
ruta <- "/Users/ecoss/Documents/RNAseq_classFEB2026/Practica_Dia3/counts/kallisto_quant"
# Listar todos los archivos .h5 dentro de las subcarpetas
archivos_h5 <- list.files(path = ruta,
                          pattern = "\\.h5$",
                          full.names = TRUE,
                          recursive = TRUE)
# Extraer el nombre SRR de la ruta (basename de la carpeta)
nombres_srr <- basename(dirname(archivos_h5))
# colocarle nombres
names(archivos_h5) <- nombres_srr
# importar archivo
txi.kallisto <- tximport(archivos_h5, type = "kallisto", txOut = TRUE)
```

En este caso lo que ocurre es que `txOut = TRUE` **le dice a** `tximport` **que no colapse los transcritos a nivel de gen**, sino que te devuelva directamente las abundancias de cada transcrito.

Número de transcritos detectados:

```{r}
all_transcripts <- rownames(txi.kallisto$abundance)
length(all_transcripts)
```

### Opción B. Importar archivo `.tsv` de `kallisto` (se cargan a nivel de genes) ⭐ Manera correcta

```{r}
# Definir la carpeta donde están los resultados
ruta <- "/Users/ecoss/Documents/RNAseq_classFEB2026/Practica_Dia3/counts/kallisto_quant"
# Listar todos los archivos .h5 dentro de las subcarpetas
archivos_tsv <- list.files(path = ruta,
                          pattern = "\\.tsv$",
                          full.names = TRUE,
                          recursive = TRUE)
# Extraer el nombre SRR de la ruta (basename de la carpeta)
nombres_srr <- basename(dirname(archivos_tsv))
# colocarle nombres
names(archivos_tsv) <- nombres_srr
# importar archivo
txi.kallisto <- tximport(archivos_tsv, type = "kallisto", tx2gene = tx2gene, ignoreTxVersion=TRUE)
names(txi.kallisto)
```

::: callout-note
### Usar `ignoreAfterBar=TRUE`

Algunas anotaciones tienen IDs como `ENST00000622028.1|extraInfo`. En ese caso, `ignoreAfterBar=TRUE` ayuda a cortar lo que está después del `|`.
:::

Seleccionar todos los genes contenidos en esta variable

```{r}
all_genes <- rownames(txi.kallisto$abundance)
length(unique(all_genes))
```

## Paso 3. Filtrar biotipos

Decidi emplear el ([Ensembldb.RData](https://github.com/NeuroGenomicsMX/SLE_Transcriptomic_analysis_tolerogenic/blob/main/data/Ensembldb.RData)). Cargando la variable annotations_ENSG en nuestro ambiente, para más información revisa el script [Ensembl_database.R](https://github.com/NeuroGenomicsMX/bulkRNA_lupus2024/blob/main/scripts/Ensembl_database.R).

```{r}
load(paste0(workdir, "/data/Ensembldb.RData"))
```

Filter `txi` related with biotypes from Ensembl.

```{r}
hsapiens_biotype <- dplyr::select(annotations_ENSG, ensembl_gene_id, hgnc_symbol, gene_biotype)
hsapiens_biotype <- hsapiens_biotype %>% distinct() # eliminar duplciados
colnames(hsapiens_biotype)[1] <- "gene"

# This retrieves all gene identifiers from the txi object (the output of tximport), specifically from the abundance matrix.

my_biotypes <- hsapiens_biotype[hsapiens_biotype$gene %in% all_genes,]
dim(my_biotypes) # [1] 21906     2

# We filter the annotation data to keep only genes that are present in the expression data. The result is a data frame (my_biotypes) with 21,949 annotated genes.

levels(as.factor(my_biotypes$gene_biotype))
```

In this step, we filter out unwanted gene biotypes from the annotation data to focus on more relevant gene categories for our analysis:

```{r}
filtered_biotypes <- my_biotypes[!(my_biotypes$gene_biotype %in% c('rRNA', 'ribozyme')), ]
dim(filtered_biotypes) # [1] 21887     2

# Next, we extract the list of filtered gene names and ensure they match those present in the expression dataset:

filtered_genes <- filtered_biotypes$gene
filtered_genes <- all_genes[all_genes %in% filtered_genes] 
length(filtered_genes) # 77786

# This code filters the txi object to retain only the genes that passed the biotype filtering, preparing the dataset for downstream analysis:

# Remove counts of rRNA
txi_minuslast <- txi.kallisto[-length(txi.kallisto)]

filtered_txi <- lapply(txi_minuslast, function(df) {
  # print(head(rownames(df)))
  df[rownames(df) %in% filtered_genes, ]
})

names(filtered_txi) <- c("abundance" , "counts","length")
filtered_txi$countsFromAbundance <- "no"
dim(filtered_txi$counts) 
```

Renombrar nombres de filas por nombres symbolicos

```{r}
# Supongamos que filtered_txi$counts es tu matriz de expresión
counts_df <- filtered_txi$counts %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene")   # convertir rownames a columna

# Unir por la columna 'gene'
counts_annot <- counts_df %>%
  left_join(hsapiens_biotype, by = "gene")

# Ahora puedes usar hgnc_symbol como rownames
# NOTA: no se pude correr esto porque los geneID se repiten
# counts_annot <- counts_annot %>%
#   tibble::column_to_rownames("hgnc_symbol")

head(counts_annot)
```

Almacenar la salida

```{r}
# Guardar variable
save(filtered_txi, file = paste0(workdir, '/counts/kallisto_outputs/txi_filtered_v2.RData'))
# Guardar en csv
write_csv(as.data.frame(filtered_txi$counts), file =  paste0(workdir, '/counts/kallisto_outputs/txi_filtered_v2.csv'))
```

## Paso 4. Crear/Cargar la metadata

En este punto las cuentas tienen

```{r}
head(filtered_txi$counts)
```

Crear mi metadata a partir de la información recolectada de NCBI, GEO y ENA/EBI. Todas son muestras del día 7 posterior a la vacunación de la vacuna de COVID (BNT162b2 SARS-CoV-2). Si quieres repasar esta parte puedes revisar las [diapositivas](https://eveliacoss.github.io/RNAseq_classFEB2026/Presentaciones/Dia3_ImportarDatos.html#56) proporcionadas previamente.

```{r}
metadata <- data.frame(
  sample_ID = c("SRR27190676", "SRR27190684", "SRR27190700", 
          "SRR27190685", "SRR27190708", "SRR27190710"),
  sample_title = c("p18.S028.day7", "p17.S023.day7", "p12.S013.day7", 
                   "p14.S018.day7", "p1.S003.day7", "p8.S008.day7"),
  responder = c(rep("NR", 3), rep("R", 3)),
  condition = "SLE",
  day = "d7"
)

head(metadata) %>% 
  DT::datatable()
```

## Paso 5. Crear el objeto DESEq2

```{r}
dds <- DESeqDataSetFromTximport(txi = filtered_txi,
                                colData = metadata,
                                design = ~responder)

dds

dim(dds) # [1]  77786   6
```

Guardar objeto

```{r}
save(dds, file = paste0(workdir, "/counts/kallisto_outputs/dds_cleaned.RData"))
```

## Paso 6. Normalización de las cuentas con rlog (log regularizado)

```{r}
counts <- assays(dds)$counts
head(counts(dds))[1:4, 1:4]
```

We first extract the raw count matrix from the `DESeqDataSet` object using `assays(dds)$counts`.

```{r}
rlog_result <- rlog(counts, blind=FALSE)
head(rlog_result)[1:4, 1:4]
```

Almacenar counts y cuentas transformadas con `rlog`:

```{r}
# Save vsd
save(rlog_result, file = paste0(workdir, "/counts/kallisto_outputs/rlog.RData"))
# Save counts
save(counts, file = paste0(workdir, "/counts/kallisto_outputs/rawCounts.RData"))
```

The `rlog()` function is then applied to perform a regularized log transformation. Setting `blind = FALSE` means that the transformation will take into account the experimental design, which is recommended when the goal is visualization and downstream comparisons. Unlike a simple log, the regularization stabilizes variance across the range of counts, making low-count genes less noisy and improving interpretability in PCA, clustering, and heatmaps.

## Paso 7. Melting the Normalized Count Matrix for Visualization

To facilitate downstream plotting and statistical analysis (e.g., violin plots, boxplots, or faceted visualizations), we reshape the normalized expression data from wide to long format.

```{r}
# vsd all info
rlog_all <- rlog(dds, blind=FALSE)
melted_norm_counts <- data.frame(melt(assay(rlog_all)))
# rename columns
colnames(melted_norm_counts) <- c("gene", "sample_ID", "normalized_counts")
# join metadata
melted_norm_counts <- left_join(melted_norm_counts, metadata, by = "sample_ID")
```

Almacenamiento:

```{r}
# Save melted_norm_counts
save(melted_norm_counts, file = paste0(workdir,  "/counts/kallisto_outputs/melted_norm_counts.RData"))
# Save vsd
save(rlog_all, file = paste0(workdir,  "/counts/kallisto_outputs/rlog_all.RData"))
```

## Paso 8. Compute z-score transformation

To facilitate clustering and heatmap visualization, we compute the `z-score transformation` of the raw count matrix. Z-scores help standardize expression values across genes, highlighting relative expression patterns.

```{r}
# filter out genes that have low variance across samples
var_non_zero <- apply(counts, 1, var) !=0 
filtered_counts <- counts[var_non_zero, ]
  
zscores <- t(scale(t(filtered_counts)))
dim(zscores) #[1] 18551    99

zscore_mat <- as.matrix(zscores)
```

```{r}
# Save melted_norm_counts with z-scores
save(zscore_mat, file = paste0(workdir, "/counts/kallisto_outputs/zscore_mat.RData"))
```

## Paso 9. Visualización de PCA

Principal Component Analysis (PCA) of Normalized Expression Data

After applying rlog, we perform Principal Component Analysis (PCA) to explore the global structure of the data and detect potential clustering patterns among the samples.

```{r}
mat <- as.matrix(rlog_result) # col = samples, rows  = genes
pc <- prcomp(t(mat)) # col = genes , rows = samples
```

-   The VST-normalized matrix is converted to a standard matrix format.
-   We transpose the matrix (`t(mat)`) because `prcomp()` expects observations (samples) as rows and variables (genes) as columns.

We extract the first two principal components (PC1 and PC2) and combine them with selected metadata columns, such as sample group and clinical information. And the sample names are added as a new column, and the grouping variable is converted to a **factor** to facilitate plotting and analysis.

```{r}
# PCA analysis
pc$sample <- metadata$sample_ID # rename sample

# This calculates the percentage of variance explained by each principal component, useful for labeling axes in PCA plots.
variance_explained <- pc$sdev^2 / sum(pc$sdev^2) * 100
```

¿Cuántos PCAs calculo?

```{r}
head(pc$x)
```

```{r}
df <- data.frame(
  Componente = seq_along(variance_explained),
  Varianza = variance_explained,
  Etiqueta = paste0("PC", seq_along(variance_explained), ": ", round(variance_explained, 2), "%")
)

# Graficar con ggplot2 y ggrepel
ggplot(df, aes(x = Componente, y = Varianza)) +
  geom_point(color = "steelblue", size = 3) +
  geom_line(color = "steelblue") +
  geom_text_repel(data = df[1:4, ], aes(label = Etiqueta),
                  size = 5, box.padding = unit(0.35, "lines"),
                  point.padding = unit(0.3, "lines"),
                  color = "darkred", fontface = "bold") +
  labs(title = "Scree plot: varianza explicada por PCA",
       x = "Componente principal",
       y = "Porcentaje de varianza explicada") +
  theme_minimal()
```

### Color by `responder`

This PCA data frame (`pca_df`) is ready for visualization, allowing us to assess whether samples cluster according to biological or experimental groups.

```{r}
# join information
pca_df <- cbind(pc$x[,1:2], metadata[,c(2:4)]) %>% as.data.frame()
pca_df$sample <- rownames(pca_df)

# Chec information
levels(as.factor(pca_df$responder))
# [1]  "NR" "R"
```

Determinar paleta de colores

```{r}
responder_colors <- c(
  "NR"      = "black",
  "R"  = "palevioletred"
)
```

Visualizacion de los datos

```{r}
p <- ggplot(pca_df,aes(x = PC1, y = PC2, color = responder)) +
  stat_ellipse() +
  geom_point(aes(shape = condition), size = 3.5) + 
  scale_shape_manual(values = 17) +
  scale_color_manual(values = responder_colors) +
  labs(x = sprintf('PC1 (%.2f%%)', variance_explained[1]),
       y = sprintf('PC2 (%.2f%%)',variance_explained[2]),
       color = "Respuesta",
       shape = "Condición") +
  theme_light() +
  theme(plot.background = element_rect(fill = "white"),
        text = element_text(size = 12),
        axis.title = element_text(size = 12),
        legend.title = element_text(size = 12),
        legend.text = element_text(size = 12),
        legend.key.size = unit(1.5, "lines"))
p
```

Guardar figura con buena calidad

```{r}
# ggsave(p, file = paste0(workdir,'figures/PCA.png'), dpi = 300)
ggsave(p, file = here("Practica_Dia3/figures","PCA.png"), dpi = 300)
```

## Paso 10. Corrección de batch effect

```{r}
mat <- assay(rlog_all)
mm <- model.matrix(~responder, colData(rlog_all))
mat <- limma::removeBatchEffect(mat, batch=rlog_all$batch, design=mm)
assay(rlog_all) <- mat
plotPCA(rlog_all, intgroup= c("responder", "condition"))
```

The `design` argument is necessary to avoiding removing variation associated with the treatment conditions. See `?removeBatchEffect` in the *limma* package for details.

```{r}
pcaData <- plotPCA(rlog_all, intgroup= c("responder", "condition"), returnData=TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))
ggplot(pcaData, aes(PC1, PC2, color=responder, shape=condition)) +
  geom_point(size=3) +
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) + 
  coord_fixed()
```

## Paso 11. Conclusiones y estrategias a realizar

Las seis muestras analizadas corresponden a pacientes con lupus (SLE) que presentan respuesta a una vacuna. Tras aplicar la corrección por *batch effect,* el análisis de componentes principales (PCA) muestra que no existe una separación clara entre los grupos de pacientes con respuesta (R) y pacientes sin respuesta (NR).

Este resultado refleja tanto una **limitación estadística** como una **variabilidad biológica inherente**: con solo tres réplicas por condición, el PCA tiende a capturar principalmente la variabilidad individual de cada paciente, más que las diferencias entre grupos. En consecuencia, la falta de separación no necesariamente implica ausencia de diferencias biológicas, sino que puede deberse al tamaño reducido de la cohorte y a la heterogeneidad propia de la enfermedad.

::: callout-note
## 🔎 ¿Por qué ocurre?

-   **Número de muestras**: con solo 3 por grupo, la variabilidad biológica de cada paciente pesa mucho.

-   **PCA es exploratorio**: no está diseñado para forzar separación de grupos, sino para mostrar la variabilidad global.

-   **Efecto de outliers**: una sola muestra atípica puede “arrastrar” la nube de puntos.
:::

### ¿Qué estrategias podemos realizar?

1.  **Aumentar el número de muestras**

    -   Es la forma más sólida de ganar poder estadístico y ver si la separación es consistente.
    -   Con más pacientes, el PCA refleja mejor la tendencia global y no tanto el ruido individual.

2.  **Complementar con análisis supervisados**

    -   Usa **DESeq2** o **edgeR** para pruebas de expresión diferencial, que sí modelan condiciones.
    -   PCA no es prueba de hipótesis, solo exploración.

3.  **Seleccionar genes más informativos**

    -   Haz PCA con los genes más variables o con los que resultan significativos en DE.
    -   Esto puede mejorar la separación visual.

4.  **Otras visualizaciones**

    -   Prueba **t-SNE** o **UMAP**, que a veces capturan mejor la estructura de grupos pequeños.
    -   Haz clustering jerárquico para ver si las muestras se agrupan como esperas.
