Load data
Fig. 2b. Astrocytes saline vs LPS
Fig. 2c. Genes statistically enriched in each cluster
- Subset of genes used for plotting
Fig. 3a. Astrocytes saline vs LPS

Load data

so.astro <- readRDS(file.path(wd, objects.dir, "SO_astrocytes_clustered.rds"))

# 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

so.astro

## An object of class Seurat 
## 14491 features across 83353 samples within 2 assays 
## Active assay: integrated (2000 features)
##  1 other assay present: RNA
##  2 dimensional reductions calculated: pca, tsne

astro.markers <- readRDS(file.path(wd, objects.dir, "astro.markers.rds"))
head(astro.markers)

Fig. 2b. Astrocytes saline vs LPS

DimPlot(so.astro, reduction = "tsne", group.by = "ident", split.by = "group_id") +
  scale_color_manual(values = cluster.cols)

Fig. 2c. Genes statistically enriched in each cluster

Subset of genes used for plotting

astro.markers.subset <- astro.markers %>% 
  group_by(cluster) %>%
  slice_max(n = 10, order_by = avg_logFC)
astro.markers.subset

astro.markers.subset %>% 
  select(gene, cluster) %>% 
  group_by(gene) %>%
  slice_min(n = 1, order_by = cluster) %>% 
  ungroup() %>% 
  arrange(cluster)

# 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.data <- so.astro@[email protected][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.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), ]

 # calculated via logFC
genes.clusters.cols.lfc <- astro.markers %>%
  mutate(gene.full = gene,
         gene = word(gene.full, 2, sep = "\\.")) %>%
  filter(gene %in% rownames(cluster_means)) %>%
  select(gene, gene.full, cluster, avg_logFC) %>%
  arrange(gene, desc(avg_logFC), cluster) %>%
  group_by(gene) %>%
  slice_max(n = 1, order_by = avg_logFC) %>%
  mutate(avg_logFC = NULL) %>%
  left_join(., enframe(cluster.cols, name = "cluster", value = "colors")) %>% 
  as.data.frame()
rownames(genes.clusters.cols.lfc) <- genes.clusters.cols$gene
genes.clusters.cols.lfc <- genes.clusters.cols.lfc[rownames(cluster_means), ]
 
 # calculated via scaled expression values
genes.clusters.cols.scaled <- cluster_means %>% 
  as.data.frame() %>% 
  rownames_to_column("gene") %>% 
  pivot_longer(names_to = "cluster", cols = as.character(clusters), values_to = "expression") %>% 
  group_by(gene) %>% 
  slice_max(n = 1, order_by = expression) %>% 
  select(-expression) %>% 
  left_join(., enframe(cluster.cols, name = "cluster", value = "colors")) %>% 
  as.data.frame()
rownames(genes.clusters.cols.scaled) <- genes.clusters.cols.scaled$gene
genes.clusters.cols.scaled <- genes.clusters.cols.scaled[rownames(cluster_means), ]

row.anno.genes <-
  rowAnnotation("Cluster genes" = genes.clusters.cols$cluster,
                col = list("Cluster genes" = 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 = 10),
  left_annotation = row.anno.genes
)

GSM5026144.dgCM <- Read10X_h5("~/data/spatial/SpaceRanger_outputs/GSM5026144_S7_filtered_feature_bc_matrix.h5")

## Warning in sparseMatrix(i = indices[] + 1, p = indptr[], x = as.numeric(x =
## counts[]), : 'giveCsparse' is deprecated; setting repr="T" for you

## 'as(<dgTMatrix>, "dgCMatrix")' is deprecated.
## Use 'as(., "CsparseMatrix")' instead.
## See help("Deprecated") and help("Matrix-deprecated").

GSM5026144_S7_Seurat <- CreateSeuratObject(counts = GSM5026144.dgCM, project = "GSM5026144_S7_Spatial")

GSM5026144_S7_Seurat <- SCTransform(object = GSM5026144_S7_Seurat, return.only.var.genes = F, verbose = FALSE)

## Warning: The 'show_progress' argument is deprecated as of v0.3. Use 'verbosity'
## instead. (in sctransform::vst)

## Calculating cell attributes from input UMI matrix: log_umi

## Variance stabilizing transformation of count matrix of size 18288 by 2511

## Model formula is y ~ log_umi

## Get Negative Binomial regression parameters per gene

## Using 2000 genes, 2511 cells

