IHC-XCR statistics

library(ithi.utils)
load_base_libs()

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

Colour palettes

pal_patient <- select_palette("patient")

Parameters

db_path <- snakemake@params$db

ihc_table_path <- snakemake@input$ihc_table
tiltypes <- snakemake@params$ihc_xcr_tiltypes
tcr_diversity_file <- snakemake@input$tcr_diversity
bcr_diversity_file <- snakemake@input$bcr_diversity

Metadata

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

Results

tcr_diversity <- read_xcr_diversity_file(tcr_diversity_file, db_path)
bcr_diversity <- read_xcr_diversity_file(bcr_diversity_file, db_path)

ihc_table <- fread(ihc_table_path)
ihc_table_subset <- subset(ihc_table, select = c("condensed_id", "patient_id", 
    tiltypes))

Correlations with TIL densities

diversity_cols <- c("reads", "diversity", "observedDiversity_mean", "shannonWienerIndex_mean", 
    "inverseSimpsonIndex_mean")

dat <- list(tcr = tcr_diversity, bcr = bcr_diversity)
diversity <- lapply(names(dat), function(segment) {
    x <- dat[[segment]]
    x <- subset(x, select = c("condensed_id", "patient_id", diversity_cols))
    colnames(x) <- mapvalues(colnames(x), from = diversity_cols, to = paste(segment, 
        diversity_cols, sep = "_"))
    return(x)
})
names(diversity) <- names(dat)

df <- Reduce(f = function(x, y) merge(x, y, by = c("condensed_id", "patient_id")), 
    c(diversity, list(ihc_table_subset)))

df_melted <- melt(df, id.vars = colnames(df)[!colnames(df) %in% tiltypes], measure.vars = tiltypes, 
    variable.name = "tiltype", value.name = "density")

pval_plot <- function(df_melted, variable) {
    pvals <- setNames(ddply(df_melted, .(tiltype), function(x) {
        df <- as.data.frame(x)
        corres <- cor.test(df[, "density"], df[, variable], 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"))
    
    ggplot(df_melted, aes_string(x = "density", y = variable)) + geom_point(aes(colour = patient_id)) + 
        facet_wrap(~tiltype, scales = "free") + theme_bw() + theme_Publication() + 
        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) + 
        ggtitle(variable)
}

pval_plot(df_melted, "tcr_reads")

pval_plot(df_melted, "bcr_reads")

pval_plot(df_melted, "tcr_diversity")

pval_plot(df_melted, "bcr_diversity")

pval_plot(df_melted, "tcr_shannonWienerIndex_mean")

pval_plot(df_melted, "bcr_shannonWienerIndex_mean")