空转细胞通讯分析之COMMOT|江阴雨辰互联

空转细胞通讯分析之COMMOT

空间转录组以及单细胞转录组技术，为解析细胞-细胞通讯（Cell–Cell Communication, CCC）提供了前所未有的机会。然而，如何在重建细胞通讯网络的过程中充分整合细胞间的空间位置信息和复杂的生化过程，仍然是一个重大挑战。

COMMOT（COMMunication analysis by Optimal Transport）是一种创新的分析方法，旨在在空间转录组数据中推断细胞通讯关系。它基于集体最优传输理论（collective optimal transport），能够同时考虑多个配体-受体对之间的相互竞争关系以及细胞之间的空间距离，实现更为精确和生物学上合理的通讯建模。该算法发表于2023年的Nature Methods，题为《Screening cell–cell communication in spatial transcriptomics via collective optimal transport》。

COMMOT 的主要功能包括：

支持多种类型空间数据（如原始空间转录组或配对成像数据推断的空间单细胞转录组）；
模拟多个配体–受体对在空间中的传递与分布；
推断信号传导的空间方向性；
借助集成树模型识别下游受信号调控的关键基因；
提供全面的可视化工具，包括通讯热图、信号流向图和基因调控网络图等。

下面我们学习一下官网的参考文档（）：

首先是环境的部署：

代码语言：javascript代码运行次数：0运行复制

#构建分析环境commot
mamba create -n commot python=3.8 -y
conda activate commot
pip install commot
#mamba install -c conda-forge pygraphviz
#pip install get-version==2.0.4

mamba install -y nb_conda_kernels ipykernel
python -m ipykernel install --user --name commot --display-name commot

可选：

代码语言：javascript代码运行次数：0运行复制

#安装官网教程，在下游分析DE genes时会用到tradeseq
mamba install rpy2
mamba install r-base==3.6.3
pip install commot[tradeSeq]

代码语言：javascript代码运行次数：0运行复制

#安装完成后，进入R程序,在里面安装tradeSeq
if (!require("BiocManager", quietly = TRUE))
   install.packages("BiocManager")
BiocManager::install("tradeSeq")

参考文档总共分为两大部分：

Basic usage
以及Visium data example (mouse brain)

一. Basic usage

代码语言：javascript代码运行次数：0运行复制

import commot as ct
import scanpy as sc
import pandas as pd
import numpy as np

import matplotlib.pyplot as plt

Step1. 加载示例数据

代码语言：javascript代码运行次数：0运行复制

adata = sc.datasets.visium_sge(sample_id='V1_Mouse_Brain_Sagittal_Posterior')
adata.var_names_make_unique()

Step2. 数据预处理

该方法假定使用的数值为非负数，并反映信号分子的丰度。

代码语言：javascript代码运行次数：0运行复制

sc.pp.normalize_total(adata, inplace=True)
sc.pp.log1p(adata)

Step3. 指定配体-受体对

这里，我们以仅包含三个 LR 对为例。

可以用同样的方式指定用户自定义的 LR 数据库，或者也可以通过函数 commot.pp.ligand_receptor_database() 获取内置的 LR 数据库。

例如，df_ligrec = ct.pp.ligand_receptor_database(database='CellChat', species='mouse')。

代码语言：javascript代码运行次数：0运行复制

LR=np.array([['Tgfb1', 'Tgfbr1_Tgfbr2', 'Tgfb_pathway'],
    ['Tgfb2', 'Tgfbr1_Tgfbr2', 'Tgfb_pathway'],
    ['Tgfb3', 'Tgfbr1_Tgfbr2', 'Tgfb_pathway'],
    ['Fgf1', 'Fgfr1', 'FGF_pathway'],
    ['Fgf1', 'Fgfr2', 'FGF_pathway']],dtype=str)
df_ligrec = pd.DataFrame(data=LR)

Step4. 构建 CCC 网络

