-
Notifications
You must be signed in to change notification settings - Fork 1
/
7f_bulk_hnsc_single_cell_validation.Rmd
129 lines (106 loc) · 3.75 KB
/
7f_bulk_hnsc_single_cell_validation.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
---
title: "Figure 6F: Enrichment of TCGA-HNSC deconvolution results in single cell data"
output:
html_document:
df_print: paged
self_contained: yes
---
```{r}
source("../R/figure_utils.R")
source("../R/setup.R")
library(Seurat)
```
```{r}
save_dir <- "../out/dualsimplex_save_tcga_hnscc_4kg_9ct_no_hk_neighbours/"
custom_cols <- c(
"#1f78b4", "#33a02c", "#e31a1c", "#ff7f00", "#6a3d9a",
"#b15928", "#a6cee3", "#b2df8a", "#fb9a99", "#fdbf6f"
)
```
```{r}
lo2 <- DualSimplexSolver$from_state(save_dir, save_there = F)
```
Markers for enrichment in single cell (100 per cell type)
```{r fig.width = 14, fig.height = 6}
lo2$plot_projected("markers", use_dims = NULL, with_legend = T, with_history = F)
```
```{r fig.width=6.5, fig.height=5}
seurat_obj <- readRDS("../data/large/GSE181919_integrated.rds")
so_malignant <- readRDS("../data/large/GSE181919_malignant_integrated.rds")
so_epithelial <- readRDS("../data/large/GSE181919_epithelial_integrated.rds")
plt <- Seurat::DimPlot(
seurat_obj,
group.by = "cell.type",
reduction = "tsne",
cols = custom_cols,
raster = TRUE
) + Seurat::NoAxes() + ggtitle("GSE181919 dataset")
show(plt)
Seurat::DimPlot(so_malignant, group.by = "cell.type", reduction = "umap")
Seurat::DimPlot(so_epithelial, group.by = "cell.type", reduction = "umap")
ggsave("../out/6f_gse181919_dataset.svg", plt, width = 6.5, height = 5, device = svg)
```
```{r fig.width=20, fig.height=15}
seurat_obj <- add_list_markers(seurat_obj, lo2$get_marker_genes())
so_malignant <- add_list_markers(so_malignant, lo2$get_marker_genes())
so_epithelial <- add_list_markers(so_epithelial, lo2$get_marker_genes())
plot_marker_enrichment(seurat_obj, lo2$get_ct_names(), limits = NULL)
plot_marker_enrichment(so_malignant, lo2$get_ct_names(), limits = NULL)
plot_marker_enrichment(so_epithelial, lo2$get_ct_names(), limits = NULL)
```
```{r fig.width=10, fig.height=15}
cancer_cts <- c("cell_type_4", "cell_type_5", "cell_type_6", "cell_type_8", "cell_type_9")
plt <- plot_marker_enrichment(
so_malignant,
cancer_cts,
limits = NULL,
ncol = 2,
ggadd = function(x, y) { x + Seurat::NoAxes() + Seurat::NoLegend() },
raster = T
)
ggsave("../out/6f_malignant_subset.svg", plt, width = 10, height = 15, device = svg)
plt
```
```{r fig.width=4.7, fig.height=4}
plt <- Seurat::DimPlot(so_malignant, group.by = "hpv", raster = T) + Seurat::NoAxes() + ggtitle("")
ggsave("../out/6f_malignant_subset_hpv.svg", plt, width = 4.7, height = 4, device = svg)
plt
```
```{r}
table(seurat_obj$cell.type)
```
```{r, fig.width=7, fig.height=4}
enrichment_map <- as.data.frame([email protected] %>% group_by(cell.type) %>% summarize(
across(paste0(lo2$get_ct_names(), 1), list(function(x) {
if ((cur_column() == "cell_type_51")) {
tc <- 1000
} else {
tc <- take_top_cells
}
mean(tail(sort(x), tc))
})), .groups = "drop"
))
rownames(enrichment_map) <- enrichment_map$cell.type
enrichment_map$cell.type <- NULL
colnames(enrichment_map) <- substr(colnames(enrichment_map), 1, nchar(colnames(enrichment_map)) - 3)
col_order <- c(
"cell_type_1", "cell_type_2", "cell_type_3", "cell_type_7", cancer_cts
)
row_order <- c(
"Macrophages", "Myocytes", "Epithelial.cells", "Fibroblasts", "Malignant.cells"
)
enrichment_map <- enrichment_map[row_order, col_order]
#enrichment_map_norm <- as.data.frame(t(apply(enrichment_map, 1, function(x) (x - min(x)) / (max(x) - min(x)))))
enrichment_map_norm <- as.data.frame(apply(enrichment_map, 2, function(x) (x - min(x)) / (max(x) - min(x))))
plt <- pheatmap::pheatmap(
enrichment_map_norm,
display_numbers = round(enrichment_map, 2),
cluster_rows = F,
cluster_cols = F,
treeheight_row = 0,
treeheight_col = 0,
legend = F
)
ggsave("../out/6f_heatmap.svg", plt, width = 7, height = 4, device = svglite::svglite)
plt
```