Yunshun Chen
| Snapshot Date: 2024-03-27 14:00:18 -0400 (Wed, 27 Mar 2024) |
| git_url: https://git.bioconductor.org/packages/edgeR |
| git_branch: devel |
| git_last_commit: 8cb63a6 |
| git_last_commit_date: 2024-03-23 07:21:34 -0400 (Sat, 23 Mar 2024) |
| Back to Multiple platform build/check report for BioC 3.19: simplified long |
|
This page was generated on 2024-03-28 11:40:40 -0400 (Thu, 28 Mar 2024).
| Hostname | OS | Arch (*) | R version | Installed pkgs |
|---|---|---|---|---|
| nebbiolo1 | Linux (Ubuntu 22.04.3 LTS) | x86_64 | R Under development (unstable) (2024-03-18 r86148) -- "Unsuffered Consequences" | 4708 |
| palomino3 | Windows Server 2022 Datacenter | x64 | R Under development (unstable) (2024-03-16 r86144 ucrt) -- "Unsuffered Consequences" | 4446 |
| lconway | macOS 12.7.1 Monterey | x86_64 | R Under development (unstable) (2024-03-18 r86148) -- "Unsuffered Consequences" | 4471 |
| kunpeng2 | Linux (openEuler 22.03 LTS-SP1) | aarch64 | R Under development (unstable) (2024-03-19 r86153) -- "Unsuffered Consequences" | 4426 |
| Click on any hostname to see more info about the system (e.g. compilers) (*) as reported by 'uname -p', except on Windows and Mac OS X | ||||
| Package 631/2270 | Hostname | OS / Arch | INSTALL | BUILD | CHECK | BUILD BIN | ||||||||
| edgeR 4.1.19 (landing page) Yunshun Chen
| nebbiolo1 | Linux (Ubuntu 22.04.3 LTS) / x86_64 | OK | OK | OK | | ||||||||
| palomino3 | Windows Server 2022 Datacenter / x64 | OK | OK | OK | OK | | ||||||||
| lconway | macOS 12.7.1 Monterey / x86_64 | OK | OK | OK | OK | | ||||||||
| kunpeng2 | Linux (openEuler 22.03 LTS-SP1) / aarch64 | OK | OK | OK | ||||||||||
|
To the developers/maintainers of the edgeR package: - Allow up to 24 hours (and sometimes 48 hours) for your latest push to git@git.bioconductor.org:packages/edgeR.git to reflect on this report. See Troubleshooting Build Report for more information. - Use the following Renviron settings to reproduce errors and warnings. - If 'R CMD check' started to fail recently on the Linux builder(s) over a missing dependency, add the missing dependency to 'Suggests:' in your DESCRIPTION file. See Renviron.bioc for more information. - See Martin Grigorov's blog post for how to debug Linux ARM64 related issues on a x86_64 host. |
| Package: edgeR |
| Version: 4.1.19 |
| Command: /home/biocbuild/R/R-4.4-devel-2024.03.20/bin/R CMD check --install=check:edgeR.install-out.txt --library=/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library --no-vignettes --timings edgeR_4.1.19.tar.gz |
| StartedAt: 2024-03-28 05:00:21 -0000 (Thu, 28 Mar 2024) |
| EndedAt: 2024-03-28 05:02:04 -0000 (Thu, 28 Mar 2024) |
| EllapsedTime: 102.2 seconds |
| RetCode: 0 |
| Status: OK |
| CheckDir: edgeR.Rcheck |
| Warnings: 0 |
##############################################################################
##############################################################################
###
### Running command:
###
### /home/biocbuild/R/R-4.4-devel-2024.03.20/bin/R CMD check --install=check:edgeR.install-out.txt --library=/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library --no-vignettes --timings edgeR_4.1.19.tar.gz
###
##############################################################################
##############################################################################
* using log directory ‘/home/biocbuild/bbs-3.19-bioc/meat/edgeR.Rcheck’
* using R Under development (unstable) (2024-03-19 r86153)
* using platform: aarch64-unknown-linux-gnu
* R was compiled by
gcc (GCC) 10.3.1
GNU Fortran (GCC) 10.3.1
* running under: openEuler 22.03 (LTS-SP1)
* using session charset: UTF-8
* using option ‘--no-vignettes’
* checking for file ‘edgeR/DESCRIPTION’ ... OK
* this is package ‘edgeR’ version ‘4.1.19’
* 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: ‘gcc (GCC) 10.3.1’
* used C++ compiler: ‘g++ (GCC) 10.3.1’
* checking installed package size ... NOTE
installed size is 11.3Mb
sub-directories of 1Mb or more:
libs 10.1Mb
* 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 loading without being on the library search path ... 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 ... NOTE
checkRd: (-1) equalizeLibSizes.Rd:29: Lost braces; missing escapes or markup?
29 | \code{equalizeLibSizes.default} returns a list with components {pseudo.counts} and {pseudo.lib.size}.
| ^
checkRd: (-1) equalizeLibSizes.Rd:29: Lost braces; missing escapes or markup?
29 | \code{equalizeLibSizes.default} returns a list with components {pseudo.counts} and {pseudo.lib.size}.
| ^
checkRd: (-1) maPlot.Rd:35: Lost braces; missing escapes or markup?
35 | \value{a plot to the current device (if \code{plot.it=TRUE}), and invisibly returns the \code{M} (logFC) and \code{A} (logConc) values used for the plot, plus identifiers \code{w} and \code{v} of genes for which \code{M} and {A} values, or just \code{M} values, respectively, were adjusted to make a nicer looking plot. }
| ^
checkRd: (-1) meanvar.Rd:47: Lost braces; missing escapes or markup?
47 | \item{group}{factor giving the experimental group or condition to which each sample (i.e. column of \code{x} or element of {y}) belongs}
| ^
checkRd: (-1) nbinomUnitDeviance.Rd:24: Lost braces; missing escapes or markup?
24 | The unit deviance of the negative binomial distribution is a measure of the distance between \code{y} and {mean}.
| ^
* checking Rd metadata ... OK
* checking Rd cross-references ... OK
* 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 ... NOTE
Note: information on .o files is not available
* checking files in ‘vignettes’ ... OK
* checking examples ... OK
Examples with CPU (user + system) or elapsed time > 5s
user system elapsed
nearestTSS 6.682 0.228 6.928
romer.DGEList 5.150 0.024 5.184
romer.DGEGLM 5.144 0.024 5.178
* 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
* DONE
Status: 3 NOTEs
See
‘/home/biocbuild/bbs-3.19-bioc/meat/edgeR.Rcheck/00check.log’
for details.
edgeR.Rcheck/00install.out
############################################################################## ############################################################################## ### ### Running command: ### ### /home/biocbuild/R/R-4.4-devel-2024.03.20/bin/R CMD INSTALL edgeR ### ############################################################################## ############################################################################## * installing to library ‘/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library’ * installing *source* package ‘edgeR’ ... ** using staged installation ** libs using C compiler: ‘gcc (GCC) 10.3.1’ using C++ compiler: ‘g++ (GCC) 10.3.1’ g++ -std=gnu++17 -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c R_add_prior_count.cpp -o R_add_prior_count.o g++ -std=gnu++17 -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c R_ave_log_cpm.cpp -o R_ave_log_cpm.o g++ -std=gnu++17 -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c R_calculate_cpm.cpp -o R_calculate_cpm.o g++ -std=gnu++17 -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c R_check_poisson_bound.cpp -o R_check_poisson_bound.o g++ -std=gnu++17 -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c R_compute_apl.cpp -o R_compute_apl.o g++ -std=gnu++17 -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c R_compute_nbdev.cpp -o R_compute_nbdev.o g++ -std=gnu++17 -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c R_exact_test_by_deviance.cpp -o R_exact_test_by_deviance.o g++ -std=gnu++17 -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c R_fit_levenberg.cpp -o R_fit_levenberg.o g++ -std=gnu++17 -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c R_fit_one_group.cpp -o R_fit_one_group.o g++ -std=gnu++17 -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c R_get_one_way_fitted.cpp -o R_get_one_way_fitted.o g++ -std=gnu++17 -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c R_initialize_levenberg.cpp -o R_initialize_levenberg.o g++ -std=gnu++17 -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c R_loess_by_col.cpp -o R_loess_by_col.o g++ -std=gnu++17 -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c R_maximize_interpolant.cpp -o R_maximize_interpolant.o gcc -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c R_process_hairpin_reads.c -o R_process_hairpin_reads.o gcc -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c R_ql_glm.c -o R_ql_glm.o g++ -std=gnu++17 -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c R_simple_good_turing.cpp -o R_simple_good_turing.o g++ -std=gnu++17 -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c add_prior.cpp -o add_prior.o g++ -std=gnu++17 -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c adj_coxreid.cpp -o adj_coxreid.o gcc -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c fmm_spline.c -o fmm_spline.o g++ -std=gnu++17 -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c glm_levenberg.cpp -o glm_levenberg.o g++ -std=gnu++17 -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c glm_one_group.cpp -o glm_one_group.o g++ -std=gnu++17 -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c init.cpp -o init.o g++ -std=gnu++17 -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c interpolator.cpp -o interpolator.o gcc -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c nbdev.c -o nbdev.o g++ -std=gnu++17 -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c objects.cpp -o objects.o gcc -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c ql_glm.c -o ql_glm.o gcc -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c ql_weights.c -o ql_weights.o gcc -I"/home/biocbuild/R/R-4.4-devel-2024.03.20/include" -DNDEBUG -I'/home/biocbuild/R/R-4.4-devel-2024.03.20/site-library/Rcpp/include' -I/usr/local/include -fPIC -g -O2 -Wall -c qr_hat.c -o qr_hat.o g++ -std=gnu++17 -shared -L/home/biocbuild/R/R-4.4-devel-2024.03.20/lib -L/usr/local/lib -o edgeR.so R_add_prior_count.o R_ave_log_cpm.o R_calculate_cpm.o R_check_poisson_bound.o R_compute_apl.o R_compute_nbdev.o R_exact_test_by_deviance.o R_fit_levenberg.o R_fit_one_group.o R_get_one_way_fitted.o R_initialize_levenberg.o R_loess_by_col.o R_maximize_interpolant.o R_process_hairpin_reads.o R_ql_glm.o R_simple_good_turing.o add_prior.o adj_coxreid.o fmm_spline.o glm_levenberg.o glm_one_group.o init.o interpolator.o nbdev.o objects.o ql_glm.o ql_weights.o qr_hat.o -lopenblas -lgfortran -lm -L/home/biocbuild/R/R-4.4-devel-2024.03.20/lib -lR installing to /home/biocbuild/R/R-4.4-devel-2024.03.20/site-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)
edgeR.Rcheck/tests/edgeR-Tests.Rout
R Under development (unstable) (2024-03-19 r86153) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: aarch64-unknown-linux-gnu
R 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 FDR
Tag.17 2.0450964 13.73726 0.01975954 0.4347099
Tag.21 -1.7265870 13.38327 0.06131012 0.6744114
Tag.6 -1.6329986 12.81479 0.12446044 0.8982100
Tag.2 4.0861092 11.54121 0.16331090 0.8982100
Tag.16 0.9324996 13.57074 0.29050785 0.9655885
Tag.20 0.8543138 13.76364 0.31736609 0.9655885
Tag.12 0.7081170 14.31389 0.37271028 0.9655885
Tag.19 -0.7976602 13.31405 0.40166354 0.9655885
Tag.3 -0.7300410 13.54155 0.42139935 0.9655885
Tag.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 FDR
Tag.17 2.0450964 13.73726 0.01975954 0.4347099
Tag.21 -1.7265870 13.38327 0.06131012 0.6744114
Tag.6 -1.6329986 12.81479 0.12446044 0.8982100
Tag.2 4.0861092 11.54121 0.16331090 0.8982100
Tag.16 0.9324996 13.57074 0.29050785 0.9655885
Tag.20 0.8543138 13.76364 0.31736609 0.9655885
Tag.12 0.7081170 14.31389 0.37271028 0.9655885
Tag.19 -0.7976602 13.31405 0.40166354 0.9655885
Tag.3 -0.7300410 13.54155 0.42139935 0.9655885
Tag.8 -0.7917906 12.86353 0.47117217 0.9655885
>
> de <- exactTest(d2)
> topTags(de)
Comparison of groups: 2-1
logFC logCPM PValue FDR
Tag.17 2.0450987 13.73726 0.01327001 0.2919403
Tag.21 -1.7265897 13.38327 0.05683886 0.6252275
Tag.6 -1.6329910 12.81479 0.11460208 0.8404152
Tag.2 4.0861092 11.54121 0.16126207 0.8869414
Tag.16 0.9324975 13.57074 0.28103256 0.9669238
Tag.20 0.8543178 13.76364 0.30234789 0.9669238
Tag.12 0.7081149 14.31389 0.37917895 0.9669238
Tag.19 -0.7976633 13.31405 0.40762735 0.9669238
Tag.3 -0.7300478 13.54155 0.40856822 0.9669238
Tag.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 FDR
Tag.17 2.0450951 13.73726 0.02427872 0.5341319
Tag.21 -1.7265927 13.38327 0.05234833 0.5758316
Tag.6 -1.6330014 12.81479 0.12846308 0.8954397
Tag.2 4.0861092 11.54121 0.16280722 0.8954397
Tag.16 0.9324887 13.57074 0.24308201 0.9711975
Tag.20 0.8543044 13.76364 0.35534649 0.9711975
Tag.19 -0.7976535 13.31405 0.38873717 0.9711975
Tag.3 -0.7300525 13.54155 0.40001438 0.9711975
Tag.12 0.7080985 14.31389 0.43530227 0.9711975
Tag.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 FDR
Tag.17 2.0450979 13.73726 0.01546795 0.3402949
Tag.21 -1.7266049 13.38327 0.03545446 0.3899990
Tag.6 -1.6329841 12.81479 0.10632987 0.7797524
Tag.2 4.0861092 11.54121 0.16057893 0.8831841
Tag.16 0.9324935 13.57074 0.26348818 0.9658389
Tag.20 0.8543140 13.76364 0.31674090 0.9658389
Tag.19 -0.7976354 13.31405 0.35564858 0.9658389
Tag.3 -0.7300593 13.54155 0.38833737 0.9658389
Tag.12 0.7081041 14.31389 0.41513004 0.9658389
Tag.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 FDR
Tag.17 2.0450979 13.73726 0.01546795 0.3402949
Tag.21 -1.7266049 13.38327 0.03545446 0.3899990
Tag.6 -1.6329841 12.81479 0.10632987 0.7797524
Tag.2 4.0861092 11.54121 0.16057893 0.8831841
Tag.16 0.9324935 13.57074 0.26348818 0.9658389
Tag.20 0.8543140 13.76364 0.31674090 0.9658389
Tag.19 -0.7976354 13.31405 0.35564858 0.9658389
Tag.3 -0.7300593 13.54155 0.38833737 0.9658389
Tag.12 0.7081041 14.31389 0.41513004 0.9658389
Tag.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) group2
Tag.1 -1.000000e+08 0.000000e+00
Tag.2 -1.000000e+08 1.000000e+08
Tag.3 2.525729e+00 -5.108256e-01
Tag.4 2.525729e+00 1.484200e-01
Tag.5 2.140066e+00 -1.941560e-01
Tag.6 2.079442e+00 -1.163151e+00
Tag.7 2.014903e+00 2.363888e-01
Tag.8 1.945910e+00 -5.596158e-01
Tag.9 1.504077e+00 2.006707e-01
Tag.10 2.302585e+00 2.623643e-01
$fitted.values
Sample1 Sample2 Sample3 Sample4
Tag.1 0.0 0.0 0.0 0.0
Tag.2 0.0 0.0 2.0 2.0
Tag.3 12.5 12.5 7.5 7.5
Tag.4 12.5 12.5 14.5 14.5
Tag.5 8.5 8.5 7.0 7.0
Tag.6 8.0 8.0 2.5 2.5
Tag.7 7.5 7.5 9.5 9.5
Tag.8 7.0 7.0 4.0 4.0
Tag.9 4.5 4.5 5.5 5.5
Tag.10 10.0 10.0 13.0 13.0
> mglmOneWay(d[1:10,],design,dispersion=0)
$coefficients
(Intercept) group2
Tag.1 -1.000000e+08 0.000000e+00
Tag.2 -1.000000e+08 1.000000e+08
Tag.3 2.525729e+00 -5.108256e-01
Tag.4 2.525729e+00 1.484200e-01
Tag.5 2.140066e+00 -1.941560e-01
Tag.6 2.079442e+00 -1.163151e+00
Tag.7 2.014903e+00 2.363888e-01
Tag.8 1.945910e+00 -5.596158e-01
Tag.9 1.504077e+00 2.006707e-01
Tag.10 2.302585e+00 2.623643e-01
$fitted.values
Sample1 Sample2 Sample3 Sample4
Tag.1 0.0 0.0 0.0 0.0
Tag.2 0.0 0.0 2.0 2.0
Tag.3 12.5 12.5 7.5 7.5
Tag.4 12.5 12.5 14.5 14.5
Tag.5 8.5 8.5 7.0 7.0
Tag.6 8.0 8.0 2.5 2.5
Tag.7 7.5 7.5 9.5 9.5
Tag.8 7.0 7.0 4.0 4.0
Tag.9 4.5 4.5 5.5 5.5
Tag.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 FDR
Tag.17 2.0450964 13.73726 6.0485417 0.01391779 0.3058698
Tag.2 4.0861092 11.54121 4.8400340 0.02780635 0.3058698
Tag.21 -1.7265870 13.38327 4.0777825 0.04345065 0.3186381
Tag.6 -1.6329986 12.81479 3.0078205 0.08286364 0.4557500
Tag.16 0.9324996 13.57074 1.3477682 0.24566867 0.8276702
Tag.20 0.8543138 13.76364 1.1890032 0.27553071 0.8276702
Tag.19 -0.7976602 13.31405 0.9279151 0.33540526 0.8276702
Tag.12 0.7081170 14.31389 0.9095513 0.34023349 0.8276702
Tag.3 -0.7300410 13.54155 0.8300307 0.36226364 0.8276702
Tag.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 FDR
Tag.17 2.0450964 13.73726 6.0485417 0.01391779 0.3058698
Tag.2 4.0861092 11.54121 4.8400340 0.02780635 0.3058698
Tag.21 -1.7265870 13.38327 4.0777825 0.04345065 0.3186381
Tag.6 -1.6329986 12.81479 3.0078205 0.08286364 0.4557500
Tag.16 0.9324996 13.57074 1.3477682 0.24566867 0.8276702
Tag.20 0.8543138 13.76364 1.1890032 0.27553071 0.8276702
Tag.19 -0.7976602 13.31405 0.9279151 0.33540526 0.8276702
Tag.12 0.7081170 14.31389 0.9095513 0.34023349 0.8276702
Tag.3 -0.7300410 13.54155 0.8300307 0.36226364 0.8276702
Tag.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 FDR
Tag.17 2.0450988 13.73727 6.8001118 0.009115216 0.2005348
Tag.2 4.0861092 11.54122 4.8594088 0.027495756 0.2872068
Tag.21 -1.7265904 13.38327 4.2537154 0.039164558 0.2872068
Tag.6 -1.6329904 12.81479 3.1763761 0.074710253 0.4109064
Tag.16 0.9324970 13.57074 1.4126709 0.234613512 0.8499599
Tag.20 0.8543183 13.76364 1.2721097 0.259371274 0.8499599
Tag.19 -0.7976614 13.31405 0.9190392 0.337727381 0.8499599
Tag.12 0.7081163 14.31389 0.9014515 0.342392806 0.8499599
Tag.3 -0.7300488 13.54155 0.8817937 0.347710872 0.8499599
Tag.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 <- calcNormFactors(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 FDR
Gene1 2.907024 13.56183 38.738512 4.845536e-10 4.845536e-07
Gene61 2.855317 10.27136 10.738307 1.049403e-03 5.247015e-01
Gene62 -2.123902 10.53174 8.818703 2.981585e-03 8.334760e-01
Gene134 -1.949073 10.53355 8.125889 4.363759e-03 8.334760e-01
Gene740 -1.610046 10.94907 8.013408 4.643227e-03 8.334760e-01
Gene354 2.022698 10.45066 7.826308 5.149118e-03 8.334760e-01
Gene5 1.856816 10.45249 7.214238 7.232750e-03 8.334760e-01
Gene746 -1.798331 10.53094 6.846262 8.882693e-03 8.334760e-01
Gene110 1.623148 10.68607 6.737984 9.438120e-03 8.334760e-01
Gene383 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) x
Gene1 -7.391745 2.0149958
Gene2 -7.318483 -0.7611895
Gene3 -6.831702 -0.1399478
Gene4 -7.480255 0.5172002
Gene5 -8.747793 1.2870467
995 more rows ...
$fitted.values
x0 x1 x2
Gene1 2.3570471 18.954454 138.2791328
Gene2 2.5138172 1.089292 0.4282107
Gene3 4.1580452 3.750528 3.0690081
Gene4 2.1012460 3.769592 6.1349937
Gene5 0.5080377 2.136398 8.1502486
995 more rows ...
$deviance
[1] 6.38037545 1.46644913 1.38532340 0.01593969 1.03894513
995 more elements ...
$iter
[1] 8 4 4 4 6
995 more elements ...
$failed
[1] FALSE FALSE FALSE FALSE FALSE
995 more elements ...
$method
[1] "levenberg"
$counts
x0 x1 x2
Gene1 0 30 110
Gene2 2 2 0
Gene3 3 6 2
Gene4 2 4 6
Gene5 1 1 9
995 more rows ...
$unshrunk.coefficients
(Intercept) x
Gene1 -7.437763 2.0412762
Gene2 -7.373370 -0.8796273
Gene3 -6.870127 -0.1465014
Gene4 -7.552642 0.5410832
Gene5 -8.972372 1.3929679
995 more rows ...
$df.residual
[1] 1 1 1 1 1
995 more elements ...
$design
(Intercept) x
1 1 0
2 1 1
3 1 2
attr(,"assign")
[1] 0 1
$offset
[,1] [,2] [,3]
[1,] 8.295172 8.338525 8.284484
attr(,"class")
[1] "CompressedMatrix"
attr(,"Dims")
[1] 5 3
attr(,"repeat.row")
[1] TRUE
attr(,"repeat.col")
[1] FALSE
995 more rows ...
$dispersion
[1] 0.1
$prior.count
[1] 0.1666667
$samples
group lib.size norm.factors
x0 1 4001 1.0008730
x1 1 4176 1.0014172
x2 1 3971 0.9977138
$AveLogCPM
[1] 13.561832 9.682757 10.447014 10.532113 10.452489
995 more elements ...
>
> calcNormFactors(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.05545 0.09511 0.11623 0.11014 0.13329 0.16861
> d2 <- estimateDisp(d, design, prior.df=20)
> summary(d2$tagwise.dispersion)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.04203 0.08586 0.11280 0.11010 0.12369 0.37408
> d2 <- estimateDisp(d, design, robust=TRUE)
> summary(d2$tagwise.dispersion)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.05545 0.09511 0.11623 0.11014 0.13329 0.16861
>
> # 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 FDR
1 -7.2381060 -0.0621100 17.19071 10.7712171 0.004582051 0.02291026
5 -1.6684268 -0.9326507 17.33529 1.7309951 0.420842115 0.90967967
2 1.2080938 1.0420198 18.24544 1.0496688 0.591653347 0.90967967
4 0.5416704 -0.1506381 17.57744 0.3958596 0.820427427 0.90967967
3 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 FDR
1 -7.2381060 7.1759960 -0.0621100 17.19071 10.7712171 0.004582051 0.02291026
5 -1.6684268 0.7357761 -0.9326507 17.33529 1.7309951 0.420842115 0.90967967
2 1.2080938 -0.1660740 1.0420198 18.24544 1.0496688 0.591653347 0.90967967
4 0.5416704 -0.6923084 -0.1506381 17.57744 0.3958596 0.820427427 0.90967967
3 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)
> 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.18263
> y$common.dispersion
[1] 0.2185181
> y$trended.dispersion[1:10]
[1] 0.2717585 0.2727238 0.2717585 0.2690969 0.2617797 0.2634692 0.2453099
[8] 0.1931524 0.1930463 0.1928675
> summary(y$trended.dispersion)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.1928 0.1932 0.2086 0.2226 0.2515 0.2728
> y$tagwise.dispersion[1:10]
[1] 0.2717585 0.2727238 0.2717585 0.2690969 0.2617797 0.1519257 0.1936500
[8] 0.1015399 0.2475146 0.1013302
> summary(y$tagwise.dispersion)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.09719 0.13460 0.19422 0.22321 0.27176 0.61463
>
> # glmQLFit
> fit <- glmQLFit(y,design,legacy=FALSE)
> fit$dispersion
[1] 0.2073113
> summary(fit$var.post)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.007627 0.024804 0.032000 0.143246 0.059155 2.589922
> fit$unit.deviance.adj[1:10,]
Sample1 Sample2 Sample3
1 0 0.000000e+00 0.000000e+00
2 0 0.000000e+00 0.000000e+00
3 0 0.000000e+00 0.000000e+00
4 0 0.000000e+00 0.000000e+00
5 0 0.000000e+00 0.000000e+00
6 0 7.715120e-02 9.346359e-02
7 0 2.951426e-01 2.111539e-01
8 0 1.640577e-02 1.518755e-02
9 0 1.039516e+00 5.865753e-01
10 0 6.813594e-05 6.790318e-05
> fit$unit.df.adj[1:10,]
Sample1 Sample2 Sample3
1 0 0.0000000 0.0000000
2 0 0.0000000 0.0000000
3 0 0.0000000 0.0000000
4 0 0.0000000 0.0000000
5 0 0.0000000 0.0000000
6 0 0.4867038 0.4881894
7 0 0.4992433 0.5021103
8 0 0.4942884 0.4952040
9 0 0.4902472 0.4914076
10 0 0.4942861 0.4952019
> summary(fit$deviance.adj)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.0778 0.5146 1.0412 1.6254 5.1916
> summary(fit$df.residual.adj)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.9759 0.9817 0.9363 0.9881 1.1005
>
> proc.time()
user system elapsed
4.599 0.089 4.683
edgeR.Rcheck/tests/edgeR-Tests.Rout.save
R version 4.3.2 (2023-10-31 ucrt) -- "Eye Holes"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)
R 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 FDR
Tag.17 2.0450964 13.73726 0.01975954 0.4347099
Tag.21 -1.7265870 13.38327 0.06131012 0.6744114
Tag.6 -1.6329986 12.81479 0.12446044 0.8982100
Tag.2 4.0861092 11.54121 0.16331090 0.8982100
Tag.16 0.9324996 13.57074 0.29050785 0.9655885
Tag.20 0.8543138 13.76364 0.31736609 0.9655885
Tag.12 0.7081170 14.31389 0.37271028 0.9655885
Tag.19 -0.7976602 13.31405 0.40166354 0.9655885
Tag.3 -0.7300410 13.54155 0.42139935 0.9655885
Tag.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 FDR
Tag.17 2.0450964 13.73726 0.01975954 0.4347099
Tag.21 -1.7265870 13.38327 0.06131012 0.6744114
Tag.6 -1.6329986 12.81479 0.12446044 0.8982100
Tag.2 4.0861092 11.54121 0.16331090 0.8982100
Tag.16 0.9324996 13.57074 0.29050785 0.9655885
Tag.20 0.8543138 13.76364 0.31736609 0.9655885
Tag.12 0.7081170 14.31389 0.37271028 0.9655885
Tag.19 -0.7976602 13.31405 0.40166354 0.9655885
Tag.3 -0.7300410 13.54155 0.42139935 0.9655885
Tag.8 -0.7917906 12.86353 0.47117217 0.9655885
>
> de <- exactTest(d2)
> topTags(de)
Comparison of groups: 2-1
logFC logCPM PValue FDR
Tag.17 2.0450987 13.73726 0.01327001 0.2919403
Tag.21 -1.7265897 13.38327 0.05683886 0.6252275
Tag.6 -1.6329910 12.81479 0.11460208 0.8404152
Tag.2 4.0861092 11.54121 0.16126207 0.8869414
Tag.16 0.9324975 13.57074 0.28103256 0.9669238
Tag.20 0.8543178 13.76364 0.30234789 0.9669238
Tag.12 0.7081149 14.31389 0.37917895 0.9669238
Tag.19 -0.7976633 13.31405 0.40762735 0.9669238
Tag.3 -0.7300478 13.54155 0.40856822 0.9669238
Tag.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 FDR
Tag.17 2.0450951 13.73726 0.02427872 0.5341319
Tag.21 -1.7265927 13.38327 0.05234833 0.5758316
Tag.6 -1.6330014 12.81479 0.12846308 0.8954397
Tag.2 4.0861092 11.54121 0.16280722 0.8954397
Tag.16 0.9324887 13.57074 0.24308201 0.9711975
Tag.20 0.8543044 13.76364 0.35534649 0.9711975
Tag.19 -0.7976535 13.31405 0.38873717 0.9711975
Tag.3 -0.7300525 13.54155 0.40001438 0.9711975
Tag.12 0.7080985 14.31389 0.43530227 0.9711975
Tag.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 FDR
Tag.17 2.0450979 13.73726 0.01546795 0.3402949
Tag.21 -1.7266049 13.38327 0.03545446 0.3899990
Tag.6 -1.6329841 12.81479 0.10632987 0.7797524
Tag.2 4.0861092 11.54121 0.16057893 0.8831841
Tag.16 0.9324935 13.57074 0.26348818 0.9658389
Tag.20 0.8543140 13.76364 0.31674090 0.9658389
Tag.19 -0.7976354 13.31405 0.35564858 0.9658389
Tag.3 -0.7300593 13.54155 0.38833737 0.9658389
Tag.12 0.7081041 14.31389 0.41513004 0.9658389
Tag.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 FDR
Tag.17 2.0450979 13.73726 0.01546795 0.3402949
Tag.21 -1.7266049 13.38327 0.03545446 0.3899990
Tag.6 -1.6329841 12.81479 0.10632987 0.7797524
Tag.2 4.0861092 11.54121 0.16057893 0.8831841
Tag.16 0.9324935 13.57074 0.26348818 0.9658389
Tag.20 0.8543140 13.76364 0.31674090 0.9658389
Tag.19 -0.7976354 13.31405 0.35564858 0.9658389
Tag.3 -0.7300593 13.54155 0.38833737 0.9658389
Tag.12 0.7081041 14.31389 0.41513004 0.9658389
Tag.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) group2
Tag.1 -1.000000e+08 0.000000e+00
Tag.2 -1.000000e+08 1.000000e+08
Tag.3 2.525729e+00 -5.108256e-01
Tag.4 2.525729e+00 1.484200e-01
Tag.5 2.140066e+00 -1.941560e-01
Tag.6 2.079442e+00 -1.163151e+00
Tag.7 2.014903e+00 2.363888e-01
Tag.8 1.945910e+00 -5.596158e-01
Tag.9 1.504077e+00 2.006707e-01
Tag.10 2.302585e+00 2.623643e-01
$fitted.values
Sample1 Sample2 Sample3 Sample4
Tag.1 0.0 0.0 0.0 0.0
Tag.2 0.0 0.0 2.0 2.0
Tag.3 12.5 12.5 7.5 7.5
Tag.4 12.5 12.5 14.5 14.5
Tag.5 8.5 8.5 7.0 7.0
Tag.6 8.0 8.0 2.5 2.5
Tag.7 7.5 7.5 9.5 9.5
Tag.8 7.0 7.0 4.0 4.0
Tag.9 4.5 4.5 5.5 5.5
Tag.10 10.0 10.0 13.0 13.0
> mglmOneWay(d[1:10,],design,dispersion=0)
$coefficients
(Intercept) group2
Tag.1 -1.000000e+08 0.000000e+00
Tag.2 -1.000000e+08 1.000000e+08
Tag.3 2.525729e+00 -5.108256e-01
Tag.4 2.525729e+00 1.484200e-01
Tag.5 2.140066e+00 -1.941560e-01
Tag.6 2.079442e+00 -1.163151e+00
Tag.7 2.014903e+00 2.363888e-01
Tag.8 1.945910e+00 -5.596158e-01
Tag.9 1.504077e+00 2.006707e-01
Tag.10 2.302585e+00 2.623643e-01
$fitted.values
Sample1 Sample2 Sample3 Sample4
Tag.1 0.0 0.0 0.0 0.0
Tag.2 0.0 0.0 2.0 2.0
Tag.3 12.5 12.5 7.5 7.5
Tag.4 12.5 12.5 14.5 14.5
Tag.5 8.5 8.5 7.0 7.0
Tag.6 8.0 8.0 2.5 2.5
Tag.7 7.5 7.5 9.5 9.5
Tag.8 7.0 7.0 4.0 4.0
Tag.9 4.5 4.5 5.5 5.5
Tag.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 FDR
Tag.17 2.0450964 13.73726 6.0485417 0.01391779 0.3058698
Tag.2 4.0861092 11.54121 4.8400340 0.02780635 0.3058698
Tag.21 -1.7265870 13.38327 4.0777825 0.04345065 0.3186381
Tag.6 -1.6329986 12.81479 3.0078205 0.08286364 0.4557500
Tag.16 0.9324996 13.57074 1.3477682 0.24566867 0.8276702
Tag.20 0.8543138 13.76364 1.1890032 0.27553071 0.8276702
Tag.19 -0.7976602 13.31405 0.9279151 0.33540526 0.8276702
Tag.12 0.7081170 14.31389 0.9095513 0.34023349 0.8276702
Tag.3 -0.7300410 13.54155 0.8300307 0.36226364 0.8276702
Tag.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 FDR
Tag.17 2.0450964 13.73726 6.0485417 0.01391779 0.3058698
Tag.2 4.0861092 11.54121 4.8400340 0.02780635 0.3058698
Tag.21 -1.7265870 13.38327 4.0777825 0.04345065 0.3186381
Tag.6 -1.6329986 12.81479 3.0078205 0.08286364 0.4557500
Tag.16 0.9324996 13.57074 1.3477682 0.24566867 0.8276702
Tag.20 0.8543138 13.76364 1.1890032 0.27553071 0.8276702
Tag.19 -0.7976602 13.31405 0.9279151 0.33540526 0.8276702
Tag.12 0.7081170 14.31389 0.9095513 0.34023349 0.8276702
Tag.3 -0.7300410 13.54155 0.8300307 0.36226364 0.8276702
Tag.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 FDR
Tag.17 2.0450988 13.73727 6.8001118 0.009115216 0.2005348
Tag.2 4.0861092 11.54122 4.8594088 0.027495756 0.2872068
Tag.21 -1.7265904 13.38327 4.2537154 0.039164558 0.2872068
Tag.6 -1.6329904 12.81479 3.1763761 0.074710253 0.4109064
Tag.16 0.9324970 13.57074 1.4126709 0.234613512 0.8499599
Tag.20 0.8543183 13.76364 1.2721097 0.259371274 0.8499599
Tag.19 -0.7976614 13.31405 0.9190392 0.337727381 0.8499599
Tag.12 0.7081163 14.31389 0.9014515 0.342392806 0.8499599
Tag.3 -0.7300488 13.54155 0.8817937 0.347710872 0.8499599
Tag.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 <- calcNormFactors(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 FDR
Gene1 2.907024 13.56183 38.738512 4.845536e-10 4.845536e-07
Gene61 2.855317 10.27136 10.738307 1.049403e-03 5.247015e-01
Gene62 -2.123902 10.53174 8.818703 2.981585e-03 8.334760e-01
Gene134 -1.949073 10.53355 8.125889 4.363759e-03 8.334760e-01
Gene740 -1.610046 10.94907 8.013408 4.643227e-03 8.334760e-01
Gene354 2.022698 10.45066 7.826308 5.149118e-03 8.334760e-01
Gene5 1.856816 10.45249 7.214238 7.232750e-03 8.334760e-01
Gene746 -1.798331 10.53094 6.846262 8.882693e-03 8.334760e-01
Gene110 1.623148 10.68607 6.737984 9.438120e-03 8.334760e-01
Gene383 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) x
Gene1 -7.391745 2.0149958
Gene2 -7.318483 -0.7611895
Gene3 -6.831702 -0.1399478
Gene4 -7.480255 0.5172002
Gene5 -8.747793 1.2870467
995 more rows ...
$fitted.values
x0 x1 x2
Gene1 2.3570471 18.954454 138.2791328
Gene2 2.5138172 1.089292 0.4282107
Gene3 4.1580452 3.750528 3.0690081
Gene4 2.1012460 3.769592 6.1349937
Gene5 0.5080377 2.136398 8.1502486
995 more rows ...
$deviance
[1] 6.38037545 1.46644913 1.38532340 0.01593969 1.03894513
995 more elements ...
$iter
[1] 8 4 4 4 6
995 more elements ...
$failed
[1] FALSE FALSE FALSE FALSE FALSE
995 more elements ...
$method
[1] "levenberg"
$counts
x0 x1 x2
Gene1 0 30 110
Gene2 2 2 0
Gene3 3 6 2
Gene4 2 4 6
Gene5 1 1 9
995 more rows ...
$unshrunk.coefficients
(Intercept) x
Gene1 -7.437763 2.0412762
Gene2 -7.373370 -0.8796273
Gene3 -6.870127 -0.1465014
Gene4 -7.552642 0.5410832
Gene5 -8.972372 1.3929679
995 more rows ...
$df.residual
[1] 1 1 1 1 1
995 more elements ...
$design
(Intercept) x
1 1 0
2 1 1
3 1 2
attr(,"assign")
[1] 0 1
$offset
[,1] [,2] [,3]
[1,] 8.295172 8.338525 8.284484
attr(,"class")
[1] "CompressedMatrix"
attr(,"Dims")
[1] 5 3
attr(,"repeat.row")
[1] TRUE
attr(,"repeat.col")
[1] FALSE
995 more rows ...
$dispersion
[1] 0.1
$prior.count
[1] 0.1666667
$samples
group lib.size norm.factors
x0 1 4001 1.0008730
x1 1 4176 1.0014172
x2 1 3971 0.9977138
$AveLogCPM
[1] 13.561832 9.682757 10.447014 10.532113 10.452489
995 more elements ...
>
> calcNormFactors(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.05545 0.09511 0.11623 0.11014 0.13329 0.16861
> d2 <- estimateDisp(d, design, prior.df=20)
> summary(d2$tagwise.dispersion)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.04203 0.08586 0.11280 0.11010 0.12369 0.37408
> d2 <- estimateDisp(d, design, robust=TRUE)
> summary(d2$tagwise.dispersion)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.05545 0.09511 0.11623 0.11014 0.13329 0.16861
>
> # 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 FDR
1 -7.2381060 -0.0621100 17.19071 10.7712171 0.004582051 0.02291026
5 -1.6684268 -0.9326507 17.33529 1.7309951 0.420842115 0.90967967
2 1.2080938 1.0420198 18.24544 1.0496688 0.591653347 0.90967967
4 0.5416704 -0.1506381 17.57744 0.3958596 0.820427427 0.90967967
3 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 FDR
1 -7.2381060 7.1759960 -0.0621100 17.19071 10.7712171 0.004582051 0.02291026
5 -1.6684268 0.7357761 -0.9326507 17.33529 1.7309951 0.420842115 0.90967967
2 1.2080938 -0.1660740 1.0420198 18.24544 1.0496688 0.591653347 0.90967967
4 0.5416704 -0.6923084 -0.1506381 17.57744 0.3958596 0.820427427 0.90967967
3 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)
> 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.18263
> y$common.dispersion
[1] 0.2185181
> y$trended.dispersion[1:10]
[1] 0.2717585 0.2727238 0.2717585 0.2690969 0.2617797 0.2634692 0.2453099
[8] 0.1931524 0.1930463 0.1928675
> summary(y$trended.dispersion)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.1928 0.1932 0.2086 0.2226 0.2515 0.2728
> y$tagwise.dispersion[1:10]
[1] 0.2717585 0.2727238 0.2717585 0.2690969 0.2617797 0.1519257 0.1936500
[8] 0.1015399 0.2475146 0.1013302
> summary(y$tagwise.dispersion)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.09719 0.13460 0.19422 0.22321 0.27176 0.61463
>
> # glmQLFit
> fit <- glmQLFit(y,design,legacy=FALSE)
> fit$dispersion
[1] 0.2073113
> summary(fit$var.post)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.007627 0.024804 0.032000 0.143246 0.059155 2.589922
> fit$unit.deviance.adj[1:10,]
Sample1 Sample2 Sample3
1 0 0.000000e+00 0.000000e+00
2 0 0.000000e+00 0.000000e+00
3 0 0.000000e+00 0.000000e+00
4 0 0.000000e+00 0.000000e+00
5 0 0.000000e+00 0.000000e+00
6 0 7.715120e-02 9.346359e-02
7 0 2.951426e-01 2.111539e-01
8 0 1.640577e-02 1.518755e-02
9 0 1.039516e+00 5.865753e-01
10 0 6.813594e-05 6.790318e-05
> fit$unit.df.adj[1:10,]
Sample1 Sample2 Sample3
1 0 0.0000000 0.0000000
2 0 0.0000000 0.0000000
3 0 0.0000000 0.0000000
4 0 0.0000000 0.0000000
5 0 0.0000000 0.0000000
6 0 0.4867038 0.4881894
7 0 0.4992433 0.5021103
8 0 0.4942884 0.4952040
9 0 0.4902472 0.4914076
10 0 0.4942861 0.4952019
> summary(fit$deviance.adj)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.0778 0.5146 1.0412 1.6254 5.1916
> summary(fit$df.residual.adj)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.9759 0.9817 0.9363 0.9881 1.1005
>
> proc.time()
user system elapsed
2.07 0.10 2.18
edgeR.Rcheck/edgeR-Ex.timings
| name | user | system | elapsed | |
| DGEList | 0.018 | 0.000 | 0.018 | |
| SE2DGEList | 0 | 0 | 0 | |
| Seurat2PB | 0 | 0 | 0 | |
| WLEB | 0.033 | 0.000 | 0.032 | |
| addPriorCount | 0.002 | 0.000 | 0.002 | |
| adjustedProfileLik | 0.006 | 0.000 | 0.006 | |
| aveLogCPM | 0.002 | 0.000 | 0.002 | |
| binomTest | 0.001 | 0.000 | 0.001 | |
| calcNormFactors | 0.01 | 0.00 | 0.01 | |
| camera.DGEGLM | 0.062 | 0.008 | 0.070 | |
| camera.DGEList | 0.123 | 0.008 | 0.131 | |
| catchSalmon | 0 | 0 | 0 | |
| cbind | 0 | 0 | 0 | |
| commonCondLogLikDerDelta | 0.002 | 0.000 | 0.003 | |
| condLogLikDerSize | 0.000 | 0.000 | 0.001 | |
| cpm | 0.005 | 0.000 | 0.005 | |
| cutWithMinN | 0.002 | 0.000 | 0.002 | |
| decidetestsDGE | 0.022 | 0.000 | 0.021 | |
| dglmStdResid | 0.013 | 0.000 | 0.013 | |
| diffSpliceDGE | 0.06 | 0.00 | 0.06 | |
| dim | 0.002 | 0.000 | 0.002 | |
| dispBinTrend | 0.439 | 0.004 | 0.443 | |
| dispCoxReid | 0.021 | 0.000 | 0.021 | |
| dispCoxReidInterpolateTagwise | 0.023 | 0.000 | 0.023 | |
| dispCoxReidSplineTrend | 0.628 | 0.000 | 0.629 | |
| dropEmptyLevels | 0.001 | 0.000 | 0.001 | |
| edgeRUsersGuide | 0.001 | 0.000 | 0.001 | |
| effectiveLibSizes | 0.006 | 0.000 | 0.006 | |
| equalizeLibSizes | 0.016 | 0.000 | 0.016 | |
| estimateCommonDisp | 0.025 | 0.000 | 0.025 | |
| estimateDisp | 0.22 | 0.00 | 0.22 | |
| estimateExonGenewisedisp | 0.018 | 0.000 | 0.017 | |
| estimateGLMCommonDisp | 0.058 | 0.000 | 0.058 | |
| estimateGLMRobustDisp | 0.575 | 0.000 | 0.577 | |
| estimateGLMTagwiseDisp | 0.132 | 0.000 | 0.132 | |
| estimateGLMTrendedDisp | 0.128 | 0.000 | 0.128 | |
| estimateTagwiseDisp | 0.035 | 0.000 | 0.034 | |
| estimateTrendedDisp | 0.298 | 0.004 | 0.302 | |
| exactTest | 0.013 | 0.000 | 0.012 | |
| expandAsMatrix | 0 | 0 | 0 | |
| filterByExpr | 0 | 0 | 0 | |
| getCounts | 0.009 | 0.000 | 0.009 | |
| getPriorN | 0.002 | 0.000 | 0.002 | |
| gini | 0.000 | 0.000 | 0.001 | |
| glmLRT | 0 | 0 | 0 | |
| glmQLFTest | 0 | 0 | 0 | |
| glmQLFit | 0.438 | 0.012 | 0.451 | |
| glmTreat | 0.025 | 0.000 | 0.026 | |
| glmfit | 0.014 | 0.000 | 0.013 | |
| goana | 0 | 0 | 0 | |
| gof | 0.009 | 0.000 | 0.008 | |
| goodTuring | 0.005 | 0.000 | 0.005 | |
| head | 0.004 | 0.000 | 0.005 | |
| loessByCol | 0.001 | 0.000 | 0.002 | |
| maPlot | 0.033 | 0.000 | 0.034 | |
| makeCompressedMatrix | 0.001 | 0.002 | 0.002 | |
| maximizeInterpolant | 0.000 | 0.001 | 0.001 | |
| maximizeQuadratic | 0.000 | 0.001 | 0.001 | |
| meanvar | 0.083 | 0.004 | 0.087 | |
| mglm | 0.004 | 0.000 | 0.004 | |
| modelMatrixMeth | 0.004 | 0.000 | 0.003 | |
| movingAverageByCol | 0 | 0 | 0 | |
| nbinomDeviance | 0 | 0 | 0 | |
| nbinomUnitDeviance | 0 | 0 | 0 | |
| nearestReftoX | 0.000 | 0.000 | 0.001 | |
| nearestTSS | 6.682 | 0.228 | 6.928 | |
| plotBCV | 0.414 | 0.000 | 0.414 | |
| plotExonUsage | 0.02 | 0.00 | 0.02 | |
| plotMDS.DGEList | 0.019 | 0.000 | 0.018 | |
| plotQLDisp | 0.465 | 0.020 | 0.486 | |
| plotSmear | 0.439 | 0.016 | 0.455 | |
| predFC | 0.053 | 0.000 | 0.053 | |
| q2qnbinom | 0.001 | 0.000 | 0.000 | |
| read10X | 0 | 0 | 0 | |
| readDGE | 0 | 0 | 0 | |
| roast.DGEGLM | 0.087 | 0.000 | 0.087 | |
| roast.DGEList | 0.150 | 0.000 | 0.151 | |
| romer.DGEGLM | 5.144 | 0.024 | 5.178 | |
| romer.DGEList | 5.150 | 0.024 | 5.184 | |
| rowsum | 0.002 | 0.003 | 0.005 | |
| scaleOffset | 0 | 0 | 0 | |
| spliceVariants | 0.014 | 0.004 | 0.017 | |
| splitIntoGroups | 0.002 | 0.000 | 0.002 | |
| subsetting | 0.017 | 0.000 | 0.018 | |
| sumTechReps | 0 | 0 | 0 | |
| systematicSubset | 0 | 0 | 0 | |
| thinCounts | 0.001 | 0.000 | 0.001 | |
| topTags | 0.021 | 0.000 | 0.021 | |
| validDGEList | 0.001 | 0.000 | 0.001 | |
| weightedCondLogLikDerDelta | 0.002 | 0.000 | 0.002 | |
| zscoreNBinom | 0 | 0 | 0 | |