Scanpy+Squidpy细胞通讯分析实战：从降维聚类到信号网络解读

单细胞测序做到降维聚类和细胞类型注释这一步，很多人就以为“分析到头了”。实际上，真正让数据开口讲故事的，是下一步——细胞通讯分析。我接触过的单细胞项目里，十个有八个最终都要回答同一个问题：这群细胞之间到底在“聊”什么？谁在发信号，谁在收信号，哪些通路在特定条件下被异常激活？Scanpy 负责把表达矩阵整理成干净可用的 AnnData 对象，Squidpy 则在这个基础上把空间邻接关系和配体-受体互作一起算出来，两者搭起来就是一套相当顺手的细胞通讯分析流水线。这篇内容面向的是已经跑过基础单细胞流程、手里有注释好的 AnnData、想进一步挖掘细胞间信号网络的人，也适合刚接触 Squidpy、想搞清楚它和 CellPhoneDB 这类工具差异的读者。我会把整个流程拆成可复现的步骤，把参数背后的逻辑讲透，再补上我实际踩过的坑。

1. 整体思路与工具选型拆解

1.1 为什么是 Scanpy 加 Squidpy 这个组合

单细胞通讯分析的工具生态其实挺杂的，CellPhoneDB、NicheNet、CellChat、iTALK 各有拥趸。我选 Scanpy + Squidpy 这套组合，核心原因有三个。第一是数据结构的统一，Scanpy 的 AnnData 是单细胞领域事实上的容器标准，Squidpy 直接吃 AnnData，不需要在工具之间反复导表格、转格式，省掉大量胶水代码。第二是空间信息的原生支持，Squidpy 从设计之初就把空间邻接图（spatial neighbors）作为一等公民，如果你手上有空间转录组数据，它能直接把“物理上挨着的细胞”纳入通讯计算，这是纯表达量工具做不到的。第三是计算轻量，Squidpy 的配体-受体分析基于表达量置换检验，不依赖复杂的先验数据库打分模型，跑起来快，结果也相对好解释。

需要说清楚的是，Squidpy 的ligrec模块本质上是基于表达量的共表达置换检验，它回答的是“这对配体-受体在两类细胞里是否显著共表达”，而不是“信号强度有多大”。这一点和 CellPhoneDB 的思路接近，但 Squidpy 的实现更贴近 Scanpy 生态，参数调节也更透明。如果你要的是通路级别的因果推断，那 NicheNet 更合适；如果你要的是快速、可解释、能出图的互作概览，Squidpy 够用。

1.2 分析流程的骨架长什么样

整个流程我习惯分成五段：数据准备与质控复查 → 细胞类型注释确认 → 配体-受体数据库加载 → 置换检验计算互作 → 结果可视化与解读。前两步是地基，很多人注释做得粗糙，后面通讯结果就会充满假阳性；第三步决定了你分析的“词典”范围；第四步是核心计算；第五步才是出图讲故事。

这里有个容易被忽略的点：通讯分析对细胞类型的分辨率非常敏感。如果你把 T 细胞和 NK 细胞合并成一个“淋巴细胞”大类，那它们之间本该存在的互作就被抹掉了；反过来，如果聚类过细，每类只有十几个细胞，置换检验的统计功效又不够。我的经验是，每类细胞至少保证 50 个以上，最好 100 以上，再去做通讯分析，结果才稳。

1.3 空间数据与非空间数据的分叉处理

这一步要提前想清楚。如果你的是普通 scRNA-seq（没有空间坐标），那 Squidpy 的ligrec只依赖细胞类型标签和表达矩阵，走的是“类别对类别”的置换检验。如果你的是空间转录组（比如 Visium 或 MERFISH），那可以额外构建空间邻接图，用sq.gr.spatial_neighbors把物理邻近关系纳入进来，做“邻域内通讯”的分析。两条路的代码前半段一样，后半段分叉，我在第 3 节会分别给出来。

2. 数据准备与核心细节解析

2.1 AnnData 对象的体检清单

在动手算通讯之前，我强烈建议先对 AnnData 做一次“体检”。Squidpy 的ligrec要求输入的对象里，细胞类型信息存在obs的某一列，表达矩阵在X或layers里，且基因名要和配体-受体数据库对得上。我通常跑这么一段检查：

import scanpy as sc import squidpy as sq adata = sc.read_h5ad("annotated_data.h5ad") # 检查细胞类型列 print(adata.obs['cell_type'].value_counts()) # 检查基因名格式（人类通常是全大写，小鼠首字母大写） print(adata.var_names[:10].tolist()) # 检查是否有原始计数（置换检验建议用归一化后的数据，但保留 counts 层） print(adata.layers.keys())

