企业官网建设流程全解析

单细胞测序做到降维聚类那一步，很多人就停住了——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 python=3.10 conda activate cellcomm pip install scanpy==1.9.6 squidpy==1.4.0 omnipath==0.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做一次"体检"，检查项如下表。这张表是我从多次返工中总结出来的，每一条都对应过真实的报错或异常结果：

基因名格式这一条特别隐蔽。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_sum=1e4) sc.pp.log1p(adata) adata.layers['lognorm'] = adata.X.copy()

我习惯把归一化后的矩阵单独存一个 layer，因为后面可能还要做别的分析会改动X，留个备份心里踏实。跑ligrec时通过layer='lognorm'指定，避免用到被改过的X。这里target_sum=1e4是 Scanpy 的常用默认值，也有用1e6的，差别不大，关键是全流程统一。如果你之前已经做过归一化，检查一下adata.X.max()，对数化后的值一般在 10 以内，如果看到几百上千，说明没对数化或者用了别的变换。

3. 跑通第一次 ligrec：参数背后的取舍

3.1 核心调用与最小可运行示例

数据收拾好之后，第一次跑建议用最小配置，先看结果长什么样，别一上来就堆参数。核心调用就一行：

import squidpy as sq sq.gr.ligrec( adata, n_perms=100, cluster_key='cell_type', use_raw=False, layer='lognorm', copy=False )

跑完之后结果存在adata.uns['ligrec']里。n_perms=100是置换次数，官方默认是 1000，我建议第一次先用 100 快速看结果，确认流程通了再加大。置换次数直接决定 p 值的精度下限——100 次置换最小 p 值约 0.01，1000 次能到 0.001。如果你要做多重检验校正后还显著，100 次往往不够，因为校正后阈值会被压得很低。我一般正式跑用 1000，赶时间用 500，低于 100 的结果基本只能当探索看。

use_raw=False配合layer='lognorm'是我推荐的组合，明确告诉函数用哪个矩阵，避免它去翻raw里可能没归一化的数据。copy=False表示结果写回adata.uns，如果你想单独拿一个结果对象做后续处理，改成copy=True接返回值。

3.2 n_perms 与 seed 的稳定性权衡

置换检验有个绕不开的问题：结果有随机性。同样的数据跑两次，p 值会有细微差别，偶尔某对互作在 0.05 边缘反复横跳。解决办法是固定随机种子：

sq.gr.ligrec( adata, n_perms=1000, cluster_key='cell_type', layer='lognorm', seed=42, copy=False )

seed参数在较新版本里支持，老版本可能没有，如果没有就通过全局np.random.seed(42)控制。固定种子后结果可复现，写方法部分时也说得清楚。但要注意，种子固定不代表结果"绝对正确"，它只是让随机过程确定化。真正影响可信度的是n_perms和样本量。我在一个只有 3 个样本的数据集上试过，n_perms 加到 5000，显著互作数量反而比 1000 时少——因为置换更充分后，一些边缘信号被正确判定为不显著。这说明加大置换次数不是"让结果更多"，而是"让判断更准"。

3.3 结果对象的结构拆解

跑完之后得知道怎么取数据。adata.uns['ligrec']是一个字典，关键字段有这几个：

• means：配体受体对在各细胞群组合上的平均表达，是个多层索引的 DataFrame
• pvalues：对应的置换检验 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 跑下来，检验的配体受体对数量轻松上千，如果不做多重检验校正，按 p<0.05 筛出来的"显著"互作里有一大半是假阳性。这是统计学的常识，但我在不少分析报告里看到有人直接用原始 p 值，结果列了几百对互作，看着热闹实则不可靠。校正方法我用 Benjamini-Hochberg（BH），控制 FDR：

