edgeR的TMM歸一化算法怎么使用

本篇內(nèi)容主要講解“edgeR的TMM歸一化算法怎么使用”，感興趣的朋友不妨來看看。本文介紹的方法操作簡單快捷，實用性強。下面就讓小編來帶大家學(xué)習(xí)“edgeR的TMM歸一化算法怎么使用”吧!

我們都知道raw count的定量方式，是無法直接在樣本間進(jìn)行比較的。所以差異分析時，都會對原始的表達(dá)量數(shù)據(jù)進(jìn)行歸一化。

造成raw count無法直接比較的因素有很多，最常見的有以下兩個因素
1.測序量
2.RNA的組成

測序量對count的影響很好理解，測序的數(shù)據(jù)量越大，對應(yīng)的reads也就越多。RNA的組成是如何影響表達(dá)量的呢？

由于RNA的組織特異性，時間特異性等因素，我們無法保證兩個樣本中表達(dá)的RNA的種類和數(shù)量完全相同。假設(shè)兩個樣本A和B, B中的RNA的種類是A的兩倍，共有的RNA表達(dá)量相同，在相同測序量的情況下，共有的RNA在A中的表達(dá)量會是B中的兩倍，由此可見，不同樣本RNA的構(gòu)成也會對檢測到的RNA表達(dá)量造成影響。

歸一化時，通常的做法是只考慮樣本間相同的RNA, 在此基礎(chǔ)上，再消除測序量的影響。

DESeq2的歸一化算法只考慮在所有樣本中表達(dá)量都大于零的基因，也是出于相同RNA構(gòu)成的考慮。edgeR采取了參照樣本的策略，首先從所有樣本中挑選一個樣本作為參照，在對其他樣本進(jìn)行歸一化時，只考慮哪些在參照樣本和待歸一化的樣本間共有的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ù)，采用該值代表每個樣本的表達(dá)
水平，選取與所有樣本第三四分位數(shù)均值相差最小的樣本作為參照樣本。

參照樣本選好之后，采用循環(huán)對每個樣本進(jìn)行歸一化。在歸一化時，重點關(guā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]

ref代表參照樣本的表達(dá)量，obs代表待歸一化樣本的表達(dá)量。通過構(gòu)建的3個統(tǒng)計量對基因進(jìn)行初步過濾，logR 其實就是樣本間的log2FD值，absE是表達(dá)量的均值，通過fin指定的過濾措施，過濾掉在任意樣本中表達(dá)量為的基因，通過absE過濾掉一部分表達(dá)量很低的基因。

在此基礎(chǔ)上，分別從頭尾再去除部分基因，代碼如下

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)

這樣做的目的是保留在樣本間表達(dá)量沒有差異的基因。最后在利用下列公式計算每個樣本的sizefactor

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

相比DESeq2的歸一化算法，TMM算法基于表達(dá)量沒有差異的基因來進(jìn)行歸一化。

到此，相信大家對“edgeR的TMM歸一化算法怎么使用”有了更深的了解，不妨來實際操作一番吧！這里是億速云網(wǎng)站，更多相關(guān)內(nèi)容可以進(jìn)入相關(guān)頻道進(jìn)行查詢，關(guān)注我們，繼續(xù)學(xué)習(xí)！