这里有个关键决策：用哪一层表达量做通讯分析。Squidpy 官方建议用归一化并 log 转换后的数据（也就是 Scanpy 流程里sc.pp.normalize_total+sc.pp.log1p之后的X），而不是原始 counts。原因是置换检验比较的是表达量的相对高低，log 转换能压缩极端值的影响，让检验更稳健。但如果你后面还想做别的分析，记得把原始 counts 存在layers['counts']里别丢。

注意：基因名格式必须和数据库一致。Squidpy 内置的ligrec数据库用的是人类基因符号（全大写），如果你分析的是小鼠数据，需要先把基因名转成人类同源格式，或者换用小鼠版本的数据库。我见过不少人直接拿小鼠数据跑，结果一个互作都出不来，排查半天才发现是大小写问题。

2.2 细胞类型注释的复核要点

通讯分析的结果质量，七成取决于注释质量。我在跑ligrec之前会做两件事。第一是看每类细胞的标记基因是否干净，如果某类细胞的 top marker 里混着别的谱系基因，说明注释可能有问题。第二是检查类别数量，前面说过每类至少 50 个细胞，低于这个数我会考虑合并或剔除。

还有一个细节：亚群命名要语义清晰。Squidpy 的结果表里，source和target两列直接用的是obs里的类别名。如果你命名成 “0”、“1”、“2”，后面解读结果时得反复对照，非常痛苦。我习惯在注释阶段就把名字定成 “CD8_T”、“NK”、“Mono” 这种一眼能懂的形式。

2.3 配体-受体数据库的选择与加载

Squidpy 的ligrec默认用的是内置的omnipath数据库，通过sq.gr.ligrec的interactions_params参数可以指定。默认情况下它会加载一套预编译的配体-受体对，覆盖常见的免疫、生长因子、趋化因子通路。如果你想用自定义数据库，可以传一个DataFrame，要求两列：source（配体基因）和target（受体基因）。

# 查看默认数据库规模 from squidpy.gr import ligrec # 默认参数下会自动加载内置数据库

我实际用下来，内置数据库对免疫和肿瘤微环境覆盖得不错，但如果你研究的是神经、发育这类特殊组织，可能需要补充自定义的配体-受体对。补充时注意基因名必须和adata.var_names完全一致，差一个字符就匹配不上。

2.4 置换检验的参数逻辑

sq.gr.ligrec的核心参数有这么几个，我逐个说清楚它们的作用和取值逻辑：

n_perms这个参数值得多说两句。置换检验的原理是：把细胞类型标签随机打乱，重新计算配体-受体的共表达统计量，重复n_perms次，得到一个零分布，再看真实统计量在零分布里的位置，算出 p 值。1000 次是常规选择，如果你要做严格的多重检验校正后仍然显著，可以加到 5000，但耗时会线性增长。我一般先用 1000 跑一遍看整体格局，锁定感兴趣的细胞对之后再针对性加次数。

threshold参数控制的是表达比例门槛。比如设为 0.1，意味着某个配体在源细胞类群里至少要有 10% 的细胞表达，才纳入检验。这个参数设太低会引入大量噪声互作，设太高又会漏掉低表达但生物学重要的信号。0.1 是我试下来比较平衡的值，如果你的数据稀疏度高，可以降到 0.05。

3. 实操过程与核心环节实现

3.1 非空间数据的完整通讯分析流程

先给非空间数据的完整代码，这是最常用的场景。假设adata已经注释好，obs['cell_type']存着细胞类型。

import scanpy as sc import squidpy as sq import pandas as pd import matplotlib.pyplot as plt # 1. 确认数据状态 adata = sc.read_h5ad("annotated_data.h5ad") # 确保 X 是 log 归一化后的表达 # 如果之前没做，补上： # sc.pp.normalize_total(adata, target_sum=1e4) # sc.pp.log1p(adata) # 2. 运行配体-受体分析 sq.gr.ligrec( adata, n_perms=1000, cluster_key="cell_type", threshold=0.1, use_raw=False, seed=42, n_jobs=4, copy=False ) # 3. 结果存在 adata.uns['cell_type_ligrec'] 里 res = adata.uns['cell_type_ligrec'] print(res.keys()) # 通常包含 means, pvalues, metadata 等

跑完之后，结果主要在adata.uns['cell_type_ligrec']里，means是各细胞对之间配体-受体的平均表达统计量，pvalues是对应的显著性。接下来是过滤和可视化。

