Immunoediting

library(ithi.utils)
load_base_libs()

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

Colour palettes

pal_patient <- select_palette("patient")

Parameters

db_path <- snakemake@params$db

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

neoediting_outdir <- snakemake@input$neoediting_outdir
epitope_dirs <- list.files(file.path(neoediting_outdir, "epitope_files"), full.names = TRUE)
spectrum_file <- file.path(neoediting_outdir, "background_spectrum", "signatures.txt")

Metadata

id_type <- "condensed_id"

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

TIL clusters

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"), to = c("E CD8+", 
    "E CD4+", "E CD20+", "E 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))

Immunoediting

read_spectrum <- function(f) {
    dat <- fread(f)
    colnames(dat) <- mapvalues(colnames(dat), from = c("mutationseq_trinucleotide_ref", 
        "mutationseq_trinucleotide_alt"), to = c("triref", "trialt"))
    return(dat)
}

read_epitope_files <- function(epitope_files, spectrum) {
    epitopes <- rbind.fill(lapply(epitope_files, function(f) {
        sample_key <- basename(dirname(f))
        vars <- fread(f)
        cbind(sample_key = sample_key, vars)
    }))
    epitopes <- merge(epitopes, spectrum, by = c("trinuc_sig"))
    epitopes$patient_id <- map_id(epitopes$sample_key, from = "condensed_id", 
        to = "patient_id", db_path)
    return(epitopes)
}

spectrum <- read_spectrum(spectrum_file)

silent_epitope_files <- Sys.glob(file.path(epitope_dirs, "silent.tsv"))
nonsynonymous_epitope_files <- Sys.glob(file.path(epitope_dirs, "nonsynonymous.tsv"))

silent_epitopes <- read_epitope_files(silent_epitope_files, spectrum)
silent_dups <- duplicated(subset(silent_epitopes, select = c(sample_key, Chromosome, 
    Start, Reference, Variant)))
silent_epitopes <- silent_epitopes[!silent_dups, ]

nonsynonymous_epitopes <- read_epitope_files(nonsynonymous_epitope_files, spectrum)

Overall, there aren’t any major differences between the fraction of neoepitopes determined from nonsynonymous mutations vs. what would be expected from synonymous mutations. In other words, there is, on average, no immunoediting occurring on neoantigens in HGSC. Do keep in mind this is complicated by the lack of expression data, so we cannot truly confirm which are neoantigens.

sum(silent_epitopes$frac_binder)/nrow(silent_epitopes)

[1] 0.7216169

sum(silent_epitopes$frac_strong_binder)/nrow(silent_epitopes)

[1] 0.3326825

sum(nonsynonymous_epitopes$frac_binder)/nrow(nonsynonymous_epitopes)

[1] 0.71962

sum(nonsynonymous_epitopes$frac_strong_binder)/nrow(nonsynonymous_epitopes)

[1] 0.3303586

Allelic fractions

hist(with(subset(nonsynonymous_epitopes, `NetMHCpan MT Percentile` > 2), alt_counts/(ref_counts + 
    alt_counts)), breaks = 30)

hist(with(subset(nonsynonymous_epitopes, `NetMHCpan MT Percentile` <= 2), alt_counts/(ref_counts + 
    alt_counts)), breaks = 30)

tab1 <- table(subset(nonsynonymous_epitopes, `NetMHCpan MT Percentile` > 2)$is_ancestral)
tab2 <- table(subset(nonsynonymous_epitopes, `NetMHCpan MT Percentile` <= 2)$is_ancestral)
hist(with(silent_epitopes, alt_counts/(ref_counts + alt_counts)), breaks = 30)

fisher.test(rbind(tab1, tab2))

    Fisher's Exact Test for Count Data

data:  rbind(tab1, tab2)
p-value = 0.0101
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.6776465 0.9518279
sample estimates:
odds ratio 
 0.8033361 

What does the peak near 0.5 mean? We might normally expect to see a peak at 0.5 as a result of heterozygosity, for tumours with high purity. Why does it only occur for the silent epitopes? Could it be a result of editing occuring for the nonsynonymous epitopes?

Nevermind, peak at 0.5 is gone. Was due to some really weird variant that had somehow passed my uniqueness filter.

Overall rates

silent_rates <- silent_epitopes %>% group_by(sample_key) %>% summarise(binders = sum(frac_binder), 
    strong_binder = sum(frac_strong_binder), total = n(), spbinder = binders/total, 
    spsbinder = strong_binder/total)
nonsynonymous_rates <- nonsynonymous_epitopes %>% group_by(sample_key) %>% summarise(ntotal = n(), 
    nepitopes = length(which(`NetMHCpan MT Percentile` <= 2)), nstrong = length(which(`NetMHCpan MT Percentile` <= 
        0.5)))


# nonsynonymous_rates <- nonsynonymous_epitopes %>% group_by(sample_key) %>%
# summarise(binders=sum(frac_binder), strong_binder=sum(frac_strong_binder),
# total=n(), npbinder=binders/total, npsbinder=strong_binder/total)

