[med-svn] [r-bioc-limma] 01/03: Imported Upstream version 3.26.1+dfsg
Andreas Tille
tille at debian.org
Mon Nov 2 10:33:04 UTC 2015
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-bioc-limma.
commit 811b74f1380d7b6e78040e0031fe6773ba98c1cc
Author: Andreas Tille <tille at debian.org>
Date: Mon Nov 2 10:58:33 2015 +0100
Imported Upstream version 3.26.1+dfsg
---
DESCRIPTION | 7 +-
NAMESPACE | 17 +++
R/avearrays.R | 5 +-
R/barcodeplot.R | 10 +-
R/decidetests.R | 2 +-
R/geneset-camera.R | 19 ++--
R/geneset-fry.R | 40 ++++---
R/geneset-roast.R | 13 ++-
R/goana.R | 2 +-
R/kegga.R | 2 +-
R/lmfit.R | 18 ++--
R/read.R | 23 +++--
R/sepchannel.R | 14 ++-
R/tricubeMovingAverage.R | 36 ++++++-
build/vignette.rds | Bin 231 -> 230 bytes
inst/doc/changelog.txt | 88 +++++++++++++---
inst/doc/intro.pdf | Bin 46191 -> 46191 bytes
man/barcodeplot.Rd | 13 +--
man/camera.Rd | 5 +-
man/controlStatus.Rd | 8 +-
man/ebayes.Rd | 4 +-
man/genas.Rd | 2 +-
man/roast.Rd | 2 +
man/romer.Rd | 6 +-
man/tmixture.Rd | 26 ++---
man/tricubeMovingAverage.Rd | 28 +++--
tests/limma-Tests.R | 4 +-
tests/limma-Tests.Rout.save | 246 ++++++++++++++++++++++++++------------------
vignettes/intro.Rnw | 49 +++++++++
29 files changed, 474 insertions(+), 215 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 306a9fa..192fe60 100755
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,12 +1,13 @@
Package: limma
-Version: 3.24.15
-Date: 2015-08-05
+Version: 3.26.1
+Date: 2015-10-30
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]
Maintainer: Gordon Smyth <smyth at wehi.edu.au>
License: GPL (>=2)
Depends: R (>= 2.3.0), methods
+Imports: grDevices, graphics, stats, utils
Suggests: affy, AnnotationDbi, BiasedUrn, Biobase, ellipse, GO.db,
illuminaio, KEGGREST, locfit, MASS, org.Hs.eg.db, splines,
statmod (>= 1.2.2), vsn
@@ -20,4 +21,4 @@ biocViews: ExonArray, GeneExpression, Transcription,
MultipleComparison, Normalization, Preprocessing,
QualityControl
NeedsCompilation: yes
-Packaged: 2015-08-06 01:23:12 UTC; biocbuild
+Packaged: 2015-10-31 00:23:54 UTC; biocbuild
diff --git a/NAMESPACE b/NAMESPACE
index e2431e7..492b181 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -5,6 +5,23 @@ exportClasses("RGList","MAList","EListRaw","EList","MArrayLM","TestResults","Lar
exportMethods("show")
import(methods)
+importFrom("grDevices", "col2rgb", "dev.cur", "dev.off", "png", "rgb")
+importFrom("graphics", "abline", "axis", "barplot", "coplot", "image",
+ "legend", "lines", "matplot", "mtext", "panel.smooth",
+ "par", "plot", "plot.new", "points", "polygon", "rect",
+ "segments", "text", "title")
+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",
+ "lowess", "median", "model.matrix", "na.exclude", "na.omit",
+ "nlminb", "nls", "nls.control", "optim", "p.adjust",
+ "pbeta", "pchisq", "pexp", "pf", "pgamma", "phyper",
+ "pnorm", "ppoints", "predict", "pt", "qchisq", "qf",
+ "qnorm", "qt", "quantile", "residuals", "rnorm", "sd",
+ "uniroot", "var", "weighted.mean")
+importFrom("utils", "read.table", "write.table")
+if( tools:::.OStype() == "windows" ) importFrom("utils", "winMenuAddItem")
S3method("[",MAList)
S3method("[",MArrayLM)
diff --git a/R/avearrays.R b/R/avearrays.R
index 7d94d94..7e769d2 100644
--- a/R/avearrays.R
+++ b/R/avearrays.R
@@ -6,11 +6,12 @@ avearrays <- function(x,ID=NULL,weights=NULL) UseMethod("avearrays")
avearrays.default <- function(x,ID=colnames(x),weights=NULL)
# Average over technical replicate columns, for matrices
# Gordon Smyth
-# Created 24 Sept 2010. Last modified 1 Dec 2010.
+# Created 24 Sept 2010. Last modified 13 Aug 2015.
{
if(is.null(x)) return(NULL)
- if(is.null(ID)) stop("No sample IDs")
+ if(is.vector(x)) stop("x must be a matrix")
x <- as.matrix(x)
+ if(is.null(ID)) stop("No sample IDs")
ID <- as.character(ID)
if(mode(x)=="character") {
d <- duplicated(ID)
diff --git a/R/barcodeplot.R b/R/barcodeplot.R
index f1b7d20..d90931d 100644
--- a/R/barcodeplot.R
+++ b/R/barcodeplot.R
@@ -1,9 +1,9 @@
## BARCODEPLOT.R
-barcodeplot <- function (statistics, index = NULL, index2 = NULL, gene.weights = NULL, weights.label = "Weight", labels = c("Up", "Down"), quantiles = c(-1, 1), col.bars = NULL, worm = TRUE, span.worm = 0.45,...)
+barcodeplot <- function (statistics, index = NULL, index2 = NULL, gene.weights = NULL, weights.label = "Weight", labels = c("Up", "Down"), quantiles = c(-1,1)*sqrt(2), col.bars = NULL, alpha = 0.4, worm = TRUE, span.worm = 0.45,...)
# Barcode plot of one or two gene sets.
# Gordon Smyth, Di Wu and Yifang Hu
-# 20 October 2008. Last revised 7 October 2014.
+# 20 October 2008. Last revised 26 October 2015.
{
# Check statistics
if(!is.vector(statistics, mode = "numeric")) stop("statistics should be a numeric vector")
@@ -183,8 +183,10 @@ barcodeplot <- function (statistics, index = NULL, index2 = NULL, gene.weights =
# Check other arguments
quantiles <- sort(quantiles)
-# color settting
- ALP <- 0.4
+# check transparency settting for vertical bars
+ ALP <- alpha
+ ALP <- min(ALP,1)
+ ALP <- max(ALP,0.1)
if (is.null(col.bars)) {
diff --git a/R/decidetests.R b/R/decidetests.R
index be759dd..be5bf97 100755
--- a/R/decidetests.R
+++ b/R/decidetests.R
@@ -89,7 +89,7 @@ classifyTestsF <- function(object,cor.matrix=NULL,df=Inf,p.value=0.01,fstat.only
# Gordon Smyth
# 20 Mar 2003. Last revised 6 June 2009.
-# Method intended for MAList objects but accept unclassed lists as well
+# Method intended for MArrayLM objects but accept unclassed lists as well
if(is.list(object)) {
if(is.null(object$t)) stop("tstat cannot be extracted from object")
if(is.null(cor.matrix) && !is.null(object$cov.coefficients)) cor.matrix <- cov2cor(object$cov.coefficients)
diff --git a/R/geneset-camera.R b/R/geneset-camera.R
index 248d513..5c9426a 100644
--- a/R/geneset-camera.R
+++ b/R/geneset-camera.R
@@ -24,7 +24,7 @@ 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,...)
# Competitive gene set test allowing for correlation between genes
# Gordon Smyth and Di Wu
-# Created 2007. Last modified 20 July 2015
+# Created 2007. Last modified 30 October 2015
{
# Issue warning if extra arguments found
dots <- names(list(...))
@@ -34,6 +34,8 @@ camera.default <- function(y,index,design=NULL,contrast=ncol(design),weights=NUL
y <- getEAWP(y)
G <- nrow(y$exprs)
n <- ncol(y$exprs)
+ ID <- rownames(y)
+ if(G<3) stop("Two few genes in dataset: need at least 3")
# Check index
if(!is.list(index)) index <- list(set1=index)
@@ -49,10 +51,9 @@ camera.default <- function(y,index,design=NULL,contrast=ncol(design),weights=NUL
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
+ if(df.residual < 1) stop("No residual df: cannot compute t-tests")
# Check weights
if(is.null(weights)) weights <- y$weights
@@ -60,14 +61,15 @@ camera.default <- function(y,index,design=NULL,contrast=ncol(design),weights=NUL
# Check inter.gene.cor
fixed.cor <- !is.null(inter.gene.cor)
-# Set for camera tests
- if(fixed.cor)
+# Set df for camera tests
+ if(fixed.cor) {
if(use.ranks)
df.camera <- Inf
else
df.camera <- G-2
- else
+ } else {
df.camera <- min(df.residual,G-2)
+ }
# Reduce to numeric expression matrix
y <- y$exprs
@@ -112,7 +114,9 @@ camera.default <- function(y,index,design=NULL,contrast=ncol(design),weights=NUL
if(QR$qr[p,p]<0) unscaledt <- -unscaledt
} else {
effects <- matrix(0,n,G)
- unscaledt <- rep(0,n)
+ colnames(effects) <- ID
+ unscaledt <- rep.int(0,G)
+ names(unscaledt) <- ID
sw <- sqrt(weights)
yw <- y*sw
for (g in 1:G) {
@@ -150,6 +154,7 @@ camera.default <- function(y,index,design=NULL,contrast=ncol(design),weights=NUL
colnames(tab) <- c("NGenes","Correlation","Down","Up","TwoSided")
for (i in 1:nsets) {
iset <- index[[i]]
+ if(is.character(iset)) iset <- which(ID %in% iset)
StatInSet <- Stat[iset]
m <- length(StatInSet)
m2 <- G-m
diff --git a/R/geneset-fry.R b/R/geneset-fry.R
index 4a6eddd..cc969fd 100644
--- a/R/geneset-fry.R
+++ b/R/geneset-fry.R
@@ -5,7 +5,7 @@ fry.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),weights=N
# 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 23 June 2015
+# Created 30 January 2015. Last modified 6 October 2015
{
# Issue warning if extra arguments found
dots <- names(list(...))
@@ -36,7 +36,6 @@ fry.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),weights=N
p <- ncol(design)
df.residual <- n-p
- df.camera <- min(df.residual,G-2)
# Check weights
if(is.null(weights)) weights <- y$weights
@@ -98,7 +97,7 @@ fry.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),weights=N
}
# Standardized residuals
- U <- effects[-(1:(p-1)),,drop=FALSE]
+ U <- t(effects[-(1:(p-1)),,drop=FALSE])
# Global statistics
nsets <- length(index)
@@ -107,24 +106,41 @@ fry.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),weights=N
PValue.Mixed <- PValue <- rep.int(0,nsets)
for (i in 1:nsets) {
iset <- index[[i]]
- USet <- U[,iset,drop=FALSE]
- NGenes[i] <- ncol(USet)
- MeanUSet <- rowMeans(USet)
+ 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)
-# MSqUSet <- rowMeans(USet^2)
-# F.stat <- MSqUSet[1L] / mean(MSqUSet[-1L])
-# PValue.Mixed[i] <- pf(F.stat,df1=1L,df2=df.residual,lower.tail=FALSE)
+
+ if(NGenes[i]==1) {
+ PValue.Mixed[i] <- PValue[i]
+ } else {
+ SVD <- svd(USet,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])
+ Frb.mean <- (sum(A) * beta.mean - A[d1]) / (A[1]-A[d1])
+ COV <- matrix(-beta.var/d,d1,d1)
+ diag(COV) <- beta.var
+ Frb.var <- (A %*% COV %*% A ) / (A[1]-A[d1])^2
+ alphaplusbeta <- Frb.mean*(1-Frb.mean)/Frb.var-1
+ alpha <- alphaplusbeta*Frb.mean
+ beta <- alphaplusbeta-alpha
+ PValue.Mixed[i] <- pbeta(Fobs,shape1=alpha,shape2=beta,lower.tail=FALSE)
+ }
}
# Add FDR
if(nsets>1) {
FDR <- p.adjust(PValue,method="BH")
-# FDR.Mixed <- p.adjust(PValue.Mixed,method="BH")
- tab <- data.frame(NGenes=NGenes,Direction=Direction,PValue=PValue,FDR=FDR)
+ FDR.Mixed <- p.adjust(PValue.Mixed,method="BH")
+ tab <- data.frame(NGenes=NGenes,Direction=Direction,PValue=PValue,FDR=FDR,PValue.Mixed=PValue.Mixed,FDR.Mixed=FDR.Mixed)
} else {
- tab <- data.frame(NGenes=NGenes,Direction=Direction,PValue=PValue)
+ tab <- data.frame(NGenes=NGenes,Direction=Direction,PValue=PValue,PValue.Mixed=PValue.Mixed)
}
rownames(tab) <- names(index)
diff --git a/R/geneset-roast.R b/R/geneset-roast.R
index c56f8b0..a705747 100644
--- a/R/geneset-roast.R
+++ b/R/geneset-roast.R
@@ -357,7 +357,7 @@ 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 23 June 2015.
+# Created 28 Jan 2010. Last revised 7 Oct 2015.
{
# Extract components from y
y <- getEAWP(y)
@@ -494,10 +494,13 @@ mroast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.st
# Sort by p-value
sort <- match.arg(sort,c("directional","mixed","none"))
if(sort=="none") return(tab)
- if(sort=="directional")
- o <- order(tab$PValue,tab$PValue.Mixed,-tab$NGenes)
- else
- o <- order(tab$PValue.Mixed,tab$PValue,-tab$NGenes)
+ 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]
}
diff --git a/R/goana.R b/R/goana.R
index 3860dd0..114e88f 100644
--- a/R/goana.R
+++ b/R/goana.R
@@ -108,7 +108,7 @@ goana.default <- function(de, universe = NULL, species = "Hs", prior.prob = NULL
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,full.length=TRUE)
+ prior.prob[o] <- tricubeMovingAverage(isDE[o],span=span)
if(plot) barcodeplot(covariate, index=(isDE==1), worm=TRUE, span.worm=span)
}
diff --git a/R/kegga.R b/R/kegga.R
index 4dd3d76..0ef60f7 100644
--- a/R/kegga.R
+++ b/R/kegga.R
@@ -108,7 +108,7 @@ kegga.default <- function(de, universe = NULL, species = "Hs", prior.prob = NULL
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,full.length=TRUE)
+ prior.prob[o] <- tricubeMovingAverage(isDE[o],span=span)
if(plot) barcodeplot(covariate, index=(isDE==1), worm=TRUE, span.worm=span)
}
diff --git a/R/lmfit.R b/R/lmfit.R
index 82904f7..c1f80ab 100644
--- a/R/lmfit.R
+++ b/R/lmfit.R
@@ -3,7 +3,7 @@
lmFit <- function(object,design=NULL,ndups=1,spacing=1,block=NULL,correlation,weights=NULL,method="ls",...)
# Fit genewise linear models
# Gordon Smyth
-# 30 June 2003. Last modified 7 July 2015.
+# 30 June 2003. Last modified 6 Oct 2015.
{
# Extract components from y
y <- getEAWP(object)
@@ -51,9 +51,8 @@ lmFit <- function(object,design=NULL,ndups=1,spacing=1,block=NULL,correlation,we
# Possible warning on missing coefs
if(NCOL(fit$coef)>1) {
- i <- is.na(fit$coef)
- i <- apply(i[,1]==i[,-1,drop=FALSE],1,all)
- n <- sum(!i)
+ n <- rowSums(is.na(fit$coef))
+ n <- sum(n>0 & n<NCOL(fit$coef))
if(n>0) warning("Partial NA coefficients for ",n," probe(s)",call.=FALSE)
}
@@ -391,7 +390,7 @@ getEAWP <- function(object)
# Given any microarray data object, extract basic information needed
# for linear modelling.
# Gordon Smyth
-# 9 March 2008. Last modified 7 July 2015.
+# 9 March 2008. Last modified 18 August 2015.
{
y <- list()
@@ -431,13 +430,20 @@ getEAWP <- function(object)
}
if(length(object at maA)) y$Amean <- rowMeans(object at maA,na.rm=TRUE)
} else {
+ if(is(object,"eSet")) {
+ if(!requireNamespace("Biobase",quietly=TRUE)) stop("Biobase package required but is not available")
+ y$exprs <- Biobase::assayDataElement(object,"exprs")
+ if(length(object at featureData@data)) y$probes <- object at featureData@data
+ y$Amean <- rowMeans(y$exprs,na.rm=TRUE)
+ if("weights" %in% Biobase::assayDataElementNames(object)) y$weights <- Biobase::assayDataElement(object,"weights")
+ } else {
# Default method for matrices, data.frames, vsn objects etc.
if(is.vector(object))
y$exprs <- matrix(object,nrow=1)
else
y$exprs <- as.matrix(object)
y$Amean <- rowMeans(y$exprs,na.rm=TRUE)
- }}}}
+ }}}}}
# Check expression values are numeric
if(mode(y$exprs) != "numeric") stop("Data object doesn't contain numeric expression values")
diff --git a/R/read.R b/R/read.R
index 9144df4..cf31e7e 100755
--- a/R/read.R
+++ b/R/read.R
@@ -177,26 +177,35 @@ readSpotTypes <- function(file="SpotTypes.txt",path=NULL,sep="\t",check.names=FA
}
controlStatus <- function(types, genes, spottypecol="SpotType", regexpcol, verbose=TRUE)
-# Set status of each spot
+# Set status of each spot (row)
# Gordon Smyth
-# 19 Oct 2003. Last modified 9 Nov 2003.
+# 19 Oct 2003. Last modified 19 Sep 2015.
{
-# Check input
- if (is(genes, "RGList") || is(genes, "MAList") || is(genes, "MArrayLM")) genes <- genes$genes
+# Check types
+ if(!is.data.frame(types)) stop("types must be a data.frame")
cntypes <- colnames(types)
+ if(is.null(cntypes)) stop("types must have column names")
+
+# Check genes
+ if (is(genes,"RGList") || is(genes,"MAList") || is(genes,"MArrayLM") || is(genes,"EListRaw") || is(genes,"EList")) genes <- genes$genes
+ if(!is.data.frame(genes)) stop("genes must be a data.frame")
cngenes <- colnames(genes)
- if(is.null(cntypes) || is.null(cngenes)) stop("types and genes must have column names")
+ if(is.null(cngenes)) stop("genes must have column names")
+
+# Check have at least one row
ntypes <- nrow(types)
nspots <- nrow(genes)
if(ntypes==0 || nspots==0) return(NULL)
+
# Any undo conversion of types to factors
for (j in cntypes) {
x <- types[,j]
if(is.factor(x) && is.character(levels(x))) types[,j] <- as.character(x)
}
+
+# Check spottypecol
if(is.numeric(spottypecol)) spottypecol <- cntypes[spottypecol[1]]
- if(is.null(spottypecol) || !is.character(spottypecol) || !is.element(spottypecol,cntypes))
- stop("spottypecol not valid column of types")
+ if(is.null(spottypecol) || !is.character(spottypecol) || !is.element(spottypecol,cntypes)) stop("spottypecol not valid column of types")
# Find common columns between types and genes
if(!missing(regexpcol) && is.numeric(regexpcol)) regexpcol <- cntypes[regexpcol]
diff --git a/R/sepchannel.R b/R/sepchannel.R
index 30fc95d..1991562 100755
--- a/R/sepchannel.R
+++ b/R/sepchannel.R
@@ -4,22 +4,30 @@ lmscFit <- function(object,design,correlation)
# Fit single channel linear model for each gene to a series of microarrays
# allowing for known correlation between the channels on each spot.
# Gordon Smyth
-# 14 March 2004. Last modified 26 June 2004.
+# 14 March 2004. Last modified 1 Oct 2015.
{
-# Check input
+# Check object (assumed to a MAList)
+ if(!is.list(object)) stop("object should be an MAList")
+ if(is.null(object$M) || is.null(object$A)) stop("object must have components M and A")
M <- as.matrix(object$M)
A <- as.matrix(object$A)
- if(is.null(M) || is.null(A)) stop("object must have components M and A")
dimM <- dim(M)
dimA <- dim(A)
if(any(dimM != dimA)) stop("dimensions of M and A don't match")
if(!all(is.finite(M)) || !all(is.finite(A))) stop("Missing or infinite values found in M or A")
+ if(!is.null(object$weights)) warning("weights found in object but not used")
+
+# Check design
if(missing(design)) stop("design matrix must be specified")
narrays <- dimM[2]
ny <- 2*narrays
design <- as.matrix(design)
if(nrow(design) != ny) stop("The number of rows of the design matrix should match the number of channel intensities, i.e., twice the number of arrays")
+
+# Check correlation
if(missing(correlation)) stop("intra-spot correlation must be specified")
+ correlation <- as.vector(correlation)[1]
+ if(is.na(correlation)) stop("intra-spot correlation is NA")
if(abs(correlation) >= 1) stop("correlation must be strictly between -1 and 1")
# Dimensions
diff --git a/R/tricubeMovingAverage.R b/R/tricubeMovingAverage.R
index 0ab341b..46c51b6 100644
--- a/R/tricubeMovingAverage.R
+++ b/R/tricubeMovingAverage.R
@@ -1,18 +1,48 @@
-tricubeMovingAverage <- function(x,span=0.5,full.length=TRUE) {
+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.
+{
+# Check span
+ if(length(span)>1) {
+ warning("only first value of span used")
+ span <- span[1]
+ }
+ if(span>1) span <- 1
+ if(span<=0) return(x)
+
+# Check power
+ if(length(power)>1) {
+ warning("only first value of power used")
+ span <- span[1]
+ }
+ if(power<0) power <- 0
+
+# Convert span to width of moving window
n <- length(x)
width <- span*n
+
# 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
+
+# Tricube weights with all positive values
u <- seq(from=-1,to=1,length=width) * width / (width+1)
- tricube.weights <- (1-abs(u)^3)^3
+ tricube.weights <- (1-abs(u)^3)^power
tricube.weights <- tricube.weights/sum(tricube.weights)
- if(!full.length) return(as.vector(filter(x,tricube.weights),mode="numeric")[(hwidth+1):(n-hwidth)])
+
+# Extend x with zeros and compute weighted moving averages
z <- numeric(hwidth)
x <- as.vector(filter(c(z,x,z),tricube.weights),mode="numeric")[(hwidth+1):(n+hwidth)]
+
+# Rescale boundary values to remove influence of the outside values of zero
cw <- cumsum(tricube.weights)
x[1:hwidth] <- x[1:hwidth] / cw[(width-hwidth):(width-1)]
x[(n-hwidth+1):n] <- x[(n-hwidth+1):n] / cw[(width-1):(width-hwidth)]
+
x
}
diff --git a/build/vignette.rds b/build/vignette.rds
index fe28f8c..e4054da 100644
Binary files a/build/vignette.rds and b/build/vignette.rds differ
diff --git a/inst/doc/changelog.txt b/inst/doc/changelog.txt
index dfd4b2c..5fb5e36 100755
--- a/inst/doc/changelog.txt
+++ b/inst/doc/changelog.txt
@@ -1,8 +1,65 @@
- 5 Aug 2015: limma 3.24.15
+30 Oct 2015: limma 3.26.1
+
+- Functions needed from the grDevices, graphics, stats and utils
+ packages are now imported explicitly into the NAMESPACE.
+
+- Two changes to barcodeplot(): new argument 'alpha' to set
+ 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.
+
+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
+ 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.
+
+ 6 Oct 2015: limma 3.25.17
+
+- Checking for NA coefficients in lmFit() is now done more
+ efficiently.
+
+- bug fix for fry() when the gene set has a small number of genes.
+
+- more informative error message from lmscFit() when correlation is
+ NA.
+
+19 Sep 2015: limma 3.25.16
+
+- fry() now produces "mixed" (non-directional) p-values and FDRs as
+ well as directional p-values and FDRs.
+
+- 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.
+
+- the 'genes' argument of controlStatus() now accepts EListRaw or
+ Elist objects.
+
+- avearrays() now gives a more informative error message if input
+ argument 'x' is not a matrix.
+
+- getEAWP() now works on all eSet objects, not just ExpressionSet,
+ provided there is an data element called "exprs".
+
+ 5 Aug 2015: limma 3.25.15
- argument truncate.term for topKEGG() renamed to 'truncate.path'.
-20 July 2015: limma 3.24.14
+20 July 2015: limma 3.25.14
- camera() has a new argument inter.gene.cor. This allows a preset
inter gene correlation to be set for each test test, resulting in
@@ -20,7 +77,7 @@
- lmFit() now treats infinite expression values as NA. In the past,
this was done somewhat inconsistently.
- 7 July 2015: limma 3.24.13
+ 7 July 2015: limma 3.25.13
- New 'tolerance' argument for read.idat() to allow manifest and idat
to differ by a certain number of probes
@@ -29,13 +86,13 @@
vector data 'object' is now coerced to a matrix with one row
instead of one column
-26 June 2015: limma 3.24.12
+26 June 2015: limma 3.25.12
- 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.
-23 June 2015: limma 3.24.11
+23 June 2015: limma 3.25.11
- camera(), fry(), mroast() and romer() now give a more informative
error message when the index list is empty.
@@ -43,41 +100,41 @@
- camera() no longer gives NA results when a gene has a residual
variance exactly zero.
- 9 June 2015: limma 3.24.10
+23 June 2015: limma 3.25.10
- bug fix to write.fit() when method="global" and fit contains more
than one column.
- 8 June 2015: limma 3.24.9
+ 8 June 2015: limma 3.25.9
- plotMD() can now accept arguments 'array' or 'coef' as appropriate
as synonyms for 'column'. This may help users transitioning from
plotMA ot plotMD.
- 7 June 2015: limma 3.24.8
+ 7 June 2015: limma 3.25.8
- 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.
- 5 June 2015: limma 3.24.7
+ 5 June 2015: limma 3.25.7
- The dimnames<- method for EListRaw objects now sets rownames for
the background matrix Eb as well as for the foreground matrix E.
- 3 June 2015: limma 3.24.6
+ 3 June 2015: limma 3.25.6
- new function kegga() to conduct KEGG pathway over-representation
analyses. kegga() is an S3 generic function with a method for
MArrayLM fitted model objects.
-20 May 2015: limma 3.24.5
+20 May 2015: limma 3.25.5
- bug fix for plotWithHighlights(). Legend was not displaying
correctly when there were no background points.
-13 May 2015: limma 3.24.4
+20 May 2015: limma 3.25.4
- goana() with trend now chooses the span of the tricube smoother
based on the number of differentially expressed genes. Larger
@@ -89,18 +146,18 @@
- bug fix to plotMD.RGList.
-28 April 2015: limma 3.24.3
+28 April 2015: limma 3.25.3
- arguments hi.pch, hi.col and hi.cex of plotWithHighlights() renamed
to hl.pch, hl.col and hl.cex. ('hl' is short for 'highlight'.)
- fix bug in auROC() introduced in 3.25.2.
-25 April 2015: limma 3.24.2
+25 April 2015: limma 3.25.2
- auROC() now allows for tied stat values.
-22 April 2015: limma 3.24.1
+22 April 2015: limma 3.25.1
- Update citation for the limma review article (Ritchie et al,
NAR 2015).
@@ -110,6 +167,7 @@
- topGO() has a new argument 'truncate.term' to limit number of
characters appearing in a column.
+17 April 2015: limma 3.25.0 (Bioconductor 3.2 Developmental Branch)
17 April 2015: limma 3.24.0 (Bioconductor 3.1 Release Branch)
14 April 2015: limma 3.23.17
diff --git a/inst/doc/intro.pdf b/inst/doc/intro.pdf
index 65a3286..051de65 100644
Binary files a/inst/doc/intro.pdf and b/inst/doc/intro.pdf differ
diff --git a/man/barcodeplot.Rd b/man/barcodeplot.Rd
index 8c9174b..a15d859 100644
--- a/man/barcodeplot.Rd
+++ b/man/barcodeplot.Rd
@@ -6,8 +6,9 @@ Display the enrichment of one or two gene sets in a ranked gene list.
}
\usage{
barcodeplot(statistics, index = NULL, index2 = NULL, gene.weights = NULL,
- weights.label = "Weight", labels = c("Up","Down"), quantiles = c(-1,1),
- col.bars = NULL, worm = TRUE, span.worm=0.45, \dots)
+ weights.label = "Weight", labels = c("Up","Down"),
+ quantiles = c(-1,1)*sqrt(2), col.bars = NULL, alpha = 0.4,
+ worm = TRUE, span.worm=0.45, \dots)
}
\arguments{
\item{statistics}{numeric vector giving the values of statistics to rank genes by.}
@@ -18,14 +19,14 @@ barcodeplot(statistics, index = NULL, index2 = NULL, gene.weights = NULL,
If specified, then \code{index} and \code{index2} specify positive and negative genes respectively.
Usually used to distinguish down-regulated genes from up-regulated genes.}
\item{gene.weights}{numeric vector giving directional weights for the genes in the (first) set.
- Postive and negative weights correspond to positive and negative genes.
+ Positive and negative weights correspond to positive and negative genes.
Ignored if \code{index2} is non-null.}
\item{weights.label}{label describing the entries in \code{gene.weights}.}
\item{labels}{character vector of labels for high and low statistics. First label is associated with high statistics and is displayed at the left end of the plot. Second label is associated with low or negative statistics and is displayed at the right end of the plot.}
\item{quantiles}{numeric vector of length 2, giving cutoff values for \code{statistics} considered small or large respectively. Used to color the rectangle of the barcodeplot.}
- \item{col.bars}{character vector giving colors for the bars on the barcodeplot.
- Defaults to \code{"black"} for one set or \code{c("red","blue")} for two sets.
- Defaults to semitransparent color for the bars inside the rectangle when variable gene weights are given which is intended to distinguish the positional bars from the weighted bars and also to show the density of the genes.}
+ \item{col.bars}{character vector of colors for the vertical bars of the barcodeplot showing the ranks of the gene set members.
+ Defaults to \code{"black"} for one set or \code{c("red","blue")} for two sets.}
+ \item{alpha}{transparency for vertical bars. When \code{gene.weights} are not \code{NULL}, values \code{0<alpha<1} give semitransparent colors for the vertical bars inside the rectangle. This helps distinguish position bars from the weighted bars and also helps to show the density of the bars when there are many bars. Ignored if \code{gene.weights=NULL}.}
\item{worm}{logical, should enrichment worms be plotted?}
\item{span.worm}{loess span for enrichment worms. Larger spans give smoother worms.}
\item{\ldots}{other arguments are passed to \code{plot}.}
diff --git a/man/camera.Rd b/man/camera.Rd
index f60e23d..fbb5290 100644
--- a/man/camera.Rd
+++ b/man/camera.Rd
@@ -44,8 +44,8 @@ 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.05}.
-This produces a less conservative test.
+\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.
}
\value{
@@ -100,6 +100,7 @@ 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=0.01)
}
\keyword{gene set test}
diff --git a/man/controlStatus.Rd b/man/controlStatus.Rd
index c699703..93bc003 100755
--- a/man/controlStatus.Rd
+++ b/man/controlStatus.Rd
@@ -8,12 +8,12 @@ Determine the type (or status) of each spot in the gene list.
controlStatus(types, genes, spottypecol="SpotType", regexpcol, verbose=TRUE)
}
\arguments{
- \item{types}{dataframe containing spot type specifiers, usually input using \code{readSpotTypes}}
- \item{genes}{dataframe containing the microarray gene list, or an \code{RGList}, \code{MAList} or \code{MArrayList} containing \code{genes}}
- \item{spottypecol}{integer or name specifying column of \code{types} containing spot type names}
+ \item{types}{dataframe containing spot type specifiers, usually input using \code{readSpotTypes}.}
+ \item{genes}{dataframe containing gene annotation, or an object of class \code{RGList}, \code{MAList}, \code{EListRaw}, \code{EList} or \code{MArrayLM} from which the gene annotation can be extracted.}
+ \item{spottypecol}{integer or name specifying column of \code{types} containing spot type names.}
\item{regexpcol}{vector of integers or column names specifying columns of types containing regular expressions.
Defaults to any column names in common between \code{types} and \code{genes}.}
- \item{verbose}{logical, if \code{TRUE} then progess on pattern matching is reported to the standard output channel}
+ \item{verbose}{logical, if \code{TRUE} then progess on pattern matching is reported to the standard output channel.}
}
\details{
This function constructs a vector of status codes by searching for patterns in the gene list.
diff --git a/man/ebayes.Rd b/man/ebayes.Rd
index 2eca21d..fdd2088 100755
--- a/man/ebayes.Rd
+++ b/man/ebayes.Rd
@@ -37,7 +37,7 @@ treat(fit, lfc=0, trend=FALSE, robust=FALSE, winsor.tail.p=c(0.05,0.1))
\code{ebayes} produces an ordinary list containing the above components except for \code{F} and \code{F.p.value}.
}
\details{
-These functions is used to rank genes in order of evidence for differential expression.
+These functions are used to rank genes in order of evidence for differential expression.
They use an empirical Bayes method to shrink the probe-wise sample variances towards a common value and to augmenting the degrees of freedom for the individual variances (Smyth, 2004).
The functions accept as input argument \code{fit} a fitted model object from the functions \code{lmFit}, \code{lm.series}, \code{mrlm} or \code{gls.series}.
The fitted model object may have been processed by \code{contrasts.fit} before being passed to \code{eBayes} to convert the coefficients of the design matrix into an arbitrary number of contrasts which are to be tested equal to zero.
@@ -56,7 +56,7 @@ The F-statistics \code{F} are computed by \code{classifyTestsF} with \code{fstat
\code{eBayes} doesn't compute ordinary (unmoderated) t-statistics by default, but these can be easily extracted from
the linear model output, see the example below.
-\code{ebayes} is the earlier and leaner function, kept for background capatability, while
+\code{ebayes} is the earlier and leaner function, kept for backwards compatibility, while
\code{eBayes} is the later more object-orientated version.
The difference is that \code{ebayes} outputs only the empirical Bayes statistics whereas \code{eBayes} adds them to the fitted model object \code{fit}.
\code{eBayes} is recommended for routine use as it produces objects containing all the necessary components for downstream analysis
diff --git a/man/genas.Rd b/man/genas.Rd
index abc02f8..65e4590 100644
--- a/man/genas.Rd
+++ b/man/genas.Rd
@@ -66,7 +66,7 @@ It does not work with \code{voom} pipelines, because these involve observation w
Majewski, IJ, Ritchie, ME, Phipson, B, Corbin, J, Pakusch, M, Ebert, A, Busslinger, M, Koseki, H, Hu, Y, Smyth, GK, Alexander, WS, Hilton, DJ, and Blewitt, ME (2010).
Opposing roles of polycomb repressive complexes in hematopoietic stem and progenitor cells.
\emph{Blood} 116, 731-739.
-\url{http://bloodjournal.hematologylibrary.org/content/116/5/731}
+\url{http://www.bloodjournal.org/content/116/5/731}
Phipson, B. (2013).
\emph{Empirical Bayes modelling of expression profiles and their associations}.
diff --git a/man/roast.Rd b/man/roast.Rd
index 6fbb90d..374afc6 100644
--- a/man/roast.Rd
+++ b/man/roast.Rd
@@ -68,6 +68,8 @@ Rows correspond to the alternative hypotheses Down, Up, UpOrDown (two-sided) and
\item{FDR}{two-sided directional false discovery rate}
\item{PValue.Mixed}{non-directional p-value}
\item{FDR.Mixed}{non-directional false discovery rate}
+
+\code{fry} produces the same output format as \code{mroast} but without the columns \code{PropDown} and \code{ProbUp}.
}
\details{
diff --git a/man/romer.Rd b/man/romer.Rd
index b5ec1b8..1c85c97 100644
--- a/man/romer.Rd
+++ b/man/romer.Rd
@@ -12,7 +12,7 @@ Gene set enrichment analysis for linear models using rotation tests (ROtation te
}
\arguments{
\item{y}{numeric matrix giving log-expression values.}
- \item{index}{list of indices specifying the rows of \code{y} in the gene sets. The list can be made using \link{ids2indices}.}
+ \item{index}{list of indices specifying the rows of \code{y} in the gene sets. The list can be made using \code{\link{ids2indices}}.}
\item{design}{design matrix.}
\item{contrast}{contrast for which the test is required. Can be an integer specifying a column of \code{design}, or else a contrast vector of length equal to the number of columns of \code{design}.}
\item{array.weights}{optional numeric vector of array weights.}
@@ -36,10 +36,6 @@ Like GSEA, it is a competitive test in that the different gene sets are pitted a
Instead of permutation, it uses rotation, a parametric resampling method suitable for linear models (Langsrud, 2005; Wu et al, 2010).
\code{romer} can be used with any linear model with some level of replication.
-Curated gene sets suitable for use with \code{romer} can be downloaded from
-\url{http://bioinf.wehi.edu.au/software/MSigDB/}.
-These lists are based on the molecular signatures database from the Broad Institute, but with gene symbols converted to offical gene symbols, separately for mouse and human.
-
In the output, p-values are given for each set for three possible alternative hypotheses.
The alternative "up" means the genes in the set tend to be up-regulated, with positive t-statistics.
The alternative "down" means the genes in the set tend to be down-regulated, with negative t-statistics.
diff --git a/man/tmixture.Rd b/man/tmixture.Rd
index edf5d09..82bcef2 100755
--- a/man/tmixture.Rd
+++ b/man/tmixture.Rd
@@ -3,29 +3,29 @@
\alias{tmixture.matrix}
\title{Estimate Scale Factor in Mixture of t-Distributions}
\description{
-This function estimates the unscaled standard deviation of the log fold change for differentially expressed genes.
-It is called by the function \code{ebayes} and is not intended to be called by users.
+These functions estimate the unscaled standard deviation of the true (unobserved) log fold changes for differentially expressed genes.
+They are used by the functions \code{ebayes} and \code{eBayes} and are not intended to be called directly by users.
}
\usage{
-tmixture.vector(tstat,stdev.unscaled,df,proportion,v0.lim=NULL)
-tmixture.matrix(tstat,stdev.unscaled,df,proportion,v0.lim=NULL)
+tmixture.vector(tstat, stdev.unscaled, df, proportion, v0.lim = NULL)
+tmixture.matrix(tstat, stdev.unscaled, df, proportion, v0.lim = NULL)
}
\arguments{
- \item{tstat}{numeric vector or matrix of t-statistics}
- \item{stdev.unscaled}{numeric matrix conformal with \code{tstatf} containing the unscaled standard deviations for the coefficient estimators}
- \item{df}{numeric vector giving the degrees of freedom associated with \code{tstat}}
- \item{proportion}{assumed proportion of genes which are differentially expressed}
- \item{v0.lim}{numeric vector of length 2, assumed lower and upper limits for the estimated unscaled standard deviation}
+ \item{tstat}{numeric vector or matrix of t-statistics. \code{tmixture.vector} assumes a vector while \code{tmixture.matrix} assumes a matrix.}
+ \item{stdev.unscaled}{numeric vector or matrix, conformal with \code{tstat}, containing the unscaled standard deviations of the coefficients used to compute the t-statistics.}
+ \item{df}{numeric vector giving the degrees of freedom associated with \code{tstat}.}
+ \item{proportion}{assumed proportion of genes that are differentially expressed.}
+ \item{v0.lim}{numeric vector of length 2 giving the lower and upper limits for the estimated unscaled standard deviations.}
}
\details{
The values in each column of \code{tstat} are assumed to follow a mixture of an ordinary t-distribution, with mixing proportion \code{1-proportion}, and \code{(v0+v1)/v1} times a t-distribution, with mixing proportion \code{proportion}.
-Here \code{v1=stdev.unscaled^2} and \code{v0} is the value to be estimated.
+Here \code{v1} is \code{stdev.unscaled^2} and \code{v0} is the value to be estimated.
}
\value{
-Numeric vector of length equal to the number of columns of \code{tstat} and \code{stdev.unscaled}.
+Numeric vector, of length equal to the number of columns of \code{tstat}, containing estimated \code{v0} values.
}
\seealso{
-\code{\link{ebayes}}
+\code{\link{eBayes}}
}
\author{Gordon Smyth}
-\keyword{htest}
+\keyword{distribution}
diff --git a/man/tricubeMovingAverage.Rd b/man/tricubeMovingAverage.Rd
index e742412..1e7413b 100644
--- a/man/tricubeMovingAverage.Rd
+++ b/man/tricubeMovingAverage.Rd
@@ -1,38 +1,44 @@
\name{tricubeMovingAverage}
\alias{tricubeMovingAverage}
-\alias{tricubeMovingAverage}
\title{Moving Average Smoother With Tricube Weights}
\description{
-Apply a moving average smoother with tricube distance weights to the columns of a matrix.
+Apply a moving average smoother with tricube distance weights to a numeric vector.
}
\usage{
-tricubeMovingAverage(x,span=0.5,full.length=TRUE)
+tricubeMovingAverage(x, span=0.5, power=3)
}
\arguments{
\item{x}{numeric vector}
- \item{span}{proportion of points included in the local window}
- \item{full.length}{logical value, should output have same number of length as input?}
+ \item{span}{the smoother span. This gives the proportion of \code{x} values that contribute to each moving average. Larger values give more smoothness. Should be positive but not greater than 1.}
+ \item{power}{a positive exponent used to compute the tricube weights. \code{power=3} gives the usual tricube weights. Smaller values give more even weighting. Should be greater than 0.}
}
\details{
This function smooths a vector (considered as a time series) using a moving average with tricube weights.
-This is similar to a loess curve of degree zero, with a couple of differences:
-a continuity correction is applied when computing the neighbouring points and, when \code{full.length=TRUE}, the span halves at the end points.
+Specifically, the function computes running weighted means of \code{w} consecutive values of \code{x}, where the window width \code{w} is equal to \code{2*h+1} with \code{h = 2*floor(span*length(x)/2)}.
+The window width \code{w} is always odd so that each window has one of the original \code{x} values at its center.
+Each weighted mean uses a set of tricube weights so that values near the ends of the window receive less weight.
+
+The smoother returns a vector of the same length as input.
+At the start and end of the vector, the series is considered to be extended by missing values, and the weighted average is computed only over the observed values.
+In other words, the window width is reduced to \code{h+1} at the boundaries with asymmetric weights.
+
+The result of this function is similar to a least squares loess curve of degree zero, with a couple of differences.
+First, a continuity correction is applied when computing the distance to neighbouring points, so that exactly \code{w} points are included with positive weights in each average.
+Second, the span halves at the end points so that the smoother is more sensitive to trends at the ends.
The \code{filter} function in the stats package is called to do the low-level calculations.
This function is used by \code{\link{barcodeplot}} to compute enrichment worms.
}
\value{
-Numeric vector of smoothed values.
-If \code{full.length=TRUE}, of same length as \code{x}.
-If \code{full.length=FALSE}, has \code{width-1} fewer rows than \code{x}.
+Numeric vector of same length as \code{x} containing smoothed values.
}
\examples{
x <- rbinom(100,size=1,prob=0.5)
plot(1:100,tricubeMovingAverage(x))
}
\seealso{
-\code{\link{filter}}, \code{\link{barcodeplot}}
+\code{\link{filter}}, \code{\link{barcodeplot}}, \code{\link[edgeR:loessByCol]{loessByCol}}
}
\author{Gordon Smyth}
\keyword{smooth}
diff --git a/tests/limma-Tests.R b/tests/limma-Tests.R
index 63597d3..1bcf9d1 100755
--- a/tests/limma-Tests.R
+++ b/tests/limma-Tests.R
@@ -284,5 +284,5 @@ EB <- c("133746","1339","134","1340","134083","134111","134147","134187","134218
"1370","137075","1371","137209","1373","137362","1374","137492","1375","1376",
"137682","137695","137735","1378","137814","137868","137872","137886","137902","137964")
go <- goana(fit,FDR=0.8,geneid=EB)
-topGO(go,n=10)
-topGO(go,n=10,sort="down")
+topGO(go,n=10,truncate.term=30)
+topGO(go,n=10,truncate.term=30,sort="down")
diff --git a/tests/limma-Tests.Rout.save b/tests/limma-Tests.Rout.save
index 1ac4959..6a348fc 100755
--- a/tests/limma-Tests.Rout.save
+++ b/tests/limma-Tests.Rout.save
@@ -1,5 +1,5 @@
-R version 3.2.1 (2015-06-18) -- "World-Famous Astronaut"
+R version 3.2.2 (2015-08-14) -- "Fire Safety"
Copyright (C) 2015 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)
@@ -1122,7 +1122,127 @@ attr(,"df1")
attr(,"df2")
[1] Inf
> classifyTestsT(tstat)
-TestRe.2900320 0.3006936 0.2935101 0.3646949 0.3596385 0.3064203
+TestResults matrix
+ [,1] [,2] [,3]
+[1,] 0 1 0
+[2,] 0 0 0
+[3,] 0 0 0
+[4,] 0 0 0
+> classifyTestsP(tstat)
+TestResults matrix
+ [,1] [,2] [,3]
+[1,] 0 1 0
+[2,] 0 1 0
+[3,] 0 0 0
+[4,] 0 0 0
+>
+> ### avereps
+>
+> x <- matrix(rnorm(8*3),8,3)
+> colnames(x) <- c("S1","S2","S3")
+> rownames(x) <- c("b","a","a","c","c","b","b","b")
+> avereps(x)
+ S1 S2 S3
+b -0.2353018 0.5220094 0.2302895
+a -0.4347701 0.6453498 -0.6758914
+c 0.3482980 -0.4820695 -0.3841313
+>
+> ### roast
+>
+> y <- matrix(rnorm(100*4),100,4)
+> sigma <- sqrt(2/rchisq(100,df=7))
+> y <- y*sigma
+> design <- cbind(Intercept=1,Group=c(0,0,1,1))
+> iset1 <- 1:5
+> y[iset1,3:4] <- y[iset1,3:4]+3
+> iset2 <- 6:10
+> roast(y=y,iset1,design,contrast=2)
+ Active.Prop P.Value
+Down 0 0.996498249
+Up 1 0.004002001
+UpOrDown 1 0.008000000
+Mixed 1 0.008000000
+> roast(y=y,iset1,design,contrast=2,array.weights=c(0.5,1,0.5,1))
+ Active.Prop P.Value
+Down 0 0.99899950
+Up 1 0.00150075
+UpOrDown 1 0.00300000
+Mixed 1 0.00300000
+> w <- matrix(runif(100*4),100,4)
+> roast(y=y,iset1,design,contrast=2,weights=w)
+ Active.Prop P.Value
+Down 0 0.9994997
+Up 1 0.0010005
+UpOrDown 1 0.0020000
+Mixed 1 0.0020000
+> mroast(y=y,list(set1=iset1,set2=iset2),design,contrast=2,gene.weights=runif(100))
+ NGenes PropDown PropUp Direction PValue FDR PValue.Mixed FDR.Mixed
+set1 5 0 1 Up 0.008 0.015 0.008 0.015
+set2 5 0 0 Up 0.959 0.959 0.687 0.687
+> mroast(y=y,list(set1=iset1,set2=iset2),design,contrast=2,array.weights=c(0.5,1,0.5,1))
+ NGenes PropDown PropUp Direction PValue FDR PValue.Mixed FDR.Mixed
+set1 5 0 1 Up 0.004 0.007 0.004 0.007
+set2 5 0 0 Up 0.679 0.679 0.658 0.658
+> mroast(y=y,list(set1=iset1,set2=iset2),design,contrast=2,weights=w)
+ NGenes PropDown PropUp Direction PValue FDR PValue.Mixed FDR.Mixed
+set1 5 0.0 1 Up 0.003 0.005 0.003 0.005
+set2 5 0.2 0 Down 0.950 0.950 0.250 0.250
+> mroast(y=y,list(set1=iset1,set2=iset2),design,contrast=2,weights=w,array.weights=c(0.5,1,0.5,1))
+ 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
+>
+> ### camera
+>
+> camera(y=y,iset1,design,contrast=2,weights=c(0.5,1,0.5,1))
+ NGenes Correlation Direction PValue
+set1 5 -0.2481655 Up 0.001050253
+> camera(y=y,list(set1=iset1,set2=iset2),design,contrast=2)
+ NGenes Correlation Direction PValue FDR
+set1 5 -0.2481655 Up 0.0009047749 0.00180955
+set2 5 0.1719094 Down 0.9068364378 0.90683644
+>
+> ### with EList arg
+>
+> y <- new("EList",list(E=y))
+> roast(y=y,iset1,design,contrast=2)
+ Active.Prop P.Value
+Down 0 0.997498749
+Up 1 0.003001501
+UpOrDown 1 0.006000000
+Mixed 1 0.006000000
+> camera(y=y,iset1,design,contrast=2)
+ NGenes Correlation Direction PValue
+set1 5 -0.2481655 Up 0.0009047749
+>
+> ### eBayes with trend
+>
+> fit <- lmFit(y,design)
+> fit <- eBayes(fit,trend=TRUE)
+> topTable(fit,coef=2)
+ logFC AveExpr t P.Value adj.P.Val B
+3 3.488703 1.03931081 4.860410 0.0002436118 0.01647958 0.6722078
+2 3.729512 1.73488969 4.700998 0.0003295917 0.01647958 0.3777787
+4 2.696676 1.74060725 3.280613 0.0053915597 0.17971866 -2.3313104
+1 2.391846 1.72305203 3.009776 0.0092611288 0.23152822 -2.8478458
+5 2.387967 1.63066783 2.786529 0.0144249169 0.26573834 -3.2671364
+33 -1.492317 -0.07525287 -2.735781 0.0159443006 0.26573834 -3.3613142
+80 -1.839760 -0.32802306 -2.594532 0.0210374835 0.30053548 -3.6207072
+95 -1.907074 1.26297763 -2.462009 0.0272186263 0.33449167 -3.8598265
+39 1.366141 -0.27360750 2.409767 0.0301042507 0.33449167 -3.9527943
+70 -1.789476 0.21771869 -2.184062 0.0462410739 0.46241074 -4.3445901
+> fit$df.prior
+[1] 12.17481
+> fit$s2.prior
+ [1] 0.7108745 0.7186517 0.3976222 0.7224388 0.6531157 0.3014062 0.3169880
+ [8] 0.3149772 0.3074632 0.2917431 0.3329334 0.3378027 0.2900500 0.3031741
+ [15] 0.3221763 0.2981580 0.2897078 0.2925188 0.2924234 0.3042822 0.2923686
+ [22] 0.2897022 0.3251669 0.2929813 0.4922090 0.2902725 0.3018205 0.3029119
+ [29] 0.3030051 0.3331358 0.3259651 0.2939051 0.3077824 0.3553515 0.3139985
+ [36] 0.3181689 0.3197601 0.4687993 0.3316536 0.2897621 0.2910744 0.2907116
+ [43] 0.2907966 0.3265722 0.3240487 0.3241126 0.3003970 0.3064187 0.3645035
+ [50] 0.2994391 0.3295512 0.2901076 0.2898658 0.3086659 0.2897209 0.2982976
+ [57] 0.3043910 0.2900320 0.3006936 0.2935101 0.3646949 0.3596385 0.3064203
[64] 0.3027439 0.3076483 0.3363356 0.3504336 0.3496698 0.2897618 0.2898810
[71] 0.3182290 0.3121707 0.2945001 0.2897549 0.3579410 0.3434376 0.3037970
[78] 0.3201893 0.3048412 0.3394079 0.3516034 0.3034589 0.3120384 0.3007827
@@ -1210,103 +1330,31 @@ TestRe.2900320 0.3006936 0.2935101 0.3646949 0.3596385 0.3064203
+ "1370","137075","1371","137209","1373","137362","1374","137492","1375","1376",
+ "137682","137695","137735","1378","137814","137868","137872","137886","137902","137964")
> go <- goana(fit,FDR=0.8,geneid=EB)
-> topGO(go,n=10)
- Term Ont
-GO:0055082 cellular chemical homeostasis BP
-GO:0006897 endocytosis BP
-GO:0019725 cellular homeostasis BP
-GO:0006915 apoptotic process BP
-GO:0012501 programmed cell death BP
-GO:0008283 cell proliferation BP
-GO:0042127 regulation of cell proliferation BP
-GO:0097190 apoptotic signaling pathway BP
-GO:0006909 phagocytosis BP
-GO:0007169 transmembrane receptor protein tyrosine kinase signaling pathway BP
- N Up Down P.Up P.Down
-GO:0055082 2 0 2 1.000000000 0.004242424
-GO:0006897 3 3 0 0.005046382 1.000000000
-GO:0019725 3 0 2 1.000000000 0.012294372
-GO:0006915 4 3 1 0.017844424 0.255402330
-GO:0012501 4 3 1 0.017844424 0.255402330
-GO:0008283 4 1 2 0.553950615 0.023749721
-GO:0042127 4 1 2 0.553950615 0.023749721
-GO:0097190 2 2 0 0.030909091 1.000000000
-GO:0006909 2 2 0 0.030909091 1.000000000
-GO:0007169 2 2 0 0.030909091 1.000000000
-> topGO(go,n=10,sort="down")
- Term Ont N Up Down P.Up
-GO:0055082 cellular chemical homeostasis BP 2 0 2 1.0000000
-GO:0019725 cellular homeostasis BP 3 0 2 1.0000000
-GO:0008283 cell proliferation BP 4 1 2 0.5539506
-GO:0042127 regulation of cell proliferation BP 4 1 2 0.5539506
-GO:0048878 chemical homeostasis BP 5 1 2 0.6375849
-GO:0032502 developmental process BP 22 3 4 0.8185069
-GO:0048519 negative regulation of biological process BP 14 2 3 0.7685153
-GO:0043230 extracellular organelle CC 14 1 3 0.9502416
-GO:0065010 extracellular membrane-bounded organelle CC 14 1 3 0.9502416
-GO:0070062 extracellular vesicular exosome CC 14 1 3 0.9502416
- P.Down
-GO:0055082 0.004242424
-GO:0019725 0.012294372
-GO:0008283 0.023749721
-GO:0042127 0.023749721
-GO:0048878 0.038228009
-GO:0032502 0.040079017
-GO:0048519 0.055161144
-GO:0043230 0.055161144
-GO:0065010 0.055161144
-GO:0070062 0.055161144
+> 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
+> 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
>
> proc.time()
user system elapsed
- 2.83 0.18 3.04
- limma/vignettes/ 0000755 0001263 0001264 00000000000 12560533400 015242 5 ustar 00biocbuild phs_compbio [...]
-%\VignetteDepends{}
-%\VignetteKeywords{microarray linear model}
-%\VignettePackage{limma}
-\documentclass[12pt]{article}
-
-\textwidth=6.2in
-\textheight=8.5in
-\oddsidemargin=0.2in
-\evensidemargin=0.2in
-\headheight=0in
-\headsep=0in
-
-\begin{document}
-\title{Limma Package Introduction}
-\author{Gordon Smyth}
-\date{23 October 2004, Revised 21 October 2013}
-\maketitle
-
-Limma is an R package for the analysis of gene expression microarray data, especially the use of linear models for analysing designed experiments and the assessment of differential expression.
-Limma provides the ability to analyse comparisons between many RNA targets simultaneously in arbitrary complicated designed experiments.
-Empirical Bayesian methods are used to provide stable results even when the number of arrays is small.
-The normalization and data analysis functions are for two-color spotted microarrays.
-The linear model and differential expression functions apply to all microarray technologies including Affymetrix and other single-channel oligonucleotide platforms.
-
-The full Limma User's Guide is available as part of the online documentation.
-To reach the User's Guide you need to install the limma package.
-If you've installed the package and you're using Windows, type \texttt{library(limma)} at the R prompt then click on ``limma'' from the drop-down menu called ``Vignettes''.
-If you're not using Windows, you can type
-\begin{Schunk}
-\begin{Sinput}
-> library(limma)
-> limmaUsersGuide()
-\end{Sinput}
-\end{Schunk}
-or alternatively
-\begin{Schunk}
-\begin{Sinput}
-> help.start()
-\end{Sinput}
-\end{Schunk}
-and follow the links to the limma package help.
-
-\end{document}
-
-
-
-
-
- [...]
\ No newline at end of file
+ 3.15 0.09 3.24
diff --git a/vignettes/intro.Rnw b/vignettes/intro.Rnw
new file mode 100644
index 0000000..88f952a
--- /dev/null
+++ b/vignettes/intro.Rnw
@@ -0,0 +1,49 @@
+%\VignetteIndexEntry{Limma One Page Introduction}
+%\VignetteDepends{}
+%\VignetteKeywords{microarray linear model}
+%\VignettePackage{limma}
+\documentclass[12pt]{article}
+
+\textwidth=6.2in
+\textheight=8.5in
+\oddsidemargin=0.2in
+\evensidemargin=0.2in
+\headheight=0in
+\headsep=0in
+
+\begin{document}
+\title{Limma Package Introduction}
+\author{Gordon Smyth}
+\date{23 October 2004, Revised 21 October 2013}
+\maketitle
+
+Limma is an R package for the analysis of gene expression microarray data, especially the use of linear models for analysing designed experiments and the assessment of differential expression.
+Limma provides the ability to analyse comparisons between many RNA targets simultaneously in arbitrary complicated designed experiments.
+Empirical Bayesian methods are used to provide stable results even when the number of arrays is small.
+The normalization and data analysis functions are for two-color spotted microarrays.
+The linear model and differential expression functions apply to all microarray technologies including Affymetrix and other single-channel oligonucleotide platforms.
+
+The full Limma User's Guide is available as part of the online documentation.
+To reach the User's Guide you need to install the limma package.
+If you've installed the package and you're using Windows, type \texttt{library(limma)} at the R prompt then click on ``limma'' from the drop-down menu called ``Vignettes''.
+If you're not using Windows, you can type
+\begin{Schunk}
+\begin{Sinput}
+> library(limma)
+> limmaUsersGuide()
+\end{Sinput}
+\end{Schunk}
+or alternatively
+\begin{Schunk}
+\begin{Sinput}
+> help.start()
+\end{Sinput}
+\end{Schunk}
+and follow the links to the limma package help.
+
+\end{document}
+
+
+
+
+
--
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