1. 差异基因分析前传：从原始测序数据到表达矩阵的必经之路

很多刚接触转录组的人会有一个误解：拿到测序公司的报告，里面有现成的基因表达量表格，直接拿来做差异分析就行了。实际上，那份报告的生成过程比你想象的要复杂得多，而"筛选出差异基因表达矩阵"这件事，误差恰恰容易在最不起眼的环节悄悄埋下。

先说整条链路：下机得到的FASTQ原始数据，要先做质量评估，然后比对到参考基因组或转录组，再按基因或者转录本级别统计read数，最后才得到我们熟悉的表达矩阵（count matrix）。这个矩阵的行是基因，列是样本，每个格子里面是一个非负整数，代表这个基因在这个样本里检测到的reads数量。后续所有差异分析、富集分析、可视化，全部建立在这个矩阵之上。换句话说，表达矩阵是原料，原料出了问题，后面做出来的菜再好看也是废的。

所以第一个关键决策点就是：用什么工具比对和定量。目前主流的方案有两种：一种是STAR比对到参考基因组，再用featureCounts统计基因水平的read数；另一种是直接用Salmon或者kallisto这类伪比对工具，跳到转录组层面定量，拿到转录本水平的TPM或者read count后再用tximport汇总到基因水平。

我自己用得最多的是STAR + featureCounts这条路线，原因是它对后续差异分析的兼容性最好，尤其是做基因水平分析的时候。STAR的映射速度快，内存占用虽然在比对工具里算偏高的（一个人类基因组索引大概需要30GB左右内存），但换来的是极高的比对准确性。featureCounts可以非常方便地从BAM文件里统计出每个基因的read计数，同时支持链特异性文库的设定，这一步很容易被忽略。

这里就踩过一个大坑：链特异性参数设错了也不行。如果文库是链特异性的（比如dUTP法建库，这是目前最主流的建库方式之一），但你用featureCounts统计的时候没有加-s 2（表示read1和基因反义链对应），那统计出来的计数会有相当一部分比对到错误的链上，表现为反义链的read数显著高于正常值，最终导致基因为一个基因的正反链计数几乎一样，后续差异分析的结果自然会出现大面积失真。

还有一个容易忽略的参数是多比对reads的处理方式。featureCounts默认只统计quality score足够高的唯一比对reads，这没什么问题。但如果你的物种是多倍体或者基因组里有高度相似的重复区域，那大量reads会被丢弃，这时候就需要评估一下unique mapped reads的比例。我习惯在建完索引之后，先跑一个样本做快速检查，看一眼STAR日志里的Uniquely mapped reads percentage，如果低于75%，就要考虑是否是参考基因组版本不合适，或者是转录组里重复序列比例太高的问题。

表达矩阵拿到了，接下来还有一个逃不掉的问题：要不要做基因注释过滤。参考基因组的GTF注释文件里，有大量的非编码RNA、假基因、线粒体基因，这些基因在标准的mRNA-seq差异分析里通常是要根据研究目的过滤掉的。比如线粒体基因占比过高，往往代表细胞状态有问题；rRNA残留会导致某些核糖体蛋白基因的计数异常偏高。我会在拿到count矩阵之后，根据GTF的gene_biotype字段，只保留protein_coding和部分lncRNA，具体保留哪些取决于文章的研究对象——如果课题本身就是研究lncRNA的，那就另当别论。

另外，表达矩阵的行名问题看似不起眼，却经常在后续分析里造成不可挽回的错误。不同来源的GTF文件，基因ID有可能是Ensembl ID（如ENSG00000141510）、也有可能是Symbol（如TP53），还有可能是NCBI的Gene ID。一旦某个环节混用了ID类型，后续的注释、富集分析都会对不上号。我的做法是从头到尾统一用Ensembl ID作为唯一标识，所有转换都集中在最后做注释时才进行，这样即使SQL脚本、R脚本来回切换也不容易出错。

走到这一步，原料基本就准备好了。但真正的"差异基因分析"是从这里才开始的。

2. 差异分析工具选型：DESeq2、edgeR、limma-voom到底怎么选

很多人总爱问"哪个差异分析工具最好"，其实这个问题本身就不严谨。转录组数据是一种典型的计数数据，分布高度偏离正态，而且不同基因之间的表达量基数差异极大——有的基因平均表达量上万，有的只有几十甚至几个read。直接拿普通的t检验去处理这样的数据，假阳性会高得离谱。所以差异分析工具的核心任务，是把这种计数数据的分布特征合理地建模，从而在有限样本量下给出尽量可靠的统计推断。