rates <- merge(silent_rates, nonsynonymous_rates, by = c("sample_key"))

rates$pred_epitopes <- with(rates, spbinder * ntotal)
rates$pred_strong <- with(rates, spsbinder * ntotal)

vs. TIL clusters

rates_clusts <- merge(rates, clusts, by.x = c("sample_key"), by.y = c("condensed_id"))

rates_clusts$patient_id <- map_id(rates_clusts$sample_key, from = "condensed_id", 
    to = "patient_id", db_path)

ggplot(rates_clusts, aes(x = cluster, y = nepitopes/pred_epitopes)) + geom_boxplot() + 
    theme_bw() + theme_Publication() + ylab("Observed/expected ratio") + geom_jitter(aes(colour = factor(patient_id)), 
    position = position_jitter(width = 0.2, height = 0)) + scale_color_manual(values = pal_patient)

ggplot(rates_clusts, aes(x = cluster, y = nstrong/pred_strong)) + geom_boxplot() + 
    theme_bw() + theme_Publication() + ylab("Observed/expected ratio")

This doesn’t show a significant correlation. There is a chance we’ll see better results by looking directly at TIL densities though.

vs. TIL densities

rates_ihc <- merge(rates, ihc_table, by.x = "sample_key", by.y = "condensed_id")

Binders

ggplot(rates_ihc, aes(x = E_CD8_density, y = nepitopes/pred_epitopes)) + theme_bw() + 
    theme_Publication() + ylab("Observed/expected ratio") + geom_point(aes(colour = factor(patient_id))) + 
    scale_color_manual(values = pal_patient) + facet_wrap(~patient_id, scales = "free")

From this plot we can see a negative correlation between obs/expected ratio and epithelial CD8+ density for some patients. It’s also apparent that the patients for which we don’t see this correlation are patients with low levels of TILs (so that the values on the x-axis are very tightly spaced). We formulate the following linear mixed model:

mod <- lmer(nepitopes/pred_epitopes ~ E_CD8_density + (1 | patient_id), subset(rates_ihc, 
    !is.na(E_CD8_density)))
summary(mod)

Linear mixed model fit by REML t-tests use Satterthwaite approximations
  to degrees of freedom [lmerMod]
Formula: nepitopes/pred_epitopes ~ E_CD8_density + (1 | patient_id)
   Data: subset(rates_ihc, !is.na(E_CD8_density))

REML criterion at convergence: -63.4

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.08176 -0.40628  0.04027  0.56379  1.79161 

Random effects:
 Groups     Name        Variance Std.Dev.
 patient_id (Intercept) 0.043792 0.20926 
 Residual               0.006302 0.07939 
Number of obs: 55, groups:  patient_id, 13

Fixed effects:
                Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)    1.0059464  0.0608862 12.7700000  16.522 5.39e-10 ***
E_CD8_density -0.0012898  0.0004055 47.5800000  -3.181  0.00258 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr)
E_CD8_dnsty -0.233

anova(mod)

So we’ve demonstrated that observed/expected neoepitopes is negatively correlated with epithelial CD8+ density – consistent with the theory that immunoediting of neoepitopes occurs in samples with higher epithelial CD8+ TIL.

Strong binders

ggplot(rates_ihc, aes(x = E_CD8_density, y = nstrong/pred_strong)) + theme_bw() + 
    theme_Publication() + ylab("Observed/expected ratio") + geom_point(aes(colour = factor(patient_id))) + 
    scale_color_manual(values = pal_patient) + facet_wrap(~patient_id, scales = "free")

mod <- lmer(nstrong/pred_strong ~ E_CD8_density + (1 | patient_id), rates_ihc)
summary(mod)

Linear mixed model fit by REML t-tests use Satterthwaite approximations
  to degrees of freedom [lmerMod]
Formula: nstrong/pred_strong ~ E_CD8_density + (1 | patient_id)
   Data: rates_ihc

REML criterion at convergence: 23.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.1649 -0.6097 -0.1265  0.5142  2.4078 

Random effects:
 Groups     Name        Variance Std.Dev.
 patient_id (Intercept) 0.10230  0.3198  
 Residual               0.04031  0.2008  
Number of obs: 55, groups:  patient_id, 13

Fixed effects:
                Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)    1.0460945  0.0993728 12.5900000  10.527 1.33e-07 ***
E_CD8_density -0.0017944  0.0009616 52.8100000  -1.866   0.0676 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr)
E_CD8_dnsty -0.342

anova(mod)

This isn’t significant. One could attribute this to the lower numbers of strong epitopes used for computing the background spectrum – leading to larger uncertainty. To offset this we could get more samples to use for the background. It would be worth asking Rooney et al. to provide their spectrum which is presumably computed on a much larger set of samples.

Also it’s worth noting that we can do this analysis on TCGA too, since we don’t need noncoding mutations. So I should probably download the TCGA bams too and process them.