单细胞测序细胞通讯分析：Scanpy与Squidpy实战指南

单细胞测序做到降维聚类那一步很多人就停住了——UMAP 上花花绿绿的一堆簇看着挺漂亮但接下来呢簇和簇之间到底在聊什么谁在给谁发信号哪条通路在特定条件下被激活这些才是真正能写进文章、能讲出生物学故事的东西。细胞通讯分析就是干这个的把表达矩阵翻译成谁分泌配体、谁表达受体、哪对细胞群之间存在显著互作的关系网络。Scanpy 负责把数据整理干净、聚类注释到位Squidpy 则在这个基础上做空间邻域和配体-受体互作的统计推断。这套组合我用了不下十几个项目从肿瘤微环境到发育轨迹都跑过踩过的坑能写满一页纸。下面这份教程不是官方文档的翻译而是我实际跑通、调过参、修过报错之后沉淀下来的完整流程适合已经会基本 Scanpy 操作、想往细胞通讯方向深入的分析人员也适合刚接手单细胞项目、需要快速产出互作结果的同学。1. 先想清楚细胞通讯分析到底在算什么1.1 从表达矩阵到对话关系的翻译逻辑很多人一上来就急着调sq.gr.ligrec结果跑完对着一堆 p 值发懵不知道哪些能信、哪些是噪声。问题的根源在于没搞明白这套方法在算什么。细胞通讯分析的核心假设很朴素如果一群细胞高表达某个配体另一群细胞高表达对应的受体那么这两群细胞之间就可能存在信号传递。注意我用了可能因为这只是统计关联不是实验验证。Squidpy 的ligrec模块做的事情就是把配体-受体对ligand-receptor pairs从一个预置数据库里拿出来逐对检查两个细胞群里的表达情况然后用置换检验permutation test判断这种共表达是否显著高于随机。这里有个关键点容易被忽略它算的是群层面的平均表达不是单细胞层面。也就是说哪怕某个细胞群里只有 5% 的细胞表达这个配体只要均值够高它就可能被判为显著。这就解释了为什么有些结果看起来不太合理——因为群均值掩盖了表达的异质性。我在一个肿瘤浸润淋巴细胞的数据集上就遇到过某个趋化因子配体在 T 细胞群里均值很高但实际表达它的细胞不到 10%剩下 90% 是阴性。这种结果如果直接写进结论审稿人一问就露馅。所以后面我会专门讲怎么用表达比例做二次过滤。1.2 Scanpy 和 Squidpy 的分工边界这两个库的关系得理清楚不然代码写起来会乱。Scanpy 是单细胞分析的主力负责 QC、归一化、高变基因、PCA、邻居图、聚类、UMAP、差异表达这一整套。Squidpy 是建立在 AnnData 之上的扩展专注两件事空间转录组的邻域富集分析以及配体-受体互作分析。它的ligrec功能其实不依赖空间信息普通单细胞数据也能跑只要你的adata.obs里有细胞类型注释就行。我习惯的分工是Scanpy 把adata收拾干净确保obs里有cell_type这一列X是归一化且对数化后的表达矩阵或者layers[lognorm]然后交给 Squidpy 做互作推断。千万别在原始 count 上直接跑 ligrec这是新手最常见的错误之一后面会详细说为什么。另外 Squidpy 的版本更新比较快API 有过几次变动我下面给的代码基于较新的稳定版本如果你用的是老版本函数名可能对不上建议先pip show squidpy确认一下。1.3 什么数据适合做、什么数据别硬做不是所有单细胞数据都适合做通讯分析。我的经验判断标准有三条第一细胞类型注释要足够可靠如果聚类结果本身就很糊通讯分析只会把噪声放大第二每个细胞群的细胞数不能太少低于 20 个细胞的群建议合并或剔除否则置换检验的统计效力不够p 值波动极大第三配体-受体数据库要和你的物种、组织匹配。Squidpy 默认用的是OmniPath资源覆盖人类和小鼠比较全但如果你做的是斑马鱼或者植物就得自己准备资源文件。还有一类数据要特别小心高度同质的细胞群。比如你把一群成纤维细胞细分成了 8 个亚群它们之间的通讯分析结果往往全是弱信号因为亚群之间表达谱差异本来就小。这种情况下与其硬做通讯不如先做差异表达看看亚群间到底有没有实质区别。我在一个纤维化模型里就吃过这个亏8 个成纤维亚群跑出来几十对显著互作结果一验证发现大部分是同一套配体受体在反复出现生物学意义很有限。2. 环境与数据准备别在第一步埋雷2.1 依赖安装与版本锁定的实操建议Squidpy 的依赖链比较长尤其是omnipath和statsmodels这两个版本不匹配时报错信息往往很隐晦。我的做法是单独建一个 conda 环境把关键包版本钉死。下面这套组合我在多个项目里验证过跑得比较稳conda create -n cellcomm python3.10 conda activate cellcomm pip install scanpy1.9.6 squidpy1.4.0 omnipath0.1.7 pip install leidenalg igraph为什么强调版本因为 Squidpy 1.3 到 1.4 之间ligrec的返回结构有过调整老教程里的res[means]在新版本里可能取不到得用res.mean。另外omnipath是 Squidpy 拉取配体受体资源的底层依赖它需要联网下载数据库缓存第一次跑的时候如果网络不稳会卡住甚至报错。我的处理办法是提前手动触发一次下载把缓存落到本地后面就快了import omnipath as op interactions op.interactions.import_intercell_network() print(interactions.shape)这一步跑通说明资源拉取没问题后面ligrec就不会在数据加载上翻车。2.2 AnnData 对象的状态检查清单在正式跑通讯分析前我会对adata做一次体检检查项如下表。这张表是我从多次返工中总结出来的每一条都对应过真实的报错或异常结果检查项期望状态不满足时的后果adata.obs[cell_type]存在且无缺失值ligrec 直接报 KeyErroradata.X或指定 layer归一化对数化后的值结果全是假阳性或全零细胞类型数量2 到 30 之间太多则组合爆炸太慢每群细胞数建议 ≥ 20置换检验 p 值不稳定基因名格式与资源库一致如全大写配体受体匹配不上结果为空是否含线粒体/核糖体基因建议已过滤干扰表达均值基因名格式这一条特别隐蔽。OmniPath 里的配体受体用的是 HGNC 标准符号全大写。如果你的var_names是小写或者带版本号比如从某些流程导出的ENSG00000...匹配率会骤降跑出来的互作对少得可怜。我一般会先做一次大小写统一adata.var_names adata.var_names.str.upper() adata.var_names_make_unique()var_names_make_unique()这行别省统一大小写后可能出现重名不加的话后面索引会出问题。2.3 归一化与对数化的正确姿势前面说过不能在原始 count 上跑这里展开讲为什么。ligrec内部计算的是表达均值如果输入是原始 count均值会被测序深度和细胞大小主导——测序深的细胞所有基因值都高导致它跟谁都能通讯。归一化把每个细胞的总表达拉到同一量级对数化则压缩动态范围让高表达基因不至于一家独大。标准流程是import scanpy as sc # 假设 adata 已经完成基本 QC sc.pp.normalize_total(adata, target_sum1e4) sc.pp.log1p(adata) adata.layers[lognorm] adata.X.copy()我习惯把归一化后的矩阵单独存一个 layer因为后面可能还要做别的分析会改动X留个备份心里踏实。跑ligrec时通过layerlognorm指定避免用到被改过的X。这里target_sum1e4是 Scanpy 的常用默认值也有用1e6的差别不大关键是全流程统一。如果你之前已经做过归一化检查一下adata.X.max()对数化后的值一般在 10 以内如果看到几百上千说明没对数化或者用了别的变换。3. 跑通第一次 ligrec参数背后的取舍3.1 核心调用与最小可运行示例数据收拾好之后第一次跑建议用最小配置先看结果长什么样别一上来就堆参数。核心调用就一行import squidpy as sq sq.gr.ligrec( adata, n_perms100, cluster_keycell_type, use_rawFalse, layerlognorm, copyFalse )跑完之后结果存在adata.uns[ligrec]里。n_perms100是置换次数官方默认是 1000我建议第一次先用 100 快速看结果确认流程通了再加大。置换次数直接决定 p 值的精度下限——100 次置换最小 p 值约 0.011000 次能到 0.001。如果你要做多重检验校正后还显著100 次往往不够因为校正后阈值会被压得很低。我一般正式跑用 1000赶时间用 500低于 100 的结果基本只能当探索看。use_rawFalse配合layerlognorm是我推荐的组合明确告诉函数用哪个矩阵避免它去翻raw里可能没归一化的数据。copyFalse表示结果写回adata.uns如果你想单独拿一个结果对象做后续处理改成copyTrue接返回值。3.2 n_perms 与 seed 的稳定性权衡置换检验有个绕不开的问题结果有随机性。同样的数据跑两次p 值会有细微差别偶尔某对互作在 0.05 边缘反复横跳。解决办法是固定随机种子sq.gr.ligrec( adata, n_perms1000, cluster_keycell_type, layerlognorm, seed42, copyFalse )seed参数在较新版本里支持老版本可能没有如果没有就通过全局np.random.seed(42)控制。固定种子后结果可复现写方法部分时也说得清楚。但要注意种子固定不代表结果绝对正确它只是让随机过程确定化。真正影响可信度的是n_perms和样本量。我在一个只有 3 个样本的数据集上试过n_perms 加到 5000显著互作数量反而比 1000 时少——因为置换更充分后一些边缘信号被正确判定为不显著。这说明加大置换次数不是让结果更多而是让判断更准。3.3 结果对象的结构拆解跑完之后得知道怎么取数据。adata.uns[ligrec]是一个字典关键字段有这几个means配体受体对在各细胞群组合上的平均表达是个多层索引的 DataFramepvalues对应的置换检验 p 值metadata这次运行的参数记录包括 n_perms、seed 等取显著互作的典型写法res adata.uns[ligrec] means res[means] pvals res[pvalues] # 展平多层索引方便筛选 means_df means.copy() means_df.columns [_.join(map(str, c)) for c in means_df.columns]这里有个坑means的列索引是(cluster1, cluster2, ligand, receptor)这样的多层结构直接筛选很别扭。我习惯先展平成字符串列名再配合 p 值做联合筛选。展平后列名形如Tcell_Bcell_CCL5_CCR1虽然长但一目了然后面做热图或者导出表格都方便。另外pvalues的结构和means对齐可以直接按位置对应不用额外对齐索引。4. 结果筛选把统计输出变成可信结论4.1 多重检验校正不能省一次 ligrec 跑下来检验的配体受体对数量轻松上千如果不做多重检验校正按 p0.05 筛出来的显著互作里有一大半是假阳性。这是统计学的常识但我在不少分析报告里看到有人直接用原始 p 值结果列了几百对互作看着热闹实则不可靠。校正方法我用 Benjamini-HochbergBH控制 FDRfrom statsmodels.stats.multitest import multipletests import numpy as np flat_p pvals.values.flatten() valid ~np.isnan(flat_p) _, fdr, _, _ multipletests(flat_p[valid], methodfdr_bh) fdr_full np.full(flat_p.shape, np.nan) fdr_full[valid] fdr fdr_df pvals.copy() fdr_df.values[:] fdr_full.reshape(pvals.shape)注意pvals里可能有 NaN某些配体受体对在特定群里完全没表达无法计算multipletests不接受 NaN得先剔除再填回。校正后一般用 FDR0.05 作为阈值严格一点用 0.01。我实测下来一个中等规模数据集10 个细胞群、约 2000 对配体受体校正后能留下的显著互作通常在几十到一两百对之间如果筛出来还是几百上千八成是校正没做对或者数据本身噪声太大。4.2 表达比例过滤均值骗局的破解之道前面提过群均值会掩盖表达异质性这里给出具体解法。对每一对显著互作回头检查配体在源细胞群、受体在靶细胞群里的表达细胞比例比例太低的直接砍掉def expr_fraction(adata, cluster, gene, layerlognorm, thresh0.1): mask adata.obs[cell_type] cluster sub adata[mask] vals sub[:, gene].layers[layer] if layer else sub[:, gene].X vals np.asarray(vals.todense()).flatten() if hasattr(vals, todense) else np.asarray(vals).flatten() return (vals thresh).mean()这个函数返回某基因在某群里表达值超过阈值的细胞占比。我的经验阈值是配体在源群的表达比例 ≥ 0.1受体在靶群 ≥ 0.1低于这个就标记为低置信并从主结论里剔除。为什么用 0.1 而不是更高因为单细胞数据 dropout 严重真实表达比例会被低估0.1 是个折中。如果你数据质量好比如有生物学重复、测序深度高可以提到 0.2。这一步做完显著互作数量通常会再砍掉三到五成但剩下的每一对都经得起推敲。4.3 按细胞群组合聚合看整体格局逐对看互作容易迷失在细节里我习惯先看宏观格局哪些细胞群组合之间的互作最密集。做法是把显著互作按(源群, 靶群)聚合计数sig_mask fdr_df 0.05 sig_pairs means_df[sig_mask.values] # 按源群-靶群统计显著互作数 counts {} for col in sig_pairs.columns: parts col.split(_) src, tgt parts[0], parts[1] counts[(src, tgt)] counts.get((src, tgt), 0) 1聚合完画个热图横轴靶群、纵轴源群颜色深浅代表互作数量。这张图往往能一眼看出通讯枢纽——比如某个细胞群作为源或靶的互作特别多它可能就是微环境里的关键调控者。我在一个炎症模型里就靠这张图发现中性粒细胞是互作网络的中心节点后续实验也验证了它在招募其他免疫细胞中的作用。这种宏观视角比逐对罗列更有叙事性写文章时也更好组织逻辑。5. 可视化让互作图能直接进文章5.1 互作数量热图与 dot plot 的组合拳Squidpy 自带了一些可视化函数但默认样式偏朴素我一般自己用 matplotlib 或 seaborn 重画。互作数量热图上面已经提过这里补一个 dot plot 的画法用来展示特定配体受体对在各群间的表达和显著性import matplotlib.pyplot as plt import seaborn as sns # 选一对感兴趣的互作 target CCL5_CCR1 sub means_df[[c for c in means_df.columns if target in c]] sub_p fdr_df[[c for c in fdr_df.columns if target in c]] fig, ax plt.subplots(figsize(10, 4)) sns.heatmap(sub, cmapviridis, axax, cbar_kws{label: mean expression}) ax.set_title(fLigand-receptor pair: {target}) plt.tight_layout()dot plot 的精髓在于同时编码两个维度点的大小表示表达比例颜色表示显著性-log10 FDR。这种图在文章里很常见因为信息密度高。我通常会对 top 20 的显著互作批量生成然后挑最有生物学意义的几张放进正文其余放补充材料。挑的原则是优先选已知通路比如趋化因子、生长因子优先选跨细胞类型同类型自分泌互作说服力弱优先选在多个样本里一致出现的。5.2 网络图把互作关系画成有向图热图看数量网络图看结构。用 networkx 把显著互作画成有向图节点是细胞群边是互作边宽表示互作数量import networkx as nx G nx.DiGraph() for (src, tgt), cnt in counts.items(): if cnt 3: # 只画互作数≥3的边避免图太乱 G.add_edge(src, tgt, weightcnt) pos nx.spring_layout(G, seed42, k1.5) weights [G[u][v][weight] for u, v in G.edges()] nx.draw(G, pos, with_labelsTrue, node_colorlightblue, node_size2000, width[w*0.5 for w in weights], arrowsTrue, arrowsize20, font_size10)k参数控制节点间距图太挤就调大。边宽用互作数缩放但要注意缩放系数w*0.5是我试出来的经验值太大边会糊成一团。网络图适合展示谁指向谁的方向性比如你会发现某个群主要是信号发出者出边多另一个群主要是接收者入边多这种不对称性本身就是有意思的生物学发现。5.3 空间数据叠加Squidpy 的独门优势如果你的数据是空间转录组比如 Visium、MERFISHSquidpy 还能把互作结果叠回空间坐标看信号传递在组织切片上的空间分布。这是它区别于纯单细胞工具的核心价值。基本流程是先算空间邻域再结合配体受体表达sq.gr.spatial_neighbors(adata, coord_typegrid, n_neighs6) sq.gr.nhood_enrichment(adata, cluster_keycell_type)coord_typegrid适用于网格排列的 spot 数据n_neighs6是六邻域Visium 常用。邻域富集分析能告诉你哪些细胞群在空间上倾向于相邻再结合通讯结果就能推断空间邻近 配体受体共表达的互作更可信。我在一个肿瘤切片数据上做过这个组合分析发现某对配体受体在空间邻近的肿瘤-免疫细胞边界处显著富集而单纯看单细胞数据时这对互作并不突出——空间信息补上了关键一环。不过空间数据的通讯分析对数据质量要求更高spot 里往往是多个细胞混合解读时要留有余地。6. 那些让我返工的坑与修法6.1 结果全空或显著互作异常少第一次跑 ligrec 结果为空是最打击人的情况。我排查过几次原因基本集中在三处。第一是基因名不匹配前面说过OmniPath 用全大写 HGNC 符号如果你的var_names是 Ensembl ID 或者小写匹配率接近零。修法是统一大小写或者用sc.read_10x_mtx时确保var_namesgene_symbols。第二是矩阵没归一化原始 count 下均值差异被测序深度主导置换检验很难出显著。修法是确认layerlognorm且该 layer 确实是对数化后的值。第三是细胞群太少或太小只有两三个群、每群十几个细胞时置换空间有限p 值很难小。修法是合并小群、确保每群 ≥ 20 个细胞。还有一种隐蔽情况cluster_key指定的列里有 NaN 或者空字符串。ligrec遇到缺失的群标签会静默跳过导致结果比预期少。跑之前用adata.obs[cell_type].isna().sum()检查一下有缺失就补上或剔除。6.2 p 值全为 1 或大量 NaN 的排查链路p 值全 1 通常意味着置换检验没跑起来或者数据有问题。我的排查顺序是先看n_perms是不是设成了 0 或负数低级错误但真有人犯再看输入矩阵是不是全零或常数比如误用了某个空 layer最后检查配体受体对在数据里到底有没有表达——如果某个配体在所有细胞里都是 0那对互作的检验自然无意义返回 NaN 是正常的。大量 NaN 则往往是基因名匹配问题匹配上的基因太少大部分配体受体对找不到对应基因只能返回 NaN。修法上我建议跑之前先做个匹配率检查from omnipath.interactions import import_intercell_network net import_intercell_network() ligands set(net[net[source] ligand][genesymbol].str.upper()) receptors set(net[net[target] receptor][genesymbol].str.upper()) data_genes set(adata.var_names.str.upper()) print(ligand match:, len(ligands data_genes) / len(ligands)) print(receptor match:, len(receptors data_genes) / len(receptors))匹配率低于 0.3 就要警惕了正常人类数据应该在 0.5 以上。低于这个值先解决基因名问题再往下走否则后面全是无用功。6.3 内存溢出与运行超时的应对细胞群多、n_perms大的时候ligrec 会吃不少内存尤其是置换过程要反复打乱标签重算。我在一个 15 个细胞群、n_perms1000 的数据集上遇到过内存峰值接近 16G 的情况。应对办法有几个一是降低n_perms到 500精度损失可接受二是减少细胞群数量把生物学意义相近的亚群合并三是分批次跑先按大群跑一轮再对感兴趣的群细分跑。另外copyFalse比copyTrue省内存因为不用复制整个 AnnData。运行超时方面ligrec 是单线程的群多的时候确实慢。我的经验是 10 个群、1000 次置换大概几分钟到十几分钟如果超过半小时还没动静检查一下是不是配体受体对数量异常多比如资源文件加载了全库或者数据里基因数特别大。必要时可以先在子集上测试流程确认没问题再上全量。6.4 结果与生物学预期不符时怎么想跑出来的结果和预期对不上是最考验分析功底的时候。我的处理原则是先怀疑技术原因再考虑生物学解释。技术原因包括前面说的归一化、基因名、细胞群定义生物学原因则可能是样本状态、处理条件、细胞比例差异。有一次我在一个对照组里期待看到某条炎症通路的互作结果没出来排查半天发现是那批对照样本的细胞活力偏低某些配体表达被系统性压低。这种时候不能硬解释得回到原始数据看 QC 指标。还有一种情况是结果太完美——所有预期通路都显著一个不落。这时候反而要警惕可能是细胞群注释太粗或者数据有批次效应没处理干净导致群间差异被夸大。我一般会做个阴性对照随机打乱细胞群标签再跑一次如果打乱后还能出一堆显著互作说明当前结果的可靠性存疑。7. 从结果到故事把互作分析写进结论跑完、筛完、画完图最后一步是把结果组织成能讲通的逻辑。我的习惯是先定一个核心问题比如在这个组织里哪类细胞是信号网络的主导者然后围绕它挑证据互作数量热图证明某群是枢纽dot plot 展示关键配体受体的表达网络图呈现方向性空间叠加如果有佐证邻近关系。四张图从不同角度指向同一个结论比堆十张图更有说服力。另外提醒一点通讯分析的结果本质是假设生成不是因果证明。写结论时用提示可能关联这类措辞别写成证明了某某细胞通过某某通路调控某某细胞。真正要坐实还得靠实验验证或者更严格的因果推断方法。我在项目里通常把通讯结果作为后续实验设计的依据比如挑出 top 几对互作去做阻断实验或者共培养验证这样整个故事才闭环。最后分享一个实用习惯每次跑完 ligrec把关键参数n_perms、seed、layer、资源库版本、筛选阈值FDR、表达比例、以及最终显著互作数量记到一个日志文件里。单细胞项目周期长隔几个月回头补分析时没有这些记录根本复现不了当时的结果。这个习惯帮我省过好几次返工尤其是在合作项目里需要向别人解释这个数字是怎么来的的时候。