1 Introduction

library(ggplot2)
library(vctrs)
library(patchwork,lib.loc = "/apps/rocs/2020.08/cascadelake/software/R/4.1.2-foss-2020a/lib64/R/library")
library(Seurat)
library(dplyr)
library(tidyr)
library(readr)
library(stringr)
library(progeny)
library(SingleCellExperiment)
library(BayesSpace)
source(file = "WrapperFunction/SeuratWrappers.R")

data_directory <- 
 params$data_directory
analysis_name <- 
  params$analysis_name

# setwd(paste0(data_directory, analysis_name))
input_folder <- "IntermediaryFiles/"

sample <- "SN048_A121573_Rep1"
list_results <- readRDS(paste0(data_directory, analysis_name, input_folder, 
      "BayesSpaceResults/Pathology_Tumor_Annotations/", sample , ".rds"))
sce_subset <- list_results[[2]]
sce_subset.enhanced <-list_results[[1]]
## There are actually only 3 clusters due to the non matched neighbors.
sce_subset$spatial.cluster[sce_subset$spatial.cluster == 4] <-1 

## We convert to a Seurat object 
seurat_obj.enhanced <- Seurat::CreateSeuratObject(
  counts=SingleCellExperiment::logcounts(sce_subset.enhanced),
  assay='Spatial',meta.data=as.data.frame(colData(sce_subset.enhanced)))
seurat_obj.enhanced <- Seurat::SetIdent(seurat_obj.enhanced, value = "spatial.cluster")

seurat_obj.enhanced@assays$Spatial@scale.data <-
    seurat_obj.enhanced@assays$Spatial@data %>% as.matrix %>% t %>% scale %>% t


seurat_obj.enhanced <- 
  progeny(seurat_obj.enhanced, assay_name = "Spatial", top = 500, return_assay = TRUE, scale = FALSE)
DefaultAssay(seurat_obj.enhanced) <- "progeny"
sce_obj <- as.SingleCellExperiment(seurat_obj.enhanced, assay = "progeny")

2 Results: Pathway activity at the subspot resolution

all_pathways <- rownames(seurat_obj.enhanced)
for ( current_pathway in all_pathways){
  
  cat( "## Pathway: ", current_pathway , "\n\n")
  
  feature_plot <- 
    featurePlot(sce_obj, feature=current_pathway, diverging = TRUE,  
             platform = "Visium",is.enhanced = TRUE, color = NA)
  
  print(feature_plot)
  
  cat( "\n\n")

  violin_plots <- Seurat::VlnPlot(object = seurat_obj.enhanced, 
    features = current_pathway ,group.by = "spatial.cluster")
  
  print(violin_plots)
  
  cat( "\n\n")
}

2.1 Pathway: Androgen

2.2 Pathway: EGFR

2.3 Pathway: Estrogen

2.4 Pathway: Hypoxia

2.5 Pathway: JAK-STAT

2.6 Pathway: MAPK

2.7 Pathway: NFkB

2.8 Pathway: p53

2.9 Pathway: PI3K

2.10 Pathway: TGFb

2.11 Pathway: TNFa

2.12 Pathway: Trail

2.13 Pathway: VEGF

2.14 Pathway: WNT

3 Plot for publication

# We first need to prepare the data to later on perform a Kruskal-Wallis Test to 
# assess the significance of the results
q <- 3
palette <- RColorBrewer::brewer.pal(q, "Dark2")

subspot_cluster_df <- seurat_obj.enhanced@meta.data %>% 
  tibble::rownames_to_column(var = 'spot_id') %>% 
  dplyr::select(spot_id, spatial.cluster)

pathway_spot_df <- GetAssayData(seurat_obj.enhanced) %>% 
  as.data.frame() %>% tibble::rownames_to_column(var = 'Pathway') %>% 
  pivot_longer(!Pathway, names_to = "spot_id", values_to = "Activity") %>% 
  dplyr::inner_join(subspot_cluster_df)
current_pathway <- "TGFb"

current_pathway_spot_df <- 
  dplyr::filter(pathway_spot_df, Pathway  == current_pathway)

results_kruskal <- 
  kruskal.test(Activity  ~ spatial.cluster, data = current_pathway_spot_df)

# if (results_kruskal$p.value < 2.2e-16){
#  pvalue_label <- "P-value < 2.2e-16"
#} else { 
#  pvalue_label <-paste0("P-value = ", results_kruskal$p.value) 
#}

pvalue_label <- paste0("P-value = ", signif(results_kruskal$p.value,digits=3)) 