## Found 90 outliers - those will be ignored in fitting/regularization step

## Second step: Get residuals using fitted parameters for 18288 genes

## Computing corrected count matrix for 18288 genes

## Calculating gene attributes

## Wall clock passed: Time difference of 1.154081 mins

astrocyte_markers <- c("Aldh1l1", "Slc1a2", "S100b", "Gfap", "Myoc")

cluster_0.module <- list(c("Gm42418", "Arhgap31", "Kcnq1ot1", "Rmst", "Cit", "Gli1", "Ttll3", "Neat1", "Gm26917"))
cluster_1.module <- list(c("Mfge8", "Igfbp2", "Btbd17"))
cluster_2.module <- list(c("Dbi", "Agt", "Nnat", "Ppia", "Tmsb4x", "Scrg1", "Mt3", "Gstm1", "Nkain4", "S100a1"))
cluster_3.module <- list(c("Thrsp", "Actb", "Chchd10", "Hint1", "Mrpl41", "Mpc1", "Camk2n1", "Ndufc2", "Nupr1", "Igfbp2"))
cluster_4.module <- list(c("Lcn2", "Ifitm3", "Timp1", "B2m", "Vim", "Gfap", "Psmb8", "Bst2"))
cluster_5.module <- list(c("Crym", "Prss56", "Slc43a3", "Mgst1", "Slc1a5", "Tuba1a", "Itm2b", "Chil1", "Ptn"))
cluster_6.module <- list(c("Thbs4", "Igfbp5", "Fxyd6", "Gsn", "Fxyd1", "Cd9", "Meg3"))
cluster_7.module <- list(c("Itih3", "Sparc", "Nkx6-2", "Etnppl", "Gria1", "Spon1", "5031439G07Rik"))
cluster_8.module <- list(c("Igtp", "Gm4951", "Ifit3", "Iigp1", "Irgm1", "Tap1", "Gbp2", "Ifit1", "Isg15", "Cxcl10"))
cluster_9.module <- list(c("Zic1", "Slc38a1", "Col1a2", "Kcnk2", "Luzp2", "Mgll", "Zic4", "Hopx"))


for (gene in astrocyte_markers) {
  GSM5026144_S7_Seurat <- AddModuleScore(object = GSM5026144_S7_Seurat, 
                                         features = list(gene), 
                                         name = gene)
}

GSM5026144_S7_Seurat <- AddModuleScore(object = GSM5026144_S7_Seurat, features = cluster_0.module, name = 'Cluster0_ModuleScore')

GSM5026144_S7_Seurat <- AddModuleScore(object = GSM5026144_S7_Seurat, features = cluster_1.module, name = "Cluster1_ModuleScore")

GSM5026144_S7_Seurat <- AddModuleScore(object = GSM5026144_S7_Seurat, features = cluster_2.module, name = "Cluster2_ModuleScore")

GSM5026144_S7_Seurat <- AddModuleScore(object = GSM5026144_S7_Seurat, features = cluster_3.module, name = "Cluster3_ModuleScore")

GSM5026144_S7_Seurat <- AddModuleScore(object = GSM5026144_S7_Seurat, features = cluster_4.module, name = "Cluster4_ModuleScore")

GSM5026144_S7_Seurat <- AddModuleScore(object = GSM5026144_S7_Seurat, features = cluster_5.module, name = "Cluster5_ModuleScore")

GSM5026144_S7_Seurat <- AddModuleScore(object = GSM5026144_S7_Seurat, features = cluster_6.module, name = "Cluster6_ModuleScore")

GSM5026144_S7_Seurat <- AddModuleScore(object = GSM5026144_S7_Seurat, features = cluster_7.module, name = "Cluster7_ModuleScore")

GSM5026144_S7_Seurat <- AddModuleScore(object = GSM5026144_S7_Seurat, features = cluster_8.module, name = "Cluster8_ModuleScore")

GSM5026144_S7_Seurat <- AddModuleScore(object = GSM5026144_S7_Seurat, features = cluster_9.module, name = "Cluster9_ModuleScore")

#module_scores <- [email protected] %>%
# dplyr::select(starts_with("Cluster"))

module_scores <- [email protected] %>%
  dplyr::select(starts_with("Cluster"), starts_with("Aldh1l1"), starts_with("Slc1a2"), 
                starts_with("S100b"), starts_with("Gfap"), starts_with("Myoc"))