# 提取显著互作 means = res['means'] pvals = res['pvalues'] # 展平 MultiIndex，方便筛选 means_long = means.stack().reset_index() means_long.columns = ['source', 'target', 'lr_pair', 'mean'] pvals_long = pvals.stack().reset_index() pvals_long.columns = ['source', 'target', 'lr_pair', 'pval'] merged = means_long.merge(pvals_long, on=['source', 'target', 'lr_pair']) sig = merged[merged['pval'] < 0.05] print(f"显著互作数量: {len(sig)}")

这里有个实操细节：means和pvalues的索引结构是 MultiIndex，直接stack展平后列名会乱，我习惯手动重命名。另外，多重检验校正这一步别省，Squidpy 本身不自动做 FDR 校正，需要你自己用statsmodels或scipy补上。

from statsmodels.stats.multitest import multipletests sig = sig.copy() sig['pval_adj'] = multipletests(sig['pval'], method='fdr_bh')[1] sig_final = sig[sig['pval_adj'] < 0.05]

3.2 空间数据的邻接图构建与通讯分析

空间数据的流程多了一步：构建空间邻接图。以 Visium 为例，坐标存在adata.obsm['spatial']里。

# 构建空间邻接图 sq.gr.spatial_neighbors( adata, coord_type="grid", # Visium 是网格排列 n_neighs=6, # 每个 spot 取 6 个邻居 spatial_key="spatial" ) # 查看邻接矩阵 print(adata.obsp['spatial_connectivities'].shape) # 空间通讯分析 sq.gr.ligrec( adata, n_perms=1000, cluster_key="cell_type", threshold=0.1, seed=42, n_jobs=4 )

coord_type这个参数要选对：Visium 是规则网格，用"grid"；MERFISH、Slide-seq 这类不规则点，用"generic"，此时n_neighs表示取最近的 K 个邻居。选错了邻接图会失真，后面所有空间分析都跟着错。

空间通讯的结果解读和非空间略有不同：非空间看的是“任意两类细胞之间”的互作，空间数据可以进一步结合邻接图，看“物理相邻的细胞对”之间是否富集了特定互作。Squidpy 提供了sq.gr.nhood_enrichment来做邻域富集分析，配合通讯结果能讲出更完整的故事。

3.3 结果可视化的三种出图方式

Squidpy 自带了几种可视化函数，我常用的有三种。第一种是互作数量热图，看哪两类细胞之间互作最密集：

sq.pl.ligrec( adata, cluster_key="cell_type", source="CD8_T", target="Mono", means_range=(0.5, 3), alpha=0.001, swap_axes=True )

第二种是配体-受体对的气泡图，展示特定细胞对之间的显著互作：

sq.pl.ligrec( adata, cluster_key="cell_type", source="CD8_T", target="Mono", pvalue_threshold=0.01, means_range=(0.5, 3), swap_axes=False )

第三种是自定义热图，用 seaborn 画互作数量矩阵，适合放进论文主图：

import seaborn as sns # 统计每对细胞的显著互作数 count_matrix = sig_final.groupby(['source', 'target']).size().unstack(fill_value=0) plt.figure(figsize=(10, 8)) sns.heatmap(count_matrix, cmap='YlOrRd', annot=True, fmt='d') plt.title('Significant LR interactions') plt.tight_layout() plt.savefig('lr_heatmap.pdf', dpi=300)

提示：sq.pl.ligrec的means_range参数控制颜色映射的范围，默认值经常导致图偏暗或偏亮，建议根据实际means的分位数手动设，比如取 5% 和 95% 分位数。

3.4 结果表的整理与下游分析衔接

跑完通讯分析，结果表通常有几千行，直接看会晕。我习惯做两件事：一是按细胞对聚合，统计每对细胞的显著互作数，快速定位“通讯热点”；二是按配体-受体对聚合，看哪些通路反复出现。

# 按细胞对聚合 pair_summary = sig_final.groupby(['source', 'target']).agg( n_interactions=('lr_pair', 'count'), n_unique_lr=('lr_pair', 'nunique') ).sort_values('n_interactions', ascending=False) # 按 LR 对聚合 lr_summary = sig_final.groupby('lr_pair').agg( n_pairs=('source', 'count') ).sort_values('n_pairs', ascending=False).head(20)

这两张表是后续讲故事的核心素材。pair_summary告诉你哪两类细胞“聊得最热”，lr_summary告诉你哪些信号通路是“公共语言”。如果某个 LR 对在多个细胞对里都显著，那它很可能是这个组织里的核心通讯轴，值得重点展开。

4. 常见问题与排查技巧实录

