关于RNA-seq中的reads count标准化处理的方法汇总，请先看看这篇：
当我们在说RNA-seq reads count标准化时，其实在说什么？
本文集中讨论常用的edgeR包中三种标准化方法TMM\UQ\RLE的比较
英文原贴Normalisation methods implemented in edgeR
1.首先创建一个数据集
包含四个样品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://genomebiology.com/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://davetang.org/file/pnas_expression.txt](https://davetang.org/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.
