heart_revremod_paper/analysis/bulk/bulk_06_clusters_example-genes_visualisation.R at main · saezlab/heart_revremod_paper · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
# Copyright (c) [2025] [Ines Rivero Garcia]
# ines.rivero@uni-heidelberg.de

# Visualise mRNA-seq funcomics results

library(dplyr)
library(DESeq2)
library(ggplot2)
library(ComplexHeatmap)
setwd("/mnt/sds-hd/sd22b002/projects/ines/heart_revremod/")

# Selected genes
gene_df <- read.csv("results/bulk/mrna_time_clusters_example-genes.csv", header = TRUE)

# Load expression data
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")

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

# Normalise expression
vsd <- vst(dds, blind = TRUE)

# Extract genes of interest
norm_counts <- as.data.frame(assay(vsd))
norm_counts$EnsemblID <- rownames(norm_counts)
norm_counts_genes <- norm_counts %>% filter(EnsemblID %in% gene_df$EnsemblID)
norm_counts_genes <- norm_counts_genes %>% select(!EnsemblID) %>% as.matrix()

# Translate rownames from EnsemblID to Gene name for pretty visualisation
norm_counts_genes <- norm_counts_genes[gene_df$EnsemblID,]
rownames(norm_counts_genes) <- gene_df$Gene.name

# Reorder samples to follow "time-series" on heatmap
sample_order <- c("SO_2", "SO_5", "SO_52", "SO_53", "SDB_19", "SDB_20", "SDB_22", "SDB_24", "O_10", "O_12", "O_15", "O_17", "DB_28", "DB_29", "DB_32", "DB_35")
norm_counts_genes <- norm_counts_genes[, sample_order]

# Scale genes of interest
scaled_counts_genes <- t(scale(t(norm_counts_genes)))

# Plot heatmap
col_ha <- HeatmapAnnotation(Cluster = gene_df$Cluster,
                       col = list(Cluster = c("gene cluster 1" = "#66c2a5",
                                              "gene cluster 2" = "#fc8d62",
                                              "gene cluster 3" = "#8da0cb",
                                              "gene cluster 4" = "#e78ac3",
                                              "gene cluster 5" = "#a6d854",
                                              "gene cluster 6" = "#ffd92f")))
metadata <- metadata[sample_order,]
row_ha <- rowAnnotation(Condition = metadata$condition_simplified,
                            col = list(Condition = condition_simplified_colors))

pdf("results/bulk/mrna_time_clusters_example-genes_heatmap.pdf",
    height = 3, width = 11)
Heatmap(t(scaled_counts_genes),
        name = "scaled\nexpression",
        cluster_rows = TRUE,
        cluster_columns=FALSE,
        show_row_names = FALSE,
        column_split = gene_df$Cluster,
        col = circlize::colorRamp2(c(-2, 0, 2), c("blue", "white", "red")),
        left_annotation = row_ha,
        top_annotation = col_ha
        )
dev.off()