diff --git a/code/03_build_sce/08_cluster_annotation.R b/code/03_build_sce/08_cluster_annotation.R index 8083fda..985c70d 100644 --- a/code/03_build_sce/08_cluster_annotation.R +++ b/code/03_build_sce/08_cluster_annotation.R @@ -303,16 +303,16 @@ my_plotMarkers( #### Assign annotation #### order_cell_types <- function(cell_types) { - cell_types <- unique(cell_types) - neun_mask <- grepl("Excit|Inhib", cell_types) - drop_mask <- grepl("Ambiguous|drop", cell_types) - - neun_ct <- cell_types[neun_mask] - glia_ct <- cell_types[!(neun_mask|drop_mask)] - drop_ct <- cell_types[drop_mask] - - ordered_cell_types <- c(sort(glia_ct), sort(neun_ct), sort(drop_ct)) - return(ordered_cell_types) + cell_types <- unique(cell_types) + neun_mask <- grepl("Excit|Inhib", cell_types) + drop_mask <- grepl("Ambiguous|drop", cell_types) + + neun_ct <- cell_types[neun_mask] + glia_ct <- cell_types[!(neun_mask | drop_mask)] + drop_ct <- cell_types[drop_mask] + + ordered_cell_types <- c(sort(glia_ct), sort(neun_ct), sort(drop_ct)) + return(ordered_cell_types) } ## mbkm annotation @@ -328,23 +328,25 @@ anno_k2 <- anno_k %>% ct = paste0(broad, "_", f_pad_zero(row_number(), width = 2)), n = n() ) %>% - ungroup() %>% - mutate(cellType = ifelse(n > 1, ct, broad), - cellType = factor(cellType, levels = order_cell_types(cellType)), - broad = factor(broad, levels = order_cell_types(broad))) + ungroup() %>% + mutate( + cellType = ifelse(n > 1, ct, broad), + cellType = factor(cellType, levels = order_cell_types(cellType)), + broad = factor(broad, levels = order_cell_types(broad)) + ) levels(anno_k2$broad) -# [1] "Astro" "EndoMural" "MicroOligo" "Oligo" "OPC" "Excit" "Inhib" "drop" +# [1] "Astro" "EndoMural" "MicroOligo" "Oligo" "OPC" "Excit" "Inhib" "drop" levels(anno_k2$cellType) -# [1] "Astro" "EndoMural" "MicroOligo" "Oligo_01" "Oligo_02" "OPC" "Excit_01" "Excit_02" -# [9] "Excit_03" "Excit_04" "Excit_05" "Excit_06" "Excit_07" "Excit_08" "Excit_09" "Excit_10" -# [17] "Excit_11" "Excit_12" "Excit_13" "Inhib_01" "Inhib_02" "Inhib_03" "Inhib_04" "Inhib_05" -# [25] "Inhib_06" "drop_01" "drop_02" "drop_03" "drop_04" +# [1] "Astro" "EndoMural" "MicroOligo" "Oligo_01" "Oligo_02" "OPC" "Excit_01" "Excit_02" +# [9] "Excit_03" "Excit_04" "Excit_05" "Excit_06" "Excit_07" "Excit_08" "Excit_09" "Excit_10" +# [17] "Excit_11" "Excit_12" "Excit_13" "Inhib_01" "Inhib_02" "Inhib_03" "Inhib_04" "Inhib_05" +# [25] "Inhib_06" "drop_01" "drop_02" "drop_03" "drop_04" ## HC annotation anno_hc <- read.csv(here("processed-data", "03_build_sce", "DLPFC_HC_anno.csv")) table(anno_hc$broad) -# Ambiguous Astro EndoMural Excit Inhib Micro Oligo OPC +# Ambiguous Astro EndoMural Excit Inhib Micro Oligo OPC # 1 1 2 15 6 1 3 1 anno_hc2 <- anno_hc %>% @@ -354,17 +356,19 @@ anno_hc2 <- anno_hc %>% n = n() ) %>% ungroup() %>% - mutate(cellType = ifelse(n > 1, ct, broad), - cellType = factor(cellType, levels = order_cell_types(cellType)), - broad = factor(broad, levels = order_cell_types(broad))) + mutate( + cellType = ifelse(n > 1, ct, broad), + cellType = factor(cellType, levels = order_cell_types(cellType)), + broad = factor(broad, levels = order_cell_types(broad)) + ) levels(anno_hc2$broad) # [1] "Astro" "EndoMural" "Micro" "Oligo" "OPC" "Excit" "Inhib" "Ambiguous" levels(anno_hc2$cellType) -# [1] "Astro" "EndoMural_01" "EndoMural_02" "Micro" "Oligo_01" "Oligo_02" "Oligo_03" -# [8] "OPC" "Excit_01" "Excit_02" "Excit_03" "Excit_04" "Excit_05" "Excit_06" -# [15] "Excit_07" "Excit_08" "Excit_09" "Excit_10" "Excit_11" "Excit_12" "Excit_13" -# [22] "Excit_14" "Excit_15" "Inhib_01" "Inhib_02" "Inhib_03" "Inhib_04" "Inhib_05" +# [1] "Astro" "EndoMural_01" "EndoMural_02" "Micro" "Oligo_01" "Oligo_02" "Oligo_03" +# [8] "OPC" "Excit_01" "Excit_02" "Excit_03" "Excit_04" "Excit_05" "Excit_06" +# [15] "Excit_07" "Excit_08" "Excit_09" "Excit_10" "Excit_11" "Excit_12" "Excit_13" +# [22] "Excit_14" "Excit_15" "Inhib_01" "Inhib_02" "Inhib_03" "Inhib_04" "Inhib_05" # [29] "Inhib_06" "Ambiguous" ## Assign to sce @@ -372,43 +376,43 @@ sce$cellType_broad_k <- anno_k2$broad[match(sce$kmeans, anno_k$cluster)] sce$cellType_k <- anno_k2$cellType[match(sce$kmeans, anno_k2$cluster)] table(sce$cellType_broad_k) -# Astro EndoMural MicroOligo Oligo OPC Excit Inhib drop +# Astro EndoMural MicroOligo Oligo OPC Excit Inhib drop # 3557 1330 5541 33716 1791 21233 10413 23 table(sce$cellType_k) -# Astro EndoMural MicroOligo Oligo_01 Oligo_02 OPC Excit_01 Excit_02 Excit_03 Excit_04 -# 3557 1330 5541 28085 5631 1791 5879 1058 548 1188 -# Excit_05 Excit_06 Excit_07 Excit_08 Excit_09 Excit_10 Excit_11 Excit_12 Excit_13 Inhib_01 -# 82 633 3198 2851 896 1581 1732 1292 295 132 -# Inhib_02 Inhib_03 Inhib_04 Inhib_05 Inhib_06 drop_01 drop_02 drop_03 drop_04 -# 461 2242 6018 1311 249 11 3 2 7 +# Astro EndoMural MicroOligo Oligo_01 Oligo_02 OPC Excit_01 Excit_02 Excit_03 Excit_04 +# 3557 1330 5541 28085 5631 1791 5879 1058 548 1188 +# Excit_05 Excit_06 Excit_07 Excit_08 Excit_09 Excit_10 Excit_11 Excit_12 Excit_13 Inhib_01 +# 82 633 3198 2851 896 1581 1732 1292 295 132 +# Inhib_02 Inhib_03 Inhib_04 Inhib_05 Inhib_06 drop_01 drop_02 drop_03 drop_04 +# 461 2242 6018 1311 249 11 3 2 7 sce$cellType_broad_hc <- anno_hc2$broad[match(sce$collapsedCluster, anno_hc$cluster)] sce$cellType_hc <- anno_hc2$cellType[match(sce$collapsedCluster, anno_hc2$cluster)] table(sce$cellType_broad_hc) -# Astro EndoMural Micro Oligo OPC Excit Inhib Ambiguous +# Astro EndoMural Micro Oligo OPC Excit Inhib Ambiguous # 3979 2157 1601 10894 1940 24809 11067 21157 table(sce$cellType_hc) -# Astro EndoMural_01 EndoMural_02 Micro Oligo_01 Oligo_02 Oligo_03 OPC Excit_01 -# 3979 446 1711 1601 1868 4732 4294 1940 7927 -# Excit_02 Excit_03 Excit_04 Excit_05 Excit_06 Excit_07 Excit_08 Excit_09 Excit_10 -# 2487 1309 2171 2532 329 334 1463 2561 1079 -# Excit_11 Excit_12 Excit_13 Excit_14 Excit_15 Inhib_01 Inhib_02 Inhib_03 Inhib_04 -# 482 420 1567 82 66 5366 1267 1310 565 -# Inhib_05 Inhib_06 Ambiguous -# 1192 1367 21157 +# Astro EndoMural_01 EndoMural_02 Micro Oligo_01 Oligo_02 Oligo_03 OPC Excit_01 +# 3979 446 1711 1601 1868 4732 4294 1940 7927 +# Excit_02 Excit_03 Excit_04 Excit_05 Excit_06 Excit_07 Excit_08 Excit_09 Excit_10 +# 2487 1309 2171 2532 329 334 1463 2561 1079 +# Excit_11 Excit_12 Excit_13 Excit_14 Excit_15 Inhib_01 Inhib_02 Inhib_03 Inhib_04 +# 482 420 1567 82 66 5366 1267 1310 565 +# Inhib_05 Inhib_06 Ambiguous +# 1192 1367 21157 ## Check out prop (prop_k <- 100 * round(table(sce$cellType_broad_k) / ncol(sce), 3)) -# Astro EndoMural MicroOligo Oligo OPC Excit Inhib drop -# 4.6 1.7 7.1 43.4 2.3 27.4 13.4 0.0 +# Astro EndoMural MicroOligo Oligo OPC Excit Inhib drop +# 4.6 1.7 7.1 43.4 2.3 27.4 13.4 0.0 (prop_hc <- 100 * round(table(sce$cellType_broad_hc) / ncol(sce), 3)) -# Astro EndoMural Micro Oligo OPC Excit Inhib Ambiguous -# 5.1 2.8 2.1 14.0 2.5 32.0 14.3 27.3 +# Astro EndoMural Micro Oligo OPC Excit Inhib Ambiguous +# 5.1 2.8 2.1 14.0 2.5 32.0 14.3 27.3 as.data.frame(prop_k) %>% full_join(as.data.frame(prop_hc) %>% rename(Freq_HC = Freq)) # Var1 Freq Freq_HC diff --git a/code/03_build_sce/09_cluster_annotation_explore.R b/code/03_build_sce/09_cluster_annotation_explore.R index c64ba75..c73adc6 100644 --- a/code/03_build_sce/09_cluster_annotation_explore.R +++ b/code/03_build_sce/09_cluster_annotation_explore.R @@ -12,7 +12,7 @@ library("patchwork") # Plotting set up my_theme <- theme_bw() + - theme(text = element_text(size = 15)) + theme(text = element_text(size = 15)) plot_dir <- here("plots", "03_build_sce", "09_cluster_annotation_explore") if (!dir.exists(plot_dir)) dir.create(plot_dir) @@ -31,46 +31,46 @@ cell_type_colors_broad <- metadata(sce)$cell_type_colors_broad # scale_color_manual(values = cell_type_colors[levels(sce$cellType_hc)], drop = TRUE) + # my_theme + # coord_equal() -# +# # ggsave( # TSNE_HC_cellTypes + # guides(colour = guide_legend(override.aes = list(size = 2, alpha = 1))), # filename = here(plot_dir, "TSNE_HC-29_cellType.png"), width = 10 # ) -# -# +# +# # TSNE_HC_cellTypes <- ggcells(sce, mapping = aes(x = TSNE.1, y = TSNE.2)) + # geom_point(aes(colour = cellType_broad_hc), size = 0.2, alpha = 0.3) + # scale_color_manual(values = cell_type_colors_broad[levels(sce$cellType_broad_hc)], drop = TRUE) + # my_theme + # coord_equal() -# +# # ggsave( # TSNE_HC_cellTypes + # guides(colour = guide_legend(override.aes = list(size = 2, alpha = 1))), # filename = here(plot_dir, "TSNE_HC_broad_cellType.png"), width = 10 # ) -# -# +# +# # ## UMAP # UMAP_HC_cellTypes <- ggcells(sce, mapping = aes(x = UMAP.1, y = UMAP.2, colour = cellType_hc)) + # geom_point(size = 0.2, alpha = 0.3) + # scale_color_manual(values = cell_type_colors[levels(sce$cellType_hc)], drop = TRUE) + # my_theme + # coord_equal() -# +# # ggsave( # UMAP_HC_cellTypes + # guides(colour = guide_legend(override.aes = list(size = 2, alpha = 1))), # filename = here(plot_dir, "UMAP_HC_cellType.png"), width = 10 # ) -# +# # UMAP_HC_cellTypes <- ggcells(sce, mapping = aes(x = UMAP.1, y = UMAP.2)) + # geom_point(aes(colour = cellType_broad_hc), size = 0.2, alpha = 0.3) + # scale_color_manual(values = cell_type_colors_broad[levels(sce$cellType_broad_hc)], drop = TRUE) + # my_theme + # coord_equal() -# +# # ggsave( # UMAP_HC_cellTypes + # guides(colour = guide_legend(override.aes = list(size = 2, alpha = 1))), @@ -80,15 +80,15 @@ cell_type_colors_broad <- metadata(sce)$cell_type_colors_broad #### Plot KM TSNE + UMAP #### ## jsut plot mbkm here TSNE_km_cellTypes <- ggcells(sce, mapping = aes(x = TSNE.1, y = TSNE.2, colour = cellType_k)) + - geom_point(size = 0.2, alpha = 0.3) + - scale_color_manual(values = cell_type_colors[levels(sce$cellType_k)], drop = TRUE) + - my_theme + - coord_equal() + geom_point(size = 0.2, alpha = 0.3) + + scale_color_manual(values = cell_type_colors[levels(sce$cellType_k)], drop = TRUE) + + my_theme + + coord_equal() ggsave( - TSNE_km_cellTypes + - guides(colour = guide_legend(override.aes = list(size = 2, alpha = 1))), - filename = here(plot_dir, "TSNE_mbkm-29_cellType.png"), width = 10 + TSNE_km_cellTypes + + guides(colour = guide_legend(override.aes = list(size = 2, alpha = 1))), + filename = here(plot_dir, "TSNE_mbkm-29_cellType.png"), width = 10 ) @@ -98,19 +98,19 @@ source("my_plotMarkers.R") load(here("processed-data", "03_build_sce", "markers.mathys.tran.Rdata"), verbose = TRUE) my_plotMarkers( - sce = sce, - marker_list = markers.mathys.tran, - cat = "cellType_hc", - fill_colors = cell_type_colors, - pdf_fn = here(plot_dir, "markers_mathys_hc_29_ct.pdf") + sce = sce, + marker_list = markers.mathys.tran, + cat = "cellType_hc", + fill_colors = cell_type_colors, + pdf_fn = here(plot_dir, "markers_mathys_hc_29_ct.pdf") ) my_plotMarkers( - sce = sce, - marker_list = markers.mathys.tran, - cat = "cellType_k", - fill_colors = cell_type_colors, - pdf_fn = here(plot_dir, "markers_mathys_mb_kmeans_29.pdf") + sce = sce, + marker_list = markers.mathys.tran, + cat = "cellType_k", + fill_colors = cell_type_colors, + pdf_fn = here(plot_dir, "markers_mathys_mb_kmeans_29.pdf") ) @@ -122,22 +122,22 @@ cluster_compare <- table(sce$kmeans, sce$collapsedCluster) cluster_compare_prop <- sweep(cluster_compare, 2, colSums(cluster_compare), `/`) hc_anno <- colData(sce) |> - as.data.frame() |> - count(collapsedCluster, broad = cellType_broad_hc) |> - tibble::column_to_rownames("collapsedCluster") + as.data.frame() |> + count(collapsedCluster, broad = cellType_broad_hc) |> + tibble::column_to_rownames("collapsedCluster") km_anno <- colData(sce) |> - as.data.frame() |> - count(kmeans, broad = cellType_broad_k) |> - tibble::column_to_rownames("kmeans") + as.data.frame() |> + count(kmeans, broad = cellType_broad_k) |> + tibble::column_to_rownames("kmeans") broad_pallet <- cell_type_colors_broad[unique(c(hc_anno$broad, km_anno$broad))] png(here(plot_dir, "heatmap_cluster_compare.png"), height = 800, width = 800) pheatmap(cluster_compare_prop, - annotation_col = hc_anno, - annotation_row = km_anno, - annotation_colors = list(broad = broad_pallet) + annotation_col = hc_anno, + annotation_row = km_anno, + annotation_colors = list(broad = broad_pallet) ) dev.off() @@ -146,18 +146,18 @@ jacc.mat <- linkClustersMatrix(sce$kmeans, sce$collapsedCluster) png(here(plot_dir, "heatmap_cluster_compare_jacc.png"), height = 800, width = 800) pheatmap(jacc.mat, - annotation_col = hc_anno, - annotation_row = km_anno, - annotation_colors = list(broad = broad_pallet) + annotation_col = hc_anno, + annotation_row = km_anno, + annotation_colors = list(broad = broad_pallet) ) dev.off() ## best cluster pairing best <- max.col(jacc.mat, ties.method = "first") DataFrame( - Cluster = rownames(jacc.mat), - Corresponding = colnames(jacc.mat)[best], - Index = jacc.mat[cbind(seq_len(nrow(jacc.mat)), best)] + Cluster = rownames(jacc.mat), + Corresponding = colnames(jacc.mat)[best], + Index = jacc.mat[cbind(seq_len(nrow(jacc.mat)), best)] ) ## compare cell types @@ -165,20 +165,20 @@ ct_tab <- table(sce$cellType_k, sce$cellType_hc) ct_tab_prop <- sweep(ct_tab, 2, colSums(ct_tab), `/`) hc_ct_anno <- table(sce$cellType_hc) |> - as.data.frame() |> - rename(ct = Var1, n = Freq) + as.data.frame() |> + rename(ct = Var1, n = Freq) rownames(hc_ct_anno) <- hc_ct_anno$ct km_ct_anno <- table(sce$cellType_k) |> - as.data.frame() |> - rename(ct = Var1, n = Freq) + as.data.frame() |> + rename(ct = Var1, n = Freq) rownames(km_ct_anno) <- km_ct_anno$ct png(here(plot_dir, "heatmap_cellType_compare.png"), height = 800, width = 800) pheatmap(ct_tab_prop, - annotation_col = hc_ct_anno, - annotation_row = km_ct_anno, - annotation_colors = list(ct = cell_type_colors) + annotation_col = hc_ct_anno, + annotation_row = km_ct_anno, + annotation_colors = list(ct = cell_type_colors) ) dev.off() @@ -186,9 +186,9 @@ jacc.mat.ct <- linkClustersMatrix(sce$cellType_k, sce$cellType_hc) png(here(plot_dir, "heatmap_cellType_compare_jacc.png"), height = 800, width = 800) pheatmap(jacc.mat.ct, - annotation_col = hc_ct_anno, - annotation_row = km_ct_anno, - annotation_colors = list(ct = cell_type_colors) + annotation_col = hc_ct_anno, + annotation_row = km_ct_anno, + annotation_colors = list(ct = cell_type_colors) ) dev.off() @@ -198,20 +198,20 @@ dev.off() ctb_tab_prop <- sweep(ctb_tab, 2, colSums(ctb_tab), `/`) hc_broad_anno <- table(sce$cellType_broad_hc) |> - as.data.frame() |> - rename(broad = Var1, n = Freq) + as.data.frame() |> + rename(broad = Var1, n = Freq) rownames(hc_broad_anno) <- hc_broad_anno$broad km_broad_anno <- table(sce$cellType_broad_k) |> - as.data.frame() |> - rename(broad = Var1, n = Freq) + as.data.frame() |> + rename(broad = Var1, n = Freq) rownames(km_broad_anno) <- km_broad_anno$broad png(here(plot_dir, "heatmap_cellType_broad_compare.png"), height = 800, width = 800) pheatmap(ctb_tab_prop, - annotation_col = hc_broad_anno, - annotation_row = km_broad_anno, - annotation_colors = list(broad = broad_pallet) + annotation_col = hc_broad_anno, + annotation_row = km_broad_anno, + annotation_colors = list(broad = broad_pallet) ) dev.off() @@ -219,9 +219,9 @@ jacc.mat.broad <- linkClustersMatrix(sce$cellType_broad_k, sce$cellType_broad_hc png(here(plot_dir, "heatmap_cellType_broad_compare_jacc.png"), height = 800, width = 800) pheatmap(jacc.mat.broad, - annotation_col = hc_broad_anno, - annotation_row = km_broad_anno, - annotation_colors = list(broad = broad_pallet) + annotation_col = hc_broad_anno, + annotation_row = km_broad_anno, + annotation_colors = list(broad = broad_pallet) ) dev.off() @@ -234,29 +234,29 @@ pairwiseRand(sce$cellType_broad_k, sce$cellType_broad_hc, mode = "index") #### Explore UMI #### UMI_ct_k <- ggcells(sce, mapping = aes(x = cellType_k, y = sum, fill = cellType_k)) + - geom_boxplot() + - scale_fill_manual(values = cell_type_colors[levels(sce$cellType_k)], drop = TRUE) + - my_theme + - theme( - legend.position = "None", axis.title.x = element_blank(), - axis.text.x = element_text(angle = 90, hjust = 1) - ) + geom_boxplot() + + scale_fill_manual(values = cell_type_colors[levels(sce$cellType_k)], drop = TRUE) + + my_theme + + theme( + legend.position = "None", axis.title.x = element_blank(), + axis.text.x = element_text(angle = 90, hjust = 1) + ) ggsave(UMI_ct_k, filename = here(plot_dir, "qc_UMI_mbkm-29_cellType.png"), width = 10) UMI_ct_hc <- ggcells(sce, mapping = aes( - x = reorder(cellType_hc, sum, FUN = median), - y = sum, fill = cellType_hc + x = reorder(cellType_hc, sum, FUN = median), + y = sum, fill = cellType_hc )) + - geom_boxplot() + - scale_fill_manual(values = cell_type_colors) + - my_theme + - theme( - legend.position = "None", axis.title.x = element_blank(), - axis.text.x = element_text(angle = 45, hjust = 1) - ) + - labs(y = "Total UMIs") + geom_boxplot() + + scale_fill_manual(values = cell_type_colors) + + my_theme + + theme( + legend.position = "None", axis.title.x = element_blank(), + axis.text.x = element_text(angle = 45, hjust = 1) + ) + + labs(y = "Total UMIs") ggsave(UMI_ct_hc, filename = here(plot_dir, "qc_UMI_HC_cellType.png"), height = 3) ggsave(UMI_ct_hc, filename = here(plot_dir, "qc_UMI_HC_cellType.pdf"), height = 3) @@ -264,49 +264,49 @@ ggsave(UMI_ct_hc, filename = here(plot_dir, "qc_UMI_HC_cellType.pdf"), height = #### Explore doublet scores dbs_ct_k <- ggcells(sce, mapping = aes(x = cellType_k, y = doubletScore, fill = cellType_k)) + - geom_boxplot() + - scale_fill_manual(values = cell_type_colors) + - my_theme + - theme( - legend.position = "None", axis.title.x = element_blank(), - axis.text.x = element_text(angle = 90, hjust = 1) - ) + geom_boxplot() + + scale_fill_manual(values = cell_type_colors) + + my_theme + + theme( + legend.position = "None", axis.title.x = element_blank(), + axis.text.x = element_text(angle = 90, hjust = 1) + ) ggsave(dbs_ct_k, filename = here(plot_dir, "qc_doubletScore_mbkm-29_cellType.png"), width = 10) dbs_ct_hc <- ggcells(sce, - mapping = aes( - x = reorder(cellType_hc, doubletScore, FUN = median), - y = doubletScore, - fill = cellType_hc - ) + mapping = aes( + x = reorder(cellType_hc, doubletScore, FUN = median), + y = doubletScore, + fill = cellType_hc + ) ) + - geom_boxplot() + - scale_fill_manual(values = cell_type_colors) + - my_theme + - theme( - legend.position = "None", axis.title.x = element_blank(), - axis.text.x = element_text(angle = 45, hjust = 1) - ) + - labs(y = "Doublet Score") + geom_boxplot() + + scale_fill_manual(values = cell_type_colors) + + my_theme + + theme( + legend.position = "None", axis.title.x = element_blank(), + axis.text.x = element_text(angle = 45, hjust = 1) + ) + + labs(y = "Doublet Score") ggsave(dbs_ct_hc, filename = here(plot_dir, "qc_doubletScore_HC_cellType.png"), height = 3) ggsave(dbs_ct_hc, filename = here(plot_dir, "qc_doubletScore_HC_cellType.pdf"), height = 3) #### Explore Mito rate #### mito_ct_hc <- ggcells(sce, mapping = aes( - x = reorder(cellType_hc, subsets_Mito_percent, FUN = median), - y = subsets_Mito_percent, fill = cellType_hc + x = reorder(cellType_hc, subsets_Mito_percent, FUN = median), + y = subsets_Mito_percent, fill = cellType_hc )) + - geom_boxplot() + - scale_fill_manual(values = cell_type_colors) + - my_theme + - theme( - legend.position = "None", axis.title.x = element_blank(), - axis.text.x = element_text(angle = 45, hjust = 1) - ) + - labs(y = "Percent Mito") + geom_boxplot() + + scale_fill_manual(values = cell_type_colors) + + my_theme + + theme( + legend.position = "None", axis.title.x = element_blank(), + axis.text.x = element_text(angle = 45, hjust = 1) + ) + + labs(y = "Percent Mito") ggsave(mito_ct_hc, filename = here(plot_dir, "qc_mitoRate_HC_cellType.png"), height = 3) ggsave(mito_ct_hc, filename = here(plot_dir, "qc_mitoRate_HC_cellType.pdf"), height = 3) @@ -319,9 +319,9 @@ cellType.idx <- splitit(sce$cellType_hc) # sapply(c("Excit", "Inhib", "MSN"), function(x){grep(x, names(cellType.idx))}) sapply(cellType.idx, function(x) { - quantile(sce$doubletScore[x]) + quantile(sce$doubletScore[x]) })[, order(sapply(cellType.idx, function(x) { - quantile(sce$doubletScore[x])["50%"] + quantile(sce$doubletScore[x])["50%"] }))] # Excit_15 Micro Oligo_03 Excit_09 Astro OPC Excit_10 Oligo_02 Excit_13 Excit_11 EndoMural_02 # 0% 0.0086180 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000 0.0000000 0.000000 0.050708 0.000000 @@ -343,9 +343,9 @@ sapply(cellType.idx, function(x) { # 100% 12.068628 13.872384 9.825768 8.841504 6.457070 8.640852 5.808176 4.430430 sapply(cellType.idx, function(x) { - quantile(sce$sum[x]) + quantile(sce$sum[x]) })[, order(sapply(cellType.idx, function(x) { - quantile(sce$sum[x])["50%"] + quantile(sce$sum[x])["50%"] }))] # Excit_15 drop Excit_09 Oligo_03 Oligo_01 Micro EndoMural_02 Astro Excit_13 Oligo_02 Excit_10 Excit_05 OPC # 0% 678 220 462 359 300.00 350 407.0 276.0 648.0 968.0 533 292.00 397.00 @@ -367,9 +367,9 @@ sapply(cellType.idx, function(x) { # 100% 179194.0 203794.0 188090.0 152333 296099.0 sapply(cellType.idx, function(x) { - quantile(sce$subsets_Mito_percent[x]) + quantile(sce$subsets_Mito_percent[x]) })[, order(sapply(cellType.idx, function(x) { - quantile(sce$subsets_Mito_percent[x])["50%"] + quantile(sce$subsets_Mito_percent[x])["50%"] }))] # Excit_06 Excit_08 Excit_02 Excit_07 Excit_01 Inhib_04 Excit_03 Inhib_05 Excit_04 Inhib_06 # 0% 0.002438608 0.001886792 0.002068851 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.003583716 0.0000000 @@ -416,22 +416,22 @@ table(unfactor(sce$prelimCluster)[sce$cellType_hc == "Excit_13"]) # 446 249 58 162 28 462 121 41 my_plotMarkers( - sce = sce[, sce$cellType_hc == "Excit_13"], - marker_list = markers.mathys.tran[c("astrocyte", "excit_neuron")], - cat = "prelimCluster", - pdf_fn = here(plot_dir, "markers_Excit_13_check.pdf") + sce = sce[, sce$cellType_hc == "Excit_13"], + marker_list = markers.mathys.tran[c("astrocyte", "excit_neuron")], + cat = "prelimCluster", + pdf_fn = here(plot_dir, "markers_Excit_13_check.pdf") ) dbs_excit13 <- ggcells(sce[, sce$cellType_hc == "Excit_13"], - mapping = aes(x = prelimCluster, y = doubletScore, fill = prelimCluster) + mapping = aes(x = prelimCluster, y = doubletScore, fill = prelimCluster) ) + - geom_boxplot() + - geom_hline(yintercept = 5, color = "red", linetype = "dashed") + - theme_bw() + - theme( - legend.position = "None", axis.title.x = element_blank(), - axis.text.x = element_text(angle = 90, hjust = 1) - ) + geom_boxplot() + + geom_hline(yintercept = 5, color = "red", linetype = "dashed") + + theme_bw() + + theme( + legend.position = "None", axis.title.x = element_blank(), + axis.text.x = element_text(angle = 90, hjust = 1) + ) ggsave(dbs_excit13, filename = here(plot_dir, "qc_Excit_13_doublet_scores.png")) @@ -479,7 +479,7 @@ ggsave(dbs_excit13, filename = here(plot_dir, "qc_Excit_13_doublet_scores.png")) # # ggsave(prelim_doublets, filename = here(plot_dir, "prelim_doublets.png"), width = 30) -sgejobs::job_single('09_cluster_annotation_explore', create_shell = TRUE, queue= 'bluejay', memory = '5G', command = "Rscript 09_cluster_annotation_explore.R") +sgejobs::job_single("09_cluster_annotation_explore", create_shell = TRUE, queue = "bluejay", memory = "5G", command = "Rscript 09_cluster_annotation_explore.R") ## Reproducibility information print("Reproducibility information:") diff --git a/code/03_build_sce/09_find_markers.R b/code/03_build_sce/09_find_markers.R index 1de514d..c588f9f 100644 --- a/code/03_build_sce/09_find_markers.R +++ b/code/03_build_sce/09_find_markers.R @@ -1,4 +1,3 @@ - library("SingleCellExperiment") library("scater") library("here") diff --git a/code/03_build_sce/09_find_markers_layer.R b/code/03_build_sce/09_find_markers_layer.R index 00ebe0b..2d325e5 100644 --- a/code/03_build_sce/09_find_markers_layer.R +++ b/code/03_build_sce/09_find_markers_layer.R @@ -1,4 +1,3 @@ - library("SingleCellExperiment") library("scater") library("here") diff --git a/code/03_build_sce/11_advanced_marker_search.R b/code/03_build_sce/11_advanced_marker_search.R index 604e85f..9278777 100644 --- a/code/03_build_sce/11_advanced_marker_search.R +++ b/code/03_build_sce/11_advanced_marker_search.R @@ -1,4 +1,3 @@ - library("SingleCellExperiment") library("scater") library("here") diff --git a/code/03_build_sce/cell_colors.R b/code/03_build_sce/cell_colors.R index e3a78aa..98d134a 100644 --- a/code/03_build_sce/cell_colors.R +++ b/code/03_build_sce/cell_colors.R @@ -113,13 +113,13 @@ expand_cell_colors <- function(cell_colors, cell_types, split = "_") { cell_types <- readLines(here("processed-data", "03_build_sce", "cell_types.txt")) cell_type_colors <- expand_cell_colors(cell_type_colors_broad, cell_types) (cell_type_colors <- cell_type_colors[cell_types]) -# Astro EndoMural EndoMural_01 EndoMural_02 Micro MicroOligo Oligo_01 Oligo_02 Oligo_03 -# "#3BB273" "#FF56AF" "#FF56AF" "#FFAAD7" "#663894" "#AB0091" "#E07000" "#EA9F55" "#F4CFAA" -# OPC Excit_01 Excit_02 Excit_03 Excit_04 Excit_05 Excit_06 Excit_07 Excit_08 -# "#D2B037" "#247FBC" "#3287C0" "#4190C4" "#4F98C9" "#5EA1CD" "#6DA9D2" "#7BB2D6" "#8ABADB" -# Excit_09 Excit_10 Excit_11 Excit_12 Excit_13 Excit_14 Excit_15 Inhib_01 Inhib_02 -# "#98C3DF" "#A7CBE4" "#B5D4E8" "#C4DCED" "#D3E5F1" "#E1EDF6" "#F0F6FA" "#E94F37" "#EC6C58" -# Inhib_03 Inhib_04 Inhib_05 Inhib_06 Ambiguous drop_01 drop_02 drop_03 drop_04 +# Astro EndoMural EndoMural_01 EndoMural_02 Micro MicroOligo Oligo_01 Oligo_02 Oligo_03 +# "#3BB273" "#FF56AF" "#FF56AF" "#FFAAD7" "#663894" "#AB0091" "#E07000" "#EA9F55" "#F4CFAA" +# OPC Excit_01 Excit_02 Excit_03 Excit_04 Excit_05 Excit_06 Excit_07 Excit_08 +# "#D2B037" "#247FBC" "#3287C0" "#4190C4" "#4F98C9" "#5EA1CD" "#6DA9D2" "#7BB2D6" "#8ABADB" +# Excit_09 Excit_10 Excit_11 Excit_12 Excit_13 Excit_14 Excit_15 Inhib_01 Inhib_02 +# "#98C3DF" "#A7CBE4" "#B5D4E8" "#C4DCED" "#D3E5F1" "#E1EDF6" "#F0F6FA" "#E94F37" "#EC6C58" +# Inhib_03 Inhib_04 Inhib_05 Inhib_06 Ambiguous drop_01 drop_02 drop_03 drop_04 # "#F08979" "#F3A79B" "#F7C4BC" "#FBE1DD" "#A9998F" "#000000" "#3F3F3F" "#7F7F7F" "#BFBFBF" ## plot previews of pallets @@ -171,9 +171,9 @@ dev.off() cell_type_colors_layer <- c(cell_type_colors_layer, Excit_ambig = "#A0A7A7") cell_type_colors_layer <- cell_type_colors_layer[cell_types_layer_all] -# Astro EndoMural Micro Oligo OPC Excit_L2/3 Excit_L3 Excit_L3/4/5 Excit_L4 -# "#3BB273" "#FF56AF" "#663894" "#E07000" "#D2B037" "#D0D1E6" "#A6BDDB" "#74A9CF" "#3690C0" -# Excit_L5 Excit_L5/6 Excit_L6 Excit_ambig Inhib +# Astro EndoMural Micro Oligo OPC Excit_L2/3 Excit_L3 Excit_L3/4/5 Excit_L4 +# "#3BB273" "#FF56AF" "#663894" "#E07000" "#D2B037" "#D0D1E6" "#A6BDDB" "#74A9CF" "#3690C0" +# Excit_L5 Excit_L5/6 Excit_L6 Excit_ambig Inhib # "#0570B0" "#045A8D" "#023858" "#A0A7A7" "#E94F37" ## Save diff --git a/code/03_build_sce/custom_plotExpression.R b/code/03_build_sce/custom_plotExpression.R index 92551e7..f889ef0 100644 --- a/code/03_build_sce/custom_plotExpression.R +++ b/code/03_build_sce/custom_plotExpression.R @@ -1,4 +1,3 @@ - #' Plot SCE Expression of marker genes over clusters #' #' @param sce A SingleCellExperiment object containing expression values diff --git a/code/03_build_sce/my_plotMarkers.R b/code/03_build_sce/my_plotMarkers.R index 2f0f6ac..e40d26e 100644 --- a/code/03_build_sce/my_plotMarkers.R +++ b/code/03_build_sce/my_plotMarkers.R @@ -1,38 +1,37 @@ - #' Plot list of marker genes #' -#' @param sce -#' @param marker_list -#' @param assay -#' @param cat -#' @param fill_colors -#' @param pdf_fn +#' @param sce +#' @param marker_list +#' @param assay +#' @param cat +#' @param fill_colors +#' @param pdf_fn #' #' @return #' @export #' #' @examples my_plotMarkers <- function(sce, marker_list, assay = "logcounts", cat = "cellType", fill_colors = NULL, pdf_fn) { - message("plotting: ", pdf_fn) - pdf(pdf_fn, height = 6, width = 8) - - for (i in 1:length(marker_list)) { - message(names(marker_list)[[i]]) - - m <- marker_list[[i]] - markers <- m[m %in% rownames(sce)] - - # if(m != markers) message("Missing...",paste(m[!m %in% markers], collapse = ", ")) - if(!identical(m, markers)) message("Missing markers...") - print( - custom_plotExpression(sce, - genes = markers, - title = names(marker_list)[[i]], - cat = cat, - fill_colors = fill_colors - ) - ) - } - - dev.off() + message("plotting: ", pdf_fn) + pdf(pdf_fn, height = 6, width = 8) + + for (i in 1:length(marker_list)) { + message(names(marker_list)[[i]]) + + m <- marker_list[[i]] + markers <- m[m %in% rownames(sce)] + + # if(m != markers) message("Missing...",paste(m[!m %in% markers], collapse = ", ")) + if (!identical(m, markers)) message("Missing markers...") + print( + custom_plotExpression(sce, + genes = markers, + title = names(marker_list)[[i]], + cat = cat, + fill_colors = fill_colors + ) + ) + } + + dev.off() } diff --git a/code/03_build_sce/prep_markers.R b/code/03_build_sce/prep_markers.R index 404262c..3d919dd 100644 --- a/code/03_build_sce/prep_markers.R +++ b/code/03_build_sce/prep_markers.R @@ -1,4 +1,3 @@ - markers.mathys.tran <- list( "neuron" = c("SYT1", "SNAP25", "GRIN1", "CAMK2A", "NRGN"), "excit_neuron" = c("SLC17A7", "SLC17A6", "SLC17A8"), # a7 coritical a6 sub-subcortical diff --git a/code/03_build_sce/utils.R b/code/03_build_sce/utils.R index 9206c5e..db4bf33 100644 --- a/code/03_build_sce/utils.R +++ b/code/03_build_sce/utils.R @@ -1,4 +1,3 @@ - library("patchwork") #' Plot Reduced Dim of SCE #' wrapper function for ggcells adding facet wrap to selected category @@ -46,9 +45,10 @@ plot_reducedDim_qc <- function(sce, type = "TSNE", color_by = "sum", title = "") } -plotExpressionCustom <- function(sce, features, features_name, anno_name = "cellType", - point_alpha = 0.2, point_size = 0.7, ncol = 2, xlab = NULL, - exprs_values = "logcounts", scales = "fixed") { +plotExpressionCustom <- function( + sce, features, features_name, anno_name = "cellType", + point_alpha = 0.2, point_size = 0.7, ncol = 2, xlab = NULL, + exprs_values = "logcounts", scales = "fixed") { scater::plotExpression(sce, exprs_values = exprs_values, features = features, diff --git a/code/05_explore_sce/00_pub_metrics.R b/code/05_explore_sce/00_pub_metrics.R index ec04231..30c7f97 100644 --- a/code/05_explore_sce/00_pub_metrics.R +++ b/code/05_explore_sce/00_pub_metrics.R @@ -1,4 +1,3 @@ - library("SingleCellExperiment") library("here") library("tidyverse") @@ -55,97 +54,99 @@ summary(zero_sums) #### Supp Table export #### -list.files(here("processed-data","03_build_sce")) -list.files(here("processed-data","02_cellranger_metrics")) +list.files(here("processed-data", "03_build_sce")) +list.files(here("processed-data", "02_cellranger_metrics")) -sample_info <- read.csv(here("processed-data","03_build_sce","sample_info.csv")) |> - rename(Donor = subject, Position = region) +sample_info <- read.csv(here("processed-data", "03_build_sce", "sample_info.csv")) |> + rename(Donor = subject, Position = region) -droplet_summary <- read.csv(here("processed-data","03_build_sce","drop_summary.csv")) |> - select(Sample, total_droplets = total_n, droplet_lower_cutoff = lower_cutoff, pre_QC_n = non_empty) +droplet_summary <- read.csv(here("processed-data", "03_build_sce", "drop_summary.csv")) |> + select(Sample, total_droplets = total_n, droplet_lower_cutoff = lower_cutoff, pre_QC_n = non_empty) -cell_ranger_info <-read.csv(here("processed-data","02_cellranger_metrics","cellranger_metrics.csv")) |> - select(file_id = X, Number.of.Reads, Mean.Reads.per.Cell, Median.UMI.Counts.per.Cell, Median.Genes.per.Cell) #mean unique molecular indices - only median in this table +cell_ranger_info <- read.csv(here("processed-data", "02_cellranger_metrics", "cellranger_metrics.csv")) |> + select(file_id = X, Number.of.Reads, Mean.Reads.per.Cell, Median.UMI.Counts.per.Cell, Median.Genes.per.Cell) # mean unique molecular indices - only median in this table colnames(sample_info_all) -# [1] "Sample" "file_id" -# [3] "region" "subject" -# [5] "round" "n" -# [7] "n_high_mito" "n_low_sum" -# [9] "n_low_detect" "n_discard_auto" -# [11] "Estimated.Number.of.Cells" "Mean.Reads.per.Cell" -# [13] "Median.Genes.per.Cell" "Number.of.Reads" -# [15] "Valid.Barcodes" "Sequencing.Saturation" -# [17] "Q30.Bases.in.Barcode" "Q30.Bases.in.RNA.Read" -# [19] "Q30.Bases.in.UMI" "Reads.Mapped.to.Genome" +# [1] "Sample" "file_id" +# [3] "region" "subject" +# [5] "round" "n" +# [7] "n_high_mito" "n_low_sum" +# [9] "n_low_detect" "n_discard_auto" +# [11] "Estimated.Number.of.Cells" "Mean.Reads.per.Cell" +# [13] "Median.Genes.per.Cell" "Number.of.Reads" +# [15] "Valid.Barcodes" "Sequencing.Saturation" +# [17] "Q30.Bases.in.Barcode" "Q30.Bases.in.RNA.Read" +# [19] "Q30.Bases.in.UMI" "Reads.Mapped.to.Genome" # [21] "Reads.Mapped.Confidently.to.Genome" "Reads.Mapped.Confidently.to.Intergenic.Regions" -# [23] "Reads.Mapped.Confidently.to.Intronic.Regions" "Reads.Mapped.Confidently.to.Exonic.Regions" -# [25] "Reads.Mapped.Confidently.to.Transcriptome" "Reads.Mapped.Antisense.to.Gene" -# [27] "Fraction.Reads.in.Cells" "Total.Genes.Detected" +# [23] "Reads.Mapped.Confidently.to.Intronic.Regions" "Reads.Mapped.Confidently.to.Exonic.Regions" +# [25] "Reads.Mapped.Confidently.to.Transcriptome" "Reads.Mapped.Antisense.to.Gene" +# [27] "Fraction.Reads.in.Cells" "Total.Genes.Detected" # [29] "Median.UMI.Counts.per.Cell" ## from manuscript # droplets were sequenced -sample_info_all <- sample_info |> - left_join(droplet_summary) |> - left_join(cell_ranger_info) |> - select(Sample, file_id, Donor, Position, round, ## sample info - Number.of.Reads, Mean.Reads.per.Cell, Median.UMI.Counts.per.Cell, Median.Genes.per.Cell, ## cell ranger - total_droplets, droplet_lower_cutoff, pre_QC_n, ## empty drops - n_high_mito, n_low_sum, n_low_detect, n_discard_auto) |> ## QC - mutate(n = pre_QC_n - n_discard_auto) - -sample_info_all |> - select(Median.Genes.per.Cell)|> - summary() - -# Sample file_id Donor Position round Number.of.Reads -# Length:19 Length:19 Length:19 Length:19 Length:19 Min. :162821623 -# Class :character Class :character Class :character Class :character Class :character 1st Qu.:275981256 -# Mode :character Mode :character Mode :character Mode :character Mode :character Median :323114123 -# Mean :351032590 -# 3rd Qu.:483252608 -# Max. :541289607 +sample_info_all <- sample_info |> + left_join(droplet_summary) |> + left_join(cell_ranger_info) |> + select( + Sample, file_id, Donor, Position, round, ## sample info + Number.of.Reads, Mean.Reads.per.Cell, Median.UMI.Counts.per.Cell, Median.Genes.per.Cell, ## cell ranger + total_droplets, droplet_lower_cutoff, pre_QC_n, ## empty drops + n_high_mito, n_low_sum, n_low_detect, n_discard_auto + ) |> ## QC + mutate(n = pre_QC_n - n_discard_auto) + +sample_info_all |> + select(Median.Genes.per.Cell) |> + summary() + +# Sample file_id Donor Position round Number.of.Reads +# Length:19 Length:19 Length:19 Length:19 Length:19 Min. :162821623 +# Class :character Class :character Class :character Class :character Class :character 1st Qu.:275981256 +# Mode :character Mode :character Mode :character Mode :character Mode :character Median :323114123 +# Mean :351032590 +# 3rd Qu.:483252608 +# Max. :541289607 # Mean.Reads.per.Cell Median.UMI.Counts.per.Cell Median.Genes.per.Cell total_droplets droplet_lower_cutoff -# Min. : 38062 Min. : 1323 Min. : 985 Min. : 604278 Min. :219.0 -# 1st Qu.: 47784 1st Qu.: 2144 1st Qu.:1434 1st Qu.: 944401 1st Qu.:276.5 -# Median : 55811 Median : 3334 Median :1959 Median :1495673 Median :390.0 -# Mean : 65253 Mean : 5253 Mean :2415 Mean :1400306 Mean :441.3 -# 3rd Qu.: 75312 3rd Qu.: 7366 3rd Qu.:3259 3rd Qu.:1821096 3rd Qu.:581.0 -# Max. :136799 Max. :14627 Max. :5146 Max. :2209670 Max. :910.0 -# pre_QC_n n_high_mito n_low_sum n_low_detect n_discard_auto n -# Min. :2989 Min. : 0.0 Min. : 0.00 Min. : 0.00 Min. : 0.0 Min. :2661 -# 1st Qu.:3682 1st Qu.: 26.5 1st Qu.: 0.00 1st Qu.: 0.00 1st Qu.: 31.5 1st Qu.:3156 -# Median :4134 Median : 258.0 Median : 0.00 Median : 0.00 Median : 328.0 Median :4067 -# Mean :4461 Mean : 345.2 Mean : 24.32 Mean : 65.42 Mean : 376.4 Mean :4084 -# 3rd Qu.:5270 3rd Qu.: 565.0 3rd Qu.: 0.00 3rd Qu.: 70.50 3rd Qu.: 635.0 3rd Qu.:4716 +# Min. : 38062 Min. : 1323 Min. : 985 Min. : 604278 Min. :219.0 +# 1st Qu.: 47784 1st Qu.: 2144 1st Qu.:1434 1st Qu.: 944401 1st Qu.:276.5 +# Median : 55811 Median : 3334 Median :1959 Median :1495673 Median :390.0 +# Mean : 65253 Mean : 5253 Mean :2415 Mean :1400306 Mean :441.3 +# 3rd Qu.: 75312 3rd Qu.: 7366 3rd Qu.:3259 3rd Qu.:1821096 3rd Qu.:581.0 +# Max. :136799 Max. :14627 Max. :5146 Max. :2209670 Max. :910.0 +# pre_QC_n n_high_mito n_low_sum n_low_detect n_discard_auto n +# Min. :2989 Min. : 0.0 Min. : 0.00 Min. : 0.00 Min. : 0.0 Min. :2661 +# 1st Qu.:3682 1st Qu.: 26.5 1st Qu.: 0.00 1st Qu.: 0.00 1st Qu.: 31.5 1st Qu.:3156 +# Median :4134 Median : 258.0 Median : 0.00 Median : 0.00 Median : 328.0 Median :4067 +# Mean :4461 Mean : 345.2 Mean : 24.32 Mean : 65.42 Mean : 376.4 Mean :4084 +# 3rd Qu.:5270 3rd Qu.: 565.0 3rd Qu.: 0.00 3rd Qu.: 70.50 3rd Qu.: 635.0 3rd Qu.:4716 # Max. :6269 Max. :1229.0 Max. :192.00 Max. :412.00 Max. :1244.0 Max. :5911 library(xlsx) key <- data.frame( - colnames = colnames(sample_info_all) + colnames = colnames(sample_info_all) ) ## Clear file and write key annotation_xlsx <- - here("processed-data","05_explore_sce", "00_pub_metrics", "sn_sample_info.xlsx") + here("processed-data", "05_explore_sce", "00_pub_metrics", "sn_sample_info.xlsx") write.xlsx( - key, - file = annotation_xlsx, - sheetName = "Key", - append = FALSE, - row.names = FALSE + key, + file = annotation_xlsx, + sheetName = "Key", + append = FALSE, + row.names = FALSE ) ## write annotations write.xlsx( - sample_info_all, - file = annotation_xlsx, - sheetName = paste0("snRNA-seq_sample_info"), - append = TRUE, - row.names = FALSE + sample_info_all, + file = annotation_xlsx, + sheetName = paste0("snRNA-seq_sample_info"), + append = TRUE, + row.names = FALSE ) diff --git a/code/05_explore_sce/01_reduced_dim_plots.R b/code/05_explore_sce/01_reduced_dim_plots.R index 7f073f0..1824f74 100644 --- a/code/05_explore_sce/01_reduced_dim_plots.R +++ b/code/05_explore_sce/01_reduced_dim_plots.R @@ -28,44 +28,48 @@ levels(sce$cellType_hc) cell_type_colors <- metadata(sce)$cell_type_colors[levels(sce$cellType_hc)] TSNE_cellTypes_hc <- ggcells(sce, mapping = aes(x = TSNE.1, y = TSNE.2, colour = cellType_hc)) + - geom_point(size = 0.2, alpha = 0.3) + - scale_color_manual(values = cell_type_colors) + - my_theme + - coord_equal() + - labs(x = "TSNE Dimension 1", y = "TSNE Dimension 2") + geom_point(size = 0.2, alpha = 0.3) + + scale_color_manual(values = cell_type_colors) + + my_theme + + coord_equal() + + labs(x = "TSNE Dimension 1", y = "TSNE Dimension 2") -ggsave(TSNE_cellTypes_hc + - guides(colour = guide_legend(override.aes = list(size = 2, alpha = 1))), - filename = here(plot_dir, "TSNE_cellType-Ambig.png"), - width = 9, height = 6) +ggsave( + TSNE_cellTypes_hc + + guides(colour = guide_legend(override.aes = list(size = 2, alpha = 1))), + filename = here(plot_dir, "TSNE_cellType-Ambig.png"), + width = 9, height = 6 +) -ggsave(TSNE_cellTypes_hc + - guides(colour = guide_legend(override.aes = list(size = 2, alpha = 1))), - filename = here(plot_dir, "TSNE_cellType-Ambig.pdf"), - width = 9, height = 6) +ggsave( + TSNE_cellTypes_hc + + guides(colour = guide_legend(override.aes = list(size = 2, alpha = 1))), + filename = here(plot_dir, "TSNE_cellType-Ambig.pdf"), + width = 9, height = 6 +) ## No legends TSNE_cellTypes_hc_no_legend <- TSNE_cellTypes_hc + theme(legend.position = "None") TSNE_cellTypes_hc_facet <- TSNE_cellTypes_hc_no_legend + - facet_wrap(~cellType_hc) + facet_wrap(~cellType_hc) ## Add full + facet for cell types ggsave( - TSNE_cellTypes_hc_no_legend + - TSNE_cellTypes_hc_facet + - theme(axis.title.y = element_blank()), - filename = here(plot_dir, "TSNE_cellType_full_facet-Ambig.png"), - width = 13 + TSNE_cellTypes_hc_no_legend + + TSNE_cellTypes_hc_facet + + theme(axis.title.y = element_blank()), + filename = here(plot_dir, "TSNE_cellType_full_facet-Ambig.png"), + width = 13 ) ggsave( - TSNE_cellTypes_hc_no_legend + - TSNE_cellTypes_hc_facet + - theme(axis.title.y = element_blank()), - filename = here(plot_dir, "TSNE_cellType_full_facet-Ambig.pdf"), - width = 13 + TSNE_cellTypes_hc_no_legend + + TSNE_cellTypes_hc_facet + + theme(axis.title.y = element_blank()), + filename = here(plot_dir, "TSNE_cellType_full_facet-Ambig.pdf"), + width = 13 ) #### Exclude Ambiguous nuc #### @@ -218,16 +222,16 @@ ggsave( ) ## no legend -TSNE_cellType_layer_no_legend <- TSNE_cellType_layer + - theme(legend.position = "None") +TSNE_cellType_layer_no_legend <- TSNE_cellType_layer + + theme(legend.position = "None") ggsave( - TSNE_cellType_layer_no_legend + - TSNE_cellType_layer_no_legend+ - facet_wrap(~cellType_layer) + - theme(axis.title.y = element_blank()), - filename = here(plot_dir, "TSNE_cellType_layer_ambig_facet.png"), - width = 13 + TSNE_cellType_layer_no_legend + + TSNE_cellType_layer_no_legend + + facet_wrap(~cellType_layer) + + theme(axis.title.y = element_blank()), + filename = here(plot_dir, "TSNE_cellType_layer_ambig_facet.png"), + width = 13 ) ## Drop Ambig @@ -237,57 +241,64 @@ sce$cellType_hc <- droplevels(sce$cellType_hc) cell_type_colors_layer <- metadata(sce)$cell_type_colors_layer[levels(sce$cellType_layer)] TSNE_cellType_layer <- ggcells(sce, mapping = aes(x = TSNE.1, y = TSNE.2, colour = cellType_layer)) + - geom_point(size = 0.2, alpha = 0.3) + - scale_color_manual(values = cell_type_colors_layer) + - my_theme + - coord_equal() + - labs(x = "TSNE Dimension 1", y = "TSNE Dimension 2") + geom_point(size = 0.2, alpha = 0.3) + + scale_color_manual(values = cell_type_colors_layer) + + my_theme + + coord_equal() + + labs(x = "TSNE Dimension 1", y = "TSNE Dimension 2") ggsave( - TSNE_cellType_layer + - guides(colour = guide_legend(override.aes = list(size = 2, alpha = 1))), - filename = here(plot_dir, "TSNE_cellType_layer.png"), width = 9 + TSNE_cellType_layer + + guides(colour = guide_legend(override.aes = list(size = 2, alpha = 1))), + filename = here(plot_dir, "TSNE_cellType_layer.png"), width = 9 ) ## no legend -TSNE_cellType_layer_no_legend <- TSNE_cellType_layer + - theme(legend.position = "None") +TSNE_cellType_layer_no_legend <- TSNE_cellType_layer + + theme(legend.position = "None") ggsave( - TSNE_cellType_layer_no_legend, - filename = here(plot_dir, "TSNE_cellType_layer_no_legend.png"), width = 9 + TSNE_cellType_layer_no_legend, + filename = here(plot_dir, "TSNE_cellType_layer_no_legend.png"), width = 9 ) ## Add lables -tsne_postions <- tibble::tibble(TSNE_1 = reducedDims(sce)$TSNE[,1], - TSNE_2 = reducedDims(sce)$TSNE[,2], - cellType_layer = sce$cellType_layer) |> - # dplyr::filter((TSNE_2 > 0 & grepl("Excit", cellType_layer)) | !grepl("Excit", cellType_layer)) |> - dplyr::group_by(cellType_layer) |> - dplyr::summarise(TSNE.1 = median(TSNE_1), - TSNE.2 = median(TSNE_2), - ) -##geom label +tsne_postions <- tibble::tibble( + TSNE_1 = reducedDims(sce)$TSNE[, 1], + TSNE_2 = reducedDims(sce)$TSNE[, 2], + cellType_layer = sce$cellType_layer +) |> + # dplyr::filter((TSNE_2 > 0 & grepl("Excit", cellType_layer)) | !grepl("Excit", cellType_layer)) |> + dplyr::group_by(cellType_layer) |> + dplyr::summarise( + TSNE.1 = median(TSNE_1), + TSNE.2 = median(TSNE_2), + ) +## geom label TSNE_cellType_layer_label <- TSNE_cellType_layer_no_legend + - geom_label(data = tsne_postions, - aes(label = cellType_layer), - size = 5) + geom_label( + data = tsne_postions, + aes(label = cellType_layer), + size = 5 + ) ggsave( - TSNE_cellType_layer_label, - filename = here(plot_dir, "TSNE_cellType_layer_label.png") + TSNE_cellType_layer_label, + filename = here(plot_dir, "TSNE_cellType_layer_label.png") ) -##geom_text +## geom_text TSNE_cellType_layer_text <- TSNE_cellType_layer_no_legend + - geom_text(data = tsne_postions, - aes(label = cellType_layer), - color = "black", - size = 5.5) + geom_text( + data = tsne_postions, + aes(label = cellType_layer), + color = "black", + size = 5.5 + ) ggsave( - TSNE_cellType_layer_text, - filename = here(plot_dir, "TSNE_cellType_layer_text.png") + TSNE_cellType_layer_text, + filename = here(plot_dir, "TSNE_cellType_layer_text.png") ) diff --git a/code/05_explore_sce/02_cellType_prop.R b/code/05_explore_sce/02_cellType_prop.R index 37c8eef..10f8e47 100644 --- a/code/05_explore_sce/02_cellType_prop.R +++ b/code/05_explore_sce/02_cellType_prop.R @@ -1,4 +1,3 @@ - library("SingleCellExperiment") library("DeconvoBuddies") library("tidyverse") @@ -17,10 +16,10 @@ load(here("processed-data", "sce", "sce_DLPFC.Rdata"), verbose = TRUE) ## get proportions before dropping ambig prop_ambig <- as.data.frame(colData(sce)) |> - group_by(cellType_hc, Sample, Position) |> - summarize(n = n()) |> - group_by(Sample) |> - mutate(prop = n / sum(n)) + group_by(cellType_hc, Sample, Position) |> + summarize(n = n()) |> + group_by(Sample) |> + mutate(prop = n / sum(n)) cell_type_colors_ambig <- metadata(sce)$cell_type_colors[levels(prop_ambig$cellType_hc)] @@ -36,7 +35,7 @@ cell_type_colors_layer <- metadata(sce)$cell_type_colors_layer[levels(sce$cellTy pd <- as.data.frame(colData(sce)) table(pd$cellType_broad_hc) -# Astro EndoMural Micro Oligo OPC Excit Inhib +# Astro EndoMural Micro Oligo OPC Excit Inhib # 3979 2157 1601 10894 1940 24809 11067 n_nuc <- pd |> @@ -222,30 +221,33 @@ ggsave(broad_prop_bar_pos, filename = here(plot_dir, "prop_bar_broad_Position.pn #### Plots ambig Plots #### prop_ambig_plus <- prop_ambig |> - mutate(ambig = "Pre-drop") |> - bind_rows(prop_all |> mutate(ambig = "Post-drop")) + mutate(ambig = "Pre-drop") |> + bind_rows(prop_all |> mutate(ambig = "Post-drop")) prop_ambig_bar <- ggplot(data = prop_ambig_plus, aes(x = Sample, y = prop, fill = cellType_hc)) + - geom_bar(stat = "identity") + - geom_text(aes(label = ifelse(prop > 0.02, format(round(prop, 3), 3), ""), - color = as.character(cellType_hc) == "ambig"), - # geom_text(aes(label = ifelse(prop > 0.02, cellType, "")), - size = 2.5, - position = position_stack(vjust = 0.5) - # color = "gray35" - ) + - scale_fill_manual(values = c(cell_type_colors_ambig)) + - scale_color_manual(values = c(`TRUE` = "white", `FALSE` = "black")) + - labs(y = "Proportion", fill = "Cell Type") + - facet_grid(fct_rev(ambig) ~ Position, scales = "free", space = "free") + - # my_theme + - theme_bw()+ - theme(axis.text.x = element_text(angle = 45, hjust = 1)) + - guides(color = FALSE, fill=guide_legend(ncol =1)) - -ggsave(prop_ambig_bar, filename = here(plot_dir,"prop_bar_ambig_Position.png"), width = 11, height = 9) -ggsave(prop_ambig_bar, filename = here(plot_dir,"prop_bar_ambig_Position.pdf"), width = 11, height = 9) + geom_bar(stat = "identity") + + geom_text( + aes( + label = ifelse(prop > 0.02, format(round(prop, 3), 3), ""), + color = as.character(cellType_hc) == "ambig" + ), + # geom_text(aes(label = ifelse(prop > 0.02, cellType, "")), + size = 2.5, + position = position_stack(vjust = 0.5) + # color = "gray35" + ) + + scale_fill_manual(values = c(cell_type_colors_ambig)) + + scale_color_manual(values = c(`TRUE` = "white", `FALSE` = "black")) + + labs(y = "Proportion", fill = "Cell Type") + + facet_grid(fct_rev(ambig) ~ Position, scales = "free", space = "free") + + # my_theme + + theme_bw() + + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + + guides(color = FALSE, fill = guide_legend(ncol = 1)) + +ggsave(prop_ambig_bar, filename = here(plot_dir, "prop_bar_ambig_Position.png"), width = 11, height = 9) +ggsave(prop_ambig_bar, filename = here(plot_dir, "prop_bar_ambig_Position.pdf"), width = 11, height = 9) #### Layer Annotation Proportions #### table(pd$cellType_layer) @@ -305,42 +307,55 @@ ggsave(layer_prop_bar_all, filename = here(plot_dir, "prop_bar_layer_ALL.png"), #### compare hc v layer #### prop_compare <- pd |> - group_by(cellType_layer) |> - summarize(prop = n()/nrow(pd)) |> - mutate(Annotation = "Layer", - cellType = as.character(cellType_layer)) |> - bind_rows(pd |> - group_by(cellType_hc, cellType_layer) |> - summarize(prop = n()/nrow(pd)) |> - mutate(Annotation = "hc", - cellType = as.character(cellType_hc))) |> - replace_na(list(cellType = "Exclude")) |> - group_by(cellType_layer) |> - mutate(order = as.numeric(cellType_layer), - order2 = row_number()) |> - replace_na(list(order = 12.5)) |> ## cheat Excluded cells to be between Excit & Inhib - ungroup() |> - mutate(order = order + (order2*.1), - cellType = factor(cellType, levels=unique(cellType[order(order)]))) - -prop_compare |> filter(Annotation == "hc") |> arrange(cellType_layer)|> print(n=29) + group_by(cellType_layer) |> + summarize(prop = n() / nrow(pd)) |> + mutate( + Annotation = "Layer", + cellType = as.character(cellType_layer) + ) |> + bind_rows(pd |> + group_by(cellType_hc, cellType_layer) |> + summarize(prop = n() / nrow(pd)) |> + mutate( + Annotation = "hc", + cellType = as.character(cellType_hc) + )) |> + replace_na(list(cellType = "Exclude")) |> + group_by(cellType_layer) |> + mutate( + order = as.numeric(cellType_layer), + order2 = row_number() + ) |> + replace_na(list(order = 12.5)) |> ## cheat Excluded cells to be between Excit & Inhib + ungroup() |> + mutate( + order = order + (order2 * .1), + cellType = factor(cellType, levels = unique(cellType[order(order)])) + ) + +prop_compare |> + filter(Annotation == "hc") |> + arrange(cellType_layer) |> + print(n = 29) prop_compare |> filter(is.na(cellType_layer)) prop_compare_bar <- ggplot(data = prop_compare, aes(x = Annotation, y = prop, fill = cellType)) + - geom_bar(stat = "identity", width = .99) + - geom_text(aes(label = ifelse(prop > 0.02, format(round(prop, 3), 3), "")), - # geom_text(aes(label = ifelse(prop > 0.02, cellType, "")), - size = 3, - position = position_stack(vjust = 0.5) - ) + - scale_fill_manual(values = c(cell_type_colors, cell_type_colors_layer)) + - labs(y = "Proportion", fill = "Cell Type") + - my_theme + - scale_y_continuous(limits = c(0, 1), expand = expansion(mult = c(0.01,0.01))) + ## reduce white space in boarder - theme(panel.grid = element_blank(), + geom_bar(stat = "identity", width = .99) + + geom_text(aes(label = ifelse(prop > 0.02, format(round(prop, 3), 3), "")), + # geom_text(aes(label = ifelse(prop > 0.02, cellType, "")), + size = 3, + position = position_stack(vjust = 0.5) + ) + + scale_fill_manual(values = c(cell_type_colors, cell_type_colors_layer)) + + labs(y = "Proportion", fill = "Cell Type") + + my_theme + + scale_y_continuous(limits = c(0, 1), expand = expansion(mult = c(0.01, 0.01))) + ## reduce white space in boarder + theme( + panel.grid = element_blank(), # panel.border = element_blank(), - legend.position = "None") + legend.position = "None" + ) ggsave(prop_compare_bar, filename = here(plot_dir, "prop_compare_bar.png"), width = 2, height = 8) ggsave(prop_compare_bar, filename = here(plot_dir, "prop_compare_bar.pdf"), width = 2, height = 8) @@ -348,7 +363,7 @@ ggsave(prop_compare_bar, filename = here(plot_dir, "prop_compare_bar.pdf"), widt # prop_bar_fineVlayer <- prop_bar_all + theme(legend.position = "None") + labs(y = "Mean Proption - hc annotation k=29") + # layer_prop_bar_all + labs(y = "Mean Proption - layer annotation k=13") -# +# # ggsave(prop_bar_fineVlayer, filename = here(plot_dir, "prop_bar_fineVlayer_ALL.png")) # ggsave(prop_bar_fineVlayer, filename = here(plot_dir, "prop_bar_fineVlayer_ALL.pdf")) @@ -378,8 +393,10 @@ layer_prop_bar_position <- ggplot(data = prop_layer, aes(x = Sample, y = prop, f scale_fill_manual(values = metadata(sce)$cell_type_colors_layer) + labs(y = "Cell Type Proportion", fill = "Cell Type") + my_theme + - theme(axis.text.x = element_text(angle = 45, hjust = 1), - legend.position = "None") + theme( + axis.text.x = element_text(angle = 45, hjust = 1), + legend.position = "None" + ) ggsave(layer_prop_bar_position, filename = here(plot_dir, "prop_bar_layer_position.png"), width = 10) ggsave(layer_prop_bar_position, filename = here(plot_dir, "prop_bar_layer_position.pdf"), width = 10) diff --git a/code/05_explore_sce/04_3D_UMAP.R b/code/05_explore_sce/04_3D_UMAP.R index f0fe124..7118fed 100644 --- a/code/05_explore_sce/04_3D_UMAP.R +++ b/code/05_explore_sce/04_3D_UMAP.R @@ -1,4 +1,3 @@ - library("SingleCellExperiment") library("scater") # library("tidyverse") diff --git a/code/05_explore_sce/06_explore_azimuth_annotations.R b/code/05_explore_sce/06_explore_azimuth_annotations.R index e2625b1..7555368 100644 --- a/code/05_explore_sce/06_explore_azimuth_annotations.R +++ b/code/05_explore_sce/06_explore_azimuth_annotations.R @@ -119,53 +119,63 @@ dev.off() ## layer colors from spatialLIBD + intermediate layers libd_intermediate_layer_colors <- c( - "L1/2" = "#BF3889", - "L2/3" = "#50DDAC", - "L3/4" = "#8278B0", - "L3/4/5" = "#7D80CA", #Excit 3/4/5 - "L4/5" = "#BD8339", - "L5/6" = "#FFB300", - "L6/WM" = "#7A3D00", - "L1/WM" = "#650136", # Glia cells? - # "WM/L1" = "#650136", # Glia cells? - "No Assigment" = "gray" + "L1/2" = "#BF3889", + "L2/3" = "#50DDAC", + "L3/4" = "#8278B0", + "L3/4/5" = "#7D80CA", # Excit 3/4/5 + "L4/5" = "#BD8339", + "L5/6" = "#FFB300", + "L6/WM" = "#7A3D00", + "L1/WM" = "#650136", # Glia cells? + # "WM/L1" = "#650136", # Glia cells? + "No Assigment" = "gray" ) libd_intermediate_layer_colors <- - c( - spatialLIBD::libd_layer_colors, - libd_intermediate_layer_colors - ) + c( + spatialLIBD::libd_layer_colors, + libd_intermediate_layer_colors + ) names(libd_intermediate_layer_colors) <- - gsub("ayer", "", names(libd_intermediate_layer_colors)) + gsub("ayer", "", names(libd_intermediate_layer_colors)) libd_intermediate_layer_colors ## Spatial Registration -hc_annotations <- pd |> - group_by(cellType_hc, cellType_broad_hc, layer_annotation) |> - count() |> - mutate(Layer = as.character(layer_annotation), - Layer = ifelse(grepl("\\*", Layer),NA,Layer)) |> - ungroup() |> - select(cellType_hc, Layer) +hc_annotations <- pd |> + group_by(cellType_hc, cellType_broad_hc, layer_annotation) |> + count() |> + mutate( + Layer = as.character(layer_annotation), + Layer = ifelse(grepl("\\*", Layer), NA, Layer) + ) |> + ungroup() |> + select(cellType_hc, Layer) unique(hc_annotations$Layer) -row_ha <- rowAnnotation(df = as.data.frame(hc_annotations), - col = list(Layer = libd_intermediate_layer_colors, - cellType_hc = cell_type_colors), - show_legend = c(FALSE, TRUE)) +row_ha <- rowAnnotation( + df = as.data.frame(hc_annotations), + col = list( + Layer = libd_intermediate_layer_colors, + cellType_hc = cell_type_colors + ), + show_legend = c(FALSE, TRUE) +) # hc_count <- rowAnnotation(n_nuc = anno_barplot(as.numeric(table(sce$cellType_hc)))) ## azimuth annotation azimuth_anno <- azimuth_cellType_notes # az_count <- columnAnnotation(n_nuc = anno_barplot(as.numeric(table(sce$cellType_azimuth)))) -col_ha <- HeatmapAnnotation(df = azimuth_anno |> select(-azimuth), - col = list(Layer = libd_intermediate_layer_colors, - cellType_broad = cell_type_colors_broad), - annotation_name_side = "left", - show_legend = c(TRUE, FALSE)) +col_ha <- HeatmapAnnotation( + df = azimuth_anno |> select(-azimuth), + col = list( + Layer = libd_intermediate_layer_colors, + cellType_broad = cell_type_colors_broad + ), + annotation_name_side = "left", + show_legend = c(TRUE, FALSE) +) #### Complex Heatmap #### @@ -182,7 +192,7 @@ Heatmap(jacc.mat, # right_annotation = hc_count, bottom_annotation = col_ha, # bottom_annotation = az_count, - col = c("black",viridisLite::plasma(100)), + col = c("black", viridisLite::plasma(100)), na_col = "black" ) dev.off() @@ -191,30 +201,34 @@ dev.off() pd_layer <- pd |> filter(!is.na(cellType_layer)) jacc.mat.layer <- linkClustersMatrix(pd_layer$cellType_layer, pd_layer$cellType_azimuth) -layer_annotations <- pd_layer |> - group_by(cellType_layer) |> - count() |> - ## what to do if cell types have multiple layer assignments? - mutate(Layer = ifelse(grepl("Excit_L", cellType_layer),gsub("Excit_","", cellType_layer),NA)) |> - ungroup() |> - select(cellType_layer, Layer) - -row_ha_layer <- rowAnnotation(df = as.data.frame(layer_annotations), - col = list(Layer = libd_intermediate_layer_colors, - cellType_layer = metadata(sce)$cell_type_colors_layer), - show_legend = c(FALSE, TRUE)) +layer_annotations <- pd_layer |> + group_by(cellType_layer) |> + count() |> + ## what to do if cell types have multiple layer assignments? + mutate(Layer = ifelse(grepl("Excit_L", cellType_layer), gsub("Excit_", "", cellType_layer), NA)) |> + ungroup() |> + select(cellType_layer, Layer) + +row_ha_layer <- rowAnnotation( + df = as.data.frame(layer_annotations), + col = list( + Layer = libd_intermediate_layer_colors, + cellType_layer = metadata(sce)$cell_type_colors_layer + ), + show_legend = c(FALSE, TRUE) +) jacc.mat.layer <- jacc.mat.layer[layer_annotations$cellType_layer, azimuth_anno$azimuth] pdf(here(plot_dir, "azimuth_v_layer_annotation_complex.pdf"), height = 8, width = 12) # png(here(plot_dir, "azimuth_v_hc_annotation_complex.png"), height = 800, width = 800) Heatmap(jacc.mat.layer, - name = "Correspondence", - right_annotation = row_ha_layer, - # right_annotation = hc_count, - bottom_annotation = col_ha, - # bottom_annotation = az_count, - col = c("black",viridisLite::plasma(100)), + name = "Correspondence", + right_annotation = row_ha_layer, + # right_annotation = hc_count, + bottom_annotation = col_ha, + # bottom_annotation = az_count, + col = c("black", viridisLite::plasma(100)), ) dev.off() diff --git a/code/05_explore_sce/07_convert_SCE_HDF5.R b/code/05_explore_sce/07_convert_SCE_HDF5.R index a9af6c1..7c78844 100644 --- a/code/05_explore_sce/07_convert_SCE_HDF5.R +++ b/code/05_explore_sce/07_convert_SCE_HDF5.R @@ -5,9 +5,9 @@ library("sessioninfo") ## load data load(here("processed-data", "sce", "sce_DLPFC.Rdata"), verbose = TRUE) sce <- HDF5Array::loadHDF5SummarizedExperiment(here("processed-data", "sce", "sce_DLPFC")) -rhdf5::h5ls(here("processed-data", "sce", "sce_DLPFC","assays.h5")) -sce_h5 <- rhdf5::h5read(here("processed-data", "sce", "sce_DLPFC","assays.h5"), "assay001") -sce <- zellkonverter::readH5AD(here("processed-data", "sce", "sce_DLPFC","assays.h5")) +rhdf5::h5ls(here("processed-data", "sce", "sce_DLPFC", "assays.h5")) +sce_h5 <- rhdf5::h5read(here("processed-data", "sce", "sce_DLPFC", "assays.h5"), "assay001") +sce <- zellkonverter::readH5AD(here("processed-data", "sce", "sce_DLPFC", "assays.h5")) names(assays(sce)) # counts binomial_deviance_residuals logcounts diff --git a/code/05_explore_sce/08_pseudobulk_cellType.R b/code/05_explore_sce/08_pseudobulk_cellType.R index ec69827..b4324c1 100644 --- a/code/05_explore_sce/08_pseudobulk_cellType.R +++ b/code/05_explore_sce/08_pseudobulk_cellType.R @@ -1,4 +1,3 @@ - library("SingleCellExperiment") library("scuttle") library("edgeR") diff --git a/code/05_explore_sce/09_SCZ_interaction_expression.R b/code/05_explore_sce/09_SCZ_interaction_expression.R index 6750e9f..f5a02dd 100644 --- a/code/05_explore_sce/09_SCZ_interaction_expression.R +++ b/code/05_explore_sce/09_SCZ_interaction_expression.R @@ -16,9 +16,9 @@ if (!dir.exists(plot_dir)) dir.create(plot_dir) load(here("processed-data", "05_explore_sce", "08_pseudobulk_cellTypes", "sce_pseudo-cellType_layer.Rdata"), verbose = TRUE) ## Check out colData table(sce_pseudo$cellType_layer) -# Astro EndoMural Micro Oligo OPC Excit_L2/3 Excit_L3 Excit_L3/4/5 Excit_L4 -# 18 17 18 19 19 1 19 17 16 -# Excit_L5 Excit_L5/6 Excit_L6 Inhib +# Astro EndoMural Micro Oligo OPC Excit_L2/3 Excit_L3 Excit_L3/4/5 Excit_L4 +# 18 17 18 19 19 1 19 17 16 +# Excit_L5 Excit_L5/6 Excit_L6 Inhib # 16 16 16 19 @@ -39,43 +39,40 @@ colnames(cat) <- c("sample_id", "cat") expression_long <- cbind(expression_long, cat) pdf(here(plot_dir, "spatialDLPFC_SCZ_interactions.pdf")) -for(i in 1:nrow(SCZ_interactions)){ - - interaction <- SCZ_interactions[i,] - interaction_title <- paste(interaction[[1]], "->", interaction[[2]]) - message(i, " ", interaction_title) - - expression_temp <- expression_long |> - filter(Var1 %in% interaction) |> - mutate(Var1 = factor(Var1, levels = interaction)) - - expression_violin <- ggplot(data = expression_temp, aes(x = cat, y = value, fill = cat)) + - geom_violin(scale = "width", alpha = 0.5) + - geom_jitter(aes(color = cat), size = .3) + - facet_wrap(~Var1, ncol = 1, scales = "free_y") + - labs( - title = interaction_title, - y = "Expression (logcounts)" - ) + - theme_bw() + - theme( - legend.position = "None", - axis.title.x = element_blank(), - axis.text.x = element_text(angle = 90, hjust = 1), - strip.text.x = element_text(face = "italic"), - text = element_text(size = 15) - ) + - stat_summary( - fun = median, - # fun.min = median, - # fun.max = median, - geom = "crossbar", - width = 0.3 - ) + - scale_color_manual(values = metadata(sce_pseudo)$cell_type_colors_layer) + - scale_fill_manual(values = metadata(sce_pseudo)$cell_type_colors_layer) - print(expression_violin) -} -dev.off() +for (i in 1:nrow(SCZ_interactions)) { + interaction <- SCZ_interactions[i, ] + interaction_title <- paste(interaction[[1]], "->", interaction[[2]]) + message(i, " ", interaction_title) + expression_temp <- expression_long |> + filter(Var1 %in% interaction) |> + mutate(Var1 = factor(Var1, levels = interaction)) + expression_violin <- ggplot(data = expression_temp, aes(x = cat, y = value, fill = cat)) + + geom_violin(scale = "width", alpha = 0.5) + + geom_jitter(aes(color = cat), size = .3) + + facet_wrap(~Var1, ncol = 1, scales = "free_y") + + labs( + title = interaction_title, + y = "Expression (logcounts)" + ) + + theme_bw() + + theme( + legend.position = "None", + axis.title.x = element_blank(), + axis.text.x = element_text(angle = 90, hjust = 1), + strip.text.x = element_text(face = "italic"), + text = element_text(size = 15) + ) + + stat_summary( + fun = median, + # fun.min = median, + # fun.max = median, + geom = "crossbar", + width = 0.3 + ) + + scale_color_manual(values = metadata(sce_pseudo)$cell_type_colors_layer) + + scale_fill_manual(values = metadata(sce_pseudo)$cell_type_colors_layer) + print(expression_violin) +} +dev.off()