from statsmodels.stats.multitest import multipletests import numpy as np flat_p = pvals.values.flatten() valid = ~np.isnan(flat_p) _, fdr, _, _ = multipletests(flat_p[valid], method='fdr_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，得先剔除再填回。校正后一般用 FDR<0.05 作为阈值，严格一点用 0.01。我实测下来，一个中等规模数据集（10 个细胞群、约 2000 对配体受体）校正后能留下的显著互作通常在几十到一两百对之间，如果筛出来还是几百上千，八成是校正没做对或者数据本身噪声太大。

4.2 表达比例过滤：均值骗局的破解之道

前面提过群均值会掩盖表达异质性，这里给出具体解法。对每一对显著互作，回头检查配体在源细胞群、受体在靶细胞群里的表达细胞比例，比例太低的直接砍掉：

def expr_fraction(adata, cluster, gene, layer='lognorm', thresh=0.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, cmap='viridis', ax=ax, cbar_kws={'label': 'mean expression'}) ax.set_title(f'Ligand-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, weight=cnt) pos = nx.spring_layout(G, seed=42, k=1.5) weights = [G[u][v]['weight'] for u, v in G.edges()] nx.draw(G, pos, with_labels=True, node_color='lightblue', node_size=2000, width=[w*0.5 for w in weights], arrows=True, arrowsize=20, font_size=10)

k参数控制节点间距，图太挤就调大。边宽用互作数缩放，但要注意缩放系数，w*0.5是我试出来的经验值，太大边会糊成一团。网络图适合展示"谁指向谁"的方向性，比如你会发现某个群主要是信号发出者（出边多），另一个群主要是接收者（入边多），这种不对称性本身就是有意思的生物学发现。

5.3 空间数据叠加：Squidpy 的独门优势

如果你的数据是空间转录组（比如 Visium、MERFISH），Squidpy 还能把互作结果叠回空间坐标，看信号传递在组织切片上的空间分布。这是它区别于纯单细胞工具的核心价值。基本流程是先算空间邻域，再结合配体受体表达：

sq.gr.spatial_neighbors(adata, coord_type='grid', n_neighs=6) sq.gr.nhood_enrichment(adata, cluster_key='cell_type')

coord_type='grid'适用于网格排列的 spot 数据，n_neighs=6是六邻域，Visium 常用。邻域富集分析能告诉你哪些细胞群在空间上倾向于相邻，再结合通讯结果，就能推断"空间邻近 + 配体受体共表达"的互作更可信。我在一个肿瘤切片数据上做过这个组合分析，发现某对配体受体在空间邻近的肿瘤-免疫细胞边界处显著富集，而单纯看单细胞数据时这对互作并不突出——空间信息补上了关键一环。不过空间数据的通讯分析对数据质量要求更高，spot 里往往是多个细胞混合，解读时要留有余地。

6. 那些让我返工的坑与修法

6.1 结果全空或显著互作异常少

第一次跑 ligrec 结果为空，是最打击人的情况。我排查过几次，原因基本集中在三处。第一是基因名不匹配，前面说过，OmniPath 用全大写 HGNC 符号，如果你的var_names是 Ensembl ID 或者小写，匹配率接近零。修法是统一大小写，或者用sc.read_10x_mtx时确保var_names='gene_symbols'。第二是矩阵没归一化，原始 count 下均值差异被测序深度主导，置换检验很难出显著。修法是确认layer='lognorm'且该 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_perms=1000 的数据集上遇到过内存峰值接近 16G 的情况。应对办法有几个：一是降低n_perms到 500，精度损失可接受；二是减少细胞群数量，把生物学意义相近的亚群合并；三是分批次跑，先按大群跑一轮，再对感兴趣的群细分跑。另外copy=False比copy=True省内存，因为不用复制整个 AnnData。

运行超时方面，ligrec 是单线程的，群多的时候确实慢。我的经验是 10 个群、1000 次置换大概几分钟到十几分钟，如果超过半小时还没动静，检查一下是不是配体受体对数量异常多（比如资源文件加载了全库），或者数据里基因数特别大。必要时可以先在子集上测试流程，确认没问题再上全量。

6.4 结果与生物学预期不符时怎么想

跑出来的结果和预期对不上，是最考验分析功底的时候。我的处理原则是：先怀疑技术原因，再考虑生物学解释。技术原因包括前面说的归一化、基因名、细胞群定义；生物学原因则可能是样本状态、处理条件、细胞比例差异。有一次我在一个对照组里期待看到某条炎症通路的互作，结果没出来，排查半天发现是那批对照样本的细胞活力偏低，某些配体表达被系统性压低。这种时候不能硬解释，得回到原始数据看 QC 指标。

还有一种情况是结果"太完美"——所有预期通路都显著，一个不落。这时候反而要警惕，可能是细胞群注释太粗或者数据有批次效应没处理干净，导致群间差异被夸大。我一般会做个阴性对照：随机打乱细胞群标签再跑一次，如果打乱后还能出一堆显著互作，说明当前结果的可靠性存疑。

7. 从结果到故事：把互作分析写进结论

跑完、筛完、画完图，最后一步是把结果组织成能讲通的逻辑。我的习惯是先定一个核心问题，比如"在这个组织里，哪类细胞是信号网络的主导者"，然后围绕它挑证据：互作数量热图证明某群是枢纽，dot plot 展示关键配体受体的表达，网络图呈现方向性，空间叠加（如果有）佐证邻近关系。四张图从不同角度指向同一个结论，比堆十张图更有说服力。

另外提醒一点：通讯分析的结果本质是假设生成，不是因果证明。写结论时用"提示""可能""关联"这类措辞，别写成"证明了某某细胞通过某某通路调控某某细胞"。真正要坐实，还得靠实验验证或者更严格的因果推断方法。我在项目里通常把通讯结果作为后续实验设计的依据，比如挑出 top 几对互作去做阻断实验或者共培养验证，这样整个故事才闭环。

最后分享一个实用习惯：每次跑完 ligrec，把关键参数（n_perms、seed、layer、资源库版本）、筛选阈值（FDR、表达比例）、以及最终显著互作数量记到一个日志文件里。单细胞项目周期长，隔几个月回头补分析时，没有这些记录根本复现不了当时的结果。这个习惯帮我省过好几次返工，尤其是在合作项目里需要向别人解释"这个数字是怎么来的"的时候。