Molecular subtypes

library(ithi.utils)
load_base_libs()
library(nestedRanksTest)

library(ithi.meta)
library(ithi.xcr)
library(ithi.ihc)
library(ithi.expr)
library(ithi.clones)
library(nanotools)

Colour palettes

pal_patient <- select_palette("patient")

Parameters

db_path <- snakemake@params$db

molecular_subtype_file <- snakemake@input$molecular_subtypes

ihc_table_path <- snakemake@input$ihc_table
tiltypes <- unlist(snakemake@params$tiltypes)

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

xcr_table_path <- snakemake@input$xcr_table

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

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

Metadata

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

Molecular subtype classification

molsubtypes <- fread(molecular_subtype_file)
molsubtypes <- setNames(molsubtypes, c("condensed_id", "subtype"))

Counts

table(molsubtypes$subtype)

C1 C2 C4 C5 
24 36 15 41 

IHC

ihc_table <- fread(ihc_table_path)
ihc_table$patient_id <- factor(ihc_table$patient_id)

ihcm <- melt(ihc_table, id.vars = c("condensed_id", "patient_id"), measure.vars = tiltypes, 
    variable.name = "tiltype", value.name = "density")

ihcm_subtype <- merge(molsubtypes, ihcm, by = c("condensed_id"))

