edgeR中三種標準化方法TMM\UQ\RLE的比較

jingda 2022-06-01 發(fā)布于江蘇

展開全文

關于RNA-seq中的reads count標準化處理的方法匯總，請先看看這篇：
當我們在說RNA-seq reads count標準化時，其實在說什么？

本文集中討論常用的edgeR包中三種標準化方法TMM\UQ\RLE的比較

英文原貼Normalisation methods implemented in edgeR

1.首先創(chuàng)建一個數據集

包含四個樣品c1,c2是正常組，p1,p2是病人。共有50個轉錄本，每個樣品內轉錄本counts的總數都是500個，前25個轉錄本在四個樣品里都有表達，其中病人轉錄本的數目(20)是對照組(10)的兩倍, 后25個轉錄本只在正常組中檢測到。

#prepare example
control_1 <- rep(10, 50)
control_2 <- rep(10, 50)
patient_1 <- c(rep(20, 25),rep(0,25))
patient_2 <- c(rep(20, 25),rep(0,25))
 
df <- data.frame(c1=control_1,
                 c2=control_2,
                 p1=patient_1,
                 p2=patient_2)
 
head(df)
  c1 c2 p1 p2
1 10 10 20 20
2 10 10 20 20
3 10 10 20 20
4 10 10 20 20
5 10 10 20 20
6 10 10 20 20
 
tail(df)
   c1 c2 p1 p2
45 10 10  0  0
46 10 10  0  0
47 10 10  0  0
48 10 10  0  0
49 10 10  0  0
50 10 10  0  0
 
#equal depth
colSums(df)
 c1  c2  p1  p2 
500 500 500 500

數據集信息詳見Robinson and Oshlack http:///2010/11/3/R25

2.如果不做標準化處理

#load library
library(edgeR)
 
#create group vector
group <- c('control','control','patient','patient')
 
#create DGEList object
d <- DGEList(counts=df, group=group)
 
#check out the DGEList object
d
An object of class "DGEList"
$counts
  c1 c2 p1 p2
1 10 10 20 20
2 10 10 20 20
3 10 10 20 20
4 10 10 20 20
5 10 10 20 20
45 more rows ...
 
$samples
     group lib.size norm.factors
c1 control      500            1
c2 control      500            1
p1 patient      500            1
p2 patient      500            1
 
d <- DGEList(counts=df, group=group)
d <- estimateCommonDisp(d)
 
#perform the DE test
de <- exactTest(d)
 
#how many differentially expressed transcripts?
table(p.adjust(de$table$PValue, method="BH")<0.05)
 
TRUE
  50

可以看到：檢測出共50個轉錄本有差異，即每個轉錄本都是差異表達的，假陽性很高。

3.TMM normalisation

TMM <- calcNormFactors(d, method="TMM")
TMM
An object of class "DGEList"
$counts
  c1 c2 p1 p2
1 10 10 20 20
2 10 10 20 20
3 10 10 20 20
4 10 10 20 20
5 10 10 20 20
45 more rows ...
 
$samples
     group lib.size norm.factors
c1 control      500    0.7071068
c2 control      500    0.7071068
p1 patient      500    1.4142136
p2 patient      500    1.4142136

我們看到對前25個轉錄本而言，正常組和病人之間沒有差異 (10/0.7071068 (~14.14) 等于 20/1.4142136 (~14.14))。因此檢測出有25個轉錄本存在差異（后25個轉錄本）

TMM <- estimateCommonDisp(TMM)
TMM <- exactTest(TMM)
table(p.adjust(TMM$table$PValue, method="BH")<0.05)
 
FALSE  TRUE
   25    25

4.RLE normalisation

RLE
An object of class "DGEList"
$counts
  c1 c2 p1 p2
1 10 10 20 20
2 10 10 20 20
3 10 10 20 20
4 10 10 20 20
5 10 10 20 20
45 more rows ...
 
$samples
     group lib.size norm.factors
c1 control      500    0.7071068
c2 control      500    0.7071068
p1 patient      500    1.4142136
p2 patient      500    1.4142136
RLE <- estimateCommonDisp(RLE)
RLE <- exactTest(RLE)
table(p.adjust(RLE$table$PValue, method="BH")<0.05)
 
FALSE  TRUE
   25    25

5.UQ normalisation

uq
An object of class "DGEList"
$counts
  c1 c2 p1 p2
1 10 10 20 20
2 10 10 20 20
3 10 10 20 20
4 10 10 20 20
5 10 10 20 20
45 more rows ...
 
$samples
     group lib.size norm.factors
c1 control      500    0.7071068
c2 control      500    0.7071068
p1 patient      500    1.4142136
p2 patient      500    1.4142136
 
uq <- estimateCommonDisp(uq)
uq <- exactTest(uq)
table(p.adjust(uq$table$PValue, method="BH")<0.05)
 
FALSE  TRUE
   25    25

因為數據比較簡單，這里三種標準化方法得到的結果一致，那么真實測序數據的情況又如何呢？

6.測試一套真實數據

my_url <-"[https:///file/pnas_expression.txt](https:///file/pnas_expression.txt)"

data <-read.table(my_url, header=TRUE, sep="\t")

dim(data)

[1] 37435     9

ensembl_ID lane1 lane2 lane3 lane4 lane5 lane6 lane8  len

1 ENSG00000215696     0     0     0     0     0     0     0  330

2 ENSG00000215700     0     0     0     0     0     0     0 2370

3 ENSG00000215699     0     0     0     0     0     0     0 1842

4 ENSG00000215784     0     0     0     0     0     0     0 2393

