做转录组分析的朋友，一定绕不开一个问题：手上有几组样本的测序数据，怎么把组间差异表达的基因筛出来，形成一张可以直接拿去做后续分析的差异基因表达矩阵。这句话听起来简单，但实际操作里涉及到的数据格式、统计模型、参数设置、结果筛选，每一步都有不少细节。这篇基于我大量实测经验的总结，完整走一遍从原始表达矩阵到差异基因表达矩阵的流程，包含工具选型逻辑、R代码、参数解释和踩坑记录，供刚开始接触差异分析的读者和想系统梳理流程的从业者参考。
1. 差异基因分析的整体思路与核心逻辑
1.1 差异分析到底在回答什么问题
差异基因分析的出发点是生物学问题：在某种处理、某个突变、某个疾病状态下，哪些基因的表达量发生了显著改变。落到数据层面，就是比较两组或多组样本之间的基因表达水平，找出差异足够大、且统计学上可信的基因集合，再把这些基因对应的表达量提取出来，形成一张清晰的矩阵表，供后续聚类、富集分析、甚至机器学习建模使用。
这里面有两个核心量需要理解清楚。一个是变化幅度，也就是差异倍数，通常用log2FoldChange表示。另一个是统计学显著性，体现为P值和经过多重假设检验校正后的调整后P值。二者缺一不可：有的基因差异倍数很大，但样本间波动也大，统计检验不显著；有的基因P值极小，但表达量差异倍数很小，生物学意义有限。所以绝大多数分析都会同时卡这两个阈值。
用生活化的例子理解：假设你想比较两个班级的数学成绩差异。差异倍数相当于两个班平均分之差，而P值相当于把班级内部成绩的波动也考虑进去后，判断这个分值差是不是只是随机波动造成的。如果一班平均分比二班高20分，但两个班内部学生成绩天差地别，那这个20分未必可信；反过来，两个班平均分只差5分，但每个学生成绩都很稳定，重复很多次都一样，那这5分反而可能是真实差异。差异基因分析背后的逻辑就是这个。
表达矩阵本身也不复杂。行是基因，列是样本，每个单元格是某种定量方式得到的表达值。比如featureCounts或者HTSeq-count输出的原始count矩阵。你拿到手的第一件事，就是把这种原始计数矩阵处理好，交给差异分析软件，最终得到一个包含基因ID、log2FoldChange、P值、调整后P值的统计结果表，再从中挑出显著基因，生成差异基因表达矩阵。
1.2 三大主流工具选型：DESeq2、edgeR、limma
目前做RNA-seq差异分析，主流工具基本稳定在三个R包：DESeq2、edgeR、limma。这三个我都实际用过，各自有不同的脾性。
先说道理。RNA-seq的原始count数据有一个显著特征：它是离散的整数计数，而且存在过度离散。什么意思呢？同一个基因在生物学重复样本之间的波动，往往比泊松分布期望的更大。如果直接用常规的t检验或者方差分析处理，容易得到太多假阳性。所以主流工具都基于负二项分布来建模，把生物学变异的离散度也估计进去。
DESeq2是目前使用最广泛的工具。它的特点是稳健、参数少、对用户友好，尤其适合样本量不大的实验设计。它内置了收缩估计逻辑，对低表达基因、离散度大的基因会做一定的收缩，避免出现极端偏离的倍数变化。默认流程只需要提供一个计数矩阵、一个分组信息，一行代码就能跑出结果，新手基本不会跑偏。
edgeR是另一个经典工具，历史更久，速度也更快。它同样基于负二项分布，但估计分散度的方法和DESeq2不同。edgeR对数据量大、样本数多的情况很友好，而且提供了灵活的实验设计矩阵函数，适合比较复杂的多因素设计。代价是它的参数和选项很多，刚接触时容易迷茫，需要仔细看懂每一步在做什么。
limma最初是为基因芯片数据分析开发的，核心优势是线性模型框架非常成熟，计算速度极快。后来通过voom转换，把RNA-seq的count数据转换为log2-CPM并估计均值-方差关系后，也可以稳定地应用于RNA-seq数据。目前limma-voom在测序数据处理中的表现也很受欢迎，尤其适合样本量足够大、有几组以上条件需要同时比较的设计。
三个工具在常规两两组别比较中，结果差异不会太大，但细节上各有取舍。我常用的建议是：如果你是第一次做差异分析，组别简单、样本量不大，优先选DESeq2；如果样本量很大，或者需要跑很多组两两比较，limma-voom的计算效率会让你省很多时间；如果实验设计比较复杂、有多个因素交叉，edgeR的广义线性模型有时更顺手。
| 工具 | 统计模型 | 优势 | 适合场景 |
|---|---|---|---|
| DESeq2 | 负二项分布+收缩估计 | 稳健、参数少、对低表达基因处理合理 | 常规两组或多组比较，新手友好 |
| edgeR | 负二项分布+经验贝叶斯 | 速度快、实验设计灵活 | 复杂设计、大批量比较 |
| limma-voom | 线性模型+voom权重 | 极快、多组比较稳定 | 样本量大、多组循环比较 |
不管选哪个，数据的预处理阶段都是相通的。下面重点说清楚输入数据怎么准备，这是整个流程里最容易出问题、也最容易被一带而过的地方。
2. 输入数据准备：从定量结果到可分析的表达矩阵
2.1 上游定量工具与标准输入格式
差异分析的前提是你手上已经有一张表达定量表。RNA-seq上游数据处理链路通常是：原始测序数据经过质控、比对到参考基因组/转录组，再用定量工具统计每个基因对应的reads数，得到表达矩阵。
常见的定量工具组合：STAR或HISAT2做比对，然后使用featureCounts或HTSeq-count统计基因计数；不比对直接定量的工具则是salmon或kallisto，它们直接将reads比对到转录本并进行定量，速度快，而且对样本间的序列偏好处理更温和。featureCounts生成的结果格式是典型的计数矩阵，行是基因，列是样本，单元格是整数。
实际操作中，我拿到手里的一般是featureCounts输出的矩阵文件，或者salmon生成的quant.sf文件，后者需要先用tximport导入，把转录本水平的定量汇总到基因水平。无论哪种来源，最终都要整理成一张标准的表达矩阵：
gene_id sample1 sample2 sample3 sample4 sample5 sample6 ENSG00000141510 1256 1098 1366 856 932 1021 ENSG00000171862 89 112 76 234 198 176这张矩阵的行名必须是唯一的基因ID，不能有重复；列名是样本名；单元格是原始整数计数。很多人习惯先做CPM、TPM标准化再去做差异分析，这是一个常见的误区。DESeq2、edgeR、limma-voom内部会自动处理文库大小差异，并自己估算合适的标准化因子。你直接传原始整数计数进R包即可，如果预先做了标准化反而会破坏统计模型的假设，导致结果失真。
与表达矩阵同步的，还有一张分组信息表，通常叫colData。它必须与表达矩阵的列一一对应。最简单的形式是两列：样本名、分组名。比如三例对照组和三例处理组。
2.2 表达矩阵检查清单与常见坑
数据准备阶段不仔细，后面每一步都可能连环报错。我花在数据检查上的时间，并不比跑差异分析本身少。这里整理一份血泪换来的检查清单。
第一，基因ID的格式和版本。不同来源的表达矩阵可能使用不同注释版本的基因ID。有的用ENSEMBL ID，有的用Symbol，有的用Entrez ID。做分析时统一用ENSEMBL ID作为行名最稳妥，因为它是稳定标识，不会有Symbol合并、缺失的问题，下游做富集分析时再转换回来即可。如果矩阵里混着不同版本的ID，第一件事就是统一。
第二，有没有重复行名。featureCounts输出时，如果注释文件有冗余，可能出现同一个基因ID出现在多行。读取到R后一定要先做去重检查。否则后面DESeq2构建对象时会直接报错。
第三，低表达基因的处理逻辑。原始矩阵中会存在大量在几乎所有样本中都没有reads或者只有个位数reads的基因。这些基因不仅没有分析价值，还会干扰多重检验校正的计算量。主流做法是跑差异分析前做一步低表达过滤。DESeq2默认流程中，会自动过滤掉那些几乎不可能成为显著差异基因的行，但如果你用的是results()提取时没有做任何过滤，也可以自己在预处理阶段设定一个阈值，比如至少有几个样本的count > 10才保留。
第四，分组因子的水平顺序。这是最容易被忽略、又最容易导致结果反过来的坑。在R中，分组列如果被识别为字符型向量，DESeq2内部会默认按照字母顺序排列水平，也就是说对照组和治疗组的顺序可能完全不在你的预期里。如果后面提取log2FoldChange时发现方向反了，需要回去看因子的水平顺序，用relevel函数把基线组设为第一个水平。
第五，样本名与分组表的匹配。这是检查清单里最基础也最重要的环节。我建议拿到数据后先做一次行列名的交叉检查，比如用all(colnames(count_matrix) == rownames(colData))验证，数学上完全一致再继续。一次粗心导致的前后错位，会直接让整个差异分析结果作废。
# 整理后的输出格式 all(rownames(colData) == colnames(count_matrix))除了这些，还有批次效应问题。如果对照组和处理组的样本是分两批测序的，批次信息往往会和分组信息混杂，比如前一批全是处理组、后一批全是对照组，这在RNA-seq分析里是很让人头疼的设计。差异分析模型如果不把批次因素放进去，很可能会把批次差异当作处理差异，导致异常的基因数量暴增。解决办法是：如果知道批次信息，在DESeq2的设计公式里同时加入批次和分组。
3. 核心实操：DESeq2筛选差异基因并生成表达矩阵
3.1 完整R代码流程
有了前面的检查和整理，现在可以进入真正的差异分析环节。我以DESeq2为主线，把从读取数据到输出差异基因表达矩阵的完整代码串一遍。
# 安装并加载包 if (!requireNamespace("BiocManager", quietly = TRUE)) { install.packages("BiocManager") } BiocManager::install("DESeq2") library(DESeq2) # 读取表达矩阵 count_matrix <- read.csv("count_matrix.csv", header = TRUE, row.names = 1, check.names = FALSE) # 读取分组信息 colData <- read.csv("colData.csv", header = TRUE, row.names = 1) colData$group <- factor(colData$group, levels = c("control", "treatment")) # 构建DESeqDataSet对象 dds <- DESeqDataSetFromMatrix( countData = count_matrix, colData = colData, design = ~ group ) # 跑差异分析主流程 dds <- DESeq(dds) # 提取结果 res <- results(dds, contrast = c("group", "treatment", "control"))这里做一个参数说明。DESeqDataSetFromMatrix的三个核心参数中，countData是原始计数矩阵，colData是分组表，design是统计模型公式。design = ~ group表示比较组间差异，此时差异分析的因子就是分组这一列，把所有组间表达差异归因于分组因素。如果数据中存在批次效应，比如测序文库批次，则写成design = ~ batch + group。这么做的意义是，将批次效应作为协变量纳入模型一起估计，让分组差异的估计不受到批次的干扰。
contrast参数指定对比的具体两组。DESeq2统计算法会先估测每个基因的均值-离散度关系，然后对所有基因分别拟合负二项广义线性模型，再对指定组合的log2FoldChange做统计检验，每一步的实际计算都藏在DESeq(dds)这一行里。运行结束后，results就得到了全基因组每个基因的统计结果。
3.2 结果表逐字段解读
跑完差异分析，res是一个DataFrame，每一行是一个基因，每一列是统计结果。很多同学看到这些列名就晕，这里逐个解释清楚。
baseMean：这个基因在所有样本中的平均表达量，这里用的是标准化后的均值。它的意义是衡量基因整体表达水平，过滤掉低表达基因时可以参考这一列。
log2FoldChange：处理组相对于对照组表达量的对数化变化倍数。正数表示处理组上调，负数表示处理组下调。举个例子，处理组表达量是对照组的2倍，log2FoldChange就是1；处理组是对照组的0.5倍，就是-1。实际筛选时，一般要求绝对值大于等于1，相当于变化倍数至少2倍。
lfcSE：log2FoldChange的标准误，反映了这个倍数估计的不确定性。它越大，说明样本间波动越大，可信度越低。
stat：检验统计量，由log2FoldChange除以标准误得到。它被用于后续的P值计算。
pvalue：原始的统计学检验P值。但这里要注意，做了几万个基因的检验，如果只用原始P值挑基因，大量假阳性会混进来。
padj：多重检验校正后的P值，即BH校正后的结果。这是筛选差异基因最关键的指标。常用阈值是小于0.05，严格一些可以设为0.01。
| 列名 | 含义 | 筛选使用时注意 |
|---|---|---|
| baseMean | 标准化后的平均表达量 | 可作为低表达过滤参考 |
| log2FoldChange | 对数化的表达变化倍数 | 正负代表上下调方向 |
| lfcSE | 倍数变化的标准误 | 越大越不可信 |
| stat | 检验统计量 | 一般不用直接筛选 |
| pvalue | 原始P值 | 有假阳性风险 |
| padj | BH校正后P值 | 标准筛选指标 |
3.3 差异基因表达矩阵的生成与导出
得到全基因组的统计表后，筛选差异基因的方法本身很简单：取padj < 0.05且|log2FoldChange| >= 1。但更关键的是如何把这个筛选结果与原始表达矩阵结合，生成一张差异基因表达矩阵。这一步直接决定后续做热图、聚类、富集分析的输入质量。
在我的实操习惯里，差异基因表达矩阵生成分三步。
第一步，从res中筛出显著上下调基因的ID。注意此时res还没有和原始表达矩阵合并，我们需要先根据行名提取。
res_df <- as.data.frame(res) res_df$gene_id <- rownames(res_df) # 筛选显著基因 up_genes <- rownames(res_df[!is.na(res_df$padj) & res_df$padj < 0.05 & res_df$log2FoldChange > 1, ]) down_genes <- rownames(res_df[!is.na(res_df$padj) & res_df$padj < 0.05 & res_df$log2FoldChange < -1, ]) deg_genes <- c(up_genes, down_genes) length(deg_genes)为什么不用is.na直接过滤掉那些padj是NA的基因？在DESeq2结果中，有一部分基因会因为表达量极低或者离散度估计问题，导致padj列是NA。这些基因无法被可靠地评估显著性，筛选时要直接排除。如果筛选时忘了去掉NA，结果集里会混入无意义的数据。
第二步，提取这些基因对应的标准化表达矩阵。DESeq2内部有一个函数vst和rlog，它们会对count做方差稳定化转换，得到的矩阵更适合做热图和聚类分析。如果只是取标准化count，则用counts(dds, normalized=TRUE)。对于下游分析来说，方差稳定化后的矩阵是更好的选择。
# 获取方差稳定化表达矩阵 vsd <- vst(dds, blind = FALSE) vsd_mat <- assay(vsd) # 提取差异基因对应的表达矩阵 deg_matrix <- vsd_mat[deg_genes, ] dim(deg_matrix)第三步，记录完整信息并导出。我一般会保存三样东西：deg_matrix本身（用于热图、聚类）、包含统计量的注释表（用于筛选依据记录）、以及用于下游GSEA等工具的排序基因列表。保存用write.csv即可，行名保留为基因ID。
# 保存差异基因表达矩阵 write.csv(deg_matrix, "deg_expression_matrix.csv", quote = FALSE) # 保存差异分析完整结果 write.csv(res_df[order(res_df$padj), ], "all_deg_results.csv", quote = FALSE)这一步完成后，你就得到了一张干净、可直接用于热图绘制和下游分析的差异基因表达矩阵。矩阵的行名是显著差异基因，列名是样本名，数值是方差稳定化后的表达量。后续无论是画热图、做GO/KEGG富集分析，还是把基因列表拉出来做蛋白互作网络，你都拥有了一个可靠的下游分析起点。
4. 结果验证与可视化：PCA、火山图与热图实操
4.1 先看整体：PCA判断分组与批次效应
在跑差异分析之前或者结果出来之后，我都会先画PCA图。PCA的核心作用是看一眼整体表达谱的维度结构：样本是否能按分组分开，是否存在明显的离群样本，是否暴露出批次结构。这一步对验证上游数据质量和判断设计效果很重要。
代码上直接用DESeq2的vst结果加ggplot2即可。
library(ggplot2) library(DESeq2) # 计算PCA pca_data <- plotPCA(vsd, intgroup = "group", returnData = TRUE) percentVar <- round(100 * attr(pca_data, "percentVar")) ggplot(pca_data, aes(x = PC1, y = PC2, color = group)) + geom_point(size = 3) + xlab(paste0("PC1: ", percentVar[1], "% variance")) + ylab(paste0("PC2: ", percentVar[2], "% variance")) + theme_classic()解读PCA图的经验法则，第一看分组是否在某个主成分方向上明显分离。常见的情况是PC1只体现了对照组和处理组之间的差异，这是好事。第二看同组样本是否聚在一起，如果有样本跑到另一组区域，就要警惕这个样本是否已经污染或者标记错误。第三看PC1和PC2的方差占比。如果PC1只有20%以下，而PC2上有30%以上，且PC2对应的不是你的分组信息，那要小心是否存在强批次效应。
如果PCA图中对照组和处理组完全混在一起，没有分离趋势，有可能是组间本身差异很小，也有可能是批次效应盖过了真实的生物学信号。这时候可以先检查colData里是否有批次信息列，然后改用design = ~ batch + group再跑一次DESeq2，重新观察PCA。
4.2 火山图：差异基因的全局视图
火山图是展示差异分析结果的经典方式。横轴是log2FoldChange，纵轴是-log10(padj)，每个点是一个基因。图上的点分布能让人一眼确认上下调基因数量和变化幅度是否合理。
# 构建画图数据 plot_df <- as.data.frame(res) plot_df$gene_id <- rownames(plot_df) plot_df$significant <- "not significant" plot_df$significant[!is.na(plot_df$padj) & plot_df$padj < 0.05 & plot_df$log2FoldChange > 1] <- "up" plot_df$significant[!is.na(plot_df$padj) & plot_df$padj < 0.05 & plot_df$log2FoldChange < -1] <- "down" ggplot(plot_df, aes(x = log2FoldChange, y = -log10(padj), color = significant)) + geom_point(size = 0.8, alpha = 0.6) + scale_color_manual(values = c("down" = "blue", "not significant" = "grey", "up" = "red")) + labs(x = "log2(Fold Change)", y = "-log10(adjusted P-value)") + theme_classic()画火山图有一个经常忽略的小细节：padj为0的情况。DESeq2计算出的padj可能存在极小的值，在R里取-log10时变成Inf，画出图后这些点会跑到图外，看起来像是横坐标附近的异常点。处理办法是先把Inf值替换成一个足够大的数字，比如统一封顶到1e-300，或者直接对这些点的纵坐标做一个截断。
4.3 热图：差异基因表达矩阵的直接呈现
热图是最直观展示差异基因表达矩阵的方式。用pheatmap包基于前面提取的deg_matrix绘制即可。
library(pheatmap) # 按组排列样本 annotation_col <- data.frame(group = colData$group) rownames(annotation_col) <- colnames(deg_matrix) # 绘制热图 pheatmap( deg_matrix, annotation_col = annotation_col, scale = "row", show_rownames = FALSE, cluster_cols = TRUE, cluster_rows = TRUE, color = colorRampPalette(c("blue", "white", "red"))(100) )关于热图的参数，scale = "row"表示对每个基因做Z-score标准化，也就是把每个基因在各样本里的表达量变换成均值为0、标准差为1的状态，这样不同表达量范围的基因才能在同一张图上比较。如果不做这个标准化，表达量高的基因占据绝对主导，低表达但存在明显差异的基因会在图上完全看不出来。
聚类时，show_rownames设置是否显示基因名。如果筛选出的基因数量多到几百上千，全部显示基因名只会让图片糊成一片。我的习惯是保留注释但不显示行名，等选出兴趣基因后再单独画一张标了基因名的精简热图。
4.4 差异基因表达矩阵的下游应用
生成deg_matrix之后，最直接的应用是用于聚类和富集分析。
在聚类层面，我们可以基于差异基因表达矩阵对各样本重新聚类，观察处理组和对照组是否形成各自的分支，也可以顺便检查是否有样本因为离群而聚到另一边。
在富集分析层面，可以拿up_genes和down_genes分别去做GO和KEGG富集，看差异基因主要集中在哪些生物学过程、分子功能或信号通路里。常用工具包括clusterProfiler包、在线数据库等。另外，如果把所有基因按照log2FoldChange从大到小排序，而不是只保留显著差异基因，则可以进一步做GSEA分析，识别那些整体趋势上富集、但单个基因变化并不足够显著的通路。差异基因表达矩阵本身不会直接参与GSEA，但基于差异分析结果产生的排序列表是GSEA的标准输入。
还有一个实际用途是用于机器学习特征筛选。差异基因表达矩阵可以直接作为分类器的输入特征，筛选出的基因可以进一步做特征选择，找到最能区分组别的标志物组合。我做模拟数据练习时，直接用top差异基因作为输入，也能达到不错的分类效果。这进一步说明，从差异分析到差异基因表达矩阵这一步，是整个转录组数据分析链条里承上启下的关键环节。
5. 常见问题与排查技巧实录
5.1 高频报错与解决思路
差异分析流程中，我遇到过不少报错，这里挑几个典型的分享。
DESeqDataSetFromMatrix报错：行名有重复。出现这种情况，多半是注释文件本身有重复基因ID。解决方式是检查count_matrix的行名，用make.unique或者在导入之前按基因ID进行去重汇总。
results提取时提示contrast不存在。这通常是因为分组列虽然在colData里，但没有被正确识别为因子，或者因子的水平里没有对应的组名。检查方法是运行colData(dds)$group看它有没有levels。
padj大量为NA。前面提到过，这通常意味着这些基因的表达量极低，或者离散度估计异常，没有足够的信息判断差异是否显著。处理方式是不做强求，直接在筛选标准中剔除NA即可。如果core信息受损，可以尝试对dds重新运行estimateDispersions，观察图形是否异常。
log2FoldChange方向反了。这个问题我踩过不少次。根因是分组因子水平顺序，DESeq2默认把第二个水平作为分子，第一个水平作为分母来计算倍数变化。如果你希望对照组是基线，处理组是变化组，需要确保对照组是因子第一个水平，或者在results函数里用contrast参数显式指定两个水平。
样本间基因名字符串格式不一致。比如一部分是ENSG开头，一部分是版本号后缀，比如ENSG00000141510.17。在整理矩阵时要去掉版本号后缀，统一为无版本号的ENSEMBL ID。
下面是高频问题速查：
| 现象 | 可能原因 | 处理方式 |
|---|---|---|
| 行名重复 | 注释文件冗余 | 按基因ID去重汇总 |
| contrast错误 | 分组不是因子 | 转为factor并检查levels |
| padj大量NA | 低表达基因 | 筛选时剔除NA |
| 倍数方向相反 | 因子水平顺序 | relevel或contrast显式指定 |
| ID带版本号 | 注释格式不一致 | 去掉小数点及后方字符 |
5.2 实操心得：新手最容易忽略的关键点
这里多说几句经验。做完差异分析不等于分析结束，下载结果表格前一定要先做一次手动抽查。我会在结果表里随机挑几个逻辑上应该显著、生物学上合理的基因，去原始表达矩阵里手动核实表达量确实是跨组变化。这个习惯可以避免最终写报告时才发现整个分析的对比方向错了。
关于生物重复数量，我强烈建议至少做3个以上。从统计功效的角度讲，2个重复也能跑DESeq2，但离散度估计会非常不稳定，结果里的大量基因padj不可信。如果有条件，4个重复更稳。RNA-seq实验设计阶段如果能在重复数上多投入样本，后续差异分析的可靠性会得到极大提升。我自己曾经为了省样本用2个重复跑过一次，结果跑出的显著基因数量非常不稳定，换一种滤波方式就变化巨大，完全没法写结论。
关于上游定量工具的一致性。差异分析要求所有样本的定量结果在同一个参考序列和同一个工具下产生。如果对照组用的是某公司A的定量结果，处理组用的是公司B的定量结果，那么两个组的计数差异很可能是技术差异而非生物学差异。这一点务必在拿到数据时检查清楚。
关于基因的log2FoldChange绝对值大小，我会特别提醒一点：实际数据里很多真实差异基因的log2FoldChange只有0.5到1左右，严格卡1倍阈值可能损失部分真实信号。如果后续富集分析提示某些信号通路很活跃但差异基因数量偏少，可以尝试放宽阈值到|log2FoldChange| >= 0.585，也就是1.5倍变化，再看看富集结果是否更加合理。
最后再分享一个写文章和复现时非常实用的习惯。跑完一次完整分析后，我会立刻保存sessionInfo()和R脚本。sessionInfo记录了所有R包的版本号。转录组分析工具版本之间常有API变化，数月后重新运行脚本可能因为包版本升级而报错。保留sessionInfo相当于给自己的分析留下了一个可精确复现的环境快照。以后无论是自己回溯结果，还是被要求复核分析，都能省下大量排查时间。
实际做数据分析和交付报告时，我还习惯把整个流程的每个参数和版本都写清楚，比如DESeq2的版本、参考注释文件的来源、筛选阈值的选择依据。一方面这是分析合规性的体现，另一方面也是避免未来被质疑结果可靠性的重要防线。生物信息分析项目的可重复性，某种意义上比分析方法本身还要重要。
