Figures 2, 4, and 5
Setup: Libraries, Objects
These libraries and objects are used in Figures 2, 4, and 5.
Libraries
# 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)scRNA-seq Astro Object
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 9scRNA-seq Markers Object
astro.markers <- readRDS("~/data/Followup/astro.markers.v1.rds") %>%
remove_rownames() %>%
mutate(gene.name = word(gene, 2, sep = "\\."))
head(astro.markers)Visium Spatial objects
Their were six spatial transcriptomics samples in the paper, three Saline and three LPS. From Extended Data Fig. 4Qc and validation of six Visium spatial transcriptomics sections from saline- and LPS-treated animals:
Sample to GEO ID Mapping:
- Saline_1 (
GSM5026144_S7) - Saline_2 (
GSM5026146_S4) - Saline_3 (
GSM5026145_S3) - LPS_1 (
GSM5026149_S6) - LPS_2 (
GSM5026148_S2) - LPS_3 (
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)Fig. 2: Astrocyte Clusters
2b. t-SNE
DimPlot(so.astro, reduction = "tsne", group.by = "ident", split.by = "group_id") + scale_color_manual(values = cluster.cols)
2c. Cluster Markers
# 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)
2d. Localization by Cluster
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
}
Fig. 4: Cluster 4
4a. t-SNE
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 9DimPlot(so.astro, reduction = "tsne", group.by = "ident", split.by = "treatment") + scale_color_manual(values = c(cluster.cols.highlight))
4b. Gap43, Timp1 t-SNE
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 Timp1so.astro %>% subset(treatment %in% "LPS") %>% FeaturePlot(., features = Fig4b.highlight.gene, min.cutoff = "q1")
4c. Dot plot
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))
4d. Spatial Transcriptomics
# 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) 
4e. Mod. Score Distributions
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.072094f. Timp1 expr. by cluster
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.0384and 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 rowsggplot() + 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.2Fig. 5: Cluster 8
5a. t-SNE
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 9DimPlot(so.astro, reduction = "tsne", group.by = "ident", split.by = "treatment") + scale_color_manual(values = c(cluster.cols.highlight))
5b. Stat1, Igtp, Tap1 t-SNE
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 Stat1so.astro %>% subset(treatment %in% "LPS") %>% FeaturePlot(., features = Fig5b.highlight.gene, min.cutoff = "q1")
5c. Dot plot
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))
5d. Spatial transcriptomics
# 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 synonymsFig5_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) 
5e. Mod. Score Distributions
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.2354145f. Igtp expr. by cluster
# 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.0185expressions.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 rowsggplot() + 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]]))
Utilized packages & versions
## 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