从标题“差异基因分析筛选出差异基因表达矩阵”来看，这应该是转录组数据分析里最常被问到、也最能决定后续分析走向的一个环节。我做生信分析这些年，带过的学生、合作过的课题组，几乎每个人都有同一个诉求：拿到测序数据之后，怎么从几万个基因里找到真正有意义的那些，然后整理成后续能直接用的表达矩阵。说白了，差异基因分析要解决的就是一个问题——在两种或者多种条件下，哪些基因的表达水平发生了不容忽视的变化。而“差异基因表达矩阵”就是这个问题的最终交付物：一个只包含显著变化基因、同时保留样本间表达量信息的二维表格。
这篇内容适合正在做转录组分析的研究生、刚转行接触生信的科研人员，也适合做了多年实验但想自己上手跑数据的实验老师。只要你的项目里出现了“对照组 vs 处理组”“肿瘤 vs 正常”这类比较，差异基因分析就绕不开。我会把完整的分析流程、参数选择的逻辑、筛选阈值的依据，以及我在实际项目中踩过的坑都写清楚，争取让你看完之后能直接开跑自己的数据。
1. 项目整体设计与思路拆解
1.1 标题背后的核心需求
“差异基因分析筛选出差异基因表达矩阵”这句话，拆开看其实是两个动作：一个是“差异基因分析”，也就是统计检验的过程；另一个是“筛选出差异基因表达矩阵”，也就是从分析结果中按阈值挑出目标基因，再回到表达量层面提取出对应的矩阵。
很多人一开始会忽略后半个动作，觉得差异分析跑完拿到一个结果表就算完事了。实际上，后续做富集分析（功能注释、通路富集）、构建预后模型、做机器学习特征筛选，甚至画热图，都需要用到筛选后的基因集合以及它们在各样本中的表达值。富集分析需要的是基因列表，热图和模型需要的是表达矩阵。这两个形态的数据，都是从差异基因分析这一步开始分叉出去的。
所以这个项目的核心需求，不只是“跑一遍统计检验”，而是构建一个规范的、可复用的分析管线的中间节点：输入是基因×样本的count矩阵，输出是两组东西——差异基因统计结果表，以及对应的标准化表达矩阵。
1.2 差异基因表达矩阵到底长什么样
先把这个最容易混淆的概念说清楚。很多初学者会把“差异基因分析结果表”和“差异基因表达矩阵”当成同一个东西，但它们格式完全不一样。
差异基因分析结果表，每一行是一个基因，每一列是统计分析指标：log2FoldChange（对数倍数变化）、p-value（显著性 p 值）、padj（多重检验校正后的 p 值），以及基因在两组中的平均表达量。这张表回答的是“基因变化有多显著、变化方向是什么”，它并没有直接给出每个样本里的表达量。
差异基因表达矩阵，则是把所有显著差异基因挑出来后，再去标准化表达矩阵里把对应的行取出来，得到一张新的矩阵：行是筛选出来的基因，列是样本，值是标准化的表达量。这张矩阵回答的是“这些基因在每一个样本里到底表达了多少”，热图、聚类、建模都用它。
我从一开始就给你一个清晰的交付概念：分析结果表用来定基因，表达矩阵用来做下游分析。两个文件分开保存，命名要规范，别混在一起用。这也是我反复跟合作方强调的“文件规范意识”，数据一多，命名混乱是最大的灾难源头。
1.3 整个分析流程的总体框架
一个完整的差异基因分析流程，从上游到下游大致是这样的：
第一步，拿到测序下机的原始数据（FASTQ）之后，做完质控、比对、定量，得到基因水平的read counts矩阵。这个矩阵通常是几万行 × 十几个样本的规模，数值是非负整数，表示每个样本中比对到每个基因上的reads数量。
第二步，对count矩阵做质量过滤。低表达基因要滤掉，样本之间要检查聚类情况，有离群样本要考虑是否剔除或展开排查。
第三步，选择差异分析工具，构建样本分组信息（condition），跑统计检验，得到每个基因的差异倍数和显著性。
第四步，设定筛选阈值。最常规的是padj < 0.05且|log2FoldChange| > 1，这个阈值组合意味着显著性和变化幅度都过关。
第五步，从标准化表达矩阵中提取差异基因对应的行，输出差异基因表达矩阵，供下游分析使用。
我特别想强调流程的“顺序感”：很多人喜欢一上来就套代码，连数据是什么样的都没搞清楚，结果跑出来的结果表一两千行，却说不清每列的含义，更说不清为什么选择DESeq2而不是edgeR。分析的本质不是跑代码，而是理解数据。我把流程放前面，代码放后面，就是希望你养成“先想明白再动手”的习惯。
2. 核心工具选型与原理剖析
2.1 三大主流差异分析工具怎么选
目前做RNA-seq差异分析的主流工具就三个：DESeq2、edgeR、limma-voom。它们的核心逻辑都是基于count数据做统计建模，但细节上各有侧重。
先看DESeq2，它在近几年的使用率非常高，尤其是在医学和生命科学领域的转录组分析中。DESeq2基于负二项分布建模，内部做了基因间方差估计的收缩（shrinkage），对小样本、低表达基因比较友好。它的另一个优势是输入格式简单，直接给count矩阵和分组信息就能跑，出错率低。
再看edgeR，同样基于负二项分布，但它最早的设计思路更偏向“数据筛选”和“快速计算”。edgeR的经典做法是先计算标准化因子（normalization factor），再用精确检验或广义线性模型做差异检验。它的计算速度比DESeq2快一点，而且在某些数据分布下检出的差异基因数量会更多一些。
limma-voom则是把原来用于芯片数据的limma框架搬到RNA-seq上来。voom先将count数据转换为log2-CPM（每百万读段数取对数），同时估计均值-方差关系，然后交给limma的线性模型流程处理。limma的优势是灵活，能处理非常复杂的实验设计，比如多因素设计、配对设计、时间序列，这在临床样本分析中特别有用。
我的建议是：新手优先学DESeq2，因为它的使用逻辑最直观，文档最完善，遇到问题最容易搜到答案。当你的实验设计比较复杂（比如含有协变量、配对样本、多分组），或者你需要的不仅仅是一对一的差异比较时，可以考虑limma-voom。edgeR则可以作为备选和对照工具，用于验证DESeq2结果的稳定性。多个工具交叉验证，是论文审稿人比较认可的做法。
2.2 差异分析背后的统计学逻辑
很多人学差异分析只记住了公式和参数，说不清背后的逻辑。我用一个比较生活化的类比来解释。
你想知道两个班级的学生成绩有没有差异，不能只看平均分，因为班级人数、学生水平波动都会影响结论。你得做统计检验，算出差异的显著程度。基因表达也是一样的道理：每个基因在两个组别里的表达量有均值差异，但每个样本之间也有波动（生物学变异）。差异分析要做的事情就是，把“组间差异”从“组内噪声”里分离出来，然后判断前者是不是显著大于后者。
count数据的分布有点特殊，它不是正态分布，而是近似负二项分布。这也是DESeq2和edgeR都基于负二项分布建模的原因。如果用普通的t检验处理count数据，会严重低估离散度，导致产生大量假阳性。
还有一个关键概念是多重检验校正。你比较的基因动辄两万多个，每一个基因做一次统计检验。如果显著性阈值为0.05，那么单凭运气就能产生约1000个假阳性（20000 × 0.05）。这显然不可接受。所以需要用Benjamini-Hochberg方法对p值做校正，得到padj（调整后p值）。通过padj < 0.05这个标准筛选出来的基因，才更有可能是真正的差异基因。
我在实际分析时，筛选基因几乎总是以padj为准，而不是原始p值。原因就在这里：不校正的结果经不起重复实验的检验。
3. 实操过程与代码实现
3.1 输入数据准备：从count矩阵到DESeqDataSet
先用一个实际项目来演示。假设我们有6个样本，3个是对照组（control），3个是处理组（treat），测序定量之后得到原始count矩阵。
count矩阵的标准格式是：第一列是基因ID（Ensembl ID或Symbol），后面每一列是一个样本，每个单元格是对应的read count数。读进R之后，要把基因ID设为行名。
# 安装并加载DESeq2（需要先安装BiocManager） if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("DESeq2") library(DESeq2) # 读取count数据，注意check.names = FALSE防止R自动修改样本名 count_data <- read.csv("count_matrix.csv", row.names = 1, check.names = FALSE) # 查看矩阵维度：应该是 基因数 × 样本数 dim(count_data)count_data就是整个分析的输入起点。这里有个细节：如果上游定量工具是featureCounts或HTSeq，输出的第一列通常是Ensembl Gene ID；如果你需要把它转换成gene symbol，可以在后续做注释，也可以在筛选结果之后再做。我更倾向于在差异分析之前就完成ID标准化，否则下游富集分析时还得反复映射。
接下来构建样本分组信息。这一步最关键，也是最容易出错的一步。colData需要是一个数据框，行名必须和count矩阵的列名完全一致，顺序也要一致。
# 构建分组信息，行名必须等于count矩阵的列名 col_data <- data.frame( condition = factor(c("control", "control", "control", "treat", "treat", "treat"), levels = c("control", "treat")), row.names = colnames(count_data) ) # 检查行名是否匹配 all(rownames(col_data) == colnames(count_data)) # 构建DESeqDataSet对象 dds <- DESeqDataSetFromMatrix(countData = count_data, colData = col_data, design = ~ condition)这里design = ~ condition的意思是，我们要比较的核心因素是condition，也就是对照组和处理组之间的差异。如果你有批次信息，后续可以在design里加入，比如design = ~ batch + condition。
3.2 低表达基因过滤：为什么不是可选项
在跑DESeq之前，建议先做一步低表达基因过滤。原因很简单：一个在所有样本里表达量都接近于0的基因，它的count值几乎全是0，没法产生有意义的统计检验结果。而且保留过多低表达基因会增大多重检验校正的负担，降低统计功效。
常用的过滤策略有两种。第一种是基于总counts数，比如保留在所有样本中counts总和大于等于10的基因。第二种是基于表达量阈值，比如保留至少在n个样本中counts > 10的基因，这里的n可以取最小样本组的大小。
# 过滤：保留在至少3个样本中counts大于等于10的基因 keep <- rowSums(counts(dds) >= 10) >= 3 dds <- dds[keep, ] # 过滤后查看剩余基因数 nrow(dds)我在实际操作中，倾向于用第二种方式，因为它考虑了每个基因在多少样本里有表达，比单纯看总和更合理。比如一个基因只在某一个样本里异常高表达，总和也能超过10，但它在其他样本里都是0，这种基因很难稳定支持结论。
3.3 DESeq2差异分析核心步骤
过滤之后，直接跑DESeq函数，它会依次执行标准化、离散度估计和差异检验，因此这个函数通常需要跑一段时间，样本量大的时候可能要几分钟到十几分钟。
# 运行差异分析 dds <- DESeq(dds) # 提取差异结果，contrast指定比较的两组 res <- results(dds, contrast = c("condition", "treat", "control")) # 转为数据框，方便查看和筛选 res_df <- as.data.frame(res) res_df$gene <- rownames(res_df)res_df里每一行是一个基因，主要列的含义是：
- baseMean：该基因在所有样本中的标准化平均表达量，相当于背景表达丰度。
- log2FoldChange：处理组相对于对照组的表达倍数变化的对数，正值表示上调，负值表示下调。
- lfcSE：log2FoldChange的标准误，反映了该估计的稳定性。
- pvalue：未校正的统计检验p值。
- padj：多重检验校正后的p值，筛选时的核心依据。
这里有一个细节：contrast = c("condition", "treat", "control")的写法很关键，它明确告诉DESeq2比较的分子和分母分别是谁。分子（treat）在log2FoldChange的正值方向上，分母（control）在负值方向上。我之前见过有人图省事直接写results(dds)，这样用的其实是design里最后一个水平作为分母，如果你没注意因子水平顺序，很容易把上下调搞反。这个坑虽小，但会直接影响后续所有解读。
3.4 筛选阈值怎么定才合理
筛选差异基因，行业惯例是padj < 0.05且|log2FoldChange| > 1。这个阈值组合的具体含义是：经过多重检验校正后依然显著的基因，同时表达变化幅度至少要翻倍或减半（log2FoldChange = 1对应表达量变为原来的2倍，-1对应变为原来的0.5倍）。
有时候样本间差异较小、或生物学重复间波动较大，严格的阈值筛出来的基因数量会很少，可能只有几十个。这时候可以适当放宽，比如padj < 0.05保持，但log2FoldChange阈值降到0.5，或者用pvalue < 0.01来替代padj。不过放宽阈值意味着假阳性风险上升，在后验证阶段要格外谨慎。
# 标准筛选：padj < 0.05 且 |log2FoldChange| > 1 sig_genes <- subset(res_df, padj < 0.05 & abs(log2FoldChange) > 1) # 统计上调/下调基因数量 table(sig_genes$log2FoldChange > 0) # 按log2FoldChange排序 sig_genes <- sig_genes[order(sig_genes$log2FoldChange, decreasing = TRUE), ]这里我还建议做一个小动作：把结果表按照log2FoldChange从大到小排序，然后人工检查排名靠前的20个基因，看看它们是否符合你的生物学预期。这比单纯看数量更有价值，因为如果排名最靠前的基因跟你的实验背景完全无关，那大概率是数据质量或分组设置出了问题。
3.5 输出差异基因表达矩阵：关键一步
筛选出差异基因之后，回到exp矩阵（标准化表达矩阵），提取差异基因对应的行，这就是最后交付的差异基因表达矩阵。
# 提取标准化表达矩阵 norm_counts <- counts(dds, normalized = TRUE) # 提取差异基因对应的表达矩阵 sig_matrix <- norm_counts[rownames(norm_counts) %in% sig_genes$gene, ] # 为了后续查看方便，把gene列作为第一列 sig_matrix <- as.data.frame(sig_matrix) sig_matrix$gene <- rownames(sig_matrix) sig_matrix <- sig_matrix[, c("gene", rownames(norm_counts))] # 保存结果表 write.csv(sig_genes, "DEGs_results.csv", quote = FALSE) # 保存差异基因表达矩阵 write.csv(sig_matrix, "DEGs_expression_matrix.csv", quote = FALSE)写完保存之后，我建议立刻打开这两个文件，抽查几个基因手动验证一下。比如A基因在treat组里counts均值是对照组的4倍，理论上它的log2FoldChange应该接近2。这个验证动作能帮你确认整个流程没有出逻辑错误。
另外提醒一句：输出的差异基因表达矩阵，值默认是DESeq2的标准化counts，它适合做热图展示，也适合大多数下游分析。但如果你想用它作为机器学习模型的输入特征，有人会先做log2转换或者z-score标准化，这取决于你后续的具体任务。不要默认一步到位，要清楚每一个操作的目的。
4. 常见问题与排查技巧实录
4.1 分组设计错了怎么办
这是我在实际项目中见过最多的问题，尤其发生在交付数据时样本顺序被打乱的情况。count矩阵的列名和分组信息对应不上，分析结果就会完全失真。我见过最惨的情况是，一个学生把对照组和处理组标签弄反了，跑出来一组显著差异基因，拿去做富集分析，得的通路全是炎症和免疫——后来才发现是因为标签反了，基因方向反了，通路的解读也反了。
排查技巧其实很简单：跑完DESeq之后，先不要急着看差异结果，先看样本聚类。用rlog或vst变换后的表达量做主成分分析（PCA），看看样本是否按分组聚在一起。
# 使用rlog变换，PCA可视化样本聚类 rld <- rlog(dds, blind = TRUE) plotPCA(rld, intgroup = "condition")如果你的处理组和对照组在PCA图上完全分不开，要么是生物学差异太小，要么就是样本分组信息出了问题。如果同组样本分散在不同象限，还要考虑批次效应或者样本污染的可能。
4.2 批次效应需要处理吗
批次效应是转录组分析中的老生常谈。不同批次提取的RNA、不同时间跑的测序、甚至不同机器运行的结果都会有系统性偏差。如果实验设计阶段没有平衡分组（比如所有对照组第一批跑，所有处理组第二批跑），批次效应就会和真正的生物学差异混淆在一起，怎么分析都会出问题。
处理方案有两种。第一种是在design中加入批次变量：design = ~ batch + condition。这样DESeq2在建模时会先把批次效应剔除掉，再评估condition的影响。第二种是对校正后的表达矩阵用ComBat_seq或removeBatchEffect等工具做批量校正。前者适用于count级别的正式分析，后者常用于下游可视化。
我的建议很直接：尽量在实验设计阶段避免批次混杂。已经发生混杂的情况下，能加批次变量就加批次变量。如果批次信息都没有记录，那就只能祈祷影响不大，但这种情况非常危险。
4.3 低表达基因到底要不要过滤
有些初学者不习惯做过滤，因为怕“丢失信息”。但过滤这一步是必要的。低表达基因的count值波动几乎全是噪声，保留它们只会增加多重检验校正的负担，降低真正差异基因的检出效率。你要知道，padj校正的目的是控制全局假阳性率，当低表达基因大量存在时，padj会更严格，很多真实差异也会被“压”下去。
不过过滤强度要适中。过度过滤会误伤部分真实的但表达量较低的基因，尤其是转录因子这类本身表达丰度就不高的调控基因。我的经验是：先用下限比较宽松的标准（比如至少在n个样本中counts >= 5）跑一遍，再用严格一点的标准（counts >= 10）跑一遍，看结果中关键基因是否稳定保留。关键基因稳定，说明结论可信。
4.4 为什么padj和pvalue差这么多
经常有人拿着结果表来问我：老师，我这个基因pvalue是0.0001，怎么padj变成0.2了？是不是软件出bug了？
不是bug。这是多重检验校正起作用的正常现象。假设你检测了两万个基因，有一个基因的原始p值是0.0001，看起来非常显著。但你要想想，在两万次检验中，哪怕没有真实差异，也有可能出现这么小的p值。Benjamini-Hochberg方法会把所有基因的p值排个序，根据你在比较中的总检验次数和排序位置重新计算显著性。所以，当一个基因的原始p值很小但它排在较后的位置（因为前面还有成百上千个p值更小的基因）时，padj就会变得很大。
反过来也能解释另一个现象：有时候两组间差异极小的基因（log2FoldChange接近0），因为p值足够小而被筛选出来，这类基因在生物学上很难被解读。这也是为什么筛选时一定要同时看padj和log2FoldChange，两者缺一不可。
4.5 差异基因数量太少或者太多怎么处理
筛选出来的差异基因数量，直接反映你的实验设计和数据质量。数量太少的情况下，大多是阈值太严、样本重复太少、生物学差异较弱。我的经验是先检查padj分布：把所有基因按padj排序，看pvalue接近0的基因有多少。如果寥寥无几，说明数据本身没有足够的分辨能力。
数量太多（比如上调下调各上万）通常也值得警惕，这往往意味着组间存在非常强的系统性差异，比如批次效应没有被剔除，或者两组样本本质上就不是同一类样本。遇到这种情况，先做PCA排查样本是否来源于同一批组织。
一个实用的中间检查指标是baseMean。差异基因的信息量与表达丰度高度相关，大多数可靠的差异基因baseMean都在几十到几千之间。如果筛出来的显著基因大量集中在baseMean < 10的范围，数据质量就有问题。这个经验其实对应了一个很朴素的道理：表达量太低的基因，测量的稳定性很差，很难被重复实验验证。
5. 实操心得与经验沉淀
先说一个我自己反复踩坑之后形成的原则：任何分析跑完，第一件事永远是手动验证，而不是画图。我会随机挑两三个显著上调的基因、两三个显著下调的基因，回到原始count矩阵手动算一下每组的平均表达量，确认方向正确、幅度合理。这一步花不到五分钟，却能挡住绝大多数低级错误。
另一个技巧是管好自己的中间文件。我习惯在一个项目目录里建立counts、results、plots、tables四个子目录，分别存放原始输入、差异分析结果、可视化图片和交付矩阵。文件名里写清楚版本和日期，比如DEGs_results_v2_20240615.csv。以前没这个习惯的时候，一个月后自己翻文件夹都分不清哪个是最终版本，更别提跟别人交接了。
关于软件版本，我也啰嗦一句：DESeq2这个包每年都在更新，不同版本的默认参数可能有细微差异。发表在论文里的分析结果，要记录好R版本和DESeq2版本，并且尽可能用固定环境跑完整套分析。严谨一点的课题组会直接用Docker或conda锁定环境，这样不管什么时候复现，结果都能保持一致。
最后再分享一个小建议：差异基因表达矩阵输出来之后，别急着丢进自动化的下游流程。先做一个最基础的热图，在热图里看看显著基因是否真的能按分组区分样本。如果热图的聚类分支和你的分组不一致，说明前面某个环节出了问题，这时候停下来排查，比带着错误结果往下走要省时间得多。差异基因分析是整个转录组项目的地基，这个地基打歪了，后面建的楼越高越危险。