目前最主流的三个工具分别是DESeq2、edgeR和limma-voom。它们在统计思路上有本质的区别。

DESeq2用的是负二项分布模型（Negative Binomial），它假设每个基因的read count符合均值与方差相关的负二项分布，然后用shrinkage（收缩估计）的方式估计基因的离散度参数。这个设计的好处是当你样本量较小、基因表达极差大的时候，它会借助所有基因的离散度信息来"借用"信息，把方差估计拉向一个合理的整体水平，从而减少假阳性。坏处是它对每个基因都要拟合两层广义线性模型，计算速度相对较慢，但目前数据量下通常几分钟到十几分钟就能跑完，不算瓶颈。

edgeR的底层模型也是负二项分布，但它的标准化方式用的是TMM（Trimmed Mean of M-values），整体思路比DESeq2更"轻量"。它更适合做样本数量较多的差异分析，也适合基因数量较大时的快速筛选。如果你的项目里样本数特别多（超过20个），或者你只想做一个快速的初步筛选，edgeR是个好选择。

limma-voom走的是另一条路。limma本身是为芯片数据设计的线性模型工具，voom算法则把RNA-seq的counts转换成log2-CPM值，同时估计每个观测值的均值-方差关系，并把这个精度权重融入到线性模型里。这套组合拳的优势是极其稳定，尤其是在做复杂实验设计（多个因素交叉）的时候，limma的线性模型框架处理起来比DESeq2和edgeR都更自然。此外limma速度极快，几千个样本跑起来也不吃力。

我自己的选型经验是这样的：

这些工具有一点是共通的，它们的输出里都包含log2FoldChange（对数倍数变化）、pvalue和padj（校正后的p值），而"差异基因表达矩阵"正是从这些结果里筛选出来的。

不过在真正运行之前，还有一个绕不开的预处理步骤：标准化。原始count矩阵的每个样本测序深度不一样，总reads数可能差两三倍，直接比较的话，测序深度大的样本看起来所有基因都"上调"了，这显然是错误结论。DESeq2内部通过median-of-ratios方法计算了size factor来自动处理测序深度差异，edgeR用TMM，limma-voom则用CPM标准化。所以在跑差异分析的时候，不需要也不应该在跑之前手动把count值转换成TPM——因为这些工具本身就对count数据做了标准化处理，你手动转了一道反而会破坏负二项模型的计数分布假设。

选择工具之后，还有一个细节值得单独说说：因子水平的设定。R语言里因子水平的顺序会影响结果的解读方向。比如你的处理组叫"Treat"，对照组叫"Ctrl"，如果R默认把Ctrl作为第一个水平，那么DESeq2默认输出的log2FoldChange就是"Treat相对Ctrl"的变化；但如果因子的水平顺序反过来，log2FoldChange的符号也会反过来，火山图里上调下调的方向就完全颠倒。这个问题我在给人debug时遇到过好几次，明明代码每一步都对，结果一查发现是对照组的水平排到后面去了。

把工具确定下来，接下来才真正进入到筛选差异基因的核心环节。

3. 筛选差异基因的核心流程：阈值设定、结果解读、可视化三板斧

我见过不少人在拿到差异分析结果表之后，随手设一个padj < 0.05、|log2FC| > 1就完事了。这个流程本身没错，但如果你不理解这两个阈值背后的含义，后续就会被各种"伪差异基因"折磨。

先说padj这个概念是怎么来的。对一万多个基因分别做假设检验，每个基因都会得到一个p值。p值的含义是"在原假设为真（基因无差异）的情况下，观察到当前这么大的差异甚至更大差异的概率"。假设你设置的显著性水平是0.05，那么即使所有基因都没差异，每一万个基因里也大约有500个会因为你运气不好而出现p<0.05的结果。这就是多重假设检验问题。

为了解决这个问题，Benjamini和Hochberg提出了一个叫BH校正的方法，通过控制FDR（False Discovery Rate，错误发现率）把所有p值按大小排序后逐步调整，得到padj。padj的含义不再是"单个基因出错概率"，而是"在所有你判定为显著的基因里，预期错误的比例"。比如你筛出了1000个padj<0.05的基因，预期其中50个其实是假阳性。这是一种全局的风险控制，比单个p值更可信。

|log2FC| > 1的意思则是倍数变化（fold change）不低于2倍。因为log2FC = 1代表表达量相差2倍，log2FC = 2代表相差4倍。设这个阈值是为了过滤掉那些统计学显著但生物学效应微小的基因。如果一个基因的padj只有0.001，但log2FC只有0.1，说明它只是测序深度高、统计功效强，真实效应可能没有实际意义。

但是请注意，这两个阈值并不是铁律。在某些场景下需要更严格：比如你做全基因组筛选，研究的又是一种表达变化幅度普遍很小的调控因子，那么用|log2FC| > 1可能一个基因都筛不出来。这时候可以放宽到|log2FC| > 0.5，甚至结合具体情况单独看padj。反过来，如果你的研究是临床标志物筛选，要求结果可重复性极高，那可以把padj设成0.01甚至更严格。筛选阈值应该服务于你的生物学问题，而不是反过来。

跑完差异分析之后，R里面通常是这样操作：

library(DESeq2) # 读取count矩阵 countData <- read.table("gene_counts.txt", header = TRUE, row.names = 1) # 样本信息，保证列名与count矩阵的列名对应 colData <- data.frame( condition = factor(c("Ctrl","Ctrl","Ctrl","Treat","Treat","Treat"), levels = c("Ctrl","Treat")) ) dds <- DESeqDataSetFromMatrix( countData = countData, colData = colData, design = ~ condition ) # 过滤：保留至少在3个样本里总read数大于10的基因 keep <- rowSums(counts(dds) >= 10) >= 3 dds <- dds[keep, ] dds <- DESeq(dds) res <- results(dds, contrast = c("condition", "Treat", "Ctrl")) res_df <- as.data.frame(res) # 添加基因ID列 res_df$gene_id <- rownames(res_df) # 筛选差异基因 sig_genes <- subset(res_df, padj < 0.05 & abs(log2FoldChange) > 1)

筛选完之后，最重要的一步是检查结果是否合理。我通常会做三个检查：

第一，差异基因的数量和上下调比例。如果上调基因数量是下调基因的十倍，先不要激动，排除一下是不是某组样本质量有问题，比如处理组细胞污染、样本标签搞反等等。

第二，已知的阳性对照基因是否出现在名单里。比如你做某种药物处理，如果这个药物公认会上调某个关键通路基因，而这个基因没在差异基因列表里，那大概率是实验环节出问题了，这时候无论统计结果多么漂亮都要回去查数据来源。

第三，主成分分析（PCA）的可视化结果。在DESeq2里用vst或者rlog变换后的数据做PCA，如果不同组的样本能明显分开，说明组间差异显著大于组内差异，差异分析的结果就相对可信；如果PCA图上两组样本混在一起，那即使跑出一堆"差异基因"，也很可能是批次效应而非真实生物学差异。

拿到差异基因列表之后，接下里就是"差异基因表达矩阵"的真正落地环节了。很多人会忽略这一环：纯粹输出一份筛选结果表还不算完，还要把这组基因在样本里的表达情况整理成方便可视化和下游分析的矩阵格式。这一步的操作逻辑是这样的：从原始的vst标准化矩阵（或者TPM矩阵）中，用差异基因的行名去子集，得到一个新的小矩阵，行是显著差异基因，列是样本。这个矩阵才是后续画热图、做聚类分析、跑富集分析的输入文件。

vsd <- vst(dds, blind = FALSE) vsd_mat <- assay(vsd) # 用差异基因子集表达矩阵 deg_exp <- vsd_mat[rownames(sig_genes), ] # 将矩阵排序并存储 deg_exp_sorted <- deg_exp[order(sig_genes$log2FoldChange, decreasing = TRUE), ] write.csv(deg_exp_sorted, "DEG_expression_matrix.csv")

这个deg_exp_sorted就是我们要的东西。它不但是热图绘制的直接输入，也是后续GSEA、GO/KEGG富集分析、蛋白互作网络构建等一系列下游分析的起点。

4. 常见坑点与排查思路：从批次效应到多重检验的完整链路

我在前面提到PCA可以用来检查样本分组，但PCA还有一个更重要的作用：发现批次效应。批次效应是转录组项目里最常见、也最致命的问题之一。简单来说，如果你的3个对照组样本在一个批次跑完、3个处理组样本在另一批次跑完，那么样本间的差异可能更多来源于"不同批次的技术误差"，而不是真实的生物学差异。

