1. 单细胞数据里的“邻里关系”为什么细胞通讯分析值得花时间做单细胞转录组的人大概都有过这样的体验花了两三天把聚类分群跑完UMAP 图上花花绿绿十几群T 细胞、B 细胞、髓系细胞、上皮细胞各自抱团看起来挺漂亮。但老板或者审稿人紧接着就会问一句——这些细胞群之间到底在“聊”什么谁在给谁发信号哪条通路在疾病组里明显变强了这就是细胞通讯分析要回答的问题。简单说它是在单细胞表达矩阵的基础上结合配体-受体相互作用的先验知识库推断不同细胞群之间可能存在的信号传递关系。Scanpy 负责把数据整理干净、降维聚类、把细胞群注释清楚Squidpy 则在这个基础上提供空间邻域分析和细胞通讯打分的能力。两者搭配基本能覆盖从“细胞是谁”到“细胞在哪”再到“细胞和谁说话”的完整链条。这套流程适合谁如果你已经跑过 Scanpy 的标准流程能独立完成质控、归一化、聚类和注释那这篇内容就是下一步的自然延伸。如果你还没接触过单细胞分析建议先把 AnnData 对象的结构、PCA 和 Leiden 聚类这些基础打牢否则后面调参数会非常痛苦。我见过不少同学一上来就想做通讯分析结果连adata.obsm里存的是什么都没搞清楚排查问题的时间远超跑流程的时间。需要提前说明的是细胞通讯分析本质上是一种基于先验知识的推断不是实验验证。它告诉你“这里可能存在信号”而不是“这里一定存在信号”。理解这一点后面看结果时才不会被一堆漂亮的圈图冲昏头脑。2. 整体思路与工具选型为什么是 Scanpy 加 Squidpy2.1 两个工具的分工边界Scanpy 是单细胞分析的主力框架数据容器 AnnData、预处理、降维、聚类、差异表达一条龙。Squidpy 是建立在 AnnData 生态之上的扩展库核心能力有两块一是空间邻域富集分析二是配体-受体相互作用分析squidpy.gr.ligrec。它的通讯分析模块直接吃 AnnData 对象和 Scanpy 的流程无缝衔接不需要在工具之间来回导表格。我对比过几种常见组合。有的方案用 R 侧的 CellPhoneDB统计严谨、数据库维护得不错但需要 Python 和 R 之间来回倒数据中间格式转换容易出错尤其是细胞注释标签和基因名大小写对不齐的时候排查起来很烦。也有用 NicheNet 的它侧重下游靶基因推断适合回答“信号进来之后引发了什么转录变化”但上手门槛更高。Squidpy 的优势在于全流程留在 Python 生态里AnnData 一路到底中间不落地文件对反复调参和可视化非常友好。当然Squidpy 的ligrec在统计检验上相对轻量置换检验的默认次数不算多复杂设计比如多条件对比需要自己多做一层封装。所以我的建议是探索性分析和快速出图用 Squidpy需要严格统计报告时再考虑用 CellPhoneDB 交叉验证。两者结论方向一致就放心往下走差异大的基因对单独拎出来看表达分布。2.2 数据库选择背后的逻辑Squidpy 的ligrec默认使用 OmniPath 资源可以通过database参数指定具体来源比如cellphonedb、consensus等。不同数据库的条目数量、置信度、组织特异性差别不小。OmniPath 的好处是整合了多个来源并做了注释坏处是条目多、噪声也相对多。我的习惯是先用consensus跑一遍看整体格局再用cellphonedb跑一遍做对照。如果两次结果里排名靠前的通路高度重合说明信号比较稳健如果差异很大就要警惕是不是数据库覆盖度导致的假象。另外ligrec支持传入自定义的相互作用表格式是两列配体、受体或者四列配体、受体、方向、来源这个口子在做特定通路聚焦时非常有用比如你只关心某个炎症轴就可以自己整理一张小表喂进去结果干净得多。2.3 分析粒度的取舍一个容易被忽略的问题是通讯分析应该在什么粒度上做用原始的 Leiden 聚类标签还是手工合并后的大类标签我的经验是两个都做但用途不同。细粒度聚类能暴露亚群之间的特异信号比如两个看似都是 T 细胞的亚群其中一个高表达某趋化因子受体通讯结果里就会体现出来。但细粒度也会带来大量低表达细胞群导致统计检验不稳定置换 p 值一片飘红。所以实操中我通常先在大类注释比如 T/NK、B/Plasma、Myeloid、Epithelial、Fibroblast上跑主分析保证每群细胞数足够经验值单群不低于 50 个细胞低于这个数结果波动会很大然后再对感兴趣的亚群单独拎出来做子集分析。这样既有全局视野又能下钻细节。3. 数据准备与预处理的关键细节3.1 从原始矩阵到干净的 AnnData假设你手上已经有一个经过基础质控的 AnnData 对象adata.X里是归一化后的表达值adata.obs里有聚类标签和细胞类型注释adata.var_names是基因名。这一步看似简单但有几个坑必须提前排掉。第一基因名格式。配体-受体数据库里的基因名通常是大写的人类符号如CCL5、CCR5如果你的数据来自小鼠需要先做同源基因映射或者确认数据库是否支持小鼠符号。Squidpy 本身不做物种转换这一步得自己处理。我一般用同源映射表把小鼠符号转成人类大写再喂给ligrec虽然会损失一部分基因但比让工具静默匹配失败强。第二表达值类型。ligrec内部会对配体和受体的表达做聚合默认取均值然后做置换检验。如果你的adata.X还是原始 counts均值会被测序深度带偏建议先做sc.pp.normalize_total加sc.pp.log1p确保是 log 归一化后的值。这一点在文档里写得不算显眼但影响很大。第三细胞类型标签的清理。adata.obs里如果有NaN、空字符串或者过于细碎的标签ligrec会报错或者产生大量无意义组合。跑之前先value_counts()看一眼把细胞数过少的群合并或者剔除。import scanpy as sc import squidpy as sq import pandas as pd # 检查细胞类型分布 print(adata.obs[cell_type].value_counts()) # 剔除细胞数过少的类型示例阈值 30 counts adata.obs[cell_type].value_counts() keep_types counts[counts 30].index adata adata[adata.obs[cell_type].isin(keep_types)].copy() # 确认表达值已归一化 print(adata.X[:5, :5])3.2 空间数据与非空间数据的区别对待如果你做的是空间转录组比如 Visium 或类似平台AnnData 里还会有obsm[spatial]存坐标obsp里可能有空间邻接图。这种情况下通讯分析可以结合空间邻域信息Squidpy 的sq.gr.spatial_neighbors能构建空间图后续ligrec虽然不直接用这个图但你可以用空间邻域富集的结果来佐证通讯结论——比如某条通路在空间上确实富集在两类细胞交界处说服力就强很多。如果是普通单细胞悬液数据没有空间坐标那就老老实实做基于表达谱的通讯推断别硬套空间分析。我见过有人把 PCA 前两个主成分当坐标喂给空间邻域函数结果出来的“邻域富集”毫无生物学意义纯属自欺欺人。3.3 高变基因与降维的复用通讯分析本身不需要重新做降维但如果你还没做过聚类需要先走一遍 Scanpy 标准流程。这里提醒一点sc.pp.highly_variable_genes选出的高变基因是给降维聚类用的ligrec用的是全基因集或者你指定的基因子集两者不要混淆。有些人为了省内存把adata裁到只剩高变基因再跑通讯分析结果配体受体大量缺失结论自然不可靠。# 标准预处理如果尚未做 sc.pp.normalize_total(adata, target_sum1e4) sc.pp.log1p(adata) sc.pp.highly_variable_genes(adata, n_top_genes2000) adata.raw adata # 保留全基因 adata adata[:, adata.var.highly_variable].copy() sc.pp.scale(adata, max_value10) sc.tl.pca(adata, n_comps30) sc.pp.neighbors(adata, n_pcs30) sc.tl.leiden(adata, resolution0.8)注意adata.raw的保留时机——在裁剪高变基因之前存一份后面做通讯分析或者画基因表达图时从raw里取全基因避免基因缺失。4. 通讯分析实操从参数设置到结果解读4.1 ligrec 的核心参数逐个拆解sq.gr.ligrec的签名不算复杂但每个参数都影响结果值得逐个说清楚。sq.gr.ligrec( adata, cluster_keycell_type, interactionsNone, # 自定义相互作用表None 则用数据库 databaseconsensus, # 数据库来源 n_perms1000, # 置换次数 threshold0.1, # 表达阈值 corr_methodfdr_bh, # 多重检验校正 use_rawTrue, # 从 raw 取表达 copyFalse, n_jobs4 )cluster_key指定用哪一列做细胞群划分必须是adata.obs里的类别型列。n_perms是置换检验次数默认 1000细胞群多、基因对多的时候计算量不小。我的经验是探索阶段用 100 到 200 快速看格局正式出图再提到 1000。threshold控制配体或受体在细胞群中的最低表达比例低于这个比例的基因对会被过滤默认 0.1 意味着至少 10% 的细胞表达才纳入。这个值调高比如 0.2结果更干净但可能漏掉低丰度信号调低则噪声增多建议根据数据质量在 0.1 到 0.25 之间试。corr_method做多重检验校正fdr_bh是常用的 Benjamini-Hochberg。注意置换 p 值和校正后的 p 值要分开看前者反映统计显著性后者控制假阳性。实际筛结果时我通常要求校正后 p 值小于 0.05同时均值表达差means里的差值也要有实际意义不能只看 p 值。4.2 结果对象的结构跑完之后adata.uns[ligrec]里会存一个结果字典包含means各细胞群对的配体受体平均表达、pvalues置换 p 值、metadata相互作用对信息等。means是一个多层索引的 DataFrame行是配体-受体对列是细胞群对比如T_cell|B_cell表示配体来自 T 细胞、受体在 B 细胞。这个方向性很重要——A|B和B|A是两回事前者是 A 发信号给 B后者反过来。res adata.uns[ligrec] means res[means] pvals res[pvalues] # 查看前几行 print(means.head()) print(pvals.head()) # 筛选显著且表达有差异的相互作用 sig_mask pvals 0.05 sig_means means[sig_mask] print(f显著相互作用对数量: {sig_mask.sum().sum()})这里有个细节means和pvals的形状一致但pvals里可能有NaN表示该组合无法计算比如某群细胞完全不表达配体。筛的时候要同时处理NaN否则布尔索引会出问题。4.3 可视化圈图、热图与点图的选择Squidpy 自带sq.pl.ligrec能画圈图chord diagram和热图。圈图适合展示整体格局一眼看出哪几个细胞群是“信号枢纽”热图适合看特定配体受体对在各群对之间的表达模式。我的习惯是先用圈图给个全局印象再用热图下钻到具体通路。sq.pl.ligrec( adata, cluster_keycell_type, source_groupsT_cell, target_groups[B_cell, Myeloid], means_range(0.5, 4), alpha0.05, swap_axesTrue )source_groups和target_groups可以限定只看特定方向避免图太乱。means_range控制颜色映射范围数据动态范围大的时候手动设一下否则颜色全挤在一端。alpha是显著性阈值超过的相互作用不画。圈图好看但信息密度有限真正写报告时我更依赖热图加表格。热图的行是配体受体对列是细胞群对颜色是平均表达星号标显著性。这种图在补充材料里很实用审稿人也能看懂。4.4 一个完整的分析片段把上面的步骤串起来一个可复现的最小流程大概是这样import scanpy as sc import squidpy as sq import matplotlib.pyplot as plt # 1. 读取已注释数据 adata sc.read_h5ad(annotated_data.h5ad) # 2. 确认表达值已归一化保留 raw if raw not in adata: adata.raw adata # 3. 清理细胞类型标签 counts adata.obs[cell_type].value_counts() adata adata[adata.obs[cell_type].isin(counts[counts 30].index)].copy() # 4. 跑通讯分析 sq.gr.ligrec( adata, cluster_keycell_type, databaseconsensus, n_perms1000, threshold0.1, corr_methodfdr_bh, use_rawTrue, n_jobs4 ) # 5. 提取结果 means adata.uns[ligrec][means] pvals adata.uns[ligrec][pvalues] # 6. 筛选显著结果并导出 sig (pvals 0.05) means.notna() sig_pairs means[sig].stack().reset_index() sig_pairs.columns [ligand_receptor, cell_pair, mean_expr] sig_pairs sig_pairs.sort_values(mean_expr, ascendingFalse) sig_pairs.to_csv(significant_interactions.csv, indexFalse) # 7. 可视化 sq.pl.ligrec( adata, cluster_keycell_type, means_range(0.5, 4), alpha0.05 ) plt.savefig(ligrec_chord.png, dpi300, bbox_inchestight)这段代码跑下来快的话十几分钟慢的话半小时到一小时取决于细胞数和置换次数。导出 CSV 那步很关键后面做通路富集或者和别的样本对比都靠它。5. 常见问题与排查实录5.1 结果全是 NaN 或者显著对少得可怜这是最高频的问题。原因通常有三类表达值没归一化、基因名不匹配、细胞群细胞数太少。排查顺序建议是先看adata.X的数值范围如果是几十上百的整数说明还是 counts先归一化再看adata.var_names前几个和数据库基因名对不对得上大小写、连字符、版本号后缀都可能是元凶最后看细胞群分布把小于 30 个细胞的群合并。还有一种隐蔽情况use_rawTrue但adata.raw不存在或者基因集不对。这时候ligrec可能静默用adata.X也可能报错取决于版本。稳妥做法是显式检查adata.raw是否存在且包含配体受体基因。5.2 置换检验跑得太慢细胞群多、n_perms高、基因对多三者相乘就是计算量。加速手段有几个降低n_perms到 100 先看格局用n_jobs开多核缩小interactions到自定义的感兴趣通路表或者先用threshold过滤掉低表达基因对。我实测下来把n_perms从 1000 降到 200时间能省一大半结论方向基本不变只是 p 值精度下降正式报告再补跑。5.3 不同数据库结果差异大前面提过consensus和cellphonedb覆盖度不同结果有差异是正常的。排查思路是对比两次结果里排名前 20 的配体受体对看重合比例。如果重合度高说明核心信号稳健如果几乎不重合检查是不是某个数据库对特定基因家族覆盖不全。另外ligrec的database参数在不同版本里支持的取值可能变化跑之前help(sq.gr.ligrec)确认一下当前版本支持哪些。5.4 圈图太乱看不清圈图的线条数量随细胞群对和显著对数增长群一多就成蜘蛛网。解决办法限定source_groups和target_groups分批画提高alpha阈值只留最显著的或者干脆放弃圈图用热图加表格。我个人在正式报告里更倾向热图信息密度高且不依赖颜色区分线条。5.5 结果和生物学预期不符通讯分析是推断不是事实。如果结果和预期差很远先别急着怀疑代码检查几个方面细胞注释是否准确注释错了通讯自然错样本分组是否混在一起跑疾病组和对照组混跑会稀释信号应该分组跑再对比数据库是否覆盖你关心的通路冷门通路可能压根不在库里。我遇到过注释把一群成纤维细胞错标成上皮细胞导致整条上皮-免疫通讯轴消失回头修正注释后结果就合理了。问题现象可能原因排查动作结果大量 NaN表达值未归一化检查adata.X数值范围补做 normalize 和 log1p显著对极少基因名不匹配对比var_names与数据库基因名格式计算超时置换次数过高降n_perms开n_jobs缩小基因对范围圈图混乱细胞群过多分批画或改用热图与预期不符注释错误或分组混淆复核注释按条件分组重跑6. 让结果更可信的几个进阶做法6.1 分组对比而不是一锅炖如果你的数据包含疾病组和对照组千万别把所有样本混在一起跑一次通讯分析就完事。正确做法是按条件分组各自跑ligrec然后对比两组里同一细胞群对的显著相互作用差异。Squidpy 本身不直接提供组间对比函数需要自己写循环或者用 AnnData 的视图分别跑再把结果拼起来做差异分析。这一步多花的时间绝对值得因为很多有意义的信号恰恰藏在组间差异里混跑会被平均掉。6.2 用空间信息做交叉验证有空间数据的话sq.gr.spatial_neighbors构建空间邻接图后sq.gr.nhood_enrichment能算细胞类型在空间上的邻域富集。如果通讯分析显示 A 群发信号给 B 群而空间邻域分析也显示 A 和 B 在空间上显著相邻两条证据互相支撑结论可信度大幅提升。这种交叉验证在写文章时是加分项审稿人很吃这一套。6.3 表达分布图辅助判断统计显著不等于生物学重要。某条相互作用 p 值很小但配体和受体的表达都只在极少数细胞里实际信号强度可能很弱。我习惯对排名靠前的几对相互作用单独画配体和受体的 UMAP 或者小提琴图看表达是否集中在特定细胞群、分布是否合理。这一步能过滤掉不少统计显著但生物学意义存疑的结果。# 查看特定配体受体的表达分布 sc.pl.umap(adata, color[CCL5, CCR5], use_rawTrue, cmapviridis)6.4 结果的可复现性管理通讯分析涉及随机置换不同运行之间结果会有细微差异。为了可复现跑之前设随机种子sq.gr.ligrec内部用 numpy 随机数可以通过np.random.seed控制并把参数完整记录在脚本里。我见过有人跑完出图过两周想复现却记不清当时用的threshold和n_perms只能重跑浪费大量时间。养成把参数写进配置字典或者脚本头部的习惯省事。7. 我个人在实际操作中的几点体会跑通讯分析这几年最大的感受是别把它当成黑盒。工具给的是候选列表真正的判断还得靠生物学背景和对数据的熟悉程度。同一份数据换个数据库、调个阈值结果排序就会变所以别迷信某一个排名第一的相互作用多看几条通路的整体趋势。另外细胞注释的质量直接决定通讯分析的上限。注释错了后面所有分析都是空中楼阁。我现在的习惯是通讯分析之前先把 marker 基因的表达图逐个过一遍确认每个群的标签站得住脚再往下走。这一步多花半小时能省掉后面几小时的困惑。最后分享一个小技巧如果你只关心某几条特定通路与其跑全库再筛不如自己整理一张小的相互作用表传给interactions参数计算量小、结果干净、解释起来也聚焦。全库扫描适合探索阶段找线索聚焦分析适合验证假设两种用法配合着来效率最高。