# Create a new data frame with a column 'V1' for the row names
module_scores_new <- data.frame(V1 = rownames(module_scores), module_scores)

# Now check the result
head(module_scores_new)

spatial_coords <- read.csv("~/data/spatial/SpaceRanger_outputs/GSM5026144_tissue_positions_list.csv.gz", header = FALSE)



data <- dplyr::inner_join(spatial_coords, module_scores_new, by = "V1")

img <- imager::load.image("~/GSM5026144_S7_tissue_hires_image.png")

img_rotated <- imager::imrotate(img, -90)  # Rotates the image 90 degrees counterclockwise

dim(img_rotated)

## [1] 2000 1865    1    3

# Create a raster grob with the rotated image
GSM5026144_S7_he_img <- rasterGrob(img_rotated, interpolate=TRUE)

# Extract the underlying matrix or array from the rasterGrob
img_data <- GSM5026144_S7_he_img$raster

# Get the dimensions
img_dims <- dim(img_data)
img_width <- img_dims[2]  # The second element represents the width
img_height <- img_dims[1]  # The first element represents the height

# Define the scale factor
scale_factor <- 0.54555374

# Adjust the coordinates
data$V5_scaled <- data$V5 * scale_factor
data$V6_scaled <- data$V6 * scale_factor

# Create the main plot

color_func <- colorRamp2(seq(0, 1, length.out = 5), c("#6256a8", "#8acda3", "#faf6c0", "#f48f54", "#a11944"))

legend_only <- ggplot(data, aes(x = V5_scaled, y = V6_scaled, fill = Myoc1)) +
  geom_point(aes(fill = Myoc1), color = NA, stroke = 0.2, shape = 21) +
  scale_fill_gradientn(colors = color_func(seq(0, 1, length.out = 100)), 
                       breaks = c(min(data$Myoc1), max(data$Myoc1)),
                       labels = c("Lo", "Hi"),
                       guide = "colourbar") +
  labs(fill = "") +
  theme_minimal()

gt <- ggplotGrob(legend_only)

## Warning: Removed 2511 rows containing missing values (`geom_point()`).

# Extract the legend as its own gtable
legend_gtable <- gtable::gtable_filter(gt, "guide-box")

# Convert the gtable legend to a grob
legend_grob <- grid::grobTree(legend_gtable)


# List to store individual plots
plots <- list()

plots$Legend <- legend_grob

# Create an image-only ggplot
image_only <- ggplot(data, aes(x = V5_scaled, y = V6_scaled, fill = Myoc1)) +
  # Use the defined image dimensions
  annotation_custom(GSM5026144_S7_he_img, 
                    xmin = 0, xmax = img_width, 
                    ymin = 0, ymax = img_height) +
  geom_point(aes(fill = Myoc1), color = NA, stroke = 0.2, shape = 21) +
  scale_fill_gradientn(colors = color_func(seq(0, 1, length.out = 100)), 
                       guide = FALSE) +  # Removes the legend
  theme_minimal() +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.x = element_blank(),
    axis.ticks.y = element_blank(),
    legend.title = element_blank()
  ) +
  # Adding the "H&E" label to the top-left corner
  annotate("text", x = min(data$V5_scaled), y = max(data$V6_scaled), 
           label = "H&E", hjust = 0, vjust = 1, fontface = "bold", size = 5)


# Add the image-only plot to the beginning of the plots list
plots[[length(plots) + 1]] <- image_only

label_list <- list(
  "Aldh1l11" = "Aldh1l1",
  "Slc1a21" = "Slc1a2",
  "S100b1" = "S100b",
  "Gfap1" = "Gfap",
  "Myoc1" = "Myoc",
  "Cluster0_ModuleScore1" = "0",
  "Cluster1_ModuleScore1" = "1",
  "Cluster2_ModuleScore1" = "2",
  "Cluster3_ModuleScore1" = "3",
  "Cluster4_ModuleScore1" = "4",
  "Cluster5_ModuleScore1" = "5",
  "Cluster6_ModuleScore1" = "6",
  "Cluster7_ModuleScore1" = "7",
  "Cluster8_ModuleScore1" = "8",
  "Cluster9_ModuleScore1" = "9"
)