使用最优传输方法为满足 500 微米空间距离约束的配体-受体对构建 CCC 网络（距离超过 500 米的细胞之间禁止耦合）。例如，Tgfb1（配体）与 Tgfbr1_Tgfbr2（受体）的点对点矩阵存储在 adata.obsp['commot-user_database-Tgfb1-Tgfbr1_Tgfbr2'] 中。每个配体-受体对的总发送或接收信号分别存储在 adata.obsm['commot-user_database-sum-sender'] 和 adata.obsm['commot-user_database-sum-receiver'] 中。

代码语言：javascript代码运行次数：0运行复制

ct.tl.spatial_communication(adata,
    database_name='user_database', df_ligrec=df_ligrec, dis_thr=500, heteromeric=True, pathway_sum=True)

Step5. 可视化

代码语言：javascript代码运行次数：0运行复制

pts = adata.obsm['spatial']
s = adata.obsm['commot-user_database-sum-sender']['s-Fgf1-Fgfr1']
r = adata.obsm['commot-user_database-sum-receiver']['r-Fgf1-Fgfr1']
fig, ax = plt.subplots(1,2, figsize=(10,4))
ax[0].scatter(pts[:,0], pts[:,1], c=s, s=5, cmap='Blues')
ax[0].set_title('Sender')
ax[1].scatter(pts[:,0], pts[:,1], c=r, s=5, cmap='Reds')
ax[1].set_title('Receiver')

image-20250321214846249

可视化信号传导方向将信号传导方向插值为向量场。

代码语言：javascript代码运行次数：0运行复制

ct.tlmunication_direction(adata, database_name='user_database', lr_pair=('Fgf1','Fgfr1'), k=5)
ct.pl.plot_cell_communication(adata, database_name='user_database', lr_pair=('Fgf1','Fgfr1'), plot_method='grid', background_legend=True,
    scale=0.00003, ndsize=8, grid_density=0.4, summary='sender', background='image', clustering='leiden', cmap='Alphabet',
    normalize_v = True, normalize_v_quantile=0.995)

image-20250321214906136

或者，根据接收信号的水平绘制方向。

代码语言：javascript代码运行次数：0运行复制

ct.pl.plot_cell_communication(adata, database_name='user_database', lr_pair=('Fgf1','Fgfr1'), plot_method='grid', background_legend=True,
    scale=0.00003, ndsize=8, grid_density=0.4, summary='receiver', background='summary', clustering='leiden', cmap='Reds',
    normalize_v = True, normalize_v_quantile=0.995)

image-20250321214921482

COMMOT 分析的结果附加在 anndata 对象中：

**adata.uns['commot-user_database-info']**：所使用的配体-受体数据库
**adata.obsm['commot-user_database-sum-sender']**：每个配体-受体对的发送信号量以及各信号通路的总和
**adata.obsm['commot-user_database-sum-receiver']**：每个配体-受体对的接收信号量以及各信号通路的总和
**adata.obsm['commot_sender_vf-user_database-Fgf1-Fgfr1']**：每个点上与发送信号相关的插值信号传导方向
**adata.obsm['commot_receiver_vf-user_database-Fgf1-Fgfr1']**：每个点上与接收信号相关的插值信号传导方向
**adata.obsp['commot-user_database-Fgf1-Fgfr1']**：一个稀疏矩阵，其中的元素表示通过 Fgf1-Fgfr1 从一个点到另一个点的 CCC 分数

代码语言：javascript代码运行次数：0运行复制

adata

代码语言：javascript代码运行次数：0运行复制

 AnnData object with n_obs × n_vars = 3355 × 32285
    obs: 'in_tissue', 'array_row', 'array_col'
    var: 'gene_ids', 'feature_types', 'genome'
    uns: 'spatial', 'log1p', 'commot-user_database-info'
    obsm: 'spatial', 'commot-user_database-sum-sender', 'commot-user_database-sum-receiver', 'commot_sender_vf-user_database-Fgf1-Fgfr1', 'commot_receiver_vf-user_database-Fgf1-Fgfr1'
    obsp: 'commot-user_database-Fgf1-Fgfr2', 'commot-user_database-Fgf1-Fgfr1', 'commot-user_database-Tgfb2-Tgfbr1_Tgfbr2', 'commot-user_database-Tgfb3-Tgfbr1_Tgfbr2', 'commot-user_database-Tgfb1-Tgfbr1_Tgfbr2', 'commot-user_database-FGF_pathway', 'commot-user_database-Tgfb_pathway', 'commot-user_database-total-total'

