bulk_01_exploration.R

# Copyright (c) [2025] [Ines Rivero Garcia]
# ines.rivero@uni-heidelberg.de

setwd("/mnt/sds-hd/sd22b002/projects/ines/heart_revremod/")

library(dplyr)
library(DESeq2)
library(ggplot2)
library(ComplexHeatmap)
library(tidyr)
library(ggpubr)

# Load data and prep metadata
mRNA <- read.csv("data/CRC_bulk_readcounts.csv", header = TRUE, sep = ";")
metadata <- read.csv("data/bulk_metadata.csv", header = TRUE, sep = ";")
metadata <- metadata %>% mutate(condition_simplified = case_when(
    condition == "O" ~ "ORAB",
    condition == "DB" ~ "DB",
    condition %in% c("SO", "SDB") ~ "Control"
    )
  )
metadata$condition <- factor(metadata$condition, levels = c("SDB", "O", "DB", "SO"))
metadata$condition_simplified <- factor(metadata$condition_simplified, levels = c("Control", "ORAB", "DB"))

condition_colors <- c("O" = "#AD1332", "DB" = "#0571B0", "SO" = "#D99AA7", "SDB" = "#92C5DE")
condition_simplified_colors <- c("ORAB" = "#AD1332", "DB" = "#0571B0", "Control" = "grey50")

# Create DESeq object
mRNA <- aggregate(mRNA[ , -1], list(geneID = mRNA[ , 1]), FUN = sum)
rownames(mRNA) <- mRNA$geneID
mRNA <- mRNA[,-1] # Remove column with "gene_ID"
keep_genes <- rowSums(mRNA > 10) >= 3
mRNA <- mRNA[keep_genes, ]
rownames(metadata) <- metadata$sample
dds <- DESeqDataSetFromMatrix(countData = as.matrix(mRNA),
                              colData = metadata, 
                              design = ~ condition_simplified)
vsd <- vst(dds, blind = TRUE)

# Exploratory analysis: PCA
vars <- rowVars(assay(vsd))
top500 <- head(order(vars, decreasing = TRUE), 500)
pca_mat <- t(assay(vsd)[top500, ])

# Save HVGs for future use
hvgs <- rownames(assay(vsd)[top500,])
saveRDS(hvgs, "results/bulk/mrna_top500_hvg.rds")

# Visualise PCA
pca_res <- prcomp(pca_mat, center = TRUE, scale. = FALSE)
pca_df <- as.data.frame(pca_res$x)
pca_df$sample <- rownames(pca_df)
pca_df <- merge(pca_df, metadata, by = "sample")
var_explained <- pca_res$sdev^2 / sum(pca_res$sdev^2)
var_df <- data.frame(
  PC = paste0("PC", 1:length(var_explained)),
  Variance = var_explained * 100
)

pca_df$condition_simplified <- factor(pca_df$condition_simplified, levels = c("DB", "ORAB", "Control"))

pdf("results/bulk/mrna_scatterplot_pca.pdf", height = 4, width = 5)
ggplot(pca_df, aes(x = PC1, y = PC2, colour =  condition_simplified)) +
  geom_point(size = 3, alpha = 0.8) +
  theme_bw() +
  scale_color_manual(values = condition_simplified_colors) +
  labs(
  x = paste("PC1 (", round(var_explained[1], 2) * 100, "%)", sep = ""),
  y = paste("PC2 (", round(var_explained[2], 2) * 100, "%)", sep = ""),
  colour = "Condition"
  )
dev.off()

pca_df$condition <- factor(pca_df$condition, levels = c("O", "DB", "SO", "SDB"))

pdf("results/bulk/mrna_scatterplot_pca_allgroups.pdf", height = 4, width = 5)
ggplot(pca_df, aes(x = PC1, y = PC2, colour =  condition)) +
  geom_point(size = 3, alpha = 0.8) +
  theme_bw() +
  scale_color_manual(values = condition_colors) +
  labs(
  x = paste("PC1 (", round(var_explained[1], 2) * 100, "%)", sep = ""),
  y = paste("PC2 (", round(var_explained[2], 2) * 100, "%)", sep = "")
  )
dev.off()

pdf("results/bulk/mrna_explained_variance_pca.pdf", height = 4, width = 7)
ggplot(var_df, aes(x = reorder(PC, -Variance), y = Variance)) +
  geom_bar(stat = "identity", fill = "grey50") +
  theme_classic(base_size = 12) +
  labs(
    title = "Explained Variance per Principal Component",
    x = "Principal Component",
    y = "Variance Explained (%)"
  ) +
  scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 10))
dev.off()

write.table(pca_df, "results/bulk/mrna_pca_coordinates.txt", col.names = TRUE, row.names = FALSE, quote = FALSE, sep = "\t")

