引言
本系列 将开展全新的转录组分析专栏,主要针对使用 DESeq2
时可能出现的问题和方法进行展开描述。想要学习更多内容可以添加文末的学习交流群或客服QQ。
计数数据转换
在进行差异表达分析时,我们会基于原始计数数据,运用上一节提到的离散分布方法。但对于像可视化、聚类这类后续分析,使用经过转换的计数数据会更有优势。
对数转换是最常见的转换方式。因为基因在某些条件下计数值可能为零,而在其他条件下又不为零,所以有人建议引入伪计数,也就是采用这种转换形式:
这里的 n 是计数值,n₀ 是一个正数。
本节中,我们探讨了两种更有理论依据且能合理选择类似 n₀ 参数的替代方法。一种基于方差稳定化转换(VST)的概念,另一种是正则化对数(rlog),它考虑了样本差异的先验信息。这两种转换都能得到在 log₂ 尺度上的数据,且这些数据已经根据文库大小等归一化因素进行了调整。
VST 和 rlog 这两种转换的核心目的是消除方差与均值之间的关联,尤其是当均值较低时,计数数据对数的方差会变得很大。它们都借助实验中整体的方差随均值变化的趋势,来转换数据,以去除这种整体趋势。需要注意的是,我们并不强求所有基因转换后的方差完全相同。在下面的图表里,你会看到转换后,相同均值的基因标准差并不完全一样,但整体趋势已经变平缓了。那些行方差高于整体趋势的基因,有助于我们将样本划分成有意义的组别。
盲方差估计
vst 和 rlog 这两个函数都有一个 blind 参数,用来决定转换时是否忽略设计公式中指定的样本信息。当 blind 为 TRUE(默认设置)时,函数会只用截距来重新估计离散度。这种设置适用于完全不考虑实验分组信息来比较样本,比如做样本质量保证(QA)时用到。
不过,如果预计很多基因(行)的计数差异是由实验设计能解释的,并且想把数据转换用于后续分析,那盲方差估计就不合适了。这种情况下,盲方差估计会让离散度估计值变大,因为它把实验设计导致的差异当成了不需要的噪声,进而使转换后的数值过度相互靠近。把 blind 设置为 FALSE,就会用已经估计好的离散度来转换,或者如果没有,就用当前设计公式来估计。不过,转换时只用到均值 - 离散度趋势线上的拟合离散度估计(也就是整个实验中离散度对均值的整体依赖关系)。所以,把 blind 设置为 FALSE,大部分情况下还是没用到样本分组的信息来进行转换。
提取转换值
这些转换函数会返回一个 DESeqTransform 类的对象,它是 RangedSummarizedExperiment 类的子类。对于大约 20 个样本,在新创建的 DESeqDataSet 上,rlog 大概要 30 秒,vst 不到 1 秒。如果使用 blind = FALSE,或者已经运行过 DESeq 函数,运行时间会缩短,因为不用重新估计离散度值了。用 assay 函数可以提取出归一化值的矩阵。
vsd <- vst(dds, blind=FALSE)
rld <- rlog(dds, blind=FALSE)
head(assay(vsd), 3)
方差稳定化转换
之前我们用参数拟合方法来估计离散度。在这种情况下,vst 函数会用方差稳定化转换的封闭形式表达式。要是用局部拟合(在 estimateDispersions 里设置 fitType="locfit" 选项),那就会用数值积分来替代。转换后的数据应该是方差大致稳定的,并且也考虑了大小因子或归一化因子的校正。对于计数较大的数据,转换后的结果是在 log₂ 尺度上的。
正则化对数转换
rlog 函数的意思是正则化对数,它通过拟合一个模型来把原始的计数数据转换到 log₂ 尺度。这个模型包含每个样本的项,还有从数据中估计出来的系数的先验分布。这和 DESeq 以及 nbinomWaldTest 用的对数倍数变化的收缩方法(有时也叫正则化或调和)是一个路子。转换后得到的数据里包含这样的元素:
qij 是个参数,和基因 i、样本 j 的预期真实片段浓度成正比(具体公式下面有),βi₀ 是个不收缩的截距,βij 是样本特有的效应,会根据整个数据集的离散度 - 均值趋势向零收缩。通常低计数的基因离散度高,所以这些基因在 rlog 转换中收缩得更厉害。
得注意的是,qij 表示在把大小因子 sj 除掉后均值 μij 的一部分,这就说明 rlog 转换本身就能考虑到测序深度的不同。要是没有先验,这个设计矩阵就会有多个解,但给非截距的 beta 加上先验就能找到唯一的解。
转换对方差的影响
下面这张图展示了用移位对数转换、正则化对数转换和方差稳定化转换后的数据标准差,和样本均值的关系。移位对数在计数少的时候标准差会变高,正则化对数稍微好点,而方差稳定化转换后的数据,标准差在整个范围内基本保持不变。
还要注意,这种图的纵轴是所有样本方差的平方根,也就是把实验条件导致的方差也算进去了。虽然从均值上看,方差平方根的平坦曲线好像是转换追求的效果,但对于那些因实验条件有很多真实差异的数据集,这可能就不合适了。
# this gives log2(n + 1)
ntd <- normTransform(dds)
library("vsn")
meanSdPlot(assay(ntd))
meanSdPlot(assay(vsd))
meanSdPlot(assay(rld))