二. Visium data example (mouse brain)

代码语言：javascript代码运行次数：0运行复制

import os
import gc
import ot
import pickle
import anndata
import scanpy as sc
import pandas as pd
import numpy as np
from scipy import sparse
from scipy.stats import spearmanr, pearsonr
from scipy.spatial import distance_matrix
import matplotlib.pyplot as plt

import commot as ct

adata = sc.datasets.visium_sge(sample_id='V1_Mouse_Brain_Sagittal_Posterior')

Step1. 数据预处理

这里需要合理反映分子丰度的非负数值。我们这里只进行最基础的数据预处理，包括归一化总计数和 log1p 变换。也可以进行更详细的预处理，例如回归细胞周期基因的影响。

代码语言：javascript代码运行次数：0运行复制

adata.var_names_make_unique()
adata.raw = adata
sc.pp.normalize_total(adata, inplace=True)
sc.pp.log1p(adata)
adata_dis500 = adata.copy()

Step2. 基本聚类

使用常见的 scRNA-seq 数据聚类方法对数据进行基本聚类。或者，可以使用专门用于空间转录组数据的工具，例如 SpaceFlow

代码语言：javascript代码运行次数：0运行复制

sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
adata = adata[:, adata.var.highly_variable]
sc.tl.pca(adata, svd_solver='arpack')
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)
sc.tl.umap(adata)
sc.tl.leiden(adata, resolution=0.4)
sc.pl.spatial(adata, color='leiden')

image-20250310221142588

Step3. 空间通信推断

这里我们使用 CellChatDB 配体-受体数据库，仅使用分泌信号的 LR 对。

代码语言：javascript代码运行次数：0运行复制

df_cellchat = ct.pp.ligand_receptor_database(species='mouse', signaling_type='Secreted Signaling', database='CellChat')#database=‘CellPhoneDB_v4.0’
#df_cellchat = ct.pp.ligand_receptor_database(species='human', signaling_type='Secreted Signaling', database='CellChat')#database=‘CellPhoneDB_v4.0’

print(df_cellchat.shape)
# (1209, 4)

接下来我们筛选 LR 对，仅保留在至少 5% 的点上同时表达配体和受体的对：

代码语言：javascript代码运行次数：0运行复制

df_cellchat_filtered = ct.pp.filter_lr_database(df_cellchat, adata_dis500, min_cell_pct=0.05)
print(df_cellchat_filtered.shape)
# (250, 4)
print(df_cellchat_filtered.head())
#       0              1     2                   3
#0  Tgfb1  Tgfbr1_Tgfbr2  TGFb  Secreted Signaling
#1  Tgfb2  Tgfbr1_Tgfbr2  TGFb  Secreted Signaling
#2  Tgfb3  Tgfbr1_Tgfbr2  TGFb  Secreted Signaling
#3  Tgfb1  Acvr1b_Tgfbr2  TGFb  Secreted Signaling
#4  Tgfb1  Acvr1c_Tgfbr2  TGFb  Secreted Signaling

现在对这 250 个配体-受体对进行空间通信推断，空间距离限制为 500。 CellChat 数据库考虑了多聚体单元。信号传导结果存储在 obsp 槽中的点对矩阵中。例如，通过 LR 对从点 i 到点 j 的信号传导得分可从adata_dis500.obsp['commot-cellchat-Wnt4-Fzd4_Lrp6'][i,j] 中获得。

代码语言：javascript代码运行次数：0运行复制

ct.tl.spatial_communication(adata_dis500,
    database_name='cellchat', df_ligrec=df_cellchat_filtered, dis_thr=500, heteromeric=True, pathway_sum=True)