排查批次效应最直观的方法是看PCA图的前两个主成分。如果PC1（主成分1）能完美分开对照组和处理组，同时PC2或PC3上还出现了一批按"上机日期"或者"送样时间"聚类的样本，那就要高度怀疑批次效应的存在。还有一种情况更隐蔽：两组样本本身生物学差异很大，比如野生型和基因敲除型，这时候第一批次的整体效应会叠加在真实差异之上，PCA图看不出来，但不代表批次效应不存在。

针对批次效应的处理，有几种常规手段。最理想的情况是实验设计时就把不同条件随机打散到不同批次，俗称"随机化"，但这在现实中经常做不到。如果已经做了，能补救的标准方案是用limma包的removeBatchEffect函数，传入已知的批次变量来校正数据；或者用RUVSeq、sva这类工具先估计隐藏的因子，再在差异分析模型里把这个因子作为协变量放进去。

我实际跑过的最典型的坑是在DESeq2的设计公式里漏掉了批次变量。比如设计公式只写了~ condition，但样本实际上分了几个批次，结果差异分析把批次效应和生物学效应混在一起，导致大量假阳性基因——尤其是那些对批次敏感的低表达基因。后来把设计公式改成~ batch + condition，差异基因数量从7000多降到了1000多，富集结果也跟文献报道更吻合了。这件事给我的教训是：任何与生物学处理无关的已知技术变量，都应该尽可能纳入设计公式。

再说一个容易被忽视的低表达基因过滤问题。很多人一上来就直接跑DESeq2，不过滤任何基因。这时候大量在所有样本里读数都是0或者只有1~2的基因也会参与离散度估计和多重检验校正。这些基因不仅会增加多重检验矫正的负担（padj会变得更保守），还会在计算离散度时引入噪音，导致真实表达、稳定表达的基因反而被误伤了。

那到底怎么过滤？通常的做法是保留至少在最小样本数（比如每组最少样本数）里，count数达到一定阈值（比如10或者1）的基因。注意，这个过滤用的count数，一定要用原始的整数count，而不是TPM或者CPM。原因是在模型假设里count服从负二项分布，低计数的基因信息量本来就少，过滤掉它们不至于损失太多真实的生物学信号。如果用TPM过滤，大量本身转录本长度很短的基因会被错误地丢掉——因为TPM对转录本长度做了归一化，短基因的TPM天然偏低，但这不代表它不表达。

接下来是多重检验校正顺序的问题。R里results()出来的padj是已经用BH方法校正过的，如果你自作聪明地对pvalue再调一次p.adjust，会把结果弄得更严格或者更宽松——这取决于你用的方法。我见过有同学在拿到res之后发现padj列全是NA，然后自己重新用p.adjust(res$pvalue, method = "bonferroni")跑了一遍，结果把差异基因筛得一个不剩，最后跑来问我为什么结果为空。一问才知道他把已经被独立筛选过一次的p值又做了一次非常保守的Bonferroni校正。这不是工具问题，而是对统计流程的理解问题。

还有一组需要单独拎出来说的基因：那些padj为NA的表达量极低基因。在DESeq2里，如果一个基因在所有样本里的count都很低，它可能没有足够的数据来估算离散度，对应的padj就会显示为NA。在筛差异基因时，千万不要把含NA的阈值判断直接拿去比较——R里NA > 1的结果是NA，这会让你的子集结果莫名其妙少掉大量的行。安全的做法是先过滤掉padj是NA的行，再做筛选。

低样本量下的p值分布偏移问题也值得一提。当每组只有2~3个样本时，单基因的p值分布往往不服从均匀分布，而是偏向0或偏向1——这是因为极度缺乏自由度导致的检验功效不足。这种情况下DSSeq2的收缩估计虽然能在一定程度上缓解，但仍要警惕：如果你在每组只有2个样本、甚至没有生物学重复的情况下硬跑差异分析，得到的p值分布会非常糟糕，很多p值极小的基因其实只是在两组样本间恰好差异大而已，没有可重复性可言。这也是为什么很多期刊现在强制要求至少3个生物学重复作为RNA-seq实验的最低门槛。

