XCR-Mapscape plots

## Modified version of mapscape. Cannot be put into the same package as some
## of the components conflict
library(xcrmapscape)
library(htmlwidgets)

library(ithi.utils)
load_base_libs()
library(ithi.meta)
library(ithi.ihc)
library(ithi.xcr)
library(ithi.figures)
library(ithi.clones)

Parameters

db_path <- snakemake@params$db

mapscape_input_directories = snakemake@input$mapscape_inputs
xcr_table_path <- snakemake@input$xcr_table
show_full_track <- (snakemake@params$track_option == "full")
full_page_mode <- (snakemake@params$page_mode == "full")
segment_type <- snakemake@params$segment_type

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

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

distance_method <- "horn"

tree_branch_data <- read_clone_tree_data(clone_tree_file, clone_branch_length_file, 
    clone_prevalence_file, db_path)
xcr_table <- read_clonotypes(xcr_table_path, duplicates = FALSE, db_path)

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

annotation_colours <- ithi.figures::get_annotation_colours()
tcr_segment_type <- "TRB"
bcr_segment_type <- "IGH"

patients <- unique(xcr_table$patient_id)

id_type <- "condensed_id"

dists <- 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 <- ithi.xcr::cross_tabulate(tcr_clonotypes, id_type = id_type)
    bcr_cross_table <- ithi.xcr::cross_tabulate(bcr_clonotypes, id_type = id_type)
    
    cross_tables <- list(tcr = tcr_cross_table, bcr = bcr_cross_table)
    
    distance_matrices <- lapply(cross_tables, function(cross_table) {
        distmat <- ithi.xcr::compute_immune_distance_matrix(cross_table, method = distance_method)
        mat <- as.matrix(distmat)
        return(mat)
    })
    return(distance_matrices)
})
tcr_dists <- lapply(dists, function(x) x$tcr)
bcr_dists <- lapply(dists, function(x) x$bcr)

xcr_dists <- tibble::tibble(patient_id = patients, tcr = tcr_dists, bcr = bcr_dists)
xcr_dists$patient_id <- factor(xcr_dists$patient_id)
xcr_dists <- xcr_dists[order(xcr_dists$patient_id), ]


trees <- lapply(tree_branch_data, function(x) x$tree)
prevalences <- plyr::rbind.fill(lapply(tree_branch_data, function(x) x$prevalence_dat))
branch_lengths <- plyr::rbind.fill(lapply(tree_branch_data, function(x) x$branch_dat))

ccfs <- ithi.clones::compute_ccf(prevalences, trees, id_type = id_type)
ccfs_labeled <- merge(ccfs, branch_lengths, by = c("label", "node", "patient_id"))

## Need non-normalized distances
clone_dists <- ithi.clones::clone_distances(ccfs_labeled, normalize = FALSE, 
    id_type = "condensed_id")

clone_dists$patient_id <- factor(clone_dists$patient_id)

total_dists <- inner_join(clone_dists, xcr_dists)


pairwise_similarities <- plyr::rbind.fill(lapply(1:nrow(total_dists), function(i) {
    clone_distmat <- total_dists[i, ]$dist_clones_weighted[[1]]
    tcr_distmat <- total_dists[i, ]$tcr[[1]]
    bcr_distmat <- total_dists[i, ]$bcr[[1]]
    
    inputs <- list(clone_distmat, tcr_distmat, bcr_distmat)
    names(inputs) <- c("clones", "tcr", "bcr")
    
    patient_id <- total_dists[i, ]$patient_id
    
    pairwise_sims <- cbind(patient_id = patient_id, ithi.supp:::compute_overlap_similarities(inputs))
    return(pairwise_sims)
}))
pairwise_similarities$patient_id <- factor(pairwise_similarities$patient_id)

dissimilarity_matrices <- dplyr::bind_rows(lapply(1:nrow(total_dists), function(i) {
    clone_distmat <- total_dists[i, ]$dist_clones_weighted[[1]]
    tcr_distmat <- total_dists[i, ]$tcr[[1]]
    bcr_distmat <- total_dists[i, ]$bcr[[1]]
    patient_id <- total_dists[i, ]$patient_id
    
    inputs <- list(clone_distmat, tcr_distmat, bcr_distmat)
    names(inputs) <- c("clones", "tcr", "bcr")
    
    matrices_subset <- ithi.supp:::compute_overlap_matrices(inputs)
    names(matrices_subset) <- names(inputs)
    
    tibble(patient_id = patient_id, clones = list(matrices_subset$clones), tcr = list(matrices_subset$tcr), 
        bcr = list(matrices_subset$bcr))
}))

segment_name <- plyr::mapvalues(segment_type, from = c("TRB", "IGH"), to = c("tcr", 
    "bcr"))

stats <- plyr::rbind.fill(lapply(1:nrow(dissimilarity_matrices), function(i) {
    row <- dissimilarity_matrices[i, ]
    xcrmat <- dissimilarity_matrices[, segment_name][[1]][[i]]
    clonemat <- dissimilarity_matrices$clones[[i]]
    
    manttest <- vegan::mantel(xcrmat, clonemat, method = "pearson")
    data.frame(patient_id = row$patient_id, p.value = manttest$signif, stat = manttest$statistic)
}))

plot_xcrmapscape <- function(clonal_prev, sample_locations, tree_edges, chord_data, 
    track_data, resize_factor, plot_title, full_page_mode, pvalue) {
    img_ref = "/shahlab/alzhang/projects/ITH_Immune/paper/miscellaneous/FF4D00-0.8.png"
    # img_ref =
    # '~/shahlab/projects/ITH_Immune/paper/miscellaneous/FF4D00-0.8.png'
    
    if (!full_page_mode) {
        xcrmapscape::xcrmapscape(clonal_prev = clonal_prev, tree_edges = tree_edges, 
            sample_locations = sample_locations, img_ref = img_ref, show_warnings = FALSE, 
            chord_data = chord_data, track_data = track_data, resize_factor = resize_factor, 
            plot_title = plot_title, phylo_stroke_width = "5px", phylo_stroke_color = "#000000", 
            legendTitleHeight = 26, centerSize = 0.6, treeRelativeSize = 0, 
            titleScaleFactor = 28, pvalue = pvalue)
    } else {
        xcrmapscape::xcrmapscape(clonal_prev = clonal_prev, tree_edges = tree_edges, 
            sample_locations = sample_locations, img_ref = img_ref, show_warnings = FALSE, 
            chord_data = chord_data, track_data = track_data, resize_factor = resize_factor, 
            plot_title = plot_title, phylo_stroke_width = "1px", phylo_stroke_color = "#9E9A9A", 
            legendTitleHeight = 16, centerSize = 0.5, treeRelativeSize = 0.33, 
            titleScaleFactor = 40, pvalue = pvalue)
    }
}

intersect_table <- compute_overlap_table(xcr_table, segment_type = segment_type, 
    prevalence_option = "clones", id_type = "condensed_id")
intersect_table$patient1 <- ithi.meta::map_id(intersect_table$sample1, from = "condensed_id", 
    to = "patient_id", db_path)
intersect_table$patient2 <- ithi.meta::map_id(intersect_table$sample2, from = "condensed_id", 
    to = "patient_id", db_path)

intersect_table <- subset(intersect_table, patient1 == patient2)

clonotype_table <- subset(xcr_table, type == segment_type)
clonotype_counts <- clonotype_table %>% group_by(condensed_id, patient_id) %>% 
    summarise(n = n()) %>% plyr::rename(c(condensed_id = "label", n = "len"))

reformat_andrew_table <- function(dat, db_path, patient_id) {
    dat$patient_id <- patient_id
    dat <- subset(dat, !sample_id %in% c("Pv2", "Rec2", "Rec3"))
    dat <- dat %>% plyr::rename(c(sample_id = "plot_id", long_sample_id = "sample_id"))
    dat_remapped <- ithi.meta::remap_sample_ids(dat, db_path)
    dat_remapped$abbrev_id <- stringr::str_extract(dat_remapped$condensed_id, 
        "[A-Za-z0-9]+$")
    dat_remapped <- dat_remapped %>% plyr::rename(c(sample_id = "old_sample_id", 
        sample_key = "long_sample_id", plot_id = "old_plot_id", abbrev_id = "sample_id"))
    if ("location_id" %in% colnames(dat_remapped)) {
        dat_remapped$location_id <- dat_remapped$sample_id
    }
    return(dat_remapped)
}

export_cydney_inputs <- function(xcr_table, patient_id, mapscape_base_dir = "/Users/alzhang/shahlab/projects/ITH_Immune/paper/results/intermediates/run2/mapscape_postprocessed_HMC", 
    outdir) {
    patient_tmp <- patient_id
    clonotypes_sub <- xcr_table[(xcr_table$type == segment_type) & (xcr_table$patient_id == 
        as.numeric(patient_tmp)), ]
    clonotypes_sub$condensed_id <- factor(clonotypes_sub$condensed_id)
    cross_table <- ithi.xcr::cross_tabulate(clonotypes_sub, id_type = "condensed_id")
    colnames(cross_table) <- str_replace(colnames(cross_table), "^[0-9]+_", 
        "")
    
    bin_cross_table <- cross_table
    bin_cross_table[bin_cross_table > 0] <- 1
    coln <- colnames(cross_table)
    
    cross_summary <- bin_cross_table %>% group_by_(.dots = coln) %>% summarise(n_shared_clonotypes = n())
    
    dir <- file.path(mapscape_base_dir, patient_id)
    
    clonal_prev = read.table(file.path(dir, "clonal_prevalence_table.tsv"), 
        sep = "\t", header = T)
    tree_edges = read.table(file.path(dir, "tree_edges.tsv"), sep = "\t", header = T)
    
    root <- tree_edges$source[!(tree_edges$source %in% tree_edges$target)]
    tree_edges <- subset(tree_edges, source != root)
    clonal_prev <- reformat_andrew_table(clonal_prev, db_path, patient_id)
    # clonal_prev <- subset(clonal_prev, long_sample_id %in% sample_ids)
    
    clonal_prev_short <- subset(clonal_prev, !duplicate, select = c(patient_id, 
        clone_id, clone_idx, condensed_id, site_id, clonal_prev))
    
    if (!dir.exists(outdir)) {
        dir.create(outdir)
        
        clone_file_output <- file.path(outdir, "clone_prevalences.tsv")
        tree_edges_output <- file.path(outdir, "tree_edges.tsv")
        intersect_output <- file.path(outdir, "intersect.tsv")
        
        write.table(clonal_prev_short, file = clone_file_output, sep = "\t", 
            quote = FALSE, row.names = FALSE, col.names = TRUE)
        write.table(tree_edges, file = tree_edges_output, sep = "\t", quote = FALSE, 
            row.names = FALSE, col.names = TRUE)
        write.table(cross_summary, file = intersect_output, sep = "\t", quote = FALSE, 
            row.names = FALSE, col.names = TRUE)
    }
}

Plots

htmltools::tagList(lapply(mapscape_input_directories, function(dir) {
    patient_id <- stringr::str_extract(dir, "[0-9]+$")
    df <- subset(intersect_table, patient1 == patient_id)
    clonal_prev = read.table(file.path(dir, "clonal_prevalence_table.tsv"), 
        sep = "\t", header = T)
    tree_edges = read.table(file.path(dir, "tree_edges.tsv"), sep = "\t", header = T)
    
    root <- tree_edges$source[!(tree_edges$source %in% tree_edges$target)]
    tree_edges <- subset(tree_edges, source != root)
    
    sample_locations = read.table(file.path(dir, "sample_locations.tsv"), sep = "\t", 
        header = T)
    sample_locations$long_sample_id <- df_as_map(clonal_prev, sample_locations$sample_id, 
        from = "sample_id", to = "long_sample_id")
    
    clonal_prev <- reformat_andrew_table(clonal_prev, db_path, patient_id)
    sample_locations <- reformat_andrew_table(sample_locations, db_path, patient_id)
    
    # df <- merge(df, subset(andrew_map, select=c(sample_key, sample_id)),
    # by.x=c('sample1'), by.y=c('sample_key')) %>%
    # plyr::rename(c('sample_id'='sample_id1')) df <- merge(df,
    # subset(andrew_map, select=c(sample_key, sample_id)), by.x=c('sample2'),
    # by.y=c('sample_key')) %>% plyr::rename(c('sample_id'='sample_id2'))
    
    clonotype_counts_sub <- clonotype_counts[as.character(clonotype_counts$patient_id) == 
        patient_id, ]
    
    patient_tmp <- patient_id
    
    clonotypes_sub <- xcr_table[(xcr_table$type == segment_type) & (xcr_table$patient_id == 
        as.numeric(patient_tmp)), ]
    clonotypes_sub$condensed_id <- factor(clonotypes_sub$condensed_id)
    cross_table <- ithi.xcr::cross_tabulate(clonotypes_sub, id_type = "condensed_id")
    
    # df$sample1 <- df$sample_id1 df$sample2 <- df$sample_id2
    
    sample_ids <- intersect(clonal_prev$long_sample_id, unique(df$sample1))
    
    if (length(sample_ids) == 0) {
        return(NULL)
    }
    
    df <- subset(df, (sample1 %in% sample_ids) & (sample2 %in% sample_ids))
    df$sample1 <- factor(df$sample1, levels = sample_ids)
    df$sample2 <- factor(df$sample2, levels = sample_ids)
    df <- df[order(df$sample1), ]
    
    ## IF YOU ONLY WANT THE INTERSECTION TO BE SHOWN, JUST SET CLONOTYPE TABLE TO
    ## NULL Otherwise set it to clonotype_counts_sub
    if (show_full_track) {
        res <- ithi.figures::format_xcrcircos_data(df, clonotype_counts = clonotype_counts_sub, 
            cross_table = cross_table)
    } else {
        res <- ithi.figures::format_xcrcircos_data(df, clonotype_counts = NULL)
    }
    chord_data <- res$chords
    track_data <- as.data.frame(res$tracks)
    
    resize_factor <- min(sum(track_data$len)/max(track_data$len) * 1/nrow(track_data), 
        0.9)
    
    clonal_prev <- subset(clonal_prev, long_sample_id %in% sample_ids)
    sample_locations <- subset(sample_locations, sample_id %in% clonal_prev$sample_id)
    
    if (!full_page_mode) {
        plot_title <- paste0("Patient ", patient_id)
    } else {
        plot_title <- ""
    }
    
    pvalue <- subset(stats, patient_id == patient_tmp)$p.value
    
    ## TODO: Extract pvalues from XCR-clone correlation
    
    print(patient_id)
    plot_xcrmapscape(clonal_prev, sample_locations, tree_edges, chord_data = chord_data, 
        track_data = track_data, resize_factor = resize_factor, plot_title = plot_title, 
        full_page_mode = full_page_mode, pvalue = round(pvalue, 3))
}))

[1] "1"
[1] "2"
[1] "3"
[1] "4"
[1] "7"
[1] "9"
[1] "10"
[1] "11"
[1] "12"
[1] "13"
[1] "14"
[1] "15"
[1] "16"
[1] "17"