企业官网建设流程全解析

1. 为什么单细胞差异表达分析没那么简单

1.1 差异表达分析到底要回答什么问题

拿到一张单细胞转录组图谱，最常被问到的不是“你聚类分了几个群”，而是“这群细胞和那群细胞到底哪里不一样”“处理组的巨噬细胞和对照组的巨噬细胞在基因表达上有哪些差异”。这些问题的答案，基本都落在差异表达基因的清单里。

Seurat包的FindMarkers函数，就是专门用来做这件事的工具。你可以把它理解为单细胞版的t检验或者Wilcoxon检验，但它面对的并不是常规的样本间比较，而是已经在降维聚类后定义好的细胞亚群之间、或者同一类细胞在不同条件之间的比较。它试图回答：给定两组细胞，各组的基因表达数据都存在一个分布，那么哪些基因在两个分布之间差异足够显著，并且差异的幅度足够大？

这个“足够”两个字，正是整个分析中最容易产生分歧的地方。不同文章、不同项目、不同生物学背景下，差异基因的标准可以差出好几倍。有人只看p_val_adj < 0.05，有人要求avg_log2FC > 1，有人还额外叠加了min.pct的过滤。其实这些都不是FindMarkers的问题，而是使用者在定义“差异”时的选择问题。

1.2 单细胞数据与bulk RNA-seq的三个本质差异

很多从bulk转录组转来做单细胞的人，最开始会用老思维去套。但单细胞差异表达分析有三个非常反直觉的地方，必须在动手前想清楚。

第一是dropout事件。单细胞测序受限于捕获效率和测序深度，很多低表达基因在某些细胞里检测到的表达量为0，但并不代表这个基因真的没有表达。这导致单细胞表达矩阵特别稀疏，0值比例很高。如果在做差异分析时把0当作“真阴性”，很容易高估差异显著性。

第二是样本量问题。单细胞数据里的“样本”不是指生物学重复的个体数，而是指参与比较的细胞数。一个cluster可以有3000个细胞，另一个条件只有50个细胞，这种数量不均衡会让统计检验的结果产生偏差。更重要的是，如果这50个细胞全部来自同一个个体，那它们并不是50个独立样本，而是1个样本里的50个观测值。统计学上的伪重复问题，在单细胞差异表达分析中非常普遍，很多人会忽略。

第三是生物学重复缺失。如果只有一个对照组样本和一个处理组样本，即便两个样本中各有5000个细胞，FindMarkers的p值也仅仅反映细胞间的差异，无法代表个体间的差异。这一点决定了后续是否能够当真，是否经得起验证。

这些特性决定了FindMarkers不是简单地套一个t.test就行，而是需要仔细设置参数、理解检验假设，甚至需要根据实验设计换用不同的统计模型。

2. FindMarkers核心参数与底层逻辑

2.1 test.use参数：五种检验方式到底怎么选

Seurat的FindMarkers函数最关键的参数就是test.use，它决定你用哪种统计检验来计算p值。默认值是wilcox，也就是Wilcoxon秩和检验。

wilcox是首选，但没有被过分神话。它的优点是几乎不做分布假设，对0值比较多的高稀疏数据相对稳健，而且运算速度快。配合Seurat内部的presto优化，在几万个细胞上跑FindAllMarkers也不会慢到让人崩溃。

如果你的数据量不大、只是想快速拿到一组候选标记基因，直接使用test.use = "wilcox"是完全没有问题的。

如果数据来自多个样本、多个个体，尤其存在配对或者重复测量设计时，可以选择MAST。MAST是一个基于混合模型的检验方法，它把表达量拆分成“是否表达（离散部分）”和“表达量高低（连续部分）”两个过程，对dropout的建模比较友好，因此也经常被用于单细胞差异表达分析。但MAST的运行速度明显比wilcox慢一个量级，在5万细胞级别跑一遍要等不少时间。

DESeq2和edgeR这类方法原本是为bulk RNA-seq设计的，内部会估算离散度、做负二项分布拟合，它们需要原始的整数count数据，而且对每个“样本”的定义通常是每个细胞或每个生物学重复。如果你拿它跑单细胞数据，建议先将细胞汇总成伪bulk样本，否则会把大量细胞的表达当作独立样本，得出极其乐观的假阳性p值。

还有bimod和roc两个选项。bimod是Seurat早期版本常用的似然比检验，roc则是基于分类准确性的AUC方法，不是真正的p值假设检验，通常用它的AUC值来辅助判断差异方向。日常推荐直接在wilcox、MAST、DESeq2里选，其他不太必要。

