DESeq2差异基因分析和批次效应移除

  • A+
所属分类:Transcriptomics

基因表达标准化

不同样品的测序量会有差异,最简单的标准化方式是计算 counts per million (CPM),即原始reads count除以总reads数乘以1,000,000。

这种计算方式的缺点是容易受到极高表达且在不同样品中存在差异表达的基因的影响;这些基因的打开或关闭会影响到细胞中总的分子数目,可能导致这些基因标准化之后就不存在表达差异了,而原本没有差异的基因标准化之后却有差异了。

RPKM、FPKM和TPM是CPM按照基因或转录本长度归一化后的表达,也会受到这一影响。

为了解决这一问题,研究者提出了其它的标准化方法。

量化因子 (size factor, SF)是由DESeq提出的。其方法是首先计算每个基因在所有样品中表达的几何平均值。每个细胞的量化因子(size factor)是所有基因与其在所有样品中的表达值的几何平均值的比值的中位数。由于几何平均值的使用,只有在所有样品中表达都不为0的基因才能用来计算。这一方法又被称为 RLE (relative log expression)

上四分位数 (upperquartile, UQ)是样品中所有基因的表达除以处于上四分位数的基因的表达值。同时为了保证表达水平的相对稳定,计算得到的上四分位数值要除以所有样品中上四分位数值的中位数。

TMM是M-值的加权截尾均值。选定一个样品为参照,其它样品中基因的表达相对于参照样品中对应基因表达倍数的log2值定义为M-值。随后去除M-值中最高和最低的30%,剩下的M值计算加权平均值。每一个非参照样品的基因表达值都乘以计算出的TMM。这个方法假设大部分基因的表达是没有差异的。

DESeq2 差异基因鉴定一步法

为了简化差异基因的运算,易生信做了脚本封装,DESeq2.sh,只需提供原始基因表达矩阵、样品分组信息表即可进行差异基因分析和鉴定。

具体使用

运行结束即可获得以下结果文件

DESeq2 差异基因鉴定分步法

安装包 DESeq2安装方法如下

A distributional assumption is needed because we want to estimate the probability of extreme events (large fold change just appearing by chance) from limited replicates. The negative binomial (a.k.a. Gamma-Poisson) is a good choice for RNA-seq experiments because

  • The observed data at gene level is inherently counts or estimated counts of fragments for each feature and
  • The spread of values among biological replicates is more than given by a simpler, one parameter distribution, the Poisson; and it seems to be captured by the NB sufficiently well

加载包

读入数据

产生DESeq数据集

DESeq2包采用DESeqDataSet存储原始的read count和中间计算的统计量。

每个DESeqDataSet对象都要有一个实验设计formula,用于对数据进行分组,以便计算表达值的离散度和估计表达倍数差异,通常格式为~ batch + conditions (为了方便后续计算,最为关注的分组信息放在最后一位)。

countData: 表达矩阵

colData: 样品分组信息表

design: 实验设计信息,conditions必须是colData中的一列

获取标准化后的数据

根据基因在不同的样本中表达变化的差异程度mad值对数据排序,差异越大的基因排位越前。

样品层级聚类分析,判断样品的相似性和组间组内差异

DESeq2差异基因分析和批次效应移除

样品PCA分析

差异基因分析

results函数提取差异基因分析结果,包含log2 fold changesp valuesadjusted p valuesconstrast用于指定比较的两组的信息,输出的log2FoldChange为log2(SampleA/SampleB)

DESeq2的原始输出结果增加样品平均表达信息,使得结果更容易理解和解析。

整体分析结果输出到文件

提取差异表达基因

绘制火山图

DESeq2差异基因分析和批次效应移除

差异基因热图

DESeq2差异基因分析和批次效应移除

关于批次效应

每个DESeqDataSet对象都要有一个实验设计formula,用于对数据进行分组,以便计算表达值的离散度和估计表达倍数差异,通常格式为~ batch + conditions (为了方便后续计算,最为关注的分组信息放在最后一位)。

如果记录了样本的批次信息,或者其它需要抹除的信息可以定义在design参数中,在下游回归的分析中会根据design formula来估计batch effect的影响,并在下游分析时减去这个影响。这是处理batch effect的推荐方式。

在模型中考虑batch effect并没有在数据矩阵中移除bacth effect,如果下游处理时,确实有需要可以使用limma包的removeBatchEffect来处理。

countData: 表达矩阵

colData: 样品分组信息表

design: 实验设计信息,batchconditions必须是colData中的一列

Just to be clear, there’s an important difference between removing a batch effect and modelling a batch effect. Including the batch in your design formula will model the batch effect in the regression step, which means that the raw data are not modified (so the batch effect is not removed), but instead the regression will estimate the size of the batch effect and subtract it out when performing all other tests. In addition, the model’s residual degrees of freedom will be reduced appropriately to reflect the fact that some degrees of freedom were “spent” modelling the batch effects. This is the preferred approach for any method that is capable of using it (this includes DESeq2). You would only remove the batch effect (e.g. using limma’s removeBatchEffect function) if you were going to do some kind of downstream analysis that can’t model the batch effects, such as training a classifier.

批次效应模拟

DESeq2差异基因分析和批次效应移除

DESeq2差异基因分析和批次效应移除

SVA(批次未记录时,寻找潜在影响因子,并矫正)

参考资料

发表评论

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen: