BCR Phylogenetics: Correlations

library(ithi.utils)
load_base_libs()

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

Colour palettes

pal_patient <- select_palette("patient")

Parameters

db_path <- snakemake@params$db

igpartition_outdir <- snakemake@input$igpartition_outdir
# igpartition_outdir <-
# '~/shahlab/pipeline_outputs/ith_immune/igpartition/run22'
ihc_table_file <- snakemake@input$ihc_table
tiltypes <- unlist(snakemake@params$tiltypes)

ith_stats_file <- snakemake@input$ith_stats_file
proportion_subclonality_file <- snakemake@input$subclonality

Description

This document describes the correlative analysis conducted using measures summarized from BCR phylogenetic data.

All measures are computed on 100 evenly-spaced samples from the posterior trace for each MCMC run. 2 MCMC runs are conducted for each BCR clonal family. 15 clonal families (randomly selected based on a reasonable set of criteria, listed below) are selected from each patient. Each clonal family can contain BCRs from multiple tumour sites.

Clonal family selection

Selection criteria:

• Estimated size of 300-500 BCRs (I say estimated, because I don’t know the exact size of a clonal family prior to the final, seeded clustering step.)

Size estimation is done by taking the clonal family size from subsampling and naive (vsearch) clustering, and multiplying it by the ratio of total library size (for all samples corresponding to the patient in question) to downsampled library size. Since subsampling is done in a stratified manner, this is correct.

Note that even though 15 clonal families were run per patient, several of them may have failed somewhere along the pipeline. This can be due to: poor estimation of family sizes (families larger than 1200 BCRs are not run for computation time issues, and empty families yield no results), cluster failures (will have to rerun a few because of this). As far as I’m aware, all easily-resolved checks/fixes have been implemented.

If we absolutely need 15 per patient, I suggest rerunning more families. This can be easily done but takes time (~ 1 week per complete run).

Definitions

For the analysis below, there will often be two sets of plots. These are named (open the code chunks to reveal):

• Mugration: Full phylogeographic model
• Nongeographic: Normal phylogenetic model with no geographic information

There are arguments to why using mugration wouldn’t be ‘fair’. Obviously, incorporating geography as a discrete trait naturally leads to clustering by geographic location – this might introduce biases into the genealogical sorting analysis below. Hence why I also do the same analysis for the nongeographic trees.

Data

ihc_table <- fread(ihc_table_file)
ihc_subset <- subset(ihc_table, select = c("condensed_id", "patient_id", tiltypes))
ihc_means <- ihc_subset %>% group_by(patient_id) %>% select(-c(condensed_id)) %>% 
    summarise_all(mean)

ith_stats <- read_ith_stats(ith_stats_file, db_path, duplicates = FALSE)
proportion_subclonality <- fread(proportion_subclonality_file)

mugration_outdir <- file.path(igpartition_outdir, "stats/all/phylogeography/bayesian_mugration_standard_nonclock/1-end/partitioned_GTR")
nongeographic_outdir <- file.path(igpartition_outdir, "stats/all/phylogenetics/bayesian_nongeographic/1-end/partitioned")
rootfreq_outdir <- file.path(igpartition_outdir, "plots/phylogeography/bayesian_mugration_standard_nonclock/1-end/partitioned_GTR")

mugration_statfiles <- list.files(mugration_outdir, pattern = ".tsv", recursive = TRUE, 
    full.names = TRUE)
nongeographic_statfiles <- list.files(nongeographic_outdir, pattern = ".tsv", 
    recursive = TRUE, full.names = TRUE)

rootfreq_files <- list.files(rootfreq_outdir, pattern = ".txt", recursive = TRUE, 
    full.names = TRUE)

Genealogical sorting analysis

read_gsi_stats <- function(f) {
    obskey <- "observed"
    dat <- fread(f)
    dat <- subset(dat, median_gsi != -1)
    
    obs <- subset(dat, permutation == obskey)
    bkg <- subset(dat, permutation != obskey)
    
    bkg_mean <- bkg %>% group_by(group) %>% summarise(gsi_m = median(median_gsi, 
        na.rm = TRUE), gsi_s = mad(median_gsi, na.rm = TRUE))
    empirical <- dat %>% group_by(group) %>% summarise(pval = 1 - length(which(median_gsi[which(permutation != 
        obskey)] < median_gsi[permutation == obskey]))/(n() - 1))
    ## Not really Gaussian but visually looks ~close
    dat_summary <- merge(obs, bkg_mean, by = c("group"))
    dat_summary$z_score <- with(dat_summary, (median_gsi - gsi_m)/gsi_s)
    
    res <- merge(dat_summary, empirical, by = c("group"))
    res <- subset(res, select = -c(permutation))
    res <- subset(res, gsi_s != 0)
    
    return(res)
}

get_gsi_table <- function(statfiles) {
    gsi_stats <- rbind.fill(lapply(statfiles[str_detect(statfiles, "gsi_values")], 
        function(f) {
            df <- read_gsi_stats(f)
            
            patient_clust_rep <- str_extract(f, "[0-9]+\\/clust[0-9]+/[0-9]+")
            strs <- strsplit(patient_clust_rep, "/")[[1]]
            patient <- strs[1]
            clust <- strs[2]
            rep <- strs[3]
            
            if (nrow(df) != 0) {
                df <- cbind(patient_id = patient, clust = clust, rep = rep, 
                  df)
                df$condensed_id <- with(df, paste(patient_id, group, sep = "_"))
            }
            return(df)
        }))
    return(gsi_stats)
}

mugration_gsi_stats <- get_gsi_table(mugration_statfiles)

nongeographic_gsi_stats <- get_gsi_table(nongeographic_statfiles)

summarize_gsi_stats <- function(gsi_stats) {
    gsi_stats_avg <- gsi_stats %>% group_by(patient_id, clust, group, condensed_id) %>% 
        select(-c(rep)) %>% summarise_all(mean)
    gsi_summary <- gsi_stats_avg %>% group_by(condensed_id, patient_id, group) %>% 
        summarise(medgsi = median(median_gsi, na.rm = TRUE), z = median(z_score, 
            na.rm = TRUE), psig = length(which(pval < 0.05))/length(which(median_gsi != 
            -1)), pval = median(pval[which(median_gsi != -1)], na.rm = TRUE), 
            num = length(which(median_gsi != -1)))
    
    gsi_summary <- subset(gsi_summary, !is.na(z))
    return(gsi_summary)
}

gsi_ihc_correlation <- function(gsi_summary, gsi_measure = "medgsi", ihc, tiltypes) {
    gsi_ihc <- merge(gsi_summary, ihc, by = c("condensed_id", "patient_id"))
    gsi_ihc_melted <- melt(gsi_ihc, id.vars = c("condensed_id", "patient_id", 
        "medgsi", "z", "psig", "pval", "num"), measure.vars = tiltypes, variable.name = "tiltype", 
        value.name = "density")
    
    ## Filter for number of observations
    gsi_ihc_melted_filtered <- subset(gsi_ihc_melted, num >= 3)
    
    pvals <- setNames(ddply(gsi_ihc_melted_filtered, .(tiltype), function(x) {
        df <- as.data.frame(x)
        corres <- cor.test(df[, "density"], df[, gsi_measure], method = "spearman")
        
        pval <- corres$p.value
        eq <- substitute(italic(P) == p, list(p = format(pval, digits = 3)))
        return(as.character(as.expression(eq)))
    }), c("tiltype", "p.value"))
    
    p <- ggplot(gsi_ihc_melted_filtered, aes_string(x = "density", y = gsi_measure)) + 
        geom_point(aes(colour = patient_id)) + theme_bw() + theme_Publication() + 
        facet_wrap(~tiltype, scales = "free") + scale_color_manual(values = pal_patient) + 
        geom_text(data = pvals, aes(x = Inf, y = Inf, label = p.value), hjust = 1.1, 
            vjust = 1.5, size = 3, parse = TRUE)
    return(p)
}

gsi_ith_correlation <- function(gsi_summary, gsi_measure = "medgsi", ith, ith_measures) {
    gsi_ith <- merge(gsi_summary, ith, by = c("condensed_id", "patient_id"))
    gsi_ith_melted <- melt(gsi_ith, id.vars = c("condensed_id", "patient_id", 
        "medgsi", "z", "psig", "pval", "num"), measure.vars = ith_measures, 
        variable.name = "ithtype", value.name = "ith")
    
    ## Filter for number of observations
    gsi_ith_melted_filtered <- subset(gsi_ith_melted, num >= 3)
    
    pvals <- setNames(ddply(gsi_ith_melted_filtered, .(ithtype), function(x) {
        df <- as.data.frame(x)
        corres <- cor.test(df[, "ith"], df[, gsi_measure], method = "spearman")
        
        pval <- corres$p.value
        eq <- substitute(italic(P) == p, list(p = format(pval, digits = 3)))
        return(as.character(as.expression(eq)))
    }), c("ithtype", "p.value"))
    
    p <- ggplot(gsi_ith_melted_filtered, aes_string(x = "ith", y = gsi_measure)) + 
        geom_point(aes(colour = patient_id)) + theme_bw() + theme_Publication() + 
        facet_wrap(~ithtype, scales = "free") + scale_color_manual(values = pal_patient) + 
        geom_text(data = pvals, aes(x = Inf, y = Inf, label = p.value), hjust = 1.1, 
            vjust = 1.5, size = 3, parse = TRUE)
    return(p)
}

mugration_gsi_summary <- summarize_gsi_stats(mugration_gsi_stats)
nongeographic_gsi_summary <- summarize_gsi_stats(nongeographic_gsi_stats)

GSI-IHC analysis

The results below show that samples with higher TIL densities have lower GSIs. Lower GSIs = less geographic exclusivity in phylogenies, i.e. less phylogenetic signal in geography, i.e. more mixing.

p <- gsi_ihc_correlation(mugration_gsi_summary, "medgsi", ihc_subset, tiltypes)
plot(p)

p <- gsi_ihc_correlation(nongeographic_gsi_summary, "medgsi", ihc_subset, tiltypes)
plot(p)

Note that this result is not significant when stratifying by patient (i.e. modeling patient as a random effect).

GSI-ITH analysis

As with all other analyses incorporating ITH this analysis is currently flawed as the ITH statistics are off.

p <- gsi_ith_correlation(mugration_gsi_summary, "medgsi", ith_stats, c("entropy", 
    "divergence"))
plot(p)

p <- gsi_ith_correlation(nongeographic_gsi_summary, "medgsi", ith_stats, c("entropy", 
    "divergence"))
plot(p)

However, to inspire some potentially misplaced hope, here’s what the analysis looks like using the naive measure ‘proportion subclonality’, which correlated very well with the old ITH measures:

proportion_subclonality <- subset(proportion_subclonality, patient_id != 15)

p <- gsi_ith_correlation(mugration_gsi_summary, "medgsi", proportion_subclonality, 
    c("proportion_subclonal", "mix"))
plot(p)

p <- gsi_ith_correlation(nongeographic_gsi_summary, "medgsi", proportion_subclonality, 
    c("proportion_subclonal", "mix"))
plot(p)

Note that patient 15 was deliberately excluded because it was the sample for which CN calls were bad (until Andrew reruns the proportion subclonality calculation I shouldn’t include it).

If this finding holds up, it would imply that more intratumoural heterogeneity = more BCR mixing within clonal families. In other words, ITH <=> heterogeneity within BCR families! Pretty cool! Beyond that, I’m not sure how to biologically interpret this.

Phylogenetic signal in readcounts

read_phylosignal_file <- function(f) {
    dat <- fread(f)
    cbind(var = c("mut_freqs", "readcount"), dat)
}

get_phylosignal_table <- function(statfiles) {
    phylosignal_stats <- rbind.fill(lapply(statfiles[str_detect(statfiles, "phylosignal")], 
        function(f) {
            df <- read_phylosignal_file(f)
            
            patient_clust_rep <- str_extract(f, "[0-9]+\\/clust[0-9]+/[0-9]+")
            strs <- strsplit(patient_clust_rep, "/")[[1]]
            patient <- strs[1]
            clust <- strs[2]
            rep <- strs[3]
            
            if (nrow(df) != 0) {
                df <- cbind(patient_id = patient, clust = clust, rep = rep, 
                  df)
            }
            return(df)
        }))
    return(phylosignal_stats)
}

