专题｜Python 分析单细胞：差异表达分析
上一篇我们学习了 CellTypist 自动注释与 scVI 潜空间。本篇进入差异表达分析：注释完成后，哪些基因在不同细胞类型或不同条件之间差异表达？这是下游通路富集和功能分析的基础。
一、为什么需要差异表达分析
注释完成后，我们已经知道每个 cluster 是什么细胞类型。接下来自然会问：
哪些基因在 T 细胞中高表达，在 B 细胞中低表达？疾病组与对照组相比，哪些基因发生了变化？这些差异基因参与哪些通路？
差异表达分析（Differential Expression Analysis，DEG）回答这些问题。单细胞中的差异分析有两个层次：
层次比较对象示例cluster 级不同细胞类型之间T cell vs B cell 哪些基因不同condition 级同一细胞类型不同条件之间疾病 T cell vs 正常 T cell
两个层次的目的不同：cluster 级找的是 marker 基因，condition 级找的是条件响应基因。
记忆：单细胞差异分析有两个层次——cluster 级找 marker，condition 级找条件响应基因。不要混为一谈。
二、rank_genes_groups 详解
1. 基本运行
sc.tl.rank_genes_groups( adata, groupby=”cell_type”, method=”wilcoxon”, use_raw=True, pts=True, key_added=”deg”,)
参数含义常用设置groupby按哪一列分组"leiden" 或 "cell_type"groups指定分析哪些组默认所有组reference参考组默认 "rest"，也可指定某一组method差异检验方法"wilcoxon"、"t-test"、"logreg"use_raw是否使用 adata.rawTrue 用全基因layer指定使用哪个 layer如 "lognorm"pts是否计算表达比例Truekey_added结果存储到 uns 的键名默认 "rank_genes_groups"n_genes保留多少基因默认全部
2. 三种方法对比
方法原理优点缺点"t-test"t 检验比较均值速度快对分布假设强，不稳健"wilcoxon"Wilcoxon 秩和检验稳健，不需分布假设速度中等"logreg"逻辑回归分类找区分性基因，可加协变量结果不直接是 p 值
推荐初学者优先使用 "wilcoxon"，它是最常用的稳健方法。
3. 提取结果到 DataFrame
deg_df = sc.get.rank_genes_groups_df( adata, group=None, key=”deg”,)print(deg_df.head(15))
输出类似：
group gene scores logfoldchanges pvals pvals_adj pct_nz_group pct_nz_reference0 T cell CD3D 15.230 3.450 1.20e-35 1.20e-33 0.95 0.021 T cell CD3E 14.870 3.120 3.45e-33 3.45e-31 0.92 0.012 T cell CD8A 13.560 2.890 8.90e-30 8.90e-28 0.88 0.033 B cell MS4A1 14.120 4.010 5.67e-32 5.67e-30 0.91 0.014 B cell CD79A 13.450 3.230 1.23e-29 1.23e-27 0.89 0.02
各列含义：
列名含义group分组名gene基因名scores得分（方法相关）logfoldchangeslog 倍数变化pvals原始 p 值pvals_adj校正后 p 值（FDR）pct_nz_group该组中表达比例pct_nz_reference参考组中表达比例
记忆：rank_genes_groups 是 Scanpy 差异分析的核心函数。结果提取到 DataFrame 后，pvals_adj 和 logfoldchanges 是筛选显著基因的两个关键列。
三、指定 cluster 之间比较
除了默认的"每组 vs 其余"，还可以指定两个特定组之间比较。
1. T cell vs B cell
sc.tl.rank_genes_groups( adata, groupby=”cell_type”, groups=[”T cell”], reference=”B cell”, method=”wilcoxon”, key_added=”T_vs_B”,)deg_T_vs_B = sc.get.rank_genes_groups_df( adata, group=”T cell”, key=”T_vs_B”,)print(deg_T_vs_B.head(20))
输出类似：
gene scores logfoldchanges pvals pvals_adj0 CD3D 18.450 4.120 5.60e-40 5.60e-381 CD3E 17.230 3.780 2.30e-36 2.30e-342 CD8A 16.560 3.450 1.10e-33 1.10e-313 TRAC 15.890 3.230 4.50e-31 4.50e-294 CD4 14.780 2.980 8.90e-28 8.90e-26
2. 可视化
sc.pl.rank_genes_groups( adata, groups=[”T cell”], n_genes=20, key=”T_vs_B”,)
这种两两比较方式更精确，适合需要知道特定两组之间差异的场景。
四、五种 DEG 可视化方法
1. rank_genes_groups 点图
sc.pl.rank_genes_groups( adata, n_genes=25, sharey=False, key=”deg”,)
参数含义n_genes每组展示前几个基因sharey是否共享 y 轴groups指定展示哪些组key对应 rank_genes_groups 的 key_added
2. Dotplot
top_genes = {}for ct in adata.obs[”cell_type”].unique(): top_genes[ct] = sc.get.rank_genes_groups_df( adata, group=ct, key=”deg” ).head(5)[”gene”].tolist()sc.pl.dotplot( adata, top_genes, groupby=”cell_type”, standard_scale=”var”, dendrogram=True,)
参数含义var_names基因字典，键为细胞类型groupby分组列standard_scale标准化方式dendrogram是否显示聚类树
3. Stacked violin
sc.pl.stacked_violin( adata, top_genes, groupby=”cell_type”, dendrogram=True,)
参数含义var_names基因字典groupby分组列dendrogram聚类树swap_axes交换行列
4. Heatmap
sc.pl.rank_genes_groups_heatmap( adata, n_genes=10, groupby=”cell_type”, key=”deg”, standard_scale=”var”, show_gene_labels=True,)
参数含义n_genes每组取多少基因groupby分组列standard_scale标准化方式show_gene_labels是否显示基因名key结果键名
5. Matrixplot
sc.pl.rank_genes_groups_matrixplot( adata, n_genes=5, groupby=”cell_type”, key=”deg”, standard_scale=”var”, cmap=”Blues”,)
参数含义n_genes每组基因数groupby分组列standard_scale标准化方式cmap颜色映射key结果键名
记忆：五种 DEG 可视化方法与注释篇相同，但这里用的是 rank_genes_groups 结果的 top 基因，而不是预定义 marker。两者结合注释更可靠。
五、Condition 级差异分析
1. 过滤到一种细胞类型
t_cells = adata[adata.obs[”cell_type”] == ”T cell”].copy()print(t_cells.obs[”condition”].value_counts())
输出类似：
Control 1600Disease 1600dtype: int64
2. 比较条件
sc.tl.rank_genes_groups( t_cells, groupby=”condition”, groups=[”Disease”], reference=”Control”, method=”wilcoxon”, key_added=”disease_vs_control”,)deg_condition = sc.get.rank_genes_groups_df( t_cells, group=”Disease”, key=”disease_vs_control”,)print(deg_condition.head(20))
输出类似：
gene scores logfoldchanges pvals pvals_adj0 ISG15 12.340 2.560 3.40e-20 3.40e-181 IFIT1 11.780 2.340 1.20e-18 1.20e-162 IFIT3 11.230 2.190 5.60e-17 5.60e-153 MX1 10.890 2.050 2.30e-15 2.30e-134 OAS1 10.450 1.920 8.90e-14 8.90e-12
3. 可视化
sc.pl.rank_genes_groups( t_cells, groups=[”Disease”], n_genes=20, key=”disease_vs_control”,)
条件级 DEG 是后续通路富集分析（GSEA/GO）的直接输入。
记忆：condition 级 DEG 必须先过滤到一种细胞类型，再比较条件。不能把所有细胞混在一起比较条件，否则差异来自细胞类型组成而非条件。
六、Pseudo-bulk 分析
1. 为什么需要 pseudo-bulk
单细胞数据把每个细胞当作独立样本做统计检验，这会夸大显著性——因为同一样本内的细胞不独立（批次效应、遗传背景相同）。正确做法是把每个样本的同一细胞类型聚合成 pseudo-bulk，再在样本级别做差异分析。
方法统计单位样本量优点缺点单细胞级每个细胞大（数千）灵敏度高夸大 p 值Pseudo-bulk每个样本×细胞类型小（几个到十几个）统计正确需要多样本
2. 使用 decoupler 生成 pseudo-bulk
import decoupler as dcpdata = dc.get_pseudobulk( adata, sample_col=”sample_id”, groups_col=”cell_type”, layer=”counts”, mode=”sum”, min_cells=10, min_counts=1000,)
参数含义常用设置sample_col样本 ID 列名"sample_id"groups_col细胞类型列名"cell_type"layer使用哪个 layer 的数据"counts"mode聚合方式"sum" 或 "mean"min_cells每个样本×类型最少细胞数10min_counts每个样本×类型最少 counts1000
输出类似：
AnnData object with n_obs × n_vars = 24 × 20000 obs: sample_id, cell_type, ...
每个 obs 行是一个样本×细胞类型的聚合体。
3. 过滤到 T 细胞并做条件比较
pdata_t = pdata[pdata.obs[”cell_type”] == ”T cell”].copy()print(pdata_t.obs[[”sample_id”, ”condition”]])
输出类似：
sample_id condition0 S01 Control1 S02 Control2 S03 Control3 S04 Disease4 S05 Disease5 S06 Disease
4. 运行差异分析
# 标准化和 log 转换sc.pp.normalize_total(pdata_t, target_sum=1e4)sc.pp.log1p(pdata_t)# 差异分析sc.tl.rank_genes_groups( pdata_t, groupby=”condition”, groups=[”Disease”], reference=”Control”, method=”t-test”, key_added=”pb_deg”,)pb_deg = sc.get.rank_genes_groups_df( pdata_t, group=”Disease”, key=”pb_deg”,)print(pb_deg.head(20))
记忆：pseudo-bulk 把每个样本×细胞类型聚合为一个统计单位，避免了把同一样本的细胞当独立重复的伪重复问题。多样本研究必须做 pseudo-bulk。
七、导出 GSEA 排序表
GSEA 需要 logFC 排序的基因列表。
1. 生成排序表
# 按 logfoldchanges 排序gsea_input = deg_condition.sort_values( ”logfoldchanges”, ascending=False)# 只保留基因名和 logFCranking = gsea_input[[”gene”, ”logfoldchanges”]].drop_duplicates( subset=”gene”)print(ranking.head(10))
输出类似：
gene logfoldchanges0 ISG15 2.5601 IFIT1 2.3402 IFIT3 2.1903 MX1 2.0504 OAS1 1.9205 IFITM1 1.8506 IFITM3 1.7807 ISG20 1.6908 RSAD2 1.6209 IFI44 1.560
2. 保存 CSV
ranking.to_csv( ”t_cell_disease_vs_control_ranking.csv”, index=False, header=False,)
这个 CSV 文件可以直接作为 GSEApy 的输入。注意：GSEA 需要从高到低的完整排序列表，不是只取显著基因。
记忆：GSEA 需要按 logFC 从高到低排序的完整基因列表，不要截断。只取显著基因做 ORA（过度表征分析），不是 GSEA。
八、多重检验校正与阈值
1. FDR 校正
rank_genes_groups 默认返回 pvals_adj，这是经过 Benjamini-Hochberg 校正的 FDR。筛选标准通常：
significant = deg_df[ (deg_df[”pvals_adj”] < 0.05) & (deg_df[”logfoldchanges”].abs() > 0.5)]print(f”显著基因数: {len(significant)}”)
输出类似：
2. 阈值选择参考
指标常用阈值说明pvals_adj< 0.05FDR 显著性logfoldchanges> 0.5 或 > 1log2 倍数变化pct_nz_group> 0.25在组中表达细胞比例pct_nz_reference< 0.25在参考组中表达比例
# 更严格的筛选strict = deg_df[ (deg_df[”pvals_adj”] < 0.01) & (deg_df[”logfoldchanges”].abs() > 1) & (deg_df[”pct_nz_group”] > 0.25)]
九、火山图
火山图同时展示 logFC 和 p 值，直观看到显著基因分布。
1. 绘制火山图
import matplotlib.pyplot as pltimport numpy as npdeg = deg_condition.copy()deg[”-log10_padj”] = -np.log10(deg[”pvals_adj”])fig, ax = plt.subplots(figsize=(8, 6))# 不显著灰色ax.scatter( deg.loc[deg[”pvals_adj”] >= 0.05, ”logfoldchanges”], deg.loc[deg[”pvals_adj”] >= 0.05, ”-log10_padj”], c=”grey”, s=5, alpha=0.5, label=”NS”)# 显著上调红色up = deg[(deg[”pvals_adj”] < 0.05) & (deg[”logfoldchanges”] > 0.5)]ax.scatter(up[”logfoldchanges”], up[”-log10_padj”], c=”red”, s=5, label=”Up”)# 显著下调蓝色down = deg[(deg[”pvals_adj”] < 0.05) & (deg[”logfoldchanges”] < -0.5)]ax.scatter(down[”logfoldchanges”], down[”-log10_padj”], c=”blue”, s=5, label=”Down”)# 标注 top 基因top_genes = up.nlargest(10, ”-log10_padj”)for _, row in top_genes.iterrows(): ax.annotate(row[”gene”], (row[”logfoldchanges”], row[”-log10_padj”]), fontsize=8, ha=”left”)ax.set_xlabel(”log2 Fold Change”)ax.set_ylabel(”-log10(FDR)”)ax.set_title(”T cell: Disease vs Control”)ax.legend()plt.tight_layout()plt.savefig(”volcano.png”, dpi=300)plt.show()
输出类似：
火山图左侧是下调基因（蓝色），右侧是上调基因（红色），灰色是不显著基因。纵轴越高表示越显著。
记忆：火山图横轴是 logFC（方向），纵轴是 -log10(FDR)（显著性）。右上角是最值得关注的上调显著基因。
十、保存结果
# 保存 DEG 表deg_df.to_csv(”deg_by_cell_type.csv”, index=False)# 保存条件级 DEGdeg_condition.to_csv(”t_cell_disease_vs_control.csv”, index=False)# 保存 pseudo-bulk DEGpb_deg.to_csv(”pseudobulk_t_cell_disease_vs_control.csv”, index=False)# 保存 GSEA 排序表ranking.to_csv(”gsea_ranking.csv”, index=False, header=False)# 保存对象adata.write_h5ad(”adata_with_deg.h5ad”)t_cells.write_h5ad(”t_cells_with_condition_deg.h5ad”)
十一、常见错误
现象常见原因回退步骤所有基因都显著统计单位是细胞而非样本，p 值被夸大改用 pseudo-bulk 分析logFC 全是 0使用了 scale 后的数据改用 use_raw=True 或 layer="lognorm"condition 级 DEG 全是细胞类型 marker没有先过滤到一种细胞类型先 adata[adata.obs["cell_type"]=="T cell"] 再比较条件pseudo-bulk 样本太少样本数不足 3 vs 3增加样本或报告局限性GSEA 结果为空排序表被截断或格式不对使用完整排序，确认格式为两列无表头火山图全是灰色p 值校正过严或 logFC 阈值过高放松 FDR 到 0.1 或 logFC 到 0.25不同方法结果差异大t-test 和 wilcoxon 假设不同优先信任 wilcoxon，报告方法敏感性
十二、完整流程代码
import scanpy as scimport pandas as pdimport numpy as npimport matplotlib.pyplot as pltimport decoupler as dc# --------------------------------------------------# 0. 前提：已完成注释，adata.obs 有 cell_type 列，# adata.obs 有 condition/sample_id 列# --------------------------------------------------# 1. cluster 级 DEG（marker 基因）sc.tl.rank_genes_groups( adata, groupby=”cell_type”, method=”wilcoxon”, use_raw=True, pts=True, key_added=”deg”,)deg_df = sc.get.rank_genes_groups_df(adata, group=None, key=”deg”)print(deg_df.head(20))# 2. 可视化sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False, key=”deg”)# 3. 提取 top 基因画 Dotplottop_genes = {}for ct in adata.obs[”cell_type”].unique(): df_ct = sc.get.rank_genes_groups_df(adata, group=ct, key=”deg”) top_genes[ct] = df_ct.head(5)[”gene”].tolist()sc.pl.dotplot(adata, top_genes, groupby=”cell_type”, standard_scale=”var”, dendrogram=True)sc.pl.stacked_violin(adata, top_genes, groupby=”cell_type”, dendrogram=True)sc.pl.rank_genes_groups_heatmap(adata, n_genes=10, groupby=”cell_type”, key=”deg”, standard_scale=”var”)sc.pl.rank_genes_groups_matrixplot(adata, n_genes=5, groupby=”cell_type”, key=”deg”, standard_scale=”var”, cmap=”Blues”)# 4. 两两比较：T cell vs B cellsc.tl.rank_genes_groups( adata, groupby=”cell_type”, groups=[”T cell”], reference=”B cell”, method=”wilcoxon”, key_added=”T_vs_B”,)deg_T_vs_B = sc.get.rank_genes_groups_df(adata, group=”T cell”, key=”T_vs_B”)print(deg_T_vs_B.head(20))# 5. condition 级 DEGt_cells = adata[adata.obs[”cell_type”] == ”T cell”].copy()sc.tl.rank_genes_groups( t_cells, groupby=”condition”, groups=[”Disease”], reference=”Control”, method=”wilcoxon”, key_added=”disease_vs_control”,)deg_condition = sc.get.rank_genes_groups_df( t_cells, group=”Disease”, key=”disease_vs_control”)print(deg_condition.head(20))# 6. 筛选显著基因significant = deg_condition[ (deg_condition[”pvals_adj”] < 0.05) & (deg_condition[”logfoldchanges”].abs() > 0.5)]print(f”显著基因数: {len(significant)}”)# 7. Pseudo-bulkpdata = dc.get_pseudobulk( adata, sample_col=”sample_id”, groups_col=”cell_type”, layer=”counts”, mode=”sum”, min_cells=10, min_counts=1000,)pdata_t = pdata[pdata.obs[”cell_type”] == ”T cell”].copy()sc.pp.normalize_total(pdata_t, target_sum=1e4)sc.pp.log1p(pdata_t)sc.tl.rank_genes_groups( pdata_t, groupby=”condition”, groups=[”Disease”], reference=”Control”, method=”t-test”, key_added=”pb_deg”,)pb_deg = sc.get.rank_genes_groups_df(pdata_t, group=”Disease”, key=”pb_deg”)print(pb_deg.head(20))# 8. GSEA 排序表gsea_input = deg_condition.sort_values(”logfoldchanges”, ascending=False)ranking = gsea_input[[”gene”, ”logfoldchanges”]].drop_duplicates(subset=”gene”)ranking.to_csv(”gsea_ranking.csv”, index=False, header=False)# 9. 火山图deg = deg_condition.copy()deg[”-log10_padj”] = -np.log10(deg[”pvals_adj”])fig, ax = plt.subplots(figsize=(8, 6))ax.scatter( deg.loc[deg[”pvals_adj”] >= 0.05, ”logfoldchanges”], deg.loc[deg[”pvals_adj”] >= 0.05, ”-log10_padj”], c=”grey”, s=5, alpha=0.5, label=”NS”)up = deg[(deg[”pvals_adj”] < 0.05) & (deg[”logfoldchanges”] > 0.5)]ax.scatter(up[”logfoldchanges”], up[”-log10_padj”], c=”red”, s=5, label=”Up”)down = deg[(deg[”pvals_adj”] < 0.05) & (deg[”logfoldchanges”] < -0.5)]ax.scatter(down[”logfoldchanges”], down[”-log10_padj”], c=”blue”, s=5, label=”Down”)top = up.nlargest(10, ”-log10_padj”)for _, row in top.iterrows(): ax.annotate(row[”gene”], (row[”logfoldchanges”], row[”-log10_padj”]), fontsize=8)ax.set_xlabel(”log2 Fold Change”)ax.set_ylabel(”-log10(FDR)”)ax.set_title(”T cell: Disease vs Control”)ax.legend()plt.tight_layout()plt.savefig(”volcano.png”, dpi=300)plt.show()# 10. 保存deg_df.to_csv(”deg_by_cell_type.csv”, index=False)deg_condition.to_csv(”t_cell_disease_vs_control.csv”, index=False)pb_deg.to_csv(”pseudobulk_t_cell_disease_vs_control.csv”, index=False)adata.write_h5ad(”adata_with_deg.h5ad”)
十三、本篇小结
差异表达分析分两个层次：cluster 级用 rank_genes_groups 找 marker，condition 级需先过滤到一种细胞类型再比较条件。多样本研究必须用 pseudo-bulk 避免伪重复。GSEA 需要按 logFC 排序的完整基因列表。多重检验校正用 pvals_adj < 0.05，配合 logFC > 0.5 筛选显著基因。火山图直观展示差异分布。
下一篇我们将学习 scCODA 样本级细胞组成分析——把细胞标签正确汇总到独立样本，在组成约束下解释相对变化。
10个月宝宝每天需要喝多少奶粉？
