支持中心

简体中文

English

简体中文

English

Appearance

单细胞空间转录组: 空间邻域富集分析（基于 CellCharter） ​

作者: SeekGene
时长: 16 分钟
字数: 3.9k 字
更新: 2026-06-05
阅读: 0 次
空间转录组 Notebooks 分析指南 共定位分析

1. 模块简介 ​

本模块基于 CellCharter算法，用于对空间转录组数据进行空间邻域富集分析。

CellCharter 通过构建细胞间的空间邻近网络，评估不同细胞群之间的空间共定位模式。实现以下核心分析目标：

  • 单样本细胞类型共定位：评估组织内两个细胞群（或状态）在空间上的邻近性是否显著高于随机预期，从而揭示细胞间的相互吸引或排斥关系。
  • 多样本差异邻域富集分析：比较不同生物学条件（如健康与疾病、不同组别）下细胞空间相互作用的显著变化。

💡 Note
本流程输入为 Seurat/Scanpy 兼容格式数据，分析结果将生成展示细胞空间互作关系的热图及相关表格。

算法原理简述：

  1. 构建网络：根据空间坐标构建节点连接网络（基于 Delaunay 三角剖分等）。2. 计算富集分数 (NE)：计算实际连接数 (Observed) 与基于节点度数的随机预期连接数 (Expected) 的比值。3. 差异分析：通过置换检验 (Permutation) 比较两组间 NE 分数的差异显著性。

2. 输入文件准备 ​

本模块需要提供包含表达矩阵和空间坐标的文件以及对应的 Metadata。

文件结构示例：

text
filtered_feature_bc_matrix/          # rds 转化矩阵 输出目录
├── matrix.mtx.gz                    # 表达矩阵
├── barcodes.tsv.gz                  # 细胞条码
├── features.tsv.gz                  # 基因特征
└── cell_locations.tsv               # 空间坐标文件

如何从 Seurat 对象导出 CellCharter 所需文件：

如果您已有 Seurat 对象（.rds），可以使用以下 R 代码片段将表达矩阵和空间坐标导出为上方所示的文件夹结构：

R
library(Seurat)
library(dplyr)
library(tibble)
library(readr)
library(Matrix)

# 1. 设定输入路径
rds_path <- "your_input.rds"
meta_path <- "your_meta.txt" 

# 2. 读取数据并整合 metadata
obj <- readRDS(rds_path)
meta <- read.delim(meta_path)
rownames(meta) <- meta$barcode
obj@meta.data <- meta

# 3. 创建输出目录
dir.create('filtered_feature_bc_matrix', showWarnings = FALSE)

# 4. 提取空间坐标并保存
spatial_data <- Embeddings(obj, reduction = 'spatial')
spatial_data <- as.data.frame(spatial_data)
colnames(spatial_data) <- c('x', 'y')
spatial_data <- rownames_to_column(spatial_data, "Barcode")
write_delim(spatial_data, 'filtered_feature_bc_matrix/cell_locations.tsv', delim = '\t')

# 5. 提取表达矩阵并保存 (MTX 格式)
writeMM(GetAssayData(obj, slot="counts"), file = 'filtered_feature_bc_matrix/matrix.mtx')

# 6. 保存基因和条码特征
Gene <- data.frame(GeneID = rownames(obj), Gene = rownames(obj))
write.table(Gene, 'filtered_feature_bc_matrix/genes.tsv', row.names = FALSE, col.names = FALSE, sep = '\t', quote = FALSE)

barcode <- data.frame(Barcode = colnames(obj))
write.table(barcode, 'filtered_feature_bc_matrix/barcodes.tsv', row.names = FALSE, col.names = FALSE, sep = '\t', quote = FALSE)

3. 数据加载与预处理 ​

python
# --- 环境准备 ---
import squidpy as sq
import cellcharter as cc
import scanpy as sc
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from itertools import combinations
import matplotlib.image as mpimg
import os
import math
output
/PROJ2/FLOAT/jinwen/apps/miniconda3/envs/cellcharter-env/lib/python3.9/site-packages/numba/core/decorators.py:246: RuntimeWarning: nopython is set for njit and is ignored
warnings.warn('nopython is set for njit and is ignored', RuntimeWarning)
/PROJ2/FLOAT/jinwen/apps/miniconda3/envs/cellcharter-env/lib/python3.9/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
from .autonotebook import tqdm as notebook_tqdm
python
# --- 输入参数配置 ---