feature_plot_TGFb <- 
    featurePlot(sce_obj, feature=current_pathway, diverging = TRUE,  
             platform = "Visium",is.enhanced = TRUE, color = NA) + 
  theme(legend.position = "left", legend.title = element_blank())
feature_plot_TGFb

current_vln_plot_TGFb <- VlnPlot(seurat_obj.enhanced, features=current_pathway, 
                cols = palette, adjust=1.25, pt.size =1) + 
  theme(axis.title.x = element_blank(), 
    plot.title = element_text(hjust = 0), axis.text.x = element_blank(), 
    legend.position = "bottom") + 
  guides(fill=guide_legend(nrow=1,byrow=TRUE, title="Cluster")) +
  labs(y = "Activity Level") + 
  geom_text(x=3, y=500, label= pvalue_label, 
  size = 4.5, family = "Arial",  colour = "darkblue")
current_vln_plot_TGFb

4 Supplementary Plots

current_pathway <- "EGFR"

current_pathway_spot_df <- 
  dplyr::filter(pathway_spot_df, Pathway  == current_pathway)

results_kruskal <- 
  kruskal.test(Activity  ~ spatial.cluster, data = current_pathway_spot_df)

# if (results_kruskal$p.value < 2.2e-16){
#  pvalue_label <- "P-value < 2.2e-16"
#} else { 
#  pvalue_label <-paste0("P-value = ", results_kruskal$p.value) 
#}
pvalue_label <- paste0("P-value = ", signif(results_kruskal$p.value,digits=3)) 


feature_plot_EGFR <- 
    featurePlot(sce_obj, feature=current_pathway, diverging = TRUE,  
             platform = "Visium",is.enhanced = TRUE, color = NA) + 
  theme(legend.position = "left", legend.title = element_blank())
feature_plot_EGFR

current_vln_plot_EGFR <- VlnPlot(seurat_obj.enhanced, features=current_pathway, 
                cols = palette, adjust=1.25, pt.size =1) + 
  theme(axis.title.x = element_blank(), 
    plot.title = element_text(hjust = 0), axis.text.x = element_blank(), 
    legend.position = "bottom") + 
  guides(fill=guide_legend(nrow=1,byrow=TRUE, title="Cluster")) +
  labs(y = "Activity Level") + 
  geom_text(x=2, y=100, label= pvalue_label, 
  size = 4.5, family = "Arial",  colour = "darkblue")
current_vln_plot_EGFR

current_pathway <- "MAPK"

current_pathway_spot_df <- 
  dplyr::filter(pathway_spot_df, Pathway  == current_pathway)

results_kruskal <- 
  kruskal.test(Activity  ~ spatial.cluster, data = current_pathway_spot_df)

# if (results_kruskal$p.value < 2.2e-16){
#  pvalue_label <- "P-value < 2.2e-16"
#} else { 
#  pvalue_label <-paste0("P-value = ", results_kruskal$p.value) 
#}
pvalue_label <- paste0("P-value = ", signif(results_kruskal$p.value,digits=3)) 


feature_plot_MAPK <- 
    featurePlot(sce_obj, feature=current_pathway, diverging = TRUE,  
             platform = "Visium",is.enhanced = TRUE, color = NA) + 
  theme(legend.position = "left", legend.title = element_blank())
feature_plot_MAPK

current_vln_plot_MAPK <- VlnPlot(seurat_obj.enhanced, features=current_pathway, 
                cols = palette, adjust=1.25, pt.size =1) + 
  theme(axis.title.x = element_blank(), 
    plot.title = element_text(hjust = 0), axis.text.x = element_blank(), 
    legend.position = "bottom") + 
  guides(fill=guide_legend(nrow=1,byrow=TRUE, title="Cluster")) +
  labs(y = "Activity Level") + 
  geom_text(x=2, y=250, label= pvalue_label, 
  size = 4.5, family = "Arial",  colour = "darkblue")
current_vln_plot_MAPK

current_pathway <- "WNT"

current_pathway_spot_df <- 
  dplyr::filter(pathway_spot_df, Pathway  == current_pathway)

results_kruskal <- 
  kruskal.test(Activity  ~ spatial.cluster, data = current_pathway_spot_df)

# if (results_kruskal$p.value < 2.2e-16){
#  pvalue_label <- "P-value < 2.2e-16"
#} else { 
#  pvalue_label <-paste0("P-value = ", results_kruskal$p.value) 
#}
pvalue_label <- paste0("P-value = ", signif(results_kruskal$p.value,digits=3)) 