stats <- setNames(ddply(ihcm_subtype, .(tiltype), function(x) {
    res <- kruskal.test(density ~ factor(subtype), data = x)
    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"))

ggplot(ihcm_subtype, aes(x = subtype, y = density)) + geom_boxplot() + geom_jitter(position = position_jitter(width = 0.2, 
    height = 0), aes(colour = patient_id)) + theme_bw() + theme_Publication() + 
    scale_color_manual(values = pal_patient) + facet_wrap(~tiltype, scales = "free", 
    nrow = 2) + geom_text(data = stats, aes(x = Inf, y = Inf, label = p.value), 
    hjust = 1.1, vjust = 1.5, size = 3, parse = TRUE)

TCR/BCR diversity

diversity_files <- list(tcr = tcr_diversity_file, bcr = bcr_diversity_file)

diversity <- lapply(diversity_files, function(f) {
    dat <- read_xcr_diversity_file(f, db_path)
    merge(dat, molsubtypes, by = c("condensed_id"))
})

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

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

xcr_table$condensed_id <- as.character(xcr_table$condensed_id)
patients <- unique(xcr_table$patient_id)

tcr_segment_type <- "TRB"
bcr_segment_type <- "IGH"
id_type <- "condensed_id"

# not actually beta diversity but still fun
betas <- rbind.fill(lapply(patients, function(patient) {
    tcr_clonotypes <- subset(xcr_table, type == tcr_segment_type & patient_id == 
        patient)
    bcr_clonotypes <- subset(xcr_table, type == bcr_segment_type & patient_id == 
        patient)
    tcr_cross_table <- cross_tabulate(tcr_clonotypes, id_type = id_type)
    bcr_cross_table <- cross_tabulate(bcr_clonotypes, id_type = id_type)
    
    cross_tables <- list(tcr = tcr_cross_table, bcr = bcr_cross_table)
    
    betavals <- rbind.fill(lapply(seq_along(cross_tables), function(i) {
        cross_table <- cross_tables[[i]]
        df <- cbind(type = names(cross_tables)[i], compute_beta_diversity(cross_table, 
            mode = "union", measure = "count"))
        return(df)
    }))
    return(cbind(patient_id = patient, betavals))
}))
betas_list <- split(betas, f = betas$type)

compute_diversity_index <- function(xcr_table, type = "rareness", rare_threshold = 1e-04) {
    if (type == "rareness") {
        PREVALENCE_TYPE <- "freq"
        
        res <- xcr_table %>% group_by_(.dots = c("condensed_id", "patient_id", 
            "type")) %>% summarise_(rareness = lazyeval::interp(~sum(var[var < 
            rare_threshold]), var = as.name(PREVALENCE_TYPE)))
        res <- split(res, f = res$type)
    } else if (type == "berger_parker") {
        res <- xcr_table %>% group_by_(.dots = c("condensed_id", "patient_id", 
            "type")) %>% summarise(berger_parker = 1/max(freq))
        res <- split(res, f = res$type)
    }
    return(res)
}

cdf_clonotype_freqs <- function(xcr_table, segment = "TRB") {
    ggplot(xcr_table, aes(x = freq)) + stat_ecdf(geom = "step") + scale_x_log10() + 
        facet_wrap(~condensed_id, ncol = 5)
}

plot_diversity_subtype <- function(diversity, col) {
    ignore <- lapply(diversity, function(df) {
        p <- ggplot(df, aes_string(x = "subtype", y = col)) + geom_boxplot() + 
            geom_jitter(position = position_jitter(width = 0.2, height = 0), 
                aes(colour = patient_id)) + theme_bw() + theme_Publication() + 
            scale_color_manual(values = pal_patient) + xlab("Subtype") + ylab("Diversity")
        print(p)
    })
}

Richness

plot_diversity_subtype(diversity, "observedDiversity_mean")

Shannon entropy

plot_diversity_subtype(diversity, "shannonWienerIndex_mean")

Divergence

plot_diversity_subtype(diversity, "inverseSimpsonIndex_mean")

D50 index

plot_diversity_subtype(diversity, "d50Index_mean")

Rareness

rareness <- compute_diversity_index(xcr_table, type = "rareness", rare_threshold = 1e-04)
rareness <- lapply(rareness, function(x) {
    x$patient_id <- as.character(x$patient_id)
    merge(x, molsubtypes, by = c("condensed_id"))
})
plot_diversity_subtype(rareness, "rareness")

Berger-Parker

dominance <- compute_diversity_index(xcr_table, type = "berger_parker")
dominance <- lapply(dominance, function(x) {
    x$patient_id <- as.character(x$patient_id)
    merge(x, molsubtypes, by = c("condensed_id"))
})
plot_diversity_subtype(dominance, "berger_parker")

Sample-specific clonotype proportion

betas_list <- lapply(betas_list, function(x) {
    x$patient_id <- as.character(x$patient_id)
    merge(x, molsubtypes, by = c("condensed_id"))
})
plot_diversity_subtype(betas_list, "prop_unique")

CDFs

cdf_clonotype_freqs(xcr_table, "TRB")

cdf_clonotype_freqs(xcr_table, "IGH")

Nanostring

exprs <- fread(nanostring_data_path)
labels <- fread(nanostring_annotations_path)

celltype_matrix <- rownames_to_column(as.data.frame(t(create_celltype_matrix(exprs, 
    labels, db_path))), var = "condensed_id")
celltype_vars <- colnames(celltype_matrix)[!colnames(celltype_matrix) == "condensed_id"]

celltype_matrix_melted <- melt(celltype_matrix, id.vars = c("condensed_id"), 
    measure.vars = celltype_vars, variable.name = "type", value.name = "expr")

celltype_subtype <- merge(celltype_matrix_melted, molsubtypes, by = c("condensed_id"))
celltype_subtype$patient_id <- factor(map_id(celltype_subtype$condensed_id, 
    from = "condensed_id", to = "patient_id", db_path))

stats <- setNames(ddply(celltype_subtype, .(type), function(x) {
    res <- kruskal.test(expr ~ factor(subtype), data = x)
    pvalue <- res$p.value
    eq <- substitute(italic(P) == p, list(p = format(pvalue, digits = 3)))
    return(as.character(as.expression(eq)))
}), c("type", "p.value"))

ggplot(celltype_subtype, aes(x = subtype, y = expr)) + geom_boxplot() + geom_jitter(position = position_jitter(width = 0.2, 
    height = 0), aes(colour = patient_id)) + theme_bw() + theme_Publication() + 
    scale_color_manual(values = pal_patient) + facet_wrap(~type, scales = "free", 
    nrow = 4) + geom_text(data = stats, aes(x = Inf, y = Inf, label = p.value), 
    hjust = 1.1, vjust = 1.5, size = 3, parse = TRUE)

pathway_matrix <- rownames_to_column(as.data.frame(t(create_pathway_matrix(exprs, 
    labels, db_path))), var = "condensed_id")
pathway_vars <- colnames(pathway_matrix)[!colnames(pathway_matrix) == "condensed_id"]

pathway_matrix_melted <- melt(pathway_matrix, id.vars = c("condensed_id"), measure.vars = pathway_vars, 
    variable.name = "type", value.name = "expr")

pathway_subtype <- merge(pathway_matrix_melted, molsubtypes, by = c("condensed_id"))
pathway_subtype$patient_id <- factor(map_id(pathway_subtype$condensed_id, from = "condensed_id", 
    to = "patient_id", db_path))

stats <- setNames(ddply(pathway_subtype, .(type), function(x) {
    res <- kruskal.test(expr ~ factor(subtype), data = x)
    pvalue <- res$p.value
    eq <- substitute(italic(P) == p, list(p = format(pvalue, digits = 3)))
    return(as.character(as.expression(eq)))
}), c("type", "p.value"))

ggplot(pathway_subtype, aes(x = subtype, y = expr)) + geom_boxplot() + geom_jitter(position = position_jitter(width = 0.2, 
    height = 0), aes(colour = patient_id)) + theme_bw() + theme_Publication() + 
    scale_color_manual(values = pal_patient) + facet_wrap(~type, scales = "free", 
    nrow = 4) + geom_text(data = stats, aes(x = Inf, y = Inf, label = p.value), 
    hjust = 1.1, vjust = 1.5, size = 3, parse = TRUE)

Relation to clonal composition

Here we’ll relate intrapatient differences in molecular subtype to differences in clonal composition. Questions: * Are tumours with different molecular subtypes more clonally different?

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 = "condensed_id")
ccfs_labeled <- merge(ccfs, branch_lengths, by = c("label", "node", "patient_id"))

## Need non-normalized distances
clone_dists <- clone_distances(ccfs_labeled, normalize = TRUE, id_type = "condensed_id")

clone_dists$patient_id <- factor(clone_dists$patient_id)

clone_dists_melted <- rbind.fill(lapply(1:nrow(clone_dists), function(i) {
    patient <- clone_dists$patient_id[i]
    distmat <- clone_dists$dist_clones_weighted[[i]]
    
    distdf <- setNames(melt(as.matrix(distmat)), c("id1", "id2", "clone_dist"))
    distdf$id1 <- factor(distdf$id1)
    distdf$id2 <- factor(distdf$id2)
    distdf <- subset(distdf, as.numeric(id1) > as.numeric(id2))
    return(distdf)
}))

df <- merge(clone_dists_melted, setNames(molsubtypes, c("id1", "subtype1")))
df <- merge(df, setNames(molsubtypes, c("id2", "subtype2")))
df$equal_subtypes <- df$subtype1 == df$subtype2
df$patient_id <- map_id(df$id1, from = "condensed_id", to = "patient_id", db_path)

ggplot(df, aes(x = equal_subtypes, y = clone_dist)) + geom_boxplot(outlier.size = -1) + 
    geom_jitter(aes(colour = patient_id), position = position_jitter(width = 0.2, 
        height = 0)) + theme_bw() + theme_Publication() + scale_colour_manual(values = pal_patient)

Making this comparison across patients is unfair, as clone distance values (even when normalized) may be very different in scale from patient to patient. It’s obvious there’s no pattern ACROSS patients, but how about WITHIN patients?

If we assume a Gaussian response variable (which is definitely not true):

summary(lmer(clone_dist ~ equal_subtypes + (1 | patient_id), df))

Linear mixed model fit by REML t-tests use Satterthwaite approximations
  to degrees of freedom [lmerMod]
Formula: clone_dist ~ equal_subtypes + (1 | patient_id)
   Data: df

REML criterion at convergence: -244.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.4180 -0.4943 -0.0424  0.5238  2.3332 

Random effects:
 Groups     Name        Variance Std.Dev.
 patient_id (Intercept) 0.02447  0.1564  
 Residual               0.01821  0.1349  
Number of obs: 252, groups:  patient_id, 14

Fixed effects:
                     Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)          0.250493   0.044280  14.090000   5.657 5.77e-05 ***
equal_subtypesTRUE  -0.005325   0.020301 242.570000  -0.262    0.793    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr)
eql_sbtTRUE -0.221

