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)

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")

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

plot_xcrmapscape <- function(clonal_prev, sample_locations, tree_edges, chord_data, 
    track_data, resize_factor) {
    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'
    
    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)
}

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

Read 19.7% of 304912 rows
Read 62.3% of 304912 rows
Read 88.6% of 304912 rows
Read 304912 rows and 18 (of 18) columns from 0.070 GB file in 00:00:05

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

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

intersect_table <- trb_intersect_table
clonotype_table <- subset(xcr_table, type == "TRB")
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)
}

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)
    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 == "TRB") & (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)
    
    print(patient_id)
    plot_xcrmapscape(clonal_prev, sample_locations, tree_edges, chord_data = chord_data, 
        track_data = track_data, resize_factor = resize_factor)
}))

[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"