Compare edgeR and DESeq2 RNA-seq results
Mark Ziemann
2020-10-06
library("edgeR")
library("DESeq2")
library("eulerr")x <- read.table("mrna_fulllen_pe_strrev_fmt.mx",row.names=1,header=TRUE)
xx <- x[which(rowMeans(x)>10),]
ss <- as.data.frame(colnames(xx))
ss$trt <- as.integer(grepl("Yap",ss[,1]))
rownames(ss) <- ss[,1]
ss[,1]=NULLdesign <- model.matrix(~ss$trt)
rownames(design) <- rownames(ss)
z <- DGEList(counts=xx)
z <- calcNormFactors(z)
z <- estimateDisp(z, design,robust=TRUE,prior.df=1)
fit <- glmFit(z, design)
lrt<-glmLRT(fit)
dge<-as.data.frame(topTags(lrt,n=Inf))
dge$dispersion<-lrt$dispersion
dge<-merge(dge,lrt$fitted.values,by='row.names')
rownames(dge)=dge$Row.names
dge$Row.names=NULL
dge<-dge[order(dge$PValue),]
head(dge,20)## logFC logCPM LR PValue
## ENSMUSG00000101968_1700027A15Rik 2.2518458 3.1219122 552.4925 3.612208e-122
## ENSMUSG00000008318_Relt 1.4374097 4.3983287 382.8638 2.961234e-85
## ENSMUSG00000006205_Htra1 -2.6002518 5.3348815 318.0284 3.893790e-71
## ENSMUSG00000013846_St3gal1 1.2370422 4.9822309 281.9756 2.786634e-63
## ENSMUSG00000030772_Dkk3 -2.2602457 2.5371639 277.5167 2.610729e-62
## ENSMUSG00000038037_Socs1 4.0317922 0.7220168 275.2180 8.273848e-62
## ENSMUSG00000026177_Slc11a1 3.8329527 0.7731070 254.9755 2.136881e-57
## ENSMUSG00000022139_Mbnl2 -0.9855745 9.8528171 242.9677 8.863872e-55
## ENSMUSG00000091514_Gm17484 -0.9339419 4.3276244 240.7327 2.722345e-54
## ENSMUSG00000090556_Olfr753-ps1 5.4676434 0.2754831 238.8882 6.872759e-54
## ENSMUSG00000026764_Kif5c 2.8727232 1.0961751 231.8811 2.317943e-52
## ENSMUSG00000053835_H2-T24 4.0337831 0.5639515 229.3520 8.253972e-52
## ENSMUSG00000085289_Gm15337 -2.3779408 3.8132843 216.6660 4.826468e-49
## ENSMUSG00000002565_Scin -1.4626953 5.2032771 212.0188 4.982374e-48
## ENSMUSG00000029826_Zc3hav1 1.1448050 3.6797174 204.1650 2.576159e-46
## ENSMUSG00000039621_Prex1 0.8815089 6.8051528 202.9579 4.724815e-46
## ENSMUSG00000072294_Klf12 -0.7125661 9.1999949 197.0681 9.112879e-45
## ENSMUSG00000020638_Cmpk2 2.3349500 0.4473717 193.3890 5.789234e-44
## ENSMUSG00000020641_Rsad2 4.8280813 0.4463170 191.9497 1.193359e-43
## ENSMUSG00000031394_Opn1mw 3.8497411 -0.6649544 191.4646 1.522796e-43
## FDR dispersion Lacz.1
## ENSMUSG00000101968_1700027A15Rik 5.922576e-118 0.235711384 196.923663
## ENSMUSG00000008318_Relt 2.427619e-81 0.017113170 744.224460
## ENSMUSG00000006205_Htra1 2.128086e-67 0.031452438 4542.614908
## ENSMUSG00000013846_St3gal1 1.142241e-59 0.007448975 1233.089436
## ENSMUSG00000030772_Dkk3 8.561102e-59 0.023213030 627.313926
## ENSMUSG00000038037_Socs1 2.260967e-58 0.008671645 12.086177
## ENSMUSG00000026177_Slc11a1 5.005185e-54 0.054470606 14.289205
## ENSMUSG00000022139_Mbnl2 1.816651e-51 0.002687166 80586.793835
## ENSMUSG00000091514_Gm17484 4.959508e-51 0.003361192 1726.851152
## ENSMUSG00000090556_Olfr753-ps1 1.126858e-50 0.012437114 3.271548
## ENSMUSG00000026764_Kif5c 3.454999e-49 0.157379901 32.977305
## ENSMUSG00000053835_H2-T24 1.127768e-48 0.004784200 10.788396
## ENSMUSG00000085289_Gm15337 6.087290e-46 0.008588910 1544.261063
## ENSMUSG00000002565_Scin 5.835071e-45 0.121419745 3544.574952
## ENSMUSG00000029826_Zc3hav1 2.815913e-43 0.011081432 522.086178
## ENSMUSG00000039621_Prex1 4.841754e-43 0.009786446 5158.935995
## ENSMUSG00000072294_Klf12 8.789103e-42 0.013870682 47907.418729
## ENSMUSG00000020638_Cmpk2 5.273349e-41 0.043563875 28.678638
## ENSMUSG00000020641_Rsad2 1.029806e-40 0.002165754 5.814343
## ENSMUSG00000031394_Opn1mw 1.248388e-40 0.015826155 4.946637
## Lacz.2 Lacz.3 Yap.1
## ENSMUSG00000101968_1700027A15Rik 151.029891 2.498441e+02 1145.89879
## ENSMUSG00000008318_Relt 570.780257 9.442242e+02 2461.58517
## ENSMUSG00000006205_Htra1 3483.942070 5.763378e+03 914.65185
## ENSMUSG00000013846_St3gal1 945.713482 1.564465e+03 3549.51269
## ENSMUSG00000030772_Dkk3 481.116146 7.958956e+02 159.78111
## ENSMUSG00000038037_Socs1 9.269449 1.533416e+01 243.66637
## ENSMUSG00000026177_Slc11a1 10.959054 1.812922e+01 250.60369
## ENSMUSG00000022139_Mbnl2 61805.749993 1.022433e+05 49697.93801
## ENSMUSG00000091514_Gm17484 1324.402243 2.190918e+03 1103.68495
## ENSMUSG00000090556_Olfr753-ps1 2.509102 4.150730e+00 183.16801
## ENSMUSG00000026764_Kif5c 25.291825 4.183948e+01 295.89352
## ENSMUSG00000053835_H2-T24 8.274121 1.368762e+01 218.04848
## ENSMUSG00000085289_Gm15337 1184.365435 1.959259e+03 362.66701
## ENSMUSG00000002565_Scin 2718.498936 4.497129e+03 1570.31031
## ENSMUSG00000029826_Zc3hav1 400.412105 6.623894e+02 1409.87538
## ENSMUSG00000039621_Prex1 3956.627298 6.545326e+03 11606.14769
## ENSMUSG00000072294_Klf12 36742.421479 6.078185e+04 35699.40021
## ENSMUSG00000020638_Cmpk2 21.994978 3.638561e+01 177.29108
## ENSMUSG00000020641_Rsad2 4.459290 7.376865e+00 205.68301
## ENSMUSG00000031394_Opn1mw 3.793805 6.275975e+00 89.08388
## Yap.2 Yap.3
## ENSMUSG00000101968_1700027A15Rik 792.08918 998.27883
## ENSMUSG00000008318_Relt 1701.54206 2144.47244
## ENSMUSG00000006205_Htra1 632.24243 796.82219
## ENSMUSG00000013846_St3gal1 2453.55928 3092.24813
## ENSMUSG00000030772_Dkk3 110.44683 139.19737
## ENSMUSG00000038037_Socs1 168.43154 212.27615
## ENSMUSG00000026177_Slc11a1 173.22688 218.31977
## ENSMUSG00000022139_Mbnl2 34353.12048 43295.62096
## ENSMUSG00000091514_Gm17484 762.90936 961.50318
## ENSMUSG00000090556_Olfr753-ps1 126.61275 159.57147
## ENSMUSG00000026764_Kif5c 204.53295 257.77515
## ENSMUSG00000053835_H2-T24 150.72347 189.95847
## ENSMUSG00000085289_Gm15337 250.68934 315.94657
## ENSMUSG00000002565_Scin 1085.45871 1368.01571
## ENSMUSG00000029826_Zc3hav1 974.55993 1228.24875
## ENSMUSG00000039621_Prex1 8022.61433 10110.99034
## ENSMUSG00000072294_Klf12 24676.79436 31100.43921
## ENSMUSG00000020638_Cmpk2 122.55039 154.45163
## ENSMUSG00000020641_Rsad2 142.17599 179.18598
## ENSMUSG00000031394_Opn1mw 61.57819 77.60768dge_edger <- dge
sig <- subset(dge_edger,FDR<0.05)
dge_edger_up <- rownames(subset(sig,logFC>0))
dge_edger_dn <- rownames(subset(sig,logFC<0))
length(dge_edger_up)## [1] 4279length(dge_edger_dn)## [1] 3755dds <- DESeqDataSetFromMatrix(countData = xx , colData = ss, design = ~ trt )## the design formula contains one or more numeric variables with integer values,
## specifying a model with increasing fold change for higher values.
## did you mean for this to be a factor? if so, first convert
## this variable to a factor using the factor() functionres <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingz<- results(res)
vsd <- vst(dds, blind=FALSE)
zz<-cbind(as.data.frame(z),assay(vsd))
dge<-as.data.frame(zz[order(zz$pvalue),])
head(dge,20)## baseMean log2FoldChange lfcSE stat
## ENSMUSG00000006205_Htra1 2674.7133 -2.6128492 0.12719949 -20.54135
## ENSMUSG00000101968_1700027A15Rik 571.3900 2.2376857 0.12127751 18.45095
## ENSMUSG00000022139_Mbnl2 61210.4162 -0.9984953 0.05958637 -16.75711
## ENSMUSG00000085289_Gm15337 930.8163 -2.3923528 0.15067219 -15.87787
## ENSMUSG00000008318_Relt 1389.1309 1.4242928 0.09041223 15.75332
## ENSMUSG00000018899_Irf1 793.2648 3.7441386 0.24688719 15.16538
## ENSMUSG00000024411_Aqp4 1286.7561 -2.1522544 0.14271365 -15.08093
## ENSMUSG00000030772_Dkk3 384.1212 -2.2735553 0.15353468 -14.80809
## ENSMUSG00000070695_Cntnap5a 741.8563 2.6803854 0.18202195 14.72562
## ENSMUSG00000034422_Parp14 1485.2801 3.5878730 0.24708537 14.52078
## ENSMUSG00000072294_Klf12 38921.1350 -0.7254631 0.05032542 -14.41544
## ENSMUSG00000002565_Scin 2439.6831 -1.4756926 0.10331695 -14.28316
## ENSMUSG00000027639_Samhd1 3339.7782 1.9389713 0.13591504 14.26605
## ENSMUSG00000013846_St3gal1 2084.0437 1.2239041 0.08730137 14.01930
## ENSMUSG00000041669_Prima1 1511.2296 -3.2894746 0.23947622 -13.73612
## ENSMUSG00000026104_Stat1 1148.5142 2.6748918 0.19483593 13.72895
## ENSMUSG00000021536_Adcy2 45823.7077 -0.9635923 0.07253328 -13.28483
## ENSMUSG00000024431_Nr3c1 24530.4980 -1.0679329 0.08087676 -13.20445
## ENSMUSG00000027366_Sppl2a 7193.0084 -0.8222866 0.06264149 -13.12687
## ENSMUSG00000039621_Prex1 7381.2387 0.8682319 0.06727827 12.90509
## pvalue padj Lacz.1 Lacz.2
## ENSMUSG00000006205_Htra1 9.197219e-94 1.507976e-89 12.060288 12.409609
## ENSMUSG00000101968_1700027A15Rik 5.123307e-76 4.200087e-72 9.101363 9.086642
## ENSMUSG00000022139_Mbnl2 5.024854e-63 2.746250e-59 16.267350 16.371911
## ENSMUSG00000085289_Gm15337 9.019623e-57 3.697143e-53 11.088197 10.812993
## ENSMUSG00000008318_Relt 6.516825e-56 2.136997e-52 10.200818 10.150501
## ENSMUSG00000018899_Irf1 5.994850e-52 1.638193e-48 8.632941 8.582452
## ENSMUSG00000024411_Aqp4 2.162062e-51 5.064167e-48 11.256371 11.167554
## ENSMUSG00000030772_Dkk3 1.298790e-49 2.661869e-46 9.917739 10.066069
## ENSMUSG00000070695_Cntnap5a 4.414258e-49 8.041798e-46 9.122993 9.134956
## ENSMUSG00000034422_Parp14 8.948261e-48 1.467157e-44 8.914011 9.094845
## ENSMUSG00000072294_Klf12 4.137956e-47 6.167811e-44 15.623479 15.547920
## ENSMUSG00000002565_Scin 2.786710e-46 3.807575e-43 11.988836 11.985923
## ENSMUSG00000027639_Samhd1 3.561702e-46 4.492129e-43 10.788720 10.727445
## ENSMUSG00000013846_St3gal1 1.187799e-44 1.391083e-41 10.756327 10.683221
## ENSMUSG00000041669_Prima1 6.169285e-43 6.743440e-40 11.006130 11.810044
## ENSMUSG00000026104_Stat1 6.811816e-43 6.980408e-40 9.349427 9.309648
## ENSMUSG00000021536_Adcy2 2.835051e-40 2.734323e-37 15.883009 15.896423
## ENSMUSG00000024431_Nr3c1 8.270505e-40 7.533511e-37 14.935861 15.049831
## ENSMUSG00000027366_Sppl2a 2.310126e-39 1.993517e-36 13.176559 13.270175
## ENSMUSG00000039621_Prex1 4.213442e-38 3.443704e-35 12.505391 12.487742
## Lacz.3 Yap.1 Yap.2 Yap.3
## ENSMUSG00000006205_Htra1 12.371211 10.047044 10.174809 10.193805
## ENSMUSG00000101968_1700027A15Rik 8.997703 10.360850 10.394358 10.357106
## ENSMUSG00000022139_Mbnl2 16.328852 15.311146 15.298335 15.385093
## ENSMUSG00000085289_Gm15337 10.885140 9.224973 9.337152 9.435957
## ENSMUSG00000008318_Relt 10.077178 11.268265 11.267546 11.185285
## ENSMUSG00000018899_Irf1 8.908633 10.605982 10.916354 11.046206
## ENSMUSG00000024411_Aqp4 11.424344 9.841911 9.684678 9.579021
## ENSMUSG00000030772_Dkk3 9.943792 8.754052 8.860593 8.835800
## ENSMUSG00000070695_Cntnap5a 8.923080 10.671070 10.534518 10.891233
## ENSMUSG00000034422_Parp14 9.374793 11.254149 11.728181 11.812121
## ENSMUSG00000072294_Klf12 15.561834 14.862865 14.848040 14.870729
## ENSMUSG00000002565_Scin 11.914426 10.853825 10.667225 10.605372
## ENSMUSG00000027639_Samhd1 10.852469 12.192260 12.549356 12.653327
## ENSMUSG00000013846_St3gal1 10.589862 11.712096 11.728693 11.648083
## ENSMUSG00000041669_Prima1 11.883628 9.344649 9.214597 9.307406
## ENSMUSG00000026104_Stat1 9.435907 10.904191 11.138865 11.538983
## ENSMUSG00000021536_Adcy2 15.909832 15.034882 14.826637 14.956856
## ENSMUSG00000024431_Nr3c1 15.120664 13.933450 13.942154 14.085416
## ENSMUSG00000027366_Sppl2a 13.237871 12.481383 12.472160 12.401433
## ENSMUSG00000039621_Prex1 12.383215 13.322225 13.240071 13.276269dge_deseq2 <- dge
sig <- subset(dge,padj<0.05)
dge_deseq2_up <- rownames(subset(sig,log2FoldChange>0))
dge_deseq2_dn <- rownames(subset(sig,log2FoldChange<0))
length(dge_deseq2_up)## [1] 3593length(dge_deseq2_dn)## [1] 3314v1 <- list("edgeR up"=dge_edger_up, "edgeR dn"=dge_edger_dn,
"DESeq2 up"=dge_deseq2_up,"DESeq2 dn"=dge_deseq2_dn)
plot(euler(v1),quantities = TRUE)
Let's look at some genes of interest.
# YAP1
# edgeR
t(dge_edger[grep("_Yap1",rownames(dge_edger)),])## ENSMUSG00000053110_Yap1
## logFC 1.638272e-01
## logCPM 7.097791e+00
## LR 1.234145e+00
## PValue 2.666025e-01
## FDR 3.670205e-01
## dispersion 4.500414e-02
## Lacz.1 8.472007e+03
## Lacz.2 6.497575e+03
## Lacz.3 1.074874e+04
## Yap.1 1.158955e+04
## Yap.2 8.011140e+03
## Yap.3 1.009653e+04# DESeq2
t(dge_deseq2[grep("_Yap1",rownames(dge_deseq2)),])## ENSMUSG00000053110_Yap1
## baseMean 9055.4131973
## log2FoldChange 0.1513948
## lfcSE 0.1341132
## stat 1.1288581
## pvalue 0.2589577
## padj 0.3858077
## Lacz.1 13.2192928
## Lacz.2 13.1443207
## Lacz.3 13.0330801
## Yap.1 13.5112230
## Yap.2 13.2750784
## Yap.3 13.0067550# Amotl2
# edgeR
t(dge_edger[grep("Amotl2",rownames(dge_edger)),])## ENSMUSG00000032531_Amotl2
## logFC -1.270718e+00
## logCPM 3.549146e+00
## LR 4.158470e+01
## PValue 1.128733e-10
## FDR 1.659795e-09
## dispersion 7.116365e-03
## Lacz.1 1.083376e+03
## Lacz.2 8.308910e+02
## Lacz.3 1.374517e+03
## Yap.1 5.482116e+02
## Yap.2 3.789449e+02
## Yap.3 4.775884e+02# DESeq2
t(dge_deseq2[grep("Amotl2",rownames(dge_deseq2)),])## ENSMUSG00000032531_Amotl2
## baseMean 7.746610e+02
## log2FoldChange -1.282856e+00
## lfcSE 1.812664e-01
## stat -7.077190e+00
## pvalue 1.471073e-12
## padj 3.605338e-11
## Lacz.1 1.056040e+01
## Lacz.2 1.054015e+01
## Lacz.3 1.050053e+01
## Yap.1 9.495169e+00
## Yap.2 9.904959e+00
## Yap.3 9.567554e+00# Glul
# edgeR
t(dge_edger[grep("Glul",rownames(dge_edger)),])## ENSMUSG00000026473_Glul
## logFC -1.072362e+00
## logCPM 7.552828e+00
## LR 6.614173e+00
## PValue 1.011703e-02
## FDR 2.356234e-02
## dispersion 4.563751e-02
## Lacz.1 1.668986e+04
## Lacz.2 1.280023e+04
## Lacz.3 2.117503e+04
## Yap.1 9.691673e+03
## Yap.2 6.699256e+03
## Yap.3 8.443147e+03# DESeq2
t(dge_deseq2[grep("Glul",rownames(dge_deseq2)),])## ENSMUSG00000026473_Glul
## baseMean 1.242085e+04
## log2FoldChange -1.087149e+00
## lfcSE 4.329426e-01
## stat -2.511069e+00
## pvalue 1.203660e-02
## padj 3.108887e-02
## Lacz.1 1.396869e+01
## Lacz.2 1.409649e+01
## Lacz.3 1.416425e+01
## Yap.1 1.209209e+01
## Yap.2 1.299335e+01
## Yap.3 1.361613e+01# Dkk3
# edgeR
t(dge_edger[grep("Dkk3",rownames(dge_edger)),])## ENSMUSG00000030772_Dkk3
## logFC -2.260246e+00
## logCPM 2.537164e+00
## LR 2.775167e+02
## PValue 2.610729e-62
## FDR 8.561102e-59
## dispersion 2.321303e-02
## Lacz.1 6.273139e+02
## Lacz.2 4.811161e+02
## Lacz.3 7.958956e+02
## Yap.1 1.597811e+02
## Yap.2 1.104468e+02
## Yap.3 1.391974e+02# DESeq2
t(dge_deseq2[grep("Dkk3",rownames(dge_deseq2)),])## ENSMUSG00000030772_Dkk3
## baseMean 3.841212e+02
## log2FoldChange -2.273555e+00
## lfcSE 1.535347e-01
## stat -1.480809e+01
## pvalue 1.298790e-49
## padj 2.661869e-46
## Lacz.1 9.917739e+00
## Lacz.2 1.006607e+01
## Lacz.3 9.943792e+00
## Yap.1 8.754052e+00
## Yap.2 8.860593e+00
## Yap.3 8.835800e+00