These libraries and objects are used in Figures 2, 4, and 5.
# Unique list of library calls
library(cowplot)
library(tidyverse)
library(grid)
library(gridExtra)
library(muscat)
library(purrr)
library(SingleCellExperiment)
library(UpSetR)
library(Seurat)
library(ComplexHeatmap)
library(circlize)
library(viridis)
library(ggridges)
library(dplyr)
library(imager)
so.astro <- readRDS("~/data/Followup/so.astro.rds")
so.astro@meta.data$ident <- as.factor(so.astro@meta.data$integrated_snn_res.0.2)
# Prepare a vector of colors for the clusters
clusters <- levels(Idents(so.astro)) %>% as.numeric()
cluster.cols <- muscat:::.cluster_colors[seq_along(clusters)]
names(cluster.cols) <- clusters
clusters
## [1] 0 1 2 3 4 5 6 7 8 9
astro.markers <- readRDS("~/data/Followup/astro.markers.v1.rds") %>%
remove_rownames() %>%
mutate(gene.name = word(gene, 2, sep = "\\."))
head(astro.markers)
Qc and validation of six Visium spatial transcriptomics sections from saline- and LPS-treated animals
:
Sample to GEO ID Mapping:
GSM5026144_S7
)GSM5026146_S4
)GSM5026145_S3
)GSM5026149_S6
)GSM5026148_S2
)GSM5026147_S1
)# Load the data into a Seurat object
## Saline_1 (GSM5026144_S7), used in Figures 2, 4, and 5
GSM5026144_S7_Seurat <- Load10X_Spatial(data.dir = "~/data/Followup/GSM5026144_S7/", filename = "GSM5026144_S7_filtered_feature_bc_matrix.h5", assay = "Spatial")
## LPS_2 (GSM5026148_S2), used in Figures 4 and 5
GSM5026148_S2_Seurat <- Load10X_Spatial(data.dir = "~/data/Followup/GSM5026148_S2/", filename = "GSM5026148_S2_filtered_feature_bc_matrix.h5", assay = "Spatial")
# Normalize the data with `SCTransform()` function
GSM5026144_S7_Seurat <- SCTransform(GSM5026144_S7_Seurat, assay = 'Spatial', return.only.var.genes = FALSE, verbose = FALSE)
GSM5026148_S2_Seurat<- SCTransform(GSM5026148_S2_Seurat, assay = 'Spatial', return.only.var.genes = FALSE, verbose = FALSE)
DimPlot(so.astro, reduction = "tsne", group.by = "ident", split.by = "group_id") + scale_color_manual(values = cluster.cols)
# Subset of genes used for plotting
astro.markers.subset <- astro.markers %>% group_by(cluster) %>% slice_max(n = 15, order_by = avg_log2FC)
astro.markers.subset
# Expression matrix
cells_to_clusters <- data.frame(cluster = so.astro$integrated_snn_res.0.2) %>% rownames_to_column("cell")
geneset <- unique(astro.markers.subset$gene)
mat.scale.ours_data <- so.astro@assays$integrated@scale.data[geneset, cells_to_clusters$cell]
cluster_indices <- split(cells_to_clusters$cell, cells_to_clusters$cluster)
cluster_means <- sapply(cluster_indices, function(indices) {Matrix::rowMeans(mat.scale.ours_data[, indices, drop = FALSE])})
rownames(cluster_means) <- word(rownames(cluster_means), 2, sep = "\\.")
# Gene cluster color bar annotation
# by cluster order
genes.clusters.cols <- astro.markers.subset %>% select(gene, cluster) %>% mutate(gene.full = gene, gene = word(gene.full, 2, sep = "\\.")) %>% group_by(gene) %>% slice_min(n = 1, order_by = cluster) %>% ungroup() %>% arrange(cluster) %>% left_join(., enframe(cluster.cols, name = "cluster", value = "colors"))
rownames(genes.clusters.cols) <- genes.clusters.cols$gene
genes.clusters.cols <- genes.clusters.cols[rownames(cluster_means),]
row.anno.genes <- rowAnnotation("Cluster" = genes.clusters.cols$cluster, col = list("Cluster" = cluster.cols), show_annotation_name = FALSE)
# Define the range and breaks
breaks <- seq(0, 1, length.out = 11) # This creates 11 points between 0 and 1
colors <- plasma(11)
# Create a color function that maps the range [0, 1] to the viridis palette
color_mapping <- colorRamp2(breaks, colors)
# Use the adjusted col argument
Heatmap(cluster_means, name = "Normal. expression", col = color_mapping, cluster_rows = FALSE, show_row_dend = FALSE, cluster_columns = FALSE, column_names_side = "bottom", row_names_side = "left", column_names_rot = 0, column_names_centered = TRUE, row_names_gp = gpar(fontsize = 8), left_annotation = row.anno.genes)
astrocyte_markers <- c("Aldh1l1", "Slc1a2", "S100b", "Gfap", "Myoc")
paper_cluster_0.module <- list(c("Gm42418", "Arhgap31", "Kcnq1ot1", "Rmst", "Cit", "Gli1", "Ttll3", "Neat1", "Gm26917"))
paper_cluster_1.module <- list(c("Mfge8", "Igfbp2", "Btbd17"))
paper_cluster_2.module <- list(c("Dbi", "Agt", "Nnat", "Ppia", "Tmsb4x", "Scrg1", "Mt3", "Gstm1", "Nkain4", "S100a1"))
paper_cluster_3.module <- list(c("Thrsp", "Actb", "Chchd10", "Hint1", "Mrpl41", "Mpc1", "Camk2n1", "Ndufc2", "Nupr1", "Igfbp2"))
paper_cluster_4.module <- list(c("Lcn2", "Ifitm3", "Timp1", "B2m", "Vim", "Gfap", "Psmb8", "Bst2"))
paper_cluster_5.module <- list(c("Crym", "Prss56", "Slc43a3", "Mgst1", "Slc1a5", "Tuba1a", "Itm2b", "Chil1", "Ptn"))
paper_cluster_6.module <- list(c("Thbs4", "Igfbp5", "Fxyd6", "Gsn", "Fxyd1", "Cd9", "Meg3"))
paper_cluster_7.module <- list(c("Itih3", "Sparc", "Nkx6-2", "Etnppl", "Gria1", "Spon1", "5031439G07Rik"))
paper_cluster_8.module <- list(c("Igtp", "Gm4951", "Ifit3", "Iigp1", "Irgm1", "Tap1", "Gbp2", "Ifit1", "Isg15", "Cxcl10"))
paper_cluster_9.module <- list(c("Zic1", "Slc38a1", "Col1a2", "Kcnk2", "Luzp2", "Mgll", "Zic4", "Hopx"))
paper_cluster_modules.all <- c("paper_cluster_0.module", "paper_cluster_1.module", "paper_cluster_2.module", "paper_cluster_3.module", "paper_cluster_4.module", "paper_cluster_5.module", "paper_cluster_6.module", "paper_cluster_7.module", "paper_cluster_8.module", "paper_cluster_9.module")
for (gene in astrocyte_markers) {
GSM5026144_S7_Seurat <- AddModuleScore(object = GSM5026144_S7_Seurat, features = list(gene), name = gene)
}
for (cluster_module_name in paper_cluster_modules.all) {
# Retrieve the actual module list using get()
cluster_module <- get(cluster_module_name)
# Add the module score to the Seurat object
GSM5026144_S7_Seurat <- AddModuleScore(object = GSM5026144_S7_Seurat, features = cluster_module, name = cluster_module_name)
}
## Warning: The following features are not present in the object: Gm4951, not
## searching for symbol synonyms
# List to store individual plots
paper_plots <- list()
paper_features_list <- list(
"Aldh1l11" = "Aldh1l1",
"Slc1a21" = "Slc1a2",
"S100b1" = "S100b",
"Gfap1" = "Gfap",
"Myoc1" = "Myoc",
"paper_cluster_0.module1" = "0 (Paper)",
"paper_cluster_1.module1" = "1 (Paper)",
"paper_cluster_2.module1" = "2 (Paper)",
"paper_cluster_3.module1" = "3 (Paper)",
"paper_cluster_4.module1" = "4 (Paper)",
"paper_cluster_5.module1" = "5 (Paper)",
"paper_cluster_6.module1" = "6 (Paper)",
"paper_cluster_7.module1" = "7 (Paper)",
"paper_cluster_8.module1" = "8 (Paper)",
"paper_cluster_9.module1" = "9 (Paper)"
)
paper_legend_only <- SpatialPlot(GSM5026144_S7_Seurat, features = "paper_cluster_9.module1", pt.size.factor = 0, image.alpha = 0) + theme(legend.text = element_blank(), legend.title = element_blank())
# Add the legend_only plot to the the plots list
paper_plots[[length(paper_plots) + 1]] <- paper_legend_only
# Create an image-only ggplot
paper_image_only <-SpatialPlot(GSM5026144_S7_Seurat, features = "paper_cluster_9.module1", pt.size.factor = 0) + theme(legend.position = "none") + annotate("text", x = Inf, y = Inf, label = "H&E", hjust = 1.1, vjust = 1.1, fontface = "italic", color = "black", size = 2)
# Add the image-only plot to the the plots list
paper_plots[[length(paper_plots) + 1]] <- paper_image_only
# Generate the plots
for (feature in names(paper_features_list)) {
# Create a spatial feature plot for each feature
p <-
SpatialFeaturePlot(GSM5026144_S7_Seurat, features = feature) + theme(legend.position = "none") + annotate(
"text",
x = Inf,
y = Inf,
label = paper_features_list[[feature]],
hjust = 1.1,
vjust = 1.1,
fontface = "italic",
color = "black",
size = 2
)
# Store the plot in the list
paper_plots[[feature]] <- p
}
With our analysis markers at 0.35 LogFC Cutoff
# Subset data for avg_logFC > 0.35
astro.cluster_markers.subset.0.35 <- astro.markers[astro.markers$avg_log2FC > 0.35, ]
# Remove ensemble ID from gene names
astro.cluster_markers.subset.0.35$gene <- sub(".*\\.", "", astro.cluster_markers.subset.0.35$gene)
# Specify the list of clusters
specified_clusters <- c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9)
# Initialize an empty list to store the cluster variables
cluster_vars <- list()
# Loop through each specified cluster
for (cluster in specified_clusters) {
# Subset data for the current cluster
cluster_data <- astro.cluster_markers.subset.0.35[astro.cluster_markers.subset.0.35$cluster == cluster, ]
# Create a comma-separated gene list with each gene in quotes
gene_list <- sprintf('"%s"', cluster_data$gene)
gene_list_string <- paste(gene_list, collapse = ", ")
# Create a list with the gene names and assign it to a named variable
variable_name <- paste0("vtp_cluster_", cluster, ".module")
variable_value <- list(eval(parse(text = paste("c(", gene_list_string, ")"))))
assign(variable_name, variable_value, envir = .GlobalEnv)
# Add the variable name to the list of cluster variables
cluster_vars <- append(cluster_vars, variable_name)
}
cluster_vars
## [[1]]
## [1] "vtp_cluster_0.module"
##
## [[2]]
## [1] "vtp_cluster_1.module"
##
## [[3]]
## [1] "vtp_cluster_2.module"
##
## [[4]]
## [1] "vtp_cluster_3.module"
##
## [[5]]
## [1] "vtp_cluster_4.module"
##
## [[6]]
## [1] "vtp_cluster_5.module"
##
## [[7]]
## [1] "vtp_cluster_6.module"
##
## [[8]]
## [1] "vtp_cluster_7.module"
##
## [[9]]
## [1] "vtp_cluster_8.module"
##
## [[10]]
## [1] "vtp_cluster_9.module"
# Loop through each cluster variable
for (cluster_var in cluster_vars) {
# Extract the cluster number from the variable name
cluster_number <- gsub("\\D", "", cluster_var)
# Create the name for the module score
module_score_name <- paste("vtp_cluster_", cluster_number, ".module", sep = "")
# Add the module score
GSM5026144_S7_Seurat <- AddModuleScore(GSM5026144_S7_Seurat, features = get(cluster_var), name = module_score_name)
}
## Warning: The following features are not present in the object: Fam212b, not
## searching for symbol synonyms
## Warning: The following features are not present in the object: 1110008F13Rik,
## Usmg5, 2410015M20Rik, Pyurf, 2010107E04Rik, 1500011K16Rik, Minos1, Fam213a,
## Tmem5, Aes, Pla2g16, D10Jhu81e, 1110004E09Rik, Fam96a, Lao1, 5930412G12Rik, 1,
## not searching for symbol synonyms
## Warning: The following features are not present in the object: Gm26772,
## Gm26870, Cd27, Gm10600, C130073E24Rik, Gm11713, Gm26909, 1810041L15Rik,
## B430010I23Rik, A330069E16Rik, not searching for symbol synonyms
## Warning: The following features are not present in the object: Fam213a, Pyurf,
## Pla2g16, 1, 1110004E09Rik, 5930412G12Rik, not searching for symbol synonyms
## Warning: The following features are not present in the object: Cxcl9, Mx1,
## Gm4951, 2410015M20Rik, Tmem5, Fam213a, F830016B08Rik, Pla2g16, 1500015O10Rik,
## 1, Fam19a5, Gm8186, 2900011O08Rik, Slc10a6, not searching for symbol synonyms
## Warning: The following features are not present in the object: Slc10a6,
## 1110004E09Rik, Cxcl9, Gm10600, 1, 1500015O10Rik, 2900011O08Rik, Cyr61, Gm8186,
## Fam19a5, not searching for symbol synonyms
## Warning: The following features are not present in the object: 1500015O10Rik,
## C130073E24Rik, Mx1, Gm26772, Gm26909, Cxcl9, Fam19a5, 1, Fam212b,
## 1500011K16Rik, 2900011O08Rik, 2410015M20Rik, Gm8186, Aes, Pyurf, Minos1, Usmg5,
## Pla2g16, Gm26870, D10Jhu81e, 2010107E04Rik, 1110004E09Rik, not searching for
## symbol synonyms
## Warning: The following features are not present in the object: Fam19a1,
## 1810041L15Rik, Cd27, not searching for symbol synonyms
## Warning: The following features are not present in the object: Gm4951,
## F830016B08Rik, Slc10a6, Cxcl9, Mx1, Slco1b2, not searching for symbol synonyms
## Warning: The following features are not present in the object: 2010107E04Rik,
## Fam213a, Usmg5, 2410015M20Rik, Minos1, Cxcl9, 1500011K16Rik, 1500015O10Rik,
## Aes, Gm8186, 2900011O08Rik, 1110008F13Rik, Lao1, Pla2g16, not searching for
## symbol synonyms
# List to store individual plots
vtp_plots <- list()
vtp_features_list <- list(
"Aldh1l11" = "Aldh1l1",
"Slc1a21" = "Slc1a2",
"S100b1" = "S100b",
"Gfap1" = "Gfap",
"Myoc1" = "Myoc",
"vtp_cluster_0.module1" = "0 (VTP)",
"vtp_cluster_1.module1" = "1 (VTP)",
"vtp_cluster_2.module1" = "2 (VTP)",
"vtp_cluster_3.module1" = "3 (VTP)",
"vtp_cluster_4.module1" = "4 (VTP)",
"vtp_cluster_5.module1" = "5 (VTP)",
"vtp_cluster_6.module1" = "6 (VTP)",
"vtp_cluster_7.module1" = "7 (VTP)",
"vtp_cluster_8.module1" = "8 (VTP)",
"vtp_cluster_9.module1" = "9 (VTP)"
)
# remove this
vtp_legend_only <- SpatialPlot(GSM5026144_S7_Seurat, features = "vtp_cluster_0.module1", pt.size.factor = 0, image.alpha = 0) + theme(legend.text = element_blank(), legend.title = element_blank())
# Add the legend_only plot to the the plots list
vtp_plots[[length(vtp_plots) + 1]] <- vtp_legend_only
# Create an image-only ggplot
vtp_image_only <- SpatialPlot(GSM5026144_S7_Seurat, features = "vtp_cluster_0.module1", pt.size.factor = 0) + theme(legend.position = "none")
# Add the image-only plot to the the plots list
vtp_plots[[length(vtp_plots) + 1]] <- vtp_image_only
# Generate the plots
for (feature in names(vtp_features_list)) {
# Create a spatial feature plot for each feature
p <-
SpatialFeaturePlot(GSM5026144_S7_Seurat, features = feature) + theme(legend.position = "none") + annotate(
"text",
x = Inf,
y = Inf,
label = vtp_features_list[[feature]],
hjust = 1.1,
vjust = 1.1,
fontface = "italic",
color = "black",
size = 2
)
# Store the plot in the list
vtp_plots[[feature]] <- p
}
LPSsuper.cluster <- 4
LPSsuper.cluster.order <- which(clusters == LPSsuper.cluster)
cluster.cols.highlight <- rep("gray70", length(clusters))
cluster.cols.highlight[LPSsuper.cluster.order] <- cluster.cols[LPSsuper.cluster.order]
names(cluster.cols.highlight) <- clusters
clusters
## [1] 0 1 2 3 4 5 6 7 8 9
DimPlot(so.astro, reduction = "tsne", group.by = "ident", split.by = "treatment") + scale_color_manual(values = c(cluster.cols.highlight))
Look up gene names with IDs
Fig4b.highlight.gene.names <- c("Gap43", "Timp1")
Fig4b.highlight.gene <- astro.markers %>% filter(gene.name %in% Fig4b.highlight.gene.names & cluster == LPSsuper.cluster) %>% print() %>% pull(gene)
## p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
## 1 0 2.839716 0.794 0.598 0 4 ENSMUSG00000047261.Gap43
## 2 0 2.667460 0.920 0.822 0 4 ENSMUSG00000001131.Timp1
## gene.name
## 1 Gap43
## 2 Timp1
so.astro %>% subset(treatment %in% "LPS") %>% FeaturePlot(., features = Fig4b.highlight.gene, min.cutoff = "q1")
Filter 30 top markers for the cluster of interest.
Fig4c.features <- astro.markers %>% filter(cluster == LPSsuper.cluster) %>% slice_max(n = 30, order_by = avg_log2FC) %>% print() %>% pull(gene)
## p_val avg_log2FC pct.1 pct.2 p_val_adj cluster
## 1 0.000000e+00 3.404696 0.704 0.589 0.000000e+00 4
## 2 0.000000e+00 3.079808 0.787 0.743 0.000000e+00 4
## 3 8.706374e-42 3.025500 0.474 0.379 1.741275e-38 4
## 4 0.000000e+00 2.896345 0.771 0.540 0.000000e+00 4
## 5 0.000000e+00 2.866512 0.710 0.680 0.000000e+00 4
## 6 3.110322e-109 2.840820 0.276 0.520 6.220645e-106 4
## 7 0.000000e+00 2.839716 0.794 0.598 0.000000e+00 4
## 8 0.000000e+00 2.760738 0.723 0.546 0.000000e+00 4
## 9 0.000000e+00 2.757122 0.759 0.448 0.000000e+00 4
## 10 0.000000e+00 2.705860 0.771 0.745 0.000000e+00 4
## 11 2.266662e-04 2.682382 0.401 0.307 4.533325e-01 4
## 12 0.000000e+00 2.667460 0.920 0.822 0.000000e+00 4
## 13 5.890018e-93 2.636812 0.565 0.479 1.178004e-89 4
## 14 1.711923e-63 2.608088 0.503 0.302 3.423847e-60 4
## 15 0.000000e+00 2.593842 0.835 0.784 0.000000e+00 4
## 16 0.000000e+00 2.576618 0.886 0.765 0.000000e+00 4
## 17 0.000000e+00 2.554800 0.786 0.641 0.000000e+00 4
## 18 0.000000e+00 2.480827 0.769 0.677 0.000000e+00 4
## 19 0.000000e+00 2.478133 0.637 0.252 0.000000e+00 4
## 20 0.000000e+00 2.467695 0.738 0.706 0.000000e+00 4
## 21 0.000000e+00 2.443389 0.695 0.550 0.000000e+00 4
## 22 0.000000e+00 2.435549 0.934 0.846 0.000000e+00 4
## 23 2.879040e-286 2.423453 0.677 0.360 5.758081e-283 4
## 24 0.000000e+00 2.407119 0.896 0.786 0.000000e+00 4
## 25 0.000000e+00 2.358614 0.702 0.622 0.000000e+00 4
## 26 0.000000e+00 2.304948 0.929 0.863 0.000000e+00 4
## 27 0.000000e+00 2.247477 0.962 0.854 0.000000e+00 4
## 28 0.000000e+00 2.206265 0.965 0.841 0.000000e+00 4
## 29 0.000000e+00 2.183335 0.813 0.628 0.000000e+00 4
## 30 1.838805e-108 2.177165 0.320 0.551 3.677610e-105 4
## gene gene.name
## 1 ENSMUSG00000029417.Cxcl9 Cxcl9
## 2 ENSMUSG00000036594.H2-Aa H2-Aa
## 3 ENSMUSG00000021469.Msx2 Msx2
## 4 ENSMUSG00000001123.Lgals9 Lgals9
## 5 ENSMUSG00000003617.Cp Cp
## 6 ENSMUSG00000021728.Emb Emb
## 7 ENSMUSG00000047261.Gap43 Gap43
## 8 ENSMUSG00000079547.H2-DMb1 H2-DMb1
## 9 ENSMUSG00000060586.H2-Eb1 H2-Eb1
## 10 ENSMUSG00000073421.H2-Ab1 H2-Ab1
## 11 ENSMUSG00000020473.Aebp1 Aebp1
## 12 ENSMUSG00000001131.Timp1 Timp1
## 13 ENSMUSG00000000753.Serpinf1 Serpinf1
## 14 ENSMUSG00000020264.Slc36a2 Slc36a2
## 15 ENSMUSG00000055172.C1ra C1ra
## 16 ENSMUSG00000067212.H2-T23 H2-T23
## 17 ENSMUSG00000028037.Ifi44 Ifi44
## 18 ENSMUSG00000027907.S100a11 S100a11
## 19 ENSMUSG00000020000.Moxd1 Moxd1
## 20 ENSMUSG00000024610.Cd74 Cd74
## 21 ENSMUSG00000021903.Galnt15 Galnt15
## 22 ENSMUSG00000024338.Psmb8 Psmb8
## 23 ENSMUSG00000067586.S1pr3 S1pr3
## 24 ENSMUSG00000046718.Bst2 Bst2
## 25 ENSMUSG00000027800.Tm4sf1 Tm4sf1
## 26 ENSMUSG00000026822.Lcn2 Lcn2
## 27 ENSMUSG00000073411.H2-D1 H2-D1
## 28 ENSMUSG00000025492.Ifitm3 Ifitm3
## 29 ENSMUSG00000022587.Ly6e Ly6e
## 30 ENSMUSG00000029304.Spp1 Spp1
# create a factor for y axis
so.astro$cluster_treatment <- factor(paste(Idents(so.astro), so.astro$treatment, sep = "_"), levels = paste0(rep(clusters, each = 2), "_", levels(so.astro$treatment)))
DotPlot(so.astro, features = Fig4c.features, assay = "RNA", group.by = "cluster_treatment") + theme(axis.text.x = element_text(angle = 90))
# Remove ensemble ID from gene names
Fig4c.features.mgi <- sub(".*\\.", "", Fig4c.features)
Fig4_cluster4_module <- list(Fig4c.features.mgi)
#Saline
GSM5026144_S7_Seurat <- AddModuleScore(GSM5026144_S7_Seurat, features = Fig4_cluster4_module, name = 'Fig4_cluster4_moduleScore')
## Warning: The following features are not present in the object: Cxcl9, not
## searching for symbol synonyms
#LPS
GSM5026148_S2_Seurat <- AddModuleScore(GSM5026148_S2_Seurat, features = Fig4_cluster4_module, name = 'Fig4_cluster4_moduleScore')
Fig4_SpatialPlot_Saline <- SpatialFeaturePlot(GSM5026144_S7_Seurat, features = 'Fig4_cluster4_moduleScore1') + theme(legend.title = element_blank()) + annotate("text", x = Inf, y = Inf, label = "Saline", hjust = 1.1, vjust = 1.1, fontface = "italic", color = "black", size = 2)
Fig4_SpatialPlot_LPS <- SpatialFeaturePlot(GSM5026148_S2_Seurat, features = 'Fig4_cluster4_moduleScore1') + theme(legend.title = element_blank()) + annotate("text", x = Inf, y = Inf, label = "LPS", hjust = 1.1, vjust = 1.1, fontface = "italic", color = "black", size = 2)
Fig_4_saline_vln <- VlnPlot(GSM5026144_S7_Seurat, features = 'Fig4_cluster4_moduleScore1')
Fig_4_lps_vln <- VlnPlot(GSM5026148_S2_Seurat, features = 'Fig4_cluster4_moduleScore1')
plot_grid(Fig_4_saline_vln, Fig_4_lps_vln, ncol = 2)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.243015 -0.079613 -0.040761 -0.034651 0.001811 0.616929
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.42083 -0.10876 -0.04137 -0.01693 0.03467 1.07209
First, construct a data frame with average expression values for cluster and treatment
# Expression data for
gene_of_interest <- "Timp1"
ID.gene_of_interest <- astro.markers %>% filter(gene.name == gene_of_interest, cluster == LPSsuper.cluster) %>% pull(gene)
expr.mat <- so.astro@assays$RNA@data
expressions.goi <- data.frame( expression = expr.mat[ID.gene_of_interest, ], cluster = Idents(so.astro), treatment = so.astro$treatment, sample = so.astro$sample_id)
expressions.goi.mean.sem <- expressions.goi %>% group_by(cluster, treatment) %>% summarise(expression.avg = mean(expression), expression.sd = sd(expression), sample.size = n()) %>% mutate(expression.sem = expression.sd / sqrt(sample.size)) %>% print()
## # A tibble: 20 × 6
## # Groups: cluster [10]
## cluster treatment expression.avg expression.sd sample.size expression.sem
## <fct> <fct> <dbl> <dbl> <int> <dbl>
## 1 0 saline 0.00151 0.0506 9366 0.000523
## 2 0 LPS 0.199 0.581 12079 0.00529
## 3 1 saline 0.00300 0.0591 7731 0.000672
## 4 1 LPS 0.299 0.622 13299 0.00539
## 5 2 saline 0.000242 0.0193 6414 0.000242
## 6 2 LPS 0.0993 0.463 14358 0.00386
## 7 3 saline 0.00225 0.0474 1699 0.00115
## 8 3 LPS 0.203 0.550 3171 0.00977
## 9 4 saline 0.00865 0.112 1485 0.00290
## 10 4 LPS 1.07 1.04 3253 0.0182
## 11 5 saline 0.00293 0.0543 1366 0.00147
## 12 5 LPS 0.220 0.616 2122 0.0134
## 13 6 saline 0 0 286 0
## 14 6 LPS 0.0912 0.487 2616 0.00952
## 15 7 saline 0.00171 0.0427 623 0.00171
## 16 7 LPS 0.243 0.657 1646 0.0162
## 17 8 saline 0 0 64 0
## 18 8 LPS 0.333 0.706 1474 0.0184
## 19 9 saline 0 0 154 0
## 20 9 LPS 0.323 0.712 344 0.0384
and the dataframe with average expression values for cluster, treatment and sample for individual dots.
expressions.goi.by.sample <- expressions.goi %>% group_by(sample, cluster, treatment) %>% summarise(expression.avg = mean(expression)) %>% print()
## # A tibble: 120 × 4
## # Groups: sample, cluster [120]
## sample cluster treatment expression.avg
## <chr> <fct> <fct> <dbl>
## 1 Female_LPS_1 0 LPS 0.304
## 2 Female_LPS_1 1 LPS 0.430
## 3 Female_LPS_1 2 LPS 0.150
## 4 Female_LPS_1 3 LPS 0.416
## 5 Female_LPS_1 4 LPS 1.28
## 6 Female_LPS_1 5 LPS 0.428
## 7 Female_LPS_1 6 LPS 0.165
## 8 Female_LPS_1 7 LPS 0.479
## 9 Female_LPS_1 8 LPS 0.414
## 10 Female_LPS_1 9 LPS 0.357
## # ℹ 110 more rows
ggplot() + geom_col(
data = expressions.goi.mean.sem,
aes(x = cluster, y = expression.avg, fill = treatment),
position = position_dodge2(width = 1.1)
) + geom_errorbar(
data = expressions.goi.mean.sem,
aes(
x = cluster,
y = expression.avg,
ymin = expression.avg - expression.sem,
ymax = expression.avg + expression.sem
),
position = position_dodge2(width = 1.1, padding = 1)
) + geom_point(
data = expressions.goi.by.sample,
aes(
x = cluster,
y = expression.avg,
fill = treatment,
group = treatment
),
position = position_jitterdodge(jitter.width = 0.1)
) + labs(fill = NULL, title = gene_of_interest) + scale_fill_manual(values = c("gray70", cluster.cols[[LPSsuper.cluster.order]]))
## R version 4.3.2 (2023-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 grid stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] imager_0.45.2 magrittr_2.0.3
## [3] ggridges_0.5.4 viridis_0.6.4
## [5] viridisLite_0.4.2 circlize_0.4.15
## [7] ComplexHeatmap_2.18.0 Seurat_5.0.1
## [9] SeuratObject_5.0.1 sp_2.1-1
## [11] UpSetR_1.4.0 SingleCellExperiment_1.24.0
## [13] SummarizedExperiment_1.32.0 Biobase_2.62.0
## [15] GenomicRanges_1.54.1 GenomeInfoDb_1.38.1
## [17] IRanges_2.36.0 S4Vectors_0.40.2
## [19] BiocGenerics_0.48.1 MatrixGenerics_1.14.0
## [21] matrixStats_1.1.0 muscat_1.16.0
## [23] gridExtra_2.3 lubridate_1.9.3
## [25] forcats_1.0.0 stringr_1.5.1
## [27] dplyr_1.1.4 purrr_1.0.2
## [29] readr_2.1.4 tidyr_1.3.0
## [31] tibble_3.2.1 ggplot2_3.4.4
## [33] tidyverse_2.0.0 cowplot_1.1.1
##
## loaded via a namespace (and not attached):
## [1] spatstat.sparse_3.0-3 bitops_1.0-7
## [3] httr_1.4.7 RColorBrewer_1.1-3
## [5] doParallel_1.0.17 numDeriv_2016.8-1.1
## [7] tools_4.3.2 sctransform_0.4.1
## [9] backports_1.4.1 utf8_1.2.4
## [11] R6_2.5.1 uwot_0.1.16
## [13] lazyeval_0.2.2 mgcv_1.9-0
## [15] GetoptLong_1.0.5 readbitmap_0.1.5
## [17] withr_2.5.2 prettyunits_1.2.0
## [19] progressr_0.14.0 cli_3.6.1
## [21] Cairo_1.6-1 spatstat.explore_3.2-5
## [23] fastDummies_1.7.3 labeling_0.4.3
## [25] sass_0.4.7 mvtnorm_1.2-3
## [27] spatstat.data_3.0-3 blme_1.0-5
## [29] pbapply_1.7-2 scater_1.30.1
## [31] parallelly_1.36.0 limma_3.58.1
## [33] rstudioapi_0.15.0 generics_0.1.3
## [35] shape_1.4.6 spatstat.random_3.2-1
## [37] gtools_3.9.5 ica_1.0-3
## [39] Matrix_1.6-3 ggbeeswarm_0.7.2
## [41] fansi_1.0.5 abind_1.4-5
## [43] lifecycle_1.0.4 yaml_2.3.7
## [45] edgeR_4.0.2 glmGamPoi_1.14.0
## [47] gplots_3.1.3 SparseArray_1.2.2
## [49] Rtsne_0.16 promises_1.2.1
## [51] crayon_1.5.2 miniUI_0.1.1.1
## [53] lattice_0.21-9 beachmat_2.18.0
## [55] pillar_1.9.0 knitr_1.45
## [57] rjson_0.2.21 boot_1.3-28.1
## [59] corpcor_1.6.10 future.apply_1.11.0
## [61] codetools_0.2-19 leiden_0.4.3.1
## [63] glue_1.6.2 data.table_1.14.8
## [65] vctrs_0.6.4 png_0.1-8
## [67] spam_2.10-0 Rdpack_2.6
## [69] gtable_0.3.4 cachem_1.0.8
## [71] xfun_0.41 rbibutils_2.2.16
## [73] S4Arrays_1.2.0 mime_0.12
## [75] survival_3.5-7 iterators_1.0.14
## [77] statmod_1.5.0 ellipsis_0.3.2
## [79] fitdistrplus_1.1-11 ROCR_1.0-11
## [81] nlme_3.1-163 pbkrtest_0.5.2
## [83] bit64_4.0.5 progress_1.2.2
## [85] EnvStats_2.8.1 RcppAnnoy_0.0.21
## [87] bslib_0.6.0 TMB_1.9.7
## [89] irlba_2.3.5.1 vipor_0.4.5
## [91] KernSmooth_2.23-22 colorspace_2.1-0
## [93] ggrastr_1.0.2 DESeq2_1.42.0
## [95] tidyselect_1.2.0 bit_4.0.5
## [97] compiler_4.3.2 BiocNeighbors_1.20.0
## [99] hdf5r_1.3.8 DelayedArray_0.28.0
## [101] plotly_4.10.3 scales_1.2.1
## [103] caTools_1.18.2 remaCor_0.0.16
## [105] lmtest_0.9-40 tiff_0.1-11
## [107] goftest_1.2-3 digest_0.6.33
## [109] spatstat.utils_3.0-4 minqa_1.2.6
## [111] variancePartition_1.32.2 rmarkdown_2.25
## [113] aod_1.3.2 XVector_0.42.0
## [115] RhpcBLASctl_0.23-42 jpeg_0.1-10
## [117] htmltools_0.5.7 pkgconfig_2.0.3
## [119] lme4_1.1-35.1 sparseMatrixStats_1.14.0
## [121] highr_0.10 fastmap_1.1.1
## [123] rlang_1.1.2 GlobalOptions_0.1.2
## [125] htmlwidgets_1.6.3 shiny_1.8.0
## [127] DelayedMatrixStats_1.24.0 farver_2.1.1
## [129] jquerylib_0.1.4 zoo_1.8-12
## [131] jsonlite_1.8.7 BiocParallel_1.36.0
## [133] BiocSingular_1.18.0 RCurl_1.98-1.13
## [135] scuttle_1.12.0 GenomeInfoDbData_1.2.11
## [137] dotCall64_1.1-0 patchwork_1.1.3
## [139] munsell_0.5.0 Rcpp_1.0.11
## [141] reticulate_1.34.0 stringi_1.8.1
## [143] zlibbioc_1.48.0 MASS_7.3-60
## [145] plyr_1.8.9 parallel_4.3.2
## [147] listenv_0.9.0 ggrepel_0.9.4
## [149] deldir_1.0-9 splines_4.3.2
## [151] tensor_1.5 hms_1.1.3
## [153] locfit_1.5-9.8 igraph_1.5.1
## [155] spatstat.geom_3.2-7 RcppHNSW_0.5.0
## [157] reshape2_1.4.4 ScaledMatrix_1.10.0
## [159] evaluate_0.23 bmp_0.3
## [161] nloptr_2.0.3 tzdb_0.4.0
## [163] foreach_1.5.2 httpuv_1.6.12
## [165] polyclip_1.10-6 RANN_2.6.1
## [167] future_1.33.0 clue_0.3-65
## [169] scattermore_1.2 rsvd_1.0.5
## [171] broom_1.0.5 xtable_1.8-4
## [173] fANCOVA_0.6-1 RSpectra_0.16-1
## [175] later_1.3.1 lmerTest_3.1-3
## [177] glmmTMB_1.1.8 beeswarm_0.4.0
## [179] cluster_2.1.4 timechange_0.2.0
## [181] globals_0.16.2
## R version 4.3.2 (2023-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 grid stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] imager_0.45.2 magrittr_2.0.3
## [3] ggridges_0.5.4 viridis_0.6.4
## [5] viridisLite_0.4.2 circlize_0.4.15
## [7] ComplexHeatmap_2.18.0 Seurat_5.0.1
## [9] SeuratObject_5.0.1 sp_2.1-1
## [11] UpSetR_1.4.0 SingleCellExperiment_1.24.0
## [13] SummarizedExperiment_1.32.0 Biobase_2.62.0
## [15] GenomicRanges_1.54.1 GenomeInfoDb_1.38.1
## [17] IRanges_2.36.0 S4Vectors_0.40.2
## [19] BiocGenerics_0.48.1 MatrixGenerics_1.14.0
## [21] matrixStats_1.1.0 muscat_1.16.0
## [23] gridExtra_2.3 lubridate_1.9.3
## [25] forcats_1.0.0 stringr_1.5.1
## [27] dplyr_1.1.4 purrr_1.0.2
## [29] readr_2.1.4 tidyr_1.3.0
## [31] tibble_3.2.1 ggplot2_3.4.4
## [33] tidyverse_2.0.0 cowplot_1.1.1
##
## loaded via a namespace (and not attached):
## [1] spatstat.sparse_3.0-3 bitops_1.0-7
## [3] httr_1.4.7 RColorBrewer_1.1-3
## [5] doParallel_1.0.17 numDeriv_2016.8-1.1
## [7] tools_4.3.2 sctransform_0.4.1
## [9] backports_1.4.1 utf8_1.2.4
## [11] R6_2.5.1 uwot_0.1.16
## [13] lazyeval_0.2.2 mgcv_1.9-0
## [15] GetoptLong_1.0.5 readbitmap_0.1.5
## [17] withr_2.5.2 prettyunits_1.2.0
## [19] progressr_0.14.0 cli_3.6.1
## [21] Cairo_1.6-1 spatstat.explore_3.2-5
## [23] fastDummies_1.7.3 labeling_0.4.3
## [25] sass_0.4.7 mvtnorm_1.2-3
## [27] spatstat.data_3.0-3 blme_1.0-5
## [29] pbapply_1.7-2 scater_1.30.1
## [31] parallelly_1.36.0 limma_3.58.1
## [33] rstudioapi_0.15.0 generics_0.1.3
## [35] shape_1.4.6 spatstat.random_3.2-1
## [37] gtools_3.9.5 ica_1.0-3
## [39] Matrix_1.6-3 ggbeeswarm_0.7.2
## [41] fansi_1.0.5 abind_1.4-5
## [43] lifecycle_1.0.4 yaml_2.3.7
## [45] edgeR_4.0.2 glmGamPoi_1.14.0
## [47] gplots_3.1.3 SparseArray_1.2.2
## [49] Rtsne_0.16 promises_1.2.1
## [51] crayon_1.5.2 miniUI_0.1.1.1
## [53] lattice_0.21-9 beachmat_2.18.0
## [55] pillar_1.9.0 knitr_1.45
## [57] rjson_0.2.21 boot_1.3-28.1
## [59] corpcor_1.6.10 future.apply_1.11.0
## [61] codetools_0.2-19 leiden_0.4.3.1
## [63] glue_1.6.2 data.table_1.14.8
## [65] vctrs_0.6.4 png_0.1-8
## [67] spam_2.10-0 Rdpack_2.6
## [69] gtable_0.3.4 cachem_1.0.8
## [71] xfun_0.41 rbibutils_2.2.16
## [73] S4Arrays_1.2.0 mime_0.12
## [75] survival_3.5-7 iterators_1.0.14
## [77] statmod_1.5.0 ellipsis_0.3.2
## [79] fitdistrplus_1.1-11 ROCR_1.0-11
## [81] nlme_3.1-163 pbkrtest_0.5.2
## [83] bit64_4.0.5 progress_1.2.2
## [85] EnvStats_2.8.1 RcppAnnoy_0.0.21
## [87] bslib_0.6.0 TMB_1.9.7
## [89] irlba_2.3.5.1 vipor_0.4.5
## [91] KernSmooth_2.23-22 colorspace_2.1-0
## [93] ggrastr_1.0.2 DESeq2_1.42.0
## [95] tidyselect_1.2.0 bit_4.0.5
## [97] compiler_4.3.2 BiocNeighbors_1.20.0
## [99] hdf5r_1.3.8 DelayedArray_0.28.0
## [101] plotly_4.10.3 scales_1.2.1
## [103] caTools_1.18.2 remaCor_0.0.16
## [105] lmtest_0.9-40 tiff_0.1-11
## [107] goftest_1.2-3 digest_0.6.33
## [109] spatstat.utils_3.0-4 minqa_1.2.6
## [111] variancePartition_1.32.2 rmarkdown_2.25
## [113] aod_1.3.2 XVector_0.42.0
## [115] RhpcBLASctl_0.23-42 jpeg_0.1-10
## [117] htmltools_0.5.7 pkgconfig_2.0.3
## [119] lme4_1.1-35.1 sparseMatrixStats_1.14.0
## [121] highr_0.10 fastmap_1.1.1
## [123] rlang_1.1.2 GlobalOptions_0.1.2
## [125] htmlwidgets_1.6.3 shiny_1.8.0
## [127] DelayedMatrixStats_1.24.0 farver_2.1.1
## [129] jquerylib_0.1.4 zoo_1.8-12
## [131] jsonlite_1.8.7 BiocParallel_1.36.0
## [133] BiocSingular_1.18.0 RCurl_1.98-1.13
## [135] scuttle_1.12.0 GenomeInfoDbData_1.2.11
## [137] dotCall64_1.1-0 patchwork_1.1.3
## [139] munsell_0.5.0 Rcpp_1.0.11
## [141] reticulate_1.34.0 stringi_1.8.1
## [143] zlibbioc_1.48.0 MASS_7.3-60
## [145] plyr_1.8.9 parallel_4.3.2
## [147] listenv_0.9.0 ggrepel_0.9.4
## [149] deldir_1.0-9 splines_4.3.2
## [151] tensor_1.5 hms_1.1.3
## [153] locfit_1.5-9.8 igraph_1.5.1
## [155] spatstat.geom_3.2-7 RcppHNSW_0.5.0
## [157] reshape2_1.4.4 ScaledMatrix_1.10.0
## [159] evaluate_0.23 bmp_0.3
## [161] nloptr_2.0.3 tzdb_0.4.0
## [163] foreach_1.5.2 httpuv_1.6.12
## [165] polyclip_1.10-6 RANN_2.6.1
## [167] future_1.33.0 clue_0.3-65
## [169] scattermore_1.2 rsvd_1.0.5
## [171] broom_1.0.5 xtable_1.8-4
## [173] fANCOVA_0.6-1 RSpectra_0.16-1
## [175] later_1.3.1 lmerTest_3.1-3
## [177] glmmTMB_1.1.8 beeswarm_0.4.0
## [179] cluster_2.1.4 timechange_0.2.0
## [181] globals_0.16.2
LPSsuper.cluster <- 8
LPSsuper.cluster.order <- which(clusters == LPSsuper.cluster)
cluster.cols.highlight <- rep("gray70", length(clusters))
cluster.cols.highlight[LPSsuper.cluster.order] <- cluster.cols[LPSsuper.cluster.order]
names(cluster.cols.highlight) <- clusters
clusters
## [1] 0 1 2 3 4 5 6 7 8 9
DimPlot(so.astro, reduction = "tsne", group.by = "ident", split.by = "treatment") + scale_color_manual(values = c(cluster.cols.highlight))
Look up gene names with IDs
Fig5b.highlight.gene.names <- c("Stat1", "Igtp", "Tap1")
Fig5b.highlight.gene <- astro.markers %>% filter(gene.name %in% Fig5b.highlight.gene.names & cluster == LPSsuper.cluster) %>% print() %>% pull(gene)
## p_val avg_log2FC pct.1 pct.2 p_val_adj cluster
## 1 0.000000e+00 5.656206 0.900 0.809 0.000000e+00 8
## 2 9.659551e-228 4.266044 0.772 0.736 1.931910e-224 8
## 3 7.928586e-150 3.657293 0.698 0.602 1.585717e-146 8
## gene gene.name
## 1 ENSMUSG00000078853.Igtp Igtp
## 2 ENSMUSG00000037321.Tap1 Tap1
## 3 ENSMUSG00000026104.Stat1 Stat1
so.astro %>% subset(treatment %in% "LPS") %>% FeaturePlot(., features = Fig5b.highlight.gene, min.cutoff = "q1")
Filter 30 top markers for the cluster of interest.
Fig5c.features <- astro.markers %>% filter(cluster == LPSsuper.cluster) %>% slice_max(n = 30, order_by = avg_log2FC) %>% print() %>% pull(gene)
## p_val avg_log2FC pct.1 pct.2 p_val_adj cluster
## 1 3.053385e-08 7.428279 0.346 0.361 6.106769e-05 8
## 2 0.000000e+00 6.370964 0.872 0.826 0.000000e+00 8
## 3 8.149065e-135 6.036104 0.672 0.668 1.629813e-131 8
## 4 1.275258e-81 5.849959 0.625 0.700 2.550516e-78 8
## 5 0.000000e+00 5.656206 0.900 0.809 0.000000e+00 8
## 6 7.795422e-12 5.228698 0.491 0.449 1.559084e-08 8
## 7 0.000000e+00 5.145975 0.925 0.856 0.000000e+00 8
## 8 5.467254e-14 5.120974 0.494 0.473 1.093451e-10 8
## 9 9.273811e-53 5.112492 0.626 0.704 1.854762e-49 8
## 10 2.221687e-101 5.033238 0.641 0.713 4.443375e-98 8
## 11 1.292782e-93 4.616974 0.652 0.759 2.585564e-90 8
## 12 2.138071e-213 4.516876 0.761 0.806 4.276141e-210 8
## 13 2.677176e-104 4.513474 0.659 0.678 5.354352e-101 8
## 14 2.840381e-34 4.454622 0.570 0.559 5.680763e-31 8
## 15 9.659551e-228 4.266044 0.772 0.736 1.931910e-224 8
## 16 4.022704e-257 4.257000 0.807 0.819 8.045407e-254 8
## 17 1.676676e-111 4.195347 0.711 0.746 3.353351e-108 8
## 18 1.498099e-24 4.192008 0.552 0.615 2.996198e-21 8
## 19 2.899336e-211 4.132627 0.745 0.714 5.798672e-208 8
## 20 2.160565e-04 4.088715 0.629 0.748 4.321130e-01 8
## 21 3.054553e-10 4.078543 0.435 0.653 6.109107e-07 8
## 22 3.466921e-07 4.074707 0.613 0.709 6.933843e-04 8
## 23 1.249664e-15 4.061852 0.651 0.749 2.499327e-12 8
## 24 1.133929e-08 3.936605 0.535 0.587 2.267858e-05 8
## 25 2.441125e-04 3.866454 0.474 0.557 4.882249e-01 8
## 26 7.385745e-69 3.845176 0.654 0.738 1.477149e-65 8
## 27 8.417197e-95 3.835184 0.638 0.516 1.683439e-91 8
## 28 5.322822e-52 3.788480 0.604 0.631 1.064564e-48 8
## 29 7.675254e-43 3.708913 0.605 0.745 1.535051e-39 8
## 30 7.928586e-150 3.657293 0.698 0.602 1.585717e-146 8
## gene gene.name
## 1 ENSMUSG00000000386.Mx1 Mx1
## 2 ENSMUSG00000054072.Iigp1 Iigp1
## 3 ENSMUSG00000078920.Ifi47 Ifi47
## 4 ENSMUSG00000079363.Gbp4 Gbp4
## 5 ENSMUSG00000078853.Igtp Igtp
## 6 ENSMUSG00000016496.Cd274 Cd274
## 7 ENSMUSG00000073555.Gm4951 Gm4951
## 8 ENSMUSG00000068245.Phf11d Phf11d
## 9 ENSMUSG00000034855.Cxcl10 Cxcl10
## 10 ENSMUSG00000069874.Irgm2 Irgm2
## 11 ENSMUSG00000104713.Gbp6 Gbp6
## 12 ENSMUSG00000090942.F830016B08Rik F830016B08Rik
## 13 ENSMUSG00000034459.Ifit1 Ifit1
## 14 ENSMUSG00000020641.Rsad2 Rsad2
## 15 ENSMUSG00000037321.Tap1 Tap1
## 16 ENSMUSG00000028270.Gbp2 Gbp2
## 17 ENSMUSG00000035692.Isg15 Isg15
## 18 ENSMUSG00000025498.Irf7 Irf7
## 19 ENSMUSG00000074896.Ifit3 Ifit3
## 20 ENSMUSG00000036594.H2-Aa H2-Aa
## 21 ENSMUSG00000028037.Ifi44 Ifi44
## 22 ENSMUSG00000024610.Cd74 Cd74
## 23 ENSMUSG00000073421.H2-Ab1 H2-Ab1
## 24 ENSMUSG00000105504.Gbp5 Gbp5
## 25 ENSMUSG00000079547.H2-DMb1 H2-DMb1
## 26 ENSMUSG00000030107.Usp18 Usp18
## 27 ENSMUSG00000040253.Gbp7 Gbp7
## 28 ENSMUSG00000062488.Ifit3b Ifit3b
## 29 ENSMUSG00000029561.Oasl2 Oasl2
## 30 ENSMUSG00000026104.Stat1 Stat1
# create a factor for y axis
so.astro$cluster_treatment <- factor(paste(Idents(so.astro), so.astro$treatment, sep = "_"), levels = paste0(rep(clusters, each = 2), "_", levels(so.astro$treatment)))
DotPlot(so.astro, features = Fig5c.features, assay = "RNA", group.by = "cluster_treatment") + theme(axis.text.x = element_text(angle = 90))
# Remove ensemble ID from gene names
Fig5c.features.mgi <- sub(".*\\.", "", Fig5c.features)
Fig5_cluster8_module <- list(Fig5c.features.mgi)
#Saline
GSM5026144_S7_Seurat <- AddModuleScore(GSM5026144_S7_Seurat, features = Fig5_cluster8_module, name = 'Fig5_cluster8_moduleScore')
## Warning: The following features are not present in the object: Mx1, Gm4951,
## F830016B08Rik, not searching for symbol synonyms
#LPS
GSM5026148_S2_Seurat <- AddModuleScore(GSM5026148_S2_Seurat, features = Fig5_cluster8_module, name = 'Fig5_cluster8_moduleScore')
## Warning: The following features are not present in the object: Mx1, Gm4951,
## F830016B08Rik, not searching for symbol synonyms
Fig5_SpatialPlot_Saline <- SpatialFeaturePlot(GSM5026144_S7_Seurat, features = 'Fig5_cluster8_moduleScore1') + theme(legend.title = element_blank()) + annotate("text", x = Inf, y = Inf, label = "Saline", hjust = 1.1, vjust = 1.1, fontface = "italic", color = "black", size = 2)
Fig5_SpatialPlot_LPS <- SpatialFeaturePlot(GSM5026148_S2_Seurat, features = 'Fig5_cluster8_moduleScore1') + theme(legend.title = element_blank()) + annotate("text", x = Inf, y = Inf, label = "LPS", hjust = 1.1, vjust = 1.1, fontface = "italic", color = "black", size = 2)
Fig_5_saline_vln <- VlnPlot(GSM5026144_S7_Seurat, features = 'Fig5_cluster8_moduleScore1')
Fig_5_lps_vln <- VlnPlot(GSM5026148_S2_Seurat, features = 'Fig5_cluster8_moduleScore1')
plot_grid(Fig_5_saline_vln, Fig_5_lps_vln, ncol = 2)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.075352 -0.025200 -0.009361 -0.003294 0.011243 0.455683
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.177668 -0.058772 -0.021851 -0.008729 0.020309 1.235414
# Expression data for
gene_of_interest <- "Igtp"
ID.gene_of_interest <- astro.markers %>% filter(gene.name == gene_of_interest, cluster == LPSsuper.cluster) %>% pull(gene)
expr.mat <- so.astro@assays$RNA@data
expressions.goi <- data.frame( expression = expr.mat[ID.gene_of_interest, ], cluster = Idents(so.astro), treatment = so.astro$treatment, sample = so.astro$sample_id)
expressions.goi.mean.sem <- expressions.goi %>% group_by(cluster, treatment) %>% summarise(expression.avg = mean(expression), expression.sd = sd(expression), sample.size = n()) %>% mutate(expression.sem = expression.sd / sqrt(sample.size)) %>% print()
## # A tibble: 20 × 6
## # Groups: cluster [10]
## cluster treatment expression.avg expression.sd sample.size expression.sem
## <fct> <fct> <dbl> <dbl> <int> <dbl>
## 1 0 saline 0.0201 0.175 9366 0.00181
## 2 0 LPS 0.0824 0.374 12079 0.00341
## 3 1 saline 0.0225 0.160 7731 0.00182
## 4 1 LPS 0.0974 0.348 13299 0.00302
## 5 2 saline 0.0119 0.154 6414 0.00192
## 6 2 LPS 0.0517 0.333 14358 0.00278
## 7 3 saline 0.0191 0.154 1699 0.00374
## 8 3 LPS 0.0971 0.362 3171 0.00642
## 9 4 saline 0.0385 0.198 1485 0.00513
## 10 4 LPS 0.145 0.437 3253 0.00766
## 11 5 saline 0.0234 0.171 1366 0.00463
## 12 5 LPS 0.0734 0.343 2122 0.00746
## 13 6 saline 0.0104 0.176 286 0.0104
## 14 6 LPS 0.0506 0.369 2616 0.00721
## 15 7 saline 0.00753 0.108 623 0.00434
## 16 7 LPS 0.0553 0.315 1646 0.00775
## 17 8 saline 0.477 0.857 64 0.107
## 18 8 LPS 1.52 1.38 1474 0.0360
## 19 9 saline 0 0 154 0
## 20 9 LPS 0.0685 0.342 344 0.0185
expressions.goi.by.sample <- expressions.goi %>% group_by(sample, cluster, treatment) %>% summarise(expression.avg = mean(expression)) %>% print()
## # A tibble: 120 × 4
## # Groups: sample, cluster [120]
## sample cluster treatment expression.avg
## <chr> <fct> <fct> <dbl>
## 1 Female_LPS_1 0 LPS 0.143
## 2 Female_LPS_1 1 LPS 0.144
## 3 Female_LPS_1 2 LPS 0.0970
## 4 Female_LPS_1 3 LPS 0.172
## 5 Female_LPS_1 4 LPS 0.213
## 6 Female_LPS_1 5 LPS 0.0973
## 7 Female_LPS_1 6 LPS 0.246
## 8 Female_LPS_1 7 LPS 0.0534
## 9 Female_LPS_1 8 LPS 2.09
## 10 Female_LPS_1 9 LPS 0.142
## # ℹ 110 more rows
ggplot() + geom_col(
data = expressions.goi.mean.sem,
aes(x = cluster, y = expression.avg, fill = treatment),
position = position_dodge2(width = 1.1)
) + geom_errorbar(
data = expressions.goi.mean.sem,
aes(
x = cluster,
y = expression.avg,
ymin = expression.avg - expression.sem,
ymax = expression.avg + expression.sem
),
position = position_dodge2(width = 1.1, padding = 1)
) + geom_point(
data = expressions.goi.by.sample,
aes(
x = cluster,
y = expression.avg,
fill = treatment,
group = treatment
),
position = position_jitterdodge(jitter.width = 0.1)
) + labs(fill = NULL, title = gene_of_interest) + scale_fill_manual(values = c("gray70", cluster.cols[[LPSsuper.cluster.order]]))
## R version 4.3.2 (2023-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 grid stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] imager_0.45.2 magrittr_2.0.3
## [3] ggridges_0.5.4 viridis_0.6.4
## [5] viridisLite_0.4.2 circlize_0.4.15
## [7] ComplexHeatmap_2.18.0 Seurat_5.0.1
## [9] SeuratObject_5.0.1 sp_2.1-1
## [11] UpSetR_1.4.0 SingleCellExperiment_1.24.0
## [13] SummarizedExperiment_1.32.0 Biobase_2.62.0
## [15] GenomicRanges_1.54.1 GenomeInfoDb_1.38.1
## [17] IRanges_2.36.0 S4Vectors_0.40.2
## [19] BiocGenerics_0.48.1 MatrixGenerics_1.14.0
## [21] matrixStats_1.1.0 muscat_1.16.0
## [23] gridExtra_2.3 lubridate_1.9.3
## [25] forcats_1.0.0 stringr_1.5.1
## [27] dplyr_1.1.4 purrr_1.0.2
## [29] readr_2.1.4 tidyr_1.3.0
## [31] tibble_3.2.1 ggplot2_3.4.4
## [33] tidyverse_2.0.0 cowplot_1.1.1
##
## loaded via a namespace (and not attached):
## [1] spatstat.sparse_3.0-3 bitops_1.0-7
## [3] httr_1.4.7 RColorBrewer_1.1-3
## [5] doParallel_1.0.17 numDeriv_2016.8-1.1
## [7] tools_4.3.2 sctransform_0.4.1
## [9] backports_1.4.1 utf8_1.2.4
## [11] R6_2.5.1 uwot_0.1.16
## [13] lazyeval_0.2.2 mgcv_1.9-0
## [15] GetoptLong_1.0.5 readbitmap_0.1.5
## [17] withr_2.5.2 prettyunits_1.2.0
## [19] progressr_0.14.0 cli_3.6.1
## [21] Cairo_1.6-1 spatstat.explore_3.2-5
## [23] fastDummies_1.7.3 labeling_0.4.3
## [25] sass_0.4.7 mvtnorm_1.2-3
## [27] spatstat.data_3.0-3 blme_1.0-5
## [29] pbapply_1.7-2 scater_1.30.1
## [31] parallelly_1.36.0 limma_3.58.1
## [33] rstudioapi_0.15.0 generics_0.1.3
## [35] shape_1.4.6 spatstat.random_3.2-1
## [37] gtools_3.9.5 ica_1.0-3
## [39] Matrix_1.6-3 ggbeeswarm_0.7.2
## [41] fansi_1.0.5 abind_1.4-5
## [43] lifecycle_1.0.4 yaml_2.3.7
## [45] edgeR_4.0.2 glmGamPoi_1.14.0
## [47] gplots_3.1.3 SparseArray_1.2.2
## [49] Rtsne_0.16 promises_1.2.1
## [51] crayon_1.5.2 miniUI_0.1.1.1
## [53] lattice_0.21-9 beachmat_2.18.0
## [55] pillar_1.9.0 knitr_1.45
## [57] rjson_0.2.21 boot_1.3-28.1
## [59] corpcor_1.6.10 future.apply_1.11.0
## [61] codetools_0.2-19 leiden_0.4.3.1
## [63] glue_1.6.2 data.table_1.14.8
## [65] vctrs_0.6.4 png_0.1-8
## [67] spam_2.10-0 Rdpack_2.6
## [69] gtable_0.3.4 cachem_1.0.8
## [71] xfun_0.41 rbibutils_2.2.16
## [73] S4Arrays_1.2.0 mime_0.12
## [75] survival_3.5-7 iterators_1.0.14
## [77] statmod_1.5.0 ellipsis_0.3.2
## [79] fitdistrplus_1.1-11 ROCR_1.0-11
## [81] nlme_3.1-163 pbkrtest_0.5.2
## [83] bit64_4.0.5 progress_1.2.2
## [85] EnvStats_2.8.1 RcppAnnoy_0.0.21
## [87] bslib_0.6.0 TMB_1.9.7
## [89] irlba_2.3.5.1 vipor_0.4.5
## [91] KernSmooth_2.23-22 colorspace_2.1-0
## [93] ggrastr_1.0.2 DESeq2_1.42.0
## [95] tidyselect_1.2.0 bit_4.0.5
## [97] compiler_4.3.2 BiocNeighbors_1.20.0
## [99] hdf5r_1.3.8 DelayedArray_0.28.0
## [101] plotly_4.10.3 scales_1.2.1
## [103] caTools_1.18.2 remaCor_0.0.16
## [105] lmtest_0.9-40 tiff_0.1-11
## [107] goftest_1.2-3 digest_0.6.33
## [109] spatstat.utils_3.0-4 minqa_1.2.6
## [111] variancePartition_1.32.2 rmarkdown_2.25
## [113] aod_1.3.2 XVector_0.42.0
## [115] RhpcBLASctl_0.23-42 jpeg_0.1-10
## [117] htmltools_0.5.7 pkgconfig_2.0.3
## [119] lme4_1.1-35.1 sparseMatrixStats_1.14.0
## [121] highr_0.10 fastmap_1.1.1
## [123] rlang_1.1.2 GlobalOptions_0.1.2
## [125] htmlwidgets_1.6.3 shiny_1.8.0
## [127] DelayedMatrixStats_1.24.0 farver_2.1.1
## [129] jquerylib_0.1.4 zoo_1.8-12
## [131] jsonlite_1.8.7 BiocParallel_1.36.0
## [133] BiocSingular_1.18.0 RCurl_1.98-1.13
## [135] scuttle_1.12.0 GenomeInfoDbData_1.2.11
## [137] dotCall64_1.1-0 patchwork_1.1.3
## [139] munsell_0.5.0 Rcpp_1.0.11
## [141] reticulate_1.34.0 stringi_1.8.1
## [143] zlibbioc_1.48.0 MASS_7.3-60
## [145] plyr_1.8.9 parallel_4.3.2
## [147] listenv_0.9.0 ggrepel_0.9.4
## [149] deldir_1.0-9 splines_4.3.2
## [151] tensor_1.5 hms_1.1.3
## [153] locfit_1.5-9.8 igraph_1.5.1
## [155] spatstat.geom_3.2-7 RcppHNSW_0.5.0
## [157] reshape2_1.4.4 ScaledMatrix_1.10.0
## [159] evaluate_0.23 bmp_0.3
## [161] nloptr_2.0.3 tzdb_0.4.0
## [163] foreach_1.5.2 httpuv_1.6.12
## [165] polyclip_1.10-6 RANN_2.6.1
## [167] future_1.33.0 clue_0.3-65
## [169] scattermore_1.2 rsvd_1.0.5
## [171] broom_1.0.5 xtable_1.8-4
## [173] fANCOVA_0.6-1 RSpectra_0.16-1
## [175] later_1.3.1 lmerTest_3.1-3
## [177] glmmTMB_1.1.8 beeswarm_0.4.0
## [179] cluster_2.1.4 timechange_0.2.0
## [181] globals_0.16.2