Movatterモバイル変換

################################################################################################################################################################## Running command:######   /Library/Frameworks/R.framework/Resources/bin/R CMD check --install=check:edgeR.install-out.txt --library=/Library/Frameworks/R.framework/Resources/library --no-vignettes --timings edgeR_4.8.1.tar.gz###############################################################################################################################################################* using log directory ‘/Users/biocbuild/bbs-3.22-bioc/meat/edgeR.Rcheck’* using R version 4.5.2 Patched (2025-11-04 r88984)* using platform: aarch64-apple-darwin20* R was compiled by    Apple clang version 16.0.0 (clang-1600.0.26.6)    GNU Fortran (GCC) 14.2.0* running under: macOS Ventura 13.7.8* using session charset: UTF-8* using option ‘--no-vignettes’* checking for file ‘edgeR/DESCRIPTION’ ... OK* this is package ‘edgeR’ version ‘4.8.1’* checking package namespace information ... OK* checking package dependencies ... OK* checking if this is a source package ... OK* checking if there is a namespace ... OK* checking for hidden files and directories ... OK* checking for portable file names ... OK* checking for sufficient/correct file permissions ... OK* checking whether package ‘edgeR’ can be installed ... OK* used C compiler: ‘Apple clang version 15.0.0 (clang-1500.0.40.1)’* used SDK: ‘MacOSX11.3.1.sdk’* checking installed package size ... OK* checking package directory ... OK* checking ‘build’ directory ... OK* checking DESCRIPTION meta-information ... OK* checking top-level files ... OK* checking for left-over files ... OK* checking index information ... OK* checking package subdirectories ... OK* checking code files for non-ASCII characters ... OK* checking R files for syntax errors ... OK* checking whether the package can be loaded ... OK* checking whether the package can be loaded with stated dependencies ... OK* checking whether the package can be unloaded cleanly ... OK* checking whether the namespace can be loaded with stated dependencies ... OK* checking whether the namespace can be unloaded cleanly ... OK* checking whether startup messages can be suppressed ... OK* checking dependencies in R code ... OK* checking S3 generic/method consistency ... OK* checking replacement functions ... OK* checking foreign function calls ... OK* checking R code for possible problems ... OK* checking Rd files ... OK* checking Rd metadata ... OK* checking Rd cross-references ... NOTENon-topic package-anchored link(s) in Rd file 'asmatrix.Rd':  ‘[limma:asmatrix]{as.matrix}’See section 'Cross-references' in the 'Writing R Extensions' manual.Found the following Rd file(s) with Rd \link{} targets missing packageanchors:  camera.Rd: ids2indices, camera  decidetestsDGE.Rd: TestResults-class, decideTests  diffSplice.Rd: diffSplice, topSplice, plotSplice  dim.Rd: 02.Classes  glmQLFit.Rd: squeezeVar  glmTreat.Rd: treat  goana.Rd: goana.default, kegga.default, goana, topGO, kegga, topKEGG  head.Rd: head.EList  normalizeBetweenArraysDGEList.Rd: normalizeCyclicLoess,    normalizeBetweenArrays  plotMD.Rd: decideTests, plotWithHighlights  roast.DGEGLM.Rd: roast, mroast  roast.DGEList.Rd: roast, mroast  romer.DGEGLM.Rd: ids2indices, romer.default, romer  romer.DGEList.Rd: ids2indices, romer.default, romer  sampleWeights.Rd: arrayWeights  topTags.Rd: topTable  voomLmFit.Rd: eBayes, voom, lmFit, voomWithQualityWeights,    duplicateCorrelation, arrayWeights, MArrayLM-class,    chooseLowessSpanPlease provide package anchors for all Rd \link{} targets not in thepackage itself and the base packages.* checking for missing documentation entries ... OK* checking for code/documentation mismatches ... OK* checking Rd \usage sections ... OK* checking Rd contents ... OK* checking for unstated dependencies in examples ... OK* checking line endings in C/C++/Fortran sources/headers ... OK* checking line endings in Makefiles ... OK* checking compilation flags in Makevars ... OK* checking for GNU extensions in Makefiles ... OK* checking for portable use of $(BLAS_LIBS) and $(LAPACK_LIBS) ... OK* checking use of PKG_*FLAGS in Makefiles ... OK* checking compiled code ... NOTENote: information on .o files is not available* checking sizes of PDF files under ‘inst/doc’ ... OK* checking files in ‘vignettes’ ... OK* checking examples ... OKExamples with CPU (user + system) or elapsed time > 5s              user system elapsednearestTSS   6.243  0.223   9.138romer.DGEGLM 3.825  0.212   5.505* checking for unstated dependencies in ‘tests’ ... OK* checking tests ...  Running ‘edgeR-Tests.R’  Comparing ‘edgeR-Tests.Rout’ to ‘edgeR-Tests.Rout.save’ ... OK OK* checking for unstated dependencies in vignettes ... OK* checking package vignettes ... OK* checking running R code from vignettes ... SKIPPED* checking re-building of vignette outputs ... SKIPPED* checking PDF version of manual ... OK* DONEStatus: 2 NOTEsSee  ‘/Users/biocbuild/bbs-3.22-bioc/meat/edgeR.Rcheck/00check.log’for details.

################################################################################################################################################################## Running command:######   /Library/Frameworks/R.framework/Resources/bin/R CMD INSTALL edgeR###############################################################################################################################################################* installing to library ‘/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library’* installing *source* package ‘edgeR’ ...** this is package ‘edgeR’ version ‘4.8.1’** using staged installation** libsusing C compiler: ‘Apple clang version 15.0.0 (clang-1500.0.40.1)’using SDK: ‘MacOSX11.3.1.sdk’clang -arch arm64 -std=gnu2x -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c R_exports.c -o R_exports.oclang -arch arm64 -std=gnu2x -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c R_process_hairpin_reads.c -o R_process_hairpin_reads.oclang -arch arm64 -std=gnu2x -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c add_prior_count.c -o add_prior_count.oclang -arch arm64 -std=gnu2x -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c clowess.c -o clowess.oclang -arch arm64 -std=gnu2x -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c compute_apl.c -o compute_apl.oclang -arch arm64 -std=gnu2x -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c compute_cpm.c -o compute_cpm.oclang -arch arm64 -std=gnu2x -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c compute_nbdev.c -o compute_nbdev.oclang -arch arm64 -std=gnu2x -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c exact_test_by_dev.c -o exact_test_by_dev.oclang -arch arm64 -std=gnu2x -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c fmm_spline.c -o fmm_spline.oclang -arch arm64 -std=gnu2x -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c glm.c -o glm.oclang -arch arm64 -std=gnu2x -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c good_turing.c -o good_turing.oclang -arch arm64 -std=gnu2x -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c init.c -o init.oclang -arch arm64 -std=gnu2x -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c interpolator.c -o interpolator.oclang -arch arm64 -std=gnu2x -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c loess_by_col.c -o loess_by_col.oclang -arch arm64 -std=gnu2x -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c object.c -o object.oclang -arch arm64 -std=gnu2x -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c ql_glm.c -o ql_glm.oclang -arch arm64 -std=gnu2x -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c ql_weights.c -o ql_weights.oclang -arch arm64 -std=gnu2x -dynamiclib -Wl,-headerpad_max_install_names -undefined dynamic_lookup -L/Library/Frameworks/R.framework/Resources/lib -L/opt/R/arm64/lib -o edgeR.so R_exports.o R_process_hairpin_reads.o add_prior_count.o clowess.o compute_apl.o compute_cpm.o compute_nbdev.o exact_test_by_dev.o fmm_spline.o glm.o good_turing.o init.o interpolator.o loess_by_col.o object.o ql_glm.o ql_weights.o -L/Library/Frameworks/R.framework/Resources/lib -lRlapack -L/Library/Frameworks/R.framework/Resources/lib -lRblas -L/opt/gfortran/lib/gcc/aarch64-apple-darwin20.0/14.2.0 -L/opt/gfortran/lib -lemutls_w -lheapt_w -lgfortran -lquadmath -F/Library/Frameworks/R.framework/.. -framework Rinstalling to /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/00LOCK-edgeR/00new/edgeR/libs** R** inst** byte-compile and prepare package for lazy loading** help*** installing help indices** building package indices** installing vignettes** testing if installed package can be loaded from temporary location** checking absolute paths in shared objects and dynamic libraries** testing if installed package can be loaded from final location** testing if installed package keeps a record of temporary installation path* DONE (edgeR)

