Analysis of CRC Spatial Transcriptomics data from 10X:Correlation Analysis of Pathway Activities and Cell Type abundances for CMS
Alberto Valdeolivas Urbelz (PS-Biomics&Path)
31 October, 2022
1 Introduction
The present script takes the the Spatial Omics data generated with 10X technology on a set of CRC samples. In particular, it explores the correlation between pathway activities and CMS cell type proportion per spot.
# library(dplyr, lib.loc = "/pstore/apps/bioinfo/R/4.0.1-foss")
library(CCA)
library(ggplot2)
library(patchwork, lib.loc = "/apps/rocs/2020.08/cascadelake/software/R/4.1.2-foss-2020a/lib64/R/library")
library(Seurat)
library(readr)
library(stringr)
library(tibble)
library(ggrepel)
library(RColorBrewer)
library(pheatmap)
source(file = "../WrapperFunction/SeuratWrappers.R")
data_directory <- params$data_directory
analysis_name <- params$analysis_name
input_folder_seurat <- "ST_LiverMetastases_pub/IntermediaryFiles/"
input_names_seurat <- "SeuratList_Clusters_Res05_Progeny.rds"
files_to_read_seurat <-
paste0(data_directory, analysis_name,input_folder_seurat,input_names_seurat)
input_folder_C2L <-
"Cell2Location/results/LocationModelLinearDependentWMultiExperiment_8experiments_36clusters_31296locations_4188genesLiverMetastasis/"
input_file_C2L <- "W_mRNA_count_q05.csv"
files_to_read_C2L <-
paste0(data_directory, analysis_name,input_folder_C2L,input_file_C2L)We read the Seurat objects containing the pathway activity.
seurat_objects <- readRDS(files_to_read_seurat)We Load Cell2Location results cotaning the abundance of the different cell types per spot.
results_STLiver_Korean <- read_csv(files_to_read_C2L)
sum_UMI_spot <- rowSums(results_STLiver_Korean[,c(2: (ncol(results_STLiver_Korean)-1))])
results_STLiver_Korean[,-1] <- results_STLiver_Korean[,-1]/sum_UMI_spot
colnames(results_STLiver_Korean) <- str_replace(colnames(results_STLiver_Korean),
pattern = "q05_nUMI_factors", replacement = "")
## We are only going to consider the CMS signatures
results_STLiver_Korean <-
results_STLiver_Korean[,colnames(results_STLiver_Korean) %in%
c("spot_id","CMS1", "CMS2", "CMS3", "CMS4")]2 Results taking into account correlation between the spots of colon untreated samples
We display the different correlation plots when this is globally computed per all the samples.
seurat_objects_CRC_Unt <-
seurat_objects[names(seurat_objects) %in%
c("ST-colon1", "ST-colon2") == TRUE]list_df_pathway <- lapply(seurat_objects_CRC_Unt, function(x){
matrix <- GetAssayData(x)
colnames(matrix) <- paste0(x$orig.ident,"_", colnames(matrix))
return(matrix %>% t() %>% as.data.frame() %>% tibble::rownames_to_column(var = "spot_id"))
})
pathway_df_allsamples <- do.call(rbind.data.frame, list_df_pathway)
merge_CL2_pathways <- pathway_df_allsamples %>%
dplyr::inner_join(results_STLiver_Korean, by ="spot_id") %>%
dplyr::arrange(spot_id)
X <- merge_CL2_pathways[,2:15]
Y <- merge_CL2_pathways[,16:19]
correl <- matcor(X, Y)
canonical_correlation2 <- cc(X, Y)
cat( "\n\n")cat( "## Pathway activities Correlation: ", "\n\n")2.1 Pathway activities Correlation:
MatrixCorrX <- correl$Xcor %>%
as.data.frame() %>%
tibble::rownames_to_column(var ="pathway") %>%
tidyr::pivot_longer(!pathway, values_drop_na = TRUE) # %>%
# dplyr::mutate(value = round(value,2))
p_cor_pathways <- ggplot(data = MatrixCorrX, aes(pathway , name , fill = value))+
geom_tile(color = "white")+
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson\nCorrelation") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 9, hjust = 1))+
xlab("Pathways") + ylab("Pathways") +
coord_fixed()
print(p_cor_pathways)
cat( "\n\n")cat( "## Cell type abundances Correlation: ", "\n\n")2.2 Cell type abundances Correlation:
MatrixCorrY <- correl$Ycor %>%
as.data.frame() %>%
tibble::rownames_to_column(var ="CellType") %>%
tidyr::pivot_longer(!CellType, values_drop_na = TRUE)
p_cor_celltype <- ggplot(MatrixCorrY, aes(name , CellType , fill = value))+
geom_tile(color = "white")+
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson\nCorrelation") +
theme_minimal()+
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 10, hjust = 1))+
xlab("Cell Type") + ylab("Cell Type") +
coord_fixed()
print(p_cor_celltype)
cat( "\n\n")cat( "## Cross Pathway-Cell type abundances Correlation: ", "\n\n")2.3 Cross Pathway-Cell type abundances Correlation:
CrossCorrelation <- correl$XYcor[1:14,15:18]
MatrixCorrXY <- CrossCorrelation %>%
as.data.frame() %>%
tibble::rownames_to_column(var ="var1") %>%
tidyr::pivot_longer(!var1, values_drop_na = TRUE) %>%
dplyr::mutate(value = round(value,2))
p_CrossCorrelation <- ggplot(data = MatrixCorrXY, aes(var1,name , fill = value))+
geom_tile(color = "white")+
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson\nCorrelation") +
theme_minimal()+
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 10, hjust = 1))+
xlab("Pathways") + ylab("Cell Type") +
coord_fixed()
print(p_CrossCorrelation)
mypalette <- rev(brewer.pal(n = 9, name ="RdBu"))
par(mar=c(10,10,10,10))
pheatmap(t(CrossCorrelation), treeheight_row=0,treeheight_col=0,
color =colorRampPalette(mypalette)(100),cluster_rows = FALSE,
fontsize = 16, cluster_cols = TRUE, breaks=seq(-1,1,length.out =
length(colorRampPalette(mypalette)(100)) + 1))
cat( "\n\n")cat( "## Variables representation on the canonical variates: ", "\n\n")2.4 Variables representation on the canonical variates:
d1 = 2
d2 = 1
df2 = data.frame(comp1 = c(canonical_correlation2$scores$corr.X.xscores[, d1],
canonical_correlation2$scores$corr.Y.xscores[, d1]),
comp2 = c(canonical_correlation2$scores$corr.X.xscores[, d2],
canonical_correlation2$scores$corr.Y.xscores[, d2]),
color = c(rep("Pathways", ncol(X)),
(rep("Cell Types", ncol(Y)))))
p_Canonival_variates <-ggplot(df2, aes(comp1, comp2), colour = color) +
geom_label_repel(aes(colour = color, label= rownames(df2)),
vjust=0,size = 3, family = "sans", direction= "both") +
theme_bw() +
theme(legend.position="none") +
theme(legend.title=element_blank())
print(p_Canonival_variates)
cat( "\n\n") 3 Results taking into account correlation between the spots of colon samples
We display the different correlation plots when this is globally computed per all the samples.
seurat_objects_CRC <-
seurat_objects[names(seurat_objects) %in%
c("ST-colon1", "ST-colon2", "ST-colon3", "ST-colon4") == TRUE]list_df_pathway <- lapply(seurat_objects_CRC, function(x){
matrix <- GetAssayData(x)
colnames(matrix) <- paste0(x$orig.ident,"_", colnames(matrix))
return(matrix %>% t() %>% as.data.frame() %>% tibble::rownames_to_column(var = "spot_id"))
})
pathway_df_allsamples <- do.call(rbind.data.frame, list_df_pathway)
merge_CL2_pathways <- pathway_df_allsamples %>%
dplyr::inner_join(results_STLiver_Korean, by ="spot_id") %>%
dplyr::arrange(spot_id)
X <- merge_CL2_pathways[,2:15]
Y <- merge_CL2_pathways[,16:19]
correl <- matcor(X, Y)
canonical_correlation2 <- cc(X, Y)
cat( "\n\n")cat( "## Pathway activities Correlation: ", "\n\n")3.1 Pathway activities Correlation:
MatrixCorrX <- correl$Xcor %>%
as.data.frame() %>%
tibble::rownames_to_column(var ="pathway") %>%
tidyr::pivot_longer(!pathway, values_drop_na = TRUE) # %>%
# dplyr::mutate(value = round(value,2))
p_cor_pathways <- ggplot(data = MatrixCorrX, aes(pathway , name , fill = value))+
geom_tile(color = "white")+
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson\nCorrelation") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 9, hjust = 1))+
xlab("Pathways") + ylab("Pathways") +
coord_fixed()
print(p_cor_pathways)
cat( "\n\n")cat( "## Cell type abundances Correlation: ", "\n\n")3.2 Cell type abundances Correlation:
MatrixCorrY <- correl$Ycor %>%
as.data.frame() %>%
tibble::rownames_to_column(var ="CellType") %>%
tidyr::pivot_longer(!CellType, values_drop_na = TRUE)
p_cor_celltype <- ggplot(MatrixCorrY, aes(name , CellType , fill = value))+
geom_tile(color = "white")+
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson\nCorrelation") +
theme_minimal()+
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 10, hjust = 1))+
xlab("Cell Type") + ylab("Cell Type") +
coord_fixed()
print(p_cor_celltype)
cat( "\n\n")cat( "## Cross Pathway-Cell type abundances Correlation: ", "\n\n")3.3 Cross Pathway-Cell type abundances Correlation:
CrossCorrelation <- correl$XYcor[1:14,15:18]
MatrixCorrXY <- CrossCorrelation %>%
as.data.frame() %>%
tibble::rownames_to_column(var ="var1") %>%
tidyr::pivot_longer(!var1, values_drop_na = TRUE) %>%
dplyr::mutate(value = round(value,2))
p_CrossCorrelation <- ggplot(data = MatrixCorrXY, aes(var1,name , fill = value))+
geom_tile(color = "white")+
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson\nCorrelation") +
theme_minimal()+
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 10, hjust = 1))+
xlab("Pathways") + ylab("Cell Type") +
coord_fixed()
print(p_CrossCorrelation)
mypalette <- rev(brewer.pal(n = 9, name ="RdBu"))
par(mar=c(10,10,10,10))
pheatmap(t(CrossCorrelation), treeheight_row=0,treeheight_col=0,
color =colorRampPalette(mypalette)(100),cluster_rows = FALSE,
fontsize = 16, cluster_cols = TRUE, breaks=seq(-1,1,length.out =
length(colorRampPalette(mypalette)(100)) + 1))
cat( "\n\n")cat( "## Variables representation on the canonical variates: ", "\n\n")3.4 Variables representation on the canonical variates:
d1 = 2
d2 = 1
df2 = data.frame(comp1 = c(canonical_correlation2$scores$corr.X.xscores[, d1],
canonical_correlation2$scores$corr.Y.xscores[, d1]),
comp2 = c(canonical_correlation2$scores$corr.X.xscores[, d2],
canonical_correlation2$scores$corr.Y.xscores[, d2]),
color = c(rep("Pathways", ncol(X)),
(rep("Cell Types", ncol(Y)))))
p_Canonival_variates <-ggplot(df2, aes(comp1, comp2), colour = color) +
geom_label_repel(aes(colour = color, label= rownames(df2)),
vjust=0,size = 3, family = "sans", direction= "both") +
theme_bw() +
theme(legend.position="none") +
theme(legend.title=element_blank())
print(p_Canonival_variates)
cat( "\n\n") 4 Results taking into account correlation between all the samples samples
We display the different correlation plots when this is globally computed per all the samples.
list_df_pathway <- lapply(seurat_objects, function(x){
matrix <- GetAssayData(x)
colnames(matrix) <- paste0(x$orig.ident,"_", colnames(matrix))
return(matrix %>% t() %>% as.data.frame() %>% tibble::rownames_to_column(var = "spot_id"))
})
pathway_df_allsamples <- do.call(rbind.data.frame, list_df_pathway)
merge_CL2_pathways <- pathway_df_allsamples %>%
dplyr::inner_join(results_STLiver_Korean, by ="spot_id") %>%
dplyr::arrange(spot_id)
X <- merge_CL2_pathways[,2:15]
Y <- merge_CL2_pathways[,16:19]
correl <- matcor(X, Y)
canonical_correlation2 <- cc(X, Y)
cat( "\n\n")cat( "## Pathway activities Correlation: ", "\n\n")4.1 Pathway activities Correlation:
MatrixCorrX <- correl$Xcor %>%
as.data.frame() %>%
tibble::rownames_to_column(var ="pathway") %>%
tidyr::pivot_longer(!pathway, values_drop_na = TRUE) # %>%
# dplyr::mutate(value = round(value,2))
p_cor_pathways <- ggplot(data = MatrixCorrX, aes(pathway , name , fill = value))+
geom_tile(color = "white")+
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson\nCorrelation") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 9, hjust = 1))+
xlab("Pathways") + ylab("Pathways") +
coord_fixed()
print(p_cor_pathways)
cat( "\n\n")cat( "## Cell type abundances Correlation: ", "\n\n")4.2 Cell type abundances Correlation:
MatrixCorrY <- correl$Ycor %>%
as.data.frame() %>%
tibble::rownames_to_column(var ="CellType") %>%
tidyr::pivot_longer(!CellType, values_drop_na = TRUE)
p_cor_celltype <- ggplot(MatrixCorrY, aes(name , CellType , fill = value))+
geom_tile(color = "white")+
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson\nCorrelation") +
theme_minimal()+
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 10, hjust = 1))+
xlab("Cell Type") + ylab("Cell Type") +
coord_fixed()
print(p_cor_celltype)
cat( "\n\n")cat( "## Cross Pathway-Cell type abundances Correlation: ", "\n\n")4.3 Cross Pathway-Cell type abundances Correlation:
CrossCorrelation <- correl$XYcor[1:14,15:18]
MatrixCorrXY <- CrossCorrelation %>%
as.data.frame() %>%
tibble::rownames_to_column(var ="var1") %>%
tidyr::pivot_longer(!var1, values_drop_na = TRUE) %>%
dplyr::mutate(value = round(value,2))
p_CrossCorrelation <- ggplot(data = MatrixCorrXY, aes(var1,name , fill = value))+
geom_tile(color = "white")+
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson\nCorrelation") +
theme_minimal()+
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 10, hjust = 1))+
xlab("Pathways") + ylab("Cell Type") +
coord_fixed()
print(p_CrossCorrelation)
mypalette <- rev(brewer.pal(n = 9, name ="RdBu"))
par(mar=c(10,10,10,10))
pheatmap(t(CrossCorrelation), treeheight_row=0,treeheight_col=0,
color =colorRampPalette(mypalette)(100),cluster_rows = FALSE,
fontsize = 16, cluster_cols = TRUE, breaks=seq(-1,1,length.out =
length(colorRampPalette(mypalette)(100)) + 1))
pheatmap(CrossCorrelation, treeheight_row=0,treeheight_col=0,
color =colorRampPalette(mypalette)(100),cluster_cols = FALSE,
fontsize = 16, cluster_rows = TRUE, breaks=seq(-1,1,length.out =
length(colorRampPalette(mypalette)(100)) + 1))
cat( "\n\n")cat( "## Variables representation on the canonical variates: ", "\n\n")4.4 Variables representation on the canonical variates:
d1 = 2
d2 = 1
df2 = data.frame(comp1 = c(canonical_correlation2$scores$corr.X.xscores[, d1],
canonical_correlation2$scores$corr.Y.xscores[, d1]),
comp2 = c(canonical_correlation2$scores$corr.X.xscores[, d2],
canonical_correlation2$scores$corr.Y.xscores[, d2]),
color = c(rep("Pathways", ncol(X)),
(rep("Cell Types", ncol(Y)))))
p_Canonival_variates <-ggplot(df2, aes(comp1, comp2), colour = color) +
geom_label_repel(aes(colour = color, label= rownames(df2)),
vjust=0,size = 3, family = "sans", direction= "both") +
theme_bw() +
theme(legend.position="none") +
theme(legend.title=element_blank())
print(p_Canonival_variates)
cat( "\n\n") 5 Conclusion
6 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] splines stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] pheatmap_1.0.12 RColorBrewer_1.1-3 ggrepel_0.9.1
## [4] tibble_3.1.7 stringr_1.4.0 readr_2.1.2
## [7] sp_1.5-0 SeuratObject_4.1.0 Seurat_4.1.0
## [10] patchwork_1.1.1 ggplot2_3.3.6 CCA_1.2.1
## [13] fields_13.3 viridis_0.6.2 viridisLite_0.4.0
## [16] spam_2.8-0 fda_6.0.3 deSolve_1.32
## [19] fds_1.8 RCurl_1.98-1.7 rainbow_3.6
## [22] pcaPP_2.0-1 MASS_7.3-57 BiocManager_1.30.18
##
## loaded via a namespace (and not attached):
## [1] plyr_1.8.7 igraph_1.3.2 lazyeval_0.2.2
## [4] listenv_0.8.0 scattermore_0.8 digest_0.6.29
## [7] htmltools_0.5.2 fansi_1.0.3 magrittr_2.0.3
## [10] tensor_1.5 cluster_2.1.3 ks_1.13.5
## [13] ROCR_1.0-11 tzdb_0.3.0 hdrcde_3.4
## [16] globals_0.15.0 matrixStats_0.62.0 vroom_1.5.7
## [19] spatstat.sparse_2.1-1 colorspace_2.1-0 xfun_0.31
## [22] dplyr_1.0.9 crayon_1.5.1 jsonlite_1.8.0
## [25] progressr_0.10.1 spatstat.data_2.2-0 survival_3.3-1
## [28] zoo_1.8-10 glue_1.6.2 polyclip_1.10-0
## [31] gtable_0.3.0 leiden_0.4.2 future.apply_1.9.0
## [34] maps_3.4.0 abind_1.4-7 scales_1.2.0
## [37] mvtnorm_1.1-3 DBI_1.1.3 spatstat.random_2.2-0
## [40] miniUI_0.1.1.1 Rcpp_1.0.8.3 xtable_1.8-6
## [43] reticulate_1.25 spatstat.core_2.4-4 bit_4.0.4
## [46] mclust_5.4.10 dotCall64_1.0-1 htmlwidgets_1.5.4
## [49] httr_1.4.3 ellipsis_0.3.2 ica_1.0-2
## [52] farver_2.1.0 pkgconfig_2.0.3 sass_0.4.1
## [55] uwot_0.1.11 deldir_1.0-6 utf8_1.2.2
## [58] labeling_0.4.2 tidyselect_1.1.2 rlang_1.0.2
## [61] reshape2_1.4.4 later_1.3.0 munsell_0.5.0
## [64] tools_4.1.2 cli_3.3.0 generics_0.1.2
## [67] ggridges_0.5.3 evaluate_0.15 fastmap_1.1.0
## [70] yaml_2.3.5 goftest_1.2-3 bit64_4.0.5
## [73] knitr_1.39 fitdistrplus_1.1-8 purrr_0.3.4
## [76] RANN_2.6.1 pbapply_1.5-0 future_1.26.1
## [79] nlme_3.1-158 mime_0.12 pracma_2.4.1
## [82] compiler_4.1.2 rstudioapi_0.13 plotly_4.10.0
## [85] png_0.1-7 spatstat.utils_2.3-1 bslib_0.3.1
## [88] stringi_1.7.6 highr_0.9 rgeos_0.5-10
## [91] lattice_0.20-45 Matrix_1.4-2 vctrs_0.4.1
## [94] pillar_1.7.0 lifecycle_1.0.1 spatstat.geom_2.4-0
## [97] lmtest_0.9-40 jquerylib_0.1.4 RcppAnnoy_0.0.19
## [100] data.table_1.14.2 cowplot_1.1.1 bitops_1.0-7
## [103] irlba_2.3.5 httpuv_1.6.5 R6_2.5.1
## [106] promises_1.2.0.1 KernSmooth_2.23-20 gridExtra_2.3
## [109] parallelly_1.32.0 codetools_0.2-18 assertthat_0.2.1
## [112] withr_2.5.0 sctransform_0.3.3 hms_1.1.1
## [115] mgcv_1.8-40 parallel_4.1.2 grid_4.1.2
## [118] rpart_4.1.16 tidyr_1.2.0 rmarkdown_2.14
## [121] Rtsne_0.16 shiny_1.7.1