2.2 min.pct、logfc.threshold与only.pos的取舍逻辑

很多新手会困惑：已经得到p值了，为什么还需要min.pct和logfc.threshold？

min.pct表达的是“这个基因在两个比较组中，至少在多少比例的细胞里检测到表达”。默认值是0.1，也就是说一个基因必须在至少一个组里有10%以上的细胞表达，才参与检验。这个参数非常实用，它把那些只在极少数细胞中表达的基因排除掉，否则很多基因因为个数太少，方差太小，p值会异常显著，实际却没有生物学意义。

我在实际项目中通常会把min.pct设为0.25或0.5，尤其是在寻找细胞类型标志基因的时候。因为一个靠谱的marker基因应当在一个亚群中有较广泛的表达，而不是两三颗细胞里的偶然信号。

logfc.threshold则是差异幅度的门槛，默认是0.25，表示两组平均表达量的log2倍数变化至少需要达到0.25，才进入最终结果。这个值设得太小，结果表会变得很长，大量生物学意义微弱的基因混进来；设得太大，又会漏掉很多真实但变化幅度不大的基因。

很多用户会把logfc.threshold = 0来获取全部检验结果，这不是不行，但你需要明白这样会拿到大量统计显著但fold change很小的基因。与其事后拼命过滤，不如一开始就想清楚自己的分析目标。如果你要的是后续做功能富集，幅度稍低的基因保留一些反而更全面；如果你要的是拿去做细胞类型注释的marker，那logfc.threshold至少得给到0.5甚至1。

only.pos参数比较简单，它只影响返回结果的方向。设成TRUE时，只返回在ident.1组中上调的基因。寻找细胞类型marker时一般会设置only.pos = TRUE，因为描述一个细胞类型，只需要它“相对于其他类型有哪些基因高表达”，极少关心它在哪些基因上低表达。

2.3 Seurat的标准化与p值背后的表达矩阵

在用FindMarkers之前，必须明确它默认使用的数据矩阵是哪个。

Seurat对象里有三个常见的表达矩阵：RNA的counts（原始UMI/read数）、data（log标准化后的数据）、scale.data（中心化并缩放到方差为1的数据）。FindMarkers默认使用的slot是data，也就是log1p(CPM/10000)变换后的结果。它不会自动使用scale.data做差异分析，原因是scale.data已经把所有基因的均值变成了0，基因之间的绝对表达差异被抹掉了，没法用来计算倍数变化。

这里的“表达量”是指每个细胞在某个基因上的log标准化值。通过比较两组细胞在这个值上的分布，计算平均值的差，得到avg_log2FC；再通过秩和检验或模型检验得到p值。所以千万不要试图把FindMarkers输出里的avg_log2FC直接拿原始count矩阵里的均值相减再取log2来验证，因为计算基准不同。

FindMarkers内部在比较前会根据min.pct过滤掉一部分低表达基因。这样既加快了计算，也避免了零膨胀对统计的干扰。如果你发现某个很感兴趣的基因没有出现在结果表里，先看看它是不是被min.pct卡掉了，再去看看它在目标群体里的表达阳性率到底有多少。

3. 实际操作：从Seurat对象到差异基因列表

3.1 准备数据与Seurat对象构建

为了演示，我们假设已经拿到了一个10X格式的单细胞转录组数据，包括三个样本：两个对照组和两个处理组。我的演示环境是R 4.3，Seurat 5，需先确保已经安装Seurat及相关包。

library(Seurat) library(dplyr) library(ggplot2) # 读取10X数据，生成对象 ctrl1 <- Read10X("data/ctrl1/filtered_feature_bc_matrix/") ctrl2 <- Read10X("data/ctrl2/filtered_feature_bc_matrix/") treat1 <- Read10X("data/treat1/filtered_feature_bc_matrix/") treat2 <- Read10X("data/treat2/filtered_feature_bc_matrix/") # 合并成表达矩阵，并打上样本标签 library(Seurat) obj <- CreateSeuratObject(counts = ctrl1, project = "ctrl1", min.cells = 3, min.features = 200) obj[["RNA"]] <- JoinLayers(obj[["RNA"]])

上面只是一个示例。实际项目中通常需要用merge或者CreateSeuratObject后逐个添加，更推荐直接使用Read10X得到列表，merge合并，再加样本标签。单细胞分析的标准流程是：质控、标准化、高变基因、PCA、UMAP、聚类、注释。等聚类和注释完成之后，我们才能定义Idents，进而执行FindMarkers。

我先给出一个简化的完整流程：

