[med-svn] [r-bioc-limma] 01/03: Imported Upstream version 3.28.10+dfsg
Gordon Ball
chronitis-guest at moszumanska.debian.org
Tue Jun 21 13:32:21 UTC 2016
This is an automated email from the git hooks/post-receive script.
chronitis-guest pushed a commit to branch master
in repository r-bioc-limma.
commit decb5cf4150bc23e7eab6f6e2091b558f36c09c6
Author: Gordon Ball <gordon at chronitis.net>
Date: Tue Jun 21 15:17:15 2016 +0200
Imported Upstream version 3.28.10+dfsg
---
DESCRIPTION | 17 +-
NAMESPACE | 4 +-
R/alias2Symbol.R | 100 +++++--
R/detectionPValues.R | 41 +++
R/diffSplice.R | 15 +-
R/dups.R | 2 +-
R/ebayes.R | 48 ++--
R/{squeezeVar.R => fitFDist.R} | 59 +---
R/fitFDistRobustly.R | 66 +++--
R/fitmixture.R | 48 ++++
R/geneset-camera.R | 11 +-
R/geneset-fry.R | 232 ++++++++-------
R/geneset-roast.R | 633 ++++++++++++++++++-----------------------
R/geneset-romer.R | 15 +-
R/goana.R | 50 +++-
R/kegga.R | 211 +++++++++-----
R/lmEffects.R | 139 +++++++++
R/lmfit.R | 39 ++-
R/loessFit.R | 2 +-
R/neqc.R | 4 +-
R/norm.R | 7 +-
R/plotrldf.R | 133 +++++++--
R/propTrueNull.R | 6 +-
R/read-maimages.R | 1 +
R/read.R | 14 +-
R/read.idat.R | 127 ++++++---
R/sepchannel.R | 2 +-
R/squeezeVar.R | 203 +++----------
R/tricubeMovingAverage.R | 15 +-
R/voom.R | 15 +-
build/vignette.rds | Bin 230 -> 229 bytes
inst/NEWS.Rd | 225 +++++++++++++++
inst/doc/changelog.txt | 234 ++++++++++++++-
inst/doc/intro.pdf | Bin 46191 -> 46191 bytes
man/03reading.Rd | 1 +
man/10GeneSetTests.Rd | 2 +-
man/alias2Symbol.Rd | 14 +-
man/camera.Rd | 40 ++-
man/detectionPValue.Rd | 58 ++++
man/dupcor.Rd | 40 ++-
man/ebayes.Rd | 14 +-
man/fitfdist.Rd | 28 +-
man/fitmixture.Rd | 66 +++++
man/fitted.MArrayLM.Rd | 5 +-
man/getEAWP.Rd | 25 +-
man/gls.series.Rd | 4 +-
man/goana.Rd | 53 +++-
man/ids2indices.Rd | 6 +-
man/lmFit.Rd | 2 +-
man/logcosh.Rd | 26 ++
man/marraylm.Rd | 1 +
man/nec.Rd | 1 +
man/plotRLDF.Rd | 70 +++--
man/propexpr.Rd | 44 ++-
man/read.idat.Rd | 77 +++--
man/read.maimages.Rd | 3 +-
man/removeext.Rd | 13 +-
man/roast.Rd | 63 ++--
man/romer.Rd | 5 +-
man/squeezeVar.Rd | 10 +-
man/targetsA2C.Rd | 6 +
man/venn.Rd | 22 +-
man/voom.Rd | 12 +-
tests/limma-Tests.R | 4 +
tests/limma-Tests.Rout.save | 81 ++++--
65 files changed, 2326 insertions(+), 1188 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 192fe60..cb99d7f 100755
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,16 +1,17 @@
Package: limma
-Version: 3.26.1
-Date: 2015-10-30
+Version: 3.28.10
+Date: 2016-06-15
Title: Linear Models for Microarray Data
Description: Data analysis, linear models and differential expression for microarray data.
-Author: Gordon Smyth [cre,aut], Matthew Ritchie [ctb], Jeremy Silver [ctb], James Wettenhall [ctb], Natalie Thorne [ctb], Davis McCarthy [ctb], Di Wu [ctb], Yifang Hu [ctb], Wei Shi [ctb], Belinda Phipson [ctb], Alicia Oshlack [ctb], Carolyn de Graaf [ctb], Mette Langaas [ctb], Egil Ferkingstad [ctb], Marcus Davy [ctb], Francois Pepin [ctb], Dongseok Choi [ctb], Aaron Lun [ctb]
+Author: Gordon Smyth [cre,aut], Yifang Hu [ctb], Matthew Ritchie [ctb], Jeremy Silver [ctb], James Wettenhall [ctb], Davis McCarthy [ctb], Di Wu [ctb], Wei Shi [ctb], Belinda Phipson [ctb], Aaron Lun [ctb], Natalie Thorne [ctb], Alicia Oshlack [ctb], Carolyn de Graaf [ctb], Yunshun Chen [ctb], Mette Langaas [ctb], Egil Ferkingstad [ctb], Marcus Davy [ctb], Francois Pepin [ctb], Dongseok Choi [ctb]
Maintainer: Gordon Smyth <smyth at wehi.edu.au>
License: GPL (>=2)
-Depends: R (>= 2.3.0), methods
-Imports: grDevices, graphics, stats, utils
+Depends: R (>= 2.3.0)
+Imports: grDevices, graphics, stats, utils, methods
Suggests: affy, AnnotationDbi, BiasedUrn, Biobase, ellipse, GO.db,
- illuminaio, KEGGREST, locfit, MASS, org.Hs.eg.db, splines,
- statmod (>= 1.2.2), vsn
+ illuminaio, locfit, MASS, org.Hs.eg.db, org.Mm.eg.db,
+ org.Dm.eg.db, org.Pt.eg.db, org.Rn.eg.db, splines, statmod (>=
+ 1.2.2), vsn
URL: http://bioinf.wehi.edu.au/limma
biocViews: ExonArray, GeneExpression, Transcription,
AlternativeSplicing, DifferentialExpression,
@@ -21,4 +22,4 @@ biocViews: ExonArray, GeneExpression, Transcription,
MultipleComparison, Normalization, Preprocessing,
QualityControl
NeedsCompilation: yes
-Packaged: 2015-10-31 00:23:54 UTC; biocbuild
+Packaged: 2016-06-16 01:27:03 UTC; biocbuild
diff --git a/NAMESPACE b/NAMESPACE
index 492b181..ea342dc 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -13,7 +13,7 @@ importFrom("graphics", "abline", "axis", "barplot", "coplot", "image",
importFrom("stats", "approx", "approxfun", "as.dist", "cmdscale",
"coef", "contr.sum", "contrasts<-", "cov", "cov2cor",
"density", "df", "dhyper", "dnorm", "filter", "fitted",
- "hat", "integrate", "lm", "lm.fit", "lm.wfit", "loess",
+ "hat", "integrate", "isoreg", "lm", "lm.fit", "lm.wfit", "loess",
"lowess", "median", "model.matrix", "na.exclude", "na.omit",
"nlminb", "nls", "nls.control", "optim", "p.adjust",
"pbeta", "pchisq", "pexp", "pf", "pgamma", "phyper",
@@ -59,6 +59,8 @@ S3method(cbind,MAList)
S3method(cbind,RGList)
S3method(cbind,EList)
S3method(cbind,EListRaw)
+S3method(detectionPValues,default)
+S3method(detectionPValues,EListRaw)
S3method(dim,MAList)
S3method(dim,MArrayLM)
S3method(dim,RGList)
diff --git a/R/alias2Symbol.R b/R/alias2Symbol.R
index 855e9fd..46ef938 100644
--- a/R/alias2Symbol.R
+++ b/R/alias2Symbol.R
@@ -4,27 +4,49 @@ alias2Symbol <- function(alias,species="Hs",expand.symbols=FALSE)
# Convert a set of alias names to official gene symbols
# via Entrez Gene identifiers
# Di Wu, Gordon Smyth and Yifang Hu
-# 4 Sep 2008. Last revised 16 Jan 2015.
+# 4 Sep 2008. Last revised 2 May 2016.
{
alias <- as.character(alias)
species <- match.arg(species,c("Dm","Hs","Mm","Rn"))
- DB <- paste("org",species,"eg","db",sep=".")
- ALIAS2EG <- paste("org",species,"egALIAS2EG",sep=".")
- SYMBOL <- paste("org",species,"egSYMBOL",sep=".")
- suppressPackageStartupMessages(require(DB,character.only=TRUE))
- if(expand.symbols)
- {
- alias <- intersect(alias,AnnotationDbi::Rkeys(get(ALIAS2EG)))
- eg <- AnnotationDbi::mappedLkeys(get(ALIAS2EG)[alias])
- AnnotationDbi::mappedRkeys(get(SYMBOL)[eg])
- }
- else
- {
- isSymbol <- alias %in% AnnotationDbi::Rkeys(get(SYMBOL))
- alias2 <- intersect(alias[!isSymbol],AnnotationDbi::Rkeys(get(ALIAS2EG)))
- eg <- AnnotationDbi::mappedLkeys(get(ALIAS2EG)[alias2])
- c(alias[isSymbol],AnnotationDbi::mappedRkeys(get(SYMBOL)[eg]))
+# Get access to required annotation functions
+ suppressPackageStartupMessages(OK <- requireNamespace("AnnotationDbi",quietly=TRUE))
+ if(!OK) stop("AnnotationDbi package required but not installed (or can't be loaded)")
+
+# Get alias to symbol mappings
+ switch(species,
+ Hs = {
+ suppressPackageStartupMessages(OK <- requireNamespace("org.Hs.eg.db",quietly=TRUE))
+ if(!OK) stop("org.Hs.eg.db package required but not installed (or can't be loaded)")
+ ALIAS2EG <- org.Hs.eg.db::org.Hs.egALIAS2EG
+ SYMBOL <- org.Hs.eg.db::org.Hs.egSYMBOL
+ }, Mm = {
+ suppressPackageStartupMessages(OK <- requireNamespace("org.Mm.eg.db",quietly=TRUE))
+ if(!OK) stop("org.Mm.eg.db package required but not installed (or can't be loaded)")
+ ALIAS2EG <- org.Mm.eg.db::org.Mm.egALIAS2EG
+ SYMBOL <- org.Mm.eg.db::org.Mm.egSYMBOL
+ }, Rn = {
+ suppressPackageStartupMessages(OK <- requireNamespace("org.Rn.eg.db",quietly=TRUE))
+ if(!OK) stop("org.Rn.eg.db package required but not installed (or can't be loaded)")
+ ALIAS2EG <- org.Rn.eg.db::org.Rn.egALIAS2EG
+ SYMBOL <- org.Rn.eg.db::org.Rn.egSYMBOL
+ }, Dm = {
+ suppressPackageStartupMessages(OK <- requireNamespace("org.Dm.eg.db",quietly=TRUE))
+ if(!OK) stop("org.Dm.eg.db package required but not installed (or can't be loaded)")
+ ALIAS2EG <- org.Dm.eg.db::org.Dm.egALIAS2EG
+ SYMBOL <- org.Dm.eg.db::org.Dm.egSYMBOL
+ }
+ )
+
+ if(expand.symbols) {
+ alias <- intersect(alias,AnnotationDbi::Rkeys(ALIAS2EG))
+ eg <- AnnotationDbi::mappedLkeys(ALIAS2EG[alias])
+ AnnotationDbi::mappedRkeys(SYMBOL[eg])
+ } else {
+ isSymbol <- alias %in% AnnotationDbi::Rkeys(SYMBOL)
+ alias2 <- intersect(alias[!isSymbol],AnnotationDbi::Rkeys(ALIAS2EG))
+ eg <- AnnotationDbi::mappedLkeys(ALIAS2EG[alias2])
+ c(alias[isSymbol],AnnotationDbi::mappedRkeys(SYMBOL[eg]))
}
}
@@ -32,31 +54,55 @@ alias2SymbolTable <- function(alias,species="Hs")
# Convert a vector of alias names to the vector of corresponding official gene symbols
# via Entrez Gene identifiers
# Di Wu, Gordon Smyth and Yifang Hu
-# Created 3 Sep 2009. Last modified 16 Jan 2015.
+# Created 3 Sep 2009. Last modified 13 April 2016.
{
alias <- as.character(alias)
species <- match.arg(species,c("Dm","Hs","Mm","Rn"))
- DB <- paste("org",species,"eg","db",sep=".")
- ALIAS2EG <- paste("org",species,"egALIAS2EG",sep=".")
- SYMBOL <- paste("org",species,"egSYMBOL",sep=".")
- suppressPackageStartupMessages(require(DB,character.only=TRUE))
- isSymbol <- alias %in% AnnotationDbi::Rkeys(get(SYMBOL))
+# Get access to required annotation functions
+ suppressPackageStartupMessages(OK <- requireNamespace("AnnotationDbi",quietly=TRUE))
+ if(!OK) stop("AnnotationDbi package required but not installed (or can't be loaded)")
+
+# Get alias to symbol mappings
+ switch(species,
+ Hs = {
+ suppressPackageStartupMessages(OK <- requireNamespace("org.Hs.eg.db",quietly=TRUE))
+ if(!OK) stop("org.Hs.eg.db package required but not installed (or can't be loaded)")
+ ALIAS2EG <- org.Hs.eg.db::org.Hs.egALIAS2EG
+ SYMBOL <- org.Hs.eg.db::org.Hs.egSYMBOL
+ }, Mm = {
+ suppressPackageStartupMessages(OK <- requireNamespace("org.Mm.eg.db",quietly=TRUE))
+ if(!OK) stop("org.Mm.eg.db package required but not installed (or can't be loaded)")
+ ALIAS2EG <- org.Mm.eg.db::org.Mm.egALIAS2EG
+ SYMBOL <- org.Mm.eg.db::org.Mm.egSYMBOL
+ }, Rn = {
+ suppressPackageStartupMessages(OK <- requireNamespace("org.Rn.eg.db",quietly=TRUE))
+ if(!OK) stop("org.Rn.eg.db package required but not installed (or can't be loaded)")
+ ALIAS2EG <- org.Rn.eg.db::org.Rn.egALIAS2EG
+ SYMBOL <- org.Rn.eg.db::org.Rn.egSYMBOL
+ }, Dm = {
+ suppressPackageStartupMessages(OK <- requireNamespace("org.Dm.eg.db",quietly=TRUE))
+ if(!OK) stop("org.Dm.eg.db package required but not installed (or can't be loaded)")
+ ALIAS2EG <- org.Dm.eg.db::org.Dm.egALIAS2EG
+ SYMBOL <- org.Dm.eg.db::org.Dm.egSYMBOL
+ }
+ )
+
+ isSymbol <- alias %in% AnnotationDbi::Rkeys(SYMBOL)
Symbol <- alias
Symbol[!isSymbol] <- NA
OtherAliases <- alias[!isSymbol]
- isAlias <- OtherAliases %in% AnnotationDbi::Rkeys(get(ALIAS2EG))
+ isAlias <- OtherAliases %in% AnnotationDbi::Rkeys(ALIAS2EG)
if(!any(isAlias)) return(Symbol)
OtherAliases <- OtherAliases[isAlias]
- AliasTbl <- AnnotationDbi::toTable(get(ALIAS2EG)[OtherAliases])
+ AliasTbl <- AnnotationDbi::toTable(ALIAS2EG[OtherAliases])
if(anyDuplicated(AliasTbl$alias_symbol)) warning("Multiple symbols ignored for one or more aliases")
- SymbolTbl <- AnnotationDbi::toTable(get(SYMBOL)[AliasTbl$gene_id])
+ SymbolTbl <- AnnotationDbi::toTable(SYMBOL[AliasTbl$gene_id])
m <- match(OtherAliases,AliasTbl$alias_symbol)
GeneID <- AliasTbl$gene_id[m]
m <- match(GeneID,SymbolTbl$gene_id)
Symbol[!isSymbol][isAlias] <- SymbolTbl$symbol[m]
Symbol
}
-
diff --git a/R/detectionPValues.R b/R/detectionPValues.R
new file mode 100644
index 0000000..a0f8a8f
--- /dev/null
+++ b/R/detectionPValues.R
@@ -0,0 +1,41 @@
+detectionPValues <- function(x,...) UseMethod("detectionPValues")
+
+detectionPValues.EListRaw <- function(x,status=NULL,...)
+# Detection p-values from negative controls for EListRaw
+# Gordon Smyth
+# Created 13 June 2016
+{
+ if(is.null(status)) status <- x$genes$Status
+ detectionPValues(x=x$E, status=status, ...)
+}
+
+detectionPValues.default <- function(x, status, negctrl="negative",...)
+# Detection p-values from negative controls
+# Gordon Smyth
+# Created 13 June 2016. Last modified 14 June 2016.
+{
+ x <- as.matrix(x)
+
+# Identify negative control rows
+ if(missing(status) || is.null(status)) stop("need to specify status vector")
+ isneg <- as.integer(status==negctrl)
+ notneg <- 1L-isneg
+ nneg <- sum(isneg)
+
+# Count the number of negative controls greater than each value
+# cs1 counts number of negative controls >= the observed value
+# cs2 counts number of negative controls > the observed value
+# By taking the average of cs1 and cs2, we count ties as 1/2
+ DetectionPValue <- x
+ cs1 <- cs2 <- rep_len(0,length.out=nrow(x))
+ for (j in 1:ncol(x)) {
+ o1 <- order(x[,j],isneg,decreasing=TRUE)
+ o2 <- order(x[,j],notneg,decreasing=TRUE)
+ cs1[o1] <- cumsum(isneg[o1])
+ cs2[o2] <- cumsum(isneg[o2])
+ DetectionPValue[,j] <- cs1+cs2
+ }
+
+# Convert to p-values
+ DetectionPValue/(2*nneg)
+}
diff --git a/R/diffSplice.R b/R/diffSplice.R
index cc1b1ca..1dd77c5 100644
--- a/R/diffSplice.R
+++ b/R/diffSplice.R
@@ -2,10 +2,13 @@ diffSplice <- function(fit,geneid,exonid=NULL,robust=FALSE,verbose=TRUE)
# Test for splicing variants between conditions
# using linear model fit of exon data.
# Gordon Smyth and Charity Law
-# Created 13 Dec 2013. Last modified 22 Sep 2014.
+# Created 13 Dec 2013. Last modified 18 Jan 2016.
{
+# Make sure there is always an annotation frame
exon.genes <- fit$genes
if(is.null(exon.genes)) exon.genes <- data.frame(ExonID=1:nrow(fit))
+
+# Get ID columns for genes and exons
if(length(geneid)==1) {
genecolname <- as.character(geneid)
geneid <- exon.genes[[genecolname]]
@@ -13,7 +16,9 @@ diffSplice <- function(fit,geneid,exonid=NULL,robust=FALSE,verbose=TRUE)
exon.genes$GeneID <- geneid
genecolname <- "GeneID"
}
- if(!is.null(exonid))
+ if(is.null(exonid)) {
+ exoncolname <- NULL
+ } else {
if(length(exonid)==1) {
exoncolname <- as.character(exonid)
exonid <- exon.genes[[exoncolname]]
@@ -21,8 +26,7 @@ diffSplice <- function(fit,geneid,exonid=NULL,robust=FALSE,verbose=TRUE)
exon.genes$ExonID <- exonid
exoncolname <- "ExonID"
}
- else
- exoncolname <- NULL
+ }
# Sort by geneid
if(is.null(exonid))
@@ -36,7 +40,7 @@ diffSplice <- function(fit,geneid,exonid=NULL,robust=FALSE,verbose=TRUE)
exon.df.residual <- fit$df.residual[o]
exon.s2 <- fit$sigma[o]^2
-# Count exons and get genewise variances
+# Count exons by gene and get genewise variances
exon.stat <- cbind(1,exon.df.residual,exon.s2)
gene.sum <- rowsum(exon.stat,geneid,reorder=FALSE)
gene.nexons <- gene.sum[,1]
@@ -102,6 +106,7 @@ diffSplice <- function(fit,geneid,exonid=NULL,robust=FALSE,verbose=TRUE)
out$gene.df.prior <- squeeze$df.prior
out$gene.df.residual <- gene.df.residual
out$gene.df.total <- gene.df.total
+ out$gene.s2 <- gene.s2[gene.keep]
out$gene.s2.post <- gene.s2.post
out$gene.F <- gene.F
out$gene.F.p.value <- gene.F.p.value
diff --git a/R/dups.R b/R/dups.R
index 90430bf..d6e8802 100755
--- a/R/dups.R
+++ b/R/dups.R
@@ -97,7 +97,7 @@ duplicateCorrelation <- function(object,design=NULL,ndups=2,spacing=1,block=NULL
design <- design %x% rep(1,ndups)
}
- if(!requireNamespace("statmod",quietly=TRUE)) stop("statmod package required but is not available")
+ if(!requireNamespace("statmod",quietly=TRUE)) stop("statmod package required but is not installed")
rho <- rep(NA,ngenes)
nafun <- function(e) NA
for (i in 1:ngenes) {
diff --git a/R/ebayes.R b/R/ebayes.R
index 314fd05..5a33cc3 100755
--- a/R/ebayes.R
+++ b/R/ebayes.R
@@ -86,11 +86,11 @@ ebayes <- function(fit,proportion=0.01,stdev.coef.lim=c(0.1,4),trend=FALSE,robus
out
}
-tmixture.matrix <- function(tstat,stdev.unscaled,df,proportion,v0.lim=NULL) {
+tmixture.matrix <- function(tstat,stdev.unscaled,df,proportion,v0.lim=NULL)
# Estimate the prior variance of the coefficients for DE genes
# Gordon Smyth
# 18 Nov 2002. Last modified 12 Dec 2003.
-
+{
tstat <- as.matrix(tstat)
stdev.unscaled <- as.matrix(stdev.unscaled)
if(any(dim(tstat) != dim(stdev.unscaled))) stop("Dims of tstat and stdev.unscaled don't match")
@@ -101,19 +101,22 @@ tmixture.matrix <- function(tstat,stdev.unscaled,df,proportion,v0.lim=NULL) {
v0
}
-tmixture.vector <- function(tstat,stdev.unscaled,df,proportion,v0.lim=NULL) {
+tmixture.vector <- function(tstat,stdev.unscaled,df,proportion,v0.lim=NULL)
# Estimate scale factor in mixture of two t-distributions
# tstat is assumed to follow sqrt(1+v0/v1)*t(df) with probability proportion and t(df) otherwise
# v1 is stdev.unscaled^2 and v0 is to be estimated
# Gordon Smyth
-# 18 Nov 2002. Last modified 13 Dec 2003.
-
- if(any(is.na(tstat))) {
+# 18 Nov 2002. Last modified 15 April 2016.
+{
+# Remove missing values
+ if(anyNA(tstat)) {
o <- !is.na(tstat)
tstat <- tstat[o]
stdev.unscaled <- stdev.unscaled[o]
df <- df[o]
}
+
+# ntarget t-statistics will be used for estimation
ngenes <- length(tstat)
ntarget <- ceiling(proportion/2*ngenes)
if(ntarget < 1) return(NA)
@@ -122,22 +125,29 @@ tmixture.vector <- function(tstat,stdev.unscaled,df,proportion,v0.lim=NULL) {
# This ensures ptarget < 1
p <- max(ntarget/ngenes,proportion)
+# Method requires that df be equal
tstat <- abs(tstat)
- if(ngenes>1)
- ttarget <- quantile(tstat,(ngenes-ntarget)/(ngenes-1))
- else
- ttarget <- tstat
- top <- (tstat >= ttarget)
- tstat <- tstat[top]
- v1 <- stdev.unscaled[top]^2
- df <- df[top]
- r <- ntarget-rank(tstat)+1
- p0 <- pt(-tstat,df=df)
- ptarget <- ( (r-0.5)/2/ngenes - (1-p)*p0 ) / p
+ MaxDF <- max(df)
+ i <- df < MaxDF
+ if(any(i)) {
+ TailP <- pt(tstat[i],df=df[i],lower.tail=FALSE,log.p=TRUE)
+ tstat[i] <- qt(TailP,df=MaxDF,lower.tail=FALSE,log.p=TRUE)
+ df[i] <- MaxDF
+ }
+
+# Select top statistics
+ o <- order(tstat,decreasing=TRUE)[1:ntarget]
+ tstat <- tstat[o]
+ v1 <- stdev.unscaled[o]^2
+
+# Compare to order statistics
+ r <- 1:ntarget
+ p0 <- 2*pt(tstat,df=MaxDF,lower.tail=FALSE)
+ ptarget <- ( (r-0.5)/ngenes - (1-p)*p0 ) / p
+ v0 <- rep.int(0,ntarget)
pos <- ptarget > p0
- v0 <- rep(0,ntarget)
if(any(pos)) {
- qtarget <- qt(ptarget[pos],df=df[pos])
+ qtarget <- qt(ptarget[pos]/2,df=MaxDF,lower.tail=FALSE)
v0[pos] <- v1[pos]*((tstat[pos]/qtarget)^2-1)
}
if(!is.null(v0.lim)) v0 <- pmin(pmax(v0,v0.lim[1]),v0.lim[2])
diff --git a/R/squeezeVar.R b/R/fitFDist.R
similarity index 66%
copy from R/squeezeVar.R
copy to R/fitFDist.R
index 387d962..db9bfa6 100644
--- a/R/squeezeVar.R
+++ b/R/fitFDist.R
@@ -1,59 +1,6 @@
-# EMPIRICAL BAYES SQUEEZING OF VARIANCES
-
-squeezeVar <- function(var, df, covariate=NULL, robust=FALSE, winsor.tail.p=c(0.05,0.1))
-# Empirical Bayes posterior variances
-# Gordon Smyth
-# 2 March 2004. Last modified 2 Dec 2013.
-{
- n <- length(var)
- if(n == 0) stop("var is empty")
- if(n == 1) return(list(var.post=var,var.prior=var,df.prior=0))
- if(length(df)==1) {
- df <- rep.int(df,n)
- } else {
- if(length(df) != n) stop("lengths differ")
- }
-
-# Estimate prior var and df
- if(robust)
- fit <- fitFDistRobustly(var, df1=df, covariate=covariate, winsor.tail.p=winsor.tail.p)
- else
- fit <- fitFDist(var, df1=df, covariate=covariate)
-
-# Prior var will be vector if robust=TRUE, otherwise scalar
- var.prior <- fit$scale
-
-# Prior df will be vector if covariate is non-NULL, otherwise scalar
- df.prior <- fit$df2.shrunk
- if(is.null(df.prior)) df.prior <- fit$df2
-
-# Check estimated prior df
- if(is.null(df.prior) || any(is.na(df.prior))) stop("Could not estimate prior df")
-
-# Squeeze the posterior variances
- df.total <- df + df.prior
- var[df==0] <- 0 # guard against missing or infinite values
- Infdf <- df.prior==Inf
- if(any(Infdf)) {
- var.post <- rep(var.prior,length.out=n)
- i <- which(!Infdf)
- if(length(i)) {
- if(is.null(covariate))
- s02 <- var.prior
- else
- s02 <- var.prior[i]
- var.post[i] <- (df[i]*var[i] + df.prior[i]*s02) / df.total[i]
- }
- } else {
- var.post <- (df*var + df.prior*var.prior) / df.total
- }
-
- list(df.prior=df.prior,var.prior=var.prior,var.post=var.post)
-}
-
fitFDist <- function(x,df1,covariate=NULL)
-# Moment estimation of the parameters of a scaled F-distribution
-# The first degrees of freedom are given
+# Moment estimation of the parameters of a scaled F-distribution.
+# The numerator degrees of freedom are given, the denominator is to be estimated.
# Gordon Smyth and Belinda Phipson
# 8 Sept 2002. Last revised 27 Oct 2012.
{
@@ -132,6 +79,7 @@ fitFDist <- function(x,df1,covariate=NULL)
list(scale=s20,df2=df2)
}
+
trigammaInverse <- function(x) {
# Solve trigamma(y) = x for y
# Gordon Smyth
@@ -189,4 +137,3 @@ trigammaInverse <- function(x) {
}
y
}
-
diff --git a/R/fitFDistRobustly.R b/R/fitFDistRobustly.R
index b2d13cf..0d82ad0 100644
--- a/R/fitFDistRobustly.R
+++ b/R/fitFDistRobustly.R
@@ -3,20 +3,17 @@ fitFDistRobustly <- function(x,df1,covariate=NULL,winsor.tail.p=c(0.05,0.1),trac
# given the first degrees of freedom, using first and second
# moments of Winsorized z-values
# Gordon Smyth and Belinda Phipson
-# 8 Sept 2002. Last revised 14 January 2015.
+# Created 7 Oct 2012. Last revised 20 November 2015.
{
# Check x
n <- length(x)
-# Eliminate case of no useful data
- if(n==0)
- if(is.null(covariate))
- return(list(scale=NaN,df2=NaN,df2.robust=NaN))
- else
- return(list(scale=numeric(0),df2=numeric(0),df2.robust=NaN))
+# Eliminate cases of no useful data
+ if(n<2) return(list(scale=NA,df2=NA,df2.shrunk=NA))
+ if(n==2) return(fitFDist(x=x,df1=df1,covariate=covariate))
# Check df1
- if(length(df1)>1) if(length(df1) != n) stop("x and df1 are different lengths")
+ if(all(length(df1) != c(1,n))) stop("x and df1 are different lengths")
# Check covariate
if(!is.null(covariate)) {
@@ -110,7 +107,7 @@ fitFDistRobustly <- function(x,df1,covariate=NULL,winsor.tail.p=c(0.05,0.1),trac
if(trace) cat("Variance of Winsorized Fisher-z",zwvar,"\n")
# Theoretical Winsorized moments
- if(!requireNamespace("statmod",quietly=TRUE)) stop("statmod package required but is not available")
+ if(!requireNamespace("statmod",quietly=TRUE)) stop("statmod package required but is not installed")
g <- statmod::gauss.quad.prob(128,dist="uniform")
linkfun <- function(x) x/(1+x)
linkinv <- function(x) x/(1-x)
@@ -189,28 +186,59 @@ fitFDistRobustly <- function(x,df1,covariate=NULL,winsor.tail.p=c(0.05,0.1),trac
r <- rank(Fstat)
EmpiricalTailProb <- (n-r+0.5)/n
ProbNotOutlier <- pmin(TailP/EmpiricalTailProb,1)
- df2.shrunk <- rep.int(df2,n)
+ ProbOutlier <- 1-ProbNotOutlier
if(any(ProbNotOutlier<1)) {
- ProbOutlier <- 1-ProbNotOutlier
+ o <- order(TailP)
+
+# Old calculation for df2.outlier
VarOutlier <- max(zresid)^2
VarOutlier <- VarOutlier-trigamma(df1/2)
if(trace) cat("VarOutlier",VarOutlier,"\n")
if(VarOutlier > 0) {
- df2.outlier <- 2*trigammaInverse(VarOutlier)
- if(trace) cat("df2.outlier",df2.outlier,"\n")
- if(df2.outlier < df2) {
- df2.shrunk <- ProbNotOutlier*df2+ProbOutlier*df2.outlier
- o <- order(TailP)
- df2.ordered <- df2.shrunk[o]
+ df2.outlier.old <- 2*trigammaInverse(VarOutlier)
+ if(trace) cat("df2.outlier.old",df2.outlier.old,"\n")
+ if(df2.outlier.old < df2) {
+ df2.shrunk.old <- ProbNotOutlier*df2+ProbOutlier*df2.outlier.old
+# Make df2.shrunk.old monotonic in TailP
+ df2.ordered <- df2.shrunk.old[o]
df2.ordered[1] <- min(df2.ordered[1],NonRobust$df2)
m <- cumsum(df2.ordered)
m <- m/(1:n)
imin <- which.min(m)
df2.ordered[1:imin] <- m[imin]
- df2.shrunk[o] <- cummax(df2.ordered)
+ df2.shrunk.old[o] <- cummax(df2.ordered)
}
}
+
+# New calculation for df2.outlier
+# Find df2.outlier to make maxFstat the median of the distribution
+# Exploit fact that log(TailP) is nearly linearly with positive 2nd deriv as a function of df2
+# Note that minTailP and NewTailP are always less than 0.5
+ df2.outlier <- log(0.5)/log(min(TailP))*df2
+# Iterate for accuracy
+ NewTailP <- pf(max(Fstat),df1=df1,df2=df2.outlier,lower.tail=FALSE)
+ df2.outlier <- log(0.5)/log(NewTailP)*df2.outlier
+ df2.shrunk <- ProbNotOutlier*df2+ProbOutlier*df2.outlier
+
+# Force df2.shrunk to be monotonic in TailP
+ o <- order(TailP)
+ df2.ordered <- df2.shrunk[o]
+ m <- cumsum(df2.ordered)
+ m <- m/(1:n)
+ imin <- which.min(m)
+ df2.ordered[1:imin] <- m[imin]
+ df2.shrunk[o] <- cummax(df2.ordered)
+
+# Use isoreg() instead. This gives similar results.
+ df2.shrunk.iso <- rep.int(df2,n)
+ o <- o[1:(n/2)]
+ df2.shrunk.iso[o] <- ProbNotOutlier[o]*df2+ProbOutlier[o]*df2.outlier
+ df2.shrunk.iso[o] <- isoreg(TailP[o],df2.shrunk.iso[o])$yf
+
+ } else {
+ df2.outlier <- df2.outlier2 <- df2
+ df2.shrunk2 <- df2.shrunk <- rep.int(df2,n)
}
- list(scale=s20,df2=df2,df2.shrunk=df2.shrunk)
+ list(scale=s20,df2=df2,tail.p.value=TailP,prob.outlier=ProbOutlier,df2.outlier=df2.outlier,df2.outlier.old=df2.outlier.old,df2.shrunk=df2.shrunk,df2.shrunk.old=df2.shrunk.old,df2.shrunk.iso=df2.shrunk.iso)
}
diff --git a/R/fitmixture.R b/R/fitmixture.R
new file mode 100644
index 0000000..0d3ba10
--- /dev/null
+++ b/R/fitmixture.R
@@ -0,0 +1,48 @@
+logcosh <- function(x)
+# Compute log(cosh(x)) without floating over or underflow
+# Gordon Smyth 8 April 2007
+{
+ y <- abs(x)-log(2)
+ i <- abs(x) < 1e-4
+ y[i] <- 0.5*x[i]^2
+ i <- !i & (abs(x) < 17)
+ y[i] <- log(cosh(x[i]))
+ y
+}
+
+fitmixture <- function(log2e,mixprop,niter=4,trace=FALSE)
+# Fit mixture model by non-linear least squares
+# Gordon Smyth 9 April 2007
+{
+ log2e <- as.matrix(log2e)
+ mixprop <- as.numeric(mixprop)
+ narrays <- ncol(log2e)
+ nprobes <- nrow(log2e)
+ y <- 2^log2e
+ start <- lm.fit(cbind(mixprop,1-mixprop),t(y))$coef
+ start <- pmax(start,1)
+ b <- as.vector(c(1,-1) %*% log(start))
+ z <- log(y)
+ pm <- matrix(1,nprobes,1) %*% matrix(2*mixprop-1,1,narrays)
+ for (i in 1:niter) {
+ mub <- logcosh(b/2)+log(1+tanh(b/2)*pm)
+ a <- rowMeans(z-mub)
+ mu <- a+mub
+ if(trace) {
+ s <- sqrt(narrays/(narrays-2)*rowMeans((z-mu)^2))
+ if(i>1) cat("stdev changes",summary(sold-s),"\n")
+ sold <- s
+ }
+ tb <- tanh(b/2)
+ dmu <- (tb+pm)/(1+tb*pm)/2
+ b <- b + rowMeans(dmu*(z-mu)) / (1e-6+rowMeans((dmu-rowMeans(dmu))^2))
+ if(trace) cat("betas",summary(b),"\n")
+ }
+ mub <- logcosh(b/2)+log(1+tanh(b/2)*pm)
+ a <- rowMeans(z-mub)
+ mu <- a+mub
+ s <- sqrt(narrays/(narrays-2)*rowMeans((z-mu)^2))
+ l2 <- log(2)
+ list(A=a/l2,M=b/l2,stdev=s/l2)
+}
+
diff --git a/R/geneset-camera.R b/R/geneset-camera.R
index 5c9426a..8c9a3a5 100644
--- a/R/geneset-camera.R
+++ b/R/geneset-camera.R
@@ -21,10 +21,10 @@ interGeneCorrelation <- function(y, design)
camera <- function(y,...) UseMethod("camera")
-camera.default <- function(y,index,design=NULL,contrast=ncol(design),weights=NULL,use.ranks=FALSE,allow.neg.cor=TRUE,inter.gene.cor=NULL,trend.var=FALSE,sort=TRUE,...)
+camera.default <- function(y,index,design=NULL,contrast=ncol(design),weights=NULL,use.ranks=FALSE,allow.neg.cor=FALSE,inter.gene.cor=0.01,trend.var=FALSE,sort=TRUE,...)
# Competitive gene set test allowing for correlation between genes
# Gordon Smyth and Di Wu
-# Created 2007. Last modified 30 October 2015
+# Created 2007. Last modified 4 June 2016.
{
# Issue warning if extra arguments found
dots <- names(list(...))
@@ -45,7 +45,7 @@ camera.default <- function(y,index,design=NULL,contrast=ncol(design),weights=NUL
# Check design
if(is.null(design)) design <- y$design
if(is.null(design))
- design <- matrix(1,n,1)
+ stop("design matrix not specified")
else {
design <- as.matrix(design)
if(mode(design) != "numeric") stop("design must be a numeric matrix")
@@ -59,7 +59,7 @@ camera.default <- function(y,index,design=NULL,contrast=ncol(design),weights=NUL
if(is.null(weights)) weights <- y$weights
# Check inter.gene.cor
- fixed.cor <- !is.null(inter.gene.cor)
+ fixed.cor <- !(is.na(inter.gene.cor) || is.null(inter.gene.cor))
# Set df for camera tests
if(fixed.cor) {
@@ -197,6 +197,9 @@ camera.default <- function(y,index,design=NULL,contrast=ncol(design),weights=NUL
tab$PValue <- tab$TwoSided
tab$Down <- tab$Up <- tab$TwoSided <- NULL
+# Remove correlation column if it was not estimated
+ if(fixed.cor) tab$Correlation <- NULL
+
# Add FDR
if(nsets>1) tab$FDR <- p.adjust(tab$PValue,method="BH")
diff --git a/R/geneset-fry.R b/R/geneset-fry.R
index cc969fd..18f5989 100644
--- a/R/geneset-fry.R
+++ b/R/geneset-fry.R
@@ -1,128 +1,146 @@
fry <- function(y,...) UseMethod("fry")
-fry.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),weights=NULL,sort=TRUE,...)
+fry.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),geneid=NULL,standardize="posterior.sd",sort="directional",...)
# Quick version of roast gene set test assuming equal variances between genes
# The up and down p-values are equivalent to those from roast with nrot=Inf
# in the special case of prior.df=Inf.
# Gordon Smyth and Goknur Giner
-# Created 30 January 2015. Last modified 6 October 2015
+# Created 30 January 2015. Last modified 11 May 2016
{
-# Issue warning if extra arguments found
- dots <- names(list(...))
- if(length(dots)) warning("Extra arguments disregarded: ",sQuote(dots))
-
-# Extract components from y
- y <- getEAWP(y)
- G <- nrow(y$exprs)
- n <- ncol(y$exprs)
+# Partial matching of extra arguments
+ Dots <- list(...)
+ PossibleArgs <- c("array.weights","weights","block","correlation","trend","robust","winsor.tail.p")
+ if(!is.null(names(Dots))) {
+ i <- pmatch(names(Dots),PossibleArgs)
+ names(Dots) <- PossibleArgs[i]
+ }
-# Check index
- if(is.null(index)) index <- list(set1=1L:G)
- if(!is.list(index)) index <- list(set1=index)
- nsets <- length(index)
- if(nsets==0) stop("index is empty")
+# Defaults for extra arguments
+ if(is.null(Dots$trend)) Dots$trend <- FALSE
+ if(is.null(Dots$robust)) Dots$robust <- FALSE
+ if(Dots$robust & is.null(Dots$winsor.tail.p)) Dots$winsor.tail.p <- c(0.05,0.1)
+ if(!is.null(Dots$block) & is.null(Dots$correlation)) stop("correlation must be set")
+
+# Covariate for trended eBayes
+ covariate <- NULL
+ if(Dots$trend) covariate <- rowMeans(as.matrix(y))
+
+# Compute effects matrix (with df.residual+1 columns)
+ Effects <- .lmEffects(y=y,design=design,contrast=contrast,
+ array.weights=Dots$array.weights,
+ weights=Dots$weights,
+ block=Dots$block,
+ correlation=Dots$correlation)
+
+# Divide out genewise standard deviations
+ standardize <- match.arg(standardize, c("none","residual.sd","posterior.sd","p2"))
+ if(!standardize=="none") {
+
+# Estimate genewise sds robustly
+ OK <- requireNamespace("statmod",quietly=TRUE)
+ if(!OK) stop("statmod package required but not installed")
+ gq <- statmod::gauss.quad.prob(128,"uniform")
+ df.residual <- ncol(Effects)-1
+ Eu2max <- sum( (df.residual+1)*gq$nodes^df.residual*qchisq(gq$nodes,df=1)*gq$weights )
+ u2max <- apply(Effects^2,1,max)
+ s2.robust <- (rowSums(Effects^2)-u2max) / (df.residual+1-Eu2max)
+
+# Empirical Bayes moderation method I: estimate hyperparameters from robust variances
+ if(standardize=="p2") {
+ sv <- squeezeVar(s2.robust, df=0.92*df.residual,
+ covariate=covariate,
+ robust=Dots$robust,
+ winsor.tail.p=Dots$winsor.tail.p)
+ s2.robust <- sv$var.post
+ }
-# Check design
- if(is.null(design)) design <- y$design
- if(is.null(design))
- design <- matrix(1,n,1)
- else {
- design <- as.matrix(design)
- if(mode(design) != "numeric") stop("design must be a numeric matrix")
- }
- if(nrow(design) != n) stop("row dimension of design matrix must match column dimension of data")
- ne <- nonEstimable(design)
- if(!is.null(ne)) cat("Coefficients not estimable:",paste(ne,collapse=" "),"\n")
-
- p <- ncol(design)
- df.residual <- n-p
-
-# Check weights
- if(is.null(weights)) weights <- y$weights
-
-# Reduce to numeric expression matrix
- y <- y$exprs
-
-# Check weights
- if(!is.null(weights)) {
- if(any(weights<=0)) stop("weights must be positive")
- if(length(weights)==n) {
- sw <- sqrt(weights)
- y <- t(t(y)*sw)
- design <- design*sw
- weights <- NULL
+# Empirical Bayes moderation method II: estimate hyperparameters from residual variances but apply squeezing to robust variances
+ if(standardize=="posterior.sd") {
+ s2 <- rowMeans(Effects[,-1,drop=FALSE]^2)
+ if(Dots$robust) {
+ fit <- fitFDistRobustly(s2, df1=df.residual, covariate=covariate, winsor.tail.p=Dots$winsor.tail.p)
+ df.prior <- fit$df2.shrunk
+ } else {
+ fit <- fitFDist(s2, df1=df.residual, covariate=covariate)
+ df.prior <- fit$df2
+ }
+ s2.robust <- .squeezeVar(s2.robust,df=0.92*df.residual,var.prior=fit$scale,df.prior=df.prior)
}
- }
- if(!is.null(weights)) {
- if(length(weights)==G) weights <- matrix(weights,G,n)
- weights <- as.matrix(weights)
- if(any( dim(weights) != dim(y) )) stop("weights not conformal with y")
- }
-# Reform design matrix so that contrast of interest is last column
- if(is.character(contrast)) {
- contrast <- which(contrast==colnames(design))
- if(length(contrast)==0) stop("coef ",contrast," not found")
- }
- if(length(contrast)==1) {
- j <- c((1:p)[-contrast], contrast)
- if(contrast<p) design <- design[,j]
- } else {
- QR <- qr(contrast)
- design <- t(qr.qty(QR,t(design)))
- if(sign(QR$qr[1,1]<0)) design[,1] <- -design[,1]
- design <- design[,c(2:p,1)]
+ Effects <- Effects/sqrt(s2.robust)
}
-# Compute effects matrix
- if(is.null(weights)) {
- QR <- qr(design)
- if(QR$rank<p) stop("design matrix is not of full rank")
- effects <- qr.qty(QR,t(y))
- unscaledt <- effects[p,]
- if(QR$qr[p,p]<0) unscaledt <- -unscaledt
+# Check geneid (can be a vector of gene IDs or an annotation column name)
+ if(is.null(geneid)) {
+ geneid <- rownames(Effects)
} else {
- effects <- matrix(0,n,G)
- unscaledt <- rep(0,n)
- sw <- sqrt(weights)
- yw <- y*sw
- for (g in 1:G) {
- xw <- design*sw[g,]
- QR <- qr(xw)
- if(QR$rank<p) stop("weighted design matrix not of full rank for gene ",g)
- effects[,g] <- qr.qty(QR,yw[g,])
- unscaledt[g] <- effects[p,g]
- if(QR$qr[p,p]<0) unscaledt[g] <- -unscaledt[g]
- }
+ geneid <- as.character(geneid)
+ if(length(geneid)==1)
+ geneid <- as.character(y$genes[,geneid])
+ else
+ if(length(geneid) != nrow(y)) stop("geneid vector should be of length nrow(y)")
}
-# Standardized residuals
- U <- t(effects[-(1:(p-1)),,drop=FALSE])
+ .fryEffects(effects=Effects,index=index,geneid=geneid,sort=sort)
+}
+
+
+.fryEffects <- function(effects,index=NULL,geneid=rownames(effects),sort=TRUE)
+# fry given the effects matrix
+# Gordon Smyth and Goknur Giner
+# Created 30 January 2015. Last modified 9 May 2016
+{
+ G <- nrow(effects)
+ neffects <- ncol(effects)
+ df.residual <- neffects-1
+
+# Check index
+ if(is.null(index)) index <- list(set1=1L:G)
+ if(is.data.frame(index) || !is.list(index)) index <- list(set1=index)
+ if(is.null(names(index))) names(index) <- paste("set",1:nsets,sep="")
# Global statistics
nsets <- length(index)
+ if(nsets==0) stop("index is empty")
NGenes <- rep.int(0L,nsets)
- Direction <- rep.int("",nsets)
- PValue.Mixed <- PValue <- rep.int(0,nsets)
+ PValue.Mixed <- t.stat <- rep.int(0,nsets)
for (i in 1:nsets) {
iset <- index[[i]]
- USet <- U[iset,,drop=FALSE]
- NGenes[i] <- nrow(USet)
- MeanUSet <- colMeans(USet)
- t.stat <- MeanUSet[1L] / sqrt(mean(MeanUSet[-1L]^2L))
- if(t.stat>0) Direction[i] <- "Up" else Direction[i] <- "Down"
- PValue[i] <- 2*pt(-abs(t.stat),df=df.residual)
-
- if(NGenes[i]==1) {
- PValue.Mixed[i] <- PValue[i]
+ if(is.data.frame(iset)) {
+ if(ncol(iset)>1 && is.numeric(iset[,2])) {
+ iw <- iset[,2]
+ iset <- iset[,1]
+ if(is.factor(iset)) iset <- as.character(iset)
+ if(is.character(iset)) {
+ if(anyDuplicated(iset)) stop("Duplicate gene ids in set ",i)
+ j <- match(geneid,iset)
+ inset <- which(!is.na(j))
+ EffectsSet <- effects[inset,,drop=FALSE]
+ iw <- iw[j[inset]]
+ } else {
+ EffectsSet <- effects[iset,,drop=FALSE]
+ }
+ MeanEffectsSet <- iw %*% EffectsSet
+ } else {
+ stop("index ",i," is a data.frame but doesn't contain gene weights")
+ }
} else {
- SVD <- svd(USet,nu=0)
+ if(is.factor(iset)) iset <- as.character(iset)
+ if(is.character(iset)) iset <- which(geneid %in% iset)
+ EffectsSet <- effects[iset,,drop=FALSE]
+ MeanEffectsSet <- colMeans(EffectsSet)
+ }
+ t.stat[i] <- MeanEffectsSet[1] / sqrt(mean(MeanEffectsSet[-1]^2))
+ NGenes[i] <- nrow(EffectsSet)
+
+ if(NGenes[i]>1) {
+ SVD <- svd(EffectsSet,nu=0)
A <- SVD$d^2
d1 <- length(A)
d <- d1-1L
beta.mean <- 1/d1
beta.var <- d/d1/d1/(d1/2+1)
- Fobs <- (sum(USet[,1]^2)-A[d1]) / (A[1]-A[d1])
+ Fobs <- (sum(EffectsSet[,1]^2)-A[d1]) / (A[1]-A[d1])
Frb.mean <- (sum(A) * beta.mean - A[d1]) / (A[1]-A[d1])
COV <- matrix(-beta.var/d,d1,d1)
diag(COV) <- beta.var
@@ -134,6 +152,11 @@ fry.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),weights=N
}
}
+ Direction <- rep.int("Up",nsets)
+ Direction[t.stat<0] <- "Down"
+ PValue <- 2*pt(-abs(t.stat),df=df.residual)
+ PValue.Mixed[NGenes==1] <- PValue[NGenes==1]
+
# Add FDR
if(nsets>1) {
FDR <- p.adjust(PValue,method="BH")
@@ -144,11 +167,14 @@ fry.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),weights=N
}
rownames(tab) <- names(index)
-# Sort by p-value
- if(sort && nsets>1) {
- o <- order(tab$PValue,-tab$NGenes)
- tab <- tab[o,]
+# Sort results
+ if(is.logical(sort)) if(sort) sort <- "directional" else sort <- "none"
+ sort <- match.arg(sort,c("directional","mixed","none"))
+ if(sort=="none") return(tab)
+ if(sort=="directional") {
+ o <- order(tab$PValue,-tab$NGenes,tab$PValue.Mixed)
+ } else {
+ o <- order(tab$PValue.Mixed,-tab$NGenes,tab$PValue)
}
-
- tab
+ tab[o,,drop=FALSE]
}
diff --git a/R/geneset-roast.R b/R/geneset-roast.R
index a705747..6bfb763 100644
--- a/R/geneset-roast.R
+++ b/R/geneset-roast.R
@@ -13,222 +13,294 @@ function(object) print(object$p.value)
roast <- function(y,...) UseMethod("roast")
-roast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.statistic="mean",gene.weights=NULL,array.weights=NULL,weights=NULL,block=NULL,correlation,var.prior=NULL,df.prior=NULL,trend.var=FALSE,nrot=999,approx.zscore=TRUE,...)
+roast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),geneid=NULL,set.statistic="mean",gene.weights=NULL,var.prior=NULL,df.prior=NULL,nrot=999,approx.zscore=TRUE,...)
# Rotation gene set testing for linear models
# Gordon Smyth and Di Wu
-# Created 24 Apr 2008. Last modified 16 July 2015.
+# Created 24 Apr 2008. Last modified 9 May 2016.
{
-# Issue warning if extra arguments found
- dots <- names(list(...))
- if(length(dots)) warning("Extra arguments disregarded: ",sQuote(dots))
-
# Check index
- if(is.list(index)) return(mroast(y=y,index=index,design=design,contrast=contrast,set.statistic=set.statistic,gene.weights=gene.weights,array.weights=array.weights,weights=weights,block=block,correlation=correlation,var.prior=var.prior,df.prior=df.prior,trend.var=trend.var,nrot=nrot))
+ if(is.list(index)) return(mroast(y=y,index=index,design=design,contrast=contrast,set.statistic=set.statistic,gene.weights=gene.weights,var.prior=var.prior,df.prior=df.prior,nrot=nrot,approx.zscore=approx.zscore,...))
+
+# Partial matching of extra arguments
+ Dots <- list(...)
+ PossibleArgs <- c("array.weights","weights","block","correlation","trend","robust","winsor.tail.p")
+ if(!is.null(names(Dots))) {
+ i <- pmatch(names(Dots),PossibleArgs)
+ names(Dots) <- PossibleArgs[i]
+ }
-# Extract components from y
- y <- getEAWP(y)
- ngenes <- nrow(y$exprs)
- n <- ncol(y$exprs)
+# Defaults for extra arguments
+ if(is.null(Dots$trend)) Dots$trend <- FALSE
+ if(is.null(Dots$robust)) Dots$robust <- FALSE
+ if(Dots$robust & is.null(Dots$winsor.tail.p)) Dots$winsor.tail.p <- c(0.05,0.1)
+ if(!is.null(Dots$block) & is.null(Dots$correlation)) stop("correlation must be set")
+
+# Covariate for trended eBayes
+ covariate <- NULL
+ if(Dots$trend) covariate <- rowMeans(as.matrix(y))
+
+# Compute effects matrix (with df.residual+1 columns)
+ Effects <- .lmEffects(y=y,design=design,contrast=contrast,
+ array.weights=Dots$array.weights,
+ weights=Dots$weights,
+ block=Dots$block,
+ correlation=Dots$correlation)
+ ngenes <- nrow(Effects)
+
+# Empirical Bayes posterior variances
+ s2 <- rowMeans(Effects[,-1L,drop=FALSE]^2)
+ df.residual <- ncol(Effects)-1L
+ if(is.null(var.prior) || is.null(df.prior)) {
+ sv <- squeezeVar(s2,df=df.residual,covariate=covariate,robust=Dots$robust,winsor.tail.p=Dots$winsor.tail.p)
+ var.prior <- sv$var.prior
+ df.prior <- sv$df.prior
+ var.post <- sv$var.post
+ } else {
+ var.post <- .squeezeVar(var=s2,df=df.residual,var.prior=var.prior,df.prior=df.prior)
+ }
-# Check index
- if(is.null(index)) index <- rep.int(TRUE,ngenes)
-
-# Check design
- if(is.null(design)) design <- y$design
- if(is.null(design))
- design <- matrix(1,n,1)
- else {
- design <- as.matrix(design)
- if(mode(design) != "numeric") stop("design must be a numeric matrix")
+# Check geneid (can be a vector of gene IDs or an annotation column name)
+ if(is.null(geneid)) {
+ geneid <- rownames(Effects)
+ } else {
+ geneid <- as.character(geneid)
+ if(length(geneid)==1)
+ geneid <- as.character(y$genes[,geneid])
+ else
+ if(length(geneid) != ngenes) stop("geneid vector should be of length nrow(y)")
}
- if(nrow(design) != n) stop("row dimension of design matrix must match column dimension of data")
- ne <- nonEstimable(design)
- if(!is.null(ne)) cat("Coefficients not estimable:",paste(ne,collapse=" "),"\n")
- p <- ncol(design)
- p0 <- p-1L
- d <- n-p
+# Subset to gene set
+ if(!is.null(index)) {
+ if(is.factor(index)) index <- as.character(index)
+ if(is.character(index)) index <- which(index %in% geneid)
+ Effects <- Effects[index,,drop=FALSE]
+ if(length(var.prior)>1) var.prior <- var.prior[index]
+ if(length(df.prior)>1) df.prior <- df.prior[index]
+ var.post <- var.post[index]
+ }
+ NGenesInSet <- nrow(Effects)
-# Check set.statistic
- set.statistic <- match.arg(set.statistic,c("mean","floormean","mean50","msq"))
+# length(gene.weights) can nrow(y) or NGenesInSet
+ lgw <- length(gene.weights)
+ if(!(lgw %in% c(0L,NGenesInSet,ngenes))) stop("length of gene.weights doesn't agree with number of genes or size of set")
+ if(lgw > NGenesInSet) gene.weights <- gene.weights[index]
-# Check var.prior and df.prior
- if(!is.null(var.prior) && var.prior<0) stop("var.prior must be non-negative")
- if(!is.null(df.prior) && df.prior<0) stop("df.prior must be non-negative")
+ .roastEffects(Effects,gene.weights=gene.weights,set.statistic=set.statistic,var.prior=var.prior,df.prior=df.prior,var.post=var.post,nrot=nrot,approx.zscore=approx.zscore)
+}
+
+
+mroast <- function(y,...) UseMethod("mroast")
-# Check array weights
- if(!is.null(array.weights)) {
- if(length(array.weights) != n) stop("Length of array.weights doesn't match number of array")
- if(any(array.weights <= 0)) stop("array.weights must be positive")
+mroast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),geneid=NULL,set.statistic="mean",gene.weights=NULL,var.prior=NULL,df.prior=NULL,nrot=999,approx.zscore=TRUE,adjust.method="BH",midp=TRUE,sort="directional",...)
+# Rotation gene set testing with multiple sets
+# Gordon Smyth and Di Wu
+# Created 28 Jan 2010. Last revised 9 May 2016.
+{
+# Partial matching of extra arguments
+ Dots <- list(...)
+ PossibleArgs <- c("array.weights","weights","block","correlation","trend","robust","winsor.tail.p")
+ if(!is.null(names(Dots))) {
+ i <- pmatch(names(Dots),PossibleArgs)
+ names(Dots) <- PossibleArgs[i]
}
-# Check observational weights
- if(is.null(weights) && is.null(array.weights)) weights <- y$weights
- if(!is.null(weights)) {
- weights <- as.matrix(weights)
- dimw <- dim(weights)
- if(dimw[1] != ngenes || dimw[2] != n) stop("weights must have same dimensions as y")
- if(any(weights <= 0)) stop("weights must be positive")
- if(!is.null(array.weights)) {
- weights <- .matvec(weights, array.weights)
- array.weights <- NULL
- }
+# Defaults for extra arguments
+ if(is.null(Dots$trend)) Dots$trend <- FALSE
+ if(is.null(Dots$robust)) Dots$robust <- FALSE
+ if(Dots$robust & is.null(Dots$winsor.tail.p)) Dots$winsor.tail.p <- c(0.05,0.1)
+ if(!is.null(Dots$block) & is.null(Dots$correlation)) stop("correlation must be set")
+
+# Covariate for trended eBayes
+ covariate <- NULL
+ if(Dots$trend) covariate <- rowMeans(as.matrix(y))
+
+# Compute effects matrix (with df.residual+1 columns)
+ Effects <- .lmEffects(y=y,design=design,contrast=contrast,
+ array.weights=Dots$array.weights,
+ weights=Dots$weights,
+ block=Dots$block,
+ correlation=Dots$correlation)
+ ngenes <- nrow(Effects)
+ df.residual <- ncol(Effects)-1L
+ geneid <- rownames(Effects)
+
+# Empirical Bayes posterior variances
+ s2 <- rowMeans(Effects[,-1L,drop=FALSE]^2)
+ df.residual <- ncol(Effects)-1L
+ if(is.null(var.prior) || is.null(df.prior)) {
+ sv <- squeezeVar(s2,df=df.residual,covariate=covariate,robust=Dots$robust,winsor.tail.p=Dots$winsor.tail.p)
+ var.prior <- sv$var.prior
+ df.prior <- sv$df.prior
+ var.post <- sv$var.post
+ } else {
+ var.post <- .squeezeVar(var=s2,df=df.residual,var.prior=var.prior,df.prior=df.prior)
}
-# Reduce to numeric expression matrix
- y <- y$exprs
+# Check index
+ if(is.null(index)) index <- rep_len(TRUE,length.out=ngenes)
+ if(!is.list(index)) index <- list(set = index)
+ nsets <- length(index)
+ if(nsets==0) stop("index is empty")
+ if(is.null(names(index))) names(index) <- paste("set",1:nsets,sep="")
+
+# Check gene.weights
+ lgw <- length(gene.weights)
+ if(lgw > 0L && lgw != ngenes) stop("gene.weights vector should be of length nrow(y)")
-# Divide out array weights
- if(!is.null(array.weights)) {
- sw <- sqrt(array.weights)
- design <- design*sw
- y <- .matvec(y,sw)
- array.weights <- NULL
+# Check geneid (can be a vector of gene IDs or an annotation column name)
+ if(is.null(geneid)) {
+ geneid <- rownames(Effects)
+ } else {
+ geneid <- as.character(geneid)
+ if(length(geneid)==1)
+ geneid <- as.character(y$genes[,geneid])
+ else
+ if(length(geneid) != ngenes) stop("geneid vector should be of length nrow(y)")
}
-# Divide out block correlation
- if(!is.null(block)) {
- block <- as.vector(block)
- if (length(block) != n) stop("Length of block does not match number of arrays")
- ub <- unique(block)
- nblocks <- length(ub)
- Z <- matrix(block,n,nblocks) == matrix(ub,n,nblocks,byrow = TRUE)
- cormatrix <- Z %*% (correlation * t(Z))
- diag(cormatrix) <- 1
- R <- chol(cormatrix)
- y <- t(backsolve(R, t(y), transpose = TRUE))
- dn <- dimnames(design)
- design <- backsolve(R, design, transpose = TRUE)
- dimnames(design) <- dn
- }
-
-# Check contrast
- if(is.character(contrast)) {
- if(length(contrast)>1L) {
- warning("using only first entry for contrast")
- contrast <- contrast[1]
+# Containers for genewise results
+ pv <- adjpv <- active <- array(0,c(nsets,4),dimnames=list(names(index),c("Down","Up","UpOrDown","Mixed")))
+ if(nsets<1L) return(pv)
+ NGenes <- rep_len(0L,length.out=nsets)
+
+# Call roast for each set
+ s20 <- var.prior
+ d0 <- df.prior
+ for(i in 1:nsets) {
+ g <- index[[i]]
+ if(is.data.frame(g)) {
+ if(ncol(g)>1 && is.numeric(g[,2])) {
+ gw <- g[,2]
+ g <- g[,1]
+ if(is.factor(g)) g <- as.character(g)
+ if(is.character(g)) {
+ if(anyDuplicated(g)) stop("Duplicate gene ids in set: ",names(index[i]))
+ j <- match(geneid,g)
+ g <- which(!is.na(j))
+ gw <- gw[j[g]]
+ }
+ } else {
+ stop("index ",names(index[i])," is a data.frame but doesn't contain gene weights")
+ }
+ } else {
+ if(is.factor(g)) g <- as.character(g)
+ if(is.character(g)) g <- which(geneid %in% g)
+ gw <- gene.weights[g]
}
- contrast <- which(contrast==colnames(design))
- if(length(contrast)==0L) stop("coef ",contrast," not found")
+ E <- Effects[g,,drop=FALSE]
+ if(length(var.prior)>1) s20 <- var.prior[g]
+ if(length(df.prior)>1) d0 <- df.prior[g]
+ s2post <- var.post[g]
+ out <- .roastEffects(E,gene.weights=gw,set.statistic=set.statistic,var.prior=s20,df.prior=d0,var.post=s2post,nrot=nrot,approx.zscore=approx.zscore)
+ pv[i,] <- out$p.value$P.Value
+ active[i,] <- out$p.value$Active.Prop
+ NGenes[i] <- out$ngenes.in.set
}
- if(all(contrast == 0)) stop("contrast all zero")
-# Reform design matrix so that contrast is last coefficient
- if(length(contrast) == 1L) {
- contrast <- as.integer(contrast)
- if(contrast < p)
- X <- cbind(design[,-contrast,drop=FALSE],design[,contrast,drop=FALSE])
- else
- X <- design
- } else {
- if(length(contrast) != p) stop("length of contrast must match column dimension of design")
- X <- contrastAsCoef(design, contrast, first=FALSE)$design
+# P-values
+ Up <- pv[,"Up"] < pv[,"Down"]
+ Direction <- rep.int("Down",nsets); Direction[Up] <- "Up"
+ TwoSidedP2 <- pv[,"UpOrDown"]
+ MixedP2 <- pv[,"Mixed"]
+ if(midp) {
+ TwoSidedP2 <- TwoSidedP2 - 1/2/(nrot+1)
+ MixedP2 <- MixedP2 - 1/2/(nrot+1)
}
- qr <- qr(X)
- signc <- sign(qr$qr[p,p])
+# Output data.frame
+ tab <- data.frame(
+ NGenes=NGenes,
+ PropDown=active[,"Down"],
+ PropUp=active[,"Up"],
+ Direction=Direction,
+ PValue=pv[,"UpOrDown"],
+ FDR=p.adjust(TwoSidedP2,method="BH"),
+ PValue.Mixed=pv[,"Mixed"],
+ FDR.Mixed=p.adjust(MixedP2,method="BH"),
+ row.names=names(index),
+ stringsAsFactors=FALSE
+ )
- if(is.null(var.prior) || is.null(df.prior)) {
-# Fit model to all genes
- if(is.null(weights)) {
- effects <- qr.qty(qr,t(y))
- } else {
- ws <- sqrt(weights)
- effects <- matrix(0,n,ngenes)
- signc <- rep.int(0,ngenes)
- for (g in 1:ngenes) {
- wX <- X*ws[g,]
- wy <- y[g,]*ws[g,]
- qrX <- qr(wX)
- signc[g] <- sign(qrX$qr[p,p])
- effects[,g] <- qr.qty(qrX,wy)
- }
- signc <- signc[index]
- }
-# Estimate global parameters s0 and d0
- s2 <- colMeans(effects[-(1:p),,drop=FALSE]^2)
- if(trend.var) covariate <- rowMeans(y,na.rm=TRUE) else covariate <- NULL
- sv <- squeezeVar(s2,df=d,covariate=covariate)
- d0 <- sv$df.prior
- s02 <- sv$var.prior
- if(trend.var) s02 <- s02[index]
- effects <- effects[,index,drop=FALSE]
- sd.post <- sqrt(sv$var.post[index])
+# Mid p-values
+ if(midp) {
+ tab$FDR <- pmax(tab$FDR, pv[,"UpOrDown"])
+ tab$FDR.Mixed <- pmax(tab$FDR.Mixed, pv[,"Mixed"])
+ }
+
+# Sort results
+ if(is.logical(sort)) if(sort) sort <- "directional" else sort <- "none"
+ sort <- match.arg(sort,c("directional","mixed","none"))
+ if(sort=="none") return(tab)
+ if(sort=="directional") {
+ Prop <- pmax(tab$PropUp,tab$PropDown)
+ o <- order(tab$PValue,-Prop,-tab$NGenes,tab$PValue.Mixed)
} else {
- d0 <- df.prior
- s02 <- var.prior
- if(length(s02)>1) {
- names(s02) <- rownames(y)
- s02 <- s02[index]
- }
- y <- y[index,,drop=FALSE]
- if(is.null(weights)) {
- effects <- qr.qty(qr,t(y))
- } else {
- ws <- sqrt(weights[index,,drop=FALSE])
- nset <- nrow(y)
- effects <- matrix(0,n,nset)
- signc <- rep.int(0,nset)
- for (g in 1:nset) {
- wX <- X*ws[g,]
- wy <- y[g,]*ws[g,]
- qrX <- qr(wX)
- signc[g] <- sign(qrX$qr[p,p])
- effects[,g] <- qr.qty(qrX,wy)
- }
- }
- s2 <- colMeans(effects[-(1:p),,drop=FALSE]^2)
- if(is.finite(d0))
- sd.post <- sqrt( (d0*s02+d*s2)/(d0+d) )
- else
- sd.post <- sqrt(s02)
+ Prop <- tab$PropUp+tab$PropDown
+ o <- order(tab$PValue.Mixed,-Prop,-tab$NGenes,tab$PValue)
}
+ tab[o,,drop=FALSE]
+}
-# From here, all results are for set only
- nset <- ncol(effects)
- if(p0>0)
- Y <- effects[-(1:p0),,drop=FALSE]
- else
- Y <- effects
- YY <- colSums(Y^2)
- B <- Y[1,]
- modt <- signc*B/sd.post
-
- statobs <- p <- rep(0,4)
- names(statobs) <- names(p) <- c("down","up","upordown","mixed")
- statrot <- array(0,c(nrot,4),dimnames=list(NULL,names(p)))
-# Convert to z-scores
- modt <- zscoreT(modt,df=d0+d,approx=approx.zscore)
+.roastEffects <- function(effects,gene.weights=NULL,set.statistic="mean",var.prior,df.prior,var.post,nrot=999,approx.zscore=TRUE)
+# Rotation gene set testing, given effects matrix for one set
+# Rows are genes. First column is primary effect. Other columns are residual effects.
+# Gordon Smyth and Di Wu
+# Created 24 Apr 2008. Last modified 7 May 2016.
+{
+ nset <- nrow(effects)
+ neffects <- ncol(effects)
+ df.residual <- neffects-1
+ df.total <- df.prior+df.residual
+
+# Observed z-statistics
+ modt <- effects[,1]/sqrt(var.post)
+ modt <- zscoreT(modt,df=df.total,approx=approx.zscore)
+
+# Estimate active proportions
+ if(is.null(gene.weights)) {
+ a1 <- mean(modt > sqrt(2))
+ a2 <- mean(modt < -sqrt(2))
+ } else {
+ s <- sign(gene.weights)
+ ss <- sum(abs(s))
+ a1 <- sum(s*modt > sqrt(2)) / ss
+ a2 <- sum(s*modt < -sqrt(2)) / ss
+ }
-# Active proportions
- if(!is.null(gene.weights)) {
- lgw <- length(gene.weights)
- if(lgw > nset && lgw==ngenes) {
- gene.weights <- gene.weights[index]
+# Rotated primary effects (neffects by nrot)
+ R <- matrix(rnorm(nrot*neffects),nrot,neffects)
+ R <- R/sqrt(rowSums(R^2))
+ Br <- tcrossprod(effects,R)
+
+# Moderated rotated variances
+ FinDf <- is.finite(df.prior)
+ if(all(FinDf)) {
+ s2r <- (rowSums(effects^2)-Br^2) / df.residual
+ s2postr <- (df.prior*var.prior+df.residual*s2r) / df.total
+ } else {
+ if(any(FinDf)) {
+ s2postr <- s2r <- (rowSums(effects^2)-Br^2) / df.residual
+ if(length(var.prior)>1L) s20 <- var.prior[FinDf] else s20 <- var.prior
+ s2postr[FinDf,] <- (df.prior[FinDf]*s20+df.residual*s2r[FinDf,]) / df.total[FinDf]
} else {
- if(lgw != nset) stop("length of gene.weights disagrees with size of set")
+ s2postr <- var.prior
}
- s <- sign(gene.weights)
- r1 <- mean(s*modt > sqrt(2))
- r2 <- mean(s*modt < -sqrt(2))
- } else {
- r1 <- mean(modt > sqrt(2))
- r2 <- mean(modt < -sqrt(2))
}
-# Random rotations
- R <- matrix(rnorm(nrot*(d+1)),nrot,d+1)
- R <- R/sqrt(rowSums(R^2))
- Br <- R %*% Y
- s2r <- (matrix(YY,nrot,nset,byrow=TRUE)-Br^2)/d
- if(is.finite(d0))
- sdr.post <- sqrt((d0*s02+d*s2r)/(d0+d))
- else
- sdr.post <- sqrt(s02)
- modtr <- signc*Br/sdr.post
- modtr <- zscoreT(modtr,df=d0+d,approx=approx.zscore)
+# Rotated z-statistics
+ modtr <- Br/sqrt(s2postr)
+ modtr <- zscoreT(modtr,df=df.total,approx=approx.zscore)
+# Setup matrices to hold output results
+ statobs <- p <- rep_len(0,length.out=4)
+ names(statobs) <- names(p) <- c("down","up","upordown","mixed")
+ statrot <- array(0,c(nrot,4),dimnames=list(NULL,names(p)))
+
+# Compute set statistics and p-values
+
+ set.statistic <- match.arg(set.statistic,c("mean","floormean","mean50","msq"))
switch(set.statistic,
"mean" = {
# Observed statistics
@@ -238,11 +310,11 @@ roast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.sta
statobs["up"] <- m
statobs["mixed"] <- mean(abs(modt))
# Simulated statistics
- if(!is.null(gene.weights)) modtr <- t(gene.weights*t(modtr))
- m <- rowMeans(modtr)
+ if(!is.null(gene.weights)) modtr <- gene.weights*modtr
+ m <- colMeans(modtr)
statrot[,"down"] <- -m
statrot[,"up"] <- m
- statrot[,"mixed"] <- rowMeans(abs(modtr))
+ statrot[,"mixed"] <- colMeans(abs(modtr))
# p-values
p["down"] <- sum(statrot[,c("down","up")] > statobs["down"])
p["up"] <- sum(statrot[,c("down","up")] > statobs["up"])
@@ -266,15 +338,15 @@ roast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.sta
# Simulated statistics
amodtr <- pmax(abs(modtr),chimed)
if(!is.null(gene.weights)) {
- amodtr <- t(gene.weights*t(amodtr))
- modtr <- t(gene.weights*t(modtr))
+ amodtr <- gene.weights*amodtr
+ modtr <- gene.weights*modtr
}
- statrot[,"down"] <- rowMeans(pmax(-modtr,0))
- statrot[,"up"] <- rowMeans(pmax(modtr,0))
+ statrot[,"down"] <- colMeans(pmax(-modtr,0))
+ statrot[,"up"] <- colMeans(pmax(modtr,0))
i <- statrot[,"up"] > statrot[,"down"]
statrot[i,"upordown"] <- statrot[i,"up"]
statrot[!i,"upordown"] <- statrot[!i,"down"]
- statrot[,"mixed"] <- rowMeans(amodtr)
+ statrot[,"mixed"] <- colMeans(amodtr)
# p-values
p["down"] <- sum(statrot[,c("down","up")] > statobs["down"])
p["up"] <- sum(statrot[,c("down","up")] > statobs["up"])
@@ -299,14 +371,14 @@ roast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.sta
s <- sort(abs(modt),partial=half2)
statobs["mixed"] <- mean(s[half2:nset])
# Simulated statistics
- if(!is.null(gene.weights)) modtr <- t(gene.weights*t(modtr))
- for (g in 1L:nrot) {
- s <- sort(modtr[g,],partial=half2)
- statrot[g,"down"] <- -mean(s[1:half1])
- statrot[g,"up"] <- mean(s[half2:nset])
- statrot[g,"upordown"] <- max(statrot[g,c("down","up")])
- s <- sort(abs(modtr[g,]),partial=half2)
- statrot[g,"mixed"] <- mean(s[half2:nset])
+ if(!is.null(gene.weights)) modtr <- gene.weights*modtr
+ for (r in 1L:nrot) {
+ s <- sort(modtr[,r],partial=half2)
+ statrot[r,"down"] <- -mean(s[1:half1])
+ statrot[r,"up"] <- mean(s[half2:nset])
+ statrot[r,"upordown"] <- max(statrot[r,c("down","up")])
+ s <- sort(abs(modtr[,r]),partial=half2)
+ statrot[r,"mixed"] <- mean(s[half2:nset])
}
# p-values
p["down"] <- sum(statrot[,c("down","up")] > statobs["down"])
@@ -330,14 +402,14 @@ roast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.sta
# Simulated statistics
if(!is.null(gene.weights)) {
gene.weights <- sqrt(abs(gene.weights))
- modtr <- t(gene.weights*t(modtr))
+ modtr <- gene.weights*modtr
}
- statrot[,"down"] <- rowMeans(pmax(-modtr,0)^2)
- statrot[,"up"] <- rowMeans(pmax(modtr,0)^2)
+ statrot[,"down"] <- colMeans(pmax(-modtr,0)^2)
+ statrot[,"up"] <- colMeans(pmax(modtr,0)^2)
i <- statrot[,"up"] > statrot[,"down"]
statrot[i,"upordown"] <- statrot[i,"up"]
statrot[!i,"upordown"] <- statrot[!i,"down"]
- statrot[,"mixed"] <- rowMeans(modtr^2)
+ statrot[,"mixed"] <- colMeans(modtr^2)
# p-values
p["down"] <- sum(statrot[,c("down","up")] > statobs["down"])
p["up"] <- sum(statrot[,c("down","up")] > statobs["up"])
@@ -347,160 +419,7 @@ roast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.sta
})
# Output
- out <- data.frame(c(r2,r1,max(r1,r2),r1+r2),p)
+ out <- data.frame(c(a2,a1,max(a1,a2),a1+a2),p)
dimnames(out) <- list(c("Down","Up","UpOrDown","Mixed"),c("Active.Prop","P.Value"))
- new("Roast",list(p.value=out,var.prior=s02,df.prior=d0,ngenes.in.set=nset))
-}
-
-mroast <- function(y,...) UseMethod("mroast")
-
-mroast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.statistic="mean",gene.weights=NULL,array.weights=NULL,weights=NULL,block=NULL,correlation,var.prior=NULL,df.prior=NULL,trend.var=FALSE,nrot=999,approx.zscore=TRUE,adjust.method="BH",midp=TRUE,sort="directional",...)
-# Rotation gene set testing with multiple sets
-# Gordon Smyth and Di Wu
-# Created 28 Jan 2010. Last revised 7 Oct 2015.
-{
-# Extract components from y
- y <- getEAWP(y)
- ngenes <- nrow(y$exprs)
- n <- ncol(y$exprs)
-
-# Check index
- if(is.null(index)) index <- rep(TRUE,ngenes)
- if(!is.list(index)) index <- list(set = index)
- nsets <- length(index)
- if(nsets==0) stop("index is empty")
- if(is.null(names(index))) names(index) <- paste("set",1:nsets,sep="")
-
-# Check design matrix
- if(is.null(design)) design <- y$design
- if(is.null(design))
- design <- matrix(1,n,1)
- else {
- design <- as.matrix(design)
- if(mode(design) != "numeric") stop("design must be a numeric matrix")
- }
- if(nrow(design) != n) stop("row dimension of design matrix must match column dimension of data")
- ne <- nonEstimable(design)
- if(!is.null(ne)) cat("Coefficients not estimable:",paste(ne,collapse=" "),"\n")
-
-# Check gene.weights
- if(!is.null(gene.weights)) if(length(gene.weights) != ngenes) stop("gene.weights must have length equal to nrow(y)")
-
-# Check array.weights
- if(!is.null(array.weights)) {
- if(length(array.weights) != n) stop("array.weights wrong length")
- if(any(array.weights <= 0)) stop("array.weights must be positive")
- }
-
-# Check weights
- if(is.null(weights) && is.null(array.weights)) weights <- y$weights
- if(!is.null(weights)) {
- weights <- as.matrix(weights)
- dimw <- dim(weights)
- if(dimw[1] != ngenes || dimw[2] != n) stop("weights must have same dimensions as y")
- if(any(weights <= 0)) stop("weights must be positive")
- if(!is.null(array.weights)) {
- weights <- .matvec(weights,array.weights)
- array.weights <- NULL
- }
- }
-
-# Reduce to numeric expression matrix
- y <- y$exprs
-
-# Divide out array.weights
- if(!is.null(array.weights)) {
- sw <- sqrt(array.weights)
- design <- design*sw
- y <- .matvec(y,sw)
- array.weights <- NULL
- }
-
-# Divide out block correlation
- if(!is.null(block)) {
- block <- as.vector(block)
- if (length(block) != n) stop("Length of block does not match number of arrays")
- ub <- unique(block)
- nblocks <- length(ub)
- Z <- matrix(block,n,nblocks) == matrix(ub,n,nblocks,byrow = TRUE)
- cormatrix <- Z %*% (correlation * t(Z))
- diag(cormatrix) <- 1
- R <- chol(cormatrix)
- y <- t(backsolve(R, t(y), transpose = TRUE))
- design <- backsolve(R, design, transpose = TRUE)
- block <- NULL
- }
-
-# Estimate var.prior and df.prior if not preset
- if(is.null(var.prior) || is.null(df.prior)) {
- fit <- lmFit(y,design=design,weights=weights)
- if(trend.var) {
- covariate <- fit$Amean
- if(is.null(covariate)) covariate <- rowMeans(y)
- } else {
- covariate=NULL
- }
- sv <- squeezeVar(fit$sigma^2,df=fit$df.residual,covariate=covariate)
- var.prior <- sv$var.prior
- df.prior <- sv$df.prior
- }
-
- pv <- adjpv <- active <- array(0,c(nsets,4),dimnames=list(names(index),c("Down","Up","UpOrDown","Mixed")))
- NGenes <- rep(0,nsets)
- if(nsets<1) return(pv)
- for(i in 1:nsets) {
- out <- roast(y=y,index=index[[i]],design=design,contrast=contrast,set.statistic=set.statistic,gene.weights=gene.weights,weights=weights,var.prior=var.prior,df.prior=df.prior,nrot=nrot,approx.zscore=approx.zscore,...)
- pv[i,] <- out$p.value$P.Value
- active[i,] <- out$p.value$Active.Prop
- NGenes[i] <- out$ngenes.in.set
- }
-
-# Use mid-p-values or ordinary p-values?
-# pv2 <- pv
-# if(midp) pv2 <- pv2-1/2/(nrot+1)
-# adjpv[,"Down"] <- p.adjust(pv2[,"Down"], method=adjust.method)
-# adjpv[,"Up"] <- p.adjust(pv2[,"Up"], method=adjust.method)
-# adjpv[,"Mixed"] <- p.adjust(pv2[,"Mixed"], method=adjust.method)
-# list(P.Value=pv, Adj.P.Value=adjpv, Active.Proportion=active)
-
-# New-style output
- Up <- pv[,"Up"] < pv[,"Down"]
- Direction <- rep.int("Down",nsets); Direction[Up] <- "Up"
- TwoSidedP2 <- pv[,"UpOrDown"]
- MixedP2 <- pv[,"Mixed"]
- if(midp) {
- TwoSidedP2 <- TwoSidedP2 - 1/2/(nrot+1)
- MixedP2 <- MixedP2 - 1/2/(nrot+1)
- }
-
- tab <- data.frame(
- NGenes=NGenes,
- PropDown=active[,"Down"],
- PropUp=active[,"Up"],
- Direction=Direction,
- PValue=pv[,"UpOrDown"],
- FDR=p.adjust(TwoSidedP2,method="BH"),
- PValue.Mixed=pv[,"Mixed"],
- FDR.Mixed=p.adjust(MixedP2,method="BH"),
- row.names=names(index),
- stringsAsFactors=FALSE
- )
-
- if(midp) {
- tab$FDR <- pmax(tab$FDR, pv[,"UpOrDown"])
- tab$FDR.Mixed <- pmax(tab$FDR.Mixed, pv[,"Mixed"])
- }
-
-# Sort by p-value
- sort <- match.arg(sort,c("directional","mixed","none"))
- if(sort=="none") return(tab)
- if(sort=="directional") {
- Prop <- pmax(tab$PropUp,tab$PropDown)
- o <- order(tab$PValue,-Prop,-tab$NGenes,tab$PValue.Mixed)
- } else {
- Prop <- tab$PropUp+tab$PropDown
- o <- order(tab$PValue.Mixed,-Prop,-tab$NGenes,tab$PValue)
- }
- tab[o,,drop=FALSE]
+ new("Roast",list(p.value=out,ngenes.in.set=nset))
}
-
diff --git a/R/geneset-romer.R b/R/geneset-romer.R
index 093b4b0..df988dc 100644
--- a/R/geneset-romer.R
+++ b/R/geneset-romer.R
@@ -1,10 +1,10 @@
## ROMER.R
romer <- function(y,...) UseMethod("romer")
-romer.default <- function(y,index,design,contrast=ncol(design),array.weights=NULL,block=NULL,correlation=NULL,set.statistic="mean",nrot=9999,shrink.resid=TRUE,...)
+romer.default <- function(y,index,design=NULL,contrast=ncol(design),array.weights=NULL,block=NULL,correlation=NULL,set.statistic="mean",nrot=9999,shrink.resid=TRUE,...)
# rotation mean-rank version of GSEA (gene set enrichment analysis) for linear models
# Gordon Smyth and Yifang Hu
-# 27 March 2009. Last modified 23 June 2015.
+# 27 March 2009. Last modified 3 May 2015.
{
# Issue warning if extra arguments found
dots <- names(list(...))
@@ -22,10 +22,17 @@ romer.default <- function(y,index,design,contrast=ncol(design),array.weights=NUL
SetSizes <- unlist(lapply(index,length))
# Check design
- design <- as.matrix(design)
+ if(is.null(design))
+ stop("design matrix not specified")
+ else {
+ design <- as.matrix(design)
+ if(mode(design) != "numeric") stop("design must be a numeric matrix")
+ }
+ if(nrow(design) != n) stop("row dimension of design matrix must match column dimension of data")
+ ne <- nonEstimable(design)
+ if(!is.null(ne)) cat("Coefficients not estimable:",paste(ne,collapse=" "),"\n")
p <- ncol(design)
if(p<2L) stop("design needs at least two columns")
- if(nrow(design) != n) stop("row dimension of design matrix must match column dimension of data")
p0 <- p-1L
d <- n-p
diff --git a/R/goana.R b/R/goana.R
index 114e88f..f0fd1e3 100644
--- a/R/goana.R
+++ b/R/goana.R
@@ -77,7 +77,7 @@ goana.MArrayLM <- function(de, coef = ncol(de), geneid = rownames(de), FDR = 0.0
goana.default <- function(de, universe = NULL, species = "Hs", prior.prob = NULL, covariate=NULL, plot=FALSE, ...)
# Gene ontology analysis of DE genes
# Gordon Smyth and Yifang Hu
-# Created 20 June 2014. Last modified 27 March 2015.
+# Created 20 June 2014. Last modified 10 April 2016.
{
# Ensure de is a list
if(!is.list(de)) de <- list(DE = de)
@@ -98,7 +98,7 @@ goana.default <- function(de, universe = NULL, species = "Hs", prior.prob = NULL
names(de) <- makeUnique(NAME)
# Select species
- species <- match.arg(species, c("Hs", "Mm", "Rn", "Dm"))
+ species <- match.arg(species, c("Hs", "Mm", "Rn", "Dm", "Pt"))
# Fit trend in DE with respect to the covariate, combining all de lists
if(!is.null(covariate)) {
@@ -112,17 +112,42 @@ goana.default <- function(de, universe = NULL, species = "Hs", prior.prob = NULL
if(plot) barcodeplot(covariate, index=(isDE==1), worm=TRUE, span.worm=span)
}
-# Load package of GO terms
- if(!suppressPackageStartupMessages(require("GO.db", character.only = TRUE))) stop("GO.db package not available")
-
-# Load species annotation package
- DB <- paste("org", species, "eg", "db", sep = ".")
- if(!suppressPackageStartupMessages(require(DB, character.only = TRUE))) stop(DB,"package not available")
+# Get access to package of GO terms
+ suppressPackageStartupMessages(OK <- requireNamespace("GO.db",quietly=TRUE))
+ if(!OK) stop("GO.db package required but not installed (or can't be loaded)")
+
+# Get access to required annotation functions
+ suppressPackageStartupMessages(OK <- requireNamespace("AnnotationDbi",quietly=TRUE))
+ if(!OK) stop("AnnotationDbi package required but not installed (or can't be loaded)")
+
+# Get gene-GOterm mappings
+ switch(species,
+ Hs = {
+ suppressPackageStartupMessages(OK <- requireNamespace("org.Hs.eg.db",quietly=TRUE))
+ if(!OK) stop("org.Hs.eg.db package required but not installed (or can't be loaded)")
+ GO2ALLEGS <- org.Hs.eg.db::org.Hs.egGO2ALLEGS
+ }, Mm = {
+ suppressPackageStartupMessages(OK <- requireNamespace("org.Mm.eg.db",quietly=TRUE))
+ if(!OK) stop("org.Mm.eg.db package required but not installed (or can't be loaded)")
+ GO2ALLEGS <- org.Mm.eg.db::org.Mm.egGO2ALLEGS
+ }, Rn = {
+ suppressPackageStartupMessages(OK <- requireNamespace("org.Rn.eg.db",quietly=TRUE))
+ if(!OK) stop("org.Rn.eg.db package required but not installed (or can't be loaded)")
+ GO2ALLEGS <- org.Rn.eg.db::org.Rn.egGO2ALLEGS
+ }, Dm = {
+ suppressPackageStartupMessages(OK <- requireNamespace("org.Dm.eg.db",quietly=TRUE))
+ if(!OK) stop("org.Dm.eg.db package required but not installed (or can't be loaded)")
+ GO2ALLEGS <- org.Dm.eg.db::org.Dm.egGO2ALLEGS
+ }, Pt = {
+ suppressPackageStartupMessages(OK <- requireNamespace("org.Pt.eg.db",quietly=TRUE))
+ if(!OK) stop("org.Pt.eg.db package required but not installed (or can't be loaded)")
+ GO2ALLEGS <- org.Pt.eg.db::org.Pt.egGO2ALLEGS
+ }
+ )
-# Get gene-GOterm mappings, and remove duplicate entries
- GO2ALLEGS <- paste("org", species, "egGO2ALLEGS", sep = ".")
+# Convert gene-GOterm mappings to data.frame and remove duplicate entries
if(is.null(universe)) {
- EG.GO <- AnnotationDbi::toTable(get(GO2ALLEGS))
+ EG.GO <- AnnotationDbi::toTable(GO2ALLEGS)
d <- duplicated(EG.GO[,c("gene_id", "go_id", "Ontology")])
EG.GO <- EG.GO[!d, ]
universe <- unique(EG.GO$gene_id)
@@ -140,7 +165,6 @@ goana.default <- function(de, universe = NULL, species = "Hs", prior.prob = NULL
}
universe <- universe[!dup]
- GO2ALLEGS <- get(GO2ALLEGS)
m <- match(AnnotationDbi::Lkeys(GO2ALLEGS),universe,0L)
universe <- universe[m]
if(!is.null(prior.prob)) prior.prob <- prior.prob[m]
@@ -185,7 +209,7 @@ goana.default <- function(de, universe = NULL, species = "Hs", prior.prob = NULL
W <- AVE.PW*(Total-S[,"N"])/(PW.ALL-S[,"N"]*AVE.PW)
# Wallenius' noncentral hypergeometric test
- if(!requireNamespace("BiasedUrn",quietly=TRUE)) stop("BiasedUrn package required but is not available")
+ if(!requireNamespace("BiasedUrn",quietly=TRUE)) stop("BiasedUrn package required but is not installed (or can't be loaded)")
for(j in 1:nsets) for(i in 1:nrow(S))
P[i,j] <- BiasedUrn::pWNCHypergeo(S[i,1+j], S[i,"N"], Total-S[i,"N"], TotalDE[[j]], W[i],lower.tail=FALSE) + BiasedUrn::dWNCHypergeo(S[i,1+j], S[i,"N"], Total-S[i,"N"], TotalDE[[j]], W[i])
S <- S[,-ncol(S)]
diff --git a/R/kegga.R b/R/kegga.R
index 0ef60f7..57d743f 100644
--- a/R/kegga.R
+++ b/R/kegga.R
@@ -74,23 +74,23 @@ kegga.MArrayLM <- function(de, coef = ncol(de), geneid = rownames(de), FDR = 0.0
kegga(de=DEGenes, universe = universe, ...)
}
-kegga.default <- function(de, universe = NULL, species = "Hs", prior.prob = NULL, covariate=NULL, plot=FALSE, ...)
+kegga.default <- function(de, universe=NULL, species="Hs", species.KEGG=NULL, convert=FALSE, gene.pathway=NULL, pathway.names = NULL, prior.prob=NULL, covariate=NULL, plot=FALSE, ...)
# KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis of DE genes
# Gordon Smyth and Yifang Hu
-# Created 18 May 2015. Modified 4 August 2015.
+# Created 18 May 2015. Modified 11 Feb 2016.
{
# Ensure de is a list
if(!is.list(de)) de <- list(DE = de)
+ nsets <- length(de)
# Stop if components of de are not vectors
if(!all(vapply(de,is.vector,TRUE))) stop("components of de should be vectors")
# Ensure gene IDs are of character mode
- de <- lapply(de, as.character)
+ for (s in 1:nsets) de[[s]] <- unique(as.character(de[[s]]))
if(!is.null(universe)) universe <- as.character(universe)
# Ensure all gene sets have unique names
- nsets <- length(de)
names(de) <- trimWhiteSpace(names(de))
NAME <- names(de)
i <- which(NAME == "" | is.na(NAME))
@@ -98,125 +98,140 @@ kegga.default <- function(de, universe = NULL, species = "Hs", prior.prob = NULL
names(de) <- makeUnique(NAME)
# Select species
- species <- match.arg(species, c("Hs", "Mm", "Rn", "Dm"))
+ if(is.null(species.KEGG)) {
+ species <- match.arg(species, c("Hs", "Mm", "Rn", "Dm", "Pt"))
+# Convert from Bioconductor to KEGG species codes
+ species.KEGG <- switch(species, "Hs"="hsa", "Dm"="dme", "Mm"="mmu", "Rn"="rno", "Pt"="ptr")
+ }
+
+# prior.prob and covariate must have same length as universe
+ if(is.null(universe)) {
+ if(!is.null(prior.prob)) stop("prior.prob ignored unless universe is specified.")
+ if(!is.null(covariate)) stop("covariate ignored unless universe is specified.")
+ prior.prob <- covariate <- NULL
+ } else {
+ lu <- length(universe)
+ if(!lu) stop("No genes in universe")
+ if(!is.null(prior.prob) && length(prior.prob)!=lu) stop("universe and prior.prob must have same length")
+ if(!is.null(covariate) && length(covariate)!=lu) stop("universe and covariate must have same length")
+ }
# Fit trend in DE with respect to the covariate, combining all de lists
if(!is.null(covariate)) {
+ if(!is.null(prior.prob)) warning("prior.prob being recomputed from covariate")
covariate <- as.numeric(covariate)
- if(length(covariate) != length(covariate)) stop("universe and covariate must have same length")
isDE <- as.numeric(universe %in% unlist(de))
o <- order(covariate)
prior.prob <- covariate
span <- approx(x=c(20,200),y=c(1,0.5),xout=sum(isDE),rule=2)$y
prior.prob[o] <- tricubeMovingAverage(isDE[o],span=span)
- if(plot) barcodeplot(covariate, index=(isDE==1), worm=TRUE, span.worm=span)
+ if(plot) barcodeplot(covariate, index=as.logical(isDE), worm=TRUE, span.worm=span, main="DE status vs covariate")
}
-# Enable REST-style online access to KEGG pathways
- if(!requireNamespace("KEGGREST",quietly=TRUE)) stop("KEGGREST required but is not available")
-
-# Convert to KEGG organism/species codes
- organism <- switch(species, "Hs"="hsa", "Dm"="dme", "Mm"="mmu", "Rn"="rno")
+# Get pathway annotation
+ if(is.null(gene.pathway))
+ EG.KEGG <- getGeneKEGGLinks(species.KEGG, convert=convert)
+ else {
+ EG.KEGG <- gene.pathway
+ d <- dim(EG.KEGG)
+ if(is.null(d)) stop("gene.pathway must be data.frame or matrix")
+ if(d[2] < 2) stop("gene.pathway must have at least 2 columns")
+ isna <- rowSums(is.na(EG.KEGG[,1:2])) > 0.5
+ EG.KEGG <- EG.KEGG[!isna,]
+ }
-# Get gene-KEGG mappings, and remove duplicate entries
+# Make sure that universe matches annotated genes
if(is.null(universe)) {
-
- path <- KEGGREST::keggLink("pathway", organism)
- path <- data.frame(kegg_id = names(path), path_id = path, stringsAsFactors = FALSE)
- if(organism == "dme") {
- EG <- KEGGREST::keggConv("dme", "ncbi-geneid")
- geneid <- names(EG)[match(path$kegg_id, EG)]
- geneid <- gsub("ncbi-geneid:", "", geneid)
- } else {
- geneid <- gsub(paste0(organism, ":"), "", path$kegg_id)
- }
- EG.KEGG <- data.frame(gene_id = geneid, path, stringsAsFactors = FALSE)
-
- universe <- unique(EG.KEGG$gene_id)
- universe <- as.character(universe)
+ universe <- as.character(unique(EG.KEGG[,1]))
} else {
+# Consolidate replicated entries in universe, averaging corresponding prior.probs
universe <- as.character(universe)
-
dup <- duplicated(universe)
if(!is.null(prior.prob)) {
- if(length(prior.prob)!=length(universe)) stop("length(prior.prob) must equal length(universe)")
prior.prob <- rowsum(prior.prob,group=universe,reorder=FALSE)
n <- rowsum(rep_len(1L,length(universe)),group=universe,reorder=FALSE)
prior.prob <- prior.prob/n
}
universe <- universe[!dup]
- path <- KEGGREST::keggLink("pathway", organism)
- path <- data.frame(kegg_id = names(path), path_id = path, stringsAsFactors = FALSE)
- if(organism == "dme") {
- EG <- KEGGREST::keggConv("dme", "ncbi-geneid")
- geneid <- names(EG)[match(path$kegg_id, EG)]
- geneid <- gsub("ncbi-geneid:", "", geneid)
- } else {
- geneid <- gsub(paste0(organism, ":"), "", path$kegg_id)
- }
- EG.KEGG <- data.frame(gene_id = geneid, path, stringsAsFactors = FALSE)
-
- m <- match(EG.KEGG$gene_id, universe)
+# Restrict universe to genes with annotation
+ m <- match(EG.KEGG[,1], universe)
+ InUni <- !is.na(m)
+ m <- m[InUni]
universe <- universe[m]
if(!is.null(prior.prob)) prior.prob <- prior.prob[m]
- EG.KEGG <- EG.KEGG[EG.KEGG$gene_id %in% universe,]
- }
-
- Total <- length(unique(EG.KEGG$gene_id))
- if(Total<1L) stop("No genes found in universe")
-# Check prior probabilities
- if(!is.null(prior.prob)) {
- if(length(prior.prob)!=length(universe)) stop("length(prior.prob) must equal length(universe)")
+# Restrict annotation to genes in universe
+ EG.KEGG <- EG.KEGG[InUni,]
}
-# Overlap with DE genes
- isDE <- lapply(de, function(x) EG.KEGG$gene_id %in% x)
- TotalDE <- lapply(isDE, function(x) length(unique(EG.KEGG$gene_id[x])))
- nDE <- length(isDE)
+ NGenes <- length(universe)
+ if(NGenes<1L) stop("No annotated genes found in universe")
- if(length(prior.prob)) {
- # Probability weight for each gene
- m <- match(EG.KEGG$gene_id, universe)
- PW2 <- list(prior.prob[m])
- X <- do.call(cbind, c(N=1, isDE, PW=PW2))
- } else
- X <- do.call(cbind, c(N=1, isDE))
+# Restrict DE genes to universe
+ nde <- rep_len(0L,nsets)
+ for (s in 1:nsets) {
+ de[[s]] <- de[[s]][de[[s]] %in% universe]
+ nde[s] <- length(de[[s]])
+ }
- group <- EG.KEGG$path_id
- S <- rowsum(X, group=group, reorder=FALSE)
+# Overlap pathways with DE genes
+ if(is.null(prior.prob)) {
+ X <- matrix(1,nrow(EG.KEGG),nsets+1)
+ colnames(X) <- c("N",names(de))
+ } else {
+ X <- matrix(1,nrow(EG.KEGG),nsets+2)
+ X[,nsets+2] <- prior.prob
+ colnames(X) <- c("N",names(de),"PP")
+ }
+ for (s in 1:nsets) X[,s+1] <- (EG.KEGG[,1] %in% de[[s]])
- P <- matrix(0, nrow = nrow(S), ncol = nsets)
+# Count #genes and #DE genes and sum prior.prob for each pathway
+ S <- rowsum(X, group=EG.KEGG[,2], reorder=FALSE)
+# Overlap tests
+ PValue <- matrix(0,nrow=nrow(S),ncol=nsets)
if(length(prior.prob)) {
# Calculate average prior prob for each set
- PW.ALL <- sum(prior.prob[universe %in% EG.KEGG$gene_id])
- AVE.PW <- S[,"PW"]/S[,"N"]
- W <- AVE.PW*(Total-S[,"N"])/(PW.ALL-S[,"N"]*AVE.PW)
+ SumPP <- sum(prior.prob)
+ AvePP.set <- S[,"PP"]/S[,"N"]
+ W <- AvePP.set*(NGenes-S[,"N"])/(SumPP-S[,"N"]*AvePP.set)
# Wallenius' noncentral hypergeometric test
if(!requireNamespace("BiasedUrn",quietly=TRUE)) stop("BiasedUrn package required but is not available")
- for(j in 1:nsets) for(i in 1:nrow(S))
- P[i,j] <- BiasedUrn::pWNCHypergeo(S[i,1+j], S[i,"N"], Total-S[i,"N"], TotalDE[[j]], W[i],lower.tail=FALSE) + BiasedUrn::dWNCHypergeo(S[i,1+j], S[i,"N"], Total-S[i,"N"], TotalDE[[j]], W[i])
+ for(j in 1:nsets) for(i in 1:nrow(S))
+ PValue[i,j] <- BiasedUrn::pWNCHypergeo(S[i,1+j], S[i,"N"], NGenes-S[i,"N"], nde[j], W[i], lower.tail=FALSE) + BiasedUrn::dWNCHypergeo(S[i,1+j], S[i,"N"], NGenes-S[i,"N"], nde[j], W[i])
+
+# Remove sum of prob column, not needed for output
S <- S[,-ncol(S)]
} else {
# Fisher's exact test
for(j in 1:nsets)
- P[,j] <- phyper(q=S[,1+j]-0.5,m=TotalDE[[j]],n=Total-TotalDE[[j]], k=S[,"N"],lower.tail=FALSE)
+ PValue[,j] <- phyper(q=S[,1+j]-0.5,m=nde[j],n=NGenes-nde[j], k=S[,"N"],lower.tail=FALSE)
+
+ }
+
+# Get list of pathway names
+ if(is.null(pathway.names))
+ PathID.PathName <- getKEGGPathwayNames(species.KEGG, remove.qualifier=TRUE)
+ else {
+ PathID.PathName <- pathway.names
+ d <- dim(PathID.PathName)
+ if(is.null(d)) stop("pathway.names must be data.frame or matrix")
+ if(d[2] < 2) stop("pathway.names must have at least 2 columns")
+ isna <- rowSums(is.na(PathID.PathName[,1:2])) > 0.5
+ PathID.PathName <- PathID.PathName[!isna,]
}
# Assemble output
g <- rownames(S)
- pathname <- KEGGREST::keggList("pathway")
- names(pathname) <- gsub("map", organism, names(pathname))
- m <- match(g, names(pathname))
- Results <- data.frame(Pathway = pathname[m], S, P, stringsAsFactors=FALSE)
+ m <- match(g, PathID.PathName[,1])
+ Results <- data.frame(Pathway = PathID.PathName[m,2], S, PValue, stringsAsFactors=FALSE)
rownames(Results) <- g
# Name p-value columns
@@ -225,10 +240,48 @@ kegga.default <- function(de, universe = NULL, species = "Hs", prior.prob = NULL
Results
}
+getGeneKEGGLinks <- function(species.KEGG="hsa", convert=FALSE)
+# Read pathway-gene mapping from KEGG website
+# Gordon Smyth
+# Created 7 Jan 2016. Last modified 11 Feb 2016.
+{
+ URL <- paste("http://rest.kegg.jp/link/pathway",species.KEGG,sep="/")
+ Path.Gene <- read.table(URL,sep="\t",quote="\"",fill=TRUE,comment.char="",stringsAsFactors=FALSE)
+ colnames(Path.Gene) <- c("GeneID","PathwayID")
+
+# Human, mouse, rat and chimp IDs are already Entrez, so don't need to convert
+ if(convert && species.KEGG %in% c("hsa","mmu","rno","ptr")) convert <- FALSE
+
+ if(convert) {
+# Convert KEGG IDs to Entrez Gene IDs
+ URL <- paste("http://rest.kegg.jp/conv",species.KEGG,"ncbi-geneid",sep="/")
+ EntrezID.KEGGGeneID <- read.table(URL,sep="\t",quote="\"",fill=TRUE,comment.char="",stringsAsFactors=FALSE)
+ m <- match(Path.Gene[,1],EntrezID.KEGGGeneID[,2])
+ Path.Gene[,1] <- EntrezID.KEGGGeneID[m,1]
+ Path.Gene[,1] <- sub("^ncbi-geneid:", "", Path.Gene[,1])
+ } else {
+ Path.Gene[,1] <- sub(paste0("^",species.KEGG,":"), "", Path.Gene[,1])
+ }
+ Path.Gene
+}
+
+getKEGGPathwayNames <- function(species.KEGG=NULL, remove.qualifier=FALSE)
+# Read pathways from KEGG website
+# Gordon Smyth
+# Created 7 Jan 2016. Last modified 8 Jan 2016.
+{
+ URL <- "http://rest.kegg.jp/list/pathway"
+ if(!is.null(species.KEGG)) URL <- paste(URL,species.KEGG,sep="/")
+ PathID.PathName <- read.table(URL,sep="\t",quote="\"",fill=TRUE,comment.char="",stringsAsFactors=FALSE)
+ colnames(PathID.PathName) <- c("PathwayID","Description")
+ if(remove.qualifier) PathID.PathName[,2] <- removeExt(PathID.PathName[,2], sep=" - ")
+ PathID.PathName
+}
+
topKEGG <- function(results, sort = NULL, number = 20L, truncate.path=NULL)
# Extract top KEGG pathways from kegga output
# Gordon Smyth and Yifang Hu
-# Created 15 May 2015. Modified 4 August 2015.
+# Created 15 May 2015. Modified 25 Jan 2016.
{
# Check results
if(!is.data.frame(results)) stop("results should be a data.frame.")
@@ -256,11 +309,11 @@ topKEGG <- function(results, sort = NULL, number = 20L, truncate.path=NULL)
# Sort by minimum p-value for specified gene lists
P.col <- 2L+nsets+isort
- if(length(P.col)==1L) {
- o <- order(results[,P.col])
- } else {
- o <- order(do.call("pmin",as.data.frame(results[,P.col,drop=FALSE])))
- }
+ if(nsets==1L)
+ P <- results[,P.col]
+ else
+ P <- do.call(pmin,results[,P.col,drop=FALSE])
+ o <- order(P,results$N,results$Pathway)
tab <- results[o[1L:number],,drop=FALSE]
# Truncate Pathway column for readability
diff --git a/R/lmEffects.R b/R/lmEffects.R
new file mode 100644
index 0000000..08171c6
--- /dev/null
+++ b/R/lmEffects.R
@@ -0,0 +1,139 @@
+.lmEffects <- function(y,design=NULL,contrast=ncol(design),array.weights=NULL,weights=NULL,block=NULL,correlation)
+# Compute matrix of effects from genewise linear models
+# Gordon Smyth
+# Created 11 Apr 2016. Last modified 14 April 2016.
+{
+# Extract components from y
+ y <- getEAWP(y)
+ ngenes <- nrow(y$exprs)
+ n <- ncol(y$exprs)
+
+# Check y
+ if(anyNA(y$exprs)) stop("All y values must be finite and non-NA")
+
+# Check design
+ if(is.null(design)) design <- y$design
+ if(is.null(design)) {
+ stop("design matrix not specified")
+ } else {
+ design <- as.matrix(design)
+ if(mode(design) != "numeric") stop("design must be a numeric matrix")
+ }
+ if(nrow(design) != n) stop("row dimension of design matrix must match column dimension of data")
+ p <- ncol(design)
+
+ if(n <= p) stop("No residual degrees of freedom")
+
+# Check contrast
+ if(is.character(contrast)) {
+ if(length(contrast)>1L) {
+ warning("using only first entry for contrast")
+ contrast <- contrast[1]
+ }
+ contrast <- which(contrast==colnames(design))
+ if(length(contrast)==0L) stop("coef ",contrast," not found")
+ }
+ if(all(contrast == 0)) stop("contrast all zero")
+
+# Reform design matrix so that contrast is last coefficient
+ if(length(contrast) == 1L) {
+ contrast <- as.integer(contrast)
+ if(contrast < p)
+ X <- cbind(design[,-contrast,drop=FALSE],design[,contrast,drop=FALSE])
+ else
+ X <- design
+ } else {
+ if(length(contrast) != p) stop("length of contrast must match column dimension of design")
+ X <- contrastAsCoef(design, contrast, first=FALSE)$design
+ }
+ CoefName <- colnames(X)[p]
+ if(is.null(CoefName)) CoefName <- "Contrast"
+
+# Check array weights
+ if(!is.null(array.weights)) {
+ if(length(array.weights) != n) stop("Length of array.weights doesn't match number of array")
+ AnyNeg <- any(array.weights <= 0)
+ if(anyNA(AnyNeg) || AnyNeg) stop("array.weights must be positive")
+ }
+
+# Check observation weights
+ if(is.null(weights)) weights <- y$weights
+ if(!is.null(weights)) {
+ weights <- as.matrix(weights)
+ dimw <- dim(weights)
+ if(dimw[1] != ngenes || dimw[2] != n) stop("weights must have same dimensions as y")
+ AnyNeg <- any(weights <= 0)
+ if(anyNA(AnyNeg) || AnyNeg) stop("weights must be positive")
+ }
+
+# Reduce to numeric expression matrix
+ y <- y$exprs
+ geneid <- rownames(y)
+ if(is.null(geneid)) geneid <- 1:ngenes
+
+# Divide out array weights
+ if(!is.null(array.weights)) {
+ ws <- sqrt(array.weights)
+ X <- X*ws
+ y <- .matvec(y,ws)
+ array.weights <- NULL
+ }
+
+# Correlation matrix
+ if(!is.null(block)) {
+ if(missing(correlation) || is.null(correlation)) stop("correlation must be set")
+ block <- as.vector(block)
+ if (length(block) != n) stop("Length of block does not match number of arrays")
+ ub <- unique(block)
+ nblocks <- length(ub)
+ Z <- matrix(block,n,nblocks) == matrix(ub,n,nblocks,byrow = TRUE)
+ cormatrix <- Z %*% (correlation * t(Z))
+ diag(cormatrix) <- 1
+ R <- chol(cormatrix)
+
+# Divide out correlations if possible
+ if(is.null(weights)) {
+ y <- t(backsolve(R, t(y), transpose = TRUE))
+ X <- backsolve(R, X, transpose = TRUE)
+ }
+ }
+
+ qrX <- qr(X)
+ if(qrX$rank < p) stop("design must be full column rank")
+
+# Compute effects for contrasts and residuals
+ if(is.null(weights)) {
+ Effects <- t(qr.qty(qrX,t(y)))
+ signc <- sign(qrX$qr[p,p])
+# Remove model effects other than the required contrast
+ if(p>1) Effects <- Effects[,p:n,drop=FALSE]
+# Preserve sign of estimated effect
+ if(signc<0) Effects[,1] <- signc*Effects[,1]
+ } else {
+ Effects <- matrix(0,ngenes,n)
+ signc <- rep.int(0,ngenes)
+ ws <- sqrt(weights)
+ for (g in 1:ngenes) {
+ wX <- X*ws[g,]
+ wy <- y[g,]*ws[g,]
+ if(!is.null(block)) {
+ wy <- backsolve(R,wy,transpose=TRUE)
+ wX <- backsolve(R,wX,transpose=TRUE)
+ }
+ qrX <- qr(wX)
+ signc[g] <- sign(qrX$qr[p,p])
+ Effects[g,] <- qr.qty(qrX,wy)
+ }
+# Remove model effects other than the required contrast
+ if(p>1) Effects <- Effects[,p:n,drop=FALSE]
+# Preserve sign of estimated effect
+ Effects[,1] <- signc*Effects[,1]
+ }
+
+# Dimension names
+ EffectNames <- p:n
+ EffectNames[1] <- CoefName
+ dimnames(Effects) <- list(geneid,c(EffectNames))
+
+ Effects
+}
diff --git a/R/lmfit.R b/R/lmfit.R
index c1f80ab..313345e 100644
--- a/R/lmfit.R
+++ b/R/lmfit.R
@@ -217,40 +217,51 @@ gls.series <- function(M,design=NULL,ndups=2,spacing=1,block=NULL,correlation=NU
# Fit linear model for each gene to a series of microarrays.
# Fit is by generalized least squares allowing for correlation between duplicate spots.
# Gordon Smyth
-# 11 May 2002. Last revised 26 June 2015.
+# 11 May 2002. Last revised 29 Dec 2015.
{
+# Check M
M <- as.matrix(M)
+ ngenes <- nrow(M)
narrays <- ncol(M)
+
+# Check design
if(is.null(design)) design <- matrix(1,narrays,1)
design <- as.matrix(design)
if(nrow(design) != narrays) stop("Number of rows of design matrix does not match number of arrays")
+ nbeta <- ncol(design)
+ coef.names <- colnames(design)
+
+# Check correlation
if(is.null(correlation)) correlation <- duplicateCorrelation(M,design=design,ndups=ndups,spacing=spacing,block=block,weights=weights,...)$consensus.correlation
+ if(abs(correlation) >= 1) stop("correlation is 1 or -1, so the model is degenerate")
+
+# Check weights
if(!is.null(weights)) {
weights[is.na(weights)] <- 0
weights <- asMatrixWeights(weights,dim(M))
M[weights < 1e-15 ] <- NA
weights[weights < 1e-15] <- NA
}
- nbeta <- ncol(design)
- coef.names <- colnames(design)
+
+# Unwrap duplicates (if necessary) and make correlation matrix
if(is.null(block)) {
+# Correlated within-array probes
if(ndups<2) {
- warning("No duplicates: correlation between duplicates set to zero")
+ warning("No duplicates (ndups<2)")
ndups <- 1
correlation <- 0
}
- if(is.null(spacing)) spacing <- 1
cormatrix <- diag(rep(correlation,len=narrays),nrow=narrays,ncol=narrays) %x% array(1,c(ndups,ndups))
+ if(is.null(spacing)) spacing <- 1
M <- unwrapdups(M,ndups=ndups,spacing=spacing)
if(!is.null(weights)) weights <- unwrapdups(weights,ndups=ndups,spacing=spacing)
design <- design %x% rep(1,ndups)
colnames(design) <- coef.names
+ ngenes <- nrow(M)
+ narrays <- ncol(M)
} else {
- if(ndups>1) {
- stop("Cannot specify ndups>2 and non-null block argument")
- } else {
- ndups <- spacing <- 1
- }
+# Correlated samples
+ ndups <- spacing <- 1
block <- as.vector(block)
if(length(block)!=narrays) stop("Length of block does not match number of arrays")
ub <- unique(block)
@@ -259,9 +270,10 @@ gls.series <- function(M,design=NULL,ndups=2,spacing=1,block=NULL,correlation=NU
cormatrix <- Z%*%(correlation*t(Z))
}
diag(cormatrix) <- 1
- ngenes <- nrow(M)
+
stdev.unscaled <- matrix(NA,ngenes,nbeta,dimnames=list(rownames(M),coef.names))
+# If weights and missing values are absent, use a fast computation
NoProbeWts <- all(is.finite(M)) && (is.null(weights) || !is.null(attr(weights,"arrayweights")))
if(NoProbeWts) {
V <- cormatrix
@@ -299,6 +311,7 @@ gls.series <- function(M,design=NULL,ndups=2,spacing=1,block=NULL,correlation=NU
return(fit)
}
+# Weights or missing values are present, to have to iterate over probes
beta <- stdev.unscaled
sigma <- rep(NA,ngenes)
df.residual <- rep(0,ngenes)
@@ -370,10 +383,10 @@ nonEstimable <- function(x)
}
}
-fitted.MArrayLM <- function(object,design=object$design,...)
+fitted.MArrayLM <- function(object,...)
# Fitted values from MArray linear model fit
# Gordon Smyth
-# 29 November 2005
+# 29 November 2005. Last modified 4 May 2016.
{
object$coefficients %*% t(object$design)
}
diff --git a/R/loessFit.R b/R/loessFit.R
index 3d18152..4691704 100644
--- a/R/loessFit.R
+++ b/R/loessFit.R
@@ -81,7 +81,7 @@ loessFit <- function(y, x, weights=NULL, span=0.3, iterations=4L, min.weight=1e-
},
"locfit" = {
# Check for locfit package
- if(!requireNamespace("locfit",quietly=TRUE)) stop("locfit required but is not available")
+ if(!requireNamespace("locfit",quietly=TRUE)) stop("locfit required but is not installed (or can't be loaded)")
# Weighted lowess with robustifying iterations
biweights <- rep(1,nobs)
for (i in 1:iterations) {
diff --git a/R/neqc.R b/R/neqc.R
index fb7a4e6..ccd5069 100644
--- a/R/neqc.R
+++ b/R/neqc.R
@@ -95,7 +95,7 @@ normexp.fit.control <- function(x,status=NULL,negctrl="negative",regular="regula
nec <- function(x,status=NULL,negctrl="negative",regular="regular",offset=16,robust=FALSE,detection.p="Detection")
# Normexp background correction aided by negative controls.
# Wei Shi and Gordon Smyth
-# Created 27 September 2010. Last modified 21 Oct 2014.
+# Created 27 September 2010. Last modified 15 Nov 2015.
{
if(is(x, "EListRaw")) {
if(!is.null(x$Eb)) {
@@ -107,7 +107,7 @@ nec <- function(x,status=NULL,negctrl="negative",regular="regular",offset=16,rob
normexp.par <- normexp.fit.control(x,status,negctrl,regular,robust)
} else {
normexp.par <- normexp.fit.detection.p(x,detection.p)
- message("Inferred negative control probe intensities were used in background correction.")
+ message("Note: inferring mean and variance of negative control probe intensities from the detection p-values.")
}
for(i in 1:ncol(x)) x$E[,i] <- normexp.signal(normexp.par[i,], x$E[,i])
x$E <- x$E + offset
diff --git a/R/norm.R b/R/norm.R
index 5c2f660..dbdf393 100755
--- a/R/norm.R
+++ b/R/norm.R
@@ -333,7 +333,12 @@ plotPrintorder <- function(object,layout,start="topleft",slide=1,method="loess",
normalizeBetweenArrays <- function(object, method=NULL, targets=NULL, cyclic.method="fast", ...) {
# Normalize between arrays
# Gordon Smyth
-# 12 Apri 2003. Last revised 21 July 2013.
+# 12 Apri 2003. Last revised 13 June 2016.
+
+ if(is.data.frame(object)) {
+ object <- as.matrix(object)
+ if(mode(object) != "numeric") stop("'object' is a data.frame and not all columns are numeric")
+ }
# Default method
if(is.null(method)) {
diff --git a/R/plotrldf.R b/R/plotrldf.R
index 2a7d236..c8fee14 100644
--- a/R/plotrldf.R
+++ b/R/plotrldf.R
@@ -1,52 +1,113 @@
# Plot regularized linear discriminant functions
-plotRLDF <- function(y,design=NULL,z=NULL,labels.y=NULL,labels.z=NULL,col.y=1,col.z=1,df.prior=5,show.dimensions=c(1,2),main=NULL,nprobes=500,...)
+plotRLDF <- function(y,design=NULL,z=NULL,labels.y=NULL,labels.z=NULL,col.y="black",col.z="black",show.dimensions=c(1,2),ndim=max(show.dimensions),nprobes=100,plot=TRUE,var.prior=NULL,df.prior=NULL,trend=FALSE,robust=FALSE,...)
# Regularized linear discriminant function
# Di Wu and Gordon Smyth
-# 29 June 2009. Last revised 30 June 2009.
+# 29 June 2009. Last revised 26 Jan 2016.
{
+# Check y
y <- as.matrix(y)
g <- nrow(y)
n <- ncol(y)
+
+# Check labels.y
+ if(is.null(labels.y))
+ labels.y <- colnames(y)
+ else {
+ if(length(labels.y) != n) stop("length(labels.y) doesn't agree with ncol(y)")
+ labels.y <- as.character(labels.y)
+ }
+ if(is.null(labels.y)) labels.y <- as.character(1:n)
+
+# Check z and labels.z
+ if(is.null(z)) {
+ labels.z <- character(0)
+ } else {
+ z <- as.matrix(z)
+ if(nrow(z) != g) stop("nrow(z) disagrees with nrow(y)")
+ if(is.null(labels.z)) labels.z <- colnames(z)
+ if(is.null(labels.z)) labels.z <- letters[1:ncol(z)]
+ if(length(labels.z) != ncol(z)) stop("length(labels.z) doesn't agree with ncol(z)")
+ labels.z <- as.character(labels.z)
+ if(!all(rownames(z)==rownames(y))) warning("y and z have different rownames - they are assumed to correspond to same probes")
+ }
+
+# Check design
if(is.null(design)) {
if(is.null(labels.y)) stop("groups not specified")
+ if(!anyDuplicated(labels.y)) stop("design not specified and all labels.y are different")
f <- as.factor(labels.y)
design <- model.matrix(~f)
+ } else {
+ design <- as.matrix(design)
+ if(nrow(design) != n) stop("nrow(design) doesn't match ncol(y)")
}
- design <- as.matrix(design)
- if(nrow(design) != n) stop("nrow(design) doesn't match ncol(y)")
- if(is.null(labels.y)) labels.y <- 1:n
+
+# Check show.dimensions
+ show.dimensions <- as.integer(show.dimensions)[1:2]
+ if(show.dimensions[1]==show.dimensions[2]) stop("show.dimensions must specify two different columns")
+ if(any(show.dimensions<1) || any(show.dimensions)>n) stop("show.dimensions must be a column number of y")
+
+# Check ndim
+ if(any(show.dimensions)>ndim) ndim <- max(show.dimensions)
+
+# Check nprobes
+ if(nprobes<1) stop("'nprobes' must be at least 1")
# Project onto between and within spaces
# Discard first column as intercept
qrd <- qr(design)
p <- qrd$rank
- if(p==n) stop("No residual degrees of freedom")
+ df.residual <- n-p
+ if(df.residual==0) stop("No residual degrees of freedom")
U <- qr.qty(qrd, t(y))
UB <- U[2:p,,drop=FALSE]
UW <- U[(p+1):n,,drop=FALSE]
+ s2 <- colMeans(UW*UW)
-# Prior variance
- s <- colMeans(UW*UW)
- s0 <- median(s)
+# Prior variance and prior df
+ if(is.null(var.prior) || is.null(df.prior)) {
+ if(trend) covariate <- rowMeans(y) else covariate <- NULL
+ sv <- squeezeVar(s2, df=df.residual, covariate=covariate, robust=robust)
+ var.prior <- sv$var.prior
+ df.prior <- sv$df.prior
+ } else {
+ if(!any(length(var.prior)==c(1,g))) stop("var.prior wrong length")
+ if(!any(length(df.prior)==c(1,g))) stop("df.prior wrong length")
+ }
+ df.prior <- pmin(df.prior, (g-1)*df.residual)
+ df.prior <- pmax(df.prior, 1)
# Select probes by moderated F
if(g>nprobes) {
- modF <- colMeans(UB*UB)/(s+df.prior*s0)
+ modF <- colMeans(UB*UB)/(s2+df.prior*var.prior)
o <- order(modF,decreasing=TRUE)
top <- o[1:nprobes]
y <- y[top,,drop=FALSE]
if(!is.null(z)) z <- z[top,,drop=FALSE]
UB <- UB[,top,drop=FALSE]
UW <- UW[,top,drop=FALSE]
+ if(length(df.prior)>1) df.prior <- df.prior[top]
+ if(length(var.prior)>1) var.prior <- var.prior[top]
g <- nprobes
+ } else {
+ top <- 1:nprobes
}
# Within group SS
W <- crossprod(UW)
-# Regularized within group SS
- Wreg <- W+diag(df.prior*s0,g,g)
+# Regularize the within-group covariance matrix
+ Wreg <- W
+ diag(Wreg) <- diag(Wreg) + df.prior*var.prior
+ df.total <- df.prior + df.residual
+ if(length(df.total)>1) {
+ df.total <- sqrt(df.total)
+ Wreg <- Wreg / df.total
+ Wreg <- (t(Wreg) / df.total)
+ } else {
+ Wreg <- Wreg / df.total
+ }
# Ratio of between to within SS
WintoB <- backsolve(chol(Wreg),t(UB),transpose=TRUE)
@@ -54,26 +115,38 @@ plotRLDF <- function(y,design=NULL,z=NULL,labels.y=NULL,labels.z=NULL,col.y=1,co
# Linear discriminant gene weights
d1 <- show.dimensions[1]
d2 <- show.dimensions[2]
- metagenes <- svd(WintoB,nu=max(d1,d2),nv=0)$u
+ SVD <- svd(WintoB,nu=ndim,nv=0)
+ metagenes <- SVD$u
- d1.y <- t(y)%*%metagenes[,d1]
- d2.y <- t(y)%*%metagenes[,d2]
- if(!is.null(z)) {
- z <- as.matrix(z)
- d1.z <- t(z)%*%metagenes[,d1]
- d2.z <- t(z)%*%metagenes[,d2]
- if(is.null(labels.z)) labels.z <- letters[1:ncol(z)]
- } else {
+ rank <- min(dim(WintoB))
+# if(ndim > rank) message("Note: only ",rank," of the discriminant functions have any predictive power")
+
+# LDF for training set
+ d.y <- t(y) %*% metagenes
+ d1.y <- d.y[,d1]
+ d2.y <- d.y[,d2]
+
+# LDF for classified set
+ if(is.null(z)) {
d1.z <- d2.z <- numeric(0)
- labels.z <- character(0)
+ } else {
+ d.z <- t(z) %*% metagenes
+ d1.z <- d.z[,d1]
+ d2.z <- d.z[,d2]
}
- lab.y <- as.character(labels.y)
- lab.z <- as.character(labels.z)
- plot(c(d1.y,d1.z),c(d2.y,d2.z),type="n",main=main,
- xlab=paste("Discriminant Function",d1),
- ylab=paste("Discriminant Function",d2))
- text(d1.y,d2.y,labels=lab.y,col=col.y,...)
- if(!is.null(z)) text(d1.z,d2.z,labels=lab.z,col=col.z,...)
- invisible(list(training=cbind(d1=d1.y,d2=d2.y),predicting=cbind(d1=d1.z,d2=d2.z),metagenes=metagenes,top=top))
+
+# Make plot
+ if(plot) {
+ plot(c(d1.y,d1.z),c(d2.y,d2.z),type="n", xlab=paste("Discriminant Function",d1), ylab=paste("Discriminant Function",d2))
+ text(d1.y,d2.y,labels=labels.y,col=col.y,...)
+ if(!is.null(z)) text(d1.z,d2.z,labels=labels.z,col=col.z,...)
+ }
+
+# Output
+ out <- list(training=d.y,top=top,metagenes=metagenes,singular.values=SVD$d,rank=rank)
+ if(!is.null(z)) out$predicting <- d.z
+ out$var.prior <- var.prior
+ out$df.prior <- df.prior
+ invisible(out)
}
diff --git a/R/propTrueNull.R b/R/propTrueNull.R
index 0c5112b..7c8b133 100644
--- a/R/propTrueNull.R
+++ b/R/propTrueNull.R
@@ -47,10 +47,12 @@ propTrueNull <- function(p, method="lfdr", nbins=20, ...)
# http://www.public.iastate.edu/~dnett/microarray/multtest.txt
# Accessed March 2012
-# Created 2 Dec 2012. Last revised 2 Dec 2012.
+# Created 2 Dec 2012. Last revised 9 May 2016.
{
bin <- c(-0.1, (1:nbins)/nbins)
- bin.counts <- tabulate(cut(p,bin))
+ bin.counts <- rep_len(0L, length.out=nbins)
+ tab <- tabulate(cut(p,bin))
+ bin.counts[1:length(tab)] <- tab
tail.means <- rev(cumsum(rev(bin.counts))/(1:nbins))
index <- which(tail.means >= bin.counts)[1]
tail.means[index]/tail.means[1]
diff --git a/R/read-maimages.R b/R/read-maimages.R
index 5fb725e..aa3ace0 100644
--- a/R/read-maimages.R
+++ b/R/read-maimages.R
@@ -176,6 +176,7 @@ read.maimages <- function(files=NULL,source="generic",path=NULL,ext=NULL,names=N
for (a in cnames) RG[[a]] <- Y
if(!is.null(wt.fun)) RG$weights <- Y
if(is.data.frame(targets)) {
+ rownames(targets) <- names
RG$targets <- targets
} else {
RG$targets <- data.frame(FileName=files,row.names=names,stringsAsFactors=FALSE)
diff --git a/R/read.R b/R/read.R
index cf31e7e..b047e8f 100755
--- a/R/read.R
+++ b/R/read.R
@@ -253,17 +253,19 @@ controlStatus <- function(types, genes, spottypecol="SpotType", regexpcol, verbo
status
}
-removeExt <- function(x) {
+removeExt <- function(x, sep=".")
# Remove any common extension from a vector of file names
# Gordon Smyth
-# 19 July 2002. Last modified 19 Jan 2005.
-
+# 19 July 2002. Last modified 8 Jan 2016.
+{
x <- as.character(x)
n <- length(x)
- if(length(grep("\\.",x)) < n) return(x)
- ext <- sub("(.*)\\.(.*)$","\\2",x)
+ if(length(grep(sep,x)) < n) return(x)
+ sep <- protectMetachar(sep)
+ RegExpr <- paste0("(.*)", sep, "(.*)$")
+ ext <- sub(RegExpr, "\\2", x)
if(all(ext[1] == ext))
- return(sub("(.*)\\.(.*)$","\\1",x))
+ return(sub(RegExpr,"\\1",x))
else
return(x)
}
diff --git a/R/read.idat.R b/R/read.idat.R
index 8e040b5..4c5f607 100644
--- a/R/read.idat.R
+++ b/R/read.idat.R
@@ -1,43 +1,88 @@
-read.idat <- function(idatfiles, bgxfile, dateinfo=FALSE, tolerance=0)
-# Read idat data from gene expression BeadChips
-# Matt Ritchie
-# Created 30 September 2013. Modified on 14 January 2015 and 26 June 2015.
+read.idat <- function(idatfiles, bgxfile, dateinfo=FALSE, annotation="Symbol", tolerance=0L, verbose=TRUE)
+# Read GenomeStudio IDAT files for Illumina gene expression BeadChips
+# Matt Ritchie and Gordon Smyth
+# Created 30 September 2013. Last modified 14 June 2016.
{
- if(!requireNamespace("illuminaio",quietly=TRUE)) stop("illuminaio package required but is not available")
- cat("Reading manifest file", bgxfile, "... ")
- bgx <- illuminaio::readBGX(bgxfile)
- cat("Done\n")
- nregprobes <-nrow(bgx$probes)
- nctrlprobes <-nrow(bgx$control)
- nprobes <- nregprobes+nctrlprobes
- nsamples <- length(idatfiles)
- elist <- new("EListRaw")
- elist$source <- "illumina"
- elist$genes <- rbind(bgx$probes[,c("Probe_Id","Array_Address_Id")], bgx$controls[,c("Probe_Id","Array_Address_Id")])
- elist$genes$Status <- "regular"
- elist$genes$Status[(nregprobes+1):nprobes] <- bgx$controls[,"Reporter_Group_Name"]
- elist$targets <- data.frame("IDATfile"=idatfiles, "DecodeInfo"=rep(NA, nsamples), "ScanInfo"=rep(NA, nsamples))
- tmp <- matrix(NA, nprobes, nsamples)
- colnames(tmp) <- idatfiles
- rownames(tmp) <- elist$genes[,"Array_Address_Id"]
- elist$E <- elist$other$STDEV <- elist$other$NumBeads <- tmp
- rm(tmp)
- for(i in 1:nsamples) {
- cat("\t", idatfiles[i], "... ")
- tmp <- illuminaio::readIDAT(idatfiles[i])
- cat("Done\n")
- ind <- match(elist$genes[,"Array_Address_Id"], tmp$Quants[,'IllumicodeBinData'])
- if(sum(is.na(ind))>tolerance)
- stop("Can't match all ids in manifest with those in idat file ", idatfiles[i], "\n", sum(is.na(ind)),
- " missing - please check that you have the right files, or consider setting \'tolerance\'=", sum(is.na(ind)))
- elist$E[!is.na(ind),i] <- round(tmp$Quants[ind[!is.na(ind)],'MeanBinData'],1) # intensity data
- elist$other$STDEV[!is.na(ind),i] <- tmp$Quants[ind[!is.na(ind)],'DevBinData'] # Bead STDEV
- elist$other$NumBeads[!is.na(ind),i] <- tmp$Quants[ind[!is.na(ind)],'NumGoodBeadsBinData'] # NumBeads
- if(dateinfo) {
- elist$targets[i,"DecodeInfo"] = paste(tmp$RunInfo[1,], sep="", collapse=" ")
- elist$targets[i,"ScanInfo"] = paste(tmp$RunInfo[2,], sep="", collapse=" ")
- }
- }
- cat("Finished reading data.\n")
- return(elist)
+# Need illuminaio package
+ OK <- requireNamespace("illuminaio",quietly=TRUE)
+ if(!OK) stop("illuminaio package required but is not installed (or can't be loaded)")
+
+# Check idatfiles
+ idatafiles <- as.character(idatfiles)
+ nsamples <- length(idatfiles)
+
+# Initialize EListRaw object
+ elist <- new("EListRaw")
+ elist$source <- "illumina"
+ elist$targets <- data.frame("IDATfile"=idatfiles,stringsAsFactors=FALSE)
+ if(dateinfo) elist$targets$DecodeInfo <- elist$targets$ScanInfo <- rep_len("", nsamples)
+
+# Read bead manifest file
+ if(verbose) cat("Reading manifest file", bgxfile, "... ")
+ bgx <- illuminaio::readBGX(bgxfile)
+ if(verbose) cat("Done\n")
+ nregprobes <- nrow(bgx$probes)
+ nctrlprobes <- nrow(bgx$control)
+ nprobes <- nregprobes+nctrlprobes
+
+# Assemble gene annotation
+ elist$genes <- rbind(bgx$probes[,c("Probe_Id","Array_Address_Id")], bgx$controls[,c("Probe_Id","Array_Address_Id")])
+
+# Set probe control status
+ elist$genes$Status <- "regular"
+ elist$genes$Status[(nregprobes+1):nprobes] <- bgx$controls[,"Reporter_Group_Name"]
+
+# Add optional annotation columns
+ if(!is.null(annotation) && !is.na(annotation)) {
+ annotation <- as.character(annotation)
+ annotation <- intersect(names(bgx$probes),annotation)
+ }
+ if(length(annotation)) {
+ ac <- annotation %in% names(bgx$controls)
+ for (i in 1:length(annotation)) {
+ elist$genes[[annotation[i]]] <- NA_character_
+ elist$genes[[annotation[i]]][1:nregprobes] <- bgx$probes[[annotation[i]]]
+ if(ac[i]) elist$genes[[annotation[i]]][(nregprobes+1L):nprobes] <- bgx$controls[[annotation[i]]]
+ }
+ }
+
+# Initalize expression matrices
+ elist$E <- matrix(0, nprobes, nsamples)
+ elist$E[] <- NA
+ colnames(elist$E) <- removeExt(idatfiles)
+ rownames(elist$E) <- elist$genes[,"Array_Address_Id"]
+ elist$other$STDEV <- elist$other$NumBeads <- elist$E
+
+# Read IDAT files
+ for(j in 1:nsamples) {
+ if(verbose) cat("\t", idatfiles[j], "... ")
+ tmp <- illuminaio::readIDAT(idatfiles[j])
+ if(verbose) cat("Done\n")
+ ind <- match(elist$genes$Array_Address_Id, tmp$Quants$IllumicodeBinData)
+
+# Check for whether values are available for all probes
+ if(anyNA(ind)) {
+ nna <- sum(is.na(ind))
+ if(nna > tolerance)
+ stop("Can't match all ids in manifest with those in idat file ", idatfiles[i], "\n", sum(is.na(ind)),
+ " missing - please check that you have the right files, or consider setting \'tolerance\'=", sum(is.na(ind)))
+ i <- which(!is.na(ind))
+ ind <- ind[i]
+ elist$E[i,j] <- tmp$Quants$MeanBinData[ind]
+ elist$other$STDEV[i,j] <- tmp$Quants$DevBinData[ind]
+ elist$other$NumBeads[i,j] <- tmp$Quants$NumGoodBeadsBinData[ind]
+ } else {
+ elist$E[,j] <- tmp$Quants$MeanBinData[ind]
+ elist$other$STDEV[,j] <- tmp$Quants$DevBinData[ind]
+ elist$other$NumBeads[,j] <- tmp$Quants$NumGoodBeadsBinData[ind]
+ }
+
+ if(dateinfo) {
+ elist$targets$DecodeInfo[j] = paste(tmp$RunInfo[1,], collapse=" ")
+ elist$targets$ScanInfo[j] = paste(tmp$RunInfo[2,], collapse=" ")
+ }
+ }
+
+ if(verbose) cat("Finished reading data.\n")
+ return(elist)
}
diff --git a/R/sepchannel.R b/R/sepchannel.R
index 1991562..cd0305c 100755
--- a/R/sepchannel.R
+++ b/R/sepchannel.R
@@ -87,7 +87,7 @@ intraspotCorrelation <- function(object,design,trim=0.15)
designA <- (Ident %x% matrix(c(0.5,0.5),1,2)) %*% design
X <- rbind(designM, designA)
Z <- diag(2) %x% rep(1,narrays)
- if(!requireNamespace("statmod",quietly=TRUE)) stop("statmod package required but is not available")
+ if(!requireNamespace("statmod",quietly=TRUE)) stop("statmod package required but is not installed")
arho <- rep(NA,ngenes)
degfre <- matrix(0,ngenes,2,dimnames=list(rownames(M),c("df.M","df.A")))
for (i in 1:ngenes) {
diff --git a/R/squeezeVar.R b/R/squeezeVar.R
index 387d962..6cf6494 100644
--- a/R/squeezeVar.R
+++ b/R/squeezeVar.R
@@ -3,190 +3,59 @@
squeezeVar <- function(var, df, covariate=NULL, robust=FALSE, winsor.tail.p=c(0.05,0.1))
# Empirical Bayes posterior variances
# Gordon Smyth
-# 2 March 2004. Last modified 2 Dec 2013.
+# Created 2 March 2004. Last modified 5 May 2016.
{
n <- length(var)
+
+# Degenerate special cases
if(n == 0) stop("var is empty")
if(n == 1) return(list(var.post=var,var.prior=var,df.prior=0))
- if(length(df)==1) {
- df <- rep.int(df,n)
- } else {
- if(length(df) != n) stop("lengths differ")
- }
-
-# Estimate prior var and df
- if(robust)
- fit <- fitFDistRobustly(var, df1=df, covariate=covariate, winsor.tail.p=winsor.tail.p)
- else
- fit <- fitFDist(var, df1=df, covariate=covariate)
-
-# Prior var will be vector if robust=TRUE, otherwise scalar
- var.prior <- fit$scale
-
-# Prior df will be vector if covariate is non-NULL, otherwise scalar
- df.prior <- fit$df2.shrunk
- if(is.null(df.prior)) df.prior <- fit$df2
-# Check estimated prior df
- if(is.null(df.prior) || any(is.na(df.prior))) stop("Could not estimate prior df")
+# When df==0, guard against missing or infinite values in var
+ if(length(df)>1) var[df==0] <- 0
-# Squeeze the posterior variances
- df.total <- df + df.prior
- var[df==0] <- 0 # guard against missing or infinite values
- Infdf <- df.prior==Inf
- if(any(Infdf)) {
- var.post <- rep(var.prior,length.out=n)
- i <- which(!Infdf)
- if(length(i)) {
- if(is.null(covariate))
- s02 <- var.prior
- else
- s02 <- var.prior[i]
- var.post[i] <- (df[i]*var[i] + df.prior[i]*s02) / df.total[i]
- }
+# Estimate hyperparameters
+ if(robust) {
+ fit <- fitFDistRobustly(var, df1=df, covariate=covariate, winsor.tail.p=winsor.tail.p)
+ df.prior <- fit$df2.shrunk
} else {
- var.post <- (df*var + df.prior*var.prior) / df.total
+ fit <- fitFDist(var, df1=df, covariate=covariate)
+ df.prior <- fit$df2
}
+ if(anyNA(df.prior)) stop("Could not estimate prior df")
- list(df.prior=df.prior,var.prior=var.prior,var.post=var.post)
+# Posterior variances
+ var.post <- .squeezeVar(var=var, df=df, var.prior=fit$scale, df.prior=df.prior)
+
+ list(df.prior=df.prior,var.prior=fit$scale,var.post=var.post)
}
-fitFDist <- function(x,df1,covariate=NULL)
-# Moment estimation of the parameters of a scaled F-distribution
-# The first degrees of freedom are given
-# Gordon Smyth and Belinda Phipson
-# 8 Sept 2002. Last revised 27 Oct 2012.
+.squeezeVar <- function(var, df, var.prior, df.prior)
+# Squeeze posterior variances given hyperparameters
+# NAs not allowed in df.prior
+# Gordon Smyth
+# Created 5 May 2016
{
-# Check covariate
- if(!is.null(covariate)) {
- if(length(covariate) != length(x)) stop("covariate and x must be of same length")
- if(any(is.na(covariate))) stop("NA covariate values not allowed")
- isfin <- is.finite(covariate)
- if(!all(isfin)) {
- if(!any(isfin))
- covariate <- sign(covariate)
- else {
- r <- range(covariate[isfin])
- covariate[covariate == -Inf] <- r[1]-1
- covariate[covariate == Inf] <- r[2]+1
- }
- }
- splinedf <- min(4,length(unique(covariate)))
- if(splinedf < 2) covariate <- NULL
- }
-# Remove missing or infinite values and zero degrees of freedom
- ok <- is.finite(x) & is.finite(df1) & (x > -1e-15) & (df1 > 1e-15)
- notallok <- !all(ok)
- if(notallok) {
- x <- x[ok]
- df1 <- df1[ok]
- if(!is.null(covariate)) {
- covariate2 <- covariate[!ok]
- covariate <- covariate[ok]
- }
- }
- n <- length(x)
- if(n==0) return(list(scale=NA,df2=NA))
-
-# Avoid exactly zero values
- x <- pmax(x,0)
- m <- median(x)
- if(m==0) {
- warning("More than half of residual variances are exactly zero: eBayes unreliable")
- m <- 1
- } else {
- if(any(x==0)) warning("Zero sample variances detected, have been offset",call.=FALSE)
- }
- x <- pmax(x, 1e-5 * m)
+ n <- length(var)
+ isfin <- is.finite(df.prior)
+ if(all(isfin)) return( (df*var + df.prior*var.prior) / (df+df.prior) )
-# Better to work on with log(F)
- z <- log(x)
- e <- z-digamma(df1/2)+log(df1/2)
+# From here, at least some df.prior are infinite
- if(is.null(covariate)) {
- emean <- mean(e)
- evar <- sum((e-emean)^2)/(n-1)
+# For infinite df.prior, return var.prior
+ if(length(var.prior) == n) {
+ var.post <- var.prior
} else {
- if(!requireNamespace("splines",quietly=TRUE)) stop("splines package required but is not available")
- design <- try(splines::ns(covariate,df=splinedf,intercept=TRUE),silent=TRUE)
- if(is(design,"try-error")) stop("Problem with covariate")
- fit <- lm.fit(design,e)
- if(notallok) {
- design2 <- predict(design,newx=covariate2)
- emean <- rep.int(0,n+length(covariate2))
- emean[ok] <- fit$fitted
- emean[!ok] <- design2 %*% fit$coefficients
- } else {
- emean <- fit$fitted
- }
- evar <- mean(fit$residuals[-(1:fit$rank)]^2)
+ var.post <- rep_len(var.prior, length.out=n)
}
- evar <- evar - mean(trigamma(df1/2))
- if(evar > 0) {
- df2 <- 2*trigammaInverse(evar)
- s20 <- exp(emean+digamma(df2/2)-log(df2/2))
- } else {
- df2 <- Inf
- s20 <- exp(emean)
- }
- list(scale=s20,df2=df2)
-}
-
-trigammaInverse <- function(x) {
-# Solve trigamma(y) = x for y
-# Gordon Smyth
-# 8 Sept 2002. Last revised 12 March 2004.
-
-# Non-numeric or zero length input
- if(!is.numeric(x)) stop("Non-numeric argument to mathematical function")
- if(length(x)==0) return(numeric(0))
-# Treat out-of-range values as special cases
- omit <- is.na(x)
- if(any(omit)) {
- y <- x
- if(any(!omit)) y[!omit] <- Recall(x[!omit])
- return(y)
- }
- omit <- (x < 0)
- if(any(omit)) {
- y <- x
- y[omit] <- NaN
- warning("NaNs produced")
- if(any(!omit)) y[!omit] <- Recall(x[!omit])
- return(y)
- }
- omit <- (x > 1e7)
- if(any(omit)) {
- y <- x
- y[omit] <- 1/sqrt(x[omit])
- if(any(!omit)) y[!omit] <- Recall(x[!omit])
- return(y)
- }
- omit <- (x < 1e-6)
- if(any(omit)) {
- y <- x
- y[omit] <- 1/x[omit]
- if(any(!omit)) y[!omit] <- Recall(x[!omit])
- return(y)
+# Maybe some df.prior are finite
+ if(any(isfin)) {
+ i <- which(isfin)
+ if(length(df)>1) df <- df[i]
+ df.prior <- df.prior[i]
+ var.post[i] <- (df*var[i] + df.prior*var.post[i]) / (df+df.prior)
}
-# Newton's method
-# 1/trigamma(y) is convex, nearly linear and strictly > y-0.5,
-# so iteration to solve 1/x = 1/trigamma is monotonically convergent
- y <- 0.5+1/x
- iter <- 0
- repeat {
- iter <- iter+1
- tri <- trigamma(y)
- dif <- tri*(1-tri/x)/psigamma(y,deriv=2)
- y <- y+dif
- if(max(-dif/y) < 1e-8) break
- if(iter > 50) {
- warning("Iteration limit exceeded")
- break
- }
- }
- y
+ var.post
}
-
diff --git a/R/tricubeMovingAverage.R b/R/tricubeMovingAverage.R
index 46c51b6..3b4b660 100644
--- a/R/tricubeMovingAverage.R
+++ b/R/tricubeMovingAverage.R
@@ -1,7 +1,7 @@
tricubeMovingAverage <- function(x,span=0.5,power=3)
-# Moving average filter for a time series with tricube weights
-# Gordon Smyth
-# Created 5 Feb 2014. Last modified 18 August 2015.
+# Moving average filter with tricube weights for a time series
+# Gordon Smyth and Yifang Hu
+# Created 5 Feb 2014. Last modified 6 Nov 2015.
{
# Check span
if(length(span)>1) {
@@ -24,11 +24,16 @@ tricubeMovingAverage <- function(x,span=0.5,power=3)
# Round width of smoothing window to nearest odd number
hwidth <- as.integer(width %/% 2L)
- if(hwidth <= 0L) return(x)
width <- 2L * hwidth + 1L
# Make sure window width can't be greater than n
- if(width>n) width <- n-1L
+ if(width>n) {
+ width <- width-2L
+ hwidth <- hwidth-1L
+ }
+
+# If width is 1 or smaller, return original series
+ if(hwidth <= 0L) return(x)
# Tricube weights with all positive values
u <- seq(from=-1,to=1,length=width) * width / (width+1)
diff --git a/R/voom.R b/R/voom.R
index fd2674b..eaa7cc0 100644
--- a/R/voom.R
+++ b/R/voom.R
@@ -1,8 +1,8 @@
-voom <- function(counts,design=NULL,lib.size=NULL,normalize.method="none",plot=FALSE,span=0.5,...)
-# Linear modelling of count data with mean-variance modelling at the observational level.
-# Creates an EList object for entry to lmFit() etc in the limma pipeline.
-# Gordon Smyth and Charity Law
-# Created 22 June 2011. Last modified 5 June 2013.
+voom <- function(counts,design=NULL,lib.size=NULL,normalize.method="none",span=0.5,plot=FALSE,save.plot=FALSE,...)
+# Linear modelling of count data with mean-variance modelling at the observational level.
+# Creates an EList object for entry to lmFit() etc in the limma pipeline.
+# Gordon Smyth and Charity Law
+# Created 22 June 2011. Last modified 27 Jan 2016.
{
out <- list()
@@ -88,6 +88,11 @@ voom <- function(counts,design=NULL,lib.size=NULL,normalize.method="none",plot=F
out$targets <- data.frame(lib.size=lib.size)
else
out$targets$lib.size <- lib.size
+ if(save.plot) {
+ out$voom.xy <- list(x=sx,y=sy,xlab="log2( count size + 0.5 )",ylab="Sqrt( standard deviation )")
+ out$voom.line <- l
+ }
+
new("EList",out)
}
diff --git a/build/vignette.rds b/build/vignette.rds
index e4054da..543d17a 100644
Binary files a/build/vignette.rds and b/build/vignette.rds differ
diff --git a/inst/NEWS.Rd b/inst/NEWS.Rd
index 01e131c..2f162b9 100644
--- a/inst/NEWS.Rd
+++ b/inst/NEWS.Rd
@@ -2,6 +2,229 @@
\title{limma News}
\encoding{UTF-8}
+
+\section{Version 3.28.0}{\itemize{
+
+\item
+ Improved capabilities and performance for fry(). fry() has two
+ new arguments, 'geneid' and 'standardize'. The index argument of
+ fry() can now be a list of data.frames containing identifiers and
+ weights for each set. The options introduced by the standardize
+ argument allows fry() to be more statistically powerful when genes have unequal
+ variances. fry() now accepts any arguments that would be suitable
+ for lmFit() or eBayes().
+ The 'sort' argument of fry() is now the same as for mroast().
+
+\item
+ roast(), mroast(), fry() and camera() now require the 'design'
+ matrix to be set. Previously the design matrix defaulted to a
+ single intercept column.
+
+\item
+ Two changes to barcodeplot(): new argument 'alpha' to set
+ semitransparency of positional bars in some circumstances;
+ the default value for 'quantiles' is increased slightly.
+
+\item
+ kegga() has several new arguments and now supports any species
+ supported by KEGG. kegga() can now accept annotation as a
+ data.frame, meaning that it can work from user-supplied annotation.
+
+\item
+ New functions getGeneKEGGLinks() and getKEGGPathwayNames() to get
+ pathway annotation from the rest.kegg website.
+
+\item
+ topKEGG() now breaks tied p-values by number of genes in pathway and
+ by name of pathway.
+
+\item
+ goana() now supports species="Pt" (chimpanzee).
+
+\item
+ plotRLDF() now produces more complete output and the output is
+ more completely documented. It can now be used as a complete LDF
+ analysis for classification based on training data. The argument
+ 'main' has been removed as an explicit argument because it and
+ other plotting arguments are better passed using the ... facility.
+ New arguments 'ndim' and 'plot' have been added. The first allows
+ all possible discriminant functions to be completed even if only
+ two are plotted. The second allows dicriminant functions to be
+ computed without a plot.
+
+ plotRLDF() also now uses squeezeVar() to estimate how by how much the
+ within-group covariance matrix is moderated. It has new arguments
+ 'trend' and 'robust' as options to the empirical Bayes estimation
+ step. It also has new argument 'var.prior' to allow the prior
+ variances to be optionally supplied. Argument 'df.prior' can now be
+ a vector of genewise values.
+
+\item
+ new argument 'save.plot' for voom().
+
+\item
+ diffSplice() now returns genewise residual variances.
+
+\item
+ removeExt() has new 'sep' argument.
+
+\item
+ Slightly improved algorithm for producing the df2.shrunk values
+ returned by fitFDistRobustly(). fitFDistRobustly() now returns
+ tail.p.value, prob.outlier and df2.outlier as well as previously
+ returned values. The minimum df.prior value returned by
+ eBayes(fit, robust=TRUE) may be slightly smaller than previously.
+
+\item
+ tmixture.vector() now handles unequal df in a better way. This
+ will be seen in slightly changed B-statistics from topTable when
+ the residual or prior df are not all equal.
+
+\item
+ If a targets data.frame is provided to read.maimages() through the
+ 'files' argument, then read.maimages() now overwrites its rownames
+ on output to agree with the column names of the RGList object.
+
+\item
+ More explicit handling of namespaces.
+ Functions needed from the grDevices, graphics, stats and utils
+ packages are now imported explicitly into the NAMESPACE, avoiding any possible masking by other packages.
+ goana(), alias2Symbol(), alias2SymbolTable() now load the name spaces of GO.db
+ the relevant organism package org.XX.eg.db instead of loading the
+ packages.
+ This keeps the user search path cleaner.
+
+\item
+ Various minor bug fixes and improvements to error messages.
+
+}}
+
+
+\section{Version 3.26.0}{\itemize{
+
+\item
+ New functions kegga() and topKEGG() to conduct KEGG pathway over-representation analyses.
+ kegga() is an S3 generic function with a method for MArrayLM fitted model objects as well a default method.
+ It uses the KEGGREST package to access up-to-date KEGG pathways using an internet connection.
+
+\item
+ goana() with trend now chooses the span of the tricube smoother
+ based on the number of differentially expressed genes. Larger
+ spans give smoother trends and are used for smaller numbers of
+ genes. Arguments 'species' and 'plot' are no longer part of
+ goana.MArrayLM() and are instead passed if present in a call to goana.default().
+
+\item
+ topGO() has a new argument 'truncate.term' to limit number of
+ characters appearing in a column of the output data.frame.
+
+\item
+ Improved power for the romer() function. The empirical Bayes
+ shrinkage of the residuals, removed on 4 April, has been
+ re-instated in a corrected form. It can optionally be turned off
+ using the 'shrink.resid' argument.
+
+\item
+ camera() has a new argument inter.gene.cor. This allows a preset
+ inter gene correlation to be set for each test test, resulting in
+ potentially less conservative tests.
+ camera() no longer gives NA results when a gene has a residual
+ variance exactly zero.
+
+\item
+ mroast() now sorts significant gene sets by proportion of genes
+ changing as well as by p-value and set size. When p-values are
+ equal, sets with a higher proportion of changing genes will now
+ rank higher in ordered results.
+ The contrast argument of roast() can now optionally be a character
+ string giving the name of a column of the design matrix.
+
+\item
+ fry() now produces "mixed" (non-directional) p-values and FDRs as
+ well as directional p-values and FDRs.
+ Bug fix for fry() when the gene set has a small number of genes.
+
+\item
+ camera(), fry(), mroast() and romer() now give more informative
+ error messages when the index list is empty.
+
+\item
+ Speedup for ids2indices().
+
+\item
+ tricubeMovingAverage() has new argument 'power' and old argument
+ 'full.length' is removed.
+ tricubeMovingAverage() now checks whether span is correctly
+ specified. It now gives identical results for any span>=1 or any
+ span<=0.
+
+\item
+ lmFit() and friends now always produce an output component 'rank'
+ that returns the column rank of the design matrix. Previously it
+ depended on the options used whether this component was given.
+
+\item
+ lmFit() now assumes that if the 'object' of expression values is a simple vector, that it represents
+ a single gene rather than a single sample.
+
+\item
+ lmFit() now always treats infinite expression values as NA. In the past,
+ this was done somewhat inconsistently.
+
+\item
+ Checking for NA coefficients in lmFit() is now done more
+ efficiently, meaning that lmFit() is now faster on large datasets.
+
+\item
+ More informative error message from lmscFit() when correlation is
+ NA.
+
+\item
+ getEAWP() now works on all eSet objects, not just ExpressionSet,
+ provided there is an data element called "exprs".
+ This allows all the linear modelling functions to handle general eSet objects.
+
+\item
+ plotMD() can now accept arguments 'array' or 'coef' as appropriate
+ as synonyms for 'column'. This may help users transitioning from
+ plotMA ot plotMD.
+
+\item
+ Arguments hi.pch, hi.col and hi.cex of plotWithHighlights() renamed
+ to hl.pch, hl.col and hl.cex. ('hl' is short for 'highlight'.)
+
+\item
+ New 'tolerance' argument for read.idat() to allow manifest and idat
+ to differ by a certain number of probes
+
+\item
+ The 'genes' argument of controlStatus() now accepts EListRaw or
+ Elist objects.
+
+\item
+ All data classes (EListRaw, EList, RGList, MAList, MArrayLM) now
+ require two arguments for subsetting. This restores the behavior
+ in limma version 3.19.10 and earlier. In limma 3.19.11 through
+ 3.25.13, subsetting by one argument y[i] was equivalent to y[i,].
+ This will now give an error message.
+
+\item
+ The dimnames<- method for EListRaw objects now sets rownames for
+ the background matrix Eb as well as for the foreground matrix E.
+
+\item
+ avearrays() now gives a more informative error message if input
+ argument 'x' is not a matrix.
+
+\item
+ auROC() now allows for tied stat values.
+
+\item
+ Bug fix to write.fit() when method="global" and fit contains more
+ than one column.
+}}
+
+
\section{Version 3.24.0}{\itemize{
\item
Limma review article (Ritchie et al, Nucleic Acids Research, 2015) is now published, and becomes the primary citation for limma.
@@ -38,6 +261,8 @@ camera(), roast(), mroast() and romer() now give a warning if
romer() is now an S3 generic function.
Speed improvements for romer() when the "mean" or "floormean" set
statistics are used.
+The empirical Bayes shrinkage of the residuals, introduced 5 October 2009, has been removed --
+this means that romer() results will revert to being somewhat more conservative.
\item
topGO() can now sort gene ontology results by several p-value
diff --git a/inst/doc/changelog.txt b/inst/doc/changelog.txt
index 5fb5e36..66c0037 100755
--- a/inst/doc/changelog.txt
+++ b/inst/doc/changelog.txt
@@ -1,4 +1,226 @@
-30 Oct 2015: limma 3.26.1
+15 Jun 2016: limma 3.28.10
+
+- new argument 'annotation' for read.idat(), to allow users to read
+ any required columns from the manifest file.
+
+15 Jun 2016: limma 3.28.9
+
+- normalizeBetweenArrays() now checks whether the input 'object' is a
+ data.frame, and converts to a matrix if possible.
+
+13 Jun 2016: limma 3.28.8
+
+- new function detectionPValues() to compute detection p-values from
+ negative control probes.
+
+12 Jun 2016: limma 3.28.7
+
+- new argument 'verbose' for read.idat().
+
+- Example code added to duplicateCorrelation.Rd.
+
+ 8 Jun 2016: limma 3.28.6
+
+- The default settings for the 'inter.gene.cor' and 'use.neg.cor'
+ arguments of camera() have been changed. camera() now uses by
+ default a preset value for the inter-gene correlation. This has the
+ effect that it tends to rank highly co-regulated,
+ biologically-interpretable more highly than before.
+
+22 May 2016: limma 3.28.5
+
+- New function fitmixture() to estimate genewise fold changes and
+ expression values from mixture experiments. Previously this
+ funtion was only available as part of the illumina package available
+ from http://bioinf.wehi.edu.au/illumina
+
+- New function logcosh() to compute log(cosh(x)) accurately without
+ floating underflow or overflow.
+
+11 May 2016: limma 3.28.4
+
+- Slight change to fry() standardize="posterior.sd" method.
+
+11 May 2016: limma 3.28.3
+
+- New flexibility for the roast() and mroast() functions, similar to
+ that previously implemented for fry(). The index vector for each
+ gene set can now be a data.frame, allowing each gene set to have
+ its own set of gene weights. The indices can now optionally be a
+ vector of gene names instead of a vector of indices. roast() and
+ mroast() now support the robust empirical Bayes option of
+ squeezeVar(). roast() and mroast() can now accept, via '...', any
+ argument that would be normally passed to lmFit() or eBayes().
+
+- bug fix to propTrueNull() when method="hist" and all the p-values
+ are less than 1/nbins.
+
+ 5 May 2016: limma 3.28.2
+
+- minor speed improvements for squeezeVar() and fry().
+
+ 4 May 2016: limma 3.28.1
+
+- unnecessary argument 'design' removed from fitted.MArrayLM.
+
+- bug fix to fry() with robust=TRUE.
+
+- bug fix to alias2Symbol() with expand.symbol=TRUE.
+
+ 4 May 2016: limma 3.28.0 (Bioconductor 3.3 Release Branch)
+
+29 April 2016: limma 3.27.20
+
+- update NEWS.Rd ready for Bioconductor 3.3 release.
+
+18 April 2016: limma 3.27.19
+
+- bug fix for fry(), which is some cases was reversing the sign of a
+ contrast effect.
+
+15 April 2016: limma 3.27.18
+
+- tmixture.vector() now handles unequal df in a better way. This
+ will be seen in slightly changed B-statistics from topTable when
+ the residual or prior df are not all equal.
+
+14 April 2016: limma 3.27.17
+
+- Improved capabilities and power for fry(). fry() has two new
+ arguments, 'geneid' and 'standardize'. The index argument of
+ fry() can now be a list of data.frames containing identifiers and
+ weights for each set. The options introduced by the standardize
+ argument allows fry() to be more powerful when genes have unequal
+ variances. fry() now accepts any arguments that would be suitable
+ for lmFit() or eBayes().
+
+13 April 2016: limma 3.27.16
+
+- alias2Symbol() and alias2SymbolTable() now load the name space of
+ the relevant organism package org.XX.eg.db instead of loading the
+ package.
+
+10 April 2016: limma 3.27.15
+
+- goana() now loads the name spaces of GO.db and org.XX.eg.db instead
+ of loading the packages. This keeps the user search path cleaner.
+
+22 March 2016: limma 3.27.14
+
+- bug fix for Active Proportions from roast() when gene.weights are
+ set.
+
+ 1 March 2016: limma 3.27.13
+
+- update references to Phipson et al (2016) article on robust
+ empirical Bayes.
+
+11 February 2016: limma 3.27.12
+
+- new argument 'save.plot' for voom().
+
+- bug fix for kegga() when 'covariate' or 'prior.prob' are not NULL
+ and when not all genes in the universe have annotation.
+
+27 January 2016: limma 3.27.11
+
+- kegga() now checks for missing values in 'gene.pathway' or
+ 'pathway.names'.
+
+26 January 2016: limma 3.27.10
+
+- kegga() gives better error message when there are no annotated genes
+ in the universe.
+
+- topKEGG() now breaks tied p-values by number of genes in pathway and
+ by name of pathway.
+
+- diffSplice() now returns genewise residual variances.
+
+- plotRLDF() now uses squeezeVar() to estimate how by how much the
+ within-group covariance matrix is moderated. It has new arguments
+ 'trend' and 'robust' as options to the empirical Bayes estimation
+ step. It also has new argument 'var.prior' to allow the prior
+ variances to be optionally supplied. Argument 'df.prior' can now be
+ a vector of genewise values.
+
+11 Jan 2016: limma 3.27.9
+
+- bug fix for kegga(), which wasn't passing the 'convert' argument
+ to getGeneKEGGLinks().
+
+ 8 Jan 2016: limma 3.27.8
+
+- goana() now supports species="Pt" (chimpanzee).
+
+- kegga() has several new arguments and now supports any species
+ supported by KEGG. kegga() can now accept annotation as a
+ data.frame, meaning that it can work from user-supplied annotation.
+
+- new functions getGeneKEGGLinks() and getKEGGPathwayNames() to get
+ pathway annotation from the rest.kegg website.
+
+- removeExt() has new 'sep' argument.
+
+29 Dec 2015: limma 3.27.7
+
+- roast(), mroast(), fry() and camera() now require the 'design'
+ matrix to be set. Previously the design matrix defaulted to a
+ single intercept column.
+
+- Bug fix to plotRLDF(), to ensure 'top' is included in output object
+ even if all the rows of data are used.
+
+- gls.series() now gives a more informative error message when the
+ correlation is equal to 1 or -1.
+
+ 1 Dec 2015: limma 3.27.6
+
+- Slightly improved algorithm for producing the df2.shrunk values
+ returned by fitFDistRobustly(). fitFDistRobustly() now returns
+ tail.p.value, prob.outlier and df2.outlier as well as previously
+ returned values. The minimum df.prior value returned by
+ eBayes(fit, robust=TRUE) may be slightly smaller than previously.
+
+- update 'sort' argument of fry() to be the same as for mroast().
+
+- edit help page for vennDiagram().
+
+15 Nov 2015: limma 3.27.5
+
+- update NEWS.Rd for Bioconductor 3.2 Release.
+
+- plotRLDF() now produces more complete output and the output is
+ more completely documented. It can now be used as a complete LDF
+ analysis for classification based on training data. The argument
+ 'main' has been removed as an explicit argument because it and
+ other plotting arguments are better passed using the ... facility.
+ New arguments 'ndim' and 'plot' have been added. The first allows
+ all possible discriminant functions to be completed even if only
+ two are plotted. The second allows dicriminant functions to be
+ computed without a plot.
+
+- If a targets data.frame is provided to read.maimages() through the
+ 'files' argument, then read.maimages() now overwrites its rownames
+ on output to agree with the column names of the RGList object.
+
+ 6 Nov 2015: limma 3.27.4
+
+- Bug fix to tricubeMovingAverage() when length of series is even
+ and span is close to 1.
+
+30 Oct 2015: limma 3.27.3
+
+- Add a check in camera() that dataset must contain as least three
+ genes.
+
+29 Oct 2015: limma 3.27.2
+
+- Fix to NAMESPACE. The function utils::winMenuAddItem(), which is
+ needed to add a GUI menu item for the limma User's Guide, is now
+ imported only conditionally in Windows.
+
+27 Oct 2015: limma 3.27.1
- Functions needed from the grDevices, graphics, stats and utils
packages are now imported explicitly into the NAMESPACE.
@@ -7,19 +229,17 @@
semitransparency of positional bars in some circumstances;
the default value for 'quantiles' is increased slightly.
-- Add a check in camera() that dataset must contain as least three
- genes.
-
-- camera() now gives explicit error when too few residual df.
-
- Bug fix to camera(), which was not accepting character values
for index when weights were non-NULL.
+- camera() now gives explicit error when too few residual df.
+
+14 Oct 2015: limma 3.27.0 (Bioconductor 3.3 Developmental Branch)
14 Oct 2015: limma 3.26.0 (Bioconductor 3.2 Release Branch)
7 Oct 2015: limma 3.25.18
-- mroast() now sorts significant gene sets up proportion of genes
+- mroast() now sorts significant gene sets by proportion of genes
changing as well as p-values and set size. When p-values are
equal, sets with a higher proportion of changing genes will now
rank higher in ordered results.
diff --git a/inst/doc/intro.pdf b/inst/doc/intro.pdf
index 051de65..c3ef58e 100644
Binary files a/inst/doc/intro.pdf and b/inst/doc/intro.pdf differ
diff --git a/man/03reading.Rd b/man/03reading.Rd
index 1eb4283..e55cacd 100644
--- a/man/03reading.Rd
+++ b/man/03reading.Rd
@@ -26,6 +26,7 @@ There are also a series of utility functions which read the header information f
and produces an \code{ElistRaw} object.
\code{\link{read.idat}} reads Illumina files in IDAT format, and produces an \code{EListRaw} object.
+\code{\link{detectionPValues}} can be used to add detection p-values.
The function \link{as.MAList} can be used to convert a \code{marrayNorm} object to an \code{MAList} object if the data was read and normalized using the marray and marrayNorm packages.
}
diff --git a/man/10GeneSetTests.Rd b/man/10GeneSetTests.Rd
index 6670298..116dfab 100644
--- a/man/10GeneSetTests.Rd
+++ b/man/10GeneSetTests.Rd
@@ -13,7 +13,7 @@ This page gives an overview of the LIMMA functions for gene set testing and path
Self-contained gene set testing for many sets.}
\item{ \code{\link{fry}} }{
- Fast version of \code{mroast} when heteroscedasticity of genes can be ignored.}
+ Fast approximation to \code{mroast}, especially useful when heteroscedasticity of genes can be ignored.}
\item{ \code{\link{camera}} }{
Competitive gene set testing.}
diff --git a/man/alias2Symbol.Rd b/man/alias2Symbol.Rd
index 29d93c5..776b3d0 100644
--- a/man/alias2Symbol.Rd
+++ b/man/alias2Symbol.Rd
@@ -13,12 +13,15 @@ alias2SymbolTable(alias, species = "Hs")
\item{alias}{character vector of gene aliases}
\item{species}{character string specifying the species.
Possible values are \code{"Dm"}, \code{"Hs"}, \code{"Mm"} or \code{"Rn"}.}
- \item{expand.symbols}{logical, should those elements of \code{alias} which are already official symbols be expanded if they are aliases for other symbols.}
+ \item{expand.symbols}{logical.
+ This affects those elements of \code{alias} that are the official gene symbol for one gene and also an alias for another gene.
+ If \code{FALSE}, then these elements will just return themselves.
+ If \code{TRUE}, then all the genes for which they are aliases will be returned.}
}
\details{
Aliases are mapped via NCBI Entrez Gene identity numbers using Bioconductor organism packages.
-Species are \code{"Dm"} for fly, \code{"Hs"} for human, \code{"Mm"} for mouse and \code{"Rn"} for rat.
-The user needs to have the appropriate Bioconductor organism package installed.
+Possible species are \code{"Dm"} for fly, \code{"Hs"} for human, \code{"Mm"} for mouse and \code{"Rn"} for rat.
+The appropriate Bioconductor organism package is assumed to be installed.
\code{alias2Symbol} maps a set of aliases to a set of symbols, without necessarily preserving order.
The output vector may be longer or shorter than the original vector, because some aliases might not be found and some aliases may map to more than one symbol.
@@ -29,7 +32,7 @@ If an alias maps to more than one symbol, then the first one found is returned.
Character vector of gene symbols.
\code{alias2SymbolTable} returns a vector of the same length and order as \code{alias}, including \code{NA} values where no gene symbol was found.
-\code{alias2Symbol} returns an unordered vector which may be longer or shorter than \code{alias}.
+\code{alias2Symbol} returns an unordered vector that may be longer or shorter than \code{alias}.
}
\author{Gordon Smyth and Yifang Hu}
\seealso{
@@ -37,7 +40,8 @@ This function is often used to assist gene set testing, see
\link{10.GeneSetTests}.
}
\examples{
-if(require("org.Hs.eg.db")) alias2Symbol(c("PUMA","NOXA","BIM"))
+alias2Symbol(c("PUMA","NOXA","BIM"), species="Hs")
+alias2Symbol("RS1", expand=TRUE)
}
\keyword{gene annotation}
diff --git a/man/camera.Rd b/man/camera.Rd
index fbb5290..ce35342 100644
--- a/man/camera.Rd
+++ b/man/camera.Rd
@@ -7,9 +7,9 @@
Test whether a set of genes is highly ranked relative to other genes in terms of differential expression, accounting for inter-gene correlation.
}
\usage{
-\S3method{camera}{default}(y, index, design, contrast=ncol(design), weights=NULL,
- use.ranks=FALSE, allow.neg.cor=TRUE, inter.gene.cor=NULL, trend.var=FALSE,
- sort=TRUE, \dots)
+\method{camera}{default}(y, index, design, contrast = ncol(design), weights = NULL,
+ use.ranks = FALSE, allow.neg.cor=FALSE, inter.gene.cor=0.01, trend.var = FALSE,
+ sort = TRUE, \dots)
interGeneCorrelation(y, design)
}
\arguments{
@@ -23,7 +23,7 @@ interGeneCorrelation(y, design)
\item{weights}{can be a numeric matrix of individual weights, of same size as \code{y}, or a numeric vector of array weights with length equal to \code{ncol(y)}, or a numeric vector of gene weights with length equal to \code{nrow(y)}.}
\item{use.ranks}{do a rank-based test (\code{TRUE}) or a parametric test (\code{FALSE})?}
\item{allow.neg.cor}{should reduced variance inflation factors be allowed for negative correlations?}
- \item{inter.gene.cor}{numeric, optional preset value for the inter gene correlation within tested sets. If not \code{NULL}, then this value will over-ride estimation of the inter gene correlations.}
+ \item{inter.gene.cor}{numeric, optional preset value for the inter-gene correlation within tested sets. If \code{NA} or \code{NULL}, then an inter-gene correlation will be estimated for each tested set.}
\item{trend.var}{logical, should an empirical Bayes trend be estimated? See \code{\link{eBayes}} for details.}
\item{sort}{logical, should the results be sorted by p-value?}
\item{\dots}{other arguments are not currently used}
@@ -44,21 +44,32 @@ This allows users to focus on differential expression for any coefficient or con
The inflation factor depends on estimated genewise correlation and the number of genes in the gene set.
By default, \code{camera} uses \code{interGeneCorrelation} to estimate the mean pair-wise correlation within each set of genes.
-\code{camera} can be used with a small preset correlation value, say \code{inter.gene.cor=0.01}.
-This produces a less conservative test that performs well in many practical situations.
+\code{camera} can alternatively be used with a preset correlation specified by \code{inter.gene.cor} that is shared by all sets.
+This usually works best with a small value, say \code{inter.gene.cor=0.01}.
+
+If \code{interGeneCorrelation=NA}, then \code{camera} will estimate the inter-gene correlation for each set.
+In this mode, \code{camera} gives rigorous error rate control for all sample sizes and all gene sets.
+However, in this mode, highly co-regulated gene sets that are biological interpretable may not always be ranked at the top of the list.
+
+With \code{interGeneCorrelation=0.01}, \code{camera} will rank biologically interpetable sets more highly.
+This gives a useful compromise between strict error rate control and interpretable gene set rankings.
}
+\note{The default settings for \code{inter.gene.cor} and \code{allow.neg.cor} were changed to the current values in limma 3.29.6.
+Previously, the default was to estimate an inter-gene correlation for each set.
+To reproduce the earlier default, use \code{allow.neg.cor=TRUE} and \code{inter.gene.cor=NA}.}
+
\value{
\code{camera} returns a data.frame with a row for each set and the following columns:
-\item{NGenes}{number of genes in set}
-\item{Correlation}{inter-gene correlation}
-\item{Direction}{direction of change (\code{"Up"} or \code{"Down"})}
-\item{PValue}{two-tailed p-value}
-\item{FDR}{Benjamini and Hochberg FDR adjusted p-value}
+\item{NGenes}{number of genes in set.}
+\item{Correlation}{inter-gene correlation (only included if the \code{inter.gene.cor} was not preset).}
+\item{Direction}{direction of change (\code{"Up"} or \code{"Down"}).}
+\item{PValue}{two-tailed p-value.}
+\item{FDR}{Benjamini and Hochberg FDR adjusted p-value.}
\code{interGeneCorrelation} returns a list with components:
-\item{vif}{variance inflation factor}
-\item{correlation}{inter-gene correlation}
+\item{vif}{variance inflation factor.}
+\item{correlation}{inter-gene correlation.}
}
\author{Di Wu and Gordon Smyth}
@@ -79,6 +90,7 @@ Analyzing gene expression data in terms of gene sets: methodological issues.
\code{\link{rankSumTestWithCorrelation}},
\code{\link{geneSetTest}},
\code{\link{roast}},
+\code{\link{fry}},
\code{\link{romer}},
\code{\link{ids2indices}}.
@@ -99,7 +111,7 @@ index2 <- 21:40
camera(y, index1, design)
camera(y, index2, design)
-camera(y, list(set1=index1,set2=index2), design)
+camera(y, list(set1=index1,set2=index2), design, inter.gene.cor=NA)
camera(y, list(set1=index1,set2=index2), design, inter.gene.cor=0.01)
}
diff --git a/man/detectionPValue.Rd b/man/detectionPValue.Rd
new file mode 100644
index 0000000..aca5353
--- /dev/null
+++ b/man/detectionPValue.Rd
@@ -0,0 +1,58 @@
+\name{detectionPValues}
+\alias{detectionPValues.default}
+\alias{detectionPValues.EListRaw}
+\alias{detectionPValues}
+
+\title{Detection P-Values from Negative Controls}
+
+\description{Compute the proportion of negative controls greater than each observed expression value.
+Particularly useful for Illumina BeadChips.}
+
+\usage{
+\method{detectionPValues}{EListRaw}(x, status = NULL, \dots)
+\method{detectionPValues}{default}(x, status, negctrl = "negative", \dots)
+}
+
+\arguments{
+ \item{x}{object of class \code{EListRaw} or a numeric \code{matrix} containing raw intensities for regular and control probes from a series of microarrays.}
+ \item{status}{character vector giving probe types. Defaults to \code{x$genes$Status} if \code{x} is an \code{EListRaw} object.}
+ \item{negctrl}{character string identifier for negative control probes.}
+ \item{\dots}{other arguments are not currently used.}
+ }
+
+\details{
+The rows of \code{x} for which \code{status == negctrl} are assumed to correspond to negative control probes.
+
+For each column of \code{x}, the detection p-values are defined as \code{(N.eq/2 + N.gt) / N.neg}, where \code{N.gt} is the number of negative controls with expression greater than the observed value, \code{N.eq} is the number of negative controls with expression equal to the observed value, and \code{N.neg} is the total number of negative controls.
+}
+
+\value{numeric matrix of same dimensions as \code{x} containing detection p-values.}
+
+\references{
+Shi, W, de Graaf, C, Kinkel, S, Achtman, A, Baldwin, T, Schofield, L, Scott, H, Hilton, D, Smyth, GK (2010).
+Estimating the proportion of microarray probes expressed in an RNA sample.
+\emph{Nucleic Acids Research} 38, 2168-2176.
+\url{http://nar.oxfordjournals.org/content/38/7/2168}
+}
+
+\author{Gordon Smyth}
+
+\seealso{
+ An overview of LIMMA functions to read expression data is given in \link{03.ReadingData}.
+
+ \code{\link{read.idat}} reads Illumina BeadChip expression data from binary IDAT files.
+
+ \code{\link{neqc}} performs normexp background correction and quantile normalization aided by control probes.
+}
+
+\examples{
+\dontrun{
+# Read Illumina binary IDAT files
+x <- read.idat(idat, bgx)
+x$genes$DectionPValue <- detectionPValues(x)
+y <- neqc(x)
+}
+}
+
+\keyword{background correction}
+\keyword{illumina beadchips}
diff --git a/man/dupcor.Rd b/man/dupcor.Rd
index 6054cc9..ee170ad 100755
--- a/man/dupcor.Rd
+++ b/man/dupcor.Rd
@@ -52,16 +52,40 @@ An overview of linear model functions in limma is given by \link{06.LinearModels
\author{Gordon Smyth}
\references{
Smyth, G. K., Michaud, J., and Scott, H. (2005). The use of within-array replicate spots for assessing differential expression in microarray experiments. \emph{Bioinformatics} 21(9), 2067-2075.
-\url{http://www.statsci.org/smyth/pubs/dupcor.pdf}
+[\url{http://bioinformatics.oxfordjournals.org/content/21/9/2067}]
+[Preprint with corrections: \url{http://www.statsci.org/smyth/pubs/dupcor.pdf}]
}
+
\examples{
-# Also see lmFit examples
+# Simulate gene expression data for 100 probes and 6 microarrays
+# Microarray are in two groups
+# First two probes are more highly expressed in second group
+# Std deviations vary between genes with prior df=4
+sd <- 0.3*sqrt(4/rchisq(100,df=4))
+y <- matrix(rnorm(100*6,sd=sd),100,6)
+rownames(y) <- paste("Gene",1:100)
+y[1:2,4:6] <- y[1:2,4:6] + 2
+design <- cbind(Grp1=1,Grp2vs1=c(0,0,0,1,1,1))
+options(digits=3)
-\dontrun{
-corfit <- duplicateCorrelation(MA, ndups=2, design)
-all.correlations <- tanh(corfit$atanh.correlations)
-boxplot(all.correlations)
-fit <- lmFit(MA, design, ndups=2, correlation=corfit$consensus)
-}
+# Fit with correlated arrays
+# Suppose each pair of arrays is a block
+block <- c(1,1,2,2,3,3)
+dupcor <- duplicateCorrelation(y,design,block=block)
+dupcor$consensus.correlation
+fit1 <- lmFit(y,design,block=block,correlation=dupcor$consensus)
+fit1 <- eBayes(fit1)
+topTable(fit1,coef=2)
+
+# Fit with duplicate probes
+# Suppose two side-by-side duplicates of each gene
+rownames(y) <- paste("Gene",rep(1:50,each=2))
+dupcor <- duplicateCorrelation(y,design,ndups=2)
+dupcor$consensus.correlation
+fit2 <- lmFit(y,design,ndups=2,correlation=dupcor$consensus)
+dim(fit2)
+fit2 <- eBayes(fit2)
+topTable(fit2,coef=2)
}
+
\keyword{multivariate}
diff --git a/man/ebayes.Rd b/man/ebayes.Rd
index fdd2088..fd51bf1 100755
--- a/man/ebayes.Rd
+++ b/man/ebayes.Rd
@@ -71,9 +71,13 @@ If \code{trend=TRUE} then an intensity-dependent trend is fitted to the prior va
Specifically, \code{squeezeVar} is called with the \code{covariate} equal to \code{Amean}, the average log2-intensity for each gene.
See \code{\link{squeezeVar}} for more details.
-If \code{robust=TRUE} then the robust empirical Bayes procedure of Phipson et al (2013) is used.
+If \code{robust=TRUE} then the robust empirical Bayes procedure of Phipson et al (2016) is used.
See \code{\link{squeezeVar}} for more details.
}
+\note{
+The algorithm used by \code{eBayes} and \code{treat} with \code{robust=TRUE} was revised slightly in limma 3.27.6.
+The minimum \code{df.prior} returned may be slightly smaller than previously.
+}
\seealso{
\code{\link{squeezeVar}}, \code{\link{fitFDist}}, \code{\link{tmixture.matrix}}.
@@ -89,10 +93,10 @@ Testing significance relative to a fold-change threshold is a TREAT.
Loennstedt, I., and Speed, T. P. (2002). Replicated microarray data. \emph{Statistica Sinica} \bold{12}, 31-46.
-Phipson, B, Lee, S, Majewski, IJ, Alexander, WS, and Smyth, GK (2013).
-Empirical Bayes in the presence of exceptional cases, with application to microarray data.
-Bioinformatics Division, Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia.
-\url{http://www.statsci.org/smyth/pubs/RobustEBayesPreprint.pdf}
+Phipson, B, Lee, S, Majewski, IJ, Alexander, WS, and Smyth, GK (2016).
+Robust hyperparameter estimation protects against hypervariable genes and improves power to detect differential expression.
+\emph{Annals of Applied Statistics} 10.
+\url{http://arxiv.org/abs/1602.08678}
Smyth, G. K. (2004). Linear models and empirical Bayes methods for assessing differential expression in microarray experiments.
\emph{Statistical Applications in Genetics and Molecular Biology}, Volume \bold{3}, Article 3.
diff --git a/man/fitfdist.Rd b/man/fitfdist.Rd
index d7a14f5..3db08b8 100755
--- a/man/fitfdist.Rd
+++ b/man/fitfdist.Rd
@@ -24,13 +24,26 @@ The parameters are estimated using the method of moments, specifically from the
\code{fitFDistRobustly} is similar to \code{fitFDist} except that it computes the moments of the Winsorized values of \code{x}, making it robust against left and right outliers.
Larger values for \code{winsor.tail.p} produce more robustness but less efficiency.
-The robust method is described by Phipson (2013).
+The robust method is described by Phipson (2013) and Phipson et al (2016).
+
+As well as estimating the F-distribution for the bulk of the cases, i.e., with outliers discounted, \code{fitFDistRobustly} also returns an estimated F-distribution with reduced df2 that might be appropriate for each outlier case.
+}
+\note{
+The algorithm used by \code{fitFDistRobustly} was revised slightly in limma 3.27.6.
+The \code{prob.outlier} value, which is the lower bound for \code{df2.shrunk}, may be slightly smaller than previously.
}
\value{
-A list containing the components
+\code{fitFDist} produces a list with the following components:
\item{scale}{scale factor for F-distribution. A vector if \code{covariate} is non-\code{NULL}, otherwise a scalar.}
- \item{df2}{the second degrees of freedom of the F-distribution.}
+ \item{df2}{the second degrees of freedom of the fitted F-distribution.}
+
+\code{fitFDistRobustly} returns the following components as well:
+ \item{tail.p.value}{right tail probability of the scaled F-distribution for each \code{x} value.}
+ \item{prob.outlier}{posterior probability that each case is an outlier relative to the scaled F-distribution with degrees of freedom \code{df1} and \code{df2}.}
+ \item{df2.outlier}{the second degrees of freedom associated with extreme outlier cases.}
+ \item{df2.shrunk}{numeric vector of values for the second degrees of freedom, with shrunk values for outliers. Most values are equal to \code{df2}, but outliers have reduced values depending on how extreme each case is. All values lie between \code{df2.outlier} and \code{df2}.}
}
+
\author{Gordon Smyth and Belinda Phipson}
\references{
@@ -45,13 +58,14 @@ Phipson, B. (2013).
PhD Thesis. University of Melbourne, Australia.
\url{http://repository.unimelb.edu.au/10187/17614}
-Phipson, B., and Smyth, G. K. (2013).
-\emph{Robust empirical Bayes estimation protetcts against hyper-variable genes and improves power to detect differential expression in RNA-seq data}.
-Bioinformatics Division, Walter and Eliza Hall Institute of Medical Research, Australia
+Phipson, B, Lee, S, Majewski, IJ, Alexander, WS, and Smyth, GK (2016).
+Robust hyperparameter estimation protects against hypervariable genes and improves power to detect differential expression.
+\emph{Annals of Applied Statistics} 10.
+\url{http://arxiv.org/abs/1602.08678}
}
\seealso{
-This function is called by \code{\link{squeezeVar}}, and hence by \code{\link{ebayes}} and \code{\link{eBayes}}.
+This function is called by \code{\link{squeezeVar}}, which in turn is called by \code{\link{ebayes}}, \code{\link{eBayes}} and \code{\link{treat}}.
This function calls \code{\link{trigammaInverse}}.
}
diff --git a/man/fitmixture.Rd b/man/fitmixture.Rd
new file mode 100644
index 0000000..5a7db5b
--- /dev/null
+++ b/man/fitmixture.Rd
@@ -0,0 +1,66 @@
+\name{fitmixture}
+\alias{fitmixture}
+\title{Fit Mixture Model by Non-Linear Least Squares}
+\description{Fit Mixture Model by Non-Linear Least Squares}
+
+\usage{
+fitmixture(log2e, mixprop, niter = 4, trace = FALSE)
+}
+
+\arguments{
+ \item{log2e}{a numeric matrix containing log2 expression values. Rows correspond to probes for genes and columns to RNA samples.}
+ \item{mixprop}{a vector of length \code{ncol(log2e)} giving the mixing proportion (between 0 and 1) for each sample.}
+ \item{niter}{integer number of iterations.}
+ \item{trace}{logical. If \code{TRUE}, summary working estimates are output from each iteration.}
+}
+
+\details{
+A mixture experiment is one in which two reference RNA sources are mixed in different proportions to create experimental samples.
+Mixture experiments have been used to evaluate genomic technologies and analysis methods (Holloway et al, 2006).
+This function uses all the data for each gene to estimate the expression level of the gene in each of two pure samples.
+
+The function fits a nonlinear mixture model to the log2 expression values for each gene.
+The expected values of \code{log2e} for each gene are assumed to be of the form
+\code{log2( mixprop*Y1 + (1-mixprop)*Y2 )}
+where \code{Y1} and \code{Y2} are the expression levels of the gene in the two reference samples being mixed.
+The \code{mixprop} values are the same for each gene but \code{Y1} and \code{Y2} are specific to the gene.
+The function returns the estimated values \code{A=0.5*log2(Y1*Y2)} and \code{M=log2(Y2/Y1)} for each gene.
+
+The nonlinear estimation algorithm implemented in \code{fitmixture} uses a nested Gauss-Newton iteration (Smyth, 1996).
+It is fully vectorized so that the estimation is done for all genes simultaneously.
+}
+
+\value{List with three components:
+\item{A}{numeric vector giving the estimated average log2 expression of the two reference samples for each gene}
+\item{M}{numeric vector giving estimated log-ratio of expression between the two reference samples for each gene}
+\item{stdev}{standard deviation of the residual term in the mixture model for each gene}
+}
+
+\author{Gordon K Smyth}
+
+\references{
+Holloway, A. J., Oshlack, A., Diyagama, D. S., Bowtell, D. D. L., and Smyth, G. K. (2006).
+Statistical analysis of an RNA titration series evaluates microarray precision and sensitivity on a whole-array basis.
+\emph{BMC Bioinformatics} 7, Article 511.
+\url{http://www.biomedcentral.com/1471-2105/7/511}
+
+Smyth, G. K. (1996).
+Partitioned algorithms for maximum likelihood and other nonlinear estimation.
+\emph{Statistics and Computing}, 6, 201-216.
+\url{http://www.statsci.org/smyth/pubs/partitio.pdf}
+}
+
+\examples{
+ngenes <- 100
+TrueY1 <- rexp(ngenes)
+TrueY2 <- rexp(ngenes)
+mixprop <- matrix(c(0,0.25,0.75,1),1,4)
+TrueExpr <- TrueY1 %*% mixprop + TrueY2 %*% (1-mixprop)
+
+log2e <- log2(TrueExpr) + matrix(rnorm(ngenes*4),ngenes,4)*0.1
+out <- fitmixture(log2e,mixprop)
+
+# Plot true vs estimated log-ratios
+plot(log2(TrueY1/TrueY2), out$M)
+}
+
diff --git a/man/fitted.MArrayLM.Rd b/man/fitted.MArrayLM.Rd
index f191974..dc482cb 100644
--- a/man/fitted.MArrayLM.Rd
+++ b/man/fitted.MArrayLM.Rd
@@ -3,12 +3,11 @@
\concept{regression}
\title{Fitted Values Method for MArrayLM Fits}
\usage{
-\S3method{fitted}{MArrayLM}(object, design = object$design, \dots)
+\S3method{fitted}{MArrayLM}(object, \dots)
}
\arguments{
\item{object}{a fitted object of class inheriting from \code{"MArrayLM"}.}
- \item{design}{numeric design matrix.}
- \item{\dots}{further arguments passed to or from other methods.}
+ \item{\dots}{other arguments are not currently used.}
}
\description{
Obtains fitted values from a fitted microarray linear model object.
diff --git a/man/getEAWP.Rd b/man/getEAWP.Rd
index 84ff6fc..2ee17b6 100644
--- a/man/getEAWP.Rd
+++ b/man/getEAWP.Rd
@@ -1,32 +1,30 @@
\name{getEAWP}
\alias{getEAWP}
-\title{Extract Basic Data from Microarray Data Objects}
+\title{Extract Basic Data from Expression Data Objects}
\description{
-Given a microarray data object of any known class, get the expression values, weights, probe annotation and A-values, which are needed for linear modelling.
+Given an expression data object of any known class, get the expression values, weights, probe annotation and A-values that are needed for linear modelling.
This function is called by the linear modelling functions in LIMMA.
}
\usage{
getEAWP(object)
}
\arguments{
- \item{object}{a microarray data object.
- An object of class \code{MAList}, \code{EList}, \code{marrayNorm}, \code{PLMset}, \code{vsn}, or any class inheriting from \code{ExpressionSet}, or any object that can be coerced to a numeric matrix.}
+ \item{object}{any matrix-like object containing log-expression values.
+ Can be an object of class \code{MAList}, \code{EList}, \code{marrayNorm}, \code{PLMset}, \code{vsn}, or any class inheriting from \code{ExpressionSet}, or any object that can be coerced to a numeric matrix.}
}
\details{
-In the case of two-color objects, the \code{Amean} is computed from the matrix of A-values.
-For single-channel objects, \code{Amean} is computed from the matrix of expression vales.
-\code{PLMset}, \code{vsn} and \code{ExpressionSet} are assumed to be single-channel for this purpose.
+Rows correspond to probes and columns to RNA samples.
-If \code{object} is a matrix, it is assumed to contain log-intensities if the values are all positive and log-ratios otherwise.
-\code{Amean} is computed in the former case but not the latter.
+In the case of two-color microarray data objects (\code{MAList} or \code{marrayNorm}), \code{Amean} is the vector of row means of the matrix of A-values.
+For other data objects, \code{Amean} is the vector of row means of the matrix of expression values.
-From April 2013, the output \code{exprs} matrix is ensured to have unique row names.
-If \code{object} has no row names, then the output row names of \code{exprs} are 1 to the number of rows.
-If \code{object} has row names but with duplicated names, then row names of \code{exprs} are set to 1 up to the number of rows and the original row names are preserved in the \code{ID} column of \code{probes}.
+From April 2013, the rownames of the output \code{exprs} matrix are required to be unique.
+If \code{object} has no row names, then the output rownames of \code{exprs} are \code{1:nrow(object)}.
+If \code{object} has row names but with duplicated names, then the rownames of \code{exprs} are set to \code{1:nrow(object)} and the original row names are preserved in the \code{ID} column of \code{probes}.
}
\value{
A list with components
-\item{exprs}{numeric matrix of log-ratios or log-intensities}
+\item{exprs}{numeric matrix of log-ratios, log-intensities or log-expression values}
\item{weights}{numeric matrix of weights}
\item{probes}{data.frame of probe-annotation}
\item{Amean}{numeric vector of average log-expression for each probe}
@@ -37,4 +35,3 @@ The other components will be \code{NULL} if not found in the input object.
\seealso{
\link{02.Classes} gives an overview of data classes used in LIMMA.
}
-\keyword{array}
diff --git a/man/gls.series.Rd b/man/gls.series.Rd
index 43ccc0b..e66b952 100755
--- a/man/gls.series.Rd
+++ b/man/gls.series.Rd
@@ -10,8 +10,8 @@ This is a utility function for \code{lmFit}.
\arguments{
\item{M}{numeric matrix containing log-ratio or log-expression values for a series of microarrays, rows correspond to genes and columns to arrays.}
\item{design}{numeric design matrix defining the linear model, with rows corresponding to arrays and columns to comparisons to be estimated. The number of rows must match the number of columns of \code{M}. Defaults to the unit vector meaning that the arrays are treated as replicates.}
- \item{ndups}{positive integer giving the number of times each gene is printed on an array. \code{nrow(M)} must be divisible by \code{ndups}.}
- \item{spacing}{the spacing between the rows of \code{M} corresponding to duplicate spots, \code{spacing=1} for consecutive spots}
+ \item{ndups}{positive integer giving the number of times each gene is printed on an array. \code{nrow(M)} must be divisible by \code{ndups}. Ignored if \code{block} is not \code{NULL}.}
+ \item{spacing}{the spacing between the rows of \code{M} corresponding to duplicate spots, \code{spacing=1} for consecutive spots. Ignored if \code{block} is not \code{NULL}.}
\item{block}{vector or factor specifying a blocking variable on the arrays.
Same length as \code{ncol(M)}.}
\item{correlation}{numeric value specifying the inter-duplicate or inter-block correlation.}
diff --git a/man/goana.Rd b/man/goana.Rd
index 4fccbfb..192f567 100644
--- a/man/goana.Rd
+++ b/man/goana.Rd
@@ -5,6 +5,8 @@
\alias{kegga}
\alias{kegga.default}
\alias{kegga.MArrayLM}
+\alias{getGeneKEGGLinks}
+\alias{getKEGGPathwayNames}
\title{Gene Ontology or KEGG Pathway Analysis}
\description{
@@ -16,8 +18,11 @@ Test for over-representation of gene ontology (GO) terms or KEGG pathways in one
\method{goana}{default}(de, universe = NULL, species = "Hs", prior.prob = NULL, covariate=NULL,
plot=FALSE, \dots)
\method{kegga}{MArrayLM}(de, coef = ncol(de), geneid = rownames(de), FDR = 0.05, trend = FALSE, \dots)
-\method{kegga}{default}(de, universe = NULL, species = "Hs", prior.prob = NULL, covariate=NULL,
- plot=FALSE, \dots)
+\method{kegga}{default}(de, universe = NULL, species = "Hs", species.KEGG = NULL, convert = FALSE,
+ gene.pathway = NULL, pathway.names = NULL,
+ prior.prob = NULL, covariate=NULL, plot=FALSE, \dots)
+getGeneKEGGLinks(species.KEGG = "hsa", convert = FALSE)
+getKEGGPathwayNames(species.KEGG = NULL, remove.qualifier = FALSE)
}
\arguments{
@@ -25,14 +30,22 @@ Test for over-representation of gene ontology (GO) terms or KEGG pathways in one
\item{coef}{column number or column name specifying for which coefficient or contrast differential expression should be assessed.}
\item{geneid}{Entrez Gene identifiers. Either a vector of length \code{nrow(de)} or the name of the column of \code{de$genes} containing the Entrez Gene IDs.}
\item{FDR}{false discovery rate cutoff for differentially expressed genes. Numeric value between 0 and 1.}
- \item{species}{species identifier. Possible values are \code{"Hs"}, \code{"Mm"}, \code{"Rn"} or \code{"Dm"}.}
+ \item{species}{two-letter Bioconductor species identifier. Possible values are \code{"Hs"}, \code{"Mm"}, \code{"Rn"}, \code{"Dm"} or \code{"Pt"}. Ignored if \code{species.KEGG} or is not \code{NULL} or if \code{gene.pathway} and \code{pathway.names} are not \code{NULL}.}
+ \item{species.KEGG}{three-letter KEGG species identifier. See \url{http://www.kegg.jp/kegg/catalog/org_list.html} or \url{http://rest.kegg.jp/list/organism} for possible values. Ignored if \code{gene.pathway} and \code{pathway.names} are not \code{NULL}.}
+ \item{convert}{if \code{TRUE} then KEGG gene identifiers will be converted to NCBI Entrez Gene identifiers. Note that KEGG IDs are the same as Entrez Gene IDs for most species anyway.}
+ \item{gene.pathway}{data.frame linking genes to pathways. First column gives gene IDs, second column gives pathway IDs. By default this is obtained automatically by \code{getGeneKEGGLinks(species.KEGG)}.}
+ \item{remove.qualifier}{if \code{TRUE}, the species qualifier will be removed from the pathway names.}
+ \item{pathway.names}{data.frame giving full names of pathways. First column gives pathway IDs, second column gives pathway names. By default this is obtained automatically using \code{getKEGGPathwayNames(species.KEGG, remove=TRUE)}.}
\item{trend}{adjust analysis for gene length or abundance?
Can be logical, or a numeric vector of covariate values, or the name of the column of \code{de$genes} containing the covariate values.
If \code{TRUE}, then \code{de$Amean} is used as the covariate.}
\item{universe}{vector specifying the set of Entrez Gene identifiers to be the background universe.
If \code{NULL} then all Entrez Gene IDs associated with any gene ontology term will be used as the universe.}
- \item{prior.prob}{optional numeric vector of the same length as \code{universe} giving the prior probability that each gene in the universe appears in a gene set. Will be computed from \code{covariate} if the latter is provided.}
- \item{covariate}{optional numeric vector of the same length as \code{universe} giving a covariate against which \code{prior.prob} should be computed.}
+ \item{prior.prob}{optional numeric vector of the same length as \code{universe} giving the prior probability that each gene in the universe appears in a gene set.
+ Will be computed from \code{covariate} if the latter is provided.
+ Ignored if \code{universe} is \code{NULL}.}
+ \item{covariate}{optional numeric vector of the same length as \code{universe} giving a covariate against which \code{prior.prob} should be computed.
+ Ignored if \code{universe} is \code{NULL}.}
\item{plot}{logical, should the \code{prior.prob} vs \code{covariate} trend be plotted?}
\item{\dots}{any other arguments in a call to the \code{MArrayLM} method are passed to the default method.}
}
@@ -43,9 +56,18 @@ The default method accepts the gene list as a vector of gene IDs,
while the \code{MArrayLM} method extracts the gene lists automatically from a linear model fit object.
\code{goana} uses annotation from the appropriate Bioconductor organism package.
-\code{kegga} uses the KEGGREST package to access KEGG pathway annotation.
+\code{goana} can be used for any of the five species specified.
+
+\code{kegga} reads KEGG pathway annotation from the KEGG website.
+Note that the species name can be provided in either Bioconductor or KEGG format.
+\code{kegga} can be used for any species supported by KEGG, of which there are more than 14,000 possibilities.
+By default, \code{kegga} obtains the KEGG annotation for the specified species from the \url{http://rest.kegg.jp} website.
+Alternatively one can supply the required pathway annotation to \code{kegga} in the form of two data.frames.
+If this is done, then an internet connection is not required.
-The default method accepts a vector \code{prior.prob} giving the prior probability that each gene in the universe appears in a gene set.
+The ability to supply data.frame annotation to \code{kegga} means that \code{kegga} can in principle be used to analyze any user-supplied gene sets.
+
+The default \code{goana} and \code{kegga} methods accept a vector \code{prior.prob} giving the prior probability that each gene in the universe appears in a gene set.
This vector can be used to correct for unwanted trends in the differential expression analysis associated with gene length, gene abundance or any other covariate.
The \code{MArrayLM} object computes the \code{prior.prob} vector automatically when \code{trend} is non-\code{NULL}.
@@ -65,6 +87,13 @@ The only methodological difference is that \code{goana} and \code{kegga} compute
While \code{tricubeMovingAverage} does not enforce monotonicity, it has the advantage of numerical stability when \code{de} contains only a small number of genes.
}
+\note{
+\code{kegga} requires an internet connection unless \code{gene.pathway} and \code{pathway.names} are both supplied.
+
+The default for \code{kegga} with \code{species="Dm"} changed from \code{convert=TRUE} to \code{convert=FALSE} in limma 3.27.8.
+Users wanting to use Entrez Gene IDs for Drosophila should set \code{convert=TRUE}, otherwise fly-base IDs are assumed.
+}
+
\value{
The \code{goana} default method produces a data frame with a row for each GO term and the following columns:
\item{Term}{GO term.}
@@ -81,8 +110,8 @@ The \code{goana} method for \code{MArrayLM} objects produces a data frame with a
\item{N}{number of genes in the GO term.}
\item{Up}{number of up-regulated differentially expressed genes.}
\item{Down}{number of down-regulated differentially expressed genes.}
- \item{P.Up}{p-value for over-representation of GO term in up-regulated genes.}
- \item{P.Down}{p-value for over-representation of GO term in down-regulated genes.}
+ \item{P.Up}{p-value for over-representation of GO term in up-regulated genes. Not adjusted for multiple testing.}
+ \item{P.Down}{p-value for over-representation of GO term in down-regulated genes. Not adjusted for multiple testing.}
The row names of the data frame give the GO term IDs.
@@ -117,7 +146,7 @@ fit <- lmFit(y, design)
fit <- eBayes(fit)
# Standard GO analysis
-go.fisher <- goana(fit)
+go.fisher <- goana(fit, species="Hs")
topGO(go.fisher, sort = "up")
topGO(go.fisher, sort = "down")
@@ -140,6 +169,10 @@ topGO(go.de, sort = "DE2")
topGO(go.de, ontology = "BP", sort = "DE3")
topGO(go.de, ontology = "CC", sort = "DE3")
topGO(go.de, ontology = "MF", sort = "DE3")
+
+## Standard KEGG analysis
+k <- kegga(fit, species="Hs")
+k <- kegga(fit, species.KEGG="hsa") # equivalent to previous
}
}
diff --git a/man/ids2indices.Rd b/man/ids2indices.Rd
index 8ca3074..69ae0fd 100644
--- a/man/ids2indices.Rd
+++ b/man/ids2indices.Rd
@@ -31,10 +31,10 @@ There is a topic page on \link{10.GeneSetTests}.
\dontrun{
-download.file("http://bioinf.wehi.edu.au/software/MSigDB/human_c2_v4.rdata",
- "human_c2_v4.rdata", mode = "wb")
+download.file("http://bioinf.wehi.edu.au/software/MSigDB/human_c2_v5.rdata",
+ "human_c2_v5.rdata", mode = "wb")
-load("human_c2_v4.rdata")
+load("human_c2_v5.rdata")
c2.indices <- ids2indices(Hs.c2, y$genes$GeneID)
camera(y, c2.indices, design)
diff --git a/man/lmFit.Rd b/man/lmFit.Rd
index a3f0da2..fe2419c 100755
--- a/man/lmFit.Rd
+++ b/man/lmFit.Rd
@@ -125,7 +125,7 @@ dupcor <- duplicateCorrelation(y,design,ndups=2)
dupcor$consensus.correlation
fit4 <- lmFit(y,design,ndups=2,correlation=dupcor$consensus)
dim(fit4)
-fit4 <- eBayes(fit3)
+fit4 <- eBayes(fit4)
topTable(fit4,coef=2)
}
\keyword{models}
diff --git a/man/logcosh.Rd b/man/logcosh.Rd
new file mode 100644
index 0000000..d64cfaa
--- /dev/null
+++ b/man/logcosh.Rd
@@ -0,0 +1,26 @@
+\name{logcosh}
+\alias{logcosh}
+\title{Logarithm of cosh}
+\description{Compute \code{log(cosh(x))} without floating overflow or underflow}
+\usage{
+logcosh(x)
+}
+\arguments{
+ \item{x}{a numeric vector or matrix.}
+}
+
+\details{
+The computation uses asymptotic expressions for very large or very small arguments.
+For intermediate arguments, \code{log(cosh(x))} is returned.
+}
+
+\value{Numeric vector or matrix of same dimensions as \code{x}.}
+
+\author{Gordon K Smyth}
+
+\examples{
+x <- c(1e-8,1e-7,1e-6,1e-5,1e-4,1,3,50,800)
+logcosh(x)
+log(cosh(x))
+}
+
diff --git a/man/marraylm.Rd b/man/marraylm.Rd
index e74e542..44ec21c 100755
--- a/man/marraylm.Rd
+++ b/man/marraylm.Rd
@@ -34,6 +34,7 @@ Objects are normally created by \code{\link{lmFit}}, and additional components a
\tabular{ll}{
{\code{s2.prior}} \tab {numeric value giving empirical Bayes estimated prior value for residual variances}\cr
{\code{df.prior}} \tab {numeric vector giving empirical Bayes estimated degrees of freedom associated with \code{s2.prior} for each gene}\cr
+ {\code{df.total}} \tab {numeric vector giving total degrees of freedom used for each gene, usually equal to \code{df.prior + df.residual}.}\cr
{\code{s2.post}} \tab {numeric vector giving posterior residual variances}\cr
{\code{var.prior}} \tab {numeric vector giving empirical Bayes estimated prior variance for each true coefficient}\cr
{\code{F}} \tab {numeric vector giving moderated F-statistics for testing all contrasts equal to zero}\cr
diff --git a/man/nec.Rd b/man/nec.Rd
index 39ea504..6d8a614 100644
--- a/man/nec.Rd
+++ b/man/nec.Rd
@@ -30,6 +30,7 @@ After background correction, an \code{offset} is added to the data.
When expression values for negative controls are not available, the \code{detection.p} argument is used instead,
In that case, these functions call \code{\link{normexp.fit.detection.p}}, which infers the negative control probe intensities from the detection p-values associated with the regular probes.
+The function outputs a message if this is done.
For more detailed descriptions of the arguments \code{x}, \code{status}, \code{negctrl}, \code{regular} and \code{detection.p}, please refer to functions \code{\link{normexp.fit.control}}, \code{\link{normexp.fit.detection.p}} and \code{\link{read.ilmn}}.
diff --git a/man/plotRLDF.Rd b/man/plotRLDF.Rd
index a280d06..e5b095f 100644
--- a/man/plotRLDF.Rd
+++ b/man/plotRLDF.Rd
@@ -2,45 +2,69 @@
\name{plotRLDF}
\alias{plotRLDF}
\description{
-Plot of regularized linear discriminant functions for microarray data.
+Plot regularized linear discriminant functions for classifying samples based on expression data.
}
\usage{
-plotRLDF(y,design=NULL,z=NULL,labels.y=NULL,labels.z=NULL,col.y=1,col.z=1,
-df.prior=5,show.dimensions=c(1,2),main=NULL,nprobes=500,\dots)}
+plotRLDF(y, design=NULL, z=NULL, labels.y=NULL, labels.z=NULL,
+ col.y="black", col.z="black", show.dimensions=c(1,2),
+ ndim=max(show.dimensions), nprobes=100, plot=TRUE,
+ var.prior=NULL, df.prior=NULL, trend=FALSE, robust=FALSE, \dots)
+}
\arguments{
- \item{y}{any data object which can be coerced to a matrix, such as \code{ExpressionSet} or \code{EList}. The training dataset.}
- \item{z}{any data object which can be coerced to a matrix, such as \code{ExpressionSet} or \code{EList}. The dataset to be classified.}
- \item{design}{the design matrix ofthe microarray experiment for \code{y}, with rows corresponding to arrays and columns to coefficients to be estimated. Defaults to the unit vector meaning that the arrays are treated as replicates.}
- \item{labels.y}{character vector of sample names or labels in \code{y}. Default is integers starting from 1.}
- \item{labels.z}{character vector of sample names or labels in \code{z}. Default is \code{letters}.}
- \item{col.y}{numeric or character vector of colors for the plotting characters of \code{y}. Default is black.}
- \item{col.z}{numeric or character vector of colors for the plotting characters of \code{z}. Default is black.}
- \item{df.prior}{prior degrees of freedom for residual variances. Used in gene selection.}
- \item{show.dimensions}{which two dimensions should be plotted, numeric vector of length two.}
- \item{main}{title of the plot.}
- \item{nprobes}{number of probes to be used for the calculations. Selected by moderated F tests.}
+ \item{y}{the training dataset. Can be any data object which can be coerced to a matrix, such as \code{ExpressionSet} or \code{EList}.}
+ \item{design}{design matrix defining the training groups to be distinguished. The first column is assumed to represent the intercept.
+ Defaults to \code{model.matrix(~factor(labels.y))}.}
+ \item{z}{the dataset to be classified. Can be any data object which can be coerced to a matrix, such as \code{ExpressionSet} or \code{EList}. Rows must correspond to rows of \code{y}.}
+ \item{labels.y}{character vector of sample names or labels in \code{y}. Defaults to \code{colnames(y)} or failing that to \code{1:n}.}
+ \item{labels.z}{character vector of sample names or labels in \code{z}. Defaults to \code{colnames(z)} or failing that to \code{letters[1:n]}.}
+ \item{col.y}{colors for the plotting \code{labels.y}.}
+ \item{col.z}{colors for the plotting \code{labels.z}.}
+ \item{show.dimensions}{integer vector of length two indicating which two discriminant functions to plot. Functions are in decreasing order of discriminatory power.}
+ \item{ndim}{number of discriminant functions to compute}
+ \item{nprobes}{number of probes to be used for the calculations. The probes will be selected by moderated F statistic.}
+ \item{plot}{logical, should a plot be created?}
+ \item{var.prior}{prior variances, for regularizing the within-group covariance matrix. By default is estimated by \code{squeezeVar}.}
+ \item{df.prior}{prior degrees of freedom for regularizing the within-group covariance matrix. By default is estimated by \code{squeezeVar}.}
+ \item{trend}{logical, should a trend be estimated for \code{var.prior}? See \code{eBayes} for details. Only used if \code{var.prior} or \code{df.prior} are \code{NULL}.}
+ \item{robust}{logical, should \code{var.prior} and \code{df.prior} be estimated robustly? See \code{eBayes} for details. Only used if \code{var.prior} or \code{df.prior} are \code{NULL}.}
\item{\dots}{any other arguments are passed to \code{plot}.}
}
+
\details{
-This function is a variation on the plot of usual linear discriminant fuction, in that the within-group covariance matrix is regularized to ensure that it is invertible, with eigenvalues bounded away from zero.
-A diagonal regulation using \code{df.prior} and the median within-group variance is used.
+The function builds discriminant functions from the training data (\code{y}) and applies them to the test data (\code{z}).
+The method is a variation on classifical linear discriminant functions (LDFs), in that the within-group covariance matrix is regularized to ensure that it is invertible, with eigenvalues bounded away from zero.
+The within-group covariance matrix is squeezed towards a diagonal matrix with empirical Bayes posterior variances as diagonal elements.
The calculations are based on a filtered list of probes.
The \code{nprobes} probes with largest moderated F statistics are used to discriminate.
-See \code{\link[graphics]{text}} for possible values for \code{col} and \code{cex}.
+The \code{ndim} argument allows all required LDFs to be computed even though only two are plotted.
}
-\value{A list containing metagene information is (invisibly) returned.
-A plot is created on the current graphics device.}
-\author{Di Wu and Gordon Smyth}
+\note{
+Default values for \code{df.prior} and \code{var.prior} were changed in limma 3.27.10.
+}
+
+\value{
+If \code{plot=TRUE} a plot is created on the current graphics device.
+A list containing the following components is (invisibly) returned:
+ \item{training}{numeric matrix with \code{ncol(y)} rows and \code{ndim} columns containing discriminant functions evaluated for the training data.}
+ \item{predicting}{numeric matrix with \code{ncol(z)} rows and \code{ndim} columns containing discriminant functions evalulated on the classification data.}
+ \item{top}{integer vector of length \code{nprobes} giving indices of probes used.}
+ \item{metagenes}{numeric matrix with \code{nprobes} rows and \code{ndim} columns containing probe weights defining each discriminant function.}
+ \item{singular.values}{singular.values showing the predictive power of each discriminant function.}
+ \item{rank}{maximum number of discriminant functions with singular.values greater than zero.}
+ \item{var.prior}{numeric vector of prior variances.}
+ \item{df.prior}{numeric vector of prior degrees of freedom.}
+}
+
+\author{Gordon Smyth, Di Wu and Yifang Hu}
\seealso{
\code{lda} in package \code{MASS}
}
\examples{
-
# Simulate gene expression data for 1000 probes and 6 microarrays.
# Samples are in two groups
# First 50 probes are differentially expressed in second group
@@ -57,7 +81,9 @@ z[1:50,1:3] <- z[1:50,1:3] - 0.2
design <- cbind(Grp1=1,Grp2vs1=c(0,0,0,1,1,1))
options(digit=3)
-plotRLDF(y,z, design=design)
+# Samples 1-6 are training set, samples a-f are test set:
+plotRLDF(y, design, z=z, col.y="black", col.z="red")
+legend("top", pch=16, col=c("black","red"), legend=c("Training","Predicted"))
}
\keyword{hplot}
diff --git a/man/propexpr.Rd b/man/propexpr.Rd
index ad28e6c..f3ce8b6 100644
--- a/man/propexpr.Rd
+++ b/man/propexpr.Rd
@@ -8,27 +8,34 @@ propexpr(x, neg.x=NULL, status=x$genes$Status, labels=c("negative","regular"))
\arguments{
\item{x}{matrix or similar object containing raw intensities for a set of arrays.}
\item{neg.x}{matrix or similar object containing raw intensities for negative control probes for the same arrays. If \code{NULL}, then negative controls must be provided in \code{x}.}
- \item{status}{character vector giving probe types.}
- \item{labels}{character vector giving probe type identifiers.}
+ \item{status}{character vector specifying control type of each probe. Only used if \code{neg.x} is \code{NULL}.}
+ \item{labels}{character vector giving the \code{status} values for negative control probes and regular (non-control) probes respectively. If of length 1, then all probes other than the negative controls are assumed to be regular. Only used if \code{neg.x} is \code{NULL}.}
}
+
\details{
-This function estimates the proportion of expressed in a microarray by utilizing the negative control probes.
+This function estimates the overall proportion of probes on each microarray that are correspond to expressed genes using the method of Shi et al (2010).
+The function is especially useful for Illumina BeadChips arrays, although it can in principle be applied to any platform with good quality negative controls.
+
+The negative controls can be supplied either as rows of \code{x} or as a separate matrix.
+If supplied as rows of \code{x}, then the negative controls are identified by the \code{status} vector.
+\code{x} might also include other types of control probes, but these will be ignored in the calculation.
+
Illumina BeadChip arrays contain 750~1600 negative control probes.
-The expression profile of these control probes can be saved to a separate file by the Illumina BeadStudio software when using it to output the expression profile for regular probes.
-The control probe profile could be re-generated if it was not generated when the regular probe profile was created by BeadStudio.
-Other microarray platforms can also use this function to estimate the proportion of expressed probes in each array, provided that they have a set of negative control probes.
-
-\code{labels} can include one or two probe type identifiers.
-Its first element should be the identifier for negative control probes (\code{negative} by default).
-If \code{labels} only contains one identifier, then it will be assumed to contain the identifier for negative control probes.
-By default, \code{regular} is the identifier for regular probes.
+If \code{read.idat} is used to read Illumina expression IDAT files, then the control probes will be populated as rows of the output \code{EListRaw} object, and the vector \code{x$genes$Status} will be set to identify control probes.
+
+Alternatively, expression values can be exported from Illumina's GenomeStudio software as tab-delimited text files.
+In this case, the control probes are usually written to a separate file from the regular probes.
}
+
\value{
Numeric vector giving the proportions of expressed probes in each array.
}
\references{
-Shi, W, de Graaf, C, Kinkel, S, Achtman, A, Baldwin, T, Schofield, L, Scott, H, Hilton, D, Smyth, GK (2010). Estimating the proportion of microarray probes expressed in an RNA sample. Nucleic Acids Research 38, 2168-2176.
+Shi, W, de Graaf, C, Kinkel, S, Achtman, A, Baldwin, T, Schofield, L, Scott, H, Hilton, D, Smyth, GK (2010).
+Estimating the proportion of microarray probes expressed in an RNA sample.
+\emph{Nucleic Acids Research} 38, 2168-2176.
+\url{http://nar.oxfordjournals.org/content/38/7/2168}
}
\author{Wei Shi and Gordon Smyth}
@@ -39,9 +46,16 @@ Description to the control probes in Illumina BeadChips can be found in \code{\l
\examples{
\dontrun{
-x <- read.ilmn(files="sample probe profile.txt",ctrlfiles="control probe profile.txt")
-propexpr(x, )
+# Read Illumina binary IDAT files
+x <- read.idat(idat, bgx)
+propexpr(x)
+
+# Read text files exported from GenomeStudio
+x <- read.ilmn(files = "sample probe profile.txt",
+ ctrlfiles = "control probe profile.txt")
+propexpr(x)
}
}
-\keyword{models}
\ No newline at end of file
+\keyword{background correction}
+\keyword{illumina beadchips}
diff --git a/man/read.idat.Rd b/man/read.idat.Rd
index a2e247d..39f66a7 100644
--- a/man/read.idat.Rd
+++ b/man/read.idat.Rd
@@ -6,64 +6,93 @@
\description{Read Illumina BeadArray data from IDAT and manifest (.bgx) files for gene expression platforms.}
\usage{
-read.idat(idatfiles, bgxfile, dateinfo=FALSE, tolerance=0)
+read.idat(idatfiles, bgxfile, dateinfo = FALSE, annotation = "Symbol",
+ tolerance = 0L, verbose = TRUE)
}
\arguments{
\item{idatfiles}{character vector specifying idat files to be read in.}
\item{bgxfile}{character string specifying bead manifest file (.bgx) to be read in.}
- \item{dateinfo}{logical. Should date and software version info be read in?}
- \item{tolerance}{numeric. The number of probe ID discrepancies allowed between the manifest and a given idat file.}
+ \item{dateinfo}{logical. Should date and software version information be read in?}
+ \item{annotation}{character vector of annotation columns to be read from the manifest file.}
+ \item{tolerance}{integer. The number of probe ID discrepancies allowed between the manifest and any of the IDAT files.}
+ \item{verbose}{logical. Should progress messages are sent to standard output?}
}
\details{
- Illumina's BeadScan/iScan software ouputs probe intensities in IDAT
- format (encrypted XML files) and probe info in a platform specific manifest file (.bgx).
+ Illumina's BeadScan/iScan software outputs probe intensities in IDAT
+ format (encrypted XML files) and uses probe information stored in a platform specific manifest file (.bgx).
These files can be processed using the low-level functions \code{readIDAT} and \code{readBGX}
- from the \code{illuminaio} package (Smith et al. 2013).
-
- The \code{read.idat} function provides a convenient way to read these files.
- into R and store them in an \code{EListRaw-class} object (similar to \code{read.ilmn},
- which imports data output by Illumina's GenomeStudio software) that can be used
- by downstream processing functions in \code{limma}.
+ from the \code{illuminaio} package (Smith et al. 2013).
+ The \code{read.idat} function provides a convenient way to read these files
+ into R and to store them in an \code{EListRaw-class} object.
+ The function serves a similar purpose to \code{\link{read.ilmn}},
+ which reads text files exported by Illumina's GenomeStudio software,
+ but it reads the IDAT files directly without any need to convert them first to text.
+
+ The function reads information on control probes as well for regular probes.
Probe types are indicated in the \code{Status} column of the \code{genes}
- component of the \code{EListRaw-class} object.
+ component of the \code{EListRaw} object.
+
+ The \code{annotation} argument specifies probe annotation columns to be extracted from the manifest file.
+ The manifest typically contains the following columns:
+ \code{"Species"}, \code{"Source"}, \code{"Search_Key"}, \code{"Transcript"},
+ \code{"ILMN_Gene"}, \code{"Source_Reference_ID"}, \code{"RefSeq_ID"},
+ \code{"Unigene_ID"}, \code{"Entrez_Gene_ID"}, \code{"GI"},
+ \code{"Accession"}, \code{"Symbol"}, \code{"Protein_Product"},
+ \code{"Probe_Id"}, \code{"Array_Address_Id"}, \code{"Probe_Type"},
+ \code{"Probe_Start"}, \code{"Probe_Sequence"}, \code{"Chromosome"},
+ \code{"Probe_Chr_Orientation"}, \code{"Probe_Coordinates"}, \code{"Cytoband"},
+ \code{"Definition"}, \code{"Ontology_Component"}, \code{"Ontology_Process"},
+ \code{"Ontology_Function"}, \code{"Synonyms"}, \code{"Obsolete_Probe_Id"}.
+ Note that \code{"Probe_Id"} and \code{"Array_Address_Id"} are always extracted and
+ do not need to included in the \code{annotation} argument.
+
+ If more than \code{tolerance} probes in the manifest cannot be found in an IDAT file then the function will return an error.
}
\value{
- An \code{EListRaw-class} object with the following components:
+ An \code{EListRaw} object with the following components:
\item{E}{ numeric matrix of raw intensities.}
- \item{other}{ list containing matrices of \code{NumBeads} and \code{STDEV} for each probe.}
- \item{genes}{ data.frame of probe annotation.}
- \item{targets}{ data.frame of sample information.}
+ \item{other$NumBeads}{ numeric matrix of same dimensions as \code{E} giving number of beads used for each intensity value.}
+ \item{other$STDEV}{ numeric matrix of same dimensions as \code{E} giving bead-level standard deviation or standard error for each intensity value.}
+ \item{genes}{ data.frame of probe annotation.
+ This includes the \code{Probe_Id} and \code{Array_Address_Id} columns extracted from the manifest file,
+ plus a \code{Status} column identifying control probes,
+ plus any other columns specified by \code{annotation}.}
+ \item{targets}{ data.frame of sample information.
+ This includes the IDAT file names plus other columns if \code{dateinfo=TRUE}.}
}
\references{
Smith ML, Baggerly KA, Bengtsson H, Ritchie ME, Hansen KD (2013).
illuminaio: An open source IDAT parsing tool. \emph{F1000 Research} 2, 264.
-\url{http://f1000research.com/articles/2-264/v1}
+\url{http://f1000research.com/articles/2-264/}
}
\author{Matt Ritchie}
\seealso{
- \code{read.ilmn} imports gene expression data output by GenomeStudio.
+ \code{\link{read.ilmn}} imports gene expression data output by GenomeStudio.
- \code{neqc} performs normexp by control background correction, log
+ \code{\link{neqc}} performs normexp by control background correction, log
transformation and quantile between-array normalization for
Illumina expression data.
- \code{propexpr} estimates the proportion of expressed probes in a microarray.
+ \code{\link{propexpr}} estimates the proportion of expressed probes in a microarray.
+
+ \code{\link{detectionPValues}} computes detection p-values from the negative controls.
}
\examples{
\dontrun{
-idatfiles = dir(pattern="idat")
-bgxfile = dir(pattern="bgx")
-data = read.idat(idatfiles, bgxfile)
+idatfiles <- dir(pattern="idat")
+bgxfile <- dir(pattern="bgx")
+x <- read.idat(idatfiles, bgxfile)
+x$other$Detection <- detectionPValues(x)
propexpr(data)
-datanorm = neqc(data)
+y <- neqc(data)
}
}
diff --git a/man/read.maimages.Rd b/man/read.maimages.Rd
index 5ff79ed..c9de379 100755
--- a/man/read.maimages.Rd
+++ b/man/read.maimages.Rd
@@ -21,7 +21,8 @@ read.imagene(files, path=NULL, ext=NULL, names=NULL, columns=NULL, other.columns
\item{path}{character string giving the directory containing the files.
The default is the current working directory.}
\item{ext}{character string giving optional extension to be added to each file name}
- \item{names}{character vector of names to be associated with each array as column name.
+ \item{names}{character vector of unique names to be associated with each array as column name.
+ Can be supplied as \code{files$Label} if \code{files} is a data.frame.
Defaults to \code{removeExt(files)}.}
\item{columns}{list, or named character vector.
For two color data, this should have fields \code{R}, \code{G}, \code{Rb} and \code{Gb} giving the column names to be used for red and green foreground and background or, in the case of Imagene data, a list with fields \code{f} and \code{b}.
diff --git a/man/removeext.Rd b/man/removeext.Rd
index 6c9a410..d4f2add 100755
--- a/man/removeext.Rd
+++ b/man/removeext.Rd
@@ -5,11 +5,19 @@
\description{Finds and removes any common extension from a vector of file names.}
\usage{
-removeExt(x)
+removeExt(x, sep=".")
}
\arguments{
\item{x}{character vector}
+\item{sep}{character string that separates the body of each character string from the extension.}
+}
+
+\details{
+This function is used for simplifying file names, or any vector of character strings, when the strings all finish with the same suffix or extension.
+If the same extension is not shared by every element of \code{x}, then it is not removed from any element.
+
+Note that \code{sep} is interpreted as a literal character string: it is not a regular expression.
}
\value{
@@ -23,6 +31,9 @@ An overview of LIMMA functions for reading data is given in \link{03.ReadingData
\examples{
x <- c("slide1.spot","slide2.spot","slide3.spot")
removeExt(x)
+
+x <- c("Harry - a name from Harry Potter","Hermione - a name from Harry Potter")
+removeExt(x, sep=" - ")
}
\author{Gordon Smyth}
diff --git a/man/roast.Rd b/man/roast.Rd
index 374afc6..fd336dc 100644
--- a/man/roast.Rd
+++ b/man/roast.Rd
@@ -13,42 +13,44 @@ Rotation gene set testing for linear models.
}
\usage{
-\S3method{fry}{default}(y, index = NULL, design = NULL, contrast = ncol(design),
- weights = NULL, sort = TRUE, \dots)
-\S3method{roast}{default}(y, index = NULL, design = NULL, contrast = ncol(design),
- set.statistic = "mean", gene.weights = NULL, array.weights = NULL, weights = NULL,
- block = NULL, correlation, var.prior = NULL, df.prior = NULL, trend.var = FALSE,
+\method{roast}{default}(y, index = NULL, design = NULL, contrast = ncol(design), geneid = NULL,
+ set.statistic = "mean", gene.weights = NULL, var.prior = NULL, df.prior = NULL,
nrot = 999, approx.zscore = TRUE, \dots)
-\S3method{mroast}{default}(y, index = NULL, design = NULL, contrast = ncol(design),
- set.statistic = "mean", gene.weights = NULL, array.weights = NULL, weights = NULL,
- block = NULL, correlation, var.prior = NULL, df.prior = NULL, trend.var = FALSE,
- nrot = 999, approx.zscore = TRUE, adjust.method = "BH", midp = TRUE,
- sort = "directional", \dots)
+\method{mroast}{default}(y, index = NULL, design = NULL, contrast = ncol(design), geneid = NULL,
+ set.statistic = "mean", gene.weights = NULL, var.prior = NULL, df.prior = NULL,
+ nrot = 999, approx.zscore = TRUE, adjust.method = "BH",
+ midp = TRUE, sort = "directional", \dots)
+\method{fry}{default}(y, index = NULL, design = NULL, contrast = ncol(design), geneid = NULL,
+ standardize = "posterior.sd", sort = "directional", \dots)
}
\arguments{
\item{y}{numeric matrix giving log-expression or log-ratio values for a series of microarrays, or any object that can coerced to a matrix including \code{ExpressionSet}, \code{MAList}, \code{EList} or \code{PLMSet} objects.
- Rows correspond to probes and columns to samples.
- If either \code{var.prior} or \code{df.prior} are null, then \code{y} should contain values for all genes on the arrays. If both prior parameters are given, then only \code{y} values for the test set are required.}
- \item{index}{index vector specifying which rows (probes) of \code{y} are in the test set. This can be a vector of indices, or a logical vector of the same length as \code{statistics}, or any vector such as \code{y[index,]} contains the values for the gene set to be tested.
- For \code{mroast}, \code{index} is a list of index vectors. The list can be made using \link{ids2indices}.}
+ Rows correspond to probes and columns to samples.
+ If either \code{var.prior} or \code{df.prior} are \code{NULL}, then \code{y} should contain values for all genes on the arrays. If both prior parameters are given, then only \code{y} values for the test set are required.}
+ \item{index}{index vector specifying which rows (probes) of \code{y} are in the test set.
+ Can be a vector of integer indices, or a logical vector of length \code{nrow(y)}, or a vector of gene IDs corresponding to entries in \code{geneid}.
+ Alternatively it can be a data.frame with the first column containing the index vector and the second column containing gene weights.
+ For \code{mroast} or \code{fry}, \code{index} is a list of index vectors or a list of data.frames. }
\item{design}{design matrix}
- \item{contrast}{contrast for which the test is required. Can be an integer specifying a column of \code{design}, or the name of a column of \code{design}, or a numeric contrast vector of length equal to the number of columns of \code{design}.}
+ \item{contrast}{contrast for which the test is required.
+ Can be an integer specifying a column of \code{design}, or the name of a column of \code{design}, or a numeric contrast vector of length equal to the number of columns of \code{design}.}
+ \item{geneid}{gene identifiers corresponding to the rows of \code{y}.
+ Can be either a vector of length \code{nrow(y)} or the name of the column of \code{y$genes} containing the gene identifiers.
+ Defaults to \code{rownames(y)}.}
\item{set.statistic}{summary set statistic. Possibilities are \code{"mean"},\code{"floormean"},\code{"mean50"} or \code{"msq"}.}
- \item{gene.weights}{optional numeric vector of weights for genes in the set. Can be positive or negative. For \code{mroast} this vector must have length equal to \code{nrow(y)}. For \code{roast}, can be of length \code{nrow(y)} or of length equal to the number of genes in the test set.}
- \item{array.weights}{optional numeric vector of array weights.}
- \item{weights}{optional matrix of observation weights. If supplied, should be of same dimensions as \code{y} and all values should be positive. If \code{y} is an \code{EList} or \code{MAList} object containing weights, then those weights will be used.}
- \item{block}{optional vector of blocks.}
- \item{correlation}{correlation between blocks.}
+ \item{gene.weights}{numeric vector of (positive or negative) probewise weights.
+ For \code{mroast} or \code{fry}, this vector must have length equal to \code{nrow(y)}.
+ For \code{roast}, can be of length \code{nrow(y)} or of length equal to the number of genes in the test set.}
\item{var.prior}{prior value for residual variances. If not provided, this is estimated from all the data using \code{squeezeVar}.}
\item{df.prior}{prior degrees of freedom for residual variances. If not provided, this is estimated using \code{squeezeVar}.}
- \item{trend.var}{logical, should a trend be estimated for \code{var.prior}? See \code{eBayes} for details. Only used if \code{var.prior} or \code{df.prior} are \code{NULL}.}
- \item{nrot}{number of rotations used to estimate the p-values.}
+ \item{nrot}{number of rotations used to compute the p-values.}
+ \item{approx.zscore}{logical, if \code{TRUE} then a fast approximation is used to convert t-statistics into z-scores prior to computing set statistics. If \code{FALSE}, z-scores will be exact.}
\item{adjust.method}{method used to adjust the p-values for multiple testing. See \code{\link{p.adjust}} for possible values.}
\item{midp}{logical, should mid-p-values be used in instead of ordinary p-values when adjusting for multiple testing?}
- \item{sort}{character, whether to sort output table by directional p-value (\code{"directional"}), non-directional p-value (\code{"mixed"}), or not at all (\code{"none"}). For \code{fry} it is logical, whether to sort or not.}
- \item{approx.zscore}{logical, if \code{TRUE} then a fast approximation is used to convert t-statistics into z-scores prior to computing set statistics. If \code{FALSE}, z-scores will be exact.}
- \item{\dots}{other arguments not currently used.}
+ \item{sort}{character, whether to sort output table by directional p-value (\code{"directional"}), non-directional p-value (\code{"mixed"}), or not at all (\code{"none"}).}
+ \item{standardize}{how to standardize for unequal probewise variances. Possibilities are \code{"residual.sd"}, \code{"posterior.sd"} or \code{"none"}.}
+ \item{\dots}{any argument that would be suitable for \code{\link{lmFit}} or \code{\link{eBayes}} can be included.}
}
\value{
@@ -84,6 +86,8 @@ The design matrix for the experiment is specified as for the \code{\link{lmFit}}
This allows users to focus on differential expression for any coefficient or contrast in a linear model.
If \code{contrast} is not specified, then the last coefficient in the linear model will be tested.
+The argument \code{index} is often made using \link{ids2indices}.
+
The argument \code{gene.weights} allows directional weights to be set for individual genes in the set.
This is often useful, because it allows each gene to be flagged as to its direction and magnitude of change based on prior experimentation.
A typical use is to make the \code{gene.weights} \code{1} or \code{-1} depending on whether the gene is up or down-regulated in the pathway under consideration.
@@ -112,8 +116,10 @@ This means that the smallest possible p-value is \code{1/(nrot+1)}.
By default, \code{mroast} reports ordinary p-values but uses mid-p-values (Routledge, 1994) at the multiple testing stage.
Mid-p-values are probably a good choice when using false discovery rates (\code{adjust.method="BH"}) but not when controlling the family-wise type I error rate (\code{adjust.method="holm"}).
-\code{fry} is a fast version of \code{mroast} in the special case that \code{df.prior} is large and \code{set.statistic="mean"}.
-In this situation, it gives the same result as \code{mroast} with an infinite number of rotations.
+\code{fry} is a fast approximation to \code{mroast}.
+In the special case that \code{df.prior} is large and \code{set.statistic="mean"}, \code{fry} gives the same result as \code{mroast} with an infinite number of rotations.
+In other circumstances, when genes have different variances, \code{fry} uses a standardization strategy to approximate the \code{mroast} results.
+Using \code{fry} may be advisable when performing tests for a large number of sets, because it is fast and because the \code{fry} p-values are not limited by the number of rotations performed.
}
\note{
@@ -160,8 +166,7 @@ y[index1,3:4] <- y[index1,3:4]+3
index2 <- 6:10
roast(y,index1,design,contrast=2)
-mroast(y,index1,design,contrast=2)
-mroast(y,list(set1=index1,set2=index2),design,contrast=2)
+fry(y,list(set1=index1,set2=index2),design,contrast=2)
}
\keyword{gene set test}
diff --git a/man/romer.Rd b/man/romer.Rd
index 1c85c97..000ad2f 100644
--- a/man/romer.Rd
+++ b/man/romer.Rd
@@ -7,8 +7,9 @@
Gene set enrichment analysis for linear models using rotation tests (ROtation testing using MEan Ranks).
}
\usage{
-\S3method{romer}{default}(y, index, design, contrast = ncol(design), array.weights = NULL, block = NULL,
- correlation, set.statistic = "mean", nrot = 9999, shrink.resid = TRUE, \dots)
+\method{romer}{default}(y, index, design = NULL, contrast = ncol(design),
+ array.weights = NULL, block = NULL, correlation,
+ set.statistic = "mean", nrot = 9999, shrink.resid = TRUE, \dots)
}
\arguments{
\item{y}{numeric matrix giving log-expression values.}
diff --git a/man/squeezeVar.Rd b/man/squeezeVar.Rd
index 6a2c7b5..9203938 100755
--- a/man/squeezeVar.Rd
+++ b/man/squeezeVar.Rd
@@ -35,7 +35,7 @@ The squeezed variances have a smaller expected mean square error to the true var
If \code{covariate} is non-null, then the scale parameter of the prior distribution is assumed to depend on the covariate.
If the covariate is average log-expression, then the effect is an intensity-dependent trend similar to that in Sartor et al (2006).
-\code{robust=TRUE} implements the robust empirical Bayes procedure of Phipson et al (2013) which allows some of the \code{var} values to be outliers.
+\code{robust=TRUE} implements the robust empirical Bayes procedure of Phipson et al (2016) which allows some of the \code{var} values to be outliers.
}
\note{
@@ -53,10 +53,10 @@ A list with components
\author{Gordon Smyth}
\references{
-Phipson, B, Lee, S, Majewski, IJ, Alexander, WS, and Smyth, GK (2013).
-Empirical Bayes in the presence of exceptional cases, with application to microarray data.
-Bioinformatics Division, Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia.
-\url{http://www.statsci.org/smyth/pubs/RobustEBayesPreprint.pdf}
+Phipson, B, Lee, S, Majewski, IJ, Alexander, WS, and Smyth, GK (2016).
+Robust hyperparameter estimation protects against hypervariable genes and improves power to detect differential expression.
+\emph{Annals of Applied Statistics} 10.
+\url{http://arxiv.org/abs/1602.08678}
Sartor MA, Tomlinson CR, Wesselkamper SC, Sivaganesan S, Leikauf GD, Medvedovic M (2006).
Intensity-based hierarchical Bayes method improves testing for differentially expressed genes in microarray experiments.
diff --git a/man/targetsA2C.Rd b/man/targetsA2C.Rd
index e14a682..602f3bd 100755
--- a/man/targetsA2C.Rd
+++ b/man/targetsA2C.Rd
@@ -29,6 +29,12 @@ The second format may be more appropriate if the data is to be analyzed in terms
data.frame with twice as many rows as \code{targets}.
Any pair of columns named by \code{channel.columns} will now be one column.
}
+\references{
+Smyth, GK, and Altman, NS (2013).
+Separate-channel analysis of two-channel microarrays: recovering inter-spot information.
+\emph{BMC Bioinformatics} 14, 165.
+\url{http://www.biomedcentral.com/1471-2105/14/165}
+}
\seealso{
\code{targetsA2C} is used by the \code{coerce} method from \code{RGList} to \code{ExpressionSet} in the convert package.
diff --git a/man/venn.Rd b/man/venn.Rd
index 7f8d6db..bf9c8c5 100755
--- a/man/venn.Rd
+++ b/man/venn.Rd
@@ -3,7 +3,7 @@
\alias{vennDiagram}
\title{Venn Diagrams}
\description{
-Compute Classification Counts and Make Venn Diagram.
+Compute classification counts and draw a Venn diagram.
}
\usage{
vennCounts(x, include="both")
@@ -11,9 +11,11 @@ vennDiagram(object, include="both", names=NULL, mar=rep(1,4), cex=c(1.5,1,0.7),
circle.col=NULL, counts.col=NULL, show.include=NULL, \dots)
}
\arguments{
- \item{x}{numeric matrix of 0's and 1's indicating significance of a test.
+ \item{x}{a \code{TestResults} matrix.
+ This is numeric matrix of 0's, 1's and -1's indicating significance of a test or membership of a set.
+ Each row corresponds to a gene and each column to a contrast or set.
Usually created by \code{\link{decideTests}}.}
- \item{object}{either a \code{TestResults} matrix from \code{decideTests} or a \code{VennCounts} object produced by \code{vennCounts}.}
+ \item{object}{either a \code{TestResults} matrix or a \code{VennCounts} object produced by \code{vennCounts}.}
\item{include}{character vector specifying whether all differentially expressed genes should be counted, or whether the counts should be restricted to genes changing in a certain direction. Choices are \code{"both"} for all differentially expressed genes, \code{"up"} for up-regulated genes only or \code{"down"} for down-regulated genes only. If \code{include=c("up","down")} then both the up and down counts will be shown. This argument is ignored if \code{object} if \code{object} is al [...]
\item{names}{character vector giving names for the sets or contrasts}
\item{mar}{numeric vector of length 4 specifying the width of the margins around the plot. This argument is passed to \code{par}.}
@@ -26,13 +28,21 @@ vennDiagram(object, include="both", names=NULL, mar=rep(1,4), cex=c(1.5,1,0.7),
}
\details{
-\code{vennDiagram} can plot up to five sets.
+Each column of \code{x} corresponds to a contrast or set, and the entries of \code{x} indicate membership of each row in each set or alternatively the significance of each row for each contrast.
+In the latter case, the entries can be negative as well as positive to indicate the direction of change.
+
\code{vennCounts} can collate intersection counts for any number of sets.
+\code{vennDiagram} can plot up to five sets.
}
\value{
-\code{vennCounts} produces a \code{VennCounts} object, which is a numeric matrix with last column \code{"Counts"} giving counts for each possible vector outcome.
-\code{vennDiagram} causes a plot to be produced on the current graphical device.
+\code{vennCounts} produces an object of class \code{"VennCounts"}.
+This contains only one slot, which is numerical matrix with \code{2^ncol{x}} rows and \code{ncol(x)+1} columns.
+Each row corresponds to a particular combination of set memberships.
+The first \code{ncol{x}} columns of output contain 1 or 0 indicating membership or not in each set.
+The last column called \code{"Counts"} gives the number of rows of \code{x} corresponding to that combination of memberships.
+
+\code{vennDiagram} produces no output but causes a plot to be produced on the current graphical device.
}
\seealso{
diff --git a/man/voom.Rd b/man/voom.Rd
index af41d09..17e3883 100644
--- a/man/voom.Rd
+++ b/man/voom.Rd
@@ -8,7 +8,7 @@ The data are then ready for linear modelling.
\usage{
voom(counts, design = NULL, lib.size = NULL, normalize.method = "none",
- plot = FALSE, span=0.5, \dots)
+ span = 0.5, plot = FALSE, save.plot = FALSE, \dots)
}
\arguments{
\item{counts}{a numeric \code{matrix} containing raw counts, or an \code{ExpressionSet} containing raw counts, or a \code{DGEList} object.}
@@ -18,8 +18,9 @@ voom(counts, design = NULL, lib.size = NULL, normalize.method = "none",
Otherwise library sizes are calculated from the columnwise counts totals.}
\item{normalize.method}{normalization method to be applied to the logCPM values.
Choices are as for the \code{method} argument of \code{normalizeBetweenArrays} when the data is single-channel.}
- \item{plot}{\code{logical}, should a plot of the mean-variance trend be displayed?}
\item{span}{width of the lowess smoothing window as a proportion.}
+ \item{plot}{logical, should a plot of the mean-variance trend be displayed?}
+ \item{save.plot}{logical, should the coordinates and line of the plot be saved in the output?}
\item{\dots}{other arguments are passed to \code{lmFit}.}
}
@@ -44,6 +45,8 @@ An \code{\link[limma:EList]{EList}} object with the following components:
\item{design}{design matrix}
\item{lib.size}{numeric vector of total normalized library sizes}
\item{genes}{dataframe of gene annotation extracted from \code{counts}}
+\item{voom.xy}{if \code{save.plot}, list containing x and y coordinates for points in mean-variance plot}
+\item{voom.line}{if \code{save.plot}, list containing coordinates of loess line in the mean-variance plot}
}
\author{Charity Law and Gordon Smyth}
@@ -56,11 +59,10 @@ Voom: precision weights unlock linear model analysis tools for RNA-seq read coun
}
\seealso{
-\code{\link{voomWithQualityWeights}}
+\code{\link{voomWithQualityWeights}}.
+\code{\link{vooma}} is similar to \code{voom} but for microarrays instead of RNA-seq.
A summary of functions for RNA-seq analysis is given in \link{11.RNAseq}.
-
-\code{\link{vooma}} is similar to \code{voom} but for microarrays instead of RNA-seq.
}
\keyword{rna-seq}
diff --git a/tests/limma-Tests.R b/tests/limma-Tests.R
index 1bcf9d1..974fff2 100755
--- a/tests/limma-Tests.R
+++ b/tests/limma-Tests.R
@@ -231,9 +231,12 @@ mroast(y=y,list(set1=iset1,set2=iset2),design,contrast=2,gene.weights=runif(100)
mroast(y=y,list(set1=iset1,set2=iset2),design,contrast=2,array.weights=c(0.5,1,0.5,1))
mroast(y=y,list(set1=iset1,set2=iset2),design,contrast=2,weights=w)
mroast(y=y,list(set1=iset1,set2=iset2),design,contrast=2,weights=w,array.weights=c(0.5,1,0.5,1))
+fry(y=y,list(set1=iset1,set2=iset2),design,contrast=2,weights=w,array.weights=c(0.5,1,0.5,1))
### camera
+camera(y=y,iset1,design,contrast=2,weights=c(0.5,1,0.5,1),allow.neg.cor=TRUE,inter.gene.cor=NA)
+camera(y=y,list(set1=iset1,set2=iset2),design,contrast=2,allow.neg.cor=TRUE,inter.gene.cor=NA)
camera(y=y,iset1,design,contrast=2,weights=c(0.5,1,0.5,1))
camera(y=y,list(set1=iset1,set2=iset2),design,contrast=2)
@@ -241,6 +244,7 @@ camera(y=y,list(set1=iset1,set2=iset2),design,contrast=2)
y <- new("EList",list(E=y))
roast(y=y,iset1,design,contrast=2)
+camera(y=y,iset1,design,contrast=2,allow.neg.cor=TRUE,inter.gene.cor=NA)
camera(y=y,iset1,design,contrast=2)
### eBayes with trend
diff --git a/tests/limma-Tests.Rout.save b/tests/limma-Tests.Rout.save
index 6a348fc..f085afd 100755
--- a/tests/limma-Tests.Rout.save
+++ b/tests/limma-Tests.Rout.save
@@ -1,6 +1,6 @@
-R version 3.2.2 (2015-08-14) -- "Fire Safety"
-Copyright (C) 2015 The R Foundation for Statistical Computing
+R version 3.3.0 (2016-05-03) -- "Supposedly Educational"
+Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
@@ -1191,16 +1191,27 @@ set2 5 0.2 0 Down 0.950 0.950 0.250 0.250
NGenes PropDown PropUp Direction PValue FDR PValue.Mixed FDR.Mixed
set1 5 0 1 Up 0.001 0.001 0.001 0.001
set2 5 0 0 Down 0.791 0.791 0.146 0.146
+> fry(y=y,list(set1=iset1,set2=iset2),design,contrast=2,weights=w,array.weights=c(0.5,1,0.5,1))
+ NGenes Direction PValue FDR PValue.Mixed FDR.Mixed
+set1 5 Up 0.0007432594 0.001486519 1.820548e-05 3.641096e-05
+set2 5 Down 0.8208140511 0.820814051 2.211837e-01 2.211837e-01
>
> ### camera
>
-> camera(y=y,iset1,design,contrast=2,weights=c(0.5,1,0.5,1))
+> camera(y=y,iset1,design,contrast=2,weights=c(0.5,1,0.5,1),allow.neg.cor=TRUE,inter.gene.cor=NA)
NGenes Correlation Direction PValue
set1 5 -0.2481655 Up 0.001050253
-> camera(y=y,list(set1=iset1,set2=iset2),design,contrast=2)
+> camera(y=y,list(set1=iset1,set2=iset2),design,contrast=2,allow.neg.cor=TRUE,inter.gene.cor=NA)
NGenes Correlation Direction PValue FDR
set1 5 -0.2481655 Up 0.0009047749 0.00180955
set2 5 0.1719094 Down 0.9068364378 0.90683644
+> camera(y=y,iset1,design,contrast=2,weights=c(0.5,1,0.5,1))
+ NGenes Direction PValue
+set1 5 Up 1.105329e-10
+> camera(y=y,list(set1=iset1,set2=iset2),design,contrast=2)
+ NGenes Direction PValue FDR
+set1 5 Up 7.334400e-12 1.466880e-11
+set2 5 Down 8.677115e-01 8.677115e-01
>
> ### with EList arg
>
@@ -1211,9 +1222,12 @@ Down 0 0.997498749
Up 1 0.003001501
UpOrDown 1 0.006000000
Mixed 1 0.006000000
-> camera(y=y,iset1,design,contrast=2)
+> camera(y=y,iset1,design,contrast=2,allow.neg.cor=TRUE,inter.gene.cor=NA)
NGenes Correlation Direction PValue
set1 5 -0.2481655 Up 0.0009047749
+> camera(y=y,iset1,design,contrast=2)
+ NGenes Direction PValue
+set1 5 Up 7.3344e-12
>
> ### eBayes with trend
>
@@ -1331,30 +1345,41 @@ set1 5 -0.2481655 Up 0.0009047749
+ "137682","137695","137735","1378","137814","137868","137872","137886","137902","137964")
> go <- goana(fit,FDR=0.8,geneid=EB)
> topGO(go,n=10,truncate.term=30)
- Term Ont N Up Down P.Up P.Down
-GO:0055082 cellular chemical homeostas... BP 2 0 2 1.000000000 0.004242424
-GO:0006897 endocytosis BP 3 3 0 0.005046382 1.000000000
-GO:0019725 cellular homeostasis BP 3 0 2 1.000000000 0.012294372
-GO:0006915 apoptotic process BP 4 3 1 0.017844424 0.255402330
-GO:0012501 programmed cell death BP 4 3 1 0.017844424 0.255402330
-GO:0008283 cell proliferation BP 4 1 2 0.553950615 0.023749721
-GO:0042127 regulation of cell prolifer... BP 4 1 2 0.553950615 0.023749721
-GO:0006909 phagocytosis BP 2 2 0 0.030909091 1.000000000
-GO:0007169 transmembrane receptor prot... BP 2 2 0 0.030909091 1.000000000
-GO:0031252 cell leading edge CC 2 2 0 0.030909091 1.000000000
+ Term Ont N Up Down P.Up
+GO:0032502 developmental process BP 24 3 6 0.868520473
+GO:0006915 apoptotic process BP 5 4 1 0.003446627
+GO:0012501 programmed cell death BP 5 4 1 0.003446627
+GO:0055082 cellular chemical homeostas... BP 2 0 2 1.000000000
+GO:0006897 endocytosis BP 3 3 0 0.005046382
+GO:0044767 single-organism development... BP 23 3 5 0.845083285
+GO:0048856 anatomical structure develo... BP 23 3 5 0.845083285
+GO:1902589 single-organism organelle o... BP 7 1 3 0.762502428
+GO:0008219 cell death BP 6 4 1 0.009129968
+GO:0016265 death BP 6 4 1 0.009129968
+ P.Down
+GO:0032502 0.0006606503
+GO:0006915 0.3096959098
+GO:0012501 0.3096959098
+GO:0055082 0.0042424242
+GO:0006897 1.0000000000
+GO:0044767 0.0066515474
+GO:0048856 0.0066515474
+GO:1902589 0.0066732856
+GO:0008219 0.3605604217
+GO:0016265 0.3605604217
> topGO(go,n=10,truncate.term=30,sort="down")
- Term Ont N Up Down P.Up P.Down
-GO:0055082 cellular chemical homeostas... BP 2 0 2 1.0000000 0.004242424
-GO:0019725 cellular homeostasis BP 3 0 2 1.0000000 0.012294372
-GO:0008283 cell proliferation BP 4 1 2 0.5539506 0.023749721
-GO:0042127 regulation of cell prolifer... BP 4 1 2 0.5539506 0.023749721
-GO:1902589 single-organism organelle o... BP 5 1 2 0.6375849 0.038228009
-GO:0048878 chemical homeostasis BP 5 1 2 0.6375849 0.038228009
-GO:0032502 developmental process BP 22 3 4 0.8185069 0.040079017
-GO:0043230 extracellular organelle CC 14 1 3 0.9502416 0.055161144
-GO:0065010 extracellular membrane-boun... CC 14 1 3 0.9502416 0.055161144
-GO:0070062 extracellular exosome CC 14 1 3 0.9502416 0.055161144
+ Term Ont N Up Down P.Up P.Down
+GO:0032502 developmental process BP 24 3 6 0.8685205 0.0006606503
+GO:0055082 cellular chemical homeostas... BP 2 0 2 1.0000000 0.0042424242
+GO:0044767 single-organism development... BP 23 3 5 0.8450833 0.0066515474
+GO:0048856 anatomical structure develo... BP 23 3 5 0.8450833 0.0066515474
+GO:1902589 single-organism organelle o... BP 7 1 3 0.7625024 0.0066732856
+GO:0007283 spermatogenesis BP 3 0 2 1.0000000 0.0122943723
+GO:0048232 male gamete generation BP 3 0 2 1.0000000 0.0122943723
+GO:0001889 liver development BP 3 0 2 1.0000000 0.0122943723
+GO:0019725 cellular homeostasis BP 3 0 2 1.0000000 0.0122943723
+GO:0035295 tube development BP 3 0 2 1.0000000 0.0122943723
>
> proc.time()
user system elapsed
- 3.15 0.09 3.24
+ 2.73 0.20 2.94
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/r-bioc-limma.git
More information about the debian-med-commit
mailing list