本文探讨了空间细胞类型共定位在疾病研究中的应用，特别是在 Spatial multi-omic map of human myocardial infarction 文章中的实践。文章指出，整体细胞类型定位揭示了生态位关系，而细胞亚群分析揭示了疾病状态转变。重点分享了MISTy工具用于共定位分析的代码，并提供了完整的分析代码链接。
作者，Evil Genius
在之前的分享中，分享了空间共定位的实现方法，但是今天回看文章的时候，又有了一点新的体会，看来温故而知新，实时温习很有作用。
我们来看看文章Spatial multi-omic map of human myocardial infarction, 2022年8月发表于nature，文章很多的代码都很经典，包括单细胞空间处理全分析代码，我放在了最后，大家可以随意拿走。
我们来看看文章对空间细胞类型共定位的运用
首先第一次的运用，是所有细胞类型之间的定位关系，识别细胞类型之间的空间生态位。
注意这里是细胞类型的整体定位，我们来看一下文章的介绍
evaluated three different neighbourhood area sizes using MISTy: (1) the importance of cell-type abundances within a spot (colocalization), (2) in the local neighbourhood (radius of 1 spot), and (3) in an extended neighbourhood that expanded to a radius of 15 spots. We observed that endothelial cells were the most predictive of the abundance of vSMCs, pericytes, adipocytes and cardiomyocytes within all spots, probably reflecting dependencies between cell types of the vasculature
初步的结论是淋巴细胞和髓细胞在免疫细胞浸润和炎症区表现出很强的依赖性（similarly captured by cell-type niche 5）。值得注意的是，观察到骨髓细胞和成纤维细胞之间存在很强的依赖性，它们在生态位4中高度共富集，这与巨噬细胞在成纤维细胞激活中和成纤维细胞在巨噬细胞吸引中所起的关键作用一致。在邻近和延伸的相邻SPOT之间，观察到与心脏血管相关的细胞(vSMCs，内皮细胞，周细胞和成纤维细胞)之间的依赖性更强，表明心肌血管网络主导心脏组织结构组织。
当然这里文章也说了一下，这是整体的情况，在不同的疾病条件下显然不可能是所有的某种细胞都发生了状态转变。所以在后续的分析中全部分析的是细胞亚群的生态位分析。
包括CMs细胞
Endo细胞
Fib
免疫细胞
也就是在整体的情况下，初步分析细胞类型的共定位情况，然后后续再分析不同疾病条件下细胞类型亚群的分布，每种亚群其实代表了不同的疾病状态转变，这种状态的转变，也会带来生态位的变化，这种空间排布的变化，解释了疾病的发生状态。
最后将整理好的共定位代码分享给大家MISTy，最后放上文章的所有代码
library(argparse)
library(tidyverse)
library(Seurat)
library(mistyR)
source("misty_utilities.R")
parser = ArgumentParser()
parser$add_argument("--rds", help="rds file,sc and sp joint",required =T)
parser$add_argument("--outdir", help="outdir,Absolute path",required =T)
args <- parser$parse_args()
rds = args$rds
outdir = args$outdir
future::plan(future::multisession)
run_colocalization <- function(slide,
assay,
useful_features,
out_label,
misty_out_alias = outdir) { ###输出目录大家自己制定
# Define assay of each view ---------------
view_assays <- list("main" = assay,
"juxta" = assay,
"para" = assay)
# Define features of each view ------------
view_features <- list("main" = useful_features,
"juxta" = useful_features,
"para" = useful_features)
# Define spatial context of each view -----
view_types <- list("main" = "intra",
"juxta" = "juxta",
"para" = "para")
# Define additional parameters (l in case of paraview,
# n of neighbors in case of juxta) --------
view_params <- list("main" = NULL,
"juxta" = 2,
"para" = 5)
misty_out <- paste0(misty_out_alias,
out_label, "_", assay)
run_misty_seurat(visium.slide = slide,
view.assays = view_assays,
view.features = view_features,
view.types = view_types,
view.params = view_params,
spot.ids = NULL,
out.alias = misty_out)
return(misty_out)
}
slide <- readRDS(rds)
DefaultAssay(slide) <- 'predictions'
useful_features <- rownames(slide) ####也可以自我设定感兴趣的细胞类型
useful_features <- useful_features[! useful_features %in% "prolif"]
mout <- run_colocalization(slide = slide,
useful_features = useful_features,
out_label = slide_id,
assay = assay,
misty_out_alias = outdir)
misty_res_slide <- collect_results(mout)
plot_folder <- paste0(mout, "/plots")
system(paste0("mkdir ", plot_folder))
pdf(file = paste0(plot_folder, "/", slide_id, "_", "summary_plots.pdf"))
mistyR::plot_improvement_stats(misty_res_slide)
mistyR::plot_view_contributions(misty_res_slide)
mistyR::plot_interaction_heatmap(misty_res_slide, "intra", cutoff = 0)
mistyR::plot_interaction_communities(misty_res_slide, "intra", cutoff = 0.5)
mistyR::plot_interaction_heatmap(misty_res_slide, "juxta_2", cutoff = 0)
mistyR::plot_interaction_communities(misty_res_slide, "juxta_2", cutoff = 0.5)
mistyR::plot_interaction_heatmap(misty_res_slide, "para_5", cutoff = 0)
mistyR::plot_interaction_communities(misty_res_slide, "para_5", cutoff = 0.5)
dev.off()
source代码如下(misty_utilities.R)：
以下内容为付费内容，定价 ¥50.00点击修改
run_misty_seurat <- function(visium.slide,
# Seurat object with spatial transcriptomics data.
view.assays,
# Named list of assays for each view.
view.features = NULL,
# Named list of features/markers to use.
# Use all by default.
view.types,
# Named list of the type of view to construct
# from the assay.
view.params,
# Named list with parameters (NULL or value)
# for each view.
spot.ids = NULL,
# spot IDs to use. Use all by default.
out.alias = "results"
# folder name for output
) {
mistyR::clear_cache()
# Extracting geometry
geometry <- GetTissueCoordinates(visium.slide,
cols = c("row", "col"), scale = NULL
)
# Extracting data
view.data <- map(view.assays,
extract_seurat_data,
geometry = geometry,
visium.slide = visium.slide
)
# Constructing and running a workflow
build_misty_pipeline(
view.data = view.data,
view.features = view.features,
view.types = view.types,
view.params = view.params,
geometry = geometry,
spot.ids = spot.ids,
out.alias = out.alias
)
}
# Extracts data from an specific assay from a Seurat object
# and aligns the IDs to the geometry
extract_seurat_data <- function(visium.slide,
assay,
geometry) {
print(assay)
data <- GetAssayData(visium.slide, assay = assay) %>%
as.matrix() %>%
t() %>%
as_tibble(rownames = NA)
return(data %>% slice(match(rownames(.), rownames(geometry))))
}
# Filters data to contain only features of interest
filter_data_features <- function(data,
features) {
if (is.null(features)) features <- colnames(data)
return(data %>% rownames_to_column() %>%
select(rowname, all_of(features)) %>% rename_with(make.names) %>%
column_to_rownames())
}
# Builds views depending on the paramaters defined
create_default_views <- function(data,
view.type,
view.param,
view.name,
spot.ids,
geometry) {
mistyR::clear_cache()
view.data.init <- create_initial_view(data)
if (!(view.type %in% c("intra", "para", "juxta"))) {
view.type <- "intra"
}
if (view.type == "intra") {
data.red <- view.data.init[["intraview"]]$data %>%
rownames_to_column() %>%
filter(rowname %in% spot.ids) %>%
select(-rowname)
} else if (view.type == "para") {
view.data.tmp <- view.data.init %>%
add_paraview(geometry, l = view.param)
data.ix <- paste0("paraview.", view.param)
data.red <- view.data.tmp[[data.ix]]$data %>%
mutate(rowname = rownames(data)) %>%
filter(rowname %in% spot.ids) %>%
select(-rowname)
} else if (view.type == "juxta") {
view.data.tmp <- view.data.init %>%
add_juxtaview(
positions = geometry,
neighbor.thr = view.param
)
data.ix <- paste0("juxtaview.", view.param)
data.red <- view.data.tmp[[data.ix]]$data %>%
mutate(rowname = rownames(data)) %>%
filter(rowname %in% spot.ids) %>%
select(-rowname)
}
if (is.null(view.param) == TRUE) {
misty.view <- create_view(
paste0(view.name),
data.red
)
} else {
misty.view <- create_view(
paste0(view.name, "_", view.param),
data.red
)
}
return(misty.view)
}
# Builds automatic MISTy workflow and runs it
build_misty_pipeline <- function(view.data,
view.features,
view.types,
view.params,
geometry,
spot.ids = NULL,
out.alias = "default") {
# Adding all spots ids in case they are not defined
if (is.null(spot.ids)) {
spot.ids <- rownames(view.data[[1]])
}
# First filter the features from the data
view.data.filt <- map2(view.data, view.features, filter_data_features)
# Create initial view
views.main <- create_initial_view(view.data.filt[[1]] %>%
rownames_to_column() %>%
filter(rowname %in% spot.ids) %>%
select(-rowname))
# Create other views
view.names <- names(view.data.filt)
all.views <- pmap(list(
view.data.filt[-1],
view.types[-1],
view.params[-1],
view.names[-1]
),
create_default_views,
spot.ids = spot.ids,
geometry = geometry
)
pline.views <- add_views(
views.main,
unlist(all.views, recursive = FALSE)
)
# Run MISTy
run_misty(pline.views, out.alias, cached = FALSE)
}
#
# Bug in collecting results
#
collect_results_v2 <- function(folders){
samples <- R.utils::getAbsolutePath(folders)
message("\nCollecting improvements")
improvements <- samples %>% furrr::future_map_dfr(function(sample) {
performance <- readr::read_table2(paste0(sample, .Platform$file.sep,
"performance.txt"), na = c("", "NA", "NaN"), col_types = readr::cols()) %>%
dplyr::distinct()
performance %>% dplyr::mutate(sample = sample, gain.RMSE = 100 *
(.data$intra.RMSE - .data$multi.RMSE)/.data$intra.RMSE,
gain.R2 = 100 * (.data$multi.R2 - .data$intra.R2),
)
}, .progress = TRUE) %>% tidyr::pivot_longer(-c(.data$sample,
.data$target), names_to = "measure")
message("\nCollecting contributions")
contributions <- samples %>% furrr::future_map_dfr(function(sample) {
coefficients <- readr::read_table2(paste0(sample, .Platform$file.sep,
"coefficients.txt"), na = c("", "NA", "NaN"), col_types = readr::cols()) %>%
dplyr::distinct()
coefficients %>% dplyr::mutate(sample = sample, .after = "target") %>%
tidyr::pivot_longer(cols = -c(.data$sample, .data$target),
names_to = "view")
}, .progress = TRUE)
improvements.stats <- improvements %>% dplyr::filter(!stringr::str_starts(.data$measure,
"p\\.")) %>% dplyr::group_by(.data$target, .data$measure) %>%
dplyr::summarise(mean = mean(.data$value), sd = stats::sd(.data$value),
cv = .data$sd/.data$mean, .groups = "drop")
contributions.stats <- dplyr::inner_join((contributions %>%
dplyr::filter(!stringr::str_starts(.data$view, "p\\.") &
.data$view != "intercept") %>% dplyr::group_by(.data$target,
.data$view) %>% dplyr::summarise(mean = mean(.data$value),
.groups = "drop_last") %>% dplyr::mutate(fraction = abs(.data$mean)/sum(abs(.data$mean))) %>%
dplyr::ungroup()), (contributions %>% dplyr::filter(stringr::str_starts(.data$view,
"p\\.") & !stringr::str_detect(.data$view, "intercept")) %>%
dplyr::group_by(.data$target, .data$view) %>% dplyr::mutate(view = stringr::str_remove(.data$view,
"^p\\.")) %>% dplyr::summarise(p.mean = mean(.data$value),
p.sd = stats::sd(.data$value), .groups = "drop")), by = c("target",
"view"))
message("\nCollecting importances")
importances <- samples %>% furrr::future_map(function(sample) {
targets <- contributions.stats %>% dplyr::pull(.data$target) %>%
unique() %>% sort()
views <- contributions.stats %>% dplyr::pull(.data$view) %>%
unique()
maps <- views %>% furrr::future_map(function(view) {
all.importances <- targets %>% purrr::map(~readr::read_csv(paste0(sample,
.Platform$file.sep, "importances_", .x, "_",
view, ".txt"), col_types = readr::cols()) %>%
dplyr::distinct() %>% dplyr::rename(feature = target))
features <- all.importances %>% purrr::map(~.x$feature) %>%
unlist() %>% unique() %>% sort()
pvalues <- contributions %>% dplyr::filter(sample ==
!!sample, view == paste0("p.", !!view)) %>% dplyr::mutate(value = 1 -
.data$value)
all.importances %>% purrr::imap_dfc(~tibble::tibble(feature = features,
zero.imp = 0) %>% dplyr::left_join(.x, by = "feature") %>%
dplyr::arrange(.data$feature) %>% dplyr::mutate(imp = scale(.data$imp)[,
1], `:=`(!!targets[.y], .data$zero.imp + (.data$imp *
(pvalues %>% dplyr::filter(target == targets[.y]) %>%
dplyr::pull(.data$value))))) %>% dplyr::select(targets[.y])) %>%
dplyr::mutate(Predictor = features)
}) %>% `names<-`(views)
}, .progress = TRUE) %>% `names<-`(samples)
message("\nAggregating")
importances.aggregated <- importances %>% purrr::reduce(function(acc,
l) {
acc %>% purrr::map2(l, ~(((.x %>% dplyr::select(-.data$Predictor)) +
(.y %>% dplyr::select(-.data$Predictor))) %>% dplyr::mutate(Predictor = .x %>%
dplyr::pull(.data$Predictor))))
}) %>% purrr::map(~.x %>% dplyr::mutate_if(is.numeric, ~./length(samples)))
return(list(improvements = improvements, improvements.stats = improvements.stats,
contributions = contributions, contributions.stats = contributions.stats,
importances = importances, importances.aggregated = importances.aggregated))
}
最后放上所有代码的链接
链接：https://pan.baidu.com/s/1eodGG5K3CsJMs0qdaWun7A?pwd=ZyfY
提取码：ZyfY
--来自百度网盘超级会员V1的分享
扫一扫
1766
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