R version 4.5.2 Patched (2025-11-04 r88984) -- "[Not] Part in a Rumble"Copyright (C) 2025 The R Foundation for Statistical ComputingPlatform: aarch64-apple-darwin20R is free software and comes with ABSOLUTELY NO WARRANTY.You are welcome to redistribute it under certain conditions.Type 'license()' or 'licence()' for distribution details.R is a collaborative project with many contributors.Type 'contributors()' for more information and'citation()' on how to cite R or R packages in publications.Type 'demo()' for some demos, 'help()' for on-line help, or'help.start()' for an HTML browser interface to help.Type 'q()' to quit R.> library(edgeR)Loading required package: limma> options(warnPartialMatchArgs=TRUE,warnPartialMatchAttr=TRUE,warnPartialMatchDollar=TRUE)> > set.seed(0); u <- runif(100)> > # generate raw counts from NB, create list object> y <- matrix(rnbinom(80,size=5,mu=10),nrow=20)> y <- rbind(0,c(0,0,2,2),y)> rownames(y) <- paste("Tag",1:nrow(y),sep=".")> d <- DGEList(counts=y,group=rep(1:2,each=2),lib.size=1001:1004)> > filterByExpr(d) Tag.1  Tag.2  Tag.3  Tag.4  Tag.5  Tag.6  Tag.7  Tag.8  Tag.9 Tag.10 Tag.11  FALSE  FALSE  FALSE   TRUE  FALSE  FALSE   TRUE  FALSE  FALSE   TRUE   TRUE Tag.12 Tag.13 Tag.14 Tag.15 Tag.16 Tag.17 Tag.18 Tag.19 Tag.20 Tag.21 Tag.22   TRUE   TRUE  FALSE  FALSE   TRUE   TRUE  FALSE  FALSE   TRUE   TRUE   TRUE > > # estimate common dispersion and find differences in expression> d <- estimateCommonDisp(d)> d$common.dispersion[1] 0.210292> de <- exactTest(d)> summary(de$table)     logFC             logCPM          PValue        Min.   :-1.7266   Min.   :10.96   Min.   :0.01976   1st Qu.:-0.4855   1st Qu.:13.21   1st Qu.:0.33120   Median : 0.2253   Median :13.37   Median :0.56514   Mean   : 0.1877   Mean   :13.26   Mean   :0.54504   3rd Qu.: 0.5258   3rd Qu.:13.70   3rd Qu.:0.81052   Max.   : 4.0861   Max.   :14.31   Max.   :1.00000  > topTags(de)Comparison of groups:  2-1             logFC   logCPM     PValue       FDRTag.17  2.0450964 13.73726 0.01975954 0.4347099Tag.21 -1.7265870 13.38327 0.06131012 0.6744114Tag.6  -1.6329986 12.81479 0.12446044 0.8982100Tag.2   4.0861092 11.54121 0.16331090 0.8982100Tag.16  0.9324996 13.57074 0.29050785 0.9655885Tag.20  0.8543138 13.76364 0.31736609 0.9655885Tag.12  0.7081170 14.31389 0.37271028 0.9655885Tag.19 -0.7976602 13.31405 0.40166354 0.9655885Tag.3  -0.7300410 13.54155 0.42139935 0.9655885Tag.8  -0.7917906 12.86353 0.47117217 0.9655885> > d2 <- estimateTagwiseDisp(d,trend="none",prior.df=20)> summary(d2$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1757  0.1896  0.1989  0.2063  0.2185  0.2677 > de <- exactTest(d2,dispersion="common")> topTags(de)Comparison of groups:  2-1             logFC   logCPM     PValue       FDRTag.17  2.0450964 13.73726 0.01975954 0.4347099Tag.21 -1.7265870 13.38327 0.06131012 0.6744114Tag.6  -1.6329986 12.81479 0.12446044 0.8982100Tag.2   4.0861092 11.54121 0.16331090 0.8982100Tag.16  0.9324996 13.57074 0.29050785 0.9655885Tag.20  0.8543138 13.76364 0.31736609 0.9655885Tag.12  0.7081170 14.31389 0.37271028 0.9655885Tag.19 -0.7976602 13.31405 0.40166354 0.9655885Tag.3  -0.7300410 13.54155 0.42139935 0.9655885Tag.8  -0.7917906 12.86353 0.47117217 0.9655885> > de <- exactTest(d2)> topTags(de)Comparison of groups:  2-1             logFC   logCPM     PValue       FDRTag.17  2.0450987 13.73726 0.01327001 0.2919403Tag.21 -1.7265897 13.38327 0.05683886 0.6252275Tag.6  -1.6329910 12.81479 0.11460208 0.8404152Tag.2   4.0861092 11.54121 0.16126207 0.8869414Tag.16  0.9324975 13.57074 0.28103256 0.9669238Tag.20  0.8543178 13.76364 0.30234789 0.9669238Tag.12  0.7081149 14.31389 0.37917895 0.9669238Tag.19 -0.7976633 13.31405 0.40762735 0.9669238Tag.3  -0.7300478 13.54155 0.40856822 0.9669238Tag.8  -0.7918243 12.86353 0.49005179 0.9669238> > d2 <- estimateTagwiseDisp(d,trend="movingave",span=0.4,prior.df=20)> summary(d2$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1005  0.1629  0.2064  0.2077  0.2585  0.3164 > de <- exactTest(d2)> topTags(de)Comparison of groups:  2-1             logFC   logCPM     PValue       FDRTag.17  2.0450951 13.73726 0.02427872 0.5341319Tag.21 -1.7265927 13.38327 0.05234833 0.5758316Tag.6  -1.6330014 12.81479 0.12846308 0.8954397Tag.2   4.0861092 11.54121 0.16280722 0.8954397Tag.16  0.9324887 13.57074 0.24308201 0.9711975Tag.20  0.8543044 13.76364 0.35534649 0.9711975Tag.19 -0.7976535 13.31405 0.38873717 0.9711975Tag.3  -0.7300525 13.54155 0.40001438 0.9711975Tag.12  0.7080985 14.31389 0.43530227 0.9711975Tag.8  -0.7918376 12.86353 0.49782701 0.9711975> > summary(exactTest(d2,rejection.region="smallp")$table$PValue)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 0.02428 0.36369 0.55662 0.54319 0.78889 1.00000 > summary(exactTest(d2,rejection.region="deviance")$table$PValue)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 0.02428 0.36369 0.55662 0.54319 0.78889 1.00000 > > d2 <- estimateTagwiseDisp(d,trend="loess",span=0.8,prior.df=20)> summary(d2$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1165  0.1449  0.1832  0.1848  0.2116  0.2825 > de <- exactTest(d2)> topTags(de)Comparison of groups:  2-1             logFC   logCPM     PValue       FDRTag.17  2.0450979 13.73726 0.01546795 0.3402949Tag.21 -1.7266049 13.38327 0.03545446 0.3899990Tag.6  -1.6329841 12.81479 0.10632987 0.7797524Tag.2   4.0861092 11.54121 0.16057893 0.8831841Tag.16  0.9324935 13.57074 0.26348818 0.9658389Tag.20  0.8543140 13.76364 0.31674090 0.9658389Tag.19 -0.7976354 13.31405 0.35564858 0.9658389Tag.3  -0.7300593 13.54155 0.38833737 0.9658389Tag.12  0.7081041 14.31389 0.41513004 0.9658389Tag.8  -0.7918152 12.86353 0.48483449 0.9658389> > d2 <- estimateTagwiseDisp(d,trend="tricube",span=0.8,prior.df=20)> summary(d2$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1165  0.1449  0.1832  0.1848  0.2116  0.2825 > de <- exactTest(d2)> topTags(de)Comparison of groups:  2-1             logFC   logCPM     PValue       FDRTag.17  2.0450979 13.73726 0.01546795 0.3402949Tag.21 -1.7266049 13.38327 0.03545446 0.3899990Tag.6  -1.6329841 12.81479 0.10632987 0.7797524Tag.2   4.0861092 11.54121 0.16057893 0.8831841Tag.16  0.9324935 13.57074 0.26348818 0.9658389Tag.20  0.8543140 13.76364 0.31674090 0.9658389Tag.19 -0.7976354 13.31405 0.35564858 0.9658389Tag.3  -0.7300593 13.54155 0.38833737 0.9658389Tag.12  0.7081041 14.31389 0.41513004 0.9658389Tag.8  -0.7918152 12.86353 0.48483449 0.9658389> > # mglmOneWay> design <- model.matrix(~group,data=d$samples)> mglmOneWay(d[1:10,],design,dispersion=0.2)$coefficients         (Intercept)        group2Tag.1  -1.000000e+08  0.000000e+00Tag.2  -1.000000e+08  1.000000e+08Tag.3   2.525729e+00 -5.108256e-01Tag.4   2.525729e+00  1.484200e-01Tag.5   2.140066e+00 -1.941560e-01Tag.6   2.079442e+00 -1.163151e+00Tag.7   2.014903e+00  2.363888e-01Tag.8   1.945910e+00 -5.596158e-01Tag.9   1.504077e+00  2.006707e-01Tag.10  2.302585e+00  2.623643e-01$fitted.values       Sample1 Sample2 Sample3 Sample4Tag.1      0.0     0.0     0.0     0.0Tag.2      0.0     0.0     2.0     2.0Tag.3     12.5    12.5     7.5     7.5Tag.4     12.5    12.5    14.5    14.5Tag.5      8.5     8.5     7.0     7.0Tag.6      8.0     8.0     2.5     2.5Tag.7      7.5     7.5     9.5     9.5Tag.8      7.0     7.0     4.0     4.0Tag.9      4.5     4.5     5.5     5.5Tag.10    10.0    10.0    13.0    13.0> mglmOneWay(d[1:10,],design,dispersion=0)$coefficients         (Intercept)        group2Tag.1  -1.000000e+08  0.000000e+00Tag.2  -1.000000e+08  1.000000e+08Tag.3   2.525729e+00 -5.108256e-01Tag.4   2.525729e+00  1.484200e-01Tag.5   2.140066e+00 -1.941560e-01Tag.6   2.079442e+00 -1.163151e+00Tag.7   2.014903e+00  2.363888e-01Tag.8   1.945910e+00 -5.596158e-01Tag.9   1.504077e+00  2.006707e-01Tag.10  2.302585e+00  2.623643e-01$fitted.values       Sample1 Sample2 Sample3 Sample4Tag.1      0.0     0.0     0.0     0.0Tag.2      0.0     0.0     2.0     2.0Tag.3     12.5    12.5     7.5     7.5Tag.4     12.5    12.5    14.5    14.5Tag.5      8.5     8.5     7.0     7.0Tag.6      8.0     8.0     2.5     2.5Tag.7      7.5     7.5     9.5     9.5Tag.8      7.0     7.0     4.0     4.0Tag.9      4.5     4.5     5.5     5.5Tag.10    10.0    10.0    13.0    13.0> > fit <- glmFit(d,design,dispersion=d$common.dispersion,prior.count=0.5/4)> lrt <- glmLRT(fit,coef=2)> topTags(lrt)Coefficient:  group2             logFC   logCPM        LR     PValue       FDRTag.17  2.0450964 13.73726 6.0485417 0.01391779 0.3058698Tag.2   4.0861092 11.54121 4.8400340 0.02780635 0.3058698Tag.21 -1.7265870 13.38327 4.0777825 0.04345065 0.3186381Tag.6  -1.6329986 12.81479 3.0078205 0.08286364 0.4557500Tag.16  0.9324996 13.57074 1.3477682 0.24566867 0.8276702Tag.20  0.8543138 13.76364 1.1890032 0.27553071 0.8276702Tag.19 -0.7976602 13.31405 0.9279151 0.33540526 0.8276702Tag.12  0.7081170 14.31389 0.9095513 0.34023349 0.8276702Tag.3  -0.7300410 13.54155 0.8300307 0.36226364 0.8276702Tag.8  -0.7917906 12.86353 0.7830377 0.37621371 0.8276702> > fit <- glmFit(d,design,dispersion=d$common.dispersion,prior.count=0.5)> summary(fit$coefficients)  (Intercept)         group2         Min.   :-7.604   Min.   :-1.13681   1st Qu.:-4.895   1st Qu.:-0.32341   Median :-4.713   Median : 0.15083   Mean   :-4.940   Mean   : 0.07817   3rd Qu.:-4.524   3rd Qu.: 0.35163   Max.   :-4.107   Max.   : 1.60864  > > fit <- glmFit(d,design,prior.count=0.5/4)> lrt <- glmLRT(fit,coef=2)> topTags(lrt)Coefficient:  group2             logFC   logCPM        LR     PValue       FDRTag.17  2.0450964 13.73726 6.0485417 0.01391779 0.3058698Tag.2   4.0861092 11.54121 4.8400340 0.02780635 0.3058698Tag.21 -1.7265870 13.38327 4.0777825 0.04345065 0.3186381Tag.6  -1.6329986 12.81479 3.0078205 0.08286364 0.4557500Tag.16  0.9324996 13.57074 1.3477682 0.24566867 0.8276702Tag.20  0.8543138 13.76364 1.1890032 0.27553071 0.8276702Tag.19 -0.7976602 13.31405 0.9279151 0.33540526 0.8276702Tag.12  0.7081170 14.31389 0.9095513 0.34023349 0.8276702Tag.3  -0.7300410 13.54155 0.8300307 0.36226364 0.8276702Tag.8  -0.7917906 12.86353 0.7830377 0.37621371 0.8276702> > dglm <- estimateGLMCommonDisp(d,design)> dglm$common.dispersion[1] 0.2033282> dglm <- estimateGLMTagwiseDisp(dglm,design,prior.df=20)> summary(dglm$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1756  0.1879  0.1998  0.2031  0.2135  0.2578 > fit <- glmFit(dglm,design,prior.count=0.5/4)> lrt <- glmLRT(fit,coef=2)> topTags(lrt)Coefficient:  group2             logFC   logCPM        LR      PValue       FDRTag.17  2.0450988 13.73727 6.8001118 0.009115216 0.2005348Tag.2   4.0861092 11.54122 4.8594088 0.027495756 0.2872068Tag.21 -1.7265904 13.38327 4.2537154 0.039164558 0.2872068Tag.6  -1.6329904 12.81479 3.1763761 0.074710253 0.4109064Tag.16  0.9324970 13.57074 1.4126709 0.234613512 0.8499599Tag.20  0.8543183 13.76364 1.2721097 0.259371274 0.8499599Tag.19 -0.7976614 13.31405 0.9190392 0.337727381 0.8499599Tag.12  0.7081163 14.31389 0.9014515 0.342392806 0.8499599Tag.3  -0.7300488 13.54155 0.8817937 0.347710872 0.8499599Tag.8  -0.7918166 12.86353 0.7356185 0.391068049 0.8603497> dglm <- estimateGLMTrendedDisp(dglm,design)> summary(dglm$trended.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1522  0.1676  0.1740  0.1887  0.2000  0.3469 > dglm <- estimateGLMTrendedDisp(dglm,design,method="power")> summary(dglm$trended.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1522  0.1676  0.1740  0.1887  0.2000  0.3469 > dglm <- estimateGLMTrendedDisp(dglm,design,method="spline")> summary(dglm$trended.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 0.09353 0.11082 0.15463 0.19006 0.23050 0.52006 > dglm <- estimateGLMTrendedDisp(dglm,design,method="bin.spline")> summary(dglm$trended.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1997  0.1997  0.1997  0.1997  0.1997  0.1997 > dglm <- estimateGLMTagwiseDisp(dglm,design,prior.df=20)> summary(dglm$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1385  0.1792  0.1964  0.1935  0.2026  0.2709 > > dglm2 <- estimateDisp(dglm, design)> summary(dglm2$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1652  0.1740  0.1821  0.1852  0.1909  0.2259 > dglm2 <- estimateDisp(dglm, design, prior.df=20)> summary(dglm2$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1527  0.1669  0.1814  0.1858  0.1951  0.2497 > dglm2 <- estimateDisp(dglm, design, robust=TRUE)> summary(dglm2$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1652  0.1735  0.1822  0.1854  0.1905  0.2280 > > # Continuous trend> nlibs <- 3> ntags <- 1000> dispersion.true <- 0.1> # Make first transcript respond to covariate x> x <- 0:2> design <- model.matrix(~x)> beta.true <- cbind(Beta1=2,Beta2=c(2,rep(0,ntags-1)))> mu.true <- 2^(beta.true %*% t(design))> # Generate count data> y <- rnbinom(ntags*nlibs,mu=mu.true,size=1/dispersion.true)> y <- matrix(y,ntags,nlibs)> colnames(y) <- c("x0","x1","x2")> rownames(y) <- paste("Gene",1:ntags,sep="")> d <- DGEList(y)> d <- normLibSizes(d,method="TMM")> fit <- glmFit(d, design, dispersion=dispersion.true, prior.count=0.5/3)> results <- glmLRT(fit, coef=2)> topTags(results)Coefficient:  x             logFC   logCPM        LR       PValue          FDRGene1    2.907024 13.56183 38.738512 4.845536e-10 4.845536e-07Gene61   2.855317 10.27136 10.738307 1.049403e-03 5.247015e-01Gene62  -2.123902 10.53174  8.818703 2.981585e-03 8.334760e-01Gene134 -1.949073 10.53355  8.125889 4.363759e-03 8.334760e-01Gene740 -1.610046 10.94907  8.013408 4.643227e-03 8.334760e-01Gene354  2.022698 10.45066  7.826308 5.149118e-03 8.334760e-01Gene5    1.856816 10.45249  7.214238 7.232750e-03 8.334760e-01Gene746 -1.798331 10.53094  6.846262 8.882693e-03 8.334760e-01Gene110  1.623148 10.68607  6.737984 9.438120e-03 8.334760e-01Gene383  1.637140 10.75412  6.687530 9.708965e-03 8.334760e-01> d1 <- estimateGLMCommonDisp(d, design, verbose=TRUE)Disp = 0.10253 , BCV = 0.3202 > glmFit(d,design,dispersion=dispersion.true, prior.count=0.5/3)An object of class "DGEGLM"$coefficients      (Intercept)          xGene1   -7.391745  2.0149958Gene2   -7.318483 -0.7611895Gene3   -6.831702 -0.1399478Gene4   -7.480255  0.5172002Gene5   -8.747793  1.2870467995 more rows ...$fitted.values             x0        x1          x2Gene1 2.3570471 18.954454 138.2791328Gene2 2.5138172  1.089292   0.4282107Gene3 4.1580452  3.750528   3.0690081Gene4 2.1012460  3.769592   6.1349937Gene5 0.5080377  2.136398   8.1502486995 more rows ...$deviance[1] 6.38037545 1.46644913 1.38532340 0.01593969 1.03894513995 more elements ...$iter[1] 8 4 4 4 6995 more elements ...$failed[1] FALSE FALSE FALSE FALSE FALSE995 more elements ...$method[1] "levenberg"$counts      x0 x1  x2Gene1  0 30 110Gene2  2  2   0Gene3  3  6   2Gene4  2  4   6Gene5  1  1   9995 more rows ...$unshrunk.coefficients      (Intercept)          xGene1   -7.437763  2.0412762Gene2   -7.373370 -0.8796273Gene3   -6.870127 -0.1465014Gene4   -7.552642  0.5410832Gene5   -8.972372  1.3929679995 more rows ...$df.residual[1] 1 1 1 1 1995 more elements ...$design  (Intercept) x1           1 02           1 13           1 2attr(,"assign")[1] 0 1$offset         [,1]     [,2]     [,3][1,] 8.295172 8.338525 8.284484attr(,"class")[1] "CompressedMatrix"attr(,"Dims")[1] 5 3attr(,"repeat.row")[1] TRUEattr(,"repeat.col")[1] FALSE995 more rows ...$dispersion[1] 0.1$prior.count[1] 0.1666667$samples   group lib.size norm.factorsx0     1     4001    1.0008730x1     1     4176    1.0014172x2     1     3971    0.9977138$AveLogCPM[1] 13.561832  9.682757 10.447014 10.532113 10.452489995 more elements ...> > normLibSizes(d$counts,method="TMMwsp")       x0        x1        x2 0.9992437 1.0077007 0.9931093 > > d2 <- estimateDisp(d, design)> summary(d2$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 0.06382 0.09689 0.11782 0.11373 0.12745 0.13507 > d2 <- estimateDisp(d, design, prior.df=20)> summary(d2$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 0.04862 0.09629 0.11283 0.11358 0.12412 0.37107 > d2 <- estimateDisp(d, design, robust=TRUE)> summary(d2$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 0.06382 0.09689 0.11782 0.11373 0.12745 0.13507 > > # Exact tests> y <- matrix(rnbinom(20,mu=10,size=3/2),5,4)> group <- factor(c(1,1,2,2))> ys <- splitIntoGroupsPseudo(y,group,pair=c(1,2))> exactTestDoubleTail(ys$y1,ys$y2,dispersion=2/3)[1] 0.1334396 0.6343568 0.7280015 0.7124912 0.3919258> > y <- matrix(rnbinom(5*7,mu=10,size=3/2),5,7)> group <- factor(c(1,1,2,2,3,3,3))> ys <- splitIntoGroupsPseudo(y,group,pair=c(1,3))> exactTestDoubleTail(ys$y1,ys$y2,dispersion=2/3)[1] 1.0000000 0.4486382 1.0000000 0.9390317 0.4591241> exactTestBetaApprox(ys$y1,ys$y2,dispersion=2/3)[1] 1.0000000 0.4492969 1.0000000 0.9421695 0.4589194> > y[1,3:4] <- 0> design <- model.matrix(~group)> fit <- glmFit(y,design,dispersion=2/3,prior.count=0.5/7)> summary(fit$coefficients)  (Intercept)         group2            group3         Min.   :-1.817   Min.   :-5.0171   Min.   :-0.64646   1st Qu.:-1.812   1st Qu.:-1.1565   1st Qu.:-0.13919   Median :-1.712   Median : 0.1994   Median :-0.10441   Mean   :-1.625   Mean   :-0.9523   Mean   :-0.04217   3rd Qu.:-1.429   3rd Qu.: 0.3755   3rd Qu.:-0.04305   Max.   :-1.356   Max.   : 0.8374   Max.   : 0.72227  > > lrt <- glmLRT(fit,contrast=cbind(c(0,1,0),c(0,0,1)))> topTags(lrt)Coefficient:  LR test on 2 degrees of freedom      logFC.1    logFC.2   logCPM         LR      PValue        FDR1 -7.2381060 -0.0621100 17.19071 10.7712171 0.004582051 0.022910265 -1.6684268 -0.9326507 17.33529  1.7309951 0.420842115 0.909679672  1.2080938  1.0420198 18.24544  1.0496688 0.591653347 0.909679674  0.5416704 -0.1506381 17.57744  0.3958596 0.820427427 0.909679673  0.2876249 -0.2008143 18.06216  0.1893255 0.909679672 0.90967967> design <- model.matrix(~0+group)> fit <- glmFit(y,design,dispersion=2/3,prior.count=0.5/7)> lrt <- glmLRT(fit,contrast=cbind(c(-1,1,0),c(0,-1,1),c(-1,0,1)))> topTags(lrt)Coefficient:  LR test on 2 degrees of freedom      logFC.1    logFC.2    logFC.3   logCPM         LR      PValue        FDR1 -7.2381060  7.1759960 -0.0621100 17.19071 10.7712171 0.004582051 0.022910265 -1.6684268  0.7357761 -0.9326507 17.33529  1.7309951 0.420842115 0.909679672  1.2080938 -0.1660740  1.0420198 18.24544  1.0496688 0.591653347 0.909679674  0.5416704 -0.6923084 -0.1506381 17.57744  0.3958596 0.820427427 0.909679673  0.2876249 -0.4884392 -0.2008143 18.06216  0.1893255 0.909679672 0.90967967> > # simple Good-Turing algorithm runs.> test1 <- 1:9> freq1 <- c(2018046, 449721, 188933, 105668, 68379, 48190, 35709, 37710, 22280)> goodTuring(rep(test1, freq1))$P0[1] 0.3814719$proportion[1] 8.035111e-08 2.272143e-07 4.060582e-07 5.773690e-07 7.516705e-07[6] 9.276808e-07 1.104759e-06 1.282549e-06 1.460837e-06$count[1] 1 2 3 4 5 6 7 8 9$n[1] 2018046  449721  188933  105668   68379   48190   35709   37710   22280$n0[1] 0> test2 <- c(312, 14491, 16401, 65124, 129797, 323321, 366051, 368599, 405261, 604962)> goodTuring(test2)$P0[1] 0$proportion [1] 0.0001362656 0.0063162959 0.0071487846 0.0283850925 0.0565733349 [6] 0.1409223124 0.1595465235 0.1606570896 0.1766365144 0.2636777866$count [1]    312  14491  16401  65124 129797 323321 366051 368599 405261 604962$n [1] 1 1 1 1 1 1 1 1 1 1$n0[1] 0> > # Dispersion estimation with fitted values equal to zero> ngenes <- 100> nsamples <- 3> y <- matrix(rnbinom(ngenes*nsamples,size=5,mu=10),ngenes,nsamples)> Group <- factor(c(1,2,2))> design <- model.matrix(~Group)> y[1:5,2:3] <- 0> y <- DGEList(counts=y,group=Group)> > fit <- glmQLFit(y,design,legacy=FALSE,keep.unit.mat=TRUE)> fit$dispersion[1] 0.1913324> summary(fit$s2.post)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.3227  0.7307  1.0467  1.1891  1.6050  3.2255 > fit$unit.deviance.adj[1:10,]   Sample1      Sample2      Sample31        0 0.000000e+00 0.000000e+002        0 0.000000e+00 0.000000e+003        0 0.000000e+00 0.000000e+004        0 0.000000e+00 0.000000e+005        0 0.000000e+00 0.000000e+006        0 8.141239e-02 9.842350e-027        0 3.047603e-01 2.190646e-018        0 1.749020e-02 1.620417e-029        0 1.098493e+00 6.237134e-0110       0 7.266233e-05 7.242688e-05> fit$unit.df.adj[1:10,]   Sample1   Sample2   Sample31        0 0.0000000 0.00000002        0 0.0000000 0.00000003        0 0.0000000 0.00000004        0 0.0000000 0.00000005        0 0.0000000 0.00000006        0 0.4865169 0.48806177        0 0.4965337 0.49943198        0 0.4945979 0.49557959        0 0.4905011 0.491729410       0 0.4945957 0.4955775> summary(fit$deviance.adj)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 0.00000 0.08228 0.53816 1.09592 1.72337 5.48746 > summary(fit$df.residual.adj)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.0000  0.9756  0.9816  0.9356  0.9882  1.0912 > > # diffSplice> GeneID <- rep(1:10,each=10)> ds <- diffSplice(fit,geneid=GeneID)Total number of exons:  100 Total number of genes:  10 Number of genes with 1 exon:  0 Mean number of exons in a gene:  10 Max number of exons in a gene:  10 > topSplice(ds,test="F")   GeneID NExons         F      P.Value          FDR1       1     10 7.2590306 2.852463e-05 0.00028524638       8     10 1.3723028 2.418392e-01 0.97052431694       4     10 1.0451997 4.278592e-01 0.97052431699       9     10 0.9638621 4.870572e-01 0.97052431695       5     10 0.9063973 5.317373e-01 0.97052431697       7     10 0.8298259 5.940627e-01 0.97052431696       6     10 0.5388506 8.349931e-01 0.97052431693       3     10 0.3743815 9.387844e-01 0.97052431692       2     10 0.3501548 9.499887e-01 0.970524316910     10     10 0.2962360 9.705243e-01 0.9705243169> topSplice(ds,test="simes")   GeneID NExons     P.Value        FDR1       1     10 0.009637957 0.096379578       8     10 0.228726421 0.985857249       9     10 0.387580728 0.985857247       7     10 0.423403151 0.985857244       4     10 0.593224387 0.985857245       5     10 0.726396305 0.985857246       6     10 0.801299346 0.985857242       2     10 0.922046804 0.9858572410     10     10 0.961722989 0.985857243       3     10 0.985857239 0.98585724> topSplice(ds,test="t")   ExonID GeneID     logFC         t     P.Value        FDR5       5      1 -4.103231 -3.607173 0.001253369 0.096379576       6      1  2.920450  3.322100 0.002596601 0.096379578       8      1  2.506536  3.279451 0.002891387 0.096379574       4      1 -3.694004 -3.149148 0.004005347 0.100133662       2      1 -3.519395 -2.962365 0.006341523 0.1268304674     74      8  2.066258  2.392110 0.022872642 0.3468721010     10      1  1.691740  2.388046 0.024281047 0.346872101       1      1 -2.755095 -2.214511 0.035519773 0.384911963       3      1 -2.755095 -2.214511 0.035519773 0.3849119685     85      9 -1.522989 -2.156323 0.038758073 0.38491196> > vfit <- voomLmFit(y,design)> ds <- diffSplice(vfit,geneid=GeneID)Total number of exons:  100 Total number of genes:  10 Number of genes with 1 exon:  0 Mean number of exons in a gene:  10 Max number of exons in a gene:  10 > topSplice(ds,test="F")   GeneID NExons         F      P.Value         FDR1       1     10 5.9148762 0.0005690819 0.0056908195       5     10 0.8456570 0.5831091010 0.9959249169       9     10 0.7383068 0.6708973693 0.9959249167       7     10 0.6718248 0.7260331196 0.9959249164       4     10 0.6611347 0.7348452566 0.9959249168       8     10 0.4939610 0.8640748170 0.9959249166       6     10 0.4207029 0.9110671096 0.9959249163       3     10 0.3260954 0.9578135984 0.9959249162       2     10 0.2599682 0.9796508664 0.99592491610     10     10 0.1651623 0.9959249163 0.995924916> topSplice(ds,test="simes")   GeneID NExons     P.Value        FDR1       1     10 0.008337696 0.083376965       5     10 0.707551263 0.977369249       9     10 0.729667397 0.977369244       4     10 0.830134473 0.977369247       7     10 0.866749023 0.977369246       6     10 0.899707598 0.977369243       3     10 0.908866165 0.977369248       8     10 0.951942595 0.977369242       2     10 0.967451764 0.9773692410     10     10 0.977369236 0.97736924> topSplice(ds,test="t")   ExonID GeneID     logFC         t     P.Value        FDR5       5      1 -4.348335 -3.800229 0.001211500 0.0833769610     10      1  2.370302  3.660050 0.001667539 0.083376964       4      1 -3.781654 -3.167928 0.005071386 0.169046202       2      1 -3.523729 -2.797191 0.011503569 0.287589238       8      1  2.844177  2.247567 0.036690931 0.733818626       6      1  3.269802  1.924552 0.069412814 0.9772240070     70      7 -1.668741 -1.786482 0.086674902 0.9772240085     85      9 -1.920152 -1.633648 0.115397947 0.9772240074     74      8  2.880866  1.482994 0.151103235 0.977224001       1      1 -2.451261 -1.489990 0.152662112 0.97722400> > y <- estimateCommonDisp(y)> y$common.dispersion[1] 0.2407907> y <- estimateGLMCommonDisp(y,design)> y$common.dispersion[1] 0.2181198> y <- estimateGLMTrendedDisp(y,design)> y$trended.dispersion[1:10] [1] 0.3398724 0.2889110 0.3398724 0.2769872 0.2494308 0.2562610 0.2368123 [8] 0.2052054 0.2024856 0.1932597> summary(y$trended.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1814  0.2065  0.2188  0.2276  0.2413  0.3399 > y <- estimateGLMTagwiseDisp(y,design)> y$tagwise.dispersion[1:10] [1] 0.4396666 0.3416568 0.4270057 0.3121297 0.2525866 0.2439067 0.2413469 [8] 0.1702664 0.2125485 0.1349552> summary(y$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1304  0.1799  0.2212  0.2323  0.2678  0.4397 > y <- estimateDisp(y,design)> y$prior.df[1] 2.18181> y$common.dispersion[1] 0.2185181> y$trended.dispersion[1:10] [1] 0.2685493 0.2668941 0.2685493 0.2624686 0.2532971 0.2554438 0.2386225 [8] 0.1977382 0.1975037 0.1967320> summary(y$trended.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1958  0.1979  0.2098  0.2216  0.2441  0.2685 > y$tagwise.dispersion[1:10] [1] 0.2685493 0.2668941 0.2685493 0.2624686 0.2532971 0.1485295 0.1890927 [8] 0.1036788 0.2510456 0.1024813> summary(y$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 0.09756 0.13374 0.19204 0.22241 0.26731 0.60911 > > # glmQLFit> fit <- glmQLFit(y,design,legacy=FALSE,keep.unit.mat=TRUE)> fit$dispersion[1] 0.1963021> summary(fit$s2.post)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.3180  0.7118  1.0289  1.1702  1.5765  3.1867 > fit$unit.deviance.adj[1:10,]   Sample1      Sample2      Sample31        0 0.0000000000 0.000000e+002        0 0.0000000000 0.000000e+003        0 0.0000000000 0.000000e+004        0 0.0000000000 0.000000e+005        0 0.0000000000 0.000000e+006        0 0.0800369771 9.682457e-027        0 0.3016851330 2.165283e-018        0 0.0171385659 1.587435e-029        0 1.0794653385 6.116829e-0110       0 0.0000711943 7.095954e-05> fit$unit.df.adj[1:10,]   Sample1   Sample2   Sample31        0 0.0000000 0.00000002        0 0.0000000 0.00000003        0 0.0000000 0.00000004        0 0.0000000 0.00000005        0 0.0000000 0.00000006        0 0.4865584 0.48808377        0 0.4973603 0.50024828        0 0.4944987 0.49545889        0 0.4904124 0.491618410       0 0.4944965 0.4954567> summary(fit$deviance.adj)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 0.00000 0.08084 0.53006 1.07828 1.69219 5.39198 > summary(fit$df.residual.adj)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.0000  0.9757  0.9814  0.9358  0.9880  1.0941 > fit <- glmQLFit(y,design,legacy=TRUE)> summary(fit$dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1958  0.1979  0.2098  0.2216  0.2441  0.2685 > summary(fit$s2.post)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 0.06182 0.60040 0.79878 0.91179 1.16628 2.59834 > > proc.time()   user  system elapsed   3.290   0.173   4.531

