# Python单细胞分析新选择：用pertpy轻松搞定milo差异丰度分析（附完整代码） 如果你是一位长期在Python生态里摸爬滚打的单细胞研究者，面对那些功能强大但深植于R语言土壤的工具，比如大名鼎鼎的`miloR`，是不是总有种“隔岸观火”的无力感？数据格式转换、环境配置、语法差异，每一步都可能成为分析路上的绊脚石。尤其是在进行差异丰度分析时，你明明知道`milo`方法能绕过传统聚类分析的局限，在KNN图上精准捕捉那些细微的细胞群体变化，却因为工具链的割裂而不得不放弃。 好消息是，这种“水土不服”的日子到头了。随着`pertpy`库的成熟与整合，原生的`milopy`功能现在已无缝融入Python的`scverse`生态。这意味着，你无需离开熟悉的`scanpy`环境，无需与`h5ad`文件格式斗智斗勇，就能在Python中完整复现甚至优化整个`milo`分析流程。今天，我们就来彻底拆解这个过程，从零开始，手把手带你用`pertpy`在Python里玩转差异丰度分析，把前沿方法真正变成你工具箱里趁手的兵器。 ## 1. 环境搭建与数据准备：打造专属分析工作流 工欲善其事，必先利其器。一个稳定、可复现的分析环境是高效科研的基石。对于单细胞分析，我强烈建议使用`conda`进行环境管理，它能有效解决包依赖冲突这个令人头疼的问题。 首先，我们创建一个独立的`conda`环境。这里选择Python 3.10，它在兼容性和稳定性上取得了很好的平衡。 ```bash conda create -n pertpy_sc python=3.10 -y conda activate pertpy_sc ``` 接下来安装核心的分析包。`pertpy`是今天的主角，它集成了包括`milo`在内的多种扰动分析工具。`scanpy`则是Python单细胞分析的事实标准，必须安装。为了绘图美观，我们也会安装`matplotlib`和`seaborn`。使用清华镜像源可以大幅加速下载过程。 ```bash pip install pertpy scanpy matplotlib seaborn -i https://pypi.tuna.tsinghua.edu.cn/simple ``` 安装完成后，在Jupyter Notebook或Python脚本中导入这些库，并设置好绘图风格。 ```python import pertpy as pt import scanpy as sc import numpy as np import pandas as pd import matplotlib.pyplot as plt import seaborn as sns sc.settings.set_figure_params(dpi=100, facecolor='white', frameon=False) sns.set_style("whitegrid") ``` 数据是分析的起点。`milo`分析要求输入数据是标准的`AnnData`对象，通常来自`scanpy`预处理后的输出。假设你已经完成了标准的质控、归一化、高变基因筛选和降维聚类，得到了一个名为`adata`的`AnnData`对象。关键点在于，你的数据中必须包含**样本标识信息**和**分组信息**，这是后续统计检验的基础。这些信息通常存储在`adata.obs`中。 例如，你的`adata.obs`里应该有类似这样的列： - `sample_id`: 每个细胞所属的生物学重复样本（如`patient_1`, `patient_2`, `sample_A_rep1`）。 - `condition`: 实验分组（如`control`, `treated`, `disease`, `healthy`）。 > **注意**：`sample_id`必须是能够唯一标识一个独立生物学样本的字段。如果你的技术重复（比如同一个样本分多个通道测序）需要合并，务必在预处理阶段完成。`milo`的统计模型是在样本层面进行的，混淆样本会导致结果无效。 ## 2. 理解核心：KNN图与邻域构建的艺术 在深入代码之前，我们有必要花点时间理解`milo`方法的精髓。传统的差异丰度分析往往依赖于事先定义的细胞聚类（clusters）。但聚类存在几个固有缺陷：边界是硬性的，忽略了细胞的连续性变化；聚类结果受分辨率和算法影响大；更重要的是，生物学上有意义的细微变化可能被淹没在宽泛的聚类中。 `milo`另辟蹊径，它不依赖离散的聚类，而是基于细胞在降维空间（如PCA、UMAP）中的**K近邻图**来定义“邻域”。简单来说，它为图上随机选取的一部分“索引细胞”建立一个邻域，这个邻域包含了该细胞及其K个最近的邻居。这些邻域是**部分重叠**的，就像 Venn 图一样，能更细腻地捕捉细胞状态的连续变化。 那么，`pertpy`中如何实现这一步呢？首先，我们需要初始化`Milo`对象，并将我们的`AnnData`数据加载进去。`pertpy`使用`MuData`对象来同时管理原始的基因表达数据（`rna` modality）和后续为`milo`分析生成的计数数据（`milo` modality）。 ```python # 初始化Milo分析对象 milo = pt.tl.Milo() # 将AnnData数据加载为MuData对象 mdata = milo.load(adata) print(mdata) ``` 你会看到输出显示这是一个包含两个`modality`的`MuData`对象：`rna`和`milo`。原始的基因表达和注释信息都完好地保存在`rna`中。 接下来是构建KNN图。这是整个流程中**至关重要且需要根据数据特性进行调优的一步**。我们使用`scanpy`的标准函数`pp.neighbors`来完成。 ```python # 构建KNN图。关键参数：n_neighbors, n_pcs sc.pp.neighbors(mdata["rna"], n_neighbors=50, n_pcs=30) ``` 这里有两个核心参数需要你仔细考量： - **`n_neighbors`**： 定义每个细胞在图中连接多少个最近邻。增大此值会使邻域更“大”、更连通，可能提高统计功效，但也可能模糊细微差异。通常设置在15到100之间，对于细胞数较多（>10万）的数据集，可以适当增大。 - **`n_pcs`**： 用于计算细胞间距离的主成分数量。应使用在`sc.pp.pca`中确定的、能够代表大部分生物异质性的PC数量。通常可以观察PCA碎石图（`sc.pl.pca_variance_ratio`）来定。 构建好KNN图后，就可以在图上“撒点”定义邻域了，对应`pertpy`的`make_nhoods`函数。 ```python # 在KNN图上定义代表性邻域 milo.make_nhoods(mdata["rna"], prop=0.2, seed=1234) ``` `prop`参数控制从所有细胞中抽样多大比例作为“索引细胞”。`prop=0.2`意味着随机抽取20%的细胞作为索引，每个索引细胞及其邻居构成一个邻域。对于中小型数据集（<3万细胞），0.1-0.2是合理的起点；对于超大型数据集，可以降低到0.05以提升计算速度。`seed`参数是为了保证结果的可重复性。 定义完邻域，务必检查一下邻域的大小分布，这直接影响后续统计检验的效力。 ```python # 检查邻域大小分布 nhood_size = np.array(mdata["rna"].obsm["nhoods"].sum(0)).ravel() plt.figure(figsize=(8,5)) plt.hist(nhood_size, bins=100, edgecolor='k', alpha=0.7) plt.xlabel('Number of cells in neighbourhood') plt.ylabel('Number of neighbourhoods') plt.axvline(x=50, color='r', linestyle='--', label='Min recommended size') plt.legend() plt.show() ``` 一个经验法则是：**平均邻域大小最好大于 `5 * N_samples`**。如果你的样本数是6，那么平均邻域大小最好在30以上。如果大部分邻域太小，你需要回头调整`make_nhoods`中的`prop`参数，或者增大`sc.pp.neighbors`中的`n_neighbors`值。 ## 3. 统计检验与结果解读：从计数到生物学发现 有了定义好的邻域，下一步就是统计每个样本在每个邻域里有多少个细胞。这步操作将我们的单细胞数据“聚合”到了邻域-样本的二维计数矩阵，为后续的广义线性模型检验做好准备。 ```python # 统计每个样本在每个邻域中的细胞数 mdata = milo.count_nhoods(mdata, sample_col="sample_id") ``` 这里的`sample_col`参数必须指定`adata.obs`中那个唯一标识生物学样本的列名。执行后，`mdata["milo"]`中会生成一个稀疏矩阵，行是邻域，列是样本，值就是细胞计数。 现在来到核心的差异丰度检验。`pertpy`在幕后使用类似于`edgeR`或`DESeq2`的负二项广义线性模型，来检验每个邻域中细胞丰度在不同实验条件间是否存在显著差异。我们需要准备好实验设计矩阵。 ```python # 准备实验设计。假设分组信息存储在adata.obs['condition']中 design_df = mdata["rna"].obs[['sample_id', 'condition']].drop_duplicates() design_df = design_df.set_index('sample_id') # 确保设计矩阵的行顺序与计数矩阵的列顺序一致 design_df = design_df.loc[mdata["milo"].uns['sample_col_order'], :] # 执行差异丰度检验 da_results = milo.da_nhoods(mdata, design='~ condition', design_df=design_df) ``` `design`参数指定了统计模型。最简单的就是`~ condition`，只考虑分组效应。如果你的实验设计更复杂，比如包含批次效应，模型可以写成`~ batch + condition`。`pertpy`会自动处理模型拟合、假设检验，并返回一个包含每个邻域统计结果的`DataFrame`。 让我们看看结果长什么样： ```python print(da_results.head()) ``` 输出会包含以下几列关键信息： - `Nhood`: 邻域ID - `logFC`: 对数倍数变化（例如，`treated`组相比`control`组）。正值表示在该邻域中，处理组的细胞更富集。 - `PValue`: 原始P值 - `SpatialFDR`: 经过空间错误发现率校正后的Q值。这是考虑了邻域间重叠性的多重检验校正，**通常我们以`SpatialFDR < 0.1`或`< 0.05`作为显著性阈值**。 - 可能还有其他统计量，如`logCPM`（平均对数计数）。 如何快速评估结果？我们可以画一个火山图。 ```python # 快速可视化：火山图 da_results['-log10(SpatialFDR)'] = -np.log10(da_results['SpatialFDR']) significant = da_results['SpatialFDR'] < 0.1 plt.figure(figsize=(7,5)) plt.scatter(da_results.loc[~significant, 'logFC'], da_results.loc[~significant, '-log10(SpatialFDR)'], c='grey', alpha=0.5, s=10, label='Not Sig') plt.scatter(da_results.loc[significant, 'logFC'], da_results.loc[significant, '-log10(SpatialFDR)'], c=da_results.loc[significant, 'logFC'], cmap='coolwarm', s=20, label='Sig (FDR<0.1)') plt.axhline(y=-np.log10(0.1), color='k', linestyle='--', linewidth=0.8, label='FDR=0.1') plt.axvline(x=0, color='k', linewidth=0.5) plt.xlabel('Log Fold Change') plt.ylabel('-log10(Spatial FDR)') plt.colorbar(label='LogFC') plt.legend() plt.tight_layout() plt.show() ``` 这张图能让你一眼看出哪些邻域发生了显著变化，以及变化的方向和幅度。 ## 4. 高级可视化与生物学注释：让结果“说话” 得到一堆统计数字还不够，我们需要将结果映射回可解释的生物学背景。`pertpy`提供了强大的可视化工具，能将抽象的统计结果与细胞的真实分布联系起来。 首先，我们可以将差异丰度结果投射回原始的UMAP图上，直观地看看到底是哪些区域的细胞丰度发生了变化。 ```python # 构建邻域图用于可视化 milo.build_nhood_graph(mdata) # 绘制差异丰度邻域图 pt.pl.plot_nhood_graph(mdata=mdata, alpha=0.1, # 显著性阈值 min_size=5, # 最小节点大小 min_logFC=0, # 显示所有显著邻域，无论logFC正负 plot_edges=True, title='Differential Abundance: Treatment vs Control', color_map='viridis') # 颜色映射，表示logFC大小 ``` 这张图里，每个节点代表一个邻域，节点大小与邻域中的细胞总数成正比，节点颜色代表logFC（例如，红色表示在处理组富集，蓝色表示在对照组富集）。节点之间的连线表示邻域之间有细胞重叠。通过这张图，你能快速定位到发生显著变化的细胞群落所在的生物学区域。 但是，我们更关心的是：这些差异丰度的邻域主要是什么细胞类型？这就需要用到细胞类型注释信息。假设你的`adata.obs`中有一个`celltype`列。 ```python # 用细胞类型注释邻域 milo.annotate_nhoods(mdata, anno_col="celltype") # 绘制蜂群图，展示不同细胞类型中差异丰度邻域的分布 pt.pl.plot_da_beeswarm(mdata, alpha=0.1, palette=["#2E86AB", "#A23B72", "#F18F01", "#C73E1D", "#6A994E"]) # 自定义颜色 ``` 蜂群图（beeswarm plot）的Y轴是不同的细胞类型，X轴是logFC。每个点代表一个邻域，颜色代表其显著性。这张图能清晰地告诉你，哪些细胞类型中包含更多在处理后发生丰度变化的细胞亚群。 有时候，一个邻域可能由多种细胞类型混合构成。`pertpy`允许我们设定一个纯度阈值，将“混合型”邻域单独标记出来。 ```python # 检查邻域的细胞类型纯度分布 plt.hist(mdata["milo"].var["nhood_annotation_frac"], bins=30, edgecolor='k') plt.xlabel("Dominant celltype fraction in neighbourhood") plt.ylabel("Number of neighbourhoods") plt.axvline(x=0.6, color='r', linestyle='--', label='60% threshold') plt.legend() plt.show() # 将纯度低于60%的邻域标记为“Mixed” import pandas as pd mdata["milo"].var["nhood_annotation"] = mdata["milo"].var["nhood_annotation"].astype('category') mdata["milo"].var["nhood_annotation"] = mdata["milo"].var["nhood_annotation"].cat.add_categories("Mixed") mdata["milo"].var.loc[mdata["milo"].var["nhood_annotation_frac"] < 0.6, "nhood_annotation"] = "Mixed" # 重新绘制蜂群图，此时图例中会多出一个“Mixed”类别 pt.pl.plot_da_beeswarm(mdata, alpha=0.1, palette=["#2E86AB", "#A23B72", "#F18F01", "#C73E1D", "#6A994E", "#5D576B"]) # 为Mixed增加一个颜色 ``` 通过这个操作，你可以把那些细胞类型组成复杂、信号可能不明确的邻域区分开来，让结果更加清晰可靠。 ## 5. 实战调优与避坑指南 纸上得来终觉浅，绝知此事要躬行。在实际项目中应用`pertpy`进行`milo`分析，有几个关键的调优点和常见陷阱需要特别注意。 **参数调优策略：** `milo`分析对`KNN`图的构建和邻域采样参数比较敏感。我个人的经验是采用一种**网格搜索**式的探索策略。你可以固定其他步骤，系统性地调整`n_neighbors` (如30, 50, 100) 和 `prop` (如0.1, 0.2, 0.3)，然后观察： 1. 邻域大小分布是否合理（避免过多极小或极大的邻域）。 2. 显著邻域的数量和`logFC`的稳定性。理想情况下，核心的、强信号的差异邻域应该在多种参数下都保持显著。 下面是一个简单的参数扫描框架： ```python param_grid = {'n_neighbors': [30, 50, 80], 'prop': [0.1, 0.15, 0.2]} results_summary = [] for k in param_grid['n_neighbors']: for p in param_grid['prop']: print(f"Testing: k={k}, prop={p}") # 1. 重新构建KNN图 (注意：需要从原始adata重新开始，或深拷贝) adata_temp = adata.copy() sc.pp.neighbors(adata_temp, n_neighbors=k, n_pcs=30) # 2. 重新进行milo分析 milo_temp = pt.tl.Milo() mdata_temp = milo_temp.load(adata_temp) milo_temp.make_nhoods(mdata_temp["rna"], prop=p, seed=42) mdata_temp = milo_temp.count_nhoods(mdata_temp, sample_col="sample_id") da_res_temp = milo_temp.da_nhoods(mdata_temp, design='~ condition', design_df=design_df) # 3. 记录结果概要 n_sig = (da_res_temp['SpatialFDR'] < 0.1).sum() median_size = np.median(np.array(mdata_temp["rna"].obsm["nhoods"].sum(0)).ravel()) results_summary.append({ 'k': k, 'prop': p, 'n_significant': n_sig, 'median_nhood_size': median_size }) # 将结果转换为DataFrame便于查看 results_df = pd.DataFrame(results_summary) print(results_df.pivot(index='k', columns='prop', values='n_significant')) ``` **处理复杂实验设计：** 现实世界的实验设计往往比简单的“case vs control”复杂。你可能需要处理**配对样本**（如同一个病人治疗前后）、**多因素设计**（如不同处理+不同时间点）或**批次效应**。`pertpy`的`design`公式可以很好地处理这些情况。 例如，对于包含批次效应的设计： ```python # 假设 adata.obs 中有 'batch' 列 design_df = mdata["rna"].obs[['sample_id', 'condition', 'batch']].drop_duplicates().set_index('sample_id') design_df = design_df.loc[mdata["milo"].uns['sample_col_order'], :] # 在模型中加入批次作为协变量 da_results = milo.da_nhoods(mdata, design='~ batch + condition', design_df=design_df) ``` 对于配对设计（重复测量），你可能需要使用混合效应模型。虽然`pertpy`的底层`statsmodels`可能支持，但设置起来更复杂，需要仔细查阅文档并可能涉及自定义模型矩阵。 **结果稳定性验证：** 一次分析的结果可能受随机抽样（`make_nhoods`中的随机种子）影响。一个稳健的做法是**多次运行**（例如，用不同的随机种子运行5-10次），然后检查那些在多次运行中都被鉴定为显著的邻域。这些“共识”邻域的结果可信度更高。 ```python # 简单示例：运行两次，比较显著邻域的交集 seeds = [1234, 5678] sig_nhoods_list = [] for s in seeds: milo_obj = pt.tl.Milo() mdata_run = milo_obj.load(adata) sc.pp.neighbors(mdata_run["rna"], n_neighbors=50, n_pcs=30) milo_obj.make_nhoods(mdata_run["rna"], prop=0.2, seed=s) mdata_run = milo_obj.count_nhoods(mdata_run, sample_col="sample_id") da_res = milo_obj.da_nhoods(mdata_run, design='~ condition', design_df=design_df) sig_nhoods = set(da_res[da_res['SpatialFDR'] < 0.1]['Nhood']) sig_nhoods_list.append(sig_nhoods) # 计算Jaccard相似指数或简单交集 consistent_sig = set.intersection(*sig_nhoods_list) print(f"Number of consistently significant neighbourhoods across {len(seeds)} runs: {len(consistent_sig)}") ``` **与下游分析衔接：** 找到差异丰度邻域只是第一步。接下来你可能想知道： - **这些邻域有什么特征基因？** `pertpy`目前可能没有直接的内置函数，但你可以基于邻域索引，从原始的`adata`中提取这些细胞的表达矩阵，然后用`scanpy`的`rank_genes_groups`功能（指定`groups`为邻域ID）进行差异表达分析。 - **这些邻域是否对应于特定的功能通路？** 将上述找到的特征基因进行通路富集分析（如GO、KEGG）。 - **如何将milo结果与细胞轨迹分析结合？** 如果数据有连续的发育或分化轨迹，你可以将差异丰度邻域映射到轨迹上（例如`PAGA`图或`diffusion map`），观察丰度变化是否沿着特定分支发生。 最后，别忘了保存你的劳动成果。你可以保存完整的`MuData`对象，也可以只保存关键的统计结果。 ```python # 保存完整的分析对象（包含原始数据和所有中间结果） mdata.write('./my_milo_analysis.h5mu') # 保存差异丰度结果表格 da_results.to_csv('./da_results.csv', index=False) # 保存用于展示的图表 plt.savefig('./da_volcano.png', dpi=300, bbox_inches='tight') ``` 走完这一整套流程，从环境配置、数据准备、核心分析到结果解读和优化，你应该已经能在Python环境中独立驾驭`milo`进行差异丰度分析了。它不再是一个遥不可及的R包，而是你`scanpy`分析流程中一个自然延伸的环节。这种整合带来的流畅体验，能让你更专注于生物学问题的挖掘，而不是在工具链的切换中消耗精力。下次当你需要探究那些细微而连续的细胞群体变化时，不妨打开`pertpy`，试试这条新的路径。