TIL density - XCR expansion

Setup

library(ithi.utils)
load_base_libs()

library(methods)
library(forcats)
library(grid)
library(gridExtra)
library(gridBase)

library(ithi.meta)
library(ithi.figures)
library(ithi.utils)
library(ithi.seq)
library(ithi.clones)
library(ithi.supp)
library(ithi.xcr)

ihc_table_path <- snakemake@input$ihc_table
xcr_table_path <- snakemake@input$xcr_table
tcr_diversity_file <- snakemake@input$tcr_diversity
bcr_diversity_file <- snakemake@input$bcr_diversity

db_path <- snakemake@params$db
total_tiltypes <- snakemake@params$total_tiltypes

annotation_colours <- ithi.figures::get_annotation_colours()

ihc_table <- fread(ihc_table_path)
xcr_table <- read_clonotypes(xcr_table_path, duplicates = FALSE, db_path = db_path)

Read 16.4% of 304822 rows
Read 49.2% of 304822 rows
Read 85.3% of 304822 rows
Read 304822 rows and 18 (of 18) columns from 0.070 GB file in 00:00:06

tcr_diversity <- read_xcr_diversity_file(tcr_diversity_file, db_path)
bcr_diversity <- read_xcr_diversity_file(bcr_diversity_file, db_path)
xcr_diversity <- ithi.supp::get_xcr_diversity(tcr_diversity_file, bcr_diversity_file, 
    db_path, xcr_table)

Analysis

This is a good point and one worth validating. As you’ll see, the results will be ‘interesting’ to this reviewer, though perhaps not to us since we’re aware that the TIL responses can be quite polyclonal.

XCR_RANK_COLS <- c("#BE2730", "#CA3339", "#E8BA40", "#EFE74C", "#90BC59", "#35A75C", 
    "#2AACED", "#2573B4", "#2A2B55", "#6B3E8F", "#D8D8D8")
names(XCR_RANK_COLS) <- paste0("Rank_", as.character(1:11))

tcr_top <- ithi.xcr::top_clonotype_table(xcr_table, xcrtype = "TRB")
bcr_top <- ithi.xcr::top_clonotype_table(xcr_table, xcrtype = "IGH")

tcr_mat <- reshape2::acast(tcr_top, formula = rank ~ condensed_id, value.var = "freq", 
    fun.aggregate = sum)
bcr_mat <- reshape2::acast(bcr_top, formula = rank ~ condensed_id, value.var = "freq", 
    fun.aggregate = sum)
tcr_mat <- tcr_mat[rev(1:nrow(tcr_mat)), ]
bcr_mat <- bcr_mat[rev(1:nrow(bcr_mat)), ]

tcr_mat_melt <- melt(as.matrix(tcr_mat)) %>% plyr::rename(c(Var1 = "rank", Var2 = "condensed_id", 
    value = "proportion"))
bcr_mat_melt <- melt(as.matrix(bcr_mat)) %>% plyr::rename(c(Var1 = "rank", Var2 = "condensed_id", 
    value = "proportion"))

# tcr_mat_fill <- fill_cols(tcr_mat, condensed_ids) bcr_mat_fill <-
# fill_cols(bcr_mat, condensed_ids)

ihc_table_subset <- subset(ihc_table, select = c("condensed_id", "patient_id", 
    total_tiltypes))
na_samples <- ihc_table_subset %>% subset(select = -c(condensed_id, patient_id)) %>% 
    apply(1, function(x) any(is.na(x)))
ihc_table_subset <- ihc_table_subset[!na_samples, ]

ihc_table_melted <- melt(ihc_table_subset, id.vars = c("condensed_id", "patient_id"), 
    measure.vars = total_tiltypes, variable.name = "tiltype", value.name = "density")
ihc_table_melted$celltype <- ifelse(ihc_table_melted$tiltype %in% c("T_CD8_density", 
    "T_CD4_density"), yes = "T_cell", no = "B_cell")

