# Load libraries
library(Seurat)
library(tidyverse)
library(circlize)
library(ComplexHeatmap)
library(paletteer)
library(viridis)
library(readxl)
library(patchwork)
library(muscat)
library(scater)10 Tumor: Plot TCR reconstruction data and T cell expression characteristics of vaccine-reactive TILs
10.1 Set up Seurat workspace
10.2 Read in clustered Seurat object
merged.18279.tumor.singlets <- readRDS("Tumor_scRNA_Part8.rds")10.3 Import and re-format T cell reactivity groupings
10.3.1 Open the scTCR barcodes in each category: Existing, Post-Nivo, Post-Vax, Reactive
p101_tils <- read.csv("sctcr_scRep_vjaa_p101_til_freq_barcodes_category_reactivity.csv")
p103_tils <- read.csv("sctcr_scRep_vjaa_p103_til_freq_barcodes_category_reactivity.csv")
p104_tils <- read.csv("sctcr_scRep_vjaa_p104_til_freq_barcodes_category_reactivity.csv")
p108_tils <- read.csv("sctcr_scRep_vjaa_p108_til_freq_barcodes_category_reactivity.csv")10.3.2 Reformat dataframes so that each scTCR barcode has its own row and filter for Post-Nivo and Post-Vax TCRs
p101_pnpv <- as_tibble(p101_tils) %>%
dplyr::filter(sctcr_category %in% c("Post-Nivolumab","Post-Vaccine")) %>%
mutate(Barcode = strsplit(as.character(barcodes), ",")) %>%
unnest(Barcode) %>%
dplyr::select(-barcodes) %>%
mutate(Barcode = str_replace(Barcode, "-1", ""))
p103_pnpv <- as_tibble(p103_tils) %>%
dplyr::filter(sctcr_category %in% c("Post-Nivolumab","Post-Vaccine")) %>%
mutate(Barcode = strsplit(as.character(barcodes), ",")) %>%
unnest(Barcode) %>%
dplyr::select(-barcodes) %>%
mutate(Barcode = str_replace(Barcode, "-1", ""))10.3.3 Merge into one dataframe
combined_pnpv <- bind_rows(p101_pnpv,p103_pnpv) %>%
dplyr::select(vjaa,Barcode,sctcr_category) %>%
dplyr::filter(Barcode %in% str_replace_all(colnames(merged.18279.tumor.singlets), "_.{1,3}mgIpi_RNA",""))
combined_pnpv# A tibble: 4,495 × 3
vjaa Barcode sctcr_category
<chr> <chr> <chr>
1 TRAV1-2.TRAJ12;CALHMDSSYKLIF_NA;NA P101_Tumor_W12_TTGGAA… Post-Nivolumab
2 TRAV13-1.TRAJ10;CAFLPGGGNKLTF_NA;NA P101_Tumor_W12_ATCACG… Post-Nivolumab
3 TRAV13-1.TRAJ49;CAASNWNTGNQFYF_NA;NA P101_Tumor_W12_CGGAGT… Post-Nivolumab
4 TRAV13-1.TRAJ6;CAAKGAGGSYIPTF_NA;NA P101_Tumor_W12_CTGAAA… Post-Nivolumab
5 TRAV13-2.TRAJ39;CAENPFIAGNMLTF_NA;NA P101_Tumor_W12_CGAGCC… Post-Nivolumab
6 TRAV17.TRAJ11;CATDAAGYSTLTF_NA;NA P101_Tumor_W12_TACTTA… Post-Nivolumab
7 TRAV19.TRAJ26;CALSEERDNYGQNFVF_NA;NA P101_Tumor_W12_CTCGTA… Post-Nivolumab
8 TRAV23DV6.TRAJ10;CAASILLTGGGNKLTF_NA;NA P101_Tumor_W12_AATCCA… Post-Nivolumab
9 TRAV26-1.TRAJ17;CIVKAAGNKLTF_NA;NA P101_Tumor_W12_TATCTC… Post-Nivolumab
10 TRAV29DV5.TRAJ33;CAASSDSNYQLIW_NA;NA P101_Tumor_W12_AACGTT… Post-Nivolumab
# ℹ 4,485 more rows10.3.4 Create pseudobulks by Patient and Timepoint, plot heatmap comparing Post-Nivo and Post-Vax clonotypes
cells.pnpv <- colnames(merged.18279.tumor.singlets)[str_replace_all(colnames(merged.18279.tumor.singlets), "_.{1,3}mgIpi_RNA","") %in% combined_pnpv$Barcode]
merged.18279.tumor.singlets.pnpv <- subset(merged.18279.tumor.singlets,
cells = cells.pnpv)
tcr_cats_pnpv <- combined_pnpv %>%
dplyr::select(-vjaa) %>%
deframe()
tcr_cats_pnpv_ordered <- tcr_cats_pnpv[str_replace_all(cells.pnpv, "_.{1,3}mgIpi_RNA","")]
merged.18279.tumor.singlets.pnpv <- AddMetaData(merged.18279.tumor.singlets.pnpv, as.character(tcr_cats_pnpv_ordered), col.name = "TCRgroup")
merged.18279.tumor.singlets.pnpv[['RNA']] <- JoinLayers(merged.18279.tumor.singlets.pnpv[['RNA']])
merged.pnpv.sce <- as.SingleCellExperiment(merged.18279.tumor.singlets.pnpv, assay="RNA")
(mergedCondition.pnpv.sce <- prepSCE(merged.pnpv.sce,
kid = "sub.cluster",
gid = "TCRgroup",
sid = "Patient",
drop = TRUE))class: SingleCellExperiment
dim: 61217 4495
metadata(1): experiment_info
assays(2): counts logcounts
rownames(61217): 5-8S-rRNA 5S-rRNA ... ZZEF1 ZZZ3
rowData names(0):
colnames(4495): P101_Tumor_W12_2.5mgIpi_RNA_TCCACACCAACACCTA
P101_Tumor_W12_2.5mgIpi_RNA_AGCTCTCTCAGCCTAA ...
P103_Tumor_W20_2.5mgIpi_RNA_AGCGTCGGTCATATGC
P103_Tumor_W20_2.5mgIpi_RNA_AAGGCAGTCTAGAGTC
colData names(3): cluster_id sample_id group_id
reducedDimNames(4): PCA UMAP.UNINTEGRATED INTEGRATED.HARMONY
UMAP.HARMONY
mainExpName: RNA
altExpNames(0):pb <- aggregateData(mergedCondition.pnpv.sce,
assay = "counts",
fun = "mean",
by = c("group_id", "sample_id"))
# Set up gene groupings for heatmap
c1 <- c("CD3E","CD4","CD8A","CD8B")
c2 <- c("CCR7","TCF7","SELL","CD28","CD27","IL7R","FAS","LEF1")
c3 <- c("NKG7","CST7","FASLG","PRF1","GZMA","GZMB","GZMH","GZMK","GZMM","GNLY","IL21","IL2","TNF","IFNG","CCL4","CCL5")
c4 <- c("CD38","ENTPD1","ITGAE","CD69","KLRG1","CD40LG","TNFRSF4","TNFRSF9","ICOS","ITGAL","ITGA1","ITGB1","CX3CR1","CXCR6","IL2RA")
c5 <- c("PDCD1","CTLA4","TIGIT","LAG3","HAVCR2","CD244","CD160","BTLA","VTCN1","TOX")
c6 <- c("MKI67","TBX1","EOMES","ZNF683","ID2","ID3","PRDM1","GATA3","FOXP3")
heat.goi <- list("T cell" = c1,
"Memory/Naive" = c2,
"Effector function" = c3,
"Phenotype" = c4,
"Inhibitory and exhaustion" = c5,
"Transcription factors" = c6
) %>%
stack() %>%
dplyr::rename("Gene" = values,"Category" = ind) %>%
as_tibble()
heat.goi# A tibble: 62 × 2
Gene Category
<chr> <fct>
1 CD3E T cell
2 CD4 T cell
3 CD8A T cell
4 CD8B T cell
5 CCR7 Memory/Naive
6 TCF7 Memory/Naive
7 SELL Memory/Naive
8 CD28 Memory/Naive
9 CD27 Memory/Naive
10 IL7R Memory/Naive
# ℹ 52 more rowspn.mat <- assays(pb)[['Post-Nivolumab']][heat.goi$Gene,]
colnames(pn.mat) <- paste0(colnames(pn.mat),"_PostNivo")
pv.mat <- assays(pb)[['Post-Vaccine']][heat.goi$Gene,]
colnames(pv.mat) <- paste0(colnames(pv.mat),"_PostVax")
ComplexHeatmap::Heatmap(cbind(pn.mat,pv.mat),
cluster_rows = F,
cluster_columns = F,
cluster_row_slices = F,
column_split = factor(c("P101","P103","P101","P103")),
column_order = c("P101_PostNivo","P101_PostVax","P103_PostNivo","P103_PostVax"),
row_split = factor(heat.goi$Category),
border = TRUE,
row_title_rot = 0,
row_title_side = "left",
row_names_side = "right",
col = circlize::colorRamp2(c(0,3),hcl_palette = "viridis"),
name = "Expression",
row_names_gp = gpar(fontsize=6,fontface="italic"),
row_title_gp = gpar(fontsize=8,fontface="bold")
)
10.3.5 Reformat dataframes so that each scTCR barcode has its own row and filter for reactive TCRs
p101_reactive <- as_tibble(p101_tils) %>%
dplyr::filter(reactive==TRUE) %>%
mutate(Barcode = strsplit(as.character(barcodes), ",")) %>%
unnest(Barcode) %>%
dplyr::select(-barcodes) %>%
mutate(Barcode = str_replace(Barcode, "-1", ""))
p103_reactive <- as_tibble(p103_tils) %>%
dplyr::filter(reactive==TRUE) %>%
mutate(Barcode = strsplit(as.character(barcodes), ",")) %>%
unnest(Barcode) %>%
dplyr::select(-barcodes) %>%
mutate(Barcode = str_replace(Barcode, "-1", ""))
p104_reactive <- as_tibble(p104_tils) %>%
dplyr::filter(reactive==TRUE) %>%
mutate(Barcode = strsplit(as.character(barcodes), ",")) %>%
unnest(Barcode) %>%
dplyr::select(-barcodes) %>%
mutate(Barcode = str_replace(Barcode, "-1", ""))
p108_reactive <- as_tibble(p108_tils) %>%
dplyr::filter(reactive==TRUE) %>%
mutate(Barcode = strsplit(as.character(barcodes), ",")) %>%
unnest(Barcode) %>%
dplyr::select(-barcodes) %>%
mutate(Barcode = str_replace(Barcode, "-1", ""))10.3.6 Merge into one dataframe
combined_reactive <- bind_rows(p101_reactive,p103_reactive,p104_reactive,p108_reactive) %>%
dplyr::select(vjaa,Barcode) %>%
dplyr::filter(Barcode %in% str_replace_all(colnames(merged.18279.tumor.singlets), "_.{1,3}mgIpi_RNA",""))
combined_reactive# A tibble: 15 × 2
vjaa Barcode
<chr> <chr>
1 TRAV21.TRAJ36;CAVRKQTGANNLFF_TRBV27.TRBJ2-1;CASSPTWGGSMVF P101_Tumor_…
2 TRAV16.TRAJ28;CALGSYQLTF_TRBV5-1.TRBJ2-7;CASSLGQGSIFSSYEQYF P101_Tumor_…
3 TRAV10.TRAJ27;CVVNTNAGKSTF_TRBV7-2.TRBJ2-1;CASSLGSGVAYNEQFF P101_Tumor_…
4 TRAV27.TRAJ22;CAVDSGSARQLTF_TRBV20-1.TRBJ2-7;CSASEGAHEQYF P103_Tumor_…
5 TRAV27.TRAJ22;CAVDSGSARQLTF_TRBV20-1.TRBJ2-7;CSASEGAHEQYF P103_Tumor_…
6 TRAV27.TRAJ22;CAVDSGSARQLTF_TRBV20-1.TRBJ2-7;CSASEGAHEQYF P103_Tumor_…
7 TRAV27.TRAJ22;CAVDSGSARQLTF_TRBV20-1.TRBJ2-7;CSASEGAHEQYF P103_Tumor_…
8 TRAV27.TRAJ22;CAVDSGSARQLTF_TRBV20-1.TRBJ2-7;CSASEGAHEQYF P103_Tumor_…
9 TRAV41.TRAJ22;CAVLLISSGSARQLTF_TRBV10-3.TRBJ2-2;CAIGPRKGGNTGELFF P104_Tumor_…
10 TRAV26-1.TRAJ23;CIGVLYNQGGKLIF_TRBV5-1.TRBJ1-2;CASSLESGQEDYGYTF P104_Tumor_…
11 TRAV9-2.TRAJ58;CALIEETSGSRLTF_TRBV5-5.TRBJ1-2;CASSTWTDYGYTF P104_Tumor_…
12 TRAV8-6.TRAJ43;CAVRDNNNDMRF_TRBV6-1.TRBJ2-1;CASRVAGGPYEQFF P104_Tumor_…
13 TRAV1-1.TRAJ3;CAVPGGGSASKIIF_TRBV19.TRBJ2-5;CASSMGPGGAETQYF P108_Tumor_…
14 TRAV38-1.TRAJ31;CAFMKQRNNARLMF_TRBV6-6.TRBJ1-2;CASSYSTRDTNYGYTF P108_Tumor_…
15 TRAV12-2.TRAJ22;CAVVRSGSARQLTF_TRBV7-6.TRBJ2-1;CASSQGLRNEQFF P108_Tumor_…10.4 Import data with numbered TCRs
numbered_tcrs_reformatted <- read.csv("reformatted_reconstructed_TCR.csv")
numbered_tcrs_original <- read_xlsx("reconstructed_TCR.xlsx")10.4.1 Check the order of the rows is the same then merge
table(numbered_tcrs_reformatted[,"CDR3A_1"] == numbered_tcrs_original[,"CDR3A_1"])
TRUE
96 # Merge dataframes so we have the TCR# and the scTCR_TRBV1 column in the same file
numbered_tcrs <- numbered_tcrs_original %>%
cbind(scTCR_TRBV_1 = numbered_tcrs_reformatted$scTCR_TRBV_1)
# Add vjaa column
numbered_tcrs <- numbered_tcrs %>%
mutate(vjaa = paste0(TRAV_1, ".", TRAJ_1, ";", CDR3A_1, "_", scTCR_TRBV_1, ".", TRBJ_1, ";", CDR3B_1)) %>%
select(`TCR#`, vjaa)
combined_reactive_numbered <- left_join(combined_reactive,numbered_tcrs,by="vjaa") %>%
mutate(TCRnum = paste0("TCR",`TCR#`)) %>%
dplyr::select(-`TCR#`)
combined_reactive_numbered# A tibble: 15 × 3
vjaa Barcode TCRnum
<chr> <chr> <chr>
1 TRAV21.TRAJ36;CAVRKQTGANNLFF_TRBV27.TRBJ2-1;CASSPTWGGSMVF P101_T… TCR17
2 TRAV16.TRAJ28;CALGSYQLTF_TRBV5-1.TRBJ2-7;CASSLGQGSIFSSYEQYF P101_T… TCR7
3 TRAV10.TRAJ27;CVVNTNAGKSTF_TRBV7-2.TRBJ2-1;CASSLGSGVAYNEQFF P101_T… TCR19
4 TRAV27.TRAJ22;CAVDSGSARQLTF_TRBV20-1.TRBJ2-7;CSASEGAHEQYF P103_T… TCR27
5 TRAV27.TRAJ22;CAVDSGSARQLTF_TRBV20-1.TRBJ2-7;CSASEGAHEQYF P103_T… TCR27
6 TRAV27.TRAJ22;CAVDSGSARQLTF_TRBV20-1.TRBJ2-7;CSASEGAHEQYF P103_T… TCR27
7 TRAV27.TRAJ22;CAVDSGSARQLTF_TRBV20-1.TRBJ2-7;CSASEGAHEQYF P103_T… TCR27
8 TRAV27.TRAJ22;CAVDSGSARQLTF_TRBV20-1.TRBJ2-7;CSASEGAHEQYF P103_T… TCR27
9 TRAV41.TRAJ22;CAVLLISSGSARQLTF_TRBV10-3.TRBJ2-2;CAIGPRKGGNTGE… P104_T… TCR53
10 TRAV26-1.TRAJ23;CIGVLYNQGGKLIF_TRBV5-1.TRBJ1-2;CASSLESGQEDYGY… P104_T… TCR40
11 TRAV9-2.TRAJ58;CALIEETSGSRLTF_TRBV5-5.TRBJ1-2;CASSTWTDYGYTF P104_T… TCR58
12 TRAV8-6.TRAJ43;CAVRDNNNDMRF_TRBV6-1.TRBJ2-1;CASRVAGGPYEQFF P104_T… TCR43
13 TRAV1-1.TRAJ3;CAVPGGGSASKIIF_TRBV19.TRBJ2-5;CASSMGPGGAETQYF P108_T… TCR88
14 TRAV38-1.TRAJ31;CAFMKQRNNARLMF_TRBV6-6.TRBJ1-2;CASSYSTRDTNYGY… P108_T… TCR87
15 TRAV12-2.TRAJ22;CAVVRSGSARQLTF_TRBV7-6.TRBJ2-1;CASSQGLRNEQFF P108_T… TCR94 10.5 Plot TCRs ranked by post-vaccine frequencies
Label with whether it was selected for reconstruction, and which of those were vaccine-reactive
10.5.1 Import scTCR object
combined_TCR <- readRDS("Skin_Tumor_scTCR.rds")10.5.2 Rank TCRs by post-vaccine frequency
reconstructed_tcr <- numbered_tcrs_reformatted %>%
mutate(TRAJ_1 = str_replace_all(TRAJ_1,";","")) %>%
mutate(vjaa = paste0(TRAV_1, ".", TRAJ_1, ";", CDR3A_1, "_", scTCR_TRBV_1, ".", TRBJ_1, ";", CDR3B_1)) %>%
mutate(Patient = paste0("P",Patient)) %>%
dplyr::select(Patient,vjaa,Reactive) %>%
distinct()
postvax_clonefreq_ranked <- bind_rows(combined_TCR) %>%
as_tibble() %>%
dplyr::filter(Site == "Tumor") %>%
dplyr::filter(Timepoint %in% c("W20", "PD")) %>%
mutate(Timepoint = str_replace_all(Timepoint, "W20|PD", "Post-vaccine")) %>%
group_by(Patient,Timepoint,vjaa) %>%
mutate(nDistinctClonesPerTimepoint = n()) %>%
ungroup() %>%
left_join(reconstructed_tcr, by=c("vjaa","Patient")) %>%
dplyr::select(Patient,Timepoint,CTaa,vjaa,nDistinctClonesPerTimepoint,Reactive) %>%
distinct() %>%
mutate(Reconstructed = case_when(
is.na(Reactive) ~ "No",
!is.na(Reactive) ~ "Yes"
)
) %>%
mutate(ReactivityCategory = case_when(
Reactive == "Yes" ~ "Reactive",
Reactive == "No" ~ "Not reactive",
Reconstructed == "No" ~ "Not tested"
)
) %>%
group_by(Patient) %>%
dplyr::arrange(Patient, -nDistinctClonesPerTimepoint) %>%
mutate(Rank = row_number())
postvax_clonefreq_ranked_list <- postvax_clonefreq_ranked %>%
group_by(Patient) %>%
group_split()
clone_lengths <- unlist(lapply(postvax_clonefreq_ranked_list,nrow))
clone_rel_lengths <- clone_lengths/min(clone_lengths)10.5.3 Plot and assemble heatmap with annotations
p101 <- ComplexHeatmap::Heatmap(as.factor(postvax_clonefreq_ranked_list[[1]]$Reconstructed),
cluster_rows = FALSE,
cluster_columns = FALSE,
border = TRUE,
name = "Reconstructed",
column_label = "P101\n",
column_names_side = "top",
column_names_rot = 0,
column_names_centered = TRUE,
column_names_gp = gpar(fontface = "bold"),
row_title = "Clones ranked by post-vaccine frequency",
row_title_gp = gpar(fontsize = 9),
col = c("gray80","blue"),
left_annotation = rowAnnotation(NumClones = anno_barplot(postvax_clonefreq_ranked_list[[1]]$nDistinctClonesPerTimepoint, axis_param = list(direction = "reverse")),width = unit(2, "cm"), show_legend = FALSE),
right_annotation = rowAnnotation(Reactive = as.factor(postvax_clonefreq_ranked_list[[1]]$Reactive),
col = list(Reactive = c("No" = "white", "Yes" = "red")),
na_col = "white",
border = TRUE,
show_legend = FALSE,
annotation_name_rot = 0
),
width = unit(2, "cm"),
show_heatmap_legend = FALSE
) %>%
draw() %>%
grid.grabExpr()
p103 <- ComplexHeatmap::Heatmap(as.factor(postvax_clonefreq_ranked_list[[2]]$Reconstructed),
cluster_rows = FALSE,
cluster_columns = FALSE,
border = TRUE,
name = "Reconstructed",
column_label = "P103\n",
column_names_side = "top",
column_names_rot = 0,
column_names_centered = TRUE,
column_names_gp = gpar(fontface = "bold"),
col = c("gray80","blue"),
left_annotation = rowAnnotation(NumClones = anno_barplot(postvax_clonefreq_ranked_list[[2]]$nDistinctClonesPerTimepoint, axis_param = list(direction = "reverse")),width = unit(2, "cm"), show_legend = FALSE),
right_annotation = rowAnnotation(Reactive = as.factor(postvax_clonefreq_ranked_list[[2]]$Reactive),
col = list(Reactive = c("No" = "white", "Yes" = "red")),
na_col = "white",
border = TRUE,
show_legend = FALSE,
annotation_name_rot = 0
),
width = unit(2, "cm"),
show_heatmap_legend = FALSE
) %>%
draw() %>%
grid.grabExpr()
p104 <- ComplexHeatmap::Heatmap(as.factor(postvax_clonefreq_ranked_list[[3]]$Reconstructed),
cluster_rows = FALSE,
cluster_columns = FALSE,
border = TRUE,
name = "Reconstructed",
column_label = "P104\n",
column_names_side = "top",
column_names_rot = 0,
column_names_centered = TRUE,
column_names_gp = gpar(fontface = "bold"),
col = c("gray80","blue"),
left_annotation = rowAnnotation(NumClones = anno_barplot(postvax_clonefreq_ranked_list[[3]]$nDistinctClonesPerTimepoint, axis_param = list(direction = "reverse")),width = unit(2, "cm"), show_legend = FALSE),
right_annotation = rowAnnotation(Reactive = as.factor(postvax_clonefreq_ranked_list[[3]]$Reactive),
col = list(Reactive = c("No" = "white", "Yes" = "red")),
na_col = "white",
border = TRUE,
show_legend = FALSE,
annotation_name_rot = 0
),
width = unit(2, "cm"),
show_heatmap_legend = FALSE
) %>%
draw() %>%
grid.grabExpr()
p108 <- ComplexHeatmap::Heatmap(as.factor(postvax_clonefreq_ranked_list[[4]]$Reconstructed),
cluster_rows = FALSE,
cluster_columns = FALSE,
border = TRUE,
name = "Reconstructed",
column_label = "P108\n",
column_names_side = "top",
column_names_rot = 0,
column_names_centered = TRUE,
column_names_gp = gpar(fontface = "bold"),
col = c("gray80","blue"),
left_annotation = rowAnnotation(NumClones = anno_barplot(postvax_clonefreq_ranked_list[[4]]$nDistinctClonesPerTimepoint, axis_param = list(direction = "reverse")),width = unit(2, "cm")),
right_annotation = rowAnnotation(Reactive = as.factor(postvax_clonefreq_ranked_list[[4]]$Reactive),
col = list(Reactive = c("No" = "white", "Yes" = "red")),
na_col = "white",
border = TRUE,
annotation_name_rot = 0
),
width = unit(2, "cm")
) %>%
draw() %>%
grid.grabExpr()
wrap_plots(p101, p103, p104, p108, ncol = 4)
10.6 Plot TILs with reactive TCRs as heatmap
10.6.1 Set up gene groupings for TIL heatmap
c1 <- c("CD3E","CD4","CD8A","CD8B")
c2 <- c("CCR7","TCF7","SELL","CD28","CD27","IL7R","FAS","LEF1")
c3 <- c("NKG7","CST7","FASLG","PRF1","GZMA","GZMB","GZMH","GZMK","GZMM","GNLY","IL21","IL2","TNF","IFNG","CCL4","CCL5")
c4 <- c("CD38","ENTPD1","ITGAE","CD69","KLRG1","CD40LG","TNFRSF4","TNFRSF9","ICOS","ITGAL","ITGA1","ITGB1","CX3CR1","CXCR6","IL2RA")
c5 <- c("PDCD1","CTLA4","TIGIT","LAG3","HAVCR2","CD244","CD160","BTLA","VTCN1","TOX")
c6 <- c("MKI67","TBX1","EOMES","ZNF683","ID2","ID3","PRDM1","GATA3","FOXP3")
heat.goi <- list("T cell" = c1,
"Memory/Naive" = c2,
"Effector function" = c3,
"Phenotype" = c4,
"Inhibitory and exhaustion" = c5,
"Transcription factors" = c6
) %>%
stack() %>%
dplyr::rename("Gene" = values,"Category" = ind) %>%
as_tibble()
heat.goi# A tibble: 62 × 2
Gene Category
<chr> <fct>
1 CD3E T cell
2 CD4 T cell
3 CD8A T cell
4 CD8B T cell
5 CCR7 Memory/Naive
6 TCF7 Memory/Naive
7 SELL Memory/Naive
8 CD28 Memory/Naive
9 CD27 Memory/Naive
10 IL7R Memory/Naive
# ℹ 52 more rows10.6.2 Add reactive TCR data as metadata label to Seurat object
reac.meta <- ifelse(str_replace_all(colnames(merged.18279.tumor.singlets), "_.{1,3}mgIpi_RNA","") %in% combined_reactive_numbered$Barcode,"Reactive","Other")
merged.18279.tumor.singlets <- AddMetaData(merged.18279.tumor.singlets, reac.meta, col.name="ReactiveTCR")10.6.3 Add subcluster annotations as metadata label to Seurat object
cluster_annot <- c(
"8" = "Macrophage",
"15" = "MonoMac",
"28" = "pDC",
"17" = "B_cells",
"13" = "Plasma_cells",
"2_1" = "CD8_TPEX",
"1_0" = "CD8_TEX",
"1_1" = "CD8_T_EM_1",
"2_0" = "CD8_T_EM_2",
"22_1" = "CD8_T_EM_3",
"22_2" = "CD8_T_EM_4",
"26_0" = "CD8_T_proliferating",
"3_0" = "CD4_T_Memory_1",
"22_0" = "CD4_T_Memory_2",
"3_1" = "CD4_T_CTL_1",
"12_0" = "CD4_T_CTL_2",
"12_1" = "CD4_Tfh",
"16_0" = "Treg",
"9_0" = "T_unknown_1",
"9_1" = "T_unknown_2",
"18_0" = "NK/T_1",
"18_1" = "NK/T_2",
"29" = "Mast",
"0" = "Melanoma_1",
"4" = "Melanoma_2",
"5" = "Melanoma_3",
"6" = "Melanoma_4",
"11" = "Melanoma_5",
"23" = "Melanoma_6",
"25" = "Melanoma_7",
"14" = "Melanoma_proliferating",
"21" = "Melanocytes",
"7" = "Endothelial_1",
"27" = "Endothelial_2",
"10" = "Fibroblasts_1",
"20" = "Fibroblasts_2",
"19" = "Pericytes/SmoothMuscle",
"24" = "Basal"
)
clustLabels <- as.data.frame(cluster_annot)[merged.18279.tumor.singlets$sub.cluster,]
merged.18279.tumor.singlets <- AddMetaData(merged.18279.tumor.singlets, list(clustLabels), col.name = "CellAnnotation")
merged.18279.tumor.singlets$CellAnnotation <- factor(x = merged.18279.tumor.singlets$CellAnnotation,
levels = as.character(cluster_annot))
# Set up broader cell classes
cluster_classes <- enframe(cluster_annot,name = "sub.cluster",value = "CellAnnotation") %>%
as_tibble() %>%
mutate(CellClass = case_when(
sub.cluster %in% c("8","15") ~ "Monocyte/\nMacrophage",
sub.cluster %in% c("28") ~ "DC",
sub.cluster %in% c("17","13") ~ "B",
sub.cluster %in% c("1_0","1_1","2_0","2_1","22_1","22_2","26_0","3_0","3_1","12_0","12_1","22_0","16_0","9_0","9_1","18_0","18_1") ~ "NK/T",
sub.cluster %in% c("29") ~ "Mast",
sub.cluster %in% c("0","4","5","6","11","23","25","14") ~ "Melanoma",
sub.cluster %in% c("21","7","27","10","20","19","24") ~ "Non-immune"
)
) %>%
dplyr::select(-CellAnnotation) %>%
column_to_rownames(var = "sub.cluster") %>%
as.data.frame()
clustClasses <- cluster_classes[merged.18279.tumor.singlets$sub.cluster,]
merged.18279.tumor.singlets <- AddMetaData(merged.18279.tumor.singlets, list(clustClasses), col.name = "CellClass")10.6.4 Subset Seurat object to NK/T clusters
merged.18279.tumor.singlets.t <- subset(merged.18279.tumor.singlets, subset = CellClass == "NK/T")10.6.5 Make metadata table of reactive TILs for heatmap annotations
meta.df <- merged.18279.tumor.singlets.t@meta.data[,c("Patient","Timepoint","sub.cluster","CellAnnotation","CellClass","functional.cluster","ReactiveTCR")]
meta.reac.df <- meta.df[meta.df$ReactiveTCR=="Reactive",]
meta.reac <- full_join(
rownames_to_column(meta.reac.df,var="Barcode") %>%
mutate(Barcode2 = str_replace_all(Barcode, "_.{1,3}mgIpi_RNA","")),
combined_reactive_numbered,
by=c("Barcode2" = "Barcode")
) %>%
dplyr::select(-Barcode2) %>%
column_to_rownames(var="Barcode") %>%
as.data.frame()
rownames_to_column(meta.reac,var="Barcode") Barcode Patient Timepoint sub.cluster
1 P101_Tumor_W20_2.5mgIpi_RNA_GCGACCACAGACAGGT P101 W20 12_1
2 P101_Tumor_W20_2.5mgIpi_RNA_AAACCTGCACCGAAAG P101 W20 3_0
3 P101_Tumor_W20_2.5mgIpi_RNA_TGGCGCAAGTGTGGCA P101 W20 22_0
4 P103_Tumor_W20_2.5mgIpi_RNA_CTAGTGAAGACGACGT P103 W20 2_1
5 P103_Tumor_W20_2.5mgIpi_RNA_TTATGCTCATCATCCC P103 W20 2_0
6 P103_Tumor_W20_2.5mgIpi_RNA_ATTGGACGTGTGTGCC P103 W20 2_0
7 P103_Tumor_W20_2.5mgIpi_RNA_CATCGGGCAGCATACT P103 W20 2_0
8 P103_Tumor_W20_2.5mgIpi_RNA_CTGCGGAAGTAGCGGT P103 W20 1_0
9 P104_Tumor_PD_2.5mgIpi_RNA_TCGCGAGCAATGTAAG P104 PD 3_0
10 P104_Tumor_PD_2.5mgIpi_RNA_CATATTCGTCCTCCAT P104 PD 12_1
11 P104_Tumor_PD_2.5mgIpi_RNA_GATCGTAAGCGACGTA P104 PD 12_1
12 P104_Tumor_PD_2.5mgIpi_RNA_ACACCGGAGCCGCCTA P104 PD 3_0
13 P108_Tumor_PD_5mgIpi_RNA_ACGGGTCCATAGAAAC P108 PD 3_0
14 P108_Tumor_PD_5mgIpi_RNA_TTTGCGCAGCTAGCCC P108 PD 1_0
15 P108_Tumor_PD_5mgIpi_RNA_ACTTGTTCAGTCAGAG P108 PD 3_0
CellAnnotation CellClass functional.cluster ReactiveTCR
1 CD4_Tfh NK/T CD4.Tfh Reactive
2 CD4_T_Memory_1 NK/T CD4.CTL_EOMES Reactive
3 CD4_T_Memory_2 NK/T CD4.Memory Reactive
4 CD8_TPEX NK/T CD8.TPEX Reactive
5 CD8_T_EM_2 NK/T CD8.TPEX Reactive
6 CD8_T_EM_2 NK/T CD8.TPEX Reactive
7 CD8_T_EM_2 NK/T CD8.TPEX Reactive
8 CD8_TEX NK/T CD8.TEX Reactive
9 CD4_T_Memory_1 NK/T CD4.Tfh Reactive
10 CD4_Tfh NK/T CD4.Tfh Reactive
11 CD4_Tfh NK/T CD4.Tfh Reactive
12 CD4_T_Memory_1 NK/T CD4.Tfh Reactive
13 CD4_T_Memory_1 NK/T CD4.Memory Reactive
14 CD8_TEX NK/T CD4.CTL_EOMES Reactive
15 CD4_T_Memory_1 NK/T CD4.Memory Reactive
vjaa TCRnum
1 TRAV16.TRAJ28;CALGSYQLTF_TRBV5-1.TRBJ2-7;CASSLGQGSIFSSYEQYF TCR7
2 TRAV21.TRAJ36;CAVRKQTGANNLFF_TRBV27.TRBJ2-1;CASSPTWGGSMVF TCR17
3 TRAV10.TRAJ27;CVVNTNAGKSTF_TRBV7-2.TRBJ2-1;CASSLGSGVAYNEQFF TCR19
4 TRAV27.TRAJ22;CAVDSGSARQLTF_TRBV20-1.TRBJ2-7;CSASEGAHEQYF TCR27
5 TRAV27.TRAJ22;CAVDSGSARQLTF_TRBV20-1.TRBJ2-7;CSASEGAHEQYF TCR27
6 TRAV27.TRAJ22;CAVDSGSARQLTF_TRBV20-1.TRBJ2-7;CSASEGAHEQYF TCR27
7 TRAV27.TRAJ22;CAVDSGSARQLTF_TRBV20-1.TRBJ2-7;CSASEGAHEQYF TCR27
8 TRAV27.TRAJ22;CAVDSGSARQLTF_TRBV20-1.TRBJ2-7;CSASEGAHEQYF TCR27
9 TRAV26-1.TRAJ23;CIGVLYNQGGKLIF_TRBV5-1.TRBJ1-2;CASSLESGQEDYGYTF TCR40
10 TRAV8-6.TRAJ43;CAVRDNNNDMRF_TRBV6-1.TRBJ2-1;CASRVAGGPYEQFF TCR43
11 TRAV41.TRAJ22;CAVLLISSGSARQLTF_TRBV10-3.TRBJ2-2;CAIGPRKGGNTGELFF TCR53
12 TRAV9-2.TRAJ58;CALIEETSGSRLTF_TRBV5-5.TRBJ1-2;CASSTWTDYGYTF TCR58
13 TRAV12-2.TRAJ22;CAVVRSGSARQLTF_TRBV7-6.TRBJ2-1;CASSQGLRNEQFF TCR94
14 TRAV38-1.TRAJ31;CAFMKQRNNARLMF_TRBV6-6.TRBJ1-2;CASSYSTRDTNYGYTF TCR87
15 TRAV1-1.TRAJ3;CAVPGGGSASKIIF_TRBV19.TRBJ2-5;CASSMGPGGAETQYF TCR8810.6.6 Set up color coding vectors for heatmap annotations
patient_cols <- c("P101" = "#59A14FFF",
"P103" = "#B07AA1FF",
"P104" = "#76B7B2FF",
"P108" = "#F6AAC9FF")
projectils_cols <- data.frame(cbind(
c(levels(factor(merged.18279.tumor.singlets$functional.cluster)),"Unknown/Other"),
c(paletteer::paletteer_d("RColorBrewer::Blues")[-c(1:2)],
"#232061FF",
paletteer::paletteer_d("RColorBrewer::Reds")[-c(1:2)],
"#CECECEFF")
)
) %>%
deframe()10.6.7 Subset data matrices to just the genes for plotting, and only reactive TILs
rna.reactive <- as.matrix(merged.18279.tumor.singlets.t@assays$RNA$data[heat.goi$Gene,rownames(meta.reac)])10.6.8 Set column ordering
reactive.order <- rownames_to_column(meta.reac,var="Cell") %>%
as_tibble() %>%
mutate(Patient = factor(Patient, c("P101","P104","P108","P103"))) %>%
dplyr::arrange(Patient) %>%
pull(Cell)10.6.9 Set up heatmap annotations
ha1 <- HeatmapAnnotation(Patient = as.factor(meta.reac$Patient),
col = list(Patient = patient_cols[as.factor(meta.reac$Patient)]), border = T)
ha2 <- HeatmapAnnotation(CellState = as.factor(meta.reac$functional.cluster),
col = list(CellState = projectils_cols[meta.reac$functional.cluster]), border = T)
ha3 <- HeatmapAnnotation(TCR = anno_text(meta.reac$TCRnum,
gp = gpar(fontsize=6),
just = "left",
location = 0))10.6.10 Capture heatmap in plotted order and print TCR info
h0 <- ComplexHeatmap::Heatmap(rna.reactive,
cluster_rows = F,
cluster_columns = F,
cluster_row_slices = F,
row_split = factor(heat.goi$Category),
column_split = factor(meta.reac$functional.cluster),
column_order = reactive.order,
top_annotation = c(ha3,ha1,ha2),
border = TRUE,
row_names_gp = gpar(fontsize=6,fontface="italic"),
row_title_gp = gpar(fontsize=8,fontface="bold"),
column_title_gp = gpar(fontsize=8,fontface="bold"),
show_column_names = F,
row_title_rot = 0,
row_title_side = "left",
row_names_side = "right",
col = circlize::colorRamp2(c(0,2), c("gray90", "red")),
name = "Expression",
column_title = "Individual TILs with Reactive TCR",
use_raster = F,
width = unit(6, "cm")
)
ht <- draw(h0)
col_order <- column_order(ht)
heatmap_barcodes_ordered <- colnames(rna.reactive)[as.numeric(unlist(col_order))]png
2 10.7 Get session info
sessionInfo()R version 4.3.2 (2023-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Rocky Linux 8.10 (Green Obsidian)
Matrix products: default
BLAS/LAPACK: /usr/lib64/libopenblasp-r0.3.15.so; LAPACK version 3.9.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
time zone: America/New_York
tzcode source: system (glibc)
attached base packages:
[1] stats4 grid stats graphics grDevices utils datasets
[8] methods base
other attached packages:
[1] scater_1.34.0 scuttle_1.12.0
[3] SingleCellExperiment_1.24.0 SummarizedExperiment_1.32.0
[5] Biobase_2.62.0 GenomicRanges_1.54.1
[7] GenomeInfoDb_1.38.8 IRanges_2.36.0
[9] S4Vectors_0.40.2 BiocGenerics_0.48.1
[11] MatrixGenerics_1.14.0 matrixStats_1.5.0
[13] muscat_1.14.0 patchwork_1.3.0
[15] readxl_1.4.5 viridis_0.6.5
[17] viridisLite_0.4.2 paletteer_1.6.0
[19] ComplexHeatmap_2.18.0 circlize_0.4.16
[21] lubridate_1.9.4 forcats_1.0.0
[23] stringr_1.5.1 dplyr_1.1.4
[25] purrr_1.0.4 readr_2.1.5
[27] tidyr_1.3.1 tibble_3.2.1
[29] ggplot2_3.5.1 tidyverse_2.0.0
[31] Seurat_5.1.0 SeuratObject_5.0.2
[33] sp_2.2-0
loaded via a namespace (and not attached):
[1] spatstat.sparse_3.1-0 bitops_1.0-9
[3] httr_1.4.7 RColorBrewer_1.1-3
[5] doParallel_1.0.17 numDeriv_2016.8-1.1
[7] backports_1.5.0 tools_4.3.2
[9] sctransform_0.4.1 utf8_1.2.4
[11] R6_2.6.1 lazyeval_0.2.2
[13] uwot_0.2.3 mgcv_1.9-3
[15] GetoptLong_1.0.5 withr_3.0.2
[17] prettyunits_1.2.0 gridExtra_2.3
[19] SeuratWrappers_0.3.19 progressr_0.15.1
[21] cli_3.6.4 Cairo_1.6-2
[23] spatstat.explore_3.2-6 fastDummies_1.7.3
[25] sandwich_3.1-1 prismatic_1.1.2
[27] mvtnorm_1.3-3 spatstat.data_3.1-6
[29] blme_1.0-5 ggridges_0.5.6
[31] pbapply_1.7-2 R.utils_2.13.0
[33] parallelly_1.43.0 limma_3.58.1
[35] generics_0.1.3 shape_1.4.6.1
[37] gtools_3.9.5 ica_1.0-3
[39] spatstat.random_3.2-2 Matrix_1.6-4
[41] ggbeeswarm_0.7.2 abind_1.4-8
[43] R.methodsS3_1.8.2 lifecycle_1.0.4
[45] multcomp_1.4-28 yaml_2.3.10
[47] edgeR_4.0.16 gplots_3.2.0
[49] SparseArray_1.2.2 Rtsne_0.17
[51] promises_1.3.2 crayon_1.5.3
[53] miniUI_0.1.1.1 lattice_0.22-7
[55] beachmat_2.18.1 cowplot_1.1.3
[57] magick_2.8.6 pillar_1.10.1
[59] knitr_1.45 rjson_0.2.23
[61] boot_1.3-31 estimability_1.5.1
[63] future.apply_1.11.3 codetools_0.2-20
[65] leiden_0.4.3.1 glue_1.8.0
[67] remotes_2.5.0 data.table_1.17.0
[69] vctrs_0.6.5 png_0.1-8
[71] spam_2.11-1 Rdpack_2.6.3
[73] cellranger_1.1.0 gtable_0.3.6
[75] rematch2_2.1.2 xfun_0.52
[77] rbibutils_2.3 S4Arrays_1.2.0
[79] mime_0.13 coda_0.19-4.1
[81] reformulas_0.4.0 survival_3.8-3
[83] iterators_1.0.14 statmod_1.5.0
[85] fitdistrplus_1.2-2 TH.data_1.1-3
[87] ROCR_1.0-11 nlme_3.1-168
[89] pbkrtest_0.5.3 EnvStats_2.8.1
[91] progress_1.2.3 RcppAnnoy_0.0.22
[93] TMB_1.9.17 irlba_2.3.5.1
[95] vipor_0.4.7 KernSmooth_2.23-26
[97] colorspace_2.1-1 nnet_7.3-20
[99] DESeq2_1.42.1 tidyselect_1.2.1
[101] emmeans_1.11.0 compiler_4.3.2
[103] BiocNeighbors_1.20.2 DelayedArray_0.28.0
[105] plotly_4.10.4 caTools_1.18.3
[107] scales_1.3.0 remaCor_0.0.18
[109] lmtest_0.9-40 digest_0.6.37
[111] goftest_1.2-3 spatstat.utils_3.1-3
[113] minqa_1.2.8 variancePartition_1.30.2
[115] rmarkdown_2.29 aod_1.3.3
[117] RhpcBLASctl_0.23-42 XVector_0.42.0
[119] htmltools_0.5.8.1 pkgconfig_2.0.3
[121] lme4_1.1-37 sparseMatrixStats_1.14.0
[123] fastmap_1.2.0 rlang_1.1.5
[125] GlobalOptions_0.1.2 htmlwidgets_1.6.4
[127] shiny_1.10.0 DelayedMatrixStats_1.24.0
[129] farver_2.1.2 zoo_1.8-13
[131] jsonlite_2.0.0 BiocParallel_1.36.0
[133] R.oo_1.27.0 BiocSingular_1.18.0
[135] RCurl_1.98-1.17 magrittr_2.0.3
[137] modeltools_0.2-23 GenomeInfoDbData_1.2.11
[139] dotCall64_1.2 munsell_0.5.1
[141] Rcpp_1.0.14 reticulate_1.35.0
[143] stringi_1.8.7 zlibbioc_1.48.2
[145] MASS_7.3-60.0.1 flexmix_2.3-20
[147] plyr_1.8.9 parallel_4.3.2
[149] listenv_0.9.1 ggrepel_0.9.6
[151] deldir_2.0-4 splines_4.3.2
[153] tensor_1.5 hms_1.1.3
[155] locfit_1.5-9.12 igraph_2.1.4
[157] spatstat.geom_3.2-8 RcppHNSW_0.6.0
[159] reshape2_1.4.4 ScaledMatrix_1.10.0
[161] evaluate_1.0.3 BiocManager_1.30.25
[163] nloptr_2.0.3 tzdb_0.5.0
[165] foreach_1.5.2 httpuv_1.6.15
[167] RANN_2.6.2 polyclip_1.10-7
[169] future_1.34.0 clue_0.3-66
[171] scattermore_1.2 rsvd_1.0.5
[173] broom_1.0.8 xtable_1.8-4
[175] RSpectra_0.16-2 later_1.4.1
[177] lmerTest_3.1-3 glmmTMB_1.1.10
[179] beeswarm_0.4.0 cluster_2.1.8.1
[181] timechange_0.3.0 globals_0.16.3