edgeR中三种标准化方法TMM\UQ\RLE的比较 edgeR中三种标准化方法TMM\UQ\R

关于RNA-seq中的reads count标准化处理的方法汇总，请先看看这篇：
当我们在说RNA-seq reads count标准化时，其实在说什么？
本文集中讨论常用的edgeR包中三种标准化方法TMM\UQ\RLE的比较
【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 1000 46 10 1000 47 10 1000 48 10 1000 49 10 1000 50 10 1000 #equal depth colSums(df) c1c2p1p2 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 control5001 c2 control5001 p1 patient5001 p2 patient5001 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 control5000.7071068 c2 control5000.7071068 p1 patient5001.4142136 p2 patient5001.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) FALSETRUE 2525

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 control5000.7071068 c2 control5000.7071068 p1 patient5001.4142136 p2 patient5001.4142136 RLE <- estimateCommonDisp(RLE) RLE <- exactTest(RLE) table(p.adjust(RLE$table$PValue, method="BH")<0.05) FALSETRUE 2525

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 control5000.7071068 c2 control5000.7071068 p1 patient5001.4142136 p2 patient5001.4142136 uq <- estimateCommonDisp(uq) uq <- exactTest(uq) table(p.adjust(uq$table$PValue, method="BH")<0.05) FALSETRUE 2525

因为数据比较简单，这里三种标准化方法得到的结果一致，那么真实测序数据的情况又如何呢？
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] 374359ensembl_ID lane1 lane2 lane3 lane4 lane5 lane6 lane8len1 ENSG0000021569600000003302 ENSG000002157000000000 23703 ENSG000002156990000000 18424 ENSG000002157840000000 23935 ENSG0000021291400000003846 ENSG00000212042000000092

准备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 ENSG000002156960000000 ENSG000002157000000000 ENSG000002156990000000 ENSG000002157840000000 ENSG000002129140000000 37430 more rows ... $samples group lib.size norm.factors lane1 Control9785761 lane2 Control11568441 lane3 Control14421691 lane4 Control14856041 lane5DHT18234601 lane6DHT18343351 lane8DHT6817431

还是先不做标准化处理

no_norm <- exactTest(no_norm) table(p.adjust(no_norm$table$PValue, method="BH")<0.05) FALSETRUE 334044031

TMM normalisation

TMM <- calcNormFactors(d, method="TMM") TMM An object of class "DGEList" $counts lane1 lane2 lane3 lane4 lane5 lane6 lane8 ENSG000002156960000000 ENSG000002157000000000 ENSG000002156990000000 ENSG000002157840000000 ENSG000002129140000000 37430 more rows ... $samples group lib.size norm.factors lane1 Control9785761.0350786 lane2 Control11568441.0379515 lane3 Control14421691.0287815 lane4 Control14856041.0222095 lane5DHT18234600.9446243 lane6DHT18343350.9412769 lane8DHT6817430.9954283 TMM <- estimateCommonDisp(TMM) TMM <- exactTest(TMM) table(p.adjust(TMM$table$PValue, method="BH")<0.05) FALSETRUE 335193916

RLE

RLE <- calcNormFactors(d, method="RLE") RLE An object of class "DGEList" $counts lane1 lane2 lane3 lane4 lane5 lane6 lane8 ENSG000002156960000000 ENSG000002157000000000 ENSG000002156990000000 ENSG000002157840000000 ENSG000002129140000000 37430 more rows ... $samples group lib.size norm.factors lane1 Control9785761.0150010 lane2 Control11568441.0236675 lane3 Control14421691.0345426 lane4 Control14856041.0399724 lane5DHT18234600.9706692 lane6DHT18343350.9734955 lane8DHT6817430.9466713 RLE <- estimateCommonDisp(RLE) RLE <- exactTest(RLE) table(p.adjust(RLE$table$PValue, method="BH")<0.05) FALSETRUE 334653970

the upper quartile method

uq <- calcNormFactors(d, method="upperquartile") uq An object of class "DGEList" $counts lane1 lane2 lane3 lane4 lane5 lane6 lane8 ENSG000002156960000000 ENSG000002157000000000 ENSG000002156990000000 ENSG000002157840000000 ENSG000002129140000000 37430 more rows ... $samples group lib.size norm.factors lane1 Control9785761.0272514 lane2 Control11568441.0222982 lane3 Control14421691.0250528 lane4 Control14856041.0348864 lane5DHT18234600.9728534 lane6DHT18343350.9670858 lane8DHT6817430.9541011 uq <- estimateCommonDisp(uq) uq <- exactTest(uq) table(p.adjust(uq$table$PValue, method="BH")<0.05) FALSETRUE 334663969

以上四种处理方法找到的差异基因取交集，可以看出不做标准化处理会得到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))

$edgeR中三种标准化方法TMM\UQ\RLE的比较$
文章图片
不做标准化会得到405个假阳性和342个假阴性的转录本
三种标准化方法找到的差异基因大部分是一致的

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

$edgeR中三种标准化方法TMM\UQ\RLE的比较$
文章图片
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.