# 质控 obj[["percent.mt"]] <- PercentageFeatureSet(obj, pattern = "^MT-") obj <- subset(obj, subset = nFeature_RNA > 200 & nFeature_RNA < 5000 & percent.mt < 20) # 标准化 obj <- NormalizeData(obj, normalization.method = "LogNormalize", scale.factor = 10000) # 高变基因与降维 obj <- FindVariableFeatures(obj, selection.method = "vst", nfeatures = 2000) obj <- ScaleData(obj) obj <- RunPCA(obj, npcs = 30) obj <- FindNeighbors(obj, dims = 1:20) obj <- FindClusters(obj, resolution = 0.5) obj <- RunUMAP(obj, dims = 1:20)

3.2 基本用法：两个指定细胞群之间找差异基因

假设我们在UMAP聚类后获得了0到10一共11个cluster。现在想比较cluster 2和cluster 5这两个亚群的转录差异。这是FindMarkers最基础的用法。

# 默认Idents是聚类结果 Idents(obj) <- "seurat_clusters" markers_cluster2_vs_5 <- FindMarkers(obj, ident.1 = "2", ident.2 = "5", test.use = "wilcox", min.pct = 0.25, logfc.threshold = 0.25, only.pos = FALSE)

这里ident.1是处理组，ident.2是对照组。返回的表格中，每个基因一行，列包括p_val、avg_log2FC、pct.1、pct.2、p_val_adj。

avg_log2FC是ident.1相对于ident.2的平均表达差异。正值表示在ident.1（这里是cluster 2）中高表达，负值表示在cluster 5中高表达。这个方向很多人搞反，我一般习惯先跑一个已知marker基因验证一下方向是否正确，再正式解读。

3.3 批量查找：FindAllMarkers遍历所有cluster

如果你刚做完聚类，还不清楚每个cluster是什么细胞类型，想一次性把每个cluster相对其他所有cluster的marker基因找出来，那么用FindAllMarkers。

all_markers <- FindAllMarkers(obj, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.5, test.use = "wilcox")

这个函数会对每个cluster都执行一次“该cluster vs 剩余所有细胞”的差异分析。输出结果包含一列cluster，表示该基因是在哪个cluster中被识别为上调marker。注意它的比较组规模很不平衡：一个500细胞的cluster在跟10000细胞的其他所有cluster比较时，统计效力会有偏差，所以p_val_adj虽然矫正了多重检验，但没有矫正细胞数量的不均衡。作为初步筛选没问题，但不能直接当成严格的组间差异。

3.4 分组比较：同一细胞类型在不同条件间的差异

更常见的医学场景是：我想知道处理组和对照组之间，某一类细胞（比如CD8+ T细胞）的转录状态发生了什么变化。

这要求先把目标细胞给分离出来，然后设置新的Idents为样本分组信息。

# 假设metadata中已有 group 列: "ctrl" 和 "treat" cd8 <- subset(obj, celltype == "CD8_T") Idents(cd8) <- "group" # 比较处理组 vs 对照组 markers_treat_vs_ctrl <- FindMarkers(cd8, ident.1 = "treat", ident.2 = "ctrl", min.pct = 0.1, logfc.threshold = 0.25, test.use = "wilcox")

注意，这里subset之后，细胞数量可能变得很少。如果只有几十个细胞，找出来的差异基因会很飘，建议至少目标组有100个以上细胞再跑。另外，如果治疗组和对照组来自多个个体，并且个体间差异很大，wilcox检验很容易把个体间的随机差异当作处理效应。此时最好换成MAST并加入个体作为随机效应，或者用Seurat::Pseudobulk先做伪bulk再跑DESeq2。关于这一点，我后面在常见问题部分会展开。

4. 结果解读与可视化

4.1 结果表的字段到底该怎么看

FindMarkers返回的是一个data.frame，行名是基因名。下面是一段示例输出：

p_val是原始p值，p_val_adj是Bonferroni校正后的p值，它直接用p值乘以检测的基因数。Seurat用的是Bonferroni校正，非常严格。虽然保守，但在差异基因数目不多的时候很好用。

avg_log2FC是两组的log2平均表达倍数变化，正值表示ident.1高表达，负值表示ident.2高表达。pct.1是基因在ident.1组中表达的细胞百分比，pct.2是ident.2组中的百分比。

一个常见误区是直接看p_val排序。p值只是统计显著性的度量，如果两组细胞数都很大，即使表达差异很小也会出现极小的p值。所以一定要先看avg_log2FC是否达到生物学上有意义的幅度，再看阳性率是否支持这个差异。

4.2 筛选差异基因的实操标准

在应用层面，我经常用一组组合条件来筛选“可靠差异基因”：

