3  Results

3.1 Reproducibility

I list below the R packages required to reproduce the analyses.

## data wrangling
library(dplyr)
library(tidyr)
library(stringr)
library(purrr)

## plotting
library(ggplot2)
library(ComplexHeatmap)
library(cowplot)
library(grid)
library(RColorBrewer)
library(igraph)
library(magick)

## auxiliary functions
source("R/utils.R")
source("R/visualisations.R")
source("R/network_exploration.R")
# today_date <- "2025-04-28"
today_date <- "2025-07-01"
## today_date <- format(Sys.Date(), "%Y-%m-%d")

## set the seed, for fixing generation of ComplexHeatmaps (dendogram clustering)
set.seed(20)

Except stated otherwise, all correlation matrices have been computed using Pearson Correlation, provided by the stats::cor function.

3.2 Compute Drug-drug Similarities

3.2.1 Drug By Drug Heatmaps

• Heatmap in Figure 3.1 is generated at the replicate level:

# Read DGEList object
dge_normalised_replicate <- readRDS(file = paste0("./data/processed/dge_normalised_replicate_",
                          today_date, ".rds"))

filename <- paste0("figures/correlation-heatmaps/sample-by-sample/heatmap_normalised_replicates_", today_date)

heatmap_object <- compute_Heatmap(dge_normalised_replicate,
                                                annotation_cols = c("Pathway", "MoA", "Compound"),
                                                matrix_title = "Pearson Correlation between biological replicates",
                                                correlation = TRUE,
                                                filename_heatmap = filename)

knitr::include_graphics(paste0(filename,".png"),
                        auto_pdf = TRUE,
                        dpi = 100)

Figure 3.1: Heatmap of normalised replication scores.

• Heatmap in Figure 3.2 is generated by averaging normalised barcode counts over replicates:

## Read DGEList object
dge_normalised_samples <- readRDS(file = paste0("./data/processed/dge_normalised_samples_",
                          today_date, ".rds"))

filename <- paste0("figures/correlation-heatmaps/sample-by-sample/heatmap_normalised_samples_", today_date)

heatmap_object <- compute_Heatmap(dge_normalised_samples,
                                                annotation_cols = c("Pathway", "MoA", "Compound"),
                                                matrix_title = "Pearson Correlation between biological replicates",
                                                correlation = TRUE,
                                                filename_heatmap = filename)

knitr::include_graphics(paste0(filename,".png"),
                        auto_pdf = TRUE,
                        dpi = 300)

Figure 3.2: Heatmap of normalised compounds by batch(aggregated over replicates) with CC.

• Heatmap at the replicates level in Figure 3.3:

# Read DGEList object
dge_binarised_replicate <- readRDS(file = paste0("./results/differential_analyses/dge_discretised_replicate_",
                          today_date, ".rds"))

filename <- paste0("figures/correlation-heatmaps/sample-by-sample/heatmap_binarised_replicates_", today_date)

heatmap_object <- compute_Heatmap(dge_binarised_replicate,
                                                annotation_cols = c("Pathway", "MoA", "Compound"),
                                                matrix_title = "Pearson Correlation between biological replicates",
                                                correlation = TRUE,
                                                filename_heatmap = filename)

knitr::include_graphics(paste0(filename,".png"),
                        auto_pdf = TRUE,
                        dpi = 100)

Figure 3.3: Heatmaps at the replicates level.

• Heatmap at the compound by batch level in Figure 3.4 (aggregated over technical replicates):

## Read DGEList object
dge_binarised_samples <- readRDS(file = paste0("./results/differential_analyses/dge_discretised_samples_",
                          today_date, ".rds"))

filename <- paste0("figures/correlation-heatmaps/sample-by-sample/heatmap_binarised_samples_", today_date)

heatmap_object <- compute_Heatmap(dge_binarised_samples,
                                                annotation_cols = c("Pathway", "MoA", "Compound"),
                                                matrix_title = "Pearson Correlation between biological replicates",
                                                correlation = TRUE,
                                                filename_heatmap = filename)

knitr::include_graphics(paste0(filename,".png"),
                        auto_pdf = TRUE,
                        dpi = 300)                           

Figure 3.4: Heatmaps at the compound by batch level.

3.2.2 Drug-Drug Networks

