suppressPackageStartupMessages({
library(SingleCellExperiment)
library(scuttle)
library(scater)
library(dplyr)
library(tidyr)
library(ggplot2)
library(cowplot)
library(ComplexHeatmap)
library(circlize)
library(ggtext)
library(scales)
library(ggh4x)
library(swissknife)
library(BiocParallel)
library(CellChat)
library(tibble)
library(patchwork)
library(ggforce)
library(ggtext)
})
source("../params.R")
source("../helpers.R")
if (!file.exists("paper_figures")) {
dir.create("paper_figures", showWarnings = FALSE, recursive = FALSE)
}Generating plots for paper (cell-cell communication)
Load packages
Load data
atlas <- readRDS("invivo_cell_atlas.rds")
cellchat <- readRDS("cellchat_menstrual_proliferative.rds")
organ <- readRDS("../scRNAseq_timecourse/organoid_scrnaseq_timecourse.rds")
# pathways of interest
poi <- c("CXCL", "ANNEXIN", "CSF", "EDN", "SEMA3", "IL6", "WNT")
sender <- "luminal early"
receivers <- setdiff(as.character(cellchat@idents), sender)Heatmaps of interaction strength by pathway
# cellchat@netP$prob give interaction probabilities
# dim 1 - sender, dim 2 - receiver, dim 3 - pathway
# luminal early as sender
plotdf <- as.data.frame(cellchat@netP$prob["luminal early", , poi]) |>
rownames_to_column("receiver") |>
pivot_longer(names_to = "pathway", values_to = "prob", -receiver) |>
group_by(pathway) |>
mutate(prob_scaled = prob / max(prob)) |>
ungroup() |>
mutate(receiver = factor(receiver,
levels = names(tissue_all_lineages_celltype_cols)[
names(tissue_all_lineages_celltype_cols) %in% receiver]
)) |>
mutate(pathway = factor(pathway, levels = poi))
(ggs <- ggplot(plotdf,
aes(x = receiver, y = pathway, fill = prob_scaled)) +
geom_tile() +
scale_fill_gradient(low = "grey98", high = "firebrick",
limits = c(0, 1), oob = scales::squish,
name = "relative\ninteraction\nprobability") +
scale_x_discrete(labels = tissue_all_lineages_celltype_labels) +
labs(x = "Receiving population", y = "Pathway",
title = "regeneration associated luminal as sender") +
theme_void() +
theme(axis.text.x = element_markdown(angle = 90, hjust = 1, vjust = 0.5),
axis.text.y = element_markdown()))
pdf("paper_figures/cellchat_pathways_luminal_early_sender.pdf",
width = 7, height = 5)
ggs
dev.off()png
2 # cellchat@netP$prob give interaction probabilities
# dim 1 - sender, dim 2 - receiver, dim 3 - pathway
# luminal early as sender
# cluster pathways
tmp <- cellchat@netP$prob["luminal early", , ]
tmp <- tmp[, colSums(tmp) > 0, drop = FALSE]
tmp <- sweep(tmp, MARGIN = 2, STATS = colSums(tmp), FUN = "/")
hcl <- colnames(tmp)[hclust(d = dist(t(tmp)), method = "complete")$order]
plotdf <- as.data.frame(cellchat@netP$prob["luminal early", , ]) |>
rownames_to_column("receiver") |>
pivot_longer(names_to = "pathway", values_to = "prob", -receiver) |>
group_by(pathway) |>
mutate(prob_scaled = prob / max(prob)) |>
ungroup() |>
filter(!is.na(prob_scaled)) |>
mutate(receiver = factor(receiver,
levels = names(tissue_all_lineages_celltype_cols)[
names(tissue_all_lineages_celltype_cols) %in% receiver]),
pathway = factor(pathway,
levels = hcl[hcl %in% pathway])
)
(ggs <- ggplot(plotdf,
aes(x = pathway, y = receiver, fill = prob_scaled)) +
geom_tile() +
scale_fill_gradient(low = "grey98", high = "firebrick",
limits = c(0, 1), oob = scales::squish,
name = "relative\ninteraction\nprobability") +
scale_y_discrete(labels = tissue_all_lineages_celltype_labels) +
labs(y = "Receiving population", x = "Pathway",
title = "regeneration associated luminal as sender") +
theme_void() +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
axis.text.y = element_markdown(hjust = 1)))
pdf("paper_figures/cellchat_pathways_luminal_early_sender_all.pdf",
width = 7, height = 6)
ggs
dev.off()png
2 png("paper_figures/cellchat_pathways_luminal_early_sender_all.png",
width = 7, height = 6, units = "in", res = 400)
ggs
dev.off()png
2 Bar plot of total interaction strength from luminal early to cell types
# cellchat@net$weight give interaction probabilities
# dim 1 - sender, dim 2 - receiver
# luminal early as sender
plotdf <- as.data.frame(cellchat@net$weight["luminal early", , drop = FALSE]) |>
rownames_to_column("sender") |>
pivot_longer(names_to = "receiver", values_to = "prob", -sender) |>
mutate(lineage = tissue_all_lineages_celltype_lineage[receiver]) |>
mutate(receiver = factor(receiver, levels = names(tissue_all_lineages_celltype_cols)))
ggc1 <- ggplot(plotdf |> dplyr::filter(lineage %in% c("Endothelial", "Immune")),
# aes(x = receiver, y = prob, fill = receiver)) +
aes(x = receiver, y = prob)) +
geom_col(fill = "grey80") +
theme_bw() +
# scale_fill_manual(values = tissue_all_lineages_celltype_cols) +
scale_x_discrete(labels = tissue_all_lineages_celltype_labels) +
facet_col(~ lineage, scales = "free_y", space = "free") +
theme(axis.text.x = element_markdown(angle = 90, hjust = 1, vjust = 0.5),
axis.text.y = element_markdown(),
legend.position = "none") +
coord_flip(ylim = c(0, max(plotdf$prob))) +
labs(x = "Receiver", y = "Summed interaction probabilities")
ggc2 <- ggplot(plotdf |> dplyr::filter(lineage %in% c("Epithelial", "Mesenchymal")),
# aes(x = receiver, y = prob, fill = receiver)) +
aes(x = receiver, y = prob)) +
geom_col(fill = "grey80") +
theme_bw() +
# scale_fill_manual(values = tissue_all_lineages_celltype_cols) +
scale_x_discrete(labels = tissue_all_lineages_celltype_labels) +
facet_col(~ lineage, scales = "free_y", space = "free") +
theme(axis.text.x = element_markdown(angle = 90, hjust = 1, vjust = 0.5),
axis.text.y = element_markdown(),
legend.position = "none") +
coord_flip(ylim = c(0, max(plotdf$prob))) +
labs(x = "Receiver", y = "Summed interaction probabilities")
ggc1 + ggc2 + plot_layout(axis_titles = "collect")
pdf("paper_figures/cellchat_celltypes_luminal_early_sender.pdf",
width = 7, height = 6)
ggc1 + ggc2 + plot_layout(axis_titles = "collect")
dev.off()png
2 png("paper_figures/cellchat_celltypes_luminal_early_sender.png",
width = 7, height = 6, units = "in", res = 400)
ggc1 + ggc2 + plot_layout(axis_titles = "collect")
dev.off()png
2 Plot expression of ligands and receptors for pathways of interest
lr <- subsetCommunication(
cellchat, slot.name = "net",
sources.use = "luminal early", targets.use = receivers,
signaling = NULL, pairLR.use = NULL, thresh = 0.05
)
allgray <- setNames(rep("gray10", length(tissue_all_lineages_celltype_cols)),
names(tissue_all_lineages_celltype_cols))
for (pw in unique(lr$pathway_name)) {
lrsub <- lr |>
filter(pathway_name == pw)
markers <- unique(c(lrsub$ligand, unlist(strsplit(lrsub$receptor, "_"))))
if (pw == "WNT") {
wntligands <- cellchat@DB$interaction |>
filter(pathway_name == "WNT") |>
pull(ligand) |>
unique()
markers <- unique(c(wntligands, markers, paste0("FZD", 1:10)))
markers <- markers[rowSums(assay(atlas, "logcounts")[markers, ]) > 0]
} else if (pw == "CXCL") {
markers <- unique(c(markers, c("CXCR1", "CXCR2")))
}
dp <- .makeDotPlot(
atlas, markers = markers, cluster_column = "celltype_annot",
cluster_order = names(tissue_all_lineages_celltype_cols),
cluster_colors = allgray,
cluster_labels = tissue_all_lineages_celltype_labels,
transpose = TRUE, split_by_ciliated = FALSE, normalize = TRUE,
min_fraction = 0.2) +
ggtitle(pw)
print(dp)
pdf(paste0("paper_figures/cellchat_pathways_ligands_receptors_", pw, ".pdf"),
height = 9, width = min(5.5, 1.25 * (0.2 * length(markers) + 3)))
print(dp)
dev.off()
}
Dotplot of WNT+Annexin+CXCL pathways
lrsubWNT <- lr |>
filter(pathway_name == "WNT")
lrsubANNEXIN <- lr |>
filter(pathway_name == "ANNEXIN")
lrsubCXCL <- lr |>
filter(pathway_name == "CXCL")
markersWNT <- unique(c(lrsubWNT$ligand, unlist(strsplit(lrsubWNT$receptor, "_"))))
markersANNEXIN <- unique(c(lrsubANNEXIN$ligand, unlist(strsplit(lrsubANNEXIN$receptor, "_"))))
markersCXCL <- unique(c(lrsubCXCL$ligand, unlist(strsplit(lrsubCXCL$receptor, "_"))))
wntligands <- cellchat@DB$interaction |>
filter(pathway_name == "WNT") |>
pull(ligand) |>
unique()
markers <- unique(c(wntligands, markersWNT, paste0("FZD", 1:10), markersANNEXIN,
markersCXCL, "CXCR1", "CXCR2"))
markers <- markers[rowSums(assay(atlas, "logcounts")[markers, ]) > 0]
dp <- .makeDotPlot(
atlas, markers = markers, cluster_column = "celltype_annot",
cluster_order = names(tissue_all_lineages_celltype_cols),
cluster_colors = allgray,
cluster_labels = tissue_all_lineages_celltype_labels,
transpose = TRUE, split_by_ciliated = FALSE, normalize = TRUE,
min_fraction = 0.2)
print(dp)
pdf(paste0("paper_figures/cellchat_pathways_ligands_receptors_WNT_ANNEXIN_CXCL.pdf"),
height = 9, width = 8)
print(dp)
dev.off()png
2 
Dotplot of all ligands in organoid timecourse
(ligands <- subsetCommunication(
cellchat, slot.name = "net",
sources.use = "luminal early", targets.use = receivers,
signaling = NULL, pairLR.use = NULL, thresh = 0.05) |>
pull(ligand) |>
unique()) [1] "GDF15" "WNT7A" "AREG" "HBEGF" "PDGFA" "PDGFB" "VEGFA"
[8] "CXCL1" "CXCL2" "CXCL3" "CXCL8" "CXCL16" "MIF" "CX3CL1"
[15] "IL6" "LIF" "CSF1" "CSF3" "TNF" "TNFSF10" "NAMPT"
[22] "MDK" "PTN" "C3" "EDN1" "EDN2" "SEMA3A" "SEMA3B"
[29] "SEMA3C" "ADM" "ANXA1" "GAS6" "GRN" "LGALS9" "RARRES2"
[36] "BAG6" (dotplot_norm_ligands <- .makeDotPlot(
dat = organ, markers = ligands,
cluster_column = "leiden_res0.125_annot",
cluster_order = sc_leiden_res0.125_annot_order,
cluster_colors = sc_leiden_res0.125_annot_colors,
cluster_labels = sc_leiden_res0.125_annot_labels,
transpose = TRUE,
split_by_ciliated = FALSE))Warning: Removed 171 rows containing missing values or values outside the scale range
(`geom_point()`).
pdf("paper_figures/dotplot_organoids_ligands_normalized.pdf", height = 6, width = 10)
dotplot_norm_ligandsWarning: Removed 171 rows containing missing values or values outside the scale range
(`geom_point()`).dev.off()png
2 png("paper_figures/dotplot_organoids_ligands_normalized.png", height = 6, width = 10,
unit = "in", res = 500)
dotplot_norm_ligandsWarning: Removed 171 rows containing missing values or values outside the scale range
(`geom_point()`).dev.off()png
2 Session info
sessionInfo()R version 4.4.2 (2024-10-31)
Platform: x86_64-pc-linux-gnu
Running under: Rocky Linux 8.10 (Green Obsidian)
Matrix products: default
BLAS/LAPACK: FlexiBLAS OPENBLAS; LAPACK version 3.12.0
locale:
[1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
[4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
[7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
[10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
time zone: Europe/Zurich
tzcode source: system (glibc)
attached base packages:
[1] grid stats4 stats graphics grDevices utils datasets
[8] methods base
other attached packages:
[1] ggforce_0.4.2 patchwork_1.3.0
[3] tibble_3.2.1 CellChat_1.6.1
[5] bigmemory_4.6.4 igraph_2.1.4
[7] BiocParallel_1.40.2 swissknife_0.43
[9] ggh4x_0.3.0 scales_1.3.0.9000
[11] ggtext_0.1.2 circlize_0.4.16
[13] ComplexHeatmap_2.22.0 cowplot_1.1.3
[15] tidyr_1.3.1 dplyr_1.1.4
[17] scater_1.34.1 ggplot2_3.5.2
[19] scuttle_1.16.0 SingleCellExperiment_1.28.1
[21] SummarizedExperiment_1.36.0 Biobase_2.66.0
[23] GenomicRanges_1.58.0 GenomeInfoDb_1.42.3
[25] IRanges_2.40.1 S4Vectors_0.44.0
[27] BiocGenerics_0.52.0 MatrixGenerics_1.18.1
[29] matrixStats_1.5.0
loaded via a namespace (and not attached):
[1] RColorBrewer_1.1-3 jsonlite_2.0.0 shape_1.4.6.1
[4] magrittr_2.0.3 ggbeeswarm_0.7.2 farver_2.1.2
[7] rmarkdown_2.29 GlobalOptions_0.1.2 fs_1.6.6
[10] zlibbioc_1.52.0 vctrs_0.6.5 rstatix_0.7.2
[13] htmltools_0.5.8.1 S4Arrays_1.6.0 usethis_3.1.0
[16] BiocNeighbors_2.0.1 broom_1.0.8 Formula_1.2-5
[19] SparseArray_1.6.2 parallelly_1.43.0 KernSmooth_2.23-24
[22] htmlwidgets_1.6.4 plyr_1.8.9 uuid_1.2-1
[25] commonmark_1.9.5 lifecycle_1.0.4 ggnetwork_0.5.13
[28] iterators_1.0.14 pkgconfig_2.0.3 rsvd_1.0.5
[31] Matrix_1.7-1 R6_2.6.1 fastmap_1.2.0
[34] GenomeInfoDbData_1.2.13 future_1.40.0 clue_0.3-66
[37] digest_0.6.37 colorspace_2.1-1 dqrng_0.4.1
[40] RSpectra_0.16-2 irlba_2.3.5.1 ggpubr_0.6.0
[43] beachmat_2.22.0 labeling_0.4.3 polyclip_1.10-7
[46] httr_1.4.7 abind_1.4-8 compiler_4.4.2
[49] rngtools_1.5.2 withr_3.0.2 doParallel_1.0.17
[52] backports_1.5.0 carData_3.0-5 viridis_0.6.5
[55] ggsignif_0.6.4 MASS_7.3-61 DelayedArray_0.32.0
[58] bluster_1.16.0 rjson_0.2.23 tools_4.4.2
[61] vipor_0.4.7 beeswarm_0.4.0 future.apply_1.11.3
[64] glue_1.8.0 gridtext_0.1.5 gridBase_0.4-7
[67] cluster_2.1.6 reshape2_1.4.4 generics_0.1.3
[70] gtable_0.3.6 sna_2.8 metapod_1.14.0
[73] car_3.1-3 BiocSingular_1.22.0 ScaledMatrix_1.14.0
[76] xml2_1.3.8 XVector_0.46.0 markdown_2.0
[79] ggrepel_0.9.6 foreach_1.5.2 pillar_1.10.2
[82] stringr_1.5.1 limma_3.62.2 tweenr_2.0.3
[85] lattice_0.22-6 FNN_1.1.4.1 tidyselect_1.2.1
[88] registry_0.5-1 locfit_1.5-9.12 pbapply_1.7-2
[91] knitr_1.50 bigmemory.sri_0.1.8 gridExtra_2.3
[94] litedown_0.7 edgeR_4.4.2 svglite_2.1.3
[97] xfun_0.52 statmod_1.5.0 NMF_0.28
[100] stringi_1.8.7 UCSC.utils_1.2.0 statnet.common_4.11.0
[103] yaml_2.3.10 evaluate_1.0.3 codetools_0.2-20
[106] BiocManager_1.30.25 cli_3.6.4 systemfonts_1.2.2
[109] reticulate_1.42.0 network_1.19.0 dichromat_2.0-0.1
[112] Rcpp_1.0.14 globals_0.17.0 coda_0.19-4.1
[115] png_0.1-8 parallel_4.4.2 ggalluvial_0.12.5
[118] scran_1.34.0 listenv_0.9.1 viridisLite_0.4.2
[121] purrr_1.0.4 crayon_1.5.3 GetoptLong_1.0.5
[124] rlang_1.1.6