for (col in names(label_list)) {
  p <- ggplot(data, aes(x = V5_scaled, y = V6_scaled, fill = !!sym(col))) +
    annotation_custom(GSM5026144_S7_he_img, 
                      xmin = 0, xmax = img_width, 
                      ymin = 0, ymax = img_height) +
    geom_point(color = "black", stroke = 0.2, shape = 21) +
    scale_fill_gradientn(colors = color_func(seq(0, 1, length.out = 100))) +
    theme_minimal() +
    theme(
      axis.title.x = element_blank(),
      axis.title.y = element_blank(),
      axis.text.x = element_blank(),
      axis.text.y = element_blank(),
      axis.ticks.x = element_blank(),
      axis.ticks.y = element_blank(),
      legend.position = "none"  
    ) +
    # Adding the label annotation to the top-left corner
    annotate("text", x = min(data$V5_scaled), y = max(data$V6_scaled), 
             label = label_list[[col]], hjust = 0, vjust = 1, fontface = "bold", size = 5)
  
  plots[[col]] <- p
}

# Arrange plots in a grid
grid.arrange(grobs = plots, ncol = 5)

Fig. 3a. Astrocytes saline vs LPS

DimPlot(so.astro, group.by = "integrated_snn_res.0.2", split.by = "treatment") +
  scale_color_manual(values = cluster.cols)

sessionInfo()