• Plot weighted undirected compound graphs with igraph::graph_from_adjacency_matrix in Figure 3.5 at the replicate level, after binarisation:

filename_drug <- paste0("./figures/drug-drug networks/replicate-normalised-drug-drug_", today_date)
n_pairwise_comparisons <- (ncol(dge_normalised_replicate) * (ncol(dge_normalised_replicate) - 1))/2
igraph_normalised_replicates <- compute_drug_networks(dge_normalised_replicate,
                                       filename_drug = filename_drug,
                                       title_graph = "Drug by drug network, after normalisation",
                                       graph_date = today_date,
                                       pval_threshold = 0.05/n_pairwise_comparisons, # basic FWER approach
                                       cor_threshold = 0.8, 
                                       export_graph = FALSE)

knitr::include_graphics(paste0(filename_drug, ".png"),
                        dpi = 200, auto_pdf = TRUE)

Figure 3.5: Drug-drug network interactions, at the replicates level.

• Plot weighted undirected compound graphs with igraph::graph_from_adjacency_matrix in Figure 3.6, averaged over replicates, after normalisation:

filename_drug <- paste0("./figures/drug-drug networks/sample-normalised-drug-drug_", today_date)
n_pairwise_comparisons <- (ncol(dge_normalised_samples) * (ncol(dge_normalised_samples) - 1))/2
igraph_normalised_samples <- compute_drug_networks(dge_normalised_samples,
                                       filename_drug = filename_drug,
                                       title_graph = "Drug by drug network, after normalisation",
                                       graph_date = today_date,
                                       pval_threshold = 0.05/n_pairwise_comparisons, # basic FWER approach
                                       cor_threshold = 0.8)

knitr::include_graphics(paste0(filename_drug, ".png"),
                        dpi = 200, auto_pdf = TRUE)

Figure 3.6: Drug-drug network interactions, at the compound level.

• Plot weighted undirected compound graphs with igraph::graph_from_adjacency_matrix in Figure 3.8 at the raplicate leval, after binarisation:

filename_drug <- paste0("./figures/drug-drug networks/replicate-binarised-drug-drug_", today_date)
n_pairwise_comparisons <- (ncol(dge_binarised_replicate) * (ncol(dge_binarised_replicate) - 1))/2
igraph_binarised_replicates <- compute_drug_networks(dge_binarised_replicate,
                                       filename_drug = filename_drug,
                                       title_graph = "Drug by drug network, after binarisation",
                                       graph_date = today_date,
                                       pval_threshold = 0.05/n_pairwise_comparisons, # basic FWER approach
                                       cor_threshold = 0.4, 
                                       export_graph = FALSE)

knitr::include_graphics(paste0(filename_drug, ".png"),
                        dpi = 200, auto_pdf = TRUE)

Figure 3.7: Drug-drug network interactions, at the replicates level.

• Plot weighted undirected compound graphs with igraph::graph_from_adjacency_matrix in Figure 3.8, averaged over replicates, after binarisation:

filename_drug <- paste0("./figures/drug-drug networks/sample-binarised-drug-drug_", today_date)
n_pairwise_comparisons <- (ncol(dge_binarised_samples) * (ncol(dge_binarised_samples) - 1))/2
igraph_binarised_samples <- compute_drug_networks(dge_binarised_samples,
                                       filename_drug = filename_drug,
                                       title_graph = "Drug by drug network, after binarisation",
                                       graph_date = today_date,
                                       pval_threshold = 0.05/n_pairwise_comparisons, # basic FWER approach
                                       cor_threshold = 0.4)

knitr::include_graphics(paste0(filename_drug, ".png"),
                        dpi = 200, auto_pdf = TRUE)

Figure 3.8: Drug-drug network interactions, at the compound level.

3.2.3 Find top hits for a given drug target

# Returns all shared neighbours for Pemetrexed
binarised_graph <- igraph_binarised_samples$igraph_object
Pemetrexed_vertices <- igraph::V(binarised_graph)[grepl("Pemetre", name)]
neigh_list <- lapply(Pemetrexed_vertices, function(v)
  igraph::neighbors(binarised_graph, v = Pemetrexed_vertices))