• p_val_adj < 0.05
• abs(avg_log2FC) > 0.5（如果要更严格，可设>1）
• pct.1 > 0.2或pct.2 > 0.2（至少在一个组中有一定表达阳性率）
• 允许pct.1与pct.2之间有明显差别，否则即使p值很小，也可能只是表达量微调，不构成类型marker

筛选代码可以这样写：

sig_markers <- markers_treat_vs_ctrl %>% filter(p_val_adj < 0.05, abs(avg_log2FC) > 0.5) %>% arrange(desc(avg_log2FC))

在细胞类型注释时，我还会再叠加一层“高表达基因特异性”的判断：一个marker基因应当在一个类型中高表达，在其他类型中低表达或不表达。此时pct.1 - pct.2的差值比单纯的平均表达量更有参考价值。

4.3 用火山图和热图把差异基因讲清楚

拿到差异基因表，下一步通常是可视化。火山图展示全局差异分布，热图展示目标基因表达模式。

library(ggplot2) library(ggrepel) # 准备带标签的数据 plot_data <- markers_treat_vs_ctrl %>% mutate(threshold = case_when( p_val_adj < 0.05 & avg_log2FC > 0.5 ~ "up", p_val_adj < 0.05 & avg_log2FC < -0.5 ~ "down", TRUE ~ "ns" )) ggplot(plot_data, aes(x = avg_log2FC, y = -log10(p_val_adj), color = threshold)) + geom_point(size = 0.5) + scale_color_manual(values = c("up" = "#E64B35", "down" = "#4DBBD5", "ns" = "grey70")) + theme_minimal() + labs(x = "log2 Fold Change (treat vs ctrl)", y = "-log10 adjusted p-value")

热图一般挑选每个cluster的前10或前20个marker基因：

top_genes <- all_markers %>% group_by(cluster) %>% slice_max(n = 10, order_by = avg_log2FC) DoHeatmap(obj, features = top_genes$gene, slot = "data")

这里提醒一下，DoHeatmap用的表达值默认是dataslot。如果细胞数量很多，画热图会比较卡。可以先对每个cluster做平均表达，再画热图，会更轻量、更容易看出梯度。

5. 常见问题与排查技巧实录

5.1 检查ident.1与ident.2的方向错误

FindMarkers的结果里，avg_log2FC的正负完全取决于ident.1和ident.2的设定。有一次我帮别人排查，对方怎么都找不到一个已知的marker基因，后来发现他ident.1写的是对照，ident.2写的是处理，结果up和down完全反了，然后他还把结果直接拿去跟文献比较，自然对不上。

建议在分析前先用levels(Idents(obj))查看当前ident水平的顺序和名称，再确认你要比较的两个组在顺序中处于什么位置。如果写错，最简单的办法是重新跑一遍，或者在代码里标记清楚注释：ident.1是处理，ident.2是对照，便于后续检查。

5.2 某些基因在部分细胞中不表达，导致p值虚高

单细胞数据0值太多，这对秩和检验的影响很大。假设一个基因在对照组几乎不表达，在处理组也只有两个细胞高表达，那Wilcoxon检验会认为这两组的分布差异极大，给出极小的p值。但这两个高表达细胞是真正的生物学信号，还是扩增伪影，无从判断。

因此，min.pct非常重要。如果你发现某个基因的两个组的阳性率分别是0.5%和1%，即便p值再小，我也不会把它列为候选基因，因为阳性率太低了，非常不稳定。最好把min.pct提高到至少0.1，或者对阳性率差额设置一个阈值。想要更稳健，还可以设置min.cells.group参数，确保每个组的表达细胞数量足够。

5.3 细胞数量悬殊时的统计偏差

比较两个cluster时，一个cluster有2000个细胞，另一个只有80个细胞，80个细胞那组的方差往往被低估，由此产生的p值会偏小。这不是Seurat独有的问题，是所有基于细胞粒度的检验都存在的。

如果你遇到这种数量悬殊的情况，可以先看看差异基因是否大多来自大细胞群一边。如果是，建议用另一种策略：将每个cluster的细胞随机下采样到较小数量，重复跑多次，看核心差异基因是否稳定出现；或者直接用伪bulk分析，规避细胞数量差异导致的假阳性。

5.4 avg_log2FC为什么和手动计算不一致

很多用户会手动用AverageExpression计算两组平均表达量，然后相减取log2，发现结果跟FindMarkers输出的avg_log2FC对不上。这是正常的。