There’s a trend toward lower clone distance between pairs of tumours with equal subtype, but this is not significant.

What if we use a nonparametric test?

tmp <- df %>% group_by(patient_id) %>% summarise(both_levels = !Reduce("&", 
    equal_subtypes) & Reduce("|", equal_subtypes))

df_filtered <- subset(df, patient_id %in% subset(tmp, both_levels)$patient_id)

res <- nestedRanksTest(clone_dist ~ !equal_subtypes | patient_id, data = df_filtered, 
    n.iter = 5000)
res

    Nested Ranks Test

data:  clone_dist by !equal_subtypes grouped by patient_id
Z = 0.00081235, p-value = 0.4866
alternative hypothesis: Z lies above bootstrapped null values
null values:
      0%       1%       5%      10%      25%      50%      75%      90% 
-0.36149 -0.25509 -0.17872 -0.14054 -0.07392 -0.00244  0.07311  0.14054 
     95%      99%     100% 
 0.17551  0.24696  0.34768 

bootstrap iterations: 5000 
group weights:
         10          11          12          13          14          15 
0.158407799 0.017059301 0.017059301 0.019496344 0.012997563 0.006498781 
         17           2           3           4           7           9 
0.019496344 0.004061738 0.604386677 0.089358245 0.043866775 0.007311129 

plot(res)

We do see a trend! We can only do this for patients with samples in both clusters. So, samples with different molecular subtypes are also more likely to have different clonal compositions. But this does not appear to be significant when I run it with the new subtyping (combined using ICGC for subtyping …).