edgeR提供的TMM归一化算法详解

edgeR提供的TMM归一化算法详解

raw count的定量方式，是无法直接在样本间进行比较的。所以差异分析时，都会对原始的表达量数据进行归一化。

造成raw count无法直接比较的因素有很多，最常见的有以下两个因素
1.测序量
2.RNA的组成

测序量对count的影响很好理解，测序的数据量越大，对应的reads也就越多。RNA的组成是如何影响表达量的呢？

由于RNA的组织特异性，时间特异性等因素，我们无法保证两个样本中表达的RNA的种类和数量完全相同。假设两个样本A和B, B中的RNA的种类是A的两倍，共有的RNA表达量相同，在相同测序量的情况下，共有的RNA在A中的表达量会是B中的两倍，由此可见，不同样本RNA的构成也会对检测到的RNA表达量造成影响。

归一化时，通常的做法是只考虑样本间相同的RNA, 在此基础上，再消除测序量的影响。

DESeq2的归一化算法只考虑在所有样本中表达量都大于零的基因，也是出于相同RNA构成的考虑。edgeR采取了参照样本的策略，首先从所有样本中挑选一个样本作为参照，在对其他样本进行归一化时，只考虑哪些在参照样本和待归一化的样本间共有的RNA。

选取参照样本的代码如下

y   <- t(t(data)/lib.size)
f75 <- apply(y,2,function(x) quantile(x,p=p0.75))
refColumn <- which.min(abs(f75-mean(f75)))

根据相对丰度，计算每个样本的第三四分位数，采用该值代表每个样本的表达水平，选取与所有样本第三四分位数均值相差最小的样本作为参照样本。
参照样本选好之后，采用循环对每个样本进行归一化。在归一化时，重点关注基因的选取。代码如下

obs <- as.numeric(obs)
ref <- as.numeric(ref)
nO <- sum(obs)
nR <- sum(ref)
logR <- log2((obs/nO)/(ref/nR))    
absE <- (log2(obs/nO) + log2(ref/nR))/2
v <- (nO-obs)/nO/obs + (nR-ref)/nR/ref
fin <- is.finite(logR) & is.finite(absE) & (absE > -1e10)
logR <- logR[fin]
absE <- absE[fin]
v <- v[fin]

ref代表参照样本的表达量，obs代表待归一化样本的表达量。通过构建的3个统计量对基因进行初步过滤，logR 其实就是样本间的log2FD值，absE是表达量的均值，通过fin指定的过滤措施，过滤掉在任意样本中表达量为的基因，通过absE过滤掉一部分表达量很低的基因。

在此基础上，分别从头尾再去除部分基因，代码如下

n <- length(logR)
loL <- floor(n * 0.3) + 1
hiL <- n + 1 - loL
loS <- floor(n * 0.05) + 1
hiS <- n + 1 - loS
keep <- (rank(logR)>=loL & rank(logR)<=hiL) & (rank(absE)>=loS & rank(absE)<=hiS)

这样做的目的是保留在样本间表达量没有差异的基因。最后在利用下列公式计算每个样本的sizefactor

f <- sum(logR[keep]/v[keep], na.rm=TRUE) / sum(1/v[keep], na.rm=TRUE)
2^f

相比DESeq2的归一化算法，TMM算法基于表达量没有差异的基因来进行归一化。