bulk_03_trajectory_clusters.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)
library(dendextend)

# 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 <- metadata %>% mutate(time = case_when(
    condition_simplified == "ORAB" ~ 2,
    condition_simplified == "DB" ~ 3,
    condition_simplified == "Control" ~ 1
    )
  )
metadata$condition <- factor(metadata$condition, levels = c("SDB", "O", "DB", "SO"))
metadata$condition_simplified <- factor(metadata$condition_simplified, levels = c("Control", "ORAB", "DB"))
metadata$time <- factor(metadata$time, levels = c(1, 2, 3))

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

geneinfo <- read.csv("data/mm39_geneannotation_dowload202507.txt", header = TRUE, sep = "\t")

# 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)

# Load DEG tables
degs_DB_vs_Control <- read.csv("results/bulk/mrna_degs_DB_vs_Control.csv", header = TRUE) %>%
    filter(padj < 0.05) %>%
    pull(EnsemblID)
degs_ORAB_vs_Control <- read.csv("results/bulk/mrna_degs_ORAB_vs_Control.csv", header = TRUE) %>%
    filter(padj < 0.05) %>%
    pull(EnsemblID)
degs_DB_vs_ORAB <- read.csv("results/bulk/mrna_degs_DB_vs_ORAB.csv", header = TRUE) %>%
    filter(padj < 0.05) %>%
    pull(EnsemblID)

degs <- unique(c(degs_DB_vs_Control, degs_ORAB_vs_Control, degs_DB_vs_ORAB))

# Get matrix of DEGs 
norm_counts <- as.data.frame(assay(vsd))
norm_counts$EnsemblID <- rownames(norm_counts)
norm_counts_degs <- norm_counts %>% filter(EnsemblID %in% degs) 

# Summarise counts 
trans_norm_mean <- norm_counts_degs %>% 
  # convert to long format
  pivot_longer(cols = SO_2:DB_35, names_to = "sample", values_to = "normcounts")  %>% 
  # join with sample info table
  full_join(metadata, by = ("sample")) %>% 
  # for each gene
  group_by(EnsemblID) %>% 
  # scale the cts column
  mutate(normcounts_scaled = (normcounts - mean(normcounts))/sd(normcounts)) %>% 
  # for each gene, strain and time
  group_by(EnsemblID, condition_simplified, time) %>%
  # calculate the mean (scaled) cts
  summarise(mean_normcounts_scaled = mean(normcounts_scaled),
            nrep = n()) %>% 
  ungroup()

# Prep matrix for clustering
hclust_matrix <- norm_counts_degs %>% 
  select(-c(EnsemblID)) %>% 
  as.matrix()
rownames(hclust_matrix) <- norm_counts_degs$EnsemblID

# Scale matrix
hclust_matrix <- hclust_matrix %>% 
  t() %>% 
  scale() %>% 
  t()

# Calculate distances
gene_dist <- dist(hclust_matrix)

# Perform hierarchical clustering
gene_hclust <- hclust(gene_dist, method = "complete")

# Try several k values
k_values <- 2:12
for(k in k_values){
    gene_cluster <- cutree(gene_hclust, k = k) %>% 
        tibble::enframe()
    colnames(gene_cluster) <- c("EnsemblID", "cluster")
    gene_cluster <- merge(gene_cluster, geneinfo[, c("Gene.stable.ID", "Gene.name")], by.x = "EnsemblID", by.y = "Gene.stable.ID", all.x = TRUE)
    trans_norm_cluster <- merge(trans_norm_mean, gene_cluster, by = "EnsemblID", all.x = TRUE)
    pdf(paste0("results/bulk/mrna_time_clusters_k", k, ".pdf"), height = 3, width = 12)
    print(ggplot(trans_norm_cluster, aes(x = time,  y = mean_normcounts_scaled)) +
              geom_line(aes(group = EnsemblID)) +
              geom_line(stat = "summary", fun = "median", colour = "brown", size = 1, aes(group = 1)) +
              facet_wrap(~ cluster, nrow = 1) +
              scale_x_discrete(labels = c("Sham", "ORAB", "DB")) + 
              theme_bw() +
              ylab("Mean scaled expression") +
              xlab(""))
    dev.off()

    dend <- gene_hclust %>% as.dendrogram
    pdf(paste0("results/bulk/mrna_time_clusters_dendogram_k", k, ".pdf"), height = 3, width = 3)
    dend %>% set("branches_k_color", k = k) %>%  plot(main = paste0("k = ", k), leaflab = 'none')
    dev.off() 
}

# Choose k = 6 as final nr clusters
k <- 6
gene_cluster <- cutree(gene_hclust, k = k) %>% 
    tibble::enframe()
colnames(gene_cluster) <- c("EnsemblID", "cluster")
gene_cluster <- merge(gene_cluster, geneinfo[, c("Gene.stable.ID", "Gene.name")], by.x = "EnsemblID", by.y = "Gene.stable.ID", all.x = TRUE)


# Test differences in avg expression per cluster across timepoints
cluster_expr <- sapply(1:6, function(x) {
  genes_k <- gene_cluster %>% filter(cluster == x) %>% pull(EnsemblID)
    colMeans(norm_counts_degs[genes_k, 1:16,]) # Last column is EnsemblID
})
cluster_expr <- as.data.frame(cluster_expr)
cluster_expr$sample <- rownames(cluster_expr)
cluster_expr<-merge(cluster_expr, metadata, by ='sample')
cluster_expr<-cluster_expr%>% select(sample, V1, V2, V3, V4, V5, V6, time)
cluster_expr<-reshape2::melt(cluster_expr, id.vars=c('sample', 'time'))
for(k in 1:6){
    cluster <- paste0("V",k)
    fit <- aov(value ~ time, data = cluster_expr[cluster_expr$variable==cluster,])
    print(summary(fit))
    print(TukeyHSD(fit))
}

# Save table of genes per cluster to do enrichment analysis
write.table(gene_cluster, "results/bulk/mrna_time_clusters.csv", col.names = TRUE, row.names = FALSE, quote = FALSE, sep = ",")

# Heatmap scaled expression genes grouped by cluster
scaled_mat <- t(scale(t(norm_counts_degs[,1:16])))  
col_anno <- HeatmapAnnotation(Condition = metadata$condition_simplified,
                              col = list(Condition = condition_simplified_colors))
gene_cluster_ensembl <- gene_cluster[, c("EnsemblID", "cluster")]
gene_cluster_ensembl <- unique(gene_cluster_ensembl)
row_anno <- rowAnnotation(Cluster = gene_cluster_ensembl$cluster,
                          col = list(Cluster = c("1" = "#66c2a5", "2" = "#fc8d62", "3" = "#8da0cb", "4" = "#e78ac3", "5" = "#a6d854", "6" = "#ffd92f")))

pdf("results/bulk/mrna_time_clusters_heatmap.pdf", height = 5, width = 4)
Heatmap(
  scaled_mat,
  name = "scaled\nexpression",
  show_row_names = FALSE,
  show_column_names = FALSE,
  cluster_rows = TRUE,
  cluster_columns = TRUE,
  bottom_annotation = col_anno,
  right_annotation = row_anno,
  row_split = gene_cluster_ensembl$cluster,   
  col = circlize::colorRamp2(c(-2, 0, 2), c("blue", "white", "red"))
)
dev.off()