do.call(igraph::intersection, neigh_list)
##  + 14/128 vertices, named, from 98b7c50:
##   [1] 5-Fluorouracil_exp070222_Conc: 6 uM_Days: 9_           
##   [2] Bortezomib_exp200921_Conc: 20 nM_Days: 9_              
##   [3] Bortezomib_exp271221_Conc: 20 nM_Days: 9_              
##   [4] Cisplatin_exp181021_Conc: 500 nM_Days: 9_              
##   [5] L-755507_exp220322_Conc: 25 uM_Days: 9_                
##   [6] MG132_exp271221_Conc: 510 nM_Days: 9_                  
##   [7] Pemetrexed_exp040821_Conc: 100 nM_Days: 9_             
##   [8] Pemetrexed_exp181021_Conc: 80 nM_Days: 9_              
##   [9] Pemetrexed_exp281022_time_course_Conc: 100 nM_Days: 19_
##  [10] Pemetrexed_exp281022_time_course_Conc: 100 nM_Days: 30_
##  + ... omitted several vertices

## Top hits normalised
cor_normalised_compound <- igraph_normalised_samples$cor_pearson
pval_normalised_compound <- igraph_normalised_samples$cor_pval
top_scores_X13271_normalised <- compute_top_hits(vertice_names = c("X13271_exp200921_Conc: 8 uM_Days: 9_", 
                                                             "X13271_exp271221_Conc: 9 uM_Days: 9_"), 
                                           cor_mat = cor_normalised_compound, 
                                           pval_mat = pval_normalised_compound)

## Top hits binarised
cor_binarised_compound <- igraph_binarised_samples$cor_pearson
pval_binarised_compound <- igraph_binarised_samples$cor_pval
# grep("X13271", colnames(cor_binarised_compound), value = TRUE)
top_scores_X13271_binarised <- compute_top_hits(vertice_names = c("X13271_exp200921_Conc: 8 uM_Days: 9_", 
                                                             "X13271_exp271221_Conc: 9 uM_Days: 9_"), 
                                           cor_mat = cor_binarised_compound, 
                                           pval_mat = pval_binarised_compound)

openxlsx::write.xlsx(list(normalised_X13271 = top_scores_X13271_normalised |> 
                            dplyr::filter(!grepl("^Contro", top_hits)), 
                          binarised_X13271 = top_scores_X13271_binarised), 
                     file = paste0("./results/top_hits/X13271_", today_date, ".xlsx"),
                     asTable = TRUE)

3.3 Compute Cell Lines Similarities

3.3.1 Cell Line by Drug Heatmap

dge_normalised_samples <- readRDS(file = paste0("./data/processed/dge_normalised_samples_",
                          today_date, ".rds"))
filename <- paste0("figures/correlation-heatmaps/barcode-by-sample/barcode_by_sample_normalised_", today_date)

heatmap_object <- compute_Heatmap(dge_normalised_samples,
                                  annotation_cols = c("Pathway", "MoA", "Compound"),
                                  matrix_title = "Barcode counts by samples, after normalisation",
                                  correlation = FALSE,
                                  filename_heatmap = filename, 
                                  size_tile = 0.3)

knitr::include_graphics(path = paste0(filename, ".png"),
                        auto_pdf = TRUE,
                        dpi = 50)                          

Figure 3.9: Heatmap cell lines by samples, averaging over replicates.

3.3.2 Cell Line By Cell Line Heatmap

threshold_cell_prevalence <- 9 ## Remove cells infrequent
selected_barcodes <- which(log1p(rowSums(dge_normalised_samples$counts)) >= threshold_cell_prevalence)
dge_normalised_samples_filtered <- dge_normalised_samples[selected_barcodes, ]

filename <- paste0("figures/correlation-heatmaps/barcode-by-barcode/heatmap_barcode_normalised_samples_", today_date)
heatmap_object <- compute_Heatmap(dge_normalised_samples_filtered,
                                  annotation_cols = c("Pathway", "MoA", "Compound"),
                                  matrix_title = "Barcode counts by samples, after normalisation",
                                  correlation = TRUE,
                                  transpose_mat = TRUE, 
                                  size_tile = 0.2,
                                  filename_heatmap = filename)

knitr::include_graphics(paste0(filename,".pdf"),
                        auto_pdf = TRUE,
                        dpi = 50)

Figure 3.10: Histogram of the distribution of cell resistance (for a given barcode ID, returns the total number of cell lines identified).