|
我們都知道raw count的定量方式，是無法直接在樣本間進行比較的。所以差異分析時，都會對原始的表達量數(shù)據(jù)進行歸一化。 在之前的文章中，我們介紹了DESeq2提供的歸一化算法，本章介紹下edgeR的TMM歸一化算法。 造成raw count無法直接比較的因素有很多，最常見的有以下兩個因素 測序量對count的影響很好理解，測序的數(shù)據(jù)量越大，對應的reads也就越多。RNA的組成是如何影響表達量的呢？ 由于RNA的組織特異性，時間特異性等因素，我們無法保證兩個樣本中表達的RNA的種類和數(shù)量完全相同。假設兩個樣本A和B, B中的RNA的種類是A的兩倍，共有的RNA表達量相同，在相同測序量的情況下，共有的RNA在A中的表達量會是B中的兩倍，由此可見，不同樣本RNA的構(gòu)成也會對檢測到的RNA表達量造成影響。 歸一化時，通常的做法是只考慮樣本間相同的RNA, 在此基礎上，再消除測序量的影響。 DESeq2的歸一化算法只考慮在所有樣本中表達量都大于零的基因，也是出于相同RNA構(gòu)成的考慮。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)))根據(jù)相對豐度，計算每個樣本的第三四分位數(shù)，采用該值代表每個樣本的表達 參照樣本選好之后，采用循環(huán)對每個樣本進行歸一化。在歸一化時，重點關注基因的選取。代碼如下 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]
在此基礎上，分別從頭尾再去除部分基因，代碼如下 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算法基于表達量沒有差異的基因來進行歸一化。 ·end· |
|
|