## meta_path：外部 metadata 文件路径
meta_path = "./rds/meta.tsv"

## col_sam：Metadata 中代表样本名称的列名
col_sam = "Sample"

## sample_name：需要分析的样本名称，多个样本用逗号分隔
sample_name = "25020230_nao_CS_expression"

## col_celltype：Metadata 中代表细胞类型的列名
col_celltype = "RNA_snn_res.0.2"

## celltypes：需要关注的细胞类型列表，用逗号分隔
celltypes = "0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19"

## col_group：Metadata 中代表实验分组的列名（用于多样本差异分析）
col_group = ""
python
# --- 参数解析 ---
sample_name=sample_name.strip().split(",")
celltypes=celltypes.strip().split(",")
python
# --- 数据加载与坐标合并 ---
adata = sc.read_10x_mtx("filtered_feature_bc_matrix/",cache=False) 

## 添加坐标信息
spatial_df = pd.read_csv("filtered_feature_bc_matrix/cell_locations.tsv",index_col=0,sep='\t')
spatial_df1 = spatial_df[spatial_df.index.isin(adata.obs_names)]
adata.obsm["spatial"] = spatial_df1.values
meta_df = pd.read_csv(meta_path,index_col=0,sep='\t')
meta_df = meta_df.reindex(adata.obs_names)
adata.obs = pd.concat([adata.obs, meta_df], axis=1)
python
# --- 数据格式转换 ---
adata.obs[col_sam] = adata.obs[col_sam].astype('category')
adata.obs[col_celltype] = adata.obs[col_celltype].astype('category')

3.1 构建空间邻域网络 ​

在进行邻域富集分析之前，首先需要基于物理空间坐标将细胞连接成一个“空间邻居网络图”。

三种常用的空间网络构建策略： 本工具提供了三种构建细胞邻居网络的方法，用户可根据数据特征（如细胞密度是否均匀）灵活选择：

  1. Delaunay 三角剖分 (Delaunay Triangulation)：
    • 设置参数：delaunay=True
    • 原理：在空间点之间连接形成连续的三角形网格，共享三角形边的细胞即被视为直接邻居。
    • 优势：对组织中细胞密度的变化具有极强的鲁棒性，无需预先猜测邻居数量或距离阈值。
  2. K-最近邻 (k-Nearest Neighbors, k-NN)：
    • 设置参数：n_neighs=K (如 n_neighs=6)
    • 原理：强制为每个细胞寻找空间距离最近的 K 个细胞作为邻居。
    • 优势：网络结构规整，每个节点的度数固定。
    • 劣势：如果组织中存在细胞稀疏区域，可能会跨越很大的空白物理空间强行连接不相关的细胞。
  3. 基于固定半径 (Radius-based)：
    • 设置参数：radius=R (如 radius=50)
    • 原理：以每个细胞为圆心，半径 R 范围内的所有细胞均被视为邻居。
    • 优势：严格基于物理距离阈值。
    • 劣势：如果组织细胞密度极不均匀，密集区的细胞会有过多邻居，而稀疏区的细胞可能没有邻居（成为孤立点）。

其他关键参数说明：

  • coord_type='generic'：表示使用通用的 2D 连续坐标（适用于 SeekSpace、Xenium 等），而不是像 Visium 那样固定的六边形/正方形网格。
  • library_key=col_sam：这是多样本分析时极其关键的参数。它指示算法必须在每个样本（切片）内部独立构建网络，严格防止将属于不同切片的细胞错误地连成邻居。
python
# --- 构建空间邻域网络 ---
# 基于空间坐标构建空间邻居图（默认使用 Delaunay 三角剖分）
sq.gr.spatial_neighbors(adata, library_key=col_sam, coord_type='generic', delaunay=True)

4. 单样本细胞类型共定位 ​

本部分主要用于评估单个组织切片内，两个细胞群（或细胞状态）之间的空间邻近性是否显著高于或低于随机分布的预期。