ids_intersect <- intersect(ihc_table_melted$condensed_id, tcr_mat_melt$condensed_id)
ihc_table_melted <- subset(ihc_table_melted, condensed_id %in% ids_intersect)
tcr_mat_melt <- subset(tcr_mat_melt, condensed_id %in% ids_intersect)
bcr_mat_melt <- subset(bcr_mat_melt, condensed_id %in% ids_intersect)

ihc_celltype_totals <- ihc_table_melted %>% group_by(condensed_id) %>% summarise(t_total = sum(density[celltype == 
    "T_cell"]), b_total = sum(density[celltype == "B_cell"]))
t_sample_order <- with(ihc_celltype_totals, condensed_id[order(t_total, decreasing = FALSE)])
b_sample_order <- with(ihc_celltype_totals, condensed_id[order(b_total, decreasing = FALSE)])

plot_stacked_xcr <- function(xcr_mat_melt, sample_order, label_x = FALSE) {
    xcr_mat_melt$condensed_id <- factor(xcr_mat_melt$condensed_id, levels = sample_order)
    p <- ggplot(xcr_mat_melt, aes(x = condensed_id, y = proportion)) + geom_bar(aes(fill = forcats::fct_rev(rank)), 
        stat = "identity") + theme_bw() + ithi.utils::theme_Publication() + 
        ithi.utils::theme_nature() + scale_fill_manual(values = XCR_RANK_COLS) + 
        scale_y_continuous(expand = c(0, 0)) + guides(fill = FALSE) + theme(axis.text.x = element_blank()) + 
        ithi.utils::ggmargins(type = "topminus")
    if (!label_x) {
        return(p + xlab(""))
    } else {
        return(p)
    }
}

p1 <- plot_stacked_xcr(tcr_mat_melt, sample_order = t_sample_order, label_x = FALSE)
p2 <- plot_stacked_xcr(bcr_mat_melt, sample_order = b_sample_order, label_x = FALSE)


celltypes <- c("T_cell", "B_cell")

create_tilplots <- function(ihc_table_melted, sample_order, selected_celltype) {
    dat <- subset(ihc_table_melted, celltype == selected_celltype)
    dat$condensed_id <- factor(dat$condensed_id, levels = sample_order)
    dat$patient_id <- ithi.meta::factor_id(dat$patient_id, type = "patient_id", 
        db_path)
    dat <- dat[order(dat$condensed_id), ]
    subtils <- dat$tiltype %>% unique
    subplots <- lapply(subtils, function(y) {
        p <- ggplot(dat %>% subset(tiltype == y), aes(x = condensed_id, y = density)) + 
            geom_bar(aes(fill = patient_id), stat = "identity") + theme_bw() + 
            theme_Publication() + theme_nature() + scale_y_continuous(expand = c(0, 
            0)) + scale_fill_manual(values = annotation_colours$patient_id) + 
            theme(axis.text.x = element_blank()) + xlab("") + ylab(y) + guides(fill = FALSE) + 
            ithi.utils::ggmargins(type = "topminus")
        return(p)
    })
    names(subplots) <- subtils
    return(subplots)
}

t_tilplots <- create_tilplots(ihc_table_melted, sample_order = t_sample_order, 
    selected_celltype = "T_cell")
b_tilplots <- create_tilplots(ihc_table_melted, sample_order = b_sample_order, 
    selected_celltype = "B_cell")
tilplots <- list(T_cell = t_tilplots, B_cell = b_tilplots)

tcr_diversity_table <- tcr_diversity %>% plyr::rename(c(observedDiversity_mean = "Unique clonotypes", 
    shannonWienerIndex_mean = "Entropy"))
bcr_diversity_table <- bcr_diversity %>% plyr::rename(c(observedDiversity_mean = "Unique clonotypes", 
    shannonWienerIndex_mean = "Entropy"))

plot_xcr_diversity <- function(xcr_diversity_table, sample_order, diversity_measures = c("Unique clonotypes", 
    "Entropy")) {
    xcr_diversity_table <- subset(xcr_diversity_table, condensed_id %in% sample_order)
    diversity_table_melted <- melt(xcr_diversity_table, id.vars = c("condensed_id", 
        "patient_id"), measure.vars = diversity_measures, variable.name = "measure", 
        value.name = "diversity")
    diversity_table_melted$condensed_id <- factor(diversity_table_melted$condensed_id, 
        levels = sample_order)
    subplots <- lapply(diversity_measures, function(divtype) {
        dat <- subset(diversity_table_melted, measure == divtype)
        dat$patient_id <- ithi.meta::factor_id(dat$patient_id, type = "patient_id", 
            db_path)
        
        ggplot(dat, aes(x = condensed_id, y = diversity)) + geom_bar(aes(fill = patient_id), 
            stat = "identity") + theme_bw() + theme_Publication() + theme_nature() + 
            scale_y_continuous(expand = c(0, 0)) + scale_fill_manual(values = annotation_colours$patient_id) + 
            theme(axis.text.x = element_blank()) + xlab("") + ylab(divtype) + 
            guides(fill = FALSE) + ithi.utils::ggmargins(type = "topminus")
    })
    names(subplots) <- diversity_measures
    return(subplots)
}

tcrdivplots <- plot_xcr_diversity(tcr_diversity_table, sample_order = t_sample_order)
bcrdivplots <- plot_xcr_diversity(bcr_diversity_table, sample_order = b_sample_order)

construct_combined_plot <- function(celltype) {
    if (celltype == "T_cell") {
        raw_gplots <- list(tilplots$T_cell$T_CD8_density, tilplots$T_cell$T_CD4_density, 
            p1, tcrdivplots$`Unique clonotypes`, tcrdivplots$Entropy)
    } else if (celltype == "B_cell") {
        raw_gplots <- list(tilplots$B_cell$T_CD20_density, tilplots$B_cell$T_Plasma_density, 
            p2, bcrdivplots$`Unique clonotypes`, bcrdivplots$Entropy)
    }
    
    grob_gplots <- lapply(raw_gplots, function(x) ggplotGrob(x))
    combined_grob <- do.call(gridExtra::rbind.gtable, grob_gplots)
    return(combined_grob)
}

t_combined_plot <- construct_combined_plot(celltype = "T_cell")
b_combined_plot <- construct_combined_plot(celltype = "B_cell")

TCR

grid.newpage()
grid.draw(t_combined_plot)

In general, we can see that higher CD8/CD4 T cell density is accompanied by an increase in the fraction of clonotype reads accounted for by rare clonotypes. In some samples/patients (e.g. patient 15) there appears to be ‘expansion’ of the most dominat clonotype(s), but if you look at the patient 15 sample with low TIL density (in the left half of the plot, near the middle) it becomes apparent that samples with lower T cell density have a higher proportion of clonotypes made up by rare clonotypes.

Therefore, it would be correct to say that the TIL responses in high-TIL density samples are polyclonal, rather than dominated by a single clone.

As you can see from the last two barplots, our diversity measures are capturing the observation that lower proportion accounted for by top 10 clonotypes = higher TCR diversity.

BCR

grid.newpage()
grid.draw(b_combined_plot)

The same isn’t as true for BCRs – we can see evidence that samples with the highest TIL densities are not necessarily polyclonal. This is in agreement with our observation that BCR clonality does not correlate with B-cell densities.

TODO: IDEALLY SHOW UNCERTAINTY IN XCR BARS TOO (using the estimates for TCR/BCR).

TODO: Relabel axis names, add colour legend for patients. Discuss if best way to show patients is with coloured bars or with a separate annotation track.

As you can see, it’s pretty obvious that the diversity measures are doing a good job capturing the distribution of top clonotypes.