Functional analysis with MISTy - pathway activity and ligand expression
Leoni Zimmermann
Heidelberg University, Heidelberg, GermanyJovan Tanevski
Heidelberg University and Heidelberg University Hospital, Heidelberg, GermanyJožef Stefan Institute, Ljubljana, Sloveniajovan.tanevski@uni-heidelberg.de
2024-03-11
Last updated: 2024-03-11
Checks: 7 0
Knit directory:
ProtocolLabRotationSaezRodriguezGroup/
This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20240306) was run prior to running
the code in the R Markdown file. Setting a seed ensures that any results
that rely on randomness, e.g. subsampling or permutations, are
reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version 9fe1532. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for
the analysis have been committed to Git prior to generating the results
(you can use wflow_publish or
wflow_git_commit). workflowr only checks the R Markdown
file, but you know if there are other scripts or data files that it
depends on. Below is the status of the Git repository when the results
were generated:
Untracked files:
Untracked: 10X_Visium_ACH005.tar.gz
Untracked: ACH005/
Untracked: bc_metadata.tsv
Untracked: data/10X_Visium_ACH005.tar.gz
Untracked: hca_p14.rds
Untracked: imc_bc_optim_zoi.RDS
Untracked: omni_resource.csv
Untracked: omnipathr-log/
Untracked: result/
Untracked: results/
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the repository in which changes were
made to the R Markdown
(analysis/FunctionalPipelinePathwayActivityLigands.Rmd) and
HTML (docs/FunctionalPipelinePathwayActivityLigands.html)
files. If you’ve configured a remote Git repository (see
?wflow_git_remote), click on the hyperlinks in the table
below to view the files as they were in that past version.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| html | a7395df | leotenshii | 2024-03-11 | Build site. |
| html | 64a10dc | leotenshii | 2024-03-11 | added html |
| Rmd | 38a9f93 | leotenshii | 2024-03-11 | Revert "changed paths to data/…" |
| Rmd | 5e9bf15 | leotenshii | 2024-03-11 | changed paths to data/… |
| Rmd | 5ec41b8 | leotenshii | 2024-03-11 | upload of vignettes |
Introduction
10X Visium captures spatially resolved transcriptomic profiles in spots containing multiple cells. In this vignette, we will use the gene expression information from Visium data to infer pathway activity and investigate spatial relationships between some pathways and ligand expression.
Load the necessary packages:
# MISTy
library(mistyR)
library(future)
#Seurat
library(Seurat)
library(SeuratObject)
# Data manipulation
library(tidyverse)
# Pathways and annotation
library(decoupleR)
library(OmnipathR)
library(progeny)
# Cleaning names
library(janitor)Get and load data
For this showcase, we use a 10X Visium spatial slide from Kuppe et al., 2022, where they created a spatial multi-omic map of human myocardial infarction. The tissue example data comes from the human heart of patient 14, which is in a later state after myocardial infarction. The Seurat object contains, among other things, the normalized and raw gene counts. First, we have to download and extract the file:
download.file("https://zenodo.org/records/6580069/files/10X_Visium_ACH005.tar.gz?download=1",
destfile = "10X_Visium_ACH005.tar.gz", method = "curl")
untar("10X_Visium_ACH005.tar.gz")The next step is to load the data, extract the normalized gene counts, names of genes expressed in at least 5% of the spots, and pixel coordinates. It is recommended to use pixel coordinates instead of row and column numbers since the rows are shifted and therefore do not express the real distance between the spots.
seurat_vs <- readRDS("ACH005/ACH005.rds")
expression <- as.matrix(GetAssayData(seurat_vs, layer = "counts", assay = "SCT"))
gene_names <- rownames(expression[(rowSums(expression > 0) / ncol(expression)) >= 0.05,])
geometry <- GetTissueCoordinates(seurat_vs, scale = NULL)Pathway activity
Now we create a Seurat object with pathway activities inferred from
PROGENy.
We delete the PROGENy assay done by Kuppe et al. and load a model matrix
with the top 1000 significant genes for each of the 15 available
pathways. We then extract the genes that are both common to the PROGENy
model and the snRNA-seq assay from the Seurat object. We compute the
weighted sum of both and scale them to infer the pathway activity. We
save the result in a Seurat assay and clean the row names to handle
problematic variables.
seurat_vs[['progeny']] <- NULL
# Matrix with important genes for each pathway
model <- get_progeny(organism = "human", top = 1000)
# Use multivariate linear model to estimate activity
est_path_act <- run_mlm(expression, model,.mor = NULL)
# Put estimated pathway activities object into the correct format
est_path_act_wide <- est_path_act %>%
pivot_wider(id_cols = condition, names_from = source, values_from = score) %>%
column_to_rownames("condition")
# Clean names
colnames(est_path_act_wide) <- est_path_act_wide %>%
clean_names(parsing_option = 0) %>%
colnames(.)
# Create a Seurat object
seurat_vs[['progeny']] <- CreateAssayObject(counts = t(est_path_act_wide))
# Format for running MISTy later
pathway_activity <- t(as.matrix(GetAssayData(seurat_vs, "progeny")))Ligands
To annotate the expressed ligands found in the tissue slide, we import an intercellular network of ligands and receptors from Omnipath. We extract the ligands that are expressed in the tissue slide and get their count data. We again clean the row names to handle problematic variables.
# Get ligands
lig_rec <- import_intercell_network(interactions_param = list(datasets = c('ligrecextra', 'omnipath', 'pathwayextra')),
transmitter_param = list(parent = 'ligand'),
receiver_param = list(parent = 'receptor'))
# Get unique ligands
ligands <- unique(lig_rec$source_genesymbol)
# Get expression of ligands in slide
slide_markers <- ligands[ligands %in% gene_names]
ligand_expr <- t(as.matrix(expression[slide_markers,])) %>% clean_names()
#clean names
rownames(seurat_vs@assays$SCT@data) <- seurat_vs@assays$SCT@data %>% clean_names(parsing_option = 0) %>% rownames(.)Visualize pathway activity
Before continuing with creating the MISTy view, we can look at the slide itself and some of the pathway activities.
#Slide
SpatialPlot(seurat_vs, alpha = 0)
| Version | Author | Date |
|---|---|---|
| 64a10dc | leotenshii | 2024-03-11 |
# Pathway activity examples
DefaultAssay(seurat_vs) <- "progeny"
SpatialFeaturePlot(seurat_vs, feature = c("mapk", "p53"), keep.scale = NULL)
| Version | Author | Date |
|---|---|---|
| 64a10dc | leotenshii | 2024-03-11 |
MISTy views
Now we need to create the MISTy views of interest. We are interested in the relationship of the pathway activity in the same spot (intraview) and the ten closest spots (paraview). Therefore we choose the family `constant` and set l to ten, which will select the ten nearest neighbors. Depending on the goal of the analysis, different families can be applied.
We are also intrigued about the relationship of ligand expression and pathway activity in the broader tissue. For this, we again create an intra- and paraview, this time for the expression of the ligands, but from this view, we only need the paraview. In the next step, we add it to the pathway activity views to achieve our intended view composition.
pathway_act_view <- create_initial_view(as_tibble(pathway_activity) ) %>%
add_paraview(geometry, l = 10, family = "constant")
Generating paraview using 10 nearest neighbors per unitligand_view <- create_initial_view(as_tibble(ligand_expr) %>% clean_names()) %>%
add_paraview(geometry, l = 10, family = "constant")
Generating paraview using 10 nearest neighbors per unitcombined_views <- pathway_act_view %>% add_views(create_view("paraview.ligand.10", ligand_view[["paraview.10"]]$data, "para.ligand.10"))Then run MISTy and collect the results:
run_misty(combined_views, "result/functional_ligand")[1] "F:/LabRotationSaez/ProtocolLabRotationSaezRodriguezGroup/result/functional_ligand"misty_results <- collect_results("result/functional_ligand/")Downstream analysis
With the collected results, we can now answer the following questions:
1. To what extent can the analyzed surrounding tissues’ activities explain the pathway activity of the spot compared to the intraview?
Here we can look at two different statistics: multi.R2 shows the total variance explained by the multiview model. gain.R2 shows the increase in explainable variance from the paraview.
misty_results %>%
plot_improvement_stats("gain.R2") %>%
plot_improvement_stats("multi.R2")Warning: Removed 14 rows containing missing values or values outside the scale range
(`geom_segment()`).
| Version | Author | Date |
|---|---|---|
| 64a10dc | leotenshii | 2024-03-11 |
Warning: Removed 14 rows containing missing values or values outside the scale range
(`geom_segment()`).
| Version | Author | Date |
|---|---|---|
| 64a10dc | leotenshii | 2024-03-11 |
The paraview particularly increases the explained variance for TGFb and PI3K. In general, the significant gain in R2 can be interpreted as the following:
“We can better explain the expression of marker X when we consider additional views other than the intrinsic view.”
To see the individual contributions of the views we can use:
misty_results %>% plot_view_contributions()
| Version | Author | Date |
|---|---|---|
| 64a10dc | leotenshii | 2024-03-11 |
2. What are the specific relations that can explain the pathway activity
The importance of the markers from each viewpoint as predictors of the spot intrinsic pathway activity can be shown individually to explain the contributions.
First, for the intrinsic view. To set an importance threshold we
apply cutoff:
misty_results %>%
plot_interaction_heatmap("intra", clean = TRUE, cutoff = 1.5)
| Version | Author | Date |
|---|---|---|
| 64a10dc | leotenshii | 2024-03-11 |
We can observe that TNFa is a significant predictor for the activity of the NFkB pathway when in the same spot. Let’s take a look at the spatial distribution of these pathway activities in the tissue slide:
SpatialFeaturePlot(seurat_vs, features = c("tnfa", "nfkb"), image.alpha = 0)
| Version | Author | Date |
|---|---|---|
| 64a10dc | leotenshii | 2024-03-11 |
We can observe a correlation between high TNFa activity and high NFkB activity.
Now we repeat this analysis with the pathway activity paraview. With
trim we display only targets with a value above 0.5 for
gain.R2.
misty_results %>%
plot_interaction_heatmap(view = "para.10",
clean = TRUE,
trim = 0.5,
trim.measure = "gain.R2",
cutoff = 1.25)
| Version | Author | Date |
|---|---|---|
| 64a10dc | leotenshii | 2024-03-11 |
From the gain.R2 we know that the paraview contributes a
lot to explaining the TGFb pathway activity. Let’s visualize it and its
most important predictor, androgen pathway activity:
SpatialFeaturePlot(seurat_vs, features = c("androgen", "tgfb"), image.alpha = 0)
| Version | Author | Date |
|---|---|---|
| 64a10dc | leotenshii | 2024-03-11 |
The plots show us an anticorrelation of these pathways.
Now we will analyze the last view, the ligand expression paraview:
misty_results %>%
plot_interaction_heatmap(view = "para.ligand.10", clean = TRUE, trim = 0.5,
trim.measure = "gain.R2", cutoff=3)
| Version | Author | Date |
|---|---|---|
| 64a10dc | leotenshii | 2024-03-11 |
The ligand SERPINF1 is a predictor of both PI3K and VEGF:
SpatialFeaturePlot(seurat_vs, features = c("pi3k"), image.alpha = 0)
| Version | Author | Date |
|---|---|---|
| 64a10dc | leotenshii | 2024-03-11 |
SpatialFeaturePlot(seurat_vs, features = c("vegf"), image.alpha = 0)
| Version | Author | Date |
|---|---|---|
| 64a10dc | leotenshii | 2024-03-11 |
DefaultAssay(seurat_vs) <- "SCT"
SpatialFeaturePlot(seurat_vs, features = c("serpinf1"), image.alpha = 0)
| Version | Author | Date |
|---|---|---|
| 64a10dc | leotenshii | 2024-03-11 |
From the slides we can see that SERPINF1 correlates positive with PI3K and negative with VEGF
See also
browseVignettes("mistyR")
Session Info
Here is the output of sessionInfo() at the point when
this document was compiled.
sessionInfo()R version 4.3.2 (2023-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)
Matrix products: default
locale:
[1] LC_COLLATE=German_Germany.utf8 LC_CTYPE=German_Germany.utf8
[3] LC_MONETARY=German_Germany.utf8 LC_NUMERIC=C
[5] LC_TIME=German_Germany.utf8
time zone: Europe/Berlin
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] distances_0.1.10 janitor_2.2.0 progeny_1.24.0 OmnipathR_3.10.1
[5] decoupleR_2.8.0 lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1
[9] dplyr_1.1.4 purrr_1.0.2 readr_2.1.5 tidyr_1.3.0
[13] tibble_3.2.1 ggplot2_3.5.0 tidyverse_2.0.0 Seurat_5.0.1
[17] SeuratObject_5.0.1 sp_2.1-2 future_1.33.1 mistyR_1.10.0
[21] workflowr_1.7.1
loaded via a namespace (and not attached):
[1] RcppAnnoy_0.0.21 splines_4.3.2 later_1.3.2
[4] filelock_1.0.3 R.oo_1.26.0 cellranger_1.1.0
[7] polyclip_1.10-6 hardhat_1.3.1 pROC_1.18.5
[10] rpart_4.1.23 fastDummies_1.7.3 lifecycle_1.0.4
[13] rprojroot_2.0.4 globals_0.16.3 processx_3.8.3
[16] lattice_0.22-5 vroom_1.6.5 MASS_7.3-60
[19] backports_1.4.1 magrittr_2.0.3 plotly_4.10.3
[22] sass_0.4.8 rmarkdown_2.25 jquerylib_0.1.4
[25] yaml_2.3.8 rlist_0.4.6.2 httpuv_1.6.13
[28] sctransform_0.4.1 spam_2.10-0 spatstat.sparse_3.0-3
[31] reticulate_1.34.0 cowplot_1.1.2 pbapply_1.7-2
[34] RColorBrewer_1.1-3 abind_1.4-5 rvest_1.0.3
[37] Rtsne_0.17 R.utils_2.12.3 nnet_7.3-19
[40] rappdirs_0.3.3 ipred_0.9-14 git2r_0.33.0
[43] lava_1.8.0 ggrepel_0.9.4 irlba_2.3.5.1
[46] listenv_0.9.1 spatstat.utils_3.0-4 goftest_1.2-3
[49] RSpectra_0.16-1 spatstat.random_3.2-2 fitdistrplus_1.1-11
[52] parallelly_1.37.1 leiden_0.4.3.1 codetools_0.2-19
[55] xml2_1.3.6 tidyselect_1.2.0 farver_2.1.1
[58] stats4_4.3.2 matrixStats_1.2.0 spatstat.explore_3.2-5
[61] jsonlite_1.8.8 caret_6.0-94 ellipsis_0.3.2
[64] progressr_0.14.0 iterators_1.0.14 ggridges_0.5.5
[67] survival_3.5-7 foreach_1.5.2 tools_4.3.2
[70] progress_1.2.3 ica_1.0-3 Rcpp_1.0.11
[73] glue_1.6.2 prodlim_2023.08.28 gridExtra_2.3
[76] ranger_0.16.0 xfun_0.41 withr_3.0.0
[79] fastmap_1.1.1 fansi_1.0.6 callr_3.7.3
[82] digest_0.6.33 timechange_0.3.0 R6_2.5.1
[85] mime_0.12 colorspace_2.1-0 scattermore_1.2
[88] tensor_1.5 spatstat.data_3.0-3 R.methodsS3_1.8.2
[91] utf8_1.2.4 generics_0.1.3 recipes_1.0.10
[94] data.table_1.15.2 class_7.3-22 ridge_3.3
[97] prettyunits_1.2.0 httr_1.4.7 htmlwidgets_1.6.4
[100] whisker_0.4.1 ModelMetrics_1.2.2.2 uwot_0.1.16
[103] pkgconfig_2.0.3 gtable_0.3.4 timeDate_4032.109
[106] lmtest_0.9-40 selectr_0.4-2 furrr_0.3.1
[109] htmltools_0.5.7 dotCall64_1.1-1 scales_1.3.0
[112] png_0.1-8 gower_1.0.1 snakecase_0.11.1
[115] knitr_1.45 rstudioapi_0.15.0 tzdb_0.4.0
[118] reshape2_1.4.4 checkmate_2.3.1 nlme_3.1-164
[121] curl_5.2.0 cachem_1.0.8 zoo_1.8-12
[124] KernSmooth_2.23-22 parallel_4.3.2 miniUI_0.1.1.1
[127] pillar_1.9.0 grid_4.3.2 logger_0.2.2
[130] vctrs_0.6.5 RANN_2.6.1 promises_1.2.1
[133] xtable_1.8-4 cluster_2.1.6 archive_1.1.7
[136] evaluate_0.23 cli_3.6.2 compiler_4.3.2
[139] rlang_1.1.2 crayon_1.5.2 future.apply_1.11.1
[142] labeling_0.4.3 ps_1.7.5 getPass_0.2-4
[145] plyr_1.8.9 fs_1.6.3 stringi_1.8.3
[148] viridisLite_0.4.2 deldir_2.0-4 assertthat_0.2.1
[151] munsell_0.5.0 lazyeval_0.2.2 spatstat.geom_3.2-7
[154] Matrix_1.6-4 RcppHNSW_0.5.0 hms_1.1.3
[157] patchwork_1.2.0.9000 bit64_4.0.5 shiny_1.8.0
[160] highr_0.10 ROCR_1.0-11 igraph_1.6.0
[163] bslib_0.6.1 bit_4.0.5 readxl_1.4.3
sessionInfo()R version 4.3.2 (2023-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)
Matrix products: default
locale:
[1] LC_COLLATE=German_Germany.utf8 LC_CTYPE=German_Germany.utf8
[3] LC_MONETARY=German_Germany.utf8 LC_NUMERIC=C
[5] LC_TIME=German_Germany.utf8
time zone: Europe/Berlin
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] distances_0.1.10 janitor_2.2.0 progeny_1.24.0 OmnipathR_3.10.1
[5] decoupleR_2.8.0 lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1
[9] dplyr_1.1.4 purrr_1.0.2 readr_2.1.5 tidyr_1.3.0
[13] tibble_3.2.1 ggplot2_3.5.0 tidyverse_2.0.0 Seurat_5.0.1
[17] SeuratObject_5.0.1 sp_2.1-2 future_1.33.1 mistyR_1.10.0
[21] workflowr_1.7.1
loaded via a namespace (and not attached):
[1] RcppAnnoy_0.0.21 splines_4.3.2 later_1.3.2
[4] filelock_1.0.3 R.oo_1.26.0 cellranger_1.1.0
[7] polyclip_1.10-6 hardhat_1.3.1 pROC_1.18.5
[10] rpart_4.1.23 fastDummies_1.7.3 lifecycle_1.0.4
[13] rprojroot_2.0.4 globals_0.16.3 processx_3.8.3
[16] lattice_0.22-5 vroom_1.6.5 MASS_7.3-60
[19] backports_1.4.1 magrittr_2.0.3 plotly_4.10.3
[22] sass_0.4.8 rmarkdown_2.25 jquerylib_0.1.4
[25] yaml_2.3.8 rlist_0.4.6.2 httpuv_1.6.13
[28] sctransform_0.4.1 spam_2.10-0 spatstat.sparse_3.0-3
[31] reticulate_1.34.0 cowplot_1.1.2 pbapply_1.7-2
[34] RColorBrewer_1.1-3 abind_1.4-5 rvest_1.0.3
[37] Rtsne_0.17 R.utils_2.12.3 nnet_7.3-19
[40] rappdirs_0.3.3 ipred_0.9-14 git2r_0.33.0
[43] lava_1.8.0 ggrepel_0.9.4 irlba_2.3.5.1
[46] listenv_0.9.1 spatstat.utils_3.0-4 goftest_1.2-3
[49] RSpectra_0.16-1 spatstat.random_3.2-2 fitdistrplus_1.1-11
[52] parallelly_1.37.1 leiden_0.4.3.1 codetools_0.2-19
[55] xml2_1.3.6 tidyselect_1.2.0 farver_2.1.1
[58] stats4_4.3.2 matrixStats_1.2.0 spatstat.explore_3.2-5
[61] jsonlite_1.8.8 caret_6.0-94 ellipsis_0.3.2
[64] progressr_0.14.0 iterators_1.0.14 ggridges_0.5.5
[67] survival_3.5-7 foreach_1.5.2 tools_4.3.2
[70] progress_1.2.3 ica_1.0-3 Rcpp_1.0.11
[73] glue_1.6.2 prodlim_2023.08.28 gridExtra_2.3
[76] ranger_0.16.0 xfun_0.41 withr_3.0.0
[79] fastmap_1.1.1 fansi_1.0.6 callr_3.7.3
[82] digest_0.6.33 timechange_0.3.0 R6_2.5.1
[85] mime_0.12 colorspace_2.1-0 scattermore_1.2
[88] tensor_1.5 spatstat.data_3.0-3 R.methodsS3_1.8.2
[91] utf8_1.2.4 generics_0.1.3 recipes_1.0.10
[94] data.table_1.15.2 class_7.3-22 ridge_3.3
[97] prettyunits_1.2.0 httr_1.4.7 htmlwidgets_1.6.4
[100] whisker_0.4.1 ModelMetrics_1.2.2.2 uwot_0.1.16
[103] pkgconfig_2.0.3 gtable_0.3.4 timeDate_4032.109
[106] lmtest_0.9-40 selectr_0.4-2 furrr_0.3.1
[109] htmltools_0.5.7 dotCall64_1.1-1 scales_1.3.0
[112] png_0.1-8 gower_1.0.1 snakecase_0.11.1
[115] knitr_1.45 rstudioapi_0.15.0 tzdb_0.4.0
[118] reshape2_1.4.4 checkmate_2.3.1 nlme_3.1-164
[121] curl_5.2.0 cachem_1.0.8 zoo_1.8-12
[124] KernSmooth_2.23-22 parallel_4.3.2 miniUI_0.1.1.1
[127] pillar_1.9.0 grid_4.3.2 logger_0.2.2
[130] vctrs_0.6.5 RANN_2.6.1 promises_1.2.1
[133] xtable_1.8-4 cluster_2.1.6 archive_1.1.7
[136] evaluate_0.23 cli_3.6.2 compiler_4.3.2
[139] rlang_1.1.2 crayon_1.5.2 future.apply_1.11.1
[142] labeling_0.4.3 ps_1.7.5 getPass_0.2-4
[145] plyr_1.8.9 fs_1.6.3 stringi_1.8.3
[148] viridisLite_0.4.2 deldir_2.0-4 assertthat_0.2.1
[151] munsell_0.5.0 lazyeval_0.2.2 spatstat.geom_3.2-7
[154] Matrix_1.6-4 RcppHNSW_0.5.0 hms_1.1.3
[157] patchwork_1.2.0.9000 bit64_4.0.5 shiny_1.8.0
[160] highr_0.10 ROCR_1.0-11 igraph_1.6.0
[163] bslib_0.6.1 bit_4.0.5 readxl_1.4.3