作者，追風(fēng)少年i

新的一周，新的開始，接上一篇，10X空間轉(zhuǎn)錄組數(shù)據(jù)分析之細(xì)胞的空間依賴性，這一篇我們來分析空間細(xì)胞類型和空間通路分布的依賴性，如下圖：

圖片.png

我們的目的就是研究空間細(xì)胞類型分布是否與生物學(xué)通路的分布存在強(qiáng)相關(guān)性。

并且我們要用代碼來實(shí)現(xiàn)一下

library(tidyverse)
library(Seurat)
library(mistyR)
source("./misty_utilities.R")

source的misty_utilities.R內(nèi)容在文章10X空間轉(zhuǎn)錄組數(shù)據(jù)分析之細(xì)胞的空間依賴性最下面

我們需要準(zhǔn)備的內(nèi)容

• 單細(xì)胞空間聯(lián)合分析的矩陣（Seurat、cell2location均可）
• 空間每個spot的通路分析結(jié)果（富集分析，或者其他方式的打分矩陣）

定義共定位函數(shù)

# Pipeline definition:
run_colocalization <- function(slide, 
                               assay, 
                               useful_features,
                               useful_features_ct,
                               out_label, 
                               misty_out_alias = "./results/tissue_structure/misty/pathway_map/pm_") {
  
  # Define assay of each view ---------------
  view_assays <- list("main" = assay,
                      "juxta" = assay,
                      "para" = assay,
                      "intra_ct" = "c2l",
                      "para_ct" = "c2l")
  # Define features of each view ------------
  view_features <- list("main" = useful_features, 
                        "juxta" = useful_features,
                        "para" = useful_features,
                        "intra_ct" = useful_features_ct,
                        "para_ct" = useful_features_ct)
  # Define spatial context of each view -----
  view_types <- list("main" = "intra", 
                     "juxta" = "juxta",
                     "para" = "para",
                     "intra_ct" = "intra",
                     "para_ct" = "para")
  # Define additional parameters (l in case of paraview,
  # n of neighbors in case of juxta) --------
  view_params <- list("main" = NULL, 
                      "juxta" = 5,
                      "para" = 15,
                      "intra_ct" = NULL,
                      "para_ct" = 15)
  
  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)
}

讀取空間的rds文件和spot打分注釋文件

slide = readRDS(spatialrds)
en_anno = read.csv(enrichment_spatial.csv)  ###這里大家需要自己分析

這里的通路富集分析推薦大家使用R包progeny，PROGENy is resource that leverages a large compendium of publicly available signaling perturbation experiments to yield a common core of pathway responsive genes for human and mouse. These, coupled with any statistical method, can be used to infer pathway activities from bulk or single-cell transcriptomics.

library(progeny)

if(species == "mouse"){
    model <- progeny::getModel(organism = "Mouse", top = top)
}else if(species == "human"){
    model <- progeny::getModel(organism = "Human", top = top)

富集分析推薦使用AUC或者GSVA，這里就不做了，直接拿著分析結(jié)果往下

progeny_scores <- scale(t(progeny_scores))  ####分?jǐn)?shù)的標(biāo)準(zhǔn)化操作

slide[['progeny']] <- CreateAssayObject(counts = t(progeny_scores))

其二是單細(xì)胞空間聯(lián)合分析的矩陣

anno = read.csv(sp_sc.csv,header = T,row.names = 1)

slide[['c2l']] <- CreateAssayObject(counts = t(anno))

數(shù)據(jù)分析

assay_label <- "progeny"

assay <- assay_label
DefaultAssay(slide) <- assay
useful_features <- rownames(slide)
useful_features <- useful_features[! useful_features %in% c("TNFa")]  ###通路特征
  
useful_features_ct <- rownames(GetAssayData(slide, assay = "c2l"))  ###細(xì)胞類型特征
useful_features_ct <- useful_features_ct[! useful_features_ct %in% "prolif"]

mout <- run_colocalization(slide = slide,
                             useful_features = useful_features,
                             useful_features_ct = useful_features_ct,
                             out_label = slide_id,
                             assay = assay,
                             misty_out_alias = "./results/tissue_structure/misty/pathway_map/pm_")
  
misty_res_slide <- collect_results(mout)

數(shù)據(jù)分析結(jié)束后進(jìn)行可視化

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_5", cutoff = 0)
mistyR::plot_interaction_communities(misty_res_slide, "juxta_5", cutoff = 0.5)
  
mistyR::plot_interaction_heatmap(misty_res_slide, "para_15", cutoff = 0)
mistyR::plot_interaction_communities(misty_res_slide, "para_15", cutoff = 0.5)
  
mistyR::plot_interaction_heatmap(misty_res_slide, "intra_ct", cutoff = 0)
  
mistyR::plot_interaction_heatmap(misty_res_slide, "para_ct_15", cutoff = 0)
  
dev.off()

最后分析得到細(xì)胞類型分布與空間分布的依賴性熱圖

圖片.png

好了，已經(jīng)分享給大家了，生活很好，有你更好，百度文庫出現(xiàn)了大量抄襲我的文章，對此我深表無奈，我寫的文章，別人掛上去賺錢，抄襲可恥，掛到百度文庫的人更可恥

這里我也感覺到很有難度了。