5 ENSG00000212914     0     0     0     0     0     0     0  384

6 ENSG00000212042     0     0     0     0     0     0     0   92

準備DGEList

rownames(d) <- data[,1]
group <- c(rep("Control",4),rep("DHT",3))
d <- DGEList(counts = d, group=group)
 
An object of class "DGEList"
$counts
                lane1 lane2 lane3 lane4 lane5 lane6 lane8
ENSG00000215696     0     0     0     0     0     0     0
ENSG00000215700     0     0     0     0     0     0     0
ENSG00000215699     0     0     0     0     0     0     0
ENSG00000215784     0     0     0     0     0     0     0
ENSG00000212914     0     0     0     0     0     0     0
37430 more rows ...
 
$samples
        group lib.size norm.factors
lane1 Control   978576            1
lane2 Control  1156844            1
lane3 Control  1442169            1
lane4 Control  1485604            1
lane5     DHT  1823460            1
lane6     DHT  1834335            1
lane8     DHT   681743            1

還是先不做標準化處理

no_norm <- exactTest(no_norm)
table(p.adjust(no_norm$table$PValue, method="BH")<0.05)
 
FALSE  TRUE
33404  4031

TMM normalisation

TMM <- calcNormFactors(d, method="TMM")
TMM
An object of class "DGEList"
$counts
                lane1 lane2 lane3 lane4 lane5 lane6 lane8
ENSG00000215696     0     0     0     0     0     0     0
ENSG00000215700     0     0     0     0     0     0     0
ENSG00000215699     0     0     0     0     0     0     0
ENSG00000215784     0     0     0     0     0     0     0
ENSG00000212914     0     0     0     0     0     0     0
37430 more rows ...
 
$samples
        group lib.size norm.factors
lane1 Control   978576    1.0350786
lane2 Control  1156844    1.0379515
lane3 Control  1442169    1.0287815
lane4 Control  1485604    1.0222095
lane5     DHT  1823460    0.9446243
lane6     DHT  1834335    0.9412769
lane8     DHT   681743    0.9954283
 
TMM <- estimateCommonDisp(TMM)
TMM <- exactTest(TMM)
table(p.adjust(TMM$table$PValue, method="BH")<0.05)
 
FALSE  TRUE
33519  3916

RLE

RLE <- calcNormFactors(d, method="RLE")
RLE
An object of class "DGEList"
$counts
                lane1 lane2 lane3 lane4 lane5 lane6 lane8
ENSG00000215696     0     0     0     0     0     0     0
ENSG00000215700     0     0     0     0     0     0     0
ENSG00000215699     0     0     0     0     0     0     0
ENSG00000215784     0     0     0     0     0     0     0
ENSG00000212914     0     0     0     0     0     0     0
37430 more rows ...
 
$samples
        group lib.size norm.factors
lane1 Control   978576    1.0150010
lane2 Control  1156844    1.0236675
lane3 Control  1442169    1.0345426
lane4 Control  1485604    1.0399724
lane5     DHT  1823460    0.9706692
lane6     DHT  1834335    0.9734955
lane8     DHT   681743    0.9466713
 
RLE <- estimateCommonDisp(RLE)
RLE <- exactTest(RLE)
table(p.adjust(RLE$table$PValue, method="BH")<0.05)
 
FALSE  TRUE
33465  3970

the upper quartile method

uq <- calcNormFactors(d, method="upperquartile")
uq
An object of class "DGEList"
$counts
                lane1 lane2 lane3 lane4 lane5 lane6 lane8
ENSG00000215696     0     0     0     0     0     0     0
ENSG00000215700     0     0     0     0     0     0     0
ENSG00000215699     0     0     0     0     0     0     0
ENSG00000215784     0     0     0     0     0     0     0
ENSG00000212914     0     0     0     0     0     0     0
37430 more rows ...
 
$samples
        group lib.size norm.factors
lane1 Control   978576    1.0272514
lane2 Control  1156844    1.0222982
lane3 Control  1442169    1.0250528
lane4 Control  1485604    1.0348864
lane5     DHT  1823460    0.9728534
lane6     DHT  1834335    0.9670858
lane8     DHT   681743    0.9541011
 
uq <- estimateCommonDisp(uq)
uq <- exactTest(uq)
table(p.adjust(uq$table$PValue, method="BH")<0.05)
 
FALSE  TRUE
33466  3969

以上四種處理方法找到的差異基因取交集，可以看出不做標準化處理會得到405個假陽性和342個假陰性的轉錄本

library(gplots)
 
get_de <- function(x, pvalue){
  my_i <- p.adjust(x$PValue, method="BH") < pvalue
  row.names(x)[my_i]
}
 
my_de_no_norm <- get_de(no_norm$table, 0.05)
my_de_tmm <- get_de(TMM$table, 0.05)
my_de_rle <- get_de(RLE$table, 0.05)
my_de_uq <- get_de(uq$table, 0.05)
 
gplots::venn(list(no_norm = my_de_no_norm, TMM = my_de_tmm, RLE = my_de_rle, UQ = my_de_uq))

不做標準化會得到405個假陽性和342個假陰性的轉錄本

三種標準化方法找到的差異基因大部分是一致的

gplots::venn(list(TMM = my_de_tmm, RLE = my_de_rle, UQ = my_de_uq))

1.png

小結

三種標準化方法效果類似，處理結果都比不做標準化要好
The normalisation factors were quite similar between all normalisation methods, which is why the results of the differential expression were quite concordant. Most methods down sized the DHT samples with a normalisation factor of less than one to account for the larger library sizes of these samples.