通过计算“邻域富集分数 (Neighborhood Enrichment Score)”，我们可以定量地描述不同细胞类型之间的空间互作模式：

  • 相互吸引 (Enrichment)：某些细胞类型倾向于在空间上聚集分布，可能暗示着它们之间存在活跃的细胞通讯或协同执行某种特定的微环境功能。
  • 相互排斥 (Depletion)：某些细胞类型在空间上刻意保持距离，可能反映了生态位的排他性或不同的发育起源。

图注说明：下图展示单样本细胞类型共定位邻域富集热图。

  • 行：源细胞类型 Ci；列：目标细胞类型 Cj。
  • 颜色 (正数)：表示细胞群 Ci 与 Cj 在空间上倾向于富集，即更有可能相邻。
  • 颜色 (负数)：表示细胞群 Ci 与 Cj 在空间上倾向于排斥，即更有可能分离。
  • 颜色 (零)：表示连接数与随机预期一致。
  • 判读：颜色越红/越深代表两类细胞在空间上越倾向于共定位聚集；颜色越蓝/越冷代表在空间上相互排斥。
python
# --- 数据格式转换 ---

for i in range(len(sample_name)):
    adata1 = adata[adata.obs[col_sam] == sample_name[i]]
    cc.gr.nhood_enrichment(
        adata1,
        cluster_key=col_celltype,
    )
    adata1 = adata1[adata1.obs[col_celltype].isin(celltypes)]
    adata1.obs[col_celltype] = adata1.obs[col_celltype].astype('category')
    #adata1.uns["FCN1_LGR5_2_nhood_enrichment"]["enrichment"] = adata1.uns["FCN1_LGR5_2_nhood_enrichment"]["enrichment"].drop(index="Others", columns="Others")

    adata1.uns[col_celltype+"_nhood_enrichment"]["enrichment"].columns = adata1.uns[col_celltype+"_nhood_enrichment"]["enrichment"].columns.astype(str)
    adata1.uns[col_celltype+"_nhood_enrichment"]["enrichment"].index = adata1.uns[col_celltype+"_nhood_enrichment"]["enrichment"].index.astype(str)
    adata1.uns[col_celltype+"_nhood_enrichment"]["enrichment"] = adata1.uns[col_celltype+"_nhood_enrichment"]["enrichment"].loc[celltypes,celltypes]
    cc.pl.nhood_enrichment(
        adata1,
        cluster_key=col_celltype,
        annotate=True,
        vmin=-0.1,
        vmax=0.1,
        figsize=(5,5),
        fontsize=7,
        save=sample_name[i]+"_"+col_celltype+"_enrichment_heatmap.png"
    )
