10X空间转录组之空间共定位的一点思考和代码

作者，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)：