4.1 一个互作都跑不出来怎么办

这是最高频的问题。我按排查顺序列一下：第一查基因名格式，人类全大写、小鼠首字母大写，和数据库对不上就全军覆没；第二查表达矩阵，如果X里全是原始 counts 没做 log 转换，表达值跨度太大，置换检验的统计量会失真；第三查细胞类型列，cluster_key指定的列名拼错、或者该列有大量 NaN，都会导致计算异常；第四查阈值，threshold设成 0.5 这种高值，大部分基因对都被过滤掉了。我遇到过最隐蔽的一次是adata.var_names里混着 Ensembl ID 和基因符号两种格式，导致匹配率极低，统一成基因符号后恢复正常。

4.2 p 值全都很小或全都很大

p 值分布异常通常指向两个方向。全都极小（比如一堆 0.001），可能是n_perms设得太低，零分布不够稳定，或者数据里存在批次效应导致某些基因在特定细胞类群里异常高表达，置换检验把它当成了强信号。全都很大，除了前面说的基因名问题，还要检查是不是细胞类型分得太细，每类细胞数太少，置换检验没有统计功效。我的处理办法是先看means的分布，如果means本身就很平，说明信号弱，不是 p 值计算的问题。

4.3 内存溢出与运行超时

ligrec在细胞类型多、n_perms大的时候很吃内存。我跑过一个 15 类细胞、5000 次置换的数据集，单机 32G 内存直接爆掉。解决办法有三个：降n_perms到 1000 先跑通；减n_jobs，并行核数开太多反而增加内存峰值；分批次跑，先按细胞大类跑，锁定感兴趣的方向后再细分。另外，copy=False让结果直接写回adata.uns，比copy=True返回新对象省内存。

4.4 结果和生物学预期对不上

这种情况我见得多了，先别急着怀疑代码。第一，通讯分析是相关性而非因果，显著的 LR 对只说明共表达，不代表信号真的传递了；第二，数据库覆盖有限，你关心的通路可能压根不在内置数据库里；第三，组织特异性，某些 LR 对在特定组织里表达低但功能重要，被threshold过滤掉了。我的建议是，把threshold调低跑一版宽松的结果，和严格版对比，看差异在哪里，往往能发现被过滤掉的“暗信号”。

4.5 常见问题速查表

5. 参数调优与结果可信度提升

5.1 置换次数与显著性的权衡

n_perms不是越大越好，它和你的检验目的挂钩。如果你只是做探索性分析、看整体格局，1000 次足够；如果结果要写进论文、做严格的 FDR 校正，建议 5000 次起步。我做过一个对比：同一份数据，1000 次置换下显著互作 320 对，5000 次下降到 280 对，差异主要集中在小效应量的互作上。这说明低置换次数会高估显著性，探索阶段可以接受，定稿阶段要加次数复核。

5.2 表达阈值的敏感性分析

threshold是另一个需要做敏感性分析的点。我习惯跑 0.05、0.1、0.2 三档，看显著互作数量的变化曲线。如果从 0.1 到 0.2 互作数断崖式下跌，说明你的数据里很多信号依赖低表达比例的细胞，这时候要谨慎，可能这些信号本身就不稳。如果三档之间变化平缓，说明结果对阈值不敏感，可信度高。

5.3 与其它工具交叉验证

单靠 Squidpy 的结果下结论有风险，我通常会拿 CellPhoneDB 或 CellChat 跑一版做交叉验证。重点看核心 LR 对是否在两个工具里都显著，如果只有 Squidpy 报显著，那这个互作要打个问号。交叉验证不需要全量对比，挑 top 20 的 LR 对核对即可，能省大量时间。

5.4 结果解读的边界意识

最后说个容易被忽视的点：通讯分析的结果是假设，不是结论。Squidpy 告诉你“这两类细胞可能通过某 LR 对交流”，但信号是否真的传递、下游通路是否激活，需要实验验证（比如空间共定位、通路富集、功能实验）。我在写报告时会把通讯结果定位成“提出假设的线索”，而不是“证明结论的证据”，这样既严谨，也给后续实验留了空间。

我个人在实际操作中的体会是，细胞通讯分析最耗时的部分从来不是跑代码，而是前期注释的打磨和后期结果的筛选解读。代码跑通只要半小时，但把几千行结果收敛成三五个有生物学意义的通讯轴，往往要花好几天。所以别急着出图，先把pair_summary和lr_summary两张表吃透，再决定讲哪个故事。另外一个小技巧：把seed固定下来，每次调参后结果可复现，排查问题时能省掉“到底是参数变了还是随机性作祟”的纠结。