20 DESeq2 差异分析(上):从原始Counts到标准化的“公正比较”
20.1 【导语】万事之源:为何要这样做?
从featureCounts等工具中获得的原始Read Counts矩阵,其数值本身携带着两个根深蒂固的“固有偏倚”。第一个是“测序深度偏倚”:一个测序文库被测得越深(总Reads数越高),其中每一个基因的Counts值也天然就越高。第二个是“基因长度偏倚”:在表达水平相同的前提下,一个越长的基因,由于其拥有更长的“渔网”,能够捕获到的测序Reads也天然就越多。
因此,如果直接使用未经处理的原始Counts值来比较不同样本间或不同基因间的表达水平,就如同在完全不考虑拳击手体重级别的情况下,去比较他们的出拳力量——其结论将是毫无科学意义的。差异表达分析的第一步,也是统计上最关键的一步,就是必须进行“标准化”(Normalization),其核心目标是科学地估算并消除这些技术性偏倚,从而确保后续所有的统计比较,都是在一个“公平”的基准之上进行。
20.2 【核心实践】从原理到决策
20.2.1 【DESeq2 工作流的起点:构建DESeqDataSet对象】
DESeq2的所有分析,都始于构建一个名为DESeqDataSet的专用数据对象。我们可以用一个严谨的临床试验来做隐喻,成功构建这个对象,你需要提供三样不可或缺的东西:
countData:这是你的原始测量数据,即我们从featureCounts获得的、未经任何修改的整数计数矩阵。
colData:这是你的样本信息表(Sample Metadata),如同临床试验中每个病人的详细记录,必须包含样本的分组信息(如'Control', 'Treated')、批次、年龄等。
design:这是你的实验设计公式,例如~ group。它是在用代码,向DESeq2明确声明:“我本次分析最核心的科学问题,是比较colData中group这一列所定义的两组或多组样本之间的差异”。
## 建立 DESeqDataSet 对象
dds <- DESeqDataSetFromMatrix(
countData = count_matrix,
colData = sample_info,
design = ~ group
)
数据结构决策:在构建时,一个最容易出错但至关重要的细节是:colData的行名,必须与countData的列名,在名称和顺序上完全、严格地一致。这是将“基因表达数据”与“样品生物学属性”进行精确关联的唯一纽带,任何不匹配都将导致DESeq2拒绝工作或产生错误的结果。
20.2.2 【DESeq2标准化的核心思想:中位数比率法】
DESeq2的标准化方法,远比简单地将Counts除以总Reads数(这种方法对高表达基因的波动非常敏感,统计上不稳定)要复杂和稳健得多。其核心思想被称为“中位数比率法”(Median of Ratios)。
简化解释:DESeq2的标准化过程大致如下: 第一步,它会计算出一个“虚拟的参考样本”,这个参考样本的表达量,是每个基因在所有真实样本中表达量的几何平均值。 第二步,对于每一个真实样本,它会计算出其所有基因表达量与这个参考样本表达量的比率。 第三步,取这些比率的中位数,作为这个样本的“表达水平缩放因子”(Size Factor)。这个因子的本质,就是对该样本整体测序深度相对于所有样本平均水平的系统性偏差的精确估计。
结果诊断:在DESeq()函数运行后,我们可以通过sizeFactors(dds)来查看计算出的每个样本的缩放因子。一个理想的数据集中,这些因子应该在1左右波动。如果某个样本的因子显著偏离1(例如,大于2或小于0.5),这便是一个明确的诊断信号,暗示该样本的测序深度与其他样本存在巨大差异,你需要在后续的质量评估中(如PCA图)对这个样本进行特别的关注。
20.2.3 【方差稳定变换 (vst / rlog):可视化的“滤镜”】
后果驱动分析:对于原始的测序Counts数据,其方差与均值存在高度的正相关关系——即表达水平越高的基因,其Counts值的绝对波动范围也越大。如果直接在这种原始或仅经过Size Factor标准化的数据上进行探索性可视化分析(如PCA或热图),那么整个图形的结构将被少数几个表达量极高的“明星基因”所主导,而数万个中低表达基因的真实模式将被完全掩盖。
决策:vst()(方差稳定变换)或rlog()(正则化对数变换)的作用,就是为可视化目的,对数据施加一个特殊的数学“滤镜”。经过这个变换后,数据的方差将与其均值基本无关。经过变换后的数据,才是用于PCA、样本聚类、热图等探索性可视化的、统计上“正确”的输入数据。
20.3 【认知升维】常见的思维陷阱与对策
20.3.1 思维陷阱一:用FPKM/RPKM/TPM值进行差异分析
FPKM, RPKM, TPM等指标,虽然在其计算中同时考虑了基因长度和测序深度,看似是“完美”的标准化值,但大量的统计学研究已经证明,它们在进行样本间的统计比较时存在系统性偏差,并不适用于像DESeq2、edgeR这类基于负二项分布的、稳健的统计模型。
其对策是,必须建立一条不可逾越的铁律:DESeq2等主流差异分析工具的输入,必须是、且只能是未经任何标准化处理的、原始的整数Read Counts。所有必要的标准化步骤,都会由DESeq2内部,以一种统计上更为严谨和稳健的方式来完成。
20.3.2 思维陷阱二:过滤低表达基因的时机错误
新手为了“净化”数据,常常会在构建DESeqDataSet对象之前,就手动地、大刀阔斧地剔除了数据集中所有低Counts的基因。这是一个极其危险的操作。
其对策是,必须理解DESeq2的中位数比率标准化方法,其准确性高度依赖于一个核心假设:数据集中大多数的基因,在不同组间,都不是差异表达的。过早、过度地进行基因过滤,可能会破坏这一假设,从而影响Size Factor估算的准确性。正确的做法是,要么只进行一个非常宽松的预过滤(例如,保留那些在至少N个样本中Counts总和大于10的基因),要么完全不作预过滤,让DESeq2在后续差异分析results()步骤中,利用其内置的、基于统计功效的“独立过滤”方法来自动完成。
20.4 【总结与拓展】构建你的思维框架
我们必须将差异表达分析前的标准化过程,视为一次精密的、基于统计模型的“仪器校准”过程。你的目标,绝非简单地让不同样本的数字看起来“差不多大”,而是要利用稳健的统计学方法,科学地估算出由测序深度这类技术因素(而非真实的生物学因素)所引入的系统性偏差,并将其从数据中移除,从而为后续的、旨在发现真实生物学信号的假设检验铺平一条公正、无偏的道路。
基于此框架,请思考一个更高级的实验设计问题:假设你的实验设计中,除了你感兴趣的处理因素(group),还包含一个你已知的、不感兴趣的“混杂因素”,例如,你的样本来自于两个不同的测序批次(batch)。你将如何在DESeq2的design公式中,体现并主动校正这个已知的批次效应,以确保你最终得到的差异表达基因列表,是真正由你的group因素所导致的,而非简单地反映了不同批次之间的技术差异?
探索生命科学前沿,提升实战技能!欢迎微信搜索并加入「生信实战圈」,获取最新技术干货、实战案例与行业动态。 点击关注,与同行一起成长!
