XCR replicate QC

library(ithi.utils)
load_base_libs()

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

Colour palettes

pal_patient <- select_palette("patient")

Parameters

db_path <- snakemake@params$db

xcr_table_path <- snakemake@input$xcr_table
prevalence_option <- snakemake@params$prevalence_option
distance_method <- snakemake@params$xcr_distance_method
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"))
duplicates <- collect(tbl(db, "duplicates"))
xcrseq <- collect(tbl(db, "xcrseq"))

TCR/BCR

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

Read 19.7% of 304822 rows
Read 52.5% 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 <- subset(xcr_table, count >= 200)

technical_duplicates <- unique(subset(xcrseq, voa_unique != voa)$voa)

duplicate_sites <- subset(duplicates, duplicate == 1)
duplicate_sites <- merge(samples, duplicate_sites, by = "voa")
duplicate_categories <- subset(duplicate_sites, select = c(patient_id, site_id))

duplicates_mapped <- merge(duplicate_categories, samples, by = c("patient_id", 
    "site_id"))
duplicates_mapped <- merge(duplicates_mapped, xcrseq, by = c("voa"))
mapped_duplicate_categories <- duplicates_mapped %>% group_by(patient_id, site_id) %>% 
    summarise(num_voa = length(unique(voa)))
mapped_categories <- subset(mapped_duplicate_categories, num_voa > 1, select = c("patient_id", 
    "site_id"))

biological_duplicates <- unique(merge(duplicates_mapped, mapped_categories, 
    by = c("patient_id", "site_id"))$voa)

tcr_segment_type <- "TRB"
bcr_segment_type <- "IGH"

id_type <- "voa_unique"

Technical duplicates

xcr_table_tdups <- subset(xcr_table, voa %in% technical_duplicates)

tcr_tdups <- subset(xcr_table_tdups, type == tcr_segment_type)
tcr_cross_table <- cross_tabulate(tcr_tdups, id_type = id_type)

bcr_tdups <- subset(xcr_table_tdups, type == bcr_segment_type)
bcr_cross_table <- cross_tabulate(bcr_tdups, id_type = id_type)

cross_tdups <- list(tcr = tcr_cross_table, bcr = bcr_cross_table)

distance_matrices <- lapply(cross_tdups, function(cross_table) {
    distmat <- compute_immune_distance_matrix(cross_table, method = distance_method)
    mat <- as.matrix(distmat)
    return(mat)
})

sims_raw <- lapply(distance_matrices, function(mat) {
    sims <- average_intrapatient_similarity(mat, filter_tissue = FALSE, id_type = id_type, 
        db_path = db_path, raw = TRUE, group = "voa", table = "xcrseq")
    sims <- setNames(subset(sims, select = -c(tissue1, tissue2, patient2)), 
        c("sample1", "sample2", "dist", "group"))
    sims$patient_id <- map_id(sims$group, from = "voa", to = "patient_id", db_path = db_path, 
        table = "samples")
    sims$condensed_id <- map_id(sims$group, from = "voa", to = "condensed_id", 
        db_path = db_path, table = "samples")
    project1 <- map_id(sims$sample1, from = "voa_unique", to = "project_id", 
        db_path = db_path, table = "xcrseq")
    project2 <- map_id(sims$sample2, from = "voa_unique", to = "project_id", 
        db_path = db_path, table = "xcrseq")
    sims$comparison_type <- c("different", "identical")[as.numeric(project1 == 
        project2) + 1]
    sims_summary <- sims %>% group_by(condensed_id, patient_id, comparison_type) %>% 
        summarise(mean_sim = 1 - mean(dist))
    return(sims_summary)
})

ignore <- lapply(list("tcr", "bcr"), function(segment) {
    dat <- sims_raw[[segment]]
    dat$patient_id <- factor(dat$patient_id)
    dat$comparison_type <- factor(dat$comparison_type)
    
    p <- ggplot(dat, aes(x = patient_id, y = mean_sim)) + geom_boxplot(aes(fill = comparison_type)) + 
        geom_jitter(position = position_jitter(width = 0.2, height = 0)) + scale_fill_manual(values = get_colour_palette(dat$comparison_type)) + 
        theme_bw() + theme_Publication() + xlab("Patient") + ylab("Pairwise repertoire similarity")
    print(p)
})

We might have good reason to suspect that duplicates conducted at different times (i.e. the original TCR/BCR-seq experiment years ago vs. the one now) might not give the same results. The RIN values for the old samples that we re-ran were quite low (3-5, compared to the usual >7). Indeed, the BCR data shows that duplicates from the same run are more similar to each other than those from different runs.

However, one case is alarming – patient 21 in the TCR data. I made an ID mapping error and missed most samples from patient 21 in the first round of processing. However, the refactoring and database construction stimulated by the newest run allowed me to recover these samples. Something fishy definitely seems to be going on here. Subsetting to clonotypes with >= 50, 200 reads did not help, indicating that the issue isn’t likely to be caused by low-level contamination. It may still be caused by high-level contamination – this is worth investigating as patient 21 has a very low total clonotype count.

Next steps: Check RIN values from the original run for patient 21’s samples. Answer: They are on the low side for IX3248 (median = 5.2), but there are several patients around the 5.4-5.7 range. So I’m tempted to say that isn’t the cause of the issue.

TODO: Also evaluate variation in diversity values.

Biological duplicates

Another class of duplicates we tested were parts of samples coming from the same tumour sample. That’s a bit confusing – what it means is that each tumour sample was split into 2 fresh frozen and 1 FFPE aliquot, and for some cases we subjected both fresh frozen aliquots to TCR/BCR-seq. I loosely refer to these cases as ‘biological duplicates’.

id_type <- "voa"

xcr_table_bdups <- subset(xcr_table, voa_unique %in% biological_duplicates)

tcr_bdups <- subset(xcr_table_bdups, type == tcr_segment_type)
tcr_cross_table <- cross_tabulate(tcr_bdups, id_type = id_type)

bcr_bdups <- subset(xcr_table_bdups, type == bcr_segment_type)
bcr_cross_table <- cross_tabulate(bcr_bdups, id_type = id_type)

cross_bdups <- list(tcr = tcr_cross_table, bcr = bcr_cross_table)

distance_matrices <- lapply(cross_bdups, function(cross_table) {
    distmat <- compute_immune_distance_matrix(cross_table, method = distance_method)
    mat <- as.matrix(distmat)
    return(mat)
})

sims_raw <- lapply(distance_matrices, function(mat) {
    sims <- average_intrapatient_similarity(mat, filter_tissue = FALSE, id_type = id_type, 
        db_path = db_path, raw = TRUE, group = "condensed_id", table = "samples")
    sims <- setNames(subset(sims, select = -c(tissue1, tissue2, patient2)), 
        c("sample1", "sample2", "dist", "condensed_id"))
    sims$patient_id <- map_id(sims$condensed_id, from = "condensed_id", to = "patient_id", 
        db_path = db_path, table = "samples")
    project1 <- map_id(sims$sample1, from = "voa_unique", to = "project_id", 
        db_path = db_path, table = "xcrseq")
    project2 <- map_id(sims$sample2, from = "voa_unique", to = "project_id", 
        db_path = db_path, table = "xcrseq")
    sims$comparison_type <- c("different", "identical")[as.numeric(project1 == 
        project2) + 1]
    sims_summary <- sims %>% group_by(condensed_id, patient_id, comparison_type) %>% 
        summarise(mean_sim = 1 - mean(dist))
    return(sims_summary)
})

ignore <- lapply(list("tcr", "bcr"), function(segment) {
    dat <- sims_raw[[segment]]
    dat$patient_id <- factor(dat$patient_id)
    dat$comparison_type <- factor(dat$comparison_type)
    
    p <- ggplot(dat, aes(x = patient_id, y = mean_sim)) + geom_boxplot(aes(fill = comparison_type)) + 
        geom_jitter(position = position_jitter(width = 0.2, height = 0)) + scale_fill_manual(values = get_colour_palette(dat$comparison_type)) + 
        theme_bw() + theme_Publication() + xlab("Patient") + ylab("Pairwise repertoire similarity")
    print(p)
})

A little surprising is how low the repertoire similarity is between biological duplicates – from aliquots that are taken from the same tumour sample! Consistent with our findings that high T-cell abundance = high levels of intrapatient similarity, patient 15 stands out as a beacon of consistency.