adata_dis500.write("./adata.h5ad")

确定信号传导方向

确定信号传导通路的空间方向，例如 PSAP 通路。点上信号发送的插值方向和信号接收的插值方向分别存储在adata_dis500.obsm['commot_sender_vf-cellchat-PSAP'] 和 adata_dis500.obsm['commot_receiver_vf-cellchat-PSAP'] 中。

代码语言：javascript代码运行次数：0运行复制

ct.tlmunication_direction(adata_dis500, database_name='cellchat', pathway_name='PSAP', k=5)
ct.pl.plot_cell_communication(adata_dis500, database_name='cellchat', pathway_name='PSAP', plot_method='grid', background_legend=True,
    scale=0.00003, ndsize=8, grid_density=0.4, summary='sender', background='image', clustering='leiden', cmap='Alphabet',
    normalize_v = True, normalize_v_quantile=0.995)
plt.savefig('sample1.PSAPmot.png',bbox_inches = 'tight')

image-20250310221310119

或者，基于 leiden 聚类可视化PSAP 信号传导通路传导情况，结果存储在adata_dis500.uns['commot_cluster-leide¸ˇn-cellchat-PSAP'] 中。

代码语言：javascript代码运行次数：0运行复制

adata_dis500.obs['leiden'] = adata.obs['leiden']
ct.tl.cluster_communication(adata_dis500, database_name='cellchat', pathway_name='PSAP', clustering='leiden',
    n_permutations=100)

我们可以将显著的结果可视化为具有自动节点嵌入的网络。

代码语言：javascript代码运行次数：0运行复制

ct.pl.plot_cluster_communication_network(adata_dis500, uns_names=['commot_cluster-leiden-cellchat-PSAP'],
    nx_node_pos=None, nx_bg_pos=False, p_value_cutoff = 5e-2, filename='PSAP_cluster.pdf', nx_node_cmap='Light24')

image-20250310222201331

或者，使用空间节点嵌入。

代码语言：javascript代码运行次数：0运行复制

ct.tl.cluster_position(adata_dis500, clustering='leiden')
ct.pl.plot_cluster_communication_network(adata_dis500, uns_names=['commot_cluster-leiden-cellchat-PSAP'], clustering='leiden',
    nx_node_pos='cluster', nx_pos_idx=np.array([0, 1]), nx_bg_pos=True, nx_bg_ndsize=0.25, p_value_cutoff=5e-2,
    filename='PSAP_cluster_spatial.pdf', nx_node_cmap='Light24')

image-20250310222236291

Step4. 下游分析

4.1 识别信号传导差异表达基因

接下来，我们进一步探究相对于推断的信号传导活性可能存在差异表达的基因。 tradeSeq 将被用于对差异表达基因进行建模。此分析类似于寻找时间性差异表达基因，只不过将伪时间变量替换为接收信号的量。需要在 adata_dis500.layers['counts'] 中提供计数数据。

参考文献： Van den Berge, Koen, et al. “Trajectory-based differential expression analysis for single-cell sequencing data.” Nature communications 11.1 (2020): 1-13.

代码语言：javascript代码运行次数：0运行复制

adata_dis500 = sc.read_h5ad("./adata.h5ad")
adata = sc.datasets.visium_sge(sample_id='V1_Mouse_Brain_Sagittal_Posterior')
adata_dis500.layers['counts'] = adata.X

寻找与 PSAP 信号传导相关的差异表达基因。

代码语言：javascript代码运行次数：0运行复制

df_deg, df_yhat = ct.tlmunication_deg_detection(adata_dis500,
    database_name = 'cellchat', pathway_name='PSAP', summary = 'receiver')

代码语言：javascript代码运行次数：0运行复制

import pickle
deg_result = {"df_deg": df_deg, "df_yhat": df_yhat}
with open('./deg_PSAP.pkl', 'wb') as handle:
    pickle.dump(deg_result, handle, protocol=pickle.HIGHEST_PROTOCOL)