Seurat的avg_log2FC并不是直接把两组平均log表达相减，而是在过去的某些版本中，使用的是所有细胞log标准化表达值的平均差（需要乘以log2因子），此外还会考虑数据的自然对数和log2的换算。不同版本的Seurat在实现上会有细节变化。因此，如果你只是想在文章里报告倍数变化，建议在方法部分明确写清楚：差异表达分析使用Seurat的FindMarkers，倍数变化为该函数输出的avg_log2FC，而不是自己另算一套。

5.5 多重检验校正与多个比较之间的一致性

FindMarkers默认对每个比较内部做Bonferroni校正。但如果你对20种细胞类型分别处理组vs对照组都跑了一次FindMarkers，每一个比较的p_val_adj只针对该比较内的几千个基因校正，没有把20种细胞类型之间的比较次数算进去。

严格来说，如果全局控制错误发现率，需要把20个比较的所有p值汇总，再做一次FDR校正。文章里如果这个工作很重要，建议额外做一步：

all_results <- lapply(celltypes, function(ct) { ... }) combined <- bind_rows(all_results, .id = "celltype") combined$p_val_adj_global <- p.adjust(combined$p_val, method = "fdr")

这样虽然更严格，但也可能漏掉一些真实差异。是否采用，取决于你是在做探索性分析还是验证性分析。如果是探索性筛选，按每个比较内校正就够了；如果差异基因是要进后续湿实验验证，建议把全局校正的阈值也纳入参考。

5.6 包含多个样本时的批次和个体效应

如果你的Seurat对象来自多个个体，却直接用FindMarkers比较“处理组”和“对照组”，很可能把个体间差异混入进来。比如处理组里面有一个个体恰好是极端表达者，它会让大量基因上调，但这些上调跟处理本身无关。

解决思路有两条。一是用Seurat::Pseudobulk函数将同一样本内的细胞表达量汇总成“伪bulk”，得到少量但真实的生物学重复样本，然后交给DESeq2或edgeR做差异分析。二是用MAST，在公式中加入个体作为随机效应。FindMarkers本身不支持这么复杂的模型，但可以用Seurat::FindMarkers(..., test.use = "MAST", latent.vars = "sample")近似处理。

后者我实际试过，速度确实慢，但能很好地把个体异质性吸收掉。如果项目组生物学重复只有2到3个，我更推荐伪bulk方法，统计推断更稳健，也更容易写进方法部分。

6. 实际操作中的一些体会

6.1 我踩过的几个坑，希望你能绕过

第一个坑是没有先看min.pct就急着解读结果。有一次跑出来一个差异基因列表，里面有一堆神经相关的基因，乍看很有故事，后来发现它们在两组中表达的细胞都只有百分之几，完全是低表达基因的噪声。加了min.pct过滤后再看，故事完全变了。

第二个坑是只跟一个参考数据库比较，不做人工验证。FindMarkers给出的统计差异不代表这个基因一定适合做marker，也不代表它在另一个数据集中同样差异。至少要在另一个独立样本中重复一次，或者用公开数据验证一下表达模式，结果才靠谱。

第三个坑是忽略了细胞类型的混杂。比较两个cluster之前，一定先确认它们的细胞类型注释是否可靠。如果cluster 3其实是增殖性T细胞的污染，而你把它当成“效应记忆样”亚群，那么找出来的差异基因很可能只是增殖相关基因，而非该亚群的真实特征。

第四个坑是用了FindAllMarkers后，直接拿所有clusters的候选marker去画热图。热图看起来有分群梯度，但其实很多基因是从“这一群细胞相对所有其他细胞”的比较中选出来的，与另一组相比可能并没那么特异。如果要展示特定两群之间的差异，还是老老实实指定ident.1和ident.2跑一次针对性的FindMarkers。

6.2 后续可以怎么延伸

差异基因列表不会停留在表格里。拿到可靠的差异基因后，我一般会先做GO和KEGG富集，看通路层面的变化趋势；再用GSVA做通路活性打分，确认差异基因是否在相应通路中一致上调或下调；如果数据里有配体-受体信息，还可以做细胞通讯分析，看差异基因是否影响了细胞间相互作用。

另外，如果把差异基因与伪时序分析结合，还能看出这些基因是在细胞状态转换的哪个阶段开始变化的。比如某基因在naive到effector的转换中期就开始上调，这比单纯在终点比较更有动态意义。

整体而言，FindMarkers只是一个入口，它的输出是后续所有生物学解释的原材料。参数怎么调、结果怎么筛，要回到你的生物学问题本身。先确认实验设计有无伪重复，再确认细胞类型有没有混杂，最后才去折腾统计参数。流程的顺序，比你用什么方法更重要。