## R version 3.6.2 (2019-12-12)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Debian GNU/Linux 10 (buster)
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.3.5.so
## 
## 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=C             
##  [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       
## 
## attached base packages:
##  [1] parallel  stats4    grid      stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] imager_0.42.1               magrittr_2.0.3             
##  [3] ggridges_0.5.4              viridis_0.5.1              
##  [5] viridisLite_0.4.2           circlize_0.4.8             
##  [7] ComplexHeatmap_2.2.0        Seurat_3.1.0               
##  [9] UpSetR_1.4.0                muscat_1.0.1               
## [11] scater_1.14.6               SingleCellExperiment_1.8.0 
## [13] SummarizedExperiment_1.16.1 DelayedArray_0.12.3        
## [15] BiocParallel_1.20.1         matrixStats_1.0.0          
## [17] Biobase_2.46.0              GenomicRanges_1.38.0       
## [19] GenomeInfoDb_1.22.1         IRanges_2.20.2             
## [21] S4Vectors_0.24.4            BiocGenerics_0.32.0        
## [23] gridExtra_2.3               forcats_1.0.0              
## [25] stringr_1.5.0               dplyr_1.1.3                
## [27] purrr_1.0.2                 readr_2.1.4                
## [29] tidyr_1.3.0                 tibble_3.2.1               
## [31] ggplot2_3.4.3               tidyverse_1.3.0            
## [33] cowplot_1.1.1              
## 
## loaded via a namespace (and not attached):
##   [1] R.methodsS3_1.8.2        acepack_1.4.2            bit64_4.0.5             
##   [4] knitr_1.43               irlba_2.3.5.1            multcomp_1.4-25         
##   [7] R.utils_2.12.2           data.table_1.14.8        rpart_4.1-15            
##  [10] RCurl_1.98-1.12          doParallel_1.0.17        generics_0.1.3          
##  [13] metap_1.8                TH.data_1.1-2            RSQLite_2.3.1           
##  [16] RANN_2.6.1               future_1.33.0            bit_4.0.5               
##  [19] tzdb_0.4.0               mutoss_0.1-13            xml2_1.3.5              
##  [22] lubridate_1.9.2          xfun_0.40                hms_1.1.3               
##  [25] jquerylib_0.1.4          evaluate_0.21            fansi_1.0.4             
##  [28] progress_1.2.2           caTools_1.18.2           dbplyr_2.3.3            
##  [31] readxl_1.4.3             igraph_1.5.1             DBI_1.1.3               
##  [34] geneplotter_1.64.0       htmlwidgets_1.6.2        backports_1.4.1         
##  [37] annotate_1.64.0          deldir_1.0-9             vctrs_0.6.3             
##  [40] ROCR_1.0-11              cachem_1.0.8             withr_2.5.0             
##  [43] checkmate_2.2.0          sctransform_0.3.5        prettyunits_1.1.1       
##  [46] mnormt_2.1.1             cluster_2.1.0            lazyeval_0.2.2          
##  [49] ape_5.7-1                crayon_1.5.2             hdf5r_1.3.0             
##  [52] genefilter_1.68.0        labeling_0.4.3           edgeR_3.28.1            
##  [55] pkgconfig_2.0.3          qqconf_1.1.1             nlme_3.1-142            
##  [58] vipor_0.4.5              blme_1.0-5               nnet_7.3-12             
##  [61] rlang_1.1.1              globals_0.16.2           lifecycle_1.0.3         
##  [64] sandwich_3.0-2           mathjaxr_1.6-0           modelr_0.1.11           
##  [67] rsvd_1.0.3               bmp_0.3                  cellranger_1.1.0        
##  [70] lmtest_0.9-40            tiff_0.1-5               Matrix_1.6-1            
##  [73] boot_1.3-23              zoo_1.8-12               reprex_2.0.2            
##  [76] base64enc_0.1-3          beeswarm_0.4.0           GlobalOptions_0.1.2     
##  [79] png_0.1-8                rjson_0.2.9              bitops_1.0-7            
##  [82] R.oo_1.25.0              KernSmooth_2.23-16       blob_1.2.4              
##  [85] DelayedMatrixStats_1.8.0 shape_1.4.6              parallelly_1.36.0       
##  [88] jpeg_0.1-10              scales_1.2.1             memoise_2.0.1           
##  [91] plyr_1.8.8               ica_1.0-3                gplots_3.1.3            
##  [94] zlibbioc_1.32.0          compiler_3.6.2           RColorBrewer_1.1-3      
##  [97] plotrix_3.8-2            clue_0.3-64              lme4_1.1-34             
## [100] DESeq2_1.26.0            fitdistrplus_1.1-11      cli_3.6.1               
## [103] XVector_0.26.0           lmerTest_3.1-3           listenv_0.9.0           
## [106] pbapply_1.7-2            TMB_1.9.6                htmlTable_2.4.1         
## [109] Formula_1.2-5            MASS_7.3-51.4            tidyselect_1.2.0        
## [112] stringi_1.7.12           highr_0.10               yaml_2.3.7              
## [115] BiocSingular_1.2.2       locfit_1.5-9.4           latticeExtra_0.6-30     
## [118] ggrepel_0.9.3            sass_0.4.7               tools_3.6.2             
## [121] timechange_0.2.0         future.apply_1.11.0      rstudioapi_0.15.0       
## [124] foreach_1.5.2            foreign_0.8-72           farver_2.1.1            
## [127] Rtsne_0.16               digest_0.6.33            Rcpp_1.0.11             
## [130] broom_1.0.5              SDMTools_1.1-221.2       RcppAnnoy_0.0.21        
## [133] readbitmap_0.1.5         httr_1.4.7               AnnotationDbi_1.48.0    
## [136] Rdpack_2.5               colorspace_2.1-0         rvest_1.0.3             
## [139] XML_3.99-0.3             fs_1.6.3                 reticulate_1.31         
## [142] splines_3.6.2            uwot_0.1.16              sn_2.1.1                
## [145] multtest_2.42.0          plotly_4.10.2            xtable_1.8-4            
## [148] jsonlite_1.8.7           nloptr_2.0.3             R6_2.5.1                
## [151] TFisher_0.2.0            Hmisc_4.4-0              pillar_1.9.0            
## [154] htmltools_0.5.6          glue_1.6.2               fastmap_1.1.1           
## [157] minqa_1.2.5              BiocNeighbors_1.4.2      codetools_0.2-16        
## [160] tsne_0.1-3.1             mvtnorm_1.2-3            utf8_1.2.3              
## [163] lattice_0.20-38          bslib_0.5.1              numDeriv_2016.8-1.1     
## [166] pbkrtest_0.4-8.6         ggbeeswarm_0.7.2         leiden_0.4.3            
## [169] colorRamps_2.3.1         gtools_3.9.4             interp_1.1-4            
## [172] survival_3.1-8           limma_3.42.2             glmmTMB_1.1.7           
## [175] rmarkdown_2.24           munsell_0.5.0            GetoptLong_1.0.5        
## [178] GenomeInfoDbData_1.2.2   iterators_1.0.14         variancePartition_1.16.1
## [181] haven_2.5.3              reshape2_1.4.4           gtable_0.3.4            
## [184] rbibutils_2.2.15

Astrocytes - cluster analysis

Load data

Fig. 2b. Astrocytes saline vs LPS

Fig. 2c. Genes statistically enriched in each cluster

Subset of genes used for plotting

Fig. 3a. Astrocytes saline vs LPS