ITH-TIL analysis

library(ithi.utils)
load_base_libs()

library(ithi.meta)
library(ithi.ihc)
library(ithi.clones)

Colour palettes

pal_patient <- select_palette("patient")

Parameters

db_path <- snakemake@params$db

ihc_table_path <- snakemake@input$ihc_table
ihc_table_slide_path <- snakemake@input$ihc_table_slide
tils_for_cluster <- unlist(snakemake@params$tils_for_cluster)
nclusts <- as.numeric(snakemake@params$nclusts)

tils_for_variability <- unlist(snakemake@params$tils_for_variability)
variability_type <- snakemake@params$variability_type

# clone_tree_dir <- snakemake@input$clone_tree_dir
clone_tree_file <- snakemake@input$clone_tree_file
clone_prevalence_file <- snakemake@input$clone_prevalence_file
clone_branch_length_file <- snakemake@input$clone_branch_length_file
clonal_measures_file <- snakemake@input$ith_stats
prevalence_threshold <- as.numeric(snakemake@params$prevalence_threshold)

proportion_subclonality_file <- snakemake@input$subclonality

Metadata

id_type <- "condensed_id"

db <- src_sqlite(db_path, create = FALSE)
duplicates <- collect(tbl(db, "duplicates"))

TIL-based clusters

We’ll use the following TIL subsets for clustering E_CD8_density, E_CD4_density, E_CD20_density, E_Plasma_density, S_CD8_density, S_CD4_density, S_CD20_density, S_Plasma_density.

ihc_table <- fread(ihc_table_path)

ihc_data <- subset(ihc_table, select = c(tils_for_cluster))
valid_rows <- apply(ihc_data, 1, function(x) all(!is.na(x)))

patient_ids <- as.data.frame(subset(ihc_table, select = c(patient_id)))
ihc_data <- as.data.frame(ihc_data[valid_rows, ])

patient_ids <- patient_ids[valid_rows, , drop = FALSE]
patient_ids$patient_id <- factor(patient_ids$patient_id)
colnames(patient_ids) <- "Patient"

rownames(ihc_data) <- rownames(patient_ids) <- as.data.frame(ihc_table)[, id_type][valid_rows]

dat_processed <- t(clip_values(scale(ihc_data)))
rownames(dat_processed) <- mapvalues(rownames(dat_processed), from = c("E_CD8_density", 
    "E_CD4_density", "E_CD20_density", "E_Plasma_density", "S_CD8_density", 
    "S_CD4_density", "S_CD20_density", "S_Plasma_density"), to = c("E CD8+", 
    "E CD4+", "E CD20+", "E Plasma", "S CD8+", "S CD4+", "S CD20+", "S Plasma"))

p <- pheatmap(dat_processed, clustering_method = "ward.D", annotation_col = patient_ids, 
    annotation_colors = list(Patient = pal_patient), cluster_rows = FALSE, show_colnames = FALSE, 
    legend = TRUE, cellwidth = 10, cellheight = 10, annotation_legend = FALSE)

clusts <- data.frame(cutree(p$tree_col, nclusts))
setDT(clusts, keep.rownames = TRUE)
clusts <- setNames(clusts, c(id_type, "cluster"))

clusts$cluster <- factor(clusts$cluster)
ihc_table_annotated <- merge(ihc_table, clusts, by = c(id_type))

Clonal diversity

tree_branch_data <- read_clone_tree_data(clone_tree_file, clone_branch_length_file, 
    clone_prevalence_file, db_path)
trees <- lapply(tree_branch_data, function(x) x$tree)
prevalences <- rbind.fill(lapply(tree_branch_data, function(x) x$prevalence_dat))
branch_lengths <- rbind.fill(lapply(tree_branch_data, function(x) x$branch_dat))

ccfs <- compute_ccf(prevalences, trees, id_type = id_type)
ccfs_labeled <- merge(ccfs, branch_lengths, by = c("label", "node", "patient_id"))

clonal_measures <- read_ith_stats(clonal_measures_file, db_path, duplicates = FALSE)

prevalence_indices <- prevalences %>% group_by(condensed_id, patient_id) %>% 
    summarise(shannon = vegan::diversity(clonal_prevalence, index = "shannon"), 
        simpson = vegan::diversity(clonal_prevalence, index = "simpson"))

prevalences_filtered <- subset(prevalences, clonal_prevalence >= prevalence_threshold)
# vafs_filtered <- subset(vafs, prevalence >= prevalence_threshold) ## Not
# exactly right but close enough
divergence <- compute_clone_divergence(prevalences_filtered, trees, id_type)
# divergence_weighted <- compute_clone_divergence(prevalences_filtered,
# vafs_filtered, clone_trees, branch_lengths, id_type, weighted=TRUE)
# colnames(divergence_weighted) <- mapvalues(colnames(divergence_weighted),
# from = 'clone_divergence', to = 'clone_divergence_weighted')
diversity <- Reduce(function(x, y) merge(x, y, by = c(id_type, "patient_id"), 
    all = TRUE), list(prevalence_indices, clonal_measures))

ITH-TIL correlations

Old measures

We’ll use proportion_subclonality as the surrogate.

proportion_subclonal <- fread(proportion_subclonality_file)

res <- plyr::join(clusts, proportion_subclonal)