feature_plot_WNT <- 
    featurePlot(sce_obj, feature=current_pathway, diverging = TRUE,  
             platform = "Visium",is.enhanced = TRUE, color = NA) + 
  theme(legend.position = "left", legend.title = element_blank())
feature_plot_WNT

current_vln_plot_WNT <- VlnPlot(seurat_obj.enhanced, features=current_pathway, 
                cols = palette, adjust=1.25, pt.size =1) + 
  theme(axis.title.x = element_blank(), 
    plot.title = element_text(hjust = 0), axis.text.x = element_blank(), 
    legend.position = "bottom") + 
  guides(fill=guide_legend(nrow=1,byrow=TRUE, title="Cluster")) +
  labs(y = "Activity Level") + 
  geom_text(x=1, y=30, label= pvalue_label, 
  size = 4.5, family = "Arial",  colour = "darkblue")
current_vln_plot_WNT

current_pathway <- "VEGF"

current_pathway_spot_df <- 
  dplyr::filter(pathway_spot_df, Pathway  == current_pathway)

results_kruskal <- 
  kruskal.test(Activity  ~ spatial.cluster, data = current_pathway_spot_df)

# if (results_kruskal$p.value < 2.2e-16){
#  pvalue_label <- "P-value < 2.2e-16"
#} else { 
#  pvalue_label <-paste0("P-value = ", results_kruskal$p.value) 
#}
pvalue_label <- paste0("P-value = ", signif(results_kruskal$p.value,digits=3)) 

feature_plot_VEGF <- 
    featurePlot(sce_obj, feature=current_pathway, diverging = TRUE,  
             platform = "Visium",is.enhanced = TRUE, color = NA) + 
  theme(legend.position = "left", legend.title = element_blank())
feature_plot_VEGF

current_vln_plot_VEGF <- VlnPlot(seurat_obj.enhanced, features=current_pathway, 
                cols = palette, adjust=1.25, pt.size =1) + 
  theme(axis.title.x = element_blank(), 
    plot.title = element_text(hjust = 0), axis.text.x = element_blank(), 
    legend.position = "bottom") + 
  guides(fill=guide_legend(nrow=1,byrow=TRUE, title="Cluster")) +
  labs(y = "Activity Level") + 
  geom_text(x=2, y=60, label= pvalue_label, 
  size = 4.5, family = "Arial",  colour = "darkblue")
current_vln_plot_VEGF

5 Session Info Details

