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)

# Tabular
library(flextable)
library(reactable)
library(htmltools)
html_bool <- knitr::is_html_output()

## 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/differential_barcoding.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 at the replicates level in Figure 3.1:

# 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.1: Heatmaps at the replicates level.

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

## Read DGEList object
dge_binarised_samples <- readRDS(file = paste0("./results/differential_analyses/dge_discretised_samples_",
                          today_date, ".rds")) |> 
  trim_constant_values()
filename <- paste0("figures/correlation-heatmaps/sample-by-sample/heatmap_binarised_samples_", today_date)

## Generate corresponding ComplexHeatmap
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.2: 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.4 at the replicate level, after binarisation:

filename_drug_plot <- paste0("./figures/drug-drug networks/replicate-binarised-drug-drug_CC_0.2_", today_date)
filename_drug_data <- paste0("./results/drug-drug-similarity-scores/replicate-binarised-drug-drug_CC_0.2_", today_date)
n_pairwise_comparisons <- (ncol(dge_binarised_replicate) * (ncol(dge_binarised_replicate) - 1))/2

# setting a CC threshold of 0.2
igraph_binarised_replicates_CC_0.2 <- compute_drug_networks(dge_binarised_replicate,
                                                     filename_drug_network_plot = filename_drug_plot, 
                                                     filename_drug_network_data = filename_drug_data,
                                                     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.2, 
                                                     export_graph = FALSE)
knitr::include_graphics(paste0(filename_drug_plot, ".png"),
                        dpi = 200, auto_pdf = TRUE)

# setting a CC threshold of 0.4
filename_drug_plot <- paste0("./figures/drug-drug networks/replicate-binarised-drug-drug_CC_0.4_", today_date)
filename_drug_data <- paste0("./results/drug-drug-similarity-scores/replicate-binarised-drug-drug_CC_0.4_", today_date)
igraph_binarised_replicates_CC_0.4 <- compute_drug_networks(dge_binarised_replicate,
                                                            filename_drug_network_plot = filename_drug_plot, 
                                                            filename_drug_network_data = filename_drug_data,
                                                            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_plot, ".png"),
                        dpi = 200, auto_pdf = TRUE)

(a) With a CC of 0.2

(b) With a CC of 0.4

Figure 3.3: Drug-drug network, after binarisation, and at the replicate level.

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

filename_drug_plot <- paste0("./figures/drug-drug networks/sample-binarised-drug-drug_CC_0.2_", today_date)
filename_drug_data <- paste0("./results/drug-drug-similarity-scores/sample-binarised-drug-drug_CC_0.2_", 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_network_plot =  filename_drug_plot,
                                                  filename_drug_network_data = filename_drug_data,
                                                  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.2)

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

filename_drug_plot <- paste0("./figures/drug-drug networks/sample-binarised-drug-drug_CC_0.4_", today_date)
filename_drug_data <- paste0("./results/drug-drug-similarity-scores/sample-binarised-drug-drug_CC_0.4_", 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_network_plot =  filename_drug_plot,
                                                  filename_drug_network_data = filename_drug_data,
                                                  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_plot, ".png"),
                        dpi = 200, auto_pdf = TRUE)

(a) With a CC of 0.2

(b) With a CC of 0.4

Figure 3.4: Drug-drug network interactions, at the compound level, and after binarisation.

3.2.3 Find top hits

# Returns all shared neighbours for Pemetrexed, in a network-centric approach
# 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)


## Top hits binarised
cor_binarised_compound <- igraph_binarised_samples$cor_pearson
pval_binarised_compound <- igraph_binarised_samples$cor_pval

pval_long <- pval_binarised_compound |> 
  tibble::as_tibble(rownames = "drug_one") |> 
  tidyr::pivot_longer(-drug_one,
                      names_to = "drug_two",
                      values_to =  "adj_pval" ) |> 
  dplyr::filter(!is.na(adj_pval))

cor_long <- cor_binarised_compound |> 
  tibble::as_tibble(rownames = "drug_one") |> 
  tidyr::pivot_longer(-drug_one,
                      names_to = "drug_two",
                      values_to =  "cor_pearson" ) |> 
  dplyr::filter(cor_pearson != 1)

top_hits <- pval_long |> 
  dplyr::inner_join(cor_long, 
                    by = join_by(drug_one, drug_two)) |> 
  dplyr::arrange(drug_one, desc(cor_pearson), adj_pval) 

# 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(binarised_top_hits = top_hits),
                     file = paste0("./results/top_hits/top_hits_", today_date, ".xlsx"),
                     asTable = TRUE)

top_hits_reactable <- reactable(top_hits,
                                columns = list(
                                  cor_pearson = colDef(
                                    minWidth = 200,  
                                    cell = function(value) {
                                      width <- paste0(max(0, value * 100 / max(top_hits$cor_pearson)), "%")
                                      bar <- div(
                                        class = "bar-chart",
                                        div(class = "bar", style = list(width = width, backgroundColor = "#3fc1c9"))
                                      )
                                      div(class = "bar-cell", span(class = "number", value), bar)
                                    }
                                  )
                                ), 
                                groupBy = "drug_one",
                                paginateSubRows = TRUE,
                                defaultColDef = colDef(
                                  cell = function(value) format(value, nsmall = 1),
                                  align = "center",
                                  headerClass = "header",
                                  minWidth = 150
                                ),
                                highlight = TRUE,
                                filterable = TRUE,
                                searchable = TRUE,
                                resizable = TRUE,
                                # to deal with the large number of observations
                                defaultPageSize = 25,
                                elementId = "top-hits-binarised-download",
                                class = "barcoding-tbl",
                                rowStyle = list(cursor = "pointer"))

top_hits_reactable <- div(
  class = "barcoding-table",
  h2(class = "barcoding-title", "Top hits for binarised values"),
  top_hits_reactable
)

htmltools::browsable(
  tagList(
    tags$button(
      tagList(fontawesome::fa("download"), "Download as CSV"),
      onclick = "Reactable.downloadDataCSV('top-hits-binarised-download', 'top-hits-binarised.csv')"
    ),
    top_hits_reactable
  )
)

Top hits for binarised values

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.5: 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.6: Histogram of the distribution of cell resistance (for a given barcode ID, returns the total number of cell lines identified).