if (nclusts == 2) {
    corres <- wilcox.test(proportion_subclonal ~ cluster, data = res)
} else {
    corres <- kruskal.test(proportion_subclonal ~ cluster, data = res)
}

pval <- corres$p.value
eq <- substitute(italic(P) == p, list(p = format(pval, digits = 3)))

ggplot(res, aes(x = cluster, y = proportion_subclonal)) + geom_boxplot(aes(fill = cluster)) + 
    theme_bw() + theme_Publication() + scale_fill_Publication() + annotate("text", 
    x = Inf, y = Inf, parse = TRUE, hjust = 1, vjust = 1, label = as.character(as.expression(eq)))

By a Kruskal-Wallis test.

If we instead choose to compare cluster 3 (high epithelial TILs) against everything else, then:

wilcox.test(proportion_subclonal ~ cluster == 3, data = res)

    Wilcoxon rank sum test

data:  proportion_subclonal by cluster == 3
W = 424, p-value = 0.005338
alternative hypothesis: true location shift is not equal to 0

New measures (in progress)

res <- merge(diversity, ihc_table_annotated, by = c(id_type, "patient_id"))

ith_measures <- c("shannon", "entropy", "simpson", "divergence", "postprocessed_divergence", 
    "proportion_subclonal", "combined_ith_normalized", "combined_ith_raw")
res_melted <- reshape2::melt(res, id.vars = c(id_type, "patient_id", "cluster"), 
    measure.vars = ith_measures, variable.name = "ith_type", value.name = "ith")
res_melted$ith_type <- factor(res_melted$ith_type)

stats <- setNames(ddply(res_melted, .(ith_type), function(x) {
    if (nclusts == 2) {
        res <- wilcox.test(ith ~ cluster, data = x)
    } else {
        res <- kruskal.test(ith ~ cluster, data = x)
    }
    pvalue <- res$p.value
    eq <- substitute(italic(P) == p, list(p = format(pvalue, digits = 3)))
    return(as.character(as.expression(eq)))
}), c("ith_type", "p.value"))

ggplot(res_melted, aes(x = cluster, y = ith)) + geom_boxplot(aes(fill = cluster)) + 
    theme_bw() + theme_Publication() + scale_fill_Publication() + facet_wrap(~ith_type, 
    scales = "free") + geom_text(data = stats, aes(x = Inf, y = Inf, label = p.value), 
    hjust = 1.1, vjust = 1.5, size = 3, parse = TRUE)

Check with Andrew as to how he computed clone divergence – my computations based on what I believe the definition is are different. Perhaps it’s a 1-minus issue?

ITH-TIL variability

ihc_table_slide <- fread(ihc_table_slide_path)

ihc_slide_data <- subset(ihc_table_slide, select = c("condensed_id", "patient_id", 
    "slide", tils_for_variability))
ihc_melted <- melt(ihc_slide_data, id.vars = c("condensed_id", "patient_id", 
    "slide"), measure.vars = tils_for_variability, variable.name = "tiltype", 
    value.name = "density")
ihc_melted$tiltype <- as.character(ihc_melted$tiltype)

intersect_samples <- intersect(unique(ihc_melted$condensed_id), unique(ccfs_labeled$condensed_id))

til_var <- compute_ihc_var(ihc_melted, ihc_table, samples = intersect_samples, 
    til_variance_option = variability_type, tils_for_variability)

clonedists <- clone_distances(ccfs_labeled, normalize = FALSE, id_type = id_type)
clonedists$patient_id <- as.character(clonedists$patient_id)

variabilities <- inner_join(til_var, subset(clonedists, select = -c(dist_clones_weighted)), 
    by = c("patient_id"))
variabilities_melted <- melt(variabilities, id.vars = c("patient_id", "mean_weighted_dist"), 
    measure.vars = tils_for_variability, variable.name = "tiltype", value.name = "tilvar")
variabilities_melted$patient_id <- factor(variabilities_melted$patient_id)

stats <- setNames(ddply(variabilities_melted, .(tiltype), function(x) {
    res <- with(x, cor.test(mean_weighted_dist, tilvar, method = "spearman"))
    pvalue <- res$p.value
    eq <- substitute(italic(P) == p, list(p = format(pvalue, digits = 3)))
    return(as.character(as.expression(eq)))
}), c("tiltype", "p.value"))

p <- ggplot(variabilities_melted, aes(x = mean_weighted_dist, y = tilvar)) + 
    geom_point(aes(colour = patient_id)) + theme_bw() + theme_Publication()
p <- p + facet_wrap(~tiltype, scales = "free", nrow = 1)
p <- p + theme(axis.text.x = element_text(size = 6, family = "Helvetica"), legend.text = element_text(size = 6, 
    family = "Helvetica"), axis.text.y = element_text(size = 6, family = "Helvetica"), 
    axis.title = element_text(size = 8, family = "Helvetica"), strip.text = element_text(size = 7, 
        family = "Helvetica"), strip.background = element_rect(fill = "white"))
p <- p + ylab("Variability") + xlab("Mean clonal distance")
p <- p + geom_text(data = stats, aes(x = Inf, y = Inf, label = p.value), hjust = 1.1, 
    vjust = 1.5, size = 3, parse = TRUE)
p <- p + scale_colour_manual(values = pal_patient) + guides(colour = guide_legend(title = "Patient"))

plot(p)