## R version 4.1.2 (2021-11-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Fedora 33 (Container Image)
## 
## Matrix products: default
## BLAS/LAPACK: /apps/rocs/2020.08/cascadelake/software/OpenBLAS/0.3.9-GCC-9.3.0/lib/libopenblas_skylakexp-r0.3.9.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=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       
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] BayesSpace_1.4.1            SingleCellExperiment_1.16.0
##  [3] SummarizedExperiment_1.24.0 Biobase_2.54.0             
##  [5] GenomicRanges_1.46.1        GenomeInfoDb_1.30.1        
##  [7] IRanges_2.28.0              S4Vectors_0.32.4           
##  [9] BiocGenerics_0.40.0         MatrixGenerics_1.6.0       
## [11] matrixStats_0.62.0          progeny_1.16.0             
## [13] stringr_1.4.0               readr_2.1.2                
## [15] tidyr_1.2.0                 dplyr_1.0.9                
## [17] sp_1.5-0                    SeuratObject_4.1.0         
## [19] Seurat_4.1.0                patchwork_1.1.1            
## [21] vctrs_0.4.1                 ggplot2_3.3.6              
## [23] BiocManager_1.30.18        
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.2                reticulate_1.25          
##   [3] tidyselect_1.1.2          RSQLite_2.2.14           
##   [5] htmlwidgets_1.5.4         grid_4.1.2               
##   [7] BiocParallel_1.28.3       Rtsne_0.16               
##   [9] miscTools_0.6-27          munsell_0.5.0            
##  [11] ScaledMatrix_1.2.0        codetools_0.2-18         
##  [13] ica_1.0-2                 xgboost_1.6.0.1          
##  [15] statmod_1.4.36            scran_1.22.1             
##  [17] future_1.26.1             miniUI_0.1.1.1           
##  [19] withr_2.5.0               spatstat.random_2.2-0    
##  [21] colorspace_2.1-0          progressr_0.10.1         
##  [23] filelock_1.0.2            highr_0.9                
##  [25] knitr_1.39                rstudioapi_0.13          
##  [27] ROCR_1.0-11               tensor_1.5               
##  [29] listenv_0.8.0             labeling_0.4.2           
##  [31] GenomeInfoDbData_1.2.7    polyclip_1.10-0          
##  [33] farver_2.1.0              bit64_4.0.5              
##  [35] rhdf5_2.38.1              coda_0.19-4              
##  [37] parallelly_1.32.0         generics_0.1.2           
##  [39] xfun_0.31                 BiocFileCache_2.2.1      
##  [41] R6_2.5.1                  ggbeeswarm_0.6.0         
##  [43] rsvd_1.0.5                locfit_1.5-9.5           
##  [45] bitops_1.0-7              rhdf5filters_1.6.0       
##  [47] spatstat.utils_2.3-1      cachem_1.0.6             
##  [49] DelayedArray_0.20.0       assertthat_0.2.1         
##  [51] promises_1.2.0.1          scales_1.2.0             
##  [53] beeswarm_0.4.0            rgeos_0.5-10             
##  [55] gtable_0.3.0              beachmat_2.10.0          
##  [57] Cairo_1.5-15              globals_0.15.0           
##  [59] goftest_1.2-3             sandwich_3.0-2           
##  [61] rlang_1.0.2               splines_4.1.2            
##  [63] lazyeval_0.2.2            spatstat.geom_2.4-0      
##  [65] yaml_2.3.5                reshape2_1.4.4           
##  [67] abind_1.4-7               httpuv_1.6.5             
##  [69] tools_4.1.2               ellipsis_0.3.2           
##  [71] spatstat.core_2.4-4       jquerylib_0.1.4          
##  [73] RColorBrewer_1.1-3        ggridges_0.5.3           
##  [75] Rcpp_1.0.8.3              plyr_1.8.7               
##  [77] sparseMatrixStats_1.6.0   zlibbioc_1.40.0          
##  [79] purrr_0.3.4               RCurl_1.98-1.7           
##  [81] rpart_4.1.16              deldir_1.0-6             
##  [83] viridis_0.6.2             pbapply_1.5-0            
##  [85] cowplot_1.1.1             zoo_1.8-10               
##  [87] ggrepel_0.9.1             cluster_2.1.3            
##  [89] magrittr_2.0.3            data.table_1.14.2        
##  [91] scattermore_0.8           lmtest_0.9-40            
##  [93] RANN_2.6.1                fitdistrplus_1.1-8       
##  [95] hms_1.1.1                 mime_0.12                
##  [97] evaluate_0.15             xtable_1.8-6             
##  [99] mclust_5.4.10             gridExtra_2.3            
## [101] compiler_4.1.2            scater_1.22.0            
## [103] tibble_3.1.7              KernSmooth_2.23-20       
## [105] crayon_1.5.1              htmltools_0.5.2          
## [107] mgcv_1.8-40               later_1.3.0              
## [109] tzdb_0.3.0                Formula_1.2-4            
## [111] DBI_1.1.2                 dbplyr_2.2.0             
## [113] MASS_7.3-57               rappdirs_0.3.3           
## [115] Matrix_1.4-2              cli_3.3.0                
## [117] metapod_1.2.0             parallel_4.1.2           
## [119] igraph_1.3.2              pkgconfig_2.0.3          
## [121] scuttle_1.4.0             plotly_4.10.0            
## [123] spatstat.sparse_2.1-1     vipor_0.4.5              
## [125] bslib_0.3.1               dqrng_0.3.0              
## [127] XVector_0.34.0            digest_0.6.29            
## [129] sctransform_0.3.3         RcppAnnoy_0.0.19         
## [131] spatstat.data_2.2-0       rmarkdown_2.14           
## [133] leiden_0.4.2              edgeR_3.36.0             
## [135] uwot_0.1.11               DelayedMatrixStats_1.16.0
## [137] maxLik_1.5-2              curl_4.3.2               
## [139] shiny_1.7.1               lifecycle_1.0.1          
## [141] nlme_3.1-158              jsonlite_1.8.0           
## [143] Rhdf5lib_1.16.0           BiocNeighbors_1.12.0     
## [145] limma_3.50.3              viridisLite_0.4.0        
## [147] fansi_1.0.3               pillar_1.7.0             
## [149] lattice_0.20-45           ggrastr_1.0.1            
## [151] fastmap_1.1.0             httr_1.4.3               
## [153] survival_3.3-1            glue_1.6.2               
## [155] png_0.1-7                 bluster_1.4.0            
## [157] bit_4.0.4                 stringi_1.7.6            
## [159] sass_0.4.1                blob_1.2.3               
## [161] BiocSingular_1.10.0       DirichletReg_0.7-1       
## [163] memoise_2.0.1             irlba_2.3.5              
## [165] future.apply_1.9.0