对下游基因进行聚类，并可视化相对于通过 PSAP 通路接收信号水平升高的表达趋势（横轴）。

代码语言：javascript代码运行次数：0运行复制

with open("./deg_PSAP.pkl", 'rb') as file:
    deg_result = pickle.load(file)
df_deg_clus, df_yhat_clus = ct.tlmunication_deg_clustering(df_deg, df_yhat, deg_clustering_res=0.4)
top_de_genes_PSAP = ct.pl.plot_communication_dependent_genes(df_deg_clus, df_yhat_clus, top_ngene_per_cluster=5,
    filename='./heatmap_deg_PSAP.pdf', font_scale=1.2, return_genes=True)

image-20250310221845550

绘制一些示例信号传导差异表达基因。

代码语言：javascript代码运行次数：0运行复制

X_sc = adata_dis500.obsm['spatial']
fig, ax = plt.subplots(1,3, figsize=(15,4))
colors = adata_dis500.obsm['commot-cellchat-sum-receiver']['r-PSAP'].values
idx = np.argsort(colors)
ax[0].scatter(X_sc[idx,0], X_sc[idx,1], c=colors[idx], cmap='coolwarm', s=10)
colors = adata_dis500[:,'Ctxn1'].X.toarray().flatten()
idx = np.argsort(colors)
ax[1].scatter(X_sc[idx,0], X_sc[idx,1], c=colors[idx], cmap='coolwarm', s=10)
colors = adata_dis500[:,'Gpr37'].X.toarray().flatten()
idx = np.argsort(colors)
ax[2].scatter(X_sc[idx,0], X_sc[idx,1], c=colors[idx], cmap='coolwarm', s=10)
ax[0].set_title('Amount of received signal')
ax[1].set_title('An example negative DE gene (Ctxn1)')
ax[2].set_title('An example positive DE gene (Gpr37)')

image-20250310221925615

4.2 进一步量化信号传导对差异表达基因的影响

上述差异表达分析显示了信号传导与下游基因表达之间的相关性。现在，我们进一步构建随机森林模型，以潜在的差异表达基因作为目标，将信号传导作为输入，通过特征重要性来量化信号传导的影响。

在每个模型中，从 top_de_genes_PSAP 中选取一个基因作为预测目标。通过 PSAP 通路接收信号的水平以及一组“背景基因”作为输入特征。与预测目标相关性最高的基因被用作“背景基因”。模型中接收信号的特征重要性反映了在考虑其他相关基因影响的情况下，该信号传导对目标基因的影响。

代码语言：javascript代码运行次数：0运行复制

df_impact_PSAP = ct.tlmunication_impact(adata_dis500, database_name='cellchat', pathway_name='PSAP',
    tree_combined=True, method='treebased_score', tree_ntrees=100, tree_repeat=100, tree_method='rf',
    ds_genes=top_de_genes_PSAP, bg_genes=500, normalize=True)

将影响评分可视化为热图。这里在 PSAP 中只有两个 LR 对，我们展示了前 30 个差异表达基因。

代码语言：javascript代码运行次数：0运行复制

ct.pl.plot_communication_impact(df_impact_PSAP, summary = 'receiver', top_ngene= 30, top_ncomm = 5, colormap='coolwarm',
    font_scale=1.2, linewidth=0, show_gene_names=True, show_comm_names=True, cluster_knn=2,
    filename = 'heatmap_impact_PSAP.pdf')

image-20250310221938600

代码语言：javascript代码运行次数：0运行复制

df_impact_PSAP

image-20250310222443988

总体而言，COMMOT 工具的使用较为简便。然而，如何将其有效应用于实际研究中，并围绕生物学问题展开深入解析，则需要借助更多相关文献的阅读与学习，例如

PMID: 39558136
PMID: 38503922
PMID: 39471809

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。原始发表：2025-03-22，如有侵权请联系 cloudcommunity@tencent 删除数据库存储database可视化数据

发布者：admin，转转请注明出处：http://www.yc00.com/web/1748181122a4743336.html

空转细胞通讯分析之COMMOT