前言
DESeq2的LRT又称之为似然比检验,用于检验跨两个以上组别评估表达变化,比方说以多个时间点作为分组,
似然比检验简介
似然比检验是用于研究你的两个统计学模型是否有差异的一种检验方式,其基本模型如下:
原假设 H0:θ = θ0;备择假设为 H1 : θ ≠(不等于) θ0(θ = θ1)
从模型中可以看到事实上 θ0 和 θ1 可以认为是代表了两个不同的模型,其含义是你有两个统计学模型,分别是 p(x; θ1) 和 p(x; θ0) 。λ 越接近 1 代表两个模型差异越小;反之,两个模型差异越大
在R里面实现LRT:
library(epicalc)
model0 <- glm(case ~ induced + spontaneous, family=binomial, data=infert)
model1 <- glm(case ~ induced, family=binomial, data=infert)
lrtest (model0, model1)
##结果
Likelihood ratio test for MLE method
Chi-squared 1 d.f. = 36.48675 , P value = 0
在R的代码里面,我们可以发现 LRT 事实上是对两个模型进行比较,来看两个模型之间的差异
对于回归模型来说:
我们对回归模型的似然值定义如下图
下面的误差项定义为真实值与预测值之间的差异
因此,不同的回归模型,每一项的误差 ei 是不一样的,因此基于误差项ei求出来的似然值也是不一样的,那么似然值越高,代表模型的误差越小,拟合效果越好
DESeq2里面的LRT
DESeq2里面也提供了LRT的检验方法,往往用于两组以上的比较,我们看官网里面的实例,改编自官网:
library("fission")
data("fission")
#以minute(时间)来建模
ddsTC <- DESeqDataSet(fission, ~ minute )
#作为比较的模型,以 1 作为对照
ddsTC <- DESeq(ddsTC, test="LRT", reduced = ~ 1)
resTC <- results(ddsTC)
resTC$symbol <- mcols(ddsTC)$symbol
head(resTC[order(resTC$padj),], 4)
结果为:
注意,这里的log2FoldChange没有实际的意义
这里以时间作为多组的分类变量,以时间为自变量来建模;参数 ~ minute 表示的是对基因表达量随时间变化来建模;而参数 ~ 1 表示的是对照,即基因表达量和时间之间没有任何线性关系
显然对于某基因来说,左图更加match,可以看出随时间变化的趋势 (时间点可以不止两个,允许有很多个时间点来建立与 gene g 表达量之间的线性模型;而LRT的思想就是两个模型进行比较,看哪个模型与你的数据更match)
这是三个时间点的情况,显然左边的线性模型相比于右边的线性模型误差更小,故左边的模型更match此数据,即该基因的表达是随时间的变化而变化
当LRT检测的p值越显著,说明这两个模型之间的差异越大,而结果中的行名代表不同的基因,p值越小,则代表接受H1(备择假设),即参数 ~ minute 建立的模型更适合你的数据,也就意味着基因表达与时间这个变量相关
将p值显著的基因挑出来,那么这些基因就是随时间变化而波动的基因了(相对于参数 ~ 1 的模型来说)