output
/PROJ2/FLOAT/jinwen/apps/miniconda3/envs/cellcharter-env/lib/python3.9/site-packages/cellcharter/gr/_nhood.py:236: ImplicitModificationWarning: Trying to modify attribute \`._uns\` of view, initializing view as actual.
adata.uns[f"{cluster_key}_nhood_enrichment"] = result
/tmp/ipykernel_1483/3764853794.py:10: ImplicitModificationWarning: Trying to modify attribute \`.obs\` of view, initializing view as actual.
adata1.obs[col_celltype] = adata1.obs[col_celltype].astype('category')

5. 多样本差异邻域富集分析 ​

当研究设计包含多个生物学重复或不同处理条件（如健康 vs. 疾病、给药前 vs. 给药后）时，我们需要评估这些特定条件是否导致了细胞空间互作模式的显著改变。

本部分通过统计学置换检验（Permutation Test），比较两组样本之间的邻域富集分数差异：

  • 识别互作变化：找出在特定疾病状态或治疗干预下，哪些细胞群体之间的空间接触被特异性地增强或削弱了。
  • 揭示微环境重塑：这种空间互作的动态变化往往是组织微环境重塑的直接证据，能够为疾病的发生机制或药物的响应机制提供关键线索。

图注说明：下图展示多样本差异邻域富集分析热图。

  • 行：源细胞类型 Ci；列：目标细胞类型 Cj。
  • 颜色 (正数)：表示在实验组 (vs 前)中，细胞群 Ci 与 Cj 在空间上倾向于富集（相互靠近）。
  • 颜色 (负数)：表示在实验组 (vs 前)中，细胞群 Ci 与 Cj 在空间上倾向于排斥（相互远离）。
  • 颜色 (零)：无明显差异。
  • 判读：用于发现特定疾病或处理条件下，特异性增强或减弱的细胞空间互作模式。
python
# --- 计算多样本差异邻域富集 ---
if len(sample_name) > 1:
    cc.gr.diff_nhood_enrichment(
        adata,
        cluster_key=col_celltype,
        condition_key=col_sam,
        library_key=col_group,
        pvalues=True,
        n_jobs=15,
        n_perms=100
    )
python
# --- 绘制差异邻域富集热图 ---
if len(sample_name) > 1:
    combinations_list = list(combinations(sample_name, 2))
    for i in range(len(combinations_list)):
        cc.pl.diff_nhood_enrichment(
            adata,
            cluster_key=col_celltype,
            condition_key=col_group,
            condition_groups=list(combinations_list[i]),
            annotate=True,
            figsize=(5,5),
            #significance=0.05,
            fontsize=5,
            save=list(combinations_list[i])[0]+"vs"+list(combinations_list[i])[1]+"_"+col_celltype+"_enrichment_heatmap.png"
        )

6. 细胞空间分布密度分析 ​

基于物理距离量化不同细胞类型相对于目标细胞的空间分布密度情况。

图注说明：下图展示细胞类型共定位密度图。

  • 横坐标 (Distance to Target)：表示某一类细胞到目标细胞（Target Cell）的空间物理距离。
  • 纵坐标 (Density)：表示细胞在不同距离区间的分布密度。
  • 颜色分组 (Phenotype)：不同的颜色曲线代表不同的细胞类型（Phenotype）。
  • 判读：曲线的峰值越靠左（距离越近），说明该细胞类型在空间上越倾向于与目标细胞共定位（聚集）；曲线越平缓或峰值越靠右，说明该细胞类型与目标细胞在空间上分布较远或呈随机分布。
python
import pandas as pd
import numpy as np
from scipy.spatial import cKDTree

def find_nearest_distance(df, 
                          id_col='Cell ID', 
                          pheno_col='Phenotype', 
                          x_col='Cell X Position', 
                          y_col='Cell Y Position'):

    results = {}
    # 获取数据中所有唯一的表型
    phenotypes = df[pheno_col].dropna().unique()
    
    # 提取所有细胞的坐标和ID
    coords_all = df[[x_col, y_col]].values
    cell_ids_all = df[id_col].values
    
    for pheno in phenotypes:
        # 筛选出特定表型的细胞
        mask_p = (df[pheno_col] == pheno)
        df_p = df[mask_p].reset_index(drop=True)
        
        dist_col_name = f'Distance to {pheno}'
        id_col_name = f'Cell ID {pheno}'
        
        n_p = len(df_p)
        
        # 如果该表型没有细胞，返回 NaN
        if n_p == 0:
            results[dist_col_name] = np.nan
            results[id_col_name] = np.nan
            continue
            
        coords_p = df_p[[x_col, y_col]].values
        ids_p = df_p[id_col].values
        
        # 构建 KDTree 以进行快速的空间近邻搜索
        tree = cKDTree(coords_p)
        
        # 如果该表型只有 1 个细胞，我们只查询最近的 1 个；否则查询最近的 2 个（为了排除自身）
        k_query = 2 if n_p > 1 else 1
        distances, indices = tree.query(coords_all, k=k_query)
        
        final_dists = np.zeros(len(df))
        final_ids = np.empty(len(df), dtype=object)
        
        if k_query == 1:
            for i in range(len(df)):
                # 如果找到的最近细胞就是它自己（ID相同），说明该表型没有"其他"细胞了
                if cell_ids_all[i] == ids_p[indices[i]]:
                    final_dists[i] = np.nan
                    final_ids[i] = pd.NA
                else:
                    final_dists[i] = distances[i]
                    final_ids[i] = ids_p[indices[i]]
        else:
            for i in range(len(df)):
                # indices[i, 0] 是最近的，indices[i, 1] 是第二近的
                # 如果最近的细胞是它自己，就取第二近的细胞
                if cell_ids_all[i] == ids_p[indices[i, 0]]:
                    final_dists[i] = distances[i, 1]
                    final_ids[i] = ids_p[indices[i, 1]]
                else:
                    final_dists[i] = distances[i, 0]
                    final_ids[i] = ids_p[indices[i, 0]]
                    
        results[dist_col_name] = final_dists
        results[id_col_name] = final_ids
        
    # 返回 DataFrame 格式，索引与原数据对齐
    return pd.DataFrame(results, index=df.index)
python
adata.obs[["Cell X Position", "Cell Y Position"]] = adata.obsm["spatial"]
position = adata.obs[["Cell X Position","Cell Y Position",col_celltype,col_sam]]
position = position.reset_index(names='Cell ID')
position.rename(columns={col_celltype: 'Phenotype'}, inplace=True)
position[["Cell X Position","Cell Y Position"]] = position[["Cell X Position","Cell Y Position"]]*0.2653
python
custom_40_colors = [
    '#1f77b4', '#aec7e8', '#ff7f0e', '#ffbb78', '#2ca02c', '#98df8a', '#d62728', '#ff9896',
    '#9467bd', '#c5b0d5', '#8c564b', '#c49c94', '#e377c2', '#f7b6d2', '#7f7f7f', '#c7c7c7',
    '#bcbd22', '#dbdb8d', '#17becf', '#9edae5', '#393b79', '#5254a3', '#6b6ecf', '#9c9ede',
    '#637939', '#8ca252', '#b5cf6b', '#cedb9c', '#8c6d31', '#bd9e39', '#e7ba52', '#e7cb94',
    '#843c39', '#ad494a', '#d6616b', '#e7969c', '#7b4173', '#a55194', '#ce6dbd', '#de9ed6'
]

# === 布局设计：一行最多排列 3 个样本的密度图 ===
max_cols = 3
num_samples = len(sample_name)

# 自动计算实际需要的列数和行数
ncols = max_cols if num_samples >= max_cols else num_samples
nrows = math.ceil(num_samples / ncols) if num_samples > 0 else 1

# 设置整体主题
sns.set_theme(style="whitegrid", rc={"axes.edgecolor": "black"})

# 动态生成画布大小：宽度为每个子图 12（留出充足的图例空间），高度为每个子图 6
fig, axes = plt.subplots(nrows=nrows, ncols=ncols, figsize=(12 * ncols, 6 * nrows))

# 统一转换为一维列表，方便在循环中按索引取用 ax
if isinstance(axes, np.ndarray):
    axes_list = axes.flatten()
else:
    axes_list = [axes]

for i in range(num_samples):
    ax = axes_list[i]
    
    # 提取当前样本数据
    sub_position = position.loc[position[col_sam] == sample_name[i]]
    distances = find_nearest_distance(sub_position)
    csd_with_distance = pd.concat([sub_position, distances], axis=1)
    
    # 绘制密度图，指定画在当前的 ax 上
    sns.kdeplot(
        data=csd_with_distance,
        x='Distance to 1',
        hue='Phenotype',
        linewidth=2,
        common_norm=False,
        clip=(0, None),
        palette=custom_40_colors,
        ax=ax
    )

    ax.set_xlim(left=0)
    ax.set_ylabel('Density', fontsize=12)
    ax.set_title(f'Sample: {sample_name[i]}', fontsize=14, fontweight='bold')

    # 使用 seaborn 专用的 move_legend 将每个子图的图例移到各自的右侧
    if ax.get_legend() is not None:
        sns.move_legend(
            ax, "upper left",
            bbox_to_anchor=(1.02, 1), 
            ncol=2, 
            title="Phenotype", 
            frameon=False
        )

# 如果子图网格数大于样本数（例如一行3个但只有4个样本），则删除多余的空白子图
for j in range(num_samples, len(axes_list)):
    fig.delaxes(axes_list[j])

# 自动调整子图间距，防止图例与相邻图重叠
plt.tight_layout()
plt.show()
output
/tmp/ipykernel_1483/54170521.py:44: UserWarning: The palette list has more values (40) than needed (20), which may not be intended.
sns.kdeplot(
python
0 条评论·0 条回复