Spatial analysis

library(ithi.utils)
load_base_libs()

library(ithi.meta)
library(ithi.ihc)
library(ithi.xcr)
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

spatial_results_file <- snakemake@input$spatial_results_file

tcr_diversity_file <- snakemake@input$tcr_diversity
bcr_diversity_file <- snakemake@input$bcr_diversity

nanostring_data_path <- snakemake@input$nanostring_data
nanostring_annotations_path <- snakemake@input$nanostring_annotations

clonal_measures_file <- snakemake@input$ith_stats

Metadata

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

id_cols <- c("voa", "patient_id", "condensed_id", "field", "tissue")
id_cols_sample <- c("voa", "patient_id", "condensed_id", "tissue")

Analysis

rename_field <- function(field) {
    voa_portion <- str_replace(field, "_HP_IM3.*$", "")
    voa_portion <- str_replace(voa_portion, "_", " ")
    field_portion <- str_extract(field, "_HP_IM3.*$")
    field_portion <- str_replace(field_portion, "__", "_[")
    field_portion <- paste0(field_portion, "]", sep = "")
    field <- paste(voa_portion, field_portion, sep = "")
    return(field)
}

fudge_cols <- function(ihc_slide_subset) {
    count_cols <- colnames(ihc_slide_subset)[str_detect(colnames(ihc_slide_subset), 
        "count")]
    ## JOIN WITH UNDERSCORES, THEN USE RESHAPE new_cols <-
    ## lapply(strsplit(str_replace_all(count_cols, '_count', ''), '_'),
    ## function(x) paste(rev(x), collapse='_')) colnames(ihc_slide_subset) <-
    ## mapvalues(colnames(ihc_slide_subset), count_cols, unlist(new_cols))
    df_melted <- melt(ihc_slide_subset, id.vars = c("voa", "field", "condensed_id"), 
        measure.vars = run1_cols, variable.name = "type", value.name = "count")
    df_melted$tissue <- mapvalues(str_extract(df_melted$type, "^(E|S)"), from = c("E", 
        "S"), to = c("epithelial", "stroma"))
    df_melted$type <- tolower(str_replace_all(df_melted$type, "(^(E|S)_|_count$)", 
        ""))
    df_cast <- dcast(df_melted, voa + field + condensed_id + tissue ~ type, 
        fun.aggregate = sum, value.var = "count")
    df_fudge <- df_cast %>% mutate(cd8_cd8 = cd8 * cd8, cd4_cd8 = cd4 * cd8, 
        cd4_cd4 = cd4 * cd4, cd20_cd4 = cd20 * cd4, cd20_cd8 = cd20 * cd8, cd20_cd20 = cd20 * 
            cd20)
    return(df_fudge)
}

spatial_data <- fread(spatial_results_file)

Read 0.0% of 1650000 rows
Read 12.1% of 1650000 rows
Read 22.4% of 1650000 rows
Read 34.5% of 1650000 rows
Read 46.7% of 1650000 rows
Read 58.2% of 1650000 rows
Read 69.7% of 1650000 rows
Read 80.6% of 1650000 rows
Read 91.5% of 1650000 rows
Read 1650000 rows and 28 (of 28) columns from 0.532 GB file in 00:00:15

if ("tissue.1" %in% colnames(spatial_data)) {
    spatial_data <- subset(spatial_data, select = -c(tissue.1))
}
spatial_data$field <- rename_field(spatial_data$field)

ihc_table_slide <- fread(ihc_table_slide_path)
ihc_table_slide$field <- str_replace(ihc_table_slide$slide, "\\.im3$", "")
ihc_table_slide <- subset(ihc_table_slide, select = -c(slide))

run1_tils <- c("CD8_count", "CD4_count", "CD20_count")
run1_cols <- c(paste0("E_", run1_tils), paste0("S_", run1_tils))
ihc_slide_subset <- subset(ihc_table_slide, select = c("voa", "field", "condensed_id", 
    run1_cols))

ihc_slide_fudge <- fudge_cols(ihc_slide_subset)

INTERACTION_VARS <- c("cd8_cd8", "cd4_cd8", "cd4_cd4", "cd20_cd4", "cd20_cd8", 
    "cd20_cd20")

ihc_slide_fudge_melted <- melt(ihc_slide_fudge, id.vars = c("voa", "field", 
    "condensed_id", "tissue"), measure.vars = INTERACTION_VARS, variable.name = "type", 
    value.name = "count")
spatial_data_melted <- melt(spatial_data, id.vars = c("voa", "field", "tissue", 
    "weight", "condensed_id", "patient_id"), measure.vars = INTERACTION_VARS, 
    variable.name = "type", value.name = "gamma")

spatial_data_merge <- merge(spatial_data_melted, ihc_slide_fudge_melted, by = c("voa", 
    "field", "tissue", "condensed_id", "type"))

THRESHOLD <- 10

spatial_data_merge <- subset(spatial_data_merge, count >= THRESHOLD)

spatial_summary_melted <- spatial_data_merge %>% group_by_(.dots = c(id_cols, 
    "type")) %>% summarise(gamma = weighted.mean(gamma, weight))

Remove super-outliers

spatial_summary_melted_filtered <- spatial_summary_melted %>% group_by(type, 
    tissue) %>% filter(!(abs(gamma - median(gamma, na.rm = TRUE)) > 10 * IQR(gamma, 
    na.rm = TRUE)))

Posterior distributions

ggplot(spatial_summary_melted_filtered, aes(x = gamma)) + geom_histogram(aes(fill = tissue), 
    bins = 30) + theme_bw() + theme_Publication() + facet_wrap(~type, scales = "free")

Correlations with TIL densities

Slide level

spatial_summary_filtered <- dcast(spatial_summary_melted_filtered, field + tissue ~ 
    type, value.var = "gamma")

spatial_ihc <- merge(ihc_table_slide, spatial_summary_filtered, by = c("field"))

IHC_VARS <- c("T_CD8_density", "T_CD4_density", "T_CD20_density", "T_Plasma_density", 
    "E_CD8_density", "E_CD4_density", "E_CD20_density", "E_Plasma_density", 
    "S_CD8_density", "S_CD4_density", "S_CD20_density", "S_Plasma_density", 
    "RUN1_E_area_pct", "RUN1_S_area_pct")

plot_correlation_heatmap <- function(dat, id_cols, VARS) {
    spatial_ihc_subset <- subset(dat, select = c(id_cols, VARS))
    spatial_ihc_subset_mat <- subset(spatial_ihc_subset, select = c(VARS))
    
    corres <- corr.test(spatial_ihc_subset_mat, adjust = "fdr", method = "spearman")
    corres.r <- corres$r
    corres.p <- corres$p
    hc <- hclust(dist(corres.r), method = "ward.D")
    
    diag(corres.r) <- NA
    diag(corres.p) <- NA
    
    corres.labels <- format(signif(corres.p, 2), 2)
    
    cols <- colorRampPalette(c("#95cbee", "#ffd73e", "#ce472e"))(100)
    pheatmap(corres.r[hc$order, hc$order], display_numbers = corres.labels, 
        color = cols, number_color = "black", cluster_rows = TRUE, cluster_cols = TRUE)
}

## Will fail due to filtering plot_correlation_heatmap(subset(spatial_ihc,
## tissue == 'epithelial'), id_cols, c(IHC_VARS, INTERACTION_VARS))

# plot_correlation_heatmap(subset(spatial_ihc, tissue == 'stroma'), id_cols,
# c(IHC_VARS, INTERACTION_VARS))

Sample level

ihc_table <- fread(ihc_table_path)

N_THRESHOLD <- 0

spatial_summary_melted_filtered_sample <- spatial_summary_melted_filtered %>% 
    group_by(voa, patient_id, condensed_id, tissue, type) %>% summarise(gamma = mean(gamma, 
    na.rm = TRUE), n = n())
spatial_summary_melted_filtered_sample <- subset(spatial_summary_melted_filtered_sample, 
    n >= N_THRESHOLD)

spatial_summary_filtered_sample <- dcast(spatial_summary_melted_filtered_sample, 
    voa + condensed_id + tissue ~ type, value.var = "gamma")

spatial_ihc_sample <- merge(ihc_table, spatial_summary_filtered_sample, by = c("voa", 
    "condensed_id"))

plot_correlation_heatmap(subset(spatial_ihc_sample, tissue == "epithelial"), 
    id_cols_sample, c(IHC_VARS, INTERACTION_VARS))

plot_correlation_heatmap(subset(spatial_ihc_sample, tissue == "stroma"), id_cols_sample, 
    c(IHC_VARS, INTERACTION_VARS))

Correlations with ITH

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

ITH_VARS <- c("entropy", "divergence", "postprocessed_divergence", "combined_ith_raw", 
    "combined_ith_normalized", "proportion_subclonal")

spatial_ith <- merge(clonal_measures, spatial_summary_filtered_sample, by = c("condensed_id"))

plot_correlation_heatmap(subset(spatial_ith, tissue == "epithelial"), id_cols_sample, 
    c(ITH_VARS, INTERACTION_VARS))

plot_correlation_heatmap(subset(spatial_ith, tissue == "stroma"), id_cols_sample, 
    c(ITH_VARS, INTERACTION_VARS))

spatial_ith_melted <- melt(spatial_ith, id.vars = colnames(spatial_ith)[!colnames(spatial_ith) %in% 
    ITH_VARS], measure.vars = ITH_VARS, variable.name = "ithtype", value.name = "ith")

spatial_ith_melted$patient_id <- factor(spatial_ith_melted$patient_id)

pvals <- setNames(ddply(spatial_ith_melted, .(tissue, ithtype), function(x) {
    df <- x
    corres <- with(df, cor.test(cd8_cd8, ith, method = "spearman"))
    
    pval <- corres$p.value
    eq <- substitute(italic(P) == p, list(p = format(pval, digits = 3)))
    return(as.character(as.expression(eq)))
}), c("tissue", "ithtype", "p.value"))

ggplot(spatial_ith_melted, aes(x = cd8_cd8, y = ith)) + geom_point(aes(colour = patient_id)) + 
    scale_colour_manual(values = pal_patient) + theme_bw() + theme_Publication() + 
    facet_wrap(tissue ~ ithtype, scales = "free") + geom_text(data = pvals, 
    aes(x = Inf, y = Inf, label = p.value), hjust = 1.1, vjust = 1.5, size = 3, 
    parse = TRUE)