最后要说的排查手段，是把自己"钉"在结果的逻辑合理性上。跑完差异分析后，脱离开统计术语，回到你的实验本身去问三个问题：这个处理在生物学上到底应该影响哪些通路？这些通路里的代表性基因在结果里出现了没有？结果里有没有哪些基因跟已知机制完全相悖？这些问题的答案，往往比任何统计指标更能帮你发现深层的数据问题。有一次我跑完一个药物浓度梯度的转录组数据，差异基因列表里竟然出现了大量的糖酵解相关基因上调，而原课题做的是免疫调控，怎么也想不通。后来重新检查才发现，处理组的细胞培养时间比对照组长了整整一天，细胞密度过高导致培养基缺氧，糖酵解通路是被环境逼出来的。这个例子说明，差异基因分析筛出来的东西，很多时候反映的是你不知道的实验偏差，这时候回到实验记录找原因，而不是急着写结论。

5. 差异基因表达矩阵的落地应用：热图、富集、实验验证一条龙

差异基因表达矩阵做出来之后，实际工作中最常见的下游使用方式有三种：可视化（热图和火山图）、通路富集分析、以及实验验证前的候选基因清单整理。每个环节的侧重点各不相同，但都可以围绕这个矩阵展开。

热图是展示差异基因表达模式最强有力的工具。绘制的时候选基因数量如果太多（比如超过1000个），整张热图会变成一团密集的马赛克，看不出任何模式。这时候我通常会过滤掉那些在各个样本间表达量波动特别小的基因，或者先按方差排序取top50~100画图。重要的一点是：热图里的颜色映射要基于标准化后的表达量（vst或者z-score），不能直接拿原始count画，因为不同样本的测序深度差异会让所有基因的颜色都偏暗或偏亮，掩盖真实的差异模式。

library(pheatmap) # 对每一行（基因）做z-score标准化 heatmap_mat <- t(scale(t(deg_exp_sorted))) pheatmap( heatmap_mat, cluster_cols = TRUE, cluster_rows = TRUE, show_rownames = FALSE, annotation_col = data.frame(condition = colData$condition, row.names = colnames(heatmap_mat)) )

热图只是"看"差异，功能富集分析才是帮助你把差异基因列表转化成一串能被审稿人理解的生物学结论的关键步骤。常规做法是使用clusterProfiler做GO和KEGG通路富集分析。

GO富集分析的输入通常是一组差异基因ID加上背景基因ID（也就是你测到的所有基因，而不是全基因组所有基因）。这一点非常关键：如果你不指定universe（背景集），富集分析会拿全部基因作为背景计算富集显著性，测序技术天生只能检测到一部分基因，这样算出来的富集p值会被严重高估，出现大量本来不显著的"伪富集"结果。

library(clusterProfiler) library(org.Hs.eg.db) # 转换ID sig_entrez <- bitr(sig_genes$gene_id, fromType = "ENSEMBL", toType = "ENTREZID", OrgDb = org.Hs.eg.db) # GO富集分析 ego <- enrichGO( gene = sig_entrez$ENTREZID, universe = background_entrez$ENTREZID, OrgDb = org.Hs.eg.db, keyType = "ENTREZID", ont = "BP", pAdjustMethod = "BH", qvalueCutoff = 0.05 )

GO富集分析虽然通用，但有个老毛病：术语体系太庞大，结果动辄几十上百条显著性通路，读起来像天书。这时候我的经验是先用simpleEnrichment或者手动的冗余去除（reduce）把相似的GO term合并，再做一层GSEA分析来辅助判断哪些通路里的基因是一致上调或一致下调的。GSEA和普通的GO富集不同，它不要求你先设定一个差异基因列表，而是直接拿所有基因的表达变化排序和foldchange信息计算通路级别的得分。这样能抓住那些单个基因变化不显著、但整个通路协同上调或者下调的生物学信号，非常适合处理那种效应"积少成多"的通路。

矩阵的下一个落地场合是实验验证，这个环节往往是整个分析链条里最容易被低估的。差异分析出来一堆基因，热图富集都很漂亮，但真要落到下一步wet lab实验，你不能把1000个基因都拿去做qPCR或者Western验证。这时候矩阵的价值就体现出来了：你要从这个矩阵里挑出那些表达变化幅度大、组内一致性高、且在已有文献里有明确机制支撑的代表性基因，而不是单纯挑log2FC最大的那几个。

什么叫组内一致性高？举个例子：某个基因在三个处理样本的log2FC分别是2.0、2.1和0.9，平均值1.67，看起来显著上调，但实际上第三个样本的数据明显偏低，可能是样本处理不均匀或者操作误差导致的。如果这个基因恰好属于一个跟课题方向高度相关的通路，你拿它去做qPCR验证，十有八九会验证失败。反过来，如果某个基因log2FC只有1.2但三个样本之间几乎无波动，它在qPCR里重复出来的概率反而更高。在矩阵里检查组内一致性，我通常会看每个基因在不同样本里的原始标准化表达值绘制箱线图，或者直接计算每组内变异系数。

还有一个实用技巧：把差异表达矩阵直接作为验证实验设计的参考。你在设计qPCR验证引物或者挑选做免疫组化的抗体时，先确认备选基因在你的矩阵里确实有足够的表达量（比如vst标准化后的均值大于某个阈值），而不是只盯着log2FC大小。表达量太低的基因在qPCR里CT值会超过35循环，验证结果会因为技术噪音而变得完全不可信。

说到这，其实"差异基因表达矩阵"的下游应用是很灵活的。有位同事做过一个有意思的尝试：他把筛选出的差异基因表达矩阵和转录因子结合位点数据库结合，通过共表达网络分析找到了一批与核心转录因子共调控的模块基因，然后把这些模块基因和临床患者的预后数据做了生存分析，最后硬是从一个普通的细胞系转录组项目里挖出了三个新的预后标志物候选基因。这个例子很好地说明，差异基因表达矩阵不仅仅是分析链路的终点，它更像一个蓄水池，后续只要你愿意，随时可以从池子里引出不同的分析支流。

6. 实操总结：我踩过的坑，希望你不用再踩一遍

最后聊几句个人经验层面的东西，不展开写大理论，只说那些我最初不理解、后来被现实教育过才明白的细节。

第一，永远不要直接拿公司给的"基因表达量"表去做差异分析。很多公司给的表格里面对应的可能是FPKM、TPM或者已经做了某种标准化的值。FPKM/TPM作为跨样本比较的表达量没问题，但它们不应该被输入到DESeq2或edgeR里做差异检验。这两个工具天然假设输入是整数count，你把小数输入进去，要么报错，要么结果畸变。如果不幸你手上只有FPKM矩阵，可以用tximport配合转录本定量结果重建count数据，或者换用limma-voom流程来处理（voom内部会对log-counts做拟合，对这类数据的容忍度更高一些）。

第二，log2FoldChange和Log2CPM是两个完全不同的东西，不要混用。差异结果表里的log2FoldChange是指两个组之间标准化后表达量的对数比率，正负号代表上下调方向；而Log2CPM是各个样本里基因的表达量归一化值，只用于展示和比较表达水平。曾经有人问我为什么某个基因log2FoldChange = 2.5但热图上颜色不深，仔细一看热图用的是原始count矩阵，样本内其他基因的表达量太高，把颜色分布整个压扁了。这个问题说大不大，说小不小，但在审稿人眼里这就是数据处理的硬伤。

第三，运行DESeq2之前，仔细检查样本列名是否和colData的行名完全一致。R里面DESeqDataSetFromMatrix要求列名与colData的DataFrame行名精确匹配，哪怕多一个空格、一个下划线，函数都会报错。这个错误信息有时候并不直观，最常见的报错是"colData must have the same samples as colnames of countData"。如果你看到这行报错，先检查一下是不是样本名没对上，而不是急着改代码。

第四，不要忽略基因注释版本的核对。同一个物种，不同版本的GTF文件，基因ID可能发生改变（gene symbol改了，Ensembl ID也可能有版本后缀，比如ENSG00000141510.17和ENSG00000141510.10是不同的东西）。在多个GTF版本之间混用时，差异基因筛选结果的可靠性会大打折扣。最稳妥的操作是在项目开始时固定一个参考基因组版本和一个GTF版本，后续所有比对、计数、注释都基于它完成，绝不中途更改。

写到这里，差异基因分析从原始数据到差异基因表达矩阵的链路差不多就完整了。这个流程里没有一步是可以用"全自动工具"跳过的"黑箱"，每一步都藏着值得仔细推敲的统计假设和生物学前提。我自己从最初只会跑DESeq2默认参数的小白，到现在每次拿到结果先做三遍"逆向检查"，中间踩过的坑远比这篇文章写的要多。如果你正在跑自己的第一批转录组数据，我的建议是别急着追求"一键汇总"，先把表达矩阵的质量检查、工具选型的理由、阈值设定的依据想清楚，后面再复杂的结果也只是水到渠成的事。