单细胞测序做到降维聚类和细胞类型注释这一步很多人就以为“分析到头了”。实际上真正让数据开口讲故事的是下一步——细胞通讯分析。我接触过的单细胞项目里十个有八个最终都要回答同一个问题这群细胞之间到底在“聊”什么谁在发信号谁在收信号哪些通路在特定条件下被异常激活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_totalsc.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置换次数1000次数越多 p 值越稳1000 是精度和耗时的平衡点threshold表达阈值0.1低于此表达比例的基因对不参与检验过滤噪声use_raw是否用 rawFalse用归一化后的 X 即可seed随机种子42保证结果可复现n_jobs并行核数4按机器配置调别开满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_sum1e4) # sc.pp.log1p(adata) # 2. 运行配体-受体分析 sq.gr.ligrec( adata, n_perms1000, cluster_keycell_type, threshold0.1, use_rawFalse, seed42, n_jobs4, copyFalse ) # 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], methodfdr_bh)[1] sig_final sig[sig[pval_adj] 0.05]3.2 空间数据的邻接图构建与通讯分析空间数据的流程多了一步构建空间邻接图。以 Visium 为例坐标存在adata.obsm[spatial]里。# 构建空间邻接图 sq.gr.spatial_neighbors( adata, coord_typegrid, # Visium 是网格排列 n_neighs6, # 每个 spot 取 6 个邻居 spatial_keyspatial ) # 查看邻接矩阵 print(adata.obsp[spatial_connectivities].shape) # 空间通讯分析 sq.gr.ligrec( adata, n_perms1000, cluster_keycell_type, threshold0.1, seed42, n_jobs4 )coord_type这个参数要选对Visium 是规则网格用gridMERFISH、Slide-seq 这类不规则点用generic此时n_neighs表示取最近的 K 个邻居。选错了邻接图会失真后面所有空间分析都跟着错。空间通讯的结果解读和非空间略有不同非空间看的是“任意两类细胞之间”的互作空间数据可以进一步结合邻接图看“物理相邻的细胞对”之间是否富集了特定互作。Squidpy 提供了sq.gr.nhood_enrichment来做邻域富集分析配合通讯结果能讲出更完整的故事。3.3 结果可视化的三种出图方式Squidpy 自带了几种可视化函数我常用的有三种。第一种是互作数量热图看哪两类细胞之间互作最密集sq.pl.ligrec( adata, cluster_keycell_type, sourceCD8_T, targetMono, means_range(0.5, 3), alpha0.001, swap_axesTrue )第二种是配体-受体对的气泡图展示特定细胞对之间的显著互作sq.pl.ligrec( adata, cluster_keycell_type, sourceCD8_T, targetMono, pvalue_threshold0.01, means_range(0.5, 3), swap_axesFalse )第三种是自定义热图用 seaborn 画互作数量矩阵适合放进论文主图import seaborn as sns # 统计每对细胞的显著互作数 count_matrix sig_final.groupby([source, target]).size().unstack(fill_value0) plt.figure(figsize(10, 8)) sns.heatmap(count_matrix, cmapYlOrRd, annotTrue, fmtd) plt.title(Significant LR interactions) plt.tight_layout() plt.savefig(lr_heatmap.pdf, dpi300)提示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, ascendingFalse) # 按 LR 对聚合 lr_summary sig_final.groupby(lr_pair).agg( n_pairs(source, count) ).sort_values(n_pairs, ascendingFalse).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并行核数开太多反而增加内存峰值分批次跑先按细胞大类跑锁定感兴趣的方向后再细分。另外copyFalse让结果直接写回adata.uns比copyTrue返回新对象省内存。4.4 结果和生物学预期对不上这种情况我见得多了先别急着怀疑代码。第一通讯分析是相关性而非因果显著的 LR 对只说明共表达不代表信号真的传递了第二数据库覆盖有限你关心的通路可能压根不在内置数据库里第三组织特异性某些 LR 对在特定组织里表达低但功能重要被threshold过滤掉了。我的建议是把threshold调低跑一版宽松的结果和严格版对比看差异在哪里往往能发现被过滤掉的“暗信号”。4.5 常见问题速查表现象可能原因排查动作零互作基因名格式不匹配统一为基因符号检查大小写零互作表达矩阵未归一化确认 X 是 log 归一化数据p 值全小批次效应或异常高表达检查批次看 means 分布p 值全大细胞数太少合并小类群保证每类 50内存溢出n_perms 或 n_jobs 过大降到 1000 和 4 核结果不可复现未设 seed固定 seed42空间图错乱coord_type 选错Visium 用 grid其他用 generic5. 参数调优与结果可信度提升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固定下来每次调参后结果可复现排查问题时能省掉“到底是参数变了还是随机性作祟”的纠结。