R version 4.5.1 (2025-06-13 ucrt) -- "Great Square Root"Copyright (C) 2025 The R Foundation for Statistical ComputingPlatform: x86_64-w64-mingw32/x64R is free software and comes with ABSOLUTELY NO WARRANTY.You are welcome to redistribute it under certain conditions.Type 'license()' or 'licence()' for distribution details.  Natural language support but running in an English localeR is a collaborative project with many contributors.Type 'contributors()' for more information and'citation()' on how to cite R or R packages in publications.Type 'demo()' for some demos, 'help()' for on-line help, or'help.start()' for an HTML browser interface to help.Type 'q()' to quit R.> library(edgeR)Loading required package: limma> options(warnPartialMatchArgs=TRUE,warnPartialMatchAttr=TRUE,warnPartialMatchDollar=TRUE)> > set.seed(0); u <- runif(100)> > # generate raw counts from NB, create list object> y <- matrix(rnbinom(80,size=5,mu=10),nrow=20)> y <- rbind(0,c(0,0,2,2),y)> rownames(y) <- paste("Tag",1:nrow(y),sep=".")> d <- DGEList(counts=y,group=rep(1:2,each=2),lib.size=1001:1004)> > filterByExpr(d) Tag.1  Tag.2  Tag.3  Tag.4  Tag.5  Tag.6  Tag.7  Tag.8  Tag.9 Tag.10 Tag.11  FALSE  FALSE  FALSE   TRUE  FALSE  FALSE   TRUE  FALSE  FALSE   TRUE   TRUE Tag.12 Tag.13 Tag.14 Tag.15 Tag.16 Tag.17 Tag.18 Tag.19 Tag.20 Tag.21 Tag.22   TRUE   TRUE  FALSE  FALSE   TRUE   TRUE  FALSE  FALSE   TRUE   TRUE   TRUE > > # estimate common dispersion and find differences in expression> d <- estimateCommonDisp(d)> d$common.dispersion[1] 0.210292> de <- exactTest(d)> summary(de$table)     logFC             logCPM          PValue        Min.   :-1.7266   Min.   :10.96   Min.   :0.01976   1st Qu.:-0.4855   1st Qu.:13.21   1st Qu.:0.33120   Median : 0.2253   Median :13.37   Median :0.56514   Mean   : 0.1877   Mean   :13.26   Mean   :0.54504   3rd Qu.: 0.5258   3rd Qu.:13.70   3rd Qu.:0.81052   Max.   : 4.0861   Max.   :14.31   Max.   :1.00000  > topTags(de)Comparison of groups:  2-1             logFC   logCPM     PValue       FDRTag.17  2.0450964 13.73726 0.01975954 0.4347099Tag.21 -1.7265870 13.38327 0.06131012 0.6744114Tag.6  -1.6329986 12.81479 0.12446044 0.8982100Tag.2   4.0861092 11.54121 0.16331090 0.8982100Tag.16  0.9324996 13.57074 0.29050785 0.9655885Tag.20  0.8543138 13.76364 0.31736609 0.9655885Tag.12  0.7081170 14.31389 0.37271028 0.9655885Tag.19 -0.7976602 13.31405 0.40166354 0.9655885Tag.3  -0.7300410 13.54155 0.42139935 0.9655885Tag.8  -0.7917906 12.86353 0.47117217 0.9655885> > d2 <- estimateTagwiseDisp(d,trend="none",prior.df=20)> summary(d2$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1757  0.1896  0.1989  0.2063  0.2185  0.2677 > de <- exactTest(d2,dispersion="common")> topTags(de)Comparison of groups:  2-1             logFC   logCPM     PValue       FDRTag.17  2.0450964 13.73726 0.01975954 0.4347099Tag.21 -1.7265870 13.38327 0.06131012 0.6744114Tag.6  -1.6329986 12.81479 0.12446044 0.8982100Tag.2   4.0861092 11.54121 0.16331090 0.8982100Tag.16  0.9324996 13.57074 0.29050785 0.9655885Tag.20  0.8543138 13.76364 0.31736609 0.9655885Tag.12  0.7081170 14.31389 0.37271028 0.9655885Tag.19 -0.7976602 13.31405 0.40166354 0.9655885Tag.3  -0.7300410 13.54155 0.42139935 0.9655885Tag.8  -0.7917906 12.86353 0.47117217 0.9655885> > de <- exactTest(d2)> topTags(de)Comparison of groups:  2-1             logFC   logCPM     PValue       FDRTag.17  2.0450987 13.73726 0.01327001 0.2919403Tag.21 -1.7265897 13.38327 0.05683886 0.6252275Tag.6  -1.6329910 12.81479 0.11460208 0.8404152Tag.2   4.0861092 11.54121 0.16126207 0.8869414Tag.16  0.9324975 13.57074 0.28103256 0.9669238Tag.20  0.8543178 13.76364 0.30234789 0.9669238Tag.12  0.7081149 14.31389 0.37917895 0.9669238Tag.19 -0.7976633 13.31405 0.40762735 0.9669238Tag.3  -0.7300478 13.54155 0.40856822 0.9669238Tag.8  -0.7918243 12.86353 0.49005179 0.9669238> > d2 <- estimateTagwiseDisp(d,trend="movingave",span=0.4,prior.df=20)> summary(d2$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1005  0.1629  0.2064  0.2077  0.2585  0.3164 > de <- exactTest(d2)> topTags(de)Comparison of groups:  2-1             logFC   logCPM     PValue       FDRTag.17  2.0450951 13.73726 0.02427872 0.5341319Tag.21 -1.7265927 13.38327 0.05234833 0.5758316Tag.6  -1.6330014 12.81479 0.12846308 0.8954397Tag.2   4.0861092 11.54121 0.16280722 0.8954397Tag.16  0.9324887 13.57074 0.24308201 0.9711975Tag.20  0.8543044 13.76364 0.35534649 0.9711975Tag.19 -0.7976535 13.31405 0.38873717 0.9711975Tag.3  -0.7300525 13.54155 0.40001438 0.9711975Tag.12  0.7080985 14.31389 0.43530227 0.9711975Tag.8  -0.7918376 12.86353 0.49782701 0.9711975> > summary(exactTest(d2,rejection.region="smallp")$table$PValue)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 0.02428 0.36369 0.55662 0.54319 0.78889 1.00000 > summary(exactTest(d2,rejection.region="deviance")$table$PValue)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 0.02428 0.36369 0.55662 0.54319 0.78889 1.00000 > > d2 <- estimateTagwiseDisp(d,trend="loess",span=0.8,prior.df=20)> summary(d2$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1165  0.1449  0.1832  0.1848  0.2116  0.2825 > de <- exactTest(d2)> topTags(de)Comparison of groups:  2-1             logFC   logCPM     PValue       FDRTag.17  2.0450979 13.73726 0.01546795 0.3402949Tag.21 -1.7266049 13.38327 0.03545446 0.3899990Tag.6  -1.6329841 12.81479 0.10632987 0.7797524Tag.2   4.0861092 11.54121 0.16057893 0.8831841Tag.16  0.9324935 13.57074 0.26348818 0.9658389Tag.20  0.8543140 13.76364 0.31674090 0.9658389Tag.19 -0.7976354 13.31405 0.35564858 0.9658389Tag.3  -0.7300593 13.54155 0.38833737 0.9658389Tag.12  0.7081041 14.31389 0.41513004 0.9658389Tag.8  -0.7918152 12.86353 0.48483449 0.9658389> > d2 <- estimateTagwiseDisp(d,trend="tricube",span=0.8,prior.df=20)> summary(d2$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1165  0.1449  0.1832  0.1848  0.2116  0.2825 > de <- exactTest(d2)> topTags(de)Comparison of groups:  2-1             logFC   logCPM     PValue       FDRTag.17  2.0450979 13.73726 0.01546795 0.3402949Tag.21 -1.7266049 13.38327 0.03545446 0.3899990Tag.6  -1.6329841 12.81479 0.10632987 0.7797524Tag.2   4.0861092 11.54121 0.16057893 0.8831841Tag.16  0.9324935 13.57074 0.26348818 0.9658389Tag.20  0.8543140 13.76364 0.31674090 0.9658389Tag.19 -0.7976354 13.31405 0.35564858 0.9658389Tag.3  -0.7300593 13.54155 0.38833737 0.9658389Tag.12  0.7081041 14.31389 0.41513004 0.9658389Tag.8  -0.7918152 12.86353 0.48483449 0.9658389> > # mglmOneWay> design <- model.matrix(~group,data=d$samples)> mglmOneWay(d[1:10,],design,dispersion=0.2)$coefficients         (Intercept)        group2Tag.1  -1.000000e+08  0.000000e+00Tag.2  -1.000000e+08  1.000000e+08Tag.3   2.525729e+00 -5.108256e-01Tag.4   2.525729e+00  1.484200e-01Tag.5   2.140066e+00 -1.941560e-01Tag.6   2.079442e+00 -1.163151e+00Tag.7   2.014903e+00  2.363888e-01Tag.8   1.945910e+00 -5.596158e-01Tag.9   1.504077e+00  2.006707e-01Tag.10  2.302585e+00  2.623643e-01$fitted.values       Sample1 Sample2 Sample3 Sample4Tag.1      0.0     0.0     0.0     0.0Tag.2      0.0     0.0     2.0     2.0Tag.3     12.5    12.5     7.5     7.5Tag.4     12.5    12.5    14.5    14.5Tag.5      8.5     8.5     7.0     7.0Tag.6      8.0     8.0     2.5     2.5Tag.7      7.5     7.5     9.5     9.5Tag.8      7.0     7.0     4.0     4.0Tag.9      4.5     4.5     5.5     5.5Tag.10    10.0    10.0    13.0    13.0> mglmOneWay(d[1:10,],design,dispersion=0)$coefficients         (Intercept)        group2Tag.1  -1.000000e+08  0.000000e+00Tag.2  -1.000000e+08  1.000000e+08Tag.3   2.525729e+00 -5.108256e-01Tag.4   2.525729e+00  1.484200e-01Tag.5   2.140066e+00 -1.941560e-01Tag.6   2.079442e+00 -1.163151e+00Tag.7   2.014903e+00  2.363888e-01Tag.8   1.945910e+00 -5.596158e-01Tag.9   1.504077e+00  2.006707e-01Tag.10  2.302585e+00  2.623643e-01$fitted.values       Sample1 Sample2 Sample3 Sample4Tag.1      0.0     0.0     0.0     0.0Tag.2      0.0     0.0     2.0     2.0Tag.3     12.5    12.5     7.5     7.5Tag.4     12.5    12.5    14.5    14.5Tag.5      8.5     8.5     7.0     7.0Tag.6      8.0     8.0     2.5     2.5Tag.7      7.5     7.5     9.5     9.5Tag.8      7.0     7.0     4.0     4.0Tag.9      4.5     4.5     5.5     5.5Tag.10    10.0    10.0    13.0    13.0> > fit <- glmFit(d,design,dispersion=d$common.dispersion,prior.count=0.5/4)> lrt <- glmLRT(fit,coef=2)> topTags(lrt)Coefficient:  group2             logFC   logCPM        LR     PValue       FDRTag.17  2.0450964 13.73726 6.0485417 0.01391779 0.3058698Tag.2   4.0861092 11.54121 4.8400340 0.02780635 0.3058698Tag.21 -1.7265870 13.38327 4.0777825 0.04345065 0.3186381Tag.6  -1.6329986 12.81479 3.0078205 0.08286364 0.4557500Tag.16  0.9324996 13.57074 1.3477682 0.24566867 0.8276702Tag.20  0.8543138 13.76364 1.1890032 0.27553071 0.8276702Tag.19 -0.7976602 13.31405 0.9279151 0.33540526 0.8276702Tag.12  0.7081170 14.31389 0.9095513 0.34023349 0.8276702Tag.3  -0.7300410 13.54155 0.8300307 0.36226364 0.8276702Tag.8  -0.7917906 12.86353 0.7830377 0.37621371 0.8276702> > fit <- glmFit(d,design,dispersion=d$common.dispersion,prior.count=0.5)> summary(fit$coefficients)  (Intercept)         group2         Min.   :-7.604   Min.   :-1.13681   1st Qu.:-4.895   1st Qu.:-0.32341   Median :-4.713   Median : 0.15083   Mean   :-4.940   Mean   : 0.07817   3rd Qu.:-4.524   3rd Qu.: 0.35163   Max.   :-4.107   Max.   : 1.60864  > > fit <- glmFit(d,design,prior.count=0.5/4)> lrt <- glmLRT(fit,coef=2)> topTags(lrt)Coefficient:  group2             logFC   logCPM        LR     PValue       FDRTag.17  2.0450964 13.73726 6.0485417 0.01391779 0.3058698Tag.2   4.0861092 11.54121 4.8400340 0.02780635 0.3058698Tag.21 -1.7265870 13.38327 4.0777825 0.04345065 0.3186381Tag.6  -1.6329986 12.81479 3.0078205 0.08286364 0.4557500Tag.16  0.9324996 13.57074 1.3477682 0.24566867 0.8276702Tag.20  0.8543138 13.76364 1.1890032 0.27553071 0.8276702Tag.19 -0.7976602 13.31405 0.9279151 0.33540526 0.8276702Tag.12  0.7081170 14.31389 0.9095513 0.34023349 0.8276702Tag.3  -0.7300410 13.54155 0.8300307 0.36226364 0.8276702Tag.8  -0.7917906 12.86353 0.7830377 0.37621371 0.8276702> > dglm <- estimateGLMCommonDisp(d,design)> dglm$common.dispersion[1] 0.2033282> dglm <- estimateGLMTagwiseDisp(dglm,design,prior.df=20)> summary(dglm$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1756  0.1879  0.1998  0.2031  0.2135  0.2578 > fit <- glmFit(dglm,design,prior.count=0.5/4)> lrt <- glmLRT(fit,coef=2)> topTags(lrt)Coefficient:  group2             logFC   logCPM        LR      PValue       FDRTag.17  2.0450988 13.73727 6.8001118 0.009115216 0.2005348Tag.2   4.0861092 11.54122 4.8594088 0.027495756 0.2872068Tag.21 -1.7265904 13.38327 4.2537154 0.039164558 0.2872068Tag.6  -1.6329904 12.81479 3.1763761 0.074710253 0.4109064Tag.16  0.9324970 13.57074 1.4126709 0.234613512 0.8499599Tag.20  0.8543183 13.76364 1.2721097 0.259371274 0.8499599Tag.19 -0.7976614 13.31405 0.9190392 0.337727381 0.8499599Tag.12  0.7081163 14.31389 0.9014515 0.342392806 0.8499599Tag.3  -0.7300488 13.54155 0.8817937 0.347710872 0.8499599Tag.8  -0.7918166 12.86353 0.7356185 0.391068049 0.8603497> dglm <- estimateGLMTrendedDisp(dglm,design)> summary(dglm$trended.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1522  0.1676  0.1740  0.1887  0.2000  0.3469 > dglm <- estimateGLMTrendedDisp(dglm,design,method="power")> summary(dglm$trended.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1522  0.1676  0.1740  0.1887  0.2000  0.3469 > dglm <- estimateGLMTrendedDisp(dglm,design,method="spline")> summary(dglm$trended.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 0.09353 0.11082 0.15463 0.19006 0.23050 0.52006 > dglm <- estimateGLMTrendedDisp(dglm,design,method="bin.spline")> summary(dglm$trended.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1997  0.1997  0.1997  0.1997  0.1997  0.1997 > dglm <- estimateGLMTagwiseDisp(dglm,design,prior.df=20)> summary(dglm$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1385  0.1792  0.1964  0.1935  0.2026  0.2709 > > dglm2 <- estimateDisp(dglm, design)> summary(dglm2$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1652  0.1740  0.1821  0.1852  0.1909  0.2259 > dglm2 <- estimateDisp(dglm, design, prior.df=20)> summary(dglm2$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1527  0.1669  0.1814  0.1858  0.1951  0.2497 > dglm2 <- estimateDisp(dglm, design, robust=TRUE)> summary(dglm2$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1652  0.1735  0.1822  0.1854  0.1905  0.2280 > > # Continuous trend> nlibs <- 3> ntags <- 1000> dispersion.true <- 0.1> # Make first transcript respond to covariate x> x <- 0:2> design <- model.matrix(~x)> beta.true <- cbind(Beta1=2,Beta2=c(2,rep(0,ntags-1)))> mu.true <- 2^(beta.true %*% t(design))> # Generate count data> y <- rnbinom(ntags*nlibs,mu=mu.true,size=1/dispersion.true)> y <- matrix(y,ntags,nlibs)> colnames(y) <- c("x0","x1","x2")> rownames(y) <- paste("Gene",1:ntags,sep="")> d <- DGEList(y)> d <- normLibSizes(d,method="TMM")> fit <- glmFit(d, design, dispersion=dispersion.true, prior.count=0.5/3)> results <- glmLRT(fit, coef=2)> topTags(results)Coefficient:  x             logFC   logCPM        LR       PValue          FDRGene1    2.907024 13.56183 38.738512 4.845536e-10 4.845536e-07Gene61   2.855317 10.27136 10.738307 1.049403e-03 5.247015e-01Gene62  -2.123902 10.53174  8.818703 2.981585e-03 8.334760e-01Gene134 -1.949073 10.53355  8.125889 4.363759e-03 8.334760e-01Gene740 -1.610046 10.94907  8.013408 4.643227e-03 8.334760e-01Gene354  2.022698 10.45066  7.826308 5.149118e-03 8.334760e-01Gene5    1.856816 10.45249  7.214238 7.232750e-03 8.334760e-01Gene746 -1.798331 10.53094  6.846262 8.882693e-03 8.334760e-01Gene110  1.623148 10.68607  6.737984 9.438120e-03 8.334760e-01Gene383  1.637140 10.75412  6.687530 9.708965e-03 8.334760e-01> d1 <- estimateGLMCommonDisp(d, design, verbose=TRUE)Disp = 0.10253 , BCV = 0.3202 > glmFit(d,design,dispersion=dispersion.true, prior.count=0.5/3)An object of class "DGEGLM"$coefficients      (Intercept)          xGene1   -7.391745  2.0149958Gene2   -7.318483 -0.7611895Gene3   -6.831702 -0.1399478Gene4   -7.480255  0.5172002Gene5   -8.747793  1.2870467995 more rows ...$fitted.values             x0        x1          x2Gene1 2.3570471 18.954454 138.2791328Gene2 2.5138172  1.089292   0.4282107Gene3 4.1580452  3.750528   3.0690081Gene4 2.1012460  3.769592   6.1349937Gene5 0.5080377  2.136398   8.1502486995 more rows ...$deviance[1] 6.38037545 1.46644913 1.38532340 0.01593969 1.03894513995 more elements ...$iter[1] 8 4 4 4 6995 more elements ...$failed[1] FALSE FALSE FALSE FALSE FALSE995 more elements ...$method[1] "levenberg"$counts      x0 x1  x2Gene1  0 30 110Gene2  2  2   0Gene3  3  6   2Gene4  2  4   6Gene5  1  1   9995 more rows ...$unshrunk.coefficients      (Intercept)          xGene1   -7.437763  2.0412762Gene2   -7.373370 -0.8796273Gene3   -6.870127 -0.1465014Gene4   -7.552642  0.5410832Gene5   -8.972372  1.3929679995 more rows ...$df.residual[1] 1 1 1 1 1995 more elements ...$design  (Intercept) x1           1 02           1 13           1 2attr(,"assign")[1] 0 1$offset         [,1]     [,2]     [,3][1,] 8.295172 8.338525 8.284484attr(,"class")[1] "CompressedMatrix"attr(,"Dims")[1] 5 3attr(,"repeat.row")[1] TRUEattr(,"repeat.col")[1] FALSE995 more rows ...$dispersion[1] 0.1$prior.count[1] 0.1666667$samples   group lib.size norm.factorsx0     1     4001    1.0008730x1     1     4176    1.0014172x2     1     3971    0.9977138$AveLogCPM[1] 13.561832  9.682757 10.447014 10.532113 10.452489995 more elements ...> > normLibSizes(d$counts,method="TMMwsp")       x0        x1        x2 0.9992437 1.0077007 0.9931093 > > d2 <- estimateDisp(d, design)> summary(d2$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 0.06382 0.09689 0.11782 0.11373 0.12745 0.13507 > d2 <- estimateDisp(d, design, prior.df=20)> summary(d2$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 0.04862 0.09629 0.11283 0.11358 0.12412 0.37107 > d2 <- estimateDisp(d, design, robust=TRUE)> summary(d2$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 0.06382 0.09689 0.11782 0.11373 0.12745 0.13507 > > # Exact tests> y <- matrix(rnbinom(20,mu=10,size=3/2),5,4)> group <- factor(c(1,1,2,2))> ys <- splitIntoGroupsPseudo(y,group,pair=c(1,2))> exactTestDoubleTail(ys$y1,ys$y2,dispersion=2/3)[1] 0.1334396 0.6343568 0.7280015 0.7124912 0.3919258> > y <- matrix(rnbinom(5*7,mu=10,size=3/2),5,7)> group <- factor(c(1,1,2,2,3,3,3))> ys <- splitIntoGroupsPseudo(y,group,pair=c(1,3))> exactTestDoubleTail(ys$y1,ys$y2,dispersion=2/3)[1] 1.0000000 0.4486382 1.0000000 0.9390317 0.4591241> exactTestBetaApprox(ys$y1,ys$y2,dispersion=2/3)[1] 1.0000000 0.4492969 1.0000000 0.9421695 0.4589194> > y[1,3:4] <- 0> design <- model.matrix(~group)> fit <- glmFit(y,design,dispersion=2/3,prior.count=0.5/7)> summary(fit$coefficients)  (Intercept)         group2            group3         Min.   :-1.817   Min.   :-5.0171   Min.   :-0.64646   1st Qu.:-1.812   1st Qu.:-1.1565   1st Qu.:-0.13919   Median :-1.712   Median : 0.1994   Median :-0.10441   Mean   :-1.625   Mean   :-0.9523   Mean   :-0.04217   3rd Qu.:-1.429   3rd Qu.: 0.3755   3rd Qu.:-0.04305   Max.   :-1.356   Max.   : 0.8374   Max.   : 0.72227  > > lrt <- glmLRT(fit,contrast=cbind(c(0,1,0),c(0,0,1)))> topTags(lrt)Coefficient:  LR test on 2 degrees of freedom      logFC.1    logFC.2   logCPM         LR      PValue        FDR1 -7.2381060 -0.0621100 17.19071 10.7712171 0.004582051 0.022910265 -1.6684268 -0.9326507 17.33529  1.7309951 0.420842115 0.909679672  1.2080938  1.0420198 18.24544  1.0496688 0.591653347 0.909679674  0.5416704 -0.1506381 17.57744  0.3958596 0.820427427 0.909679673  0.2876249 -0.2008143 18.06216  0.1893255 0.909679672 0.90967967> design <- model.matrix(~0+group)> fit <- glmFit(y,design,dispersion=2/3,prior.count=0.5/7)> lrt <- glmLRT(fit,contrast=cbind(c(-1,1,0),c(0,-1,1),c(-1,0,1)))> topTags(lrt)Coefficient:  LR test on 2 degrees of freedom      logFC.1    logFC.2    logFC.3   logCPM         LR      PValue        FDR1 -7.2381060  7.1759960 -0.0621100 17.19071 10.7712171 0.004582051 0.022910265 -1.6684268  0.7357761 -0.9326507 17.33529  1.7309951 0.420842115 0.909679672  1.2080938 -0.1660740  1.0420198 18.24544  1.0496688 0.591653347 0.909679674  0.5416704 -0.6923084 -0.1506381 17.57744  0.3958596 0.820427427 0.909679673  0.2876249 -0.4884392 -0.2008143 18.06216  0.1893255 0.909679672 0.90967967> > # simple Good-Turing algorithm runs.> test1 <- 1:9> freq1 <- c(2018046, 449721, 188933, 105668, 68379, 48190, 35709, 37710, 22280)> goodTuring(rep(test1, freq1))$P0[1] 0.3814719$proportion[1] 8.035111e-08 2.272143e-07 4.060582e-07 5.773690e-07 7.516705e-07[6] 9.276808e-07 1.104759e-06 1.282549e-06 1.460837e-06$count[1] 1 2 3 4 5 6 7 8 9$n[1] 2018046  449721  188933  105668   68379   48190   35709   37710   22280$n0[1] 0> test2 <- c(312, 14491, 16401, 65124, 129797, 323321, 366051, 368599, 405261, 604962)> goodTuring(test2)$P0[1] 0$proportion [1] 0.0001362656 0.0063162959 0.0071487846 0.0283850925 0.0565733349 [6] 0.1409223124 0.1595465235 0.1606570896 0.1766365144 0.2636777866$count [1]    312  14491  16401  65124 129797 323321 366051 368599 405261 604962$n [1] 1 1 1 1 1 1 1 1 1 1$n0[1] 0> > # Dispersion estimation with fitted values equal to zero> ngenes <- 100> nsamples <- 3> y <- matrix(rnbinom(ngenes*nsamples,size=5,mu=10),ngenes,nsamples)> Group <- factor(c(1,2,2))> design <- model.matrix(~Group)> y[1:5,2:3] <- 0> y <- DGEList(counts=y,group=Group)> > fit <- glmQLFit(y,design,legacy=FALSE,keep.unit.mat=TRUE)> fit$dispersion[1] 0.1913324> summary(fit$s2.post)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.3227  0.7307  1.0467  1.1891  1.6050  3.2255 > fit$unit.deviance.adj[1:10,]   Sample1      Sample2      Sample31        0 0.000000e+00 0.000000e+002        0 0.000000e+00 0.000000e+003        0 0.000000e+00 0.000000e+004        0 0.000000e+00 0.000000e+005        0 0.000000e+00 0.000000e+006        0 8.141239e-02 9.842350e-027        0 3.047603e-01 2.190646e-018        0 1.749020e-02 1.620417e-029        0 1.098493e+00 6.237134e-0110       0 7.266233e-05 7.242688e-05> fit$unit.df.adj[1:10,]   Sample1   Sample2   Sample31        0 0.0000000 0.00000002        0 0.0000000 0.00000003        0 0.0000000 0.00000004        0 0.0000000 0.00000005        0 0.0000000 0.00000006        0 0.4865169 0.48806177        0 0.4965337 0.49943198        0 0.4945979 0.49557959        0 0.4905011 0.491729410       0 0.4945957 0.4955775> summary(fit$deviance.adj)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 0.00000 0.08228 0.53816 1.09592 1.72337 5.48746 > summary(fit$df.residual.adj)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.0000  0.9756  0.9816  0.9356  0.9882  1.0912 > > # diffSplice> GeneID <- rep(1:10,each=10)> ds <- diffSplice(fit,geneid=GeneID)Total number of exons:  100 Total number of genes:  10 Number of genes with 1 exon:  0 Mean number of exons in a gene:  10 Max number of exons in a gene:  10 > topSplice(ds,test="F")   GeneID NExons         F      P.Value          FDR1       1     10 7.2590306 2.852463e-05 0.00028524638       8     10 1.3723028 2.418392e-01 0.97052431694       4     10 1.0451997 4.278592e-01 0.97052431699       9     10 0.9638621 4.870572e-01 0.97052431695       5     10 0.9063973 5.317373e-01 0.97052431697       7     10 0.8298259 5.940627e-01 0.97052431696       6     10 0.5388506 8.349931e-01 0.97052431693       3     10 0.3743815 9.387844e-01 0.97052431692       2     10 0.3501548 9.499887e-01 0.970524316910     10     10 0.2962360 9.705243e-01 0.9705243169> topSplice(ds,test="simes")   GeneID NExons     P.Value        FDR1       1     10 0.009637957 0.096379578       8     10 0.228726421 0.985857249       9     10 0.387580728 0.985857247       7     10 0.423403151 0.985857244       4     10 0.593224387 0.985857245       5     10 0.726396305 0.985857246       6     10 0.801299346 0.985857242       2     10 0.922046804 0.9858572410     10     10 0.961722989 0.985857243       3     10 0.985857239 0.98585724> topSplice(ds,test="t")   ExonID GeneID     logFC         t     P.Value        FDR5       5      1 -4.103231 -3.607173 0.001253369 0.096379576       6      1  2.920450  3.322100 0.002596601 0.096379578       8      1  2.506536  3.279451 0.002891387 0.096379574       4      1 -3.694004 -3.149148 0.004005347 0.100133662       2      1 -3.519395 -2.962365 0.006341523 0.1268304674     74      8  2.066258  2.392110 0.022872642 0.3468721010     10      1  1.691740  2.388046 0.024281047 0.346872101       1      1 -2.755095 -2.214511 0.035519773 0.384911963       3      1 -2.755095 -2.214511 0.035519773 0.3849119685     85      9 -1.522989 -2.156323 0.038758073 0.38491196> > vfit <- voomLmFit(y,design)> ds <- diffSplice(vfit,geneid=GeneID)Total number of exons:  100 Total number of genes:  10 Number of genes with 1 exon:  0 Mean number of exons in a gene:  10 Max number of exons in a gene:  10 > topSplice(ds,test="F")   GeneID NExons         F      P.Value         FDR1       1     10 5.9148762 0.0005690819 0.0056908195       5     10 0.8456570 0.5831091010 0.9959249169       9     10 0.7383068 0.6708973693 0.9959249167       7     10 0.6718248 0.7260331196 0.9959249164       4     10 0.6611347 0.7348452566 0.9959249168       8     10 0.4939610 0.8640748170 0.9959249166       6     10 0.4207029 0.9110671096 0.9959249163       3     10 0.3260954 0.9578135984 0.9959249162       2     10 0.2599682 0.9796508664 0.99592491610     10     10 0.1651623 0.9959249163 0.995924916> topSplice(ds,test="simes")   GeneID NExons     P.Value        FDR1       1     10 0.008337696 0.083376965       5     10 0.707551263 0.977369249       9     10 0.729667397 0.977369244       4     10 0.830134473 0.977369247       7     10 0.866749023 0.977369246       6     10 0.899707598 0.977369243       3     10 0.908866165 0.977369248       8     10 0.951942595 0.977369242       2     10 0.967451764 0.9773692410     10     10 0.977369236 0.97736924> topSplice(ds,test="t")   ExonID GeneID     logFC         t     P.Value        FDR5       5      1 -4.348335 -3.800229 0.001211500 0.0833769610     10      1  2.370302  3.660050 0.001667539 0.083376964       4      1 -3.781654 -3.167928 0.005071386 0.169046202       2      1 -3.523729 -2.797191 0.011503569 0.287589238       8      1  2.844177  2.247567 0.036690931 0.733818626       6      1  3.269802  1.924552 0.069412814 0.9772240070     70      7 -1.668741 -1.786482 0.086674902 0.9772240085     85      9 -1.920152 -1.633648 0.115397947 0.9772240074     74      8  2.880866  1.482994 0.151103235 0.977224001       1      1 -2.451261 -1.489990 0.152662112 0.97722400> > y <- estimateCommonDisp(y)> y$common.dispersion[1] 0.2407907> y <- estimateGLMCommonDisp(y,design)> y$common.dispersion[1] 0.2181198> y <- estimateGLMTrendedDisp(y,design)> y$trended.dispersion[1:10] [1] 0.3398724 0.2889110 0.3398724 0.2769872 0.2494308 0.2562610 0.2368123 [8] 0.2052054 0.2024856 0.1932597> summary(y$trended.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1814  0.2065  0.2188  0.2276  0.2413  0.3399 > y <- estimateGLMTagwiseDisp(y,design)> y$tagwise.dispersion[1:10] [1] 0.4396666 0.3416568 0.4270057 0.3121297 0.2525866 0.2439067 0.2413469 [8] 0.1702664 0.2125485 0.1349552> summary(y$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1304  0.1799  0.2212  0.2323  0.2678  0.4397 > y <- estimateDisp(y,design)> y$prior.df[1] 2.18181> y$common.dispersion[1] 0.2185181> y$trended.dispersion[1:10] [1] 0.2685493 0.2668941 0.2685493 0.2624686 0.2532971 0.2554438 0.2386225 [8] 0.1977382 0.1975037 0.1967320> summary(y$trended.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1958  0.1979  0.2098  0.2216  0.2441  0.2685 > y$tagwise.dispersion[1:10] [1] 0.2685493 0.2668941 0.2685493 0.2624686 0.2532971 0.1485295 0.1890927 [8] 0.1036788 0.2510456 0.1024813> summary(y$tagwise.dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 0.09756 0.13374 0.19204 0.22241 0.26731 0.60911 > > # glmQLFit> fit <- glmQLFit(y,design,legacy=FALSE,keep.unit.mat=TRUE)> fit$dispersion[1] 0.1963021> summary(fit$s2.post)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.3180  0.7118  1.0289  1.1702  1.5765  3.1867 > fit$unit.deviance.adj[1:10,]   Sample1      Sample2      Sample31        0 0.0000000000 0.000000e+002        0 0.0000000000 0.000000e+003        0 0.0000000000 0.000000e+004        0 0.0000000000 0.000000e+005        0 0.0000000000 0.000000e+006        0 0.0800369771 9.682457e-027        0 0.3016851330 2.165283e-018        0 0.0171385659 1.587435e-029        0 1.0794653385 6.116829e-0110       0 0.0000711943 7.095954e-05> fit$unit.df.adj[1:10,]   Sample1   Sample2   Sample31        0 0.0000000 0.00000002        0 0.0000000 0.00000003        0 0.0000000 0.00000004        0 0.0000000 0.00000005        0 0.0000000 0.00000006        0 0.4865584 0.48808377        0 0.4973603 0.50024828        0 0.4944987 0.49545889        0 0.4904124 0.491618410       0 0.4944965 0.4954567> summary(fit$deviance.adj)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 0.00000 0.08084 0.53006 1.07828 1.69219 5.39198 > summary(fit$df.residual.adj)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.0000  0.9757  0.9814  0.9358  0.9880  1.0941 > fit <- glmQLFit(y,design,legacy=TRUE)> summary(fit$dispersion)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  0.1958  0.1979  0.2098  0.2216  0.2441  0.2685 > summary(fit$s2.post)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 0.06182 0.60040 0.79878 0.91179 1.16628 2.59834 > > proc.time()   user  system elapsed    1.37    0.28    1.71