# Explore how other PCs separate data
pdf("results/bulk/mrna_pairplot_allpc_allgroups.pdf", height = 20, width = 20)
GGally::ggpairs(data = pca_df, 
                columns = c(paste0("PC", 1:16), "condition"), 
                mapping = aes(colour = condition)) + 
  theme_bw() + 
  scale_color_manual(values = condition_colors) +
  scale_fill_manual(values = condition_colors)
dev.off()

# Boxplot of weights per PC and stat test
pca_df_long <- pca_df %>% 
    select(c(starts_with("PC"), sample, condition_simplified)) %>% 
    reshape2::melt(c("sample", "condition_simplified"))

pca_df_long$condition_simplified <- factor(pca_df$condition_simplified, levels = c("ORAB", "DB", "Control"))

pdf("results/bulk/mrna_boxplot_pcs.pdf", height = 8, width = 9)
ggplot(pca_df_long, aes(x = condition_simplified, y = value, fill = condition_simplified)) +
  geom_boxplot() +
  geom_jitter() +
  facet_wrap(~ variable) +
  theme_bw(base_size = 12) +
  scale_fill_manual(values = condition_simplified_colors)
dev.off()

# Check associations between PCs and conditions
pca_pvalues <- c()
for(i in 1:16){
    aov_result <- aov(as.formula(paste(paste0("PC",i), " ~ condition_simplified")), data = pca_df)
    pca_pvalues <- c(pca_pvalues,
                     summary(aov_result)[[1]][["Pr(>F)"]][1])                   
}
pca_pvalues <- p.adjust(pca_pvalues, method = "bonferroni")          
pca_pvalues 

# Boxplot of PC1 with pvalues
pc1_anova <- aov(PC1 ~ condition_simplified, data = pca_df)
pc1_pvals <- TukeyHSD(pc1_anova) %>%
    broom::tidy() %>%
    as_tibble()
pc1_pvals <- pc1_pvals %>%
    mutate(signif_label = case_when(
        adj.p.value < 0.001 ~ "***",
        adj.p.value < 0.01  ~ "**",
        adj.p.value < 0.05  ~ "*",
        TRUE          ~ "ns"
    )) %>%
    mutate(y_position = rep(max(pca_df_long[pca_df_long$variable == "PC1", "value"]), 3)) %>%
    separate(contrast, into = c("group1", "group2"), sep = "-")

pdf("results/bulk/mrna_boxplot_pc1.pdf", height = 3, width = 4)
pca_df_long %>%
    filter(variable == "PC1") %>%
    ggplot(aes(x = condition_simplified, y = value, fill = condition_simplified)) +
        geom_boxplot() +
        geom_jitter() +
        stat_pvalue_manual(data = pc1_pvals, 
                           label = "signif_label",
                           xmin = "group1",
                           xmax = "group2",
                           y.position = "y_position",
                           tip.length = 0.01,
                           bracket.size = 0.3,
                           size = 5,
                           inherit.aes = FALSE) +
        labs(x = "Condition", y = "PC1 value", fill = "Condition") +
        theme_bw(base_size = 12) +
        scale_fill_manual(values = condition_simplified_colors)
dev.off()

# Exploratory analysis: sample clustering
dist <- dist(t(assay(vsd)))
dist.m <- as.matrix(dist)
dist.df <- dist.m %>%
  as.data.frame() %>%        
  as_tibble(rownames = "sample1") %>%
  pivot_longer(
    -sample1,
    names_to = "sample2",
    values_to = "distance"
  ) %>%
  filter(sample1 < sample2)   # keep each pair only once but exclude comparisons in which sample1 = sample2 because their distance will be 0.

# Add type of comparisons
dist.df <- dist.df %>%
  mutate(
      comparison_simplified = case_when(
          startsWith(sample1, "O_") & startsWith(sample2, "DB_") ~ "ORAB vs DB",
          startsWith(sample1, "DB_") & startsWith(sample2, "O_") ~ "ORAB vs DB",
          
          startsWith(sample1, "O_") & startsWith(sample2, "S")  ~ "ORAB vs Control",
          startsWith(sample1, "S")  & startsWith(sample2, "O_") ~ "ORAB vs Control",

          startsWith(sample1, "DB_") & startsWith(sample2, "S")  ~ "DB vs Control",
          startsWith(sample1, "S")  & startsWith(sample2, "DB_") ~ "DB vs Control",
          
          startsWith(sample1, "O_") & startsWith(sample2, "O_") ~ "ORAB vs ORAB",
          startsWith(sample1, "DB_") & startsWith(sample2, "DB_") ~ "DB vs DB",
          startsWith(sample1, "S")  & startsWith(sample2, "S")  ~ "Control vs Control",
      TRUE ~ NA_character_
    ),
      comparison =  case_when(
          startsWith(sample1, "O_") & startsWith(sample2, "DB_") ~ "ORAB vs DB",
          startsWith(sample1, "DB_") & startsWith(sample2, "O_") ~ "ORAB vs DB",
          
          startsWith(sample1, "O_") & startsWith(sample2, "SO_")  ~ "ORAB vs Sham-ORAB",
          startsWith(sample1, "SO_")  & startsWith(sample2, "O_") ~ "ORAB vs Sham-ORAB",

          startsWith(sample1, "O_") & startsWith(sample2, "SDB_") ~ "ORAB vs Sham-DB",
          startsWith(sample1, "SDB_") & startsWith(sample2, "O_") ~ "ORAB vs Sham-DB",

          startsWith(sample1, "DB_") & startsWith(sample2, "SDB")  ~ "DB vs Sham-DB",
          startsWith(sample1, "SDB_")  & startsWith(sample2, "DB_") ~ "DB vs Sham-DB",

          startsWith(sample1, "DB_") & startsWith(sample2, "SO")  ~ "DB vs Sham-ORAB",
          startsWith(sample1, "SO_")  & startsWith(sample2, "DB_") ~ "DB vs Sham-ORAB",

          startsWith(sample1, "SDB_") & startsWith(sample2, "SO")  ~ "Sham-DB vs Sham-ORAB",
          startsWith(sample1, "SO_")  & startsWith(sample2, "SDB_") ~ "Sham-DB vs Sham-ORAB",
          
          startsWith(sample1, "O_") & startsWith(sample2, "O_") ~ "ORAB vs ORAB",
          startsWith(sample1, "DB_") & startsWith(sample2, "DB_") ~ "DB vs DB",
          startsWith(sample1, "SO_")  & startsWith(sample2, "SO_")  ~ "Sham-ORAB vs Sham-ORAB",
          startsWith(sample1, "SDB_") & startsWith(sample2, "SDB_") ~ "Sham-DB vs Sham-DB",
      TRUE ~ NA_character_
    )
  )

# Test if there are differences in distances between groups
kruskal.test(distance ~ comparison_simplified, data = dist.df)
wilcox_result <- pairwise.wilcox.test(dist.df$distance, dist.df$comparison_simplified, p.adjust.method = "bonferroni")
pval_matrix <- wilcox_result$p.value
write.table(pval_matrix, "results/bulk/mrna_wilcoxtext_padj_matrix_sample_distances.tsv", col.names = TRUE, row.names = TRUE, sep = "\t", quote = FALSE)

# Clean significant padjs for nice plotting
sig_pvals_df <- pval_matrix %>%
    as.data.frame() %>%
    tibble::rownames_to_column("group1") %>%
    pivot_longer(
        cols = -group1,
        names_to = "group2",
        values_to = "p.adj"
    ) %>%
    filter(!is.na(p.adj)) %>% # keep only comparisons that exist
    filter(p.adj < 0.05) %>% # keep significant comparisons only
    mutate(
        signif_label = case_when(
            p.adj < 0.001 ~ "***",
            p.adj < 0.01  ~ "**",
            p.adj < 0.05  ~ "*"
    )
  )

# Add y-position for nice plotting
y_max <- dist.df %>%
  group_by(comparison_simplified) %>%
  summarise(max_y = max(distance))

sig_pvals_df <- sig_pvals_df %>%
  left_join(y_max, by = c("group1" = "comparison_simplified")) %>%
  left_join(y_max, by = c("group2" = "comparison_simplified")) %>%
  mutate(y_position = pmax(max_y.x, max_y.y) + 0.05 * max(dist.df$distance)) %>% # add spacing
  select(group1, group2, p.adj, signif_label, y_position)

pdf("results/bulk/mrna_boxplot_group_distances.pdf", height = 4, width = 7)
ggplot(dist.df, aes(x = reorder(comparison_simplified, -distance, FUN = median), y = distance)) + 
    geom_boxplot(alpha = 0.3, width = 0.65, fill = "grey50") + 
    geom_jitter(width = 0.2) +
    theme_bw() +
    ylim(c(0,70)) +
    stat_pvalue_manual(data = sig_pvals_df, label = "signif_label", xmin = "group1", xmax = "group2",y.position = "y_position")+
    ylab("Euclidean distance between pairs of samples") +
    xlab("Comparison") +
    theme(legend.position = "right")
dev.off()

pdf("results/bulk/mrna_boxplot_group_distances_allgroups.pdf", height = 6, width = 5)
ggplot(dist.df, aes(x = reorder(comparison, -distance, FUN = median), y = distance)) + 
    geom_boxplot(alpha = 0.3, width = 0.65, fill = "grey50") + 
    geom_jitter(width = 0.2) +
    theme_bw() +
    ylim(c(0,70)) +
    ylab("Euclidean distance between pairs of samples") +
    xlab("Comparison") +
    theme(legend.position = "right", axis.text.x = element_text(angle = 90, hjust = 1))
dev.off()