summarize_phylosignal_stats <- function(phylosignal_stats) {
    phylosignal_stats_avg <- phylosignal_stats %>% group_by(patient_id, clust, 
        var) %>% select(-c(rep)) %>% summarise_all(mean)
    phylosignal_summary <- phylosignal_stats_avg %>% group_by(patient_id, var) %>% 
        summarise(Cmean = median(stat.Cmean, na.rm = TRUE), Cmean_p = median(pvalue.Cmean, 
            na.rm = TRUE), I = median(stat.I, na.rm = TRUE), I_p = median(pvalue.I, 
            na.rm = TRUE), num = n())
    
    return(phylosignal_summary)
}

phylosignal_ihc_correlation <- function(phylosignal_summary, phylosignal_measure = "Cmean", 
    ihc, tiltypes, variable = "readcount") {
    phylosignal_summary <- subset(phylosignal_summary, var == variable)
    phylosignal_ihc <- merge(phylosignal_summary, ihc, by = c("patient_id"))
    phylosignal_ihc_melted <- melt(phylosignal_ihc, id.vars = c("patient_id", 
        "Cmean", "Cmean_p", "I", "I_p", "num"), measure.vars = tiltypes, variable.name = "tiltype", 
        value.name = "density")
    
    ## Filter for number of observations
    phylosignal_ihc_melted_filtered <- subset(phylosignal_ihc_melted, num >= 
        3)
    
    pvals <- setNames(ddply(phylosignal_ihc_melted_filtered, .(tiltype), function(x) {
        df <- as.data.frame(x)
        corres <- cor.test(df[, "density"], df[, phylosignal_measure], method = "spearman")
        
        pval <- corres$p.value
        eq <- substitute(italic(P) == p, list(p = format(pval, digits = 3)))
        return(as.character(as.expression(eq)))
    }), c("tiltype", "p.value"))
    
    p <- ggplot(phylosignal_ihc_melted_filtered, aes_string(x = "density", y = phylosignal_measure)) + 
        geom_point(aes(colour = patient_id)) + theme_bw() + theme_Publication() + 
        facet_wrap(~tiltype, scales = "free") + scale_color_manual(values = pal_patient) + 
        geom_text(data = pvals, aes(x = Inf, y = Inf, label = p.value), hjust = 1.1, 
            vjust = 1.5, size = 3, parse = TRUE)
    return(p)
}

mugration_phylosignal_stats <- get_phylosignal_table(mugration_statfiles)
nongeographic_phylosignal_stats <- get_phylosignal_table(nongeographic_statfiles)

mugration_phylosignal_summary <- summarize_phylosignal_stats(mugration_phylosignal_stats)
nongeographic_phylosignal_summary <- summarize_phylosignal_stats(nongeographic_phylosignal_stats)

Phylosignal-IHC analysis

p <- phylosignal_ihc_correlation(mugration_phylosignal_summary, "Cmean", ihc_means, 
    tiltypes, variable = "readcount")
plot(p)

p <- phylosignal_ihc_correlation(mugration_phylosignal_summary, "I", ihc_means, 
    tiltypes, variable = "readcount")
plot(p)

p <- phylosignal_ihc_correlation(mugration_phylosignal_summary, "Cmean", ihc_means, 
    tiltypes, variable = "mut_freqs")
plot(p)

p <- phylosignal_ihc_correlation(nongeographic_phylosignal_summary, "Cmean", 
    ihc_means, tiltypes, variable = "readcount")
plot(p)

p <- phylosignal_ihc_correlation(nongeographic_phylosignal_summary, "I", ihc_means, 
    tiltypes, variable = "readcount")
plot(p)

p <- phylosignal_ihc_correlation(nongeographic_phylosignal_summary, "Cmean", 
    ihc_means, tiltypes, variable = "mut_freqs")
