内容复习----10X visium多样本的共定位分析

作者，Evil Genius

假期总是很快，老家出来到当牛马的城市，总有一种恍如隔世的感觉。其实还房贷车贷也挺好的，比饥荒强。

今天本来想复习两个内容的，一个是hotspot。

另外一个是多样本的共定位分析

但是实在是不想思考了，先复习共定位分析吧。

首先10X visium的基础分析大家自己做，保存好h5ad文件，并且做好单细胞空间联合分析，拿到单细胞空间联合后的细胞矩阵文件。

接下来先要分析单样本的共定位结果

首先定义一下范围intra-view是点内（same-spot）， juxta-view为直径100um的范围， para-view为直径200um的范围。

拿到共定位分析结果文件

library(zellkonverter)
library(mistyR)
library(tidyverse)
library(SingleCellExperiment)
library(distances)

path <- '/Results/mistyR_Analysis/Objects/'

files <- grep('h5ad', list.files(path, full.names = T), value = T)

file.names <- c(样本名称列表)

计算共定位文件

counter = 1
importances <- list()
for (file in files) {
  print(file.names[counter])
  # Load AnnData Object
  sce <- readH5AD(file)

  # Get Spatial Coordinates
  coords <- as_tibble(reducedDim(sce, 'spatial'))
  colnames(coords) <- c('x', 'y')

  # Calculate radius
  geom_dist <- as.matrix(distances(as.data.frame(coords)))
  dist_nn <- apply(geom_dist, 1, function(x) (sort(x)[2]))
  juxta_radius <- ceiling(mean(dist_nn + sd(dist_nn)))

  dist_en <- apply(geom_dist, 1, function(x) (sort(x)[8]))
  paraview_radius <- ceiling(mean(dist_en + sd(dist_en)))

  # Get c2l abundances
  c2l <- as_tibble(reducedDim(sce, 'c2l'))

联合后的细胞矩阵列名，大家自己修改一下，不要c2l自带的那个列名（q05开头的）。

  colnames(c2l) <- c('Adip', 'ArtEC', 'B_cells', 'CM', 'CapEC', 'EndoEC', 'Fibroblasts', 'Fibro_activ', 'LymphEC', 'MP', 'Epi_cells', 'Ccr2_MP', 'Pericytes', 'SMC', 'T_cells', 'VeinEC')

计算共定位

  # Create misty views
  heart_views <- create_initial_view(c2l) %>%
    add_juxtaview(coords, neighbor.thr = juxta_radius, verbose = T) %>%   # Juxta View is the Hexamers
    add_paraview(coords, l = paraview_radius, verbose = T, zoi = juxta_radius)  # Para View - Juxtaview

  # Run Misty
  #new.path <- paste0(path,'tmp_misty_v2/', '_', file.names[counter] )
  new.path <- paste0(path, 'tmp_misty_v3/', '_', file.names[counter])
  dir.create(new.path)

  run_misty(heart_views, new.path)
  misty_results <- collect_results(new.path)

  df <- misty_results$importances.aggregated
  df['sample'] <- file.names[counter]
  importances[[counter]] <- df
  counter = counter + 1
}

df <- do.call(rbind, importances)
write_csv(df, '/mistyR_Analysis/importances_samples_c2l_v3.csv')

# Get improvements stats
counter <- 1
r2.values <- list()
for (file in file.names) {
  print(file)
  #new.path <- paste0(path,'tmp_misty_v2/', '_',file )
  new.path <- paste0(path, 'tmp_misty_v3/', '_', file)

  misty_results <- collect_results(new.path)

  df.r2 <- misty_results$improvements.stats
  df.r2['sample'] <- file
  r2.values[[counter]] <- df.r2
  counter <- counter + 1
}

df.r2 <- do.call(rbind, r2.values)
write_csv(df.r2, '/mistyR_Analysis/importances_samples_c2l_r2vals_v3.csv')

拿到数据之后，我们进行绘图,python绘图

import os
from tqdm import tqdm
import itertools

import anndata as ad
import scanpy as sc
import bbknn as bkn
import decoupler as dc
import pandas as pd
import numpy as np
from scipy.cluster.hierarchy import linkage, dendrogram
from scipy.stats import mannwhitneyu, shapiro
import statsmodels.api as sm
from statsmodels.formula.api import ols
from statsmodels.stats.anova import AnovaRM

import matplotlib.pyplot as plt
from matplotlib.patches import PathPatch
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
import matplotlib.lines as mlines
import seaborn as sns

import davidrPlotting
import davidrUtility
import davidrExperimental
import davidrScanpy

misty_path = os.path.join('/mistyR_Analysis/')

首先第一张图