---
title: "Figure 2 and S3-S5 -- tumor cell intrinsic phenotypic shifts"
author: "Matthew Zatzman"
execute:
warning: false
error: false
format:
html:
code-fold: true
code-tools: true
toc: true
---
Figure 2 and the supplementary figures its legend points to: S3 (cohort cell type
composition), S4 (malignant cell identification) and S5 (data integration).
Run `00_download_convert.qmd` first to build `paper_data/`.
{{< include _setup.qmd >}}
# Load data
Only the cohort object is read. S4A additionally needs the inferCNV object; the
epithelial, tumor and histotime objects are not used by any panel here.
```{r load}
# L'Ecuyer-CMRG so forked mclapply workers get reproducible streams.
RNGkind("L'Ecuyer-CMRG")
set.seed(42)
db <- load_project(objects = "cohort", cnv = TRUE)
srt <- db$srt
df_sp <- db$spec
# Sample labels get a leading star when the sample contributed no tumour cells.
srt$sample_id_has_tumor <- srt$sample_id
no_tum <- !levels(srt$sample_id_has_tumor) %in% unique(as.character(srt$sample_id[srt$sample_has_tumor]))
levels(srt$sample_id_has_tumor)[no_tum] <- paste0("*", levels(srt$sample_id_has_tumor)[no_tum])
```
# Figure 2 -- Tumor cell intrinsic phenotypic shifts accompany disease progression
## Figure 2A-D -- UMAP before batch correction
```{r fig2ad-cohort-umap-combined-nobatch, fig.width = 14, fig.height = 8}
# Patient
p1 <- DimPlot(srt, reduction = "umap_nb", group.by = "site_of_tissue_simple") +
umap_theme() +
scale_color_manual(values = colors$site_of_tissue_simple, breaks = names(colors$site_of_tissue_simple), labels = label_num_generator(srt$site_of_tissue_simple)) +
NoAxes() +
labs(title = "Tissue site") +
theme(plot.margin = margin(1, 1, 1, 1, "mm"))
p2 <- DimPlot(srt, reduction = "umap_nb", group.by = "sample_id") +
umap_theme() +
scale_color_manual(values = colors$sample_id, breaks = names(colors$sample_id), labels = label_num_generator(srt$sample_id)) +
guides(color = guide_legend(ncol = 3, override.aes = list(size = 3), keywidth = 0, keyheight = 0)) +
NoAxes() +
labs(title = "Sample") +
theme(plot.margin = margin(8, 1, 1, 1, "mm"))
p3 <- DimPlot(srt, reduction = "umap_nb", group.by = "time_point") +
umap_theme() +
NoAxes() +
scale_color_manual(values = colors$time_point, breaks = names(colors$time_point), labels = label_num_generator(srt$time_point)) +
labs(title = "Treatment state") +
theme(plot.margin = margin(1, 1, 1, 1, "mm"))
# Cell types
p4 <- DimPlot(srt, group.by = "cell_type", reduction = "umap_nb") +
umap_theme() +
NoAxes() +
scale_color_manual(values = colors$cell_type, breaks = names(colors$cell_type), labels = label_num_generator(srt$cell_type)) +
labs(title = "Cell Type") +
theme(plot.margin = margin(1, 1, 1, 1, "mm"))
pcomb <- p4 + p2 + p1 + p3 + plot_layout(ncol = 2) & theme(aspect.ratio = 1, legend.box.margin = margin(), legend.margin = margin())
pcomb
ggsave(pcomb, file = here("plots", "cohort_umap_combined_nobatch.pdf"), width = 14, height = 8)
```
## Figure 2E -- major cell type proportions across treatment states
```{r fig2e-cohort-timepoint-cell-type-prop, fig.width = 4.5, fig.height = 2.5}
p <- srt@meta.data %>%
filter(sample_has_tumor) %>%
add_count(sample_id, name = "sample_n") %>%
add_count(sample_id, cell_type, name = "cell_type_n") %>%
mutate(prop = cell_type_n / sample_n) %>%
distinct(sample_id, cell_type, time_point, sample_n, cell_type_n, prop) %>%
filter(cell_type != "Other") %>%
ggplot(aes(x = cell_type, y = prop)) +
geom_boxplot(aes(fill = time_point), outlier.size = 0.25) +
scale_fill_manual(values = colors$time_point) +
labs(x = "Cell Type", y = "Proportion", fill = "Treatment state") +
scale_y_continuous(breaks = seq(0, 1, by = 0.2)) +
guides(
x = guide_axis(cap = "both", angle = 45),
y = guide_axis(cap = "both")
) +
ggpubr::geom_pwc(tip.length = 0, symnum.args = symnum_0.1_star, label = "p.adj.signif", aes(group = time_point), hide.ns = T, vjust = 0.5, method = "wilcox_test", step.increase = 0.05) +
coord_cartesian(clip = "off") #+
# theme(legend.position = "top", plot.margin = margin(0, 1, 1, 1, "mm"), legend.title.position = "top", legend.title = element_text(hjust = 0.5), legend.box.margin = margin(), legend.box.spacing = unit(0, "mm"), legend.box.background = element_blank(), legend.background = element_blank())
p
ggsave(p, file = here("plots", "cohort_timepoint_cell_type_prop.pdf"), width = 4.5, height = 2.5)
```
## Figure 2F and 2H -- patient and treatment-state specificity scores
```{r derive-specificity-panels}
cls_order <- levels(fct_reorder(df_sp$patient_nocor$cid2, df_sp$patient_nocor$cid1_purity, median))
sp_order <- c("patient_nocor", "patient_cor", "time_point_nocor", "time_point_cor",
"site_of_tissue_simple_nocor", "site_of_tissue_simple_cor")
spec_pls <- map(sp_order, .f = function(l) {
df <- df_sp[[l]]
type <- unique(df$type)
type <- case_when(
type == "patient" ~ "Patient",
type == "time_point" ~ "Treatment state",
type == "site_of_tissue_simple" ~ "Tissue site"
)
bcor <- ifelse(unique(df$batch_cor) == "cor", "Batch Corrected", "Uncorrected")
p1 <- df %>%
filter(cid2 != "Other") %>%
# ggplot(aes(x = fct_relevel(cid2, cls_order), y = cid1_purity)) +
ggplot(aes(x = fct_reorder(cid2, cid1_purity, median), y = cid1_purity)) +
geom_violin(scale = "width", aes(fill = cid2), show.legend = F) +
scale_y_continuous(breaks = seq(0, 1, by = 0.5)) +
geom_pointrange(
stat = "summary", fatten = 0.8,
fun.min = function(z) {
quantile(z, 0.25)
},
fun.max = function(z) {
quantile(z, 0.75)
},
fun = median
) +
guides(
x = guide_axis(angle = 45, cap = "both"),
y = guide_axis(cap = "both")
) +
scale_fill_manual(values = colors$cell_type) +
labs(x = "Cell Type", y = "Specificity", title = glue("{type}: {bcor}")) +
ggpubr::stat_compare_means(comparisons = list(c("Tumor Cell", "Epithelial")), tip.length = 0.05, label.y = 1.05, method = "wilcox", size = 3, vjust = -0.5) +
coord_cartesian(clip = "off")
})
names(spec_pls) <- sp_order
```
```{r fig2f2h-sample-timepoint-specificity, fig.width = 2.85, fig.height = 3}
p1 <- spec_pls$patient_nocor + labs(title = NULL, y = "Patient\nspecificity\n(uncorrected)")
# Reviewer #1 C7: axis reads "treatment state", matching the Figure 2 legend.
p2 <- spec_pls$time_point_cor + labs(title = NULL, y = "Treatment state\nspecificity\n(corrected)")
pcomb <- p1 + p2 + plot_layout(ncol = 1, axes = "collect") & labs(x = NULL) & guides(x = guide_axis(angle = 45, cap = "both"))
pcomb
ggsave(pcomb, file = here("plots", "sample_timepoint_specificity.pdf"), width = 2.85, height = 3)
```
## Figure 2G and 2I -- UMAP after batch correction
```{r fig2g2i-cohort-umap-combined-batch-main2, fig.width = 9, fig.height = 3}
p3 <- DimPlot(srt, reduction = "umap", group.by = "time_point") +
umap_theme() +
NoAxes() +
scale_color_manual(values = colors$time_point, breaks = names(colors$time_point)) +
labs(title = "Treatment state")
# Cell types
p4 <- DimPlot(srt, group.by = "cell_type", reduction = "umap") +
umap_theme() +
NoAxes() +
scale_color_manual(values = colors$cell_type, breaks = names(colors$cell_type)) +
labs(title = "Cell Type")
#
cohort_umap_tumor_cells <- DimPlot(srt, group.by = "is_tumor_cell", cols = list("TRUE" = as.character(colors$cell_type["Tumor Cell"]), "FALSE" = "grey90")) +
labs(title = NULL) +
# theme_void() +
theme(aspect.ratio = 1, legend.key.spacing = unit(0.1, "lines")) +
umap_theme() +
NoAxes() +
NoLegend() +
labs(title = "Tumor Cells")
# cohort_umap_tumor_cells
pcomb <- wrap_plots(p4, p3, byrow = F) + plot_layout(nrow = 1) & theme(aspect.ratio = 1, plot.margin = margin(1, 1, 1, 1, "mm"))
pcomb
ggsave(pcomb, file = here("plots", "cohort_umap_combined_batch_main2.pdf"), width = 8, height = 3)
```
## Figure 2J -- tumor cell density by treatment state
```{r fig2j-tumor-cell-density, fig.width=7, fig.height=3}
library(ggtext)
df_raw <- FetchData(srt, vars = c("time_point", "is_tumor_cell", "cell_type", "umap_1", "umap_2")) %>%
dplyr::rename("UMAP1" = "umap_1", "UMAP2" = "umap_2")
df <- df_raw %>%
filter(is_tumor_cell)
# Group by timepoint and calculate scaled density per group
df <- df %>%
group_by(time_point, cell_type) %>%
mutate(density = get_density(UMAP1, UMAP2, n = 500)) %>%
mutate(density_scaled = egfrnucseq::standardize(density))
df$time_point_col <- factor(df$time_point,
levels = levels(df$time_point),
labels = glue("<b style='color:{colors$time_point[levels(df$time_point)]}'>{levels(df$time_point)}</b>")
)
df_raw_background <- df_raw %>%
filter(is_tumor_cell)
# lab_cents <- get_label_centers(srt, group_var = "RNA_cluster_labels_coarse", reduced_dim = "RNA_umap_rpca_integrated")
p <- df %>%
ungroup() %>%
arrange((density_scaled)) %>%
ggplot(aes(x = UMAP1, y = UMAP2)) +
scattermore::geom_scattermore(data = df_raw_background, color = "grey50", alpha = 0.5) +
scattermore::geom_scattermore(aes(color = density_scaled), pointsize = 2, pixels = c(512, 512)) +
# ggrastr::rasterize(geom_point(aes(color = density_scaled), size = 0.01), dpi = 600) +
facet_wrap(~time_point_col) +
scale_color_gradientn(colours = viridis::turbo(10)) +
guides(color = guide_colorbar(keywidth = 0.5, keyheight = 3)) +
coord_cartesian(clip = "on", xlim = c(-10, 6),
ylim = c(-12.1, 4)) +
theme(
strip.background = element_blank(),
strip.text.x = element_markdown(size = 16),
aspect.ratio = 1, panel.background = element_rect(fill = NA, color = "black"),
plot.title = element_text(hjust = 0.5, size = 18),
axis.line = element_blank(),
axis.ticks = element_blank(),
axis.text = element_blank()
) +
labs(x = "UMAP1", y = "UMAP2", title = NULL, color = "Tumor Cell\nDensity")
ggsave(p, file = here("plots", "tumor_cell_density.pdf"), width = 7, height = 3)
p
```
# Figure S3 -- Cohort cell type composition
## Figure S3A -- cell type marker dot plot
```{r figs3a-cell-type-markers-dotplot, fig.width = 8, fig.height = 3}
ct_markers <- markers$cell_type_markers
names(ct_markers) <- map(names(ct_markers), .f = function(ct) {
ct_col <- colors$cell_type[ct]
ct_colname <- glue("<b style='color:{ct_col}'>{ct}</b>")
return(ct_colname)
}) %>%
list_c()
library(ggtext)
srt$cell_type_col <- factor(srt$cell_type, levels = levels(srt$cell_type), labels = glue("<b style='color:{colors$cell_type[levels(srt$cell_type)]}'>{levels(srt$cell_type)}</b>"))
# srt$cell_type_col <- fct_rev(srt$cell_type_col)
pl <- DotPlotter(srt, features = ct_markers, group.by = "cell_type_col", cluster_groups = F, rotate_y_strip = F) +
theme(
strip.clip = "off",
axis.line = element_blank(),
strip.text.x = element_markdown(),
axis.text.y = element_markdown()
)
ggsave(pl, file = here("plots", "cell_type_markers_dotplot.pdf"), width = 8, height = 3)
pl
```
## Figure S3B -- cell counts and proportions per specimen
```{r figs3b-sample-cellcounts-barplot, fig.width=12, fig.height = 4}
p1 <- srt@meta.data %>%
ggplot(aes(x = sample_id_has_tumor)) +
geom_bar(aes(fill = cell_type)) +
facet_grid(. ~ time_point, scales = "free", space = "free") +
guides(x = guide_axis(angle = 90)) +
scale_fill_manual(values = colors$cell_type) +
labs(fill = "Cell Type", y = "Count", x = "Sample") +
guides(fill = guide_legend(keywidth = 0.4, keyheight = 0.2)) +
theme_clean()
p2 <- srt@meta.data %>%
ggplot(aes(x = sample_id_has_tumor)) +
geom_bar(aes(fill = cell_type), position = "fill") +
facet_grid(. ~ time_point, scales = "free", space = "free") +
guides(
x = guide_axis(angle = 90),
y = guide_axis(cap = "none")
) +
scale_fill_manual(values = colors$cell_type) +
labs(fill = "Cell Type", y = "Proportion", x = "Sample") +
guides(fill = guide_legend(keywidth = 0.4, keyheight = 0.2)) +
theme_clean()
# Final formatting
p1 <- p1 + theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) + labs(x = NULL)
p2 <- p2 + theme(strip.text.x = element_blank())
pcomb <- p1 + p2 +
plot_layout(guides = "collect", axes = "collect", ncol = 1) &
theme_clean() &
theme(plot.background = element_blank(),
strip.background = element_blank()) &
scale_y_continuous(expand = c(0, 0))
pcomb
ggsave(pcomb, file = here("plots", "sample_cellcounts_barplot.pdf"), width = 12, height = 4)
```
# Figure S4 -- Malignant cell identification
## Figure S4A -- inferred per-cell copy number heatmap
The panel below is a raster of ~100,000 cells squeezed onto the page's pixel grid.
`plots/cohort_cnv_heatmap.pdf`, written by the chunk, is the full-resolution version.
```{r, fig.width = 13, fig.height=8}
sce_full <- db$cnv
sce_full$is_tumor_cell <- srt$is_tumor_cell[match(colnames(sce_full), srt$cell_id)]
sce <- sce_full[, sce_full$sample_id %in% unique(srt$sample_id) &
sce_full$cell_type %in% c("Tumor Cell", "Epithelial")]
library(ComplexHeatmap)
# Sample 1000 tumor and 1000 normal cells from each patient
mat <- assay(sce) %>%
t() %>%
log2()
df_anno <- colData(sce) %>%
as.data.frame() %>%
dplyr::select(cell_type, time_point, sample_id)
df_anno$cell_type <- factor(df_anno$cell_type)
df_anno$time_point <- factor(df_anno$time_point, levels = c("TN", "MRD", "PD"))
df_anno$sample_id <- fct_relevel(df_anno$sample_id, levels(srt$sample_id))
df_anno <- df_anno %>%
arrange(time_point, sample_id, cell_type)
```
```{r figs4a-cohort-cnv-heatmap, fig.width = 13, fig.height=8, dpi = 600}
# 600 dpi: ComplexHeatmap rasters the body at raster_quality x 72 = 1800 dpi, and a PNG
# device keeps only what its own grid holds. The deposited PDF is the full-resolution panel.
mat <- mat[rownames(df_anno), ]
# Key by column name, never position: ComplexHeatmap matches `annotation_label`
# positionally, and `palette[factor]` indexes by integer level code, not label.
anno_cols <- list(
"cell_type" = colors$cell_type[levels(fct_drop(df_anno$cell_type))],
"time_point" = colors$time_point,
"sample_id" = colors$sample_id
)
anno_labels <- c(
"cell_type" = "Cell Type",
"time_point" = "Treatment state",
"sample_id" = "Sample"
)
stopifnot(
all(colnames(df_anno) %in% names(anno_labels)),
all(colnames(df_anno) %in% names(anno_cols)),
!anyNA(anno_cols$cell_type)
)
left_anno <- HeatmapAnnotation(
df = df_anno,
annotation_label = unname(anno_labels[colnames(df_anno)]),
which = "row",
show_legend = TRUE,
col = anno_cols[colnames(df_anno)],
annotation_legend_param = list(
"sample_id" = list(
ncol = 3
)
)
)
# Split columns by chromosome
chrs <- as.vector(gsub("chr", "", GenomeInfoDb::seqnames(SummarizedExperiment::rowRanges(sce))))
col_split <- factor(chrs, levels = unique(gtools::mixedsort(chrs)))
ht <- Heatmap(
matrix = mat,
name = "Log2Ratio",
left_annotation = left_anno,
col = logr_col_fun(breaks = c(-0.05, -0.01, 0.01, 0.05), colors = c("blue", "white", "white", "red")),
cluster_rows = F,
cluster_columns = F,
show_column_dend = F,
show_row_dend = F,
show_row_names = F,
show_column_names = F,
column_split = col_split,
row_split = df_anno$cell_type,
raster_by_magick = F, raster_device = "CairoPNG",
border = T,
raster_quality = 25
)
pdf(here("plots/cohort_cnv_heatmap.pdf"), width = 13, height = 8)
draw(ht)
dev.off()
# draw() renders to the active device only, and the pdf() above took the first one, so
# the panel needs a second draw to reach the page.
draw(ht)
```
## Figure S4B -- tumor cell proportion per specimen
```{r figs4b-tumor-cell-proportion-barplot, fig.width = 11, fig.height = 3}
df <- srt@meta.data %>%
add_count(sample_id_has_tumor, cell_type, name = "n_ct") %>%
add_count(sample_id_has_tumor, name = "n_samp") %>%
mutate(ct_prop = n_ct / n_samp) %>%
distinct(sample_id_has_tumor, cell_type, ct_prop, n_ct, n_samp, procedure_type, time_point, site_of_tissue_simple) %>%
pivot_wider(id_cols = c(sample_id_has_tumor, time_point, site_of_tissue_simple, procedure_type), names_from = cell_type, values_from = ct_prop, values_fill = 0) %>%
pivot_longer(cols = -c(sample_id_has_tumor, time_point, site_of_tissue_simple, procedure_type), names_to = "cell_type", values_to = "cell_type_prop")
p <- df %>%
filter(cell_type == "Tumor Cell") %>%
ggplot(aes(x = fct_reorder(sample_id_has_tumor, -cell_type_prop, mean), y = cell_type_prop)) +
geom_col(aes(fill = time_point), show.legend = F) +
scale_fill_manual(values = colors$time_point) +
facet_grid(. ~ time_point, scales = "free", space = "free") +
theme_clean() +
theme(strip.background = element_blank()) +
coord_cartesian(ylim = c(0, 1), clip = "off") +
scale_y_continuous(expand = c(0, 0)) +
guides(
x = guide_axis(angle = 90),
y = guide_axis(cap = "none")
) +
labs(x = "Sample", y = "Tumor Cell Proportion")
p2 <- df %>%
filter(cell_type == "Tumor Cell") %>%
ggplot(aes(x = time_point, y = cell_type_prop)) +
geom_boxplot(aes(fill = time_point), show.legend = F, outlier.size = 0.25) +
scale_fill_manual(values = colors$time_point) +
theme_clean() +
scale_y_continuous(expand = c(0, 0), breaks = seq(0, 1, by = 0.25)) +
guides(
x = guide_axis(angle = 45, cap = "none"),
y = guide_axis(cap = "both")
) +
labs(x = "Treatment state", y = "Tumor Cell Proportion") +
coord_cartesian(clip = "off", ylim = c(0, 1)) +
ggpubr::geom_pwc(tip.length = 0, label = "p={p.adj}", symnum.args = symnum_0.1_star, hide.ns = F, vjust = 0.2)
pcomb <- p + free(p2, type = "label", side = "b") + plot_layout(ncol = 2, widths = c(1, 0.1), axes = "collect") &
theme(plot.background = element_blank())
pcomb
ggsave(pcomb, file = here("plots", "tumor_cell_proportion_barplot.pdf"), width = 11, height = 3)
```
## Figure S4C-D -- tumor cellularity by treatment site
```{r figs4cd-tumor-cellularity-site-procedure, fig.width = 3.5, fig.height = 3}
p1 <- df %>%
filter(cell_type == "Tumor Cell") %>%
ggplot(aes(x = fct_relevel(site_of_tissue_simple, names(colors$site_of_tissue_simple)), y = cell_type_prop)) +
geom_boxplot(aes(fill = site_of_tissue_simple), show.legend = F) +
geom_quasirandom(alpha = 0.5) +
scale_fill_manual(values = colors$site_of_tissue_simple) +
theme_clean() +
scale_y_continuous(expand = c(0, 0), breaks = seq(0, 1, by = 0.25)) +
guides(
x = guide_axis(angle = 0, cap = "none"),
y = guide_axis(cap = "both")
) +
labs(x = "Tissue Site", y = "Tumor Cell Proportion") +
coord_cartesian(clip = "off", ylim = c(0, 1))
pcomb <- p1 + plot_layout(ncol = 1) &
theme_classic() &
ggpubr::stat_compare_means(label.x.npc = 0.1) &
guides(x = guide_axis(angle = 45))
pcomb
ggsave(pcomb, file = here("plots", "tumor_cellularity_site_procedure.pdf"), width = 3.5, height = 3)
```
# Figure S5 -- snRNA data integration
## Figure S5A-B -- specificity scores before and after integration
`spec_pls` holds all six specificity violins. Figure S5A-B is the full grid; Figure 2F/2H
picks two of the six, so the list is built once here rather than inside either panel.
```{r figs5ab-all-specificity-plots, fig.width = 8, fig.height = 5}
pcomb <- wrap_plots(spec_pls, ncol = 3, byrow = F) &
theme(
axis.title.y = element_text(hjust = 0.5),
plot.title = element_text(hjust = 0.5)
)
ggsave(pcomb, file = here("plots", "all_specificity_plots.pdf"), width = 8, height = 5)
pcomb
```
## Figure S5C-D -- batch-corrected UMAP by patient and site
```{r figs5cd-cohort-umap-combined-batch-supp, fig.width = 14, fig.height = 7}
# Patient
p1 <- DimPlot(srt, reduction = "umap", group.by = "site_of_tissue_simple", shuffle = T) +
umap_theme() +
scale_color_manual(values = colors$site_of_tissue_simple, breaks = names(colors$site_of_tissue_simple), labels = label_num_generator(srt$site_of_tissue_simple)) +
NoAxes() +
labs(title = "Tissue site")
p2 <- DimPlot(srt, reduction = "umap", group.by = "sample_id") +
umap_theme() +
scale_color_manual(values = colors$sample_id, breaks = names(colors$sample_id), labels = label_num_generator(srt$sample_id)) +
guides(color = guide_legend(ncol = 3, override.aes = list(size = 3), keywidth = 0, keyheight = 0)) +
NoAxes() +
labs(title = "Sample")
p3 <- DimPlot(srt, reduction = "umap", group.by = "time_point") +
umap_theme() +
NoAxes() +
scale_color_manual(values = colors$time_point, breaks = names(colors$time_point), labels = label_num_generator(srt$time_point)) +
labs(title = "Treatment state")
# Cell types
p4 <- DimPlot(srt, group.by = "cell_type", reduction = "umap") +
umap_theme() +
NoAxes() +
scale_color_manual(values = colors$cell_type, breaks = names(colors$cell_type), labels = label_num_generator(srt$cell_type)) +
labs(title = "Cell Type")
pcomb <- wrap_plots(p4, p2, p1, p3) & theme(aspect.ratio = 1)
pcomb
# For supplement: patient and tissue site
pcomb2 <- wrap_plots(p2, p1, byrow = F) + plot_layout(nrow = 1) & theme(aspect.ratio = 1, plot.margin = margin(1, 1, 1, 1, "mm"))
pcomb2
ggsave(pcomb2, file = here("plots", "cohort_umap_combined_batch_supp.pdf"), width = 15, height = 5)
```