plot(p)

So phylogenetic signal in readcounts correlates with CD20+ and CD4+ TIL densities. This is potentially evidence of selection occurring (clade-specific expansion) in samples with higher TIL densities, a finding that we might naively expect.

Note that I ran a few plot sets using the variable mut_freqs, which just refers to mutation frequency. We don’t expect phylogenetic signal for this parameter to correlate with any biological parameters, which is what we observe.

Phylosignal-ITH analysis

TODO.

Origin analysis

Description

Origin analysis = root frequency analysis. The goal of this analysis is to infer the geographical origin of each B-cell clonal family, equivalent to ancestral state analysis at the root node. These results can then be correlated with other known covariates.

Questions

Some really simple questions to ask include:

• Do clonal families from the same patient consistently arise from the same geographic origin?
• Is geographic origin restricted to certain anatomic sites?
• Are geographic origin/migration patterns correlated with tumour clone migration patterns?

In the previous version of this analysis, I observed a pattern where, in high-immune patients, geographic origin tended to be consistent across clonal families. The same trend was not present for low-immune clusters. A hypothetical biological explanation for this was that, in high-immune patients, B-cell clonal families respond to the same/similar antigenic threats, and therefore arose from the same site where those antigens were initially presented. Whereas, in low-immune patients, it was possible that B-cell clonal families were not tumour-reactive and were instead ‘sentinels’ patrolling each tumour site.

We check if this pattern holds up.

Results

read_root_frequencies <- function(f, clust, rep) {
    freqs <- fread(f)
    cbind(clust = clust, rep = rep, freqs)
}

get_rf_table <- function(files) {
    patient_clust_rep <- str_extract(files, "[0-9]+\\/clust[0-9]+/[0-9]+")
    strs <- strsplit(patient_clust_rep, "/")
    patient <- sapply(strs, function(x) x[1])
    clust <- sapply(strs, function(x) x[2])
    rep <- sapply(strs, function(x) x[3])
    
    info <- data.frame(patient_id = patient, clust = clust, rep = rep, tablefile = files, 
        stringsAsFactors = FALSE)
    
    freq_tables <- info %>% group_by(patient_id) %>% do(rf = rbind.fill(lapply(seq_along(.$tablefile), 
        function(i) read_root_frequencies(.$tablefile[i], .$clust[i], .$rep[i]))))
    
    freq_tables$rf <- lapply(freq_tables$rf, function(x) {
        roworder <- order(as.numeric(str_extract(x$clust, "[0-9]+")), x$rep)
        x <- x[roworder, ]
    })
    freq_tables <- freq_tables[order(as.numeric(freq_tables$patient_id)), ]
    return(freq_tables)
}

plot_rf_table <- function(rf, annotation_cols = c("clust", "rep"), title) {
    dat <- as.data.frame(rf)[, (!colnames(rf) %in% annotation_cols)]
    
    annotation_row <- subset(rf, select = c(annotation_cols))
    annotation_colours <- list(clust = get_colour_palette(annotation_row$clust), 
        rep = c(`0` = "#000000", `1` = "#FFFFFF"))
    
    pheatmap(dat, cluster_rows = FALSE, cluster_cols = FALSE, annotation_row = annotation_row, 
        annotation_colors = annotation_colours, show_rownames = FALSE, main = title)
}

rf_table <- get_rf_table(rootfreq_files)

ignore <- lapply(1:nrow(rf_table), function(i) {
    patient_id <- rf_table$patient_id[i]
    plot_title <- paste("Patient", patient_id, sep = " ")
    plot_rf_table(rf_table$rf[[i]], title = plot_title)
})

It’s pretty obvious from these plots that this pattern doesn’t hold up anymore. The only patient for which we see consistent anatomic origin is patient 2, where nearly all clonal families arise from the right ovary. The results for patient 2 corroborate our previous analysis, but no other patients do. Patient 1 and 22, which are two immune-high patients that previously followed that trend, no longer do. Likewise, patient 15 is a new (ITH-3) patient that is extremely immune-high and does not follow the trend. Therefore our hypothesis is not supported by the new data.