[med-svn] [r-bioc-limma] 01/04: Imported Upstream version 3.18.10~dfsg
Charles Plessy
plessy at moszumanska.debian.org
Thu Feb 6 10:59:32 UTC 2014
This is an automated email from the git hooks/post-receive script.
plessy pushed a commit to branch master
in repository r-bioc-limma.
commit c3b46aad9685bf85e05715137c62634a0fbd6764
Author: Charles Plessy <plessy at debian.org>
Date: Thu Feb 6 19:51:30 2014 +0900
Imported Upstream version 3.18.10~dfsg
---
DESCRIPTION | 6 +-
R/ebayes.R | 7 +-
R/geneset-roast.R | 115 ++++++++++------
R/lmfit.R | 30 ++--
R/read-ilmn.R | 7 +-
R/subsetting.R | 311 ++++++++++++++++--------------------------
R/toptable.R | 22 ++-
R/treat.R | 106 +-------------
R/venn.R | 8 +-
inst/doc/changelog.txt | 53 +++++++
inst/doc/intro.pdf | Bin 46191 -> 46191 bytes
man/camera.Rd | 5 +-
man/ebayes.Rd | 2 +-
man/normalizebetweenarrays.Rd | 2 +-
man/roast.Rd | 4 +-
man/subsetting.Rd | 31 +++--
man/voom.Rd | 8 +-
tests/limma-Tests.R | 6 +-
tests/limma-Tests.Rout.save | 58 +++++---
19 files changed, 372 insertions(+), 409 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 44d6625..acc3ea4 100755
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
Package: limma
-Version: 3.18.4
-Date: 2013/12/02
+Version: 3.18.10
+Date: 2014/01/31
Title: Linear Models for Microarray Data
Author: Gordon Smyth [cre,aut], Matthew Ritchie [ctb], Jeremy Silver [ctb], James Wettenhall [ctb], Natalie Thorne [ctb], Mette Langaas [ctb], Egil Ferkingstad [ctb], Marcus Davy [ctb], Francois Pepin [ctb], Dongseok Choi [ctb], Davis McCarthy [ctb], Di Wu [ctb], Alicia Oshlack [ctb], Carolyn de Graaf [ctb], Yifang Hu [ctb], Wei Shi [ctb], Belinda Phipson [ctb]
Maintainer: Gordon Smyth <smyth at wehi.edu.au>
@@ -14,4 +14,4 @@ URL: http://bioinf.wehi.edu.au/limma
biocViews: Microarray, OneChannel, TwoChannel, DataImport,
QualityControl, Preprocessing, Bioinformatics,
DifferentialExpression, MultipleComparisons, TimeCourse
-Packaged: 2013-12-03 04:16:44 UTC; biocbuild
+Packaged: 2014-02-01 04:15:50 UTC; biocbuild
diff --git a/R/ebayes.R b/R/ebayes.R
index 7bafa0d..314fd05 100755
--- a/R/ebayes.R
+++ b/R/ebayes.R
@@ -103,7 +103,7 @@ tmixture.matrix <- 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 (v0+v1)/v1*t(df) with probability proportion and t(df) otherwise
+# 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.
@@ -123,7 +123,10 @@ tmixture.vector <- function(tstat,stdev.unscaled,df,proportion,v0.lim=NULL) {
p <- max(ntarget/ngenes,proportion)
tstat <- abs(tstat)
- ttarget <- quantile(tstat,(ngenes-ntarget)/(ngenes-1))
+ if(ngenes>1)
+ ttarget <- quantile(tstat,(ngenes-ntarget)/(ngenes-1))
+ else
+ ttarget <- tstat
top <- (tstat >= ttarget)
tstat <- tstat[top]
v1 <- stdev.unscaled[top]^2
diff --git a/R/geneset-roast.R b/R/geneset-roast.R
index 80fcca5..85e27f0 100644
--- a/R/geneset-roast.R
+++ b/R/geneset-roast.R
@@ -37,7 +37,7 @@ roast.MAList <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.stat
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)
# Rotation gene set testing for linear models
# Gordon Smyth and Di Wu
-# Created 24 Apr 2008. Last modified 19 Sep 2013.
+# Created 24 Apr 2008. Last modified 9 Jan 2014.
{
# Check y
y <- as.matrix(y)
@@ -53,37 +53,44 @@ roast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.sta
design <- as.matrix(design)
if(nrow(design) != n) stop("row dimension of design matrix must match column dimension of data")
p <- ncol(design)
- p0 <- p-1
+ p0 <- p-1L
d <- n-p
-# Check contrast
- if(length(contrast) == 1) {
- u <- rep.int(0,p)
- u[contrast] <- 1
- contrast <- u
- }
- if(length(contrast) != p) stop("length of contrast must match column dimension of design")
- if(all(contrast==0)) stop("contrast all zero")
-
# Check set.statistic
set.statistic <- match.arg(set.statistic,c("mean","floormean","mean50","msq"))
# 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")
-
-# Check and divide out array weights
+
+# Check observation 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")
+ }
+
+# Check array weights
if(!is.null(array.weights)) {
- if(!is.null(weights)) stop("Can't specify both array.weights and observational weights")
- if(any(array.weights <= 0)) stop("array.weights must be positive")
if(length(array.weights) != n) stop("Length of array.weights doesn't match number of array")
- design <- design*sqrt(array.weights)
- y <- t(t(y)*sqrt(array.weights))
+ if(any(array.weights <= 0)) stop("array.weights must be positive")
+ if(!is.null(weights)) {
+ weights <- .matvec(weights,array.weights)
+ array.weights <- NULL
+ }
}
-# Check and divide out block correlation
+# 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)) {
- if(!is.null(weights)) stop("Can't use block with weights")
block <- as.vector(block)
if (length(block) != n) stop("Length of block does not match number of arrays")
ub <- unique(block)
@@ -96,21 +103,21 @@ roast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.sta
design <- backsolve(R, design, transpose = TRUE)
}
-# Check 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")
+# Check contrast
+ if(all(contrast==0)) stop("contrast all zero")
+
+# Reform design matrix so that contrast is last coefficient
+ if(length(contrast) == 1) {
+ 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
}
-# Reform design matrix so that contrast of interest is last column
-# qr <- qr(contrast)
-# Q <- qr.Q(qr,complete=TRUE)
-# sign1 <- sign(qr$qr[1,1])
-# Q <- cbind(Q[,-1],Q[,1])
-# X <- design %*% Q
- X <- contrastAsCoef(design, contrast, first=FALSE)$design
qr <- qr(X)
signc <- sign(qr$qr[p,p])
@@ -327,10 +334,12 @@ mroast.MAList <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.sta
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,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 19 Sep 2013.
+# Created 28 Jan 2010. Last revised 9 Jan 2014.
{
# Check y
y <- as.matrix(y)
+ ngenes <- nrow(y)
+ n <- ncol(y)
# Check index
if(is.null(index)) index <- rep(TRUE,nrow(y))
@@ -343,24 +352,50 @@ mroast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.st
design <- as.matrix(design)
# Check gene.weights
- if(!is.null(gene.weights)) if(length(gene.weights) != nrow(y)) stop("gene.weights must have length equal to nrow(y)")
+ 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) != ncol(y)) stop("array.weights wrong length")
+ if(any(array.weights<=0)) stop("array.weights must be positive")
+ }
# Check weights
if(!is.null(weights)) {
- if(!is.null(array.weights)) stop("Can't specify both array weights and observational weights")
weights <- as.matrix(weights)
if(any(dim(weights) != dim(y))) stop("weights must have same dimensions as y")
+ if(!is.null(array.weights)) {
+ weights <- .matvec(weights,array.weights)
+ array.weights <- NULL
+ }
}
-# Check array.weights
+# Divide out array.weights
if(!is.null(array.weights)) {
- if(length(array.weights) != ncol(y)) stop("array.weights wrong length")
- weights <- 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,block=block,correlation=correlation)
+ fit <- lmFit(y,design=design,weights=weights)
if(trend.var) {
covariate <- fit$Amean
if(is.null(covariate)) covariate <- rowMeans(y)
@@ -376,7 +411,7 @@ mroast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.st
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,array.weights=array.weights,weights=weights,block=block,correlation=correlation,var.prior=var.prior,df.prior=df.prior,nrot=nrot)
+ 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)
pv[i,] <- out$p.value$P.Value
active[i,] <- out$p.value$Active.Prop
NGenes[i] <- out$ngenes.in.set
diff --git a/R/lmfit.R b/R/lmfit.R
index 1ea3104..aecb915 100644
--- a/R/lmfit.R
+++ b/R/lmfit.R
@@ -110,9 +110,12 @@ lm.series <- function(M,design=NULL,ndups=1,spacing=1,weights=NULL)
fit <- lm.wfit(design, t(M), weights[1,])
fit$weights <- NULL
}
- if(fit$df.residual>0)
- fit$sigma <- sqrt(colMeans(fit$effects[(fit$rank + 1):narrays,,drop=FALSE]^2))
- else
+ if(fit$df.residual>0) {
+ if(is.matrix(fit$effects))
+ fit$sigma <- sqrt(colMeans(fit$effects[(fit$rank + 1):narrays,,drop=FALSE]^2))
+ else
+ fit$sigma <- sqrt(mean(fit$effects[(fit$rank + 1):narrays]^2))
+ } else
fit$sigma <- rep(NA,ngenes)
fit$fitted.values <- fit$residuals <- fit$effects <- NULL
fit$coefficients <- t(fit$coefficients)
@@ -147,11 +150,7 @@ lm.series <- function(M,design=NULL,ndups=1,spacing=1,weights=NULL)
beta[i,] <- out$coef
stdev.unscaled[i,est] <- sqrt(diag(chol2inv(out$qr$qr,size=out$rank)))
df.residual[i] <- out$df.residual
- if(df.residual[i] > 0)
- if(is.null(weights))
- sigma[i] <- sqrt(sum(out$residuals^2)/out$df.residual)
- else
- sigma[i] <- sqrt(sum(w*out$residuals^2)/out$df.residual)
+ if(df.residual[i] > 0) sigma[i] <- sqrt(mean(out$effects[-(1:out$rank)]^2))
}
}
@@ -161,7 +160,7 @@ lm.series <- function(M,design=NULL,ndups=1,spacing=1,weights=NULL)
est <- QR$pivot[1:QR$rank]
dimnames(cov.coef) <- list(coef.names[est],coef.names[est])
- list(coefficients=drop(beta),stdev.unscaled=drop(stdev.unscaled),sigma=sigma,df.residual=df.residual,cov.coefficients=cov.coef,pivot=QR$pivot)
+ list(coefficients=beta,stdev.unscaled=stdev.unscaled,sigma=sigma,df.residual=df.residual,cov.coefficients=cov.coef,pivot=QR$pivot)
}
mrlm <- function(M,design=NULL,ndups=1,spacing=1,weights=NULL,...)
@@ -211,7 +210,7 @@ mrlm <- function(M,design=NULL,ndups=1,spacing=1,weights=NULL,...)
cov.coef <- chol2inv(QR$qr,size=QR$rank)
est <- QR$pivot[1:QR$rank]
dimnames(cov.coef) <- list(coef.names[est],coef.names[est])
- list(coefficients=drop(beta),stdev.unscaled=drop(stdev.unscaled),sigma=sigma,df.residual=df.residual,cov.coefficients=cov.coef,pivot=QR$pivot)
+ list(coefficients=beta,stdev.unscaled=stdev.unscaled,sigma=sigma,df.residual=df.residual,cov.coefficients=cov.coef,pivot=QR$pivot)
}
gls.series <- function(M,design=NULL,ndups=2,spacing=1,block=NULL,correlation=NULL,weights=NULL,...)
@@ -276,9 +275,12 @@ gls.series <- function(M,design=NULL,ndups=2,spacing=1,block=NULL,correlation=NU
X <- backsolve(cholV,design,transpose=TRUE)
dimnames(X) <- dimnames(design)
fit <- lm.fit(X,y)
- if(fit$df.residual>0)
- fit$sigma <- sqrt(colMeans(fit$effects[-(1:fit$rank),,drop=FALSE]^2))
- else
+ if(fit$df.residual>0) {
+ if(is.matrix(fit$effects))
+ fit$sigma <- sqrt(colMeans(fit$effects[-(1:fit$rank),,drop=FALSE]^2))
+ else
+ fit$sigma <- sqrt(mean(fit$effects[-(1:fit$rank)]^2))
+ } else
fit$sigma <- rep(NA,ngenes)
fit$fitted.values <- fit$residuals <- fit$effects <- NULL
fit$coefficients <- t(fit$coefficients)
@@ -334,7 +336,7 @@ gls.series <- function(M,design=NULL,ndups=2,spacing=1,block=NULL,correlation=NU
cov.coef <- chol2inv(QR$qr,size=QR$rank)
est <- QR$pivot[1:QR$rank]
dimnames(cov.coef) <- list(coef.names[est],coef.names[est])
- list(coefficients=drop(beta),stdev.unscaled=drop(stdev.unscaled),sigma=sigma,df.residual=df.residual,ndups=ndups,spacing=spacing,block=block,correlation=correlation,cov.coefficients=cov.coef,pivot=QR$pivot)
+ list(coefficients=beta,stdev.unscaled=stdev.unscaled,sigma=sigma,df.residual=df.residual,ndups=ndups,spacing=spacing,block=block,correlation=correlation,cov.coefficients=cov.coef,pivot=QR$pivot)
}
is.fullrank <- function(x)
diff --git a/R/read-ilmn.R b/R/read-ilmn.R
index 63ed5ee..952f317 100644
--- a/R/read-ilmn.R
+++ b/R/read-ilmn.R
@@ -3,7 +3,7 @@
read.ilmn <- function(files=NULL, ctrlfiles=NULL, path=NULL, ctrlpath=NULL, probeid="Probe", annotation=c("TargetID", "SYMBOL"), expr="AVG_Signal", other.columns="Detection",sep="\t", quote="\"", verbose=TRUE, ...)
# Read one or more files of Illumina BeadStudio output
# Wei Shi and Gordon Smyth.
-# Created 15 July 2009. Last modified 27 November 2013.
+# Created 15 July 2009. Last modified 24 January 2014.
{
if(!is.null(files)){
f <- unique(files)
@@ -99,7 +99,10 @@ read.ilmn.targets <- function(targets, ...)
rownames(elist$E) <- pids
# Add probe annotation
- if(length(anncol)) elist$genes <- x[,anncol,drop=FALSE]
+ if(length(anncol)) {
+ elist$genes <- x[,anncol,drop=FALSE]
+ row.names(elist$genes) <- pids
+ }
# elist$targets <- data.frame(SampleNames=snames, stringsAsFactors=FALSE)
diff --git a/R/subsetting.R b/R/subsetting.R
index 4cb97a8..7a6bb3c 100755
--- a/R/subsetting.R
+++ b/R/subsetting.R
@@ -1,232 +1,151 @@
# SUBSET DATA SETS
-assign("[.RGList",
-function(object, i, j, ...) {
-# Subsetting for RGList objects
-# Gordon Smyth
-# 29 June 2003. Last modified 22 December 2005.
+subsetListOfArrays <- function(object,i,j,IJ,IX,I,JX)
+# Subsetting for list-like data objects
+# Gordon Smyth
+# 11 Dec 2013
+{
+# object,IJ,IX,I,JX are required arguments
- if(nargs() != 3) stop("Two subscripts required",call.=FALSE)
- oc <- names(object$other)
- if(missing(i))
- if(missing(j))
- return(object)
- else {
- object$R <- object$R[,j,drop=FALSE]
- object$G <- object$G[,j,drop=FALSE]
- object$Rb <- object$Rb[,j,drop=FALSE]
- object$Gb <- object$Gb[,j,drop=FALSE]
- object$weights <- object$weights[,j,drop=FALSE]
- object$targets <- object$targets[j,,drop=FALSE]
- for(k in oc) object$other[[k]] <- object$other[[k]][,j,drop=FALSE]
- }
- else {
+# Remove components that are absent or scalars
+ len <- vapply(object,length,0)
+ I <- intersect(I,names(len)[len>1L])
+
+ if(missing(i)) {
+ IX <- I <- character(0)
+ if(missing(j)) IJ <- character(0)
+ } else {
if(is.character(i)) {
- i <- match(i,rownames(object))
- i <- i[!is.na(i)]
+ i <- match(i, rownames(object))
+ if(any(is.na(i))) stop("Subscript not found in rownames")
}
- if(missing(j)) {
- object$R <- object$R[i,,drop=FALSE]
- object$G <- object$G[i,,drop=FALSE]
- object$Rb <- object$Rb[i,,drop=FALSE]
- object$Gb <- object$Gb[i,,drop=FALSE]
- object$weights <- object$weights[i,,drop=FALSE]
- object$genes <- object$genes[i,,drop=FALSE]
- for(k in oc) object$other[[k]] <- object$other[[k]][i,,drop=FALSE]
- } else {
- object$R <- object$R[i,j,drop=FALSE]
- object$G <- object$G[i,j,drop=FALSE]
- object$Rb <- object$Rb[i,j,drop=FALSE]
- object$Gb <- object$Gb[i,j,drop=FALSE]
- object$weights <- object$weights[i,j,drop=FALSE]
- object$genes <- object$genes[i,,drop=FALSE]
- object$targets <- object$targets[j,,drop=FALSE]
- for(k in oc) object$other[[k]] <- object$other[[k]][i,j,drop=FALSE]
+ }
+ if(missing(j)) {
+ JX <- character(0)
+ } else {
+ if(is.character(j)) {
+ j <- match(j, colnames(object))
+ if(any(is.na(j))) stop("Subscript not found in colnames")
}
}
+
+ for(a in IJ) object[[a]] <- object[[a]][i,j,drop=FALSE]
+ for(a in IX) object[[a]] <- object[[a]][i, ,drop=FALSE]
+ for(a in I ) object[[a]] <- object[[a]][i]
+ for(a in JX) object[[a]] <- object[[a]][j, ,drop=FALSE]
+
+ object
+}
+
+assign("[.RGList",
+function(object, i, j)
+# Subsetting for RGList objects
+# Gordon Smyth
+# 29 June 2003. Last modified 11 Dec 2013.
+{
+# Recognized components
+ IJ <- c("R","G","Rb","Gb","weights")
+ IX <- c("genes")
+ JX <- c("targets")
+ I <- character(0)
+ object <- subsetListOfArrays(object,i,j,IJ=IJ,IX=IX,I=I,JX=JX)
+
+ oc <- names(object$other)
+ if(!missing(i) || !missing(j)) for(a in oc) object$other[[a]] <- object$other[[a]][i,j,drop=FALSE]
+
object
})
assign("[.MAList",
-function(object, i, j, ...) {
+function(object, i, j)
# Subsetting for MAList objects
# Gordon Smyth
# 29 June 2003. Last modified 22 Dec 2005.
+{
+# Recognized components
+ IJ <- c("M","A","weights")
+ IX <- c("genes")
+ JX <- c("targets","design")
+ I <- character(0)
+ object <- subsetListOfArrays(object,i,j,IJ=IJ,IX=IX,I=I,JX=JX)
- if(nargs() != 3) stop("Two subscripts required",call.=FALSE)
- other <- names(object$other)
- if(missing(i))
- if(missing(j))
- return(object)
- else {
- object$M <- object$M[,j,drop=FALSE]
- object$A <- object$A[,j,drop=FALSE]
- object$weights <- object$weights[,j,drop=FALSE]
- object$targets <- object$targets[j,,drop=FALSE]
- if(!is.null(object$design)) {
- object$design <- as.matrix(object$design)[j,,drop=FALSE]
- if(!is.fullrank(object$design)) warning("subsetted design matrix is singular",call.=FALSE)
- }
- for(a in other) object$other[[a]] <- object$other[[a]][,j,drop=FALSE]
- }
- else {
- if(is.character(i)) {
- i <- match(i,rownames(object))
- i <- i[!is.na(i)]
- }
- if(missing(j)) {
- object$M <- object$M[i,,drop=FALSE]
- object$A <- object$A[i,,drop=FALSE]
- object$weights <- object$weights[i,,drop=FALSE]
- object$genes <- object$genes[i,,drop=FALSE]
- for(a in other) object$other[[a]] <- object$other[[a]][i,,drop=FALSE]
- } else {
- object$M <- object$M[i,j,drop=FALSE]
- object$A <- object$A[i,j,drop=FALSE]
- object$weights <- object$weights[i,j,drop=FALSE]
- object$genes <- object$genes[i,,drop=FALSE]
- object$targets <- object$targets[j,,drop=FALSE]
- if(!is.null(object$design)) {
- object$design <- as.matrix(object$design)[j,,drop=FALSE]
- if(!is.fullrank(object$design)) warning("subsetted design matrix is singular",call.=FALSE)
- }
- for(a in other) object$other[[a]] <- object$other[[a]][i,j,drop=FALSE]
- }
- }
+ if(!missing(j) && !is.null(object$design) && !is.fullrank(object$design)) warning("subsetted design matrix is singular",call.=FALSE)
+
+ oc <- names(object$other)
+ if(!missing(i) || !missing(j)) for(a in oc) object$other[[a]] <- object$other[[a]][i,j,drop=FALSE]
+
object
})
assign("[.EList",
-function(object, i, j, ...) {
+function(object, i, j)
# Subsetting for EList objects
# Gordon Smyth
-# 23 February 2009. Last modified 21 October 2010.
-
- if(nargs() != 3) stop("Two subscripts required",call.=FALSE)
- other <- names(object$other)
- if(missing(i))
- if(missing(j))
- return(object)
- else {
- object$E <- object$E[,j,drop=FALSE]
- object$Eb <- object$Eb[,j,drop=FALSE]
- object$weights <- object$weights[,j,drop=FALSE]
- object$targets <- object$targets[j,,drop=FALSE]
- if(!is.null(object$design)) {
- object$design <- as.matrix(object$design)[j,,drop=FALSE]
- if(!is.fullrank(object$design)) warning("subsetted design matrix is singular",call.=FALSE)
- }
- for(a in other) object$other[[a]] <- object$other[[a]][,j,drop=FALSE]
- }
- else {
- if(is.character(i)) {
- i <- match(i,rownames(object))
- i <- i[!is.na(i)]
- }
- if(missing(j)) {
- object$E <- object$E[i,,drop=FALSE]
- object$Eb <- object$Eb[i,,drop=FALSE]
- object$weights <- object$weights[i,,drop=FALSE]
- object$genes <- object$genes[i,,drop=FALSE]
- for(a in other) object$other[[a]] <- object$other[[a]][i,,drop=FALSE]
- } else {
- object$E <- object$E[i,j,drop=FALSE]
- object$Eb <- object$Eb[i,j,drop=FALSE]
- object$weights <- object$weights[i,j,drop=FALSE]
- object$genes <- object$genes[i,,drop=FALSE]
- object$targets <- object$targets[j,,drop=FALSE]
- if(!is.null(object$design)) {
- object$design <- as.matrix(object$design)[j,,drop=FALSE]
- if(!is.fullrank(object$design)) warning("subsetted design matrix is singular",call.=FALSE)
- }
- for(a in other) object$other[[a]] <- object$other[[a]][i,j,drop=FALSE]
- }
- }
+# 23 February 2009. Last modified 11 Dec 2013.
+{
+# Recognized components
+ IJ <- c("E","Eb","weights")
+ IX <- c("genes")
+ JX <- c("targets","design")
+ I <- character(0)
+ object <- subsetListOfArrays(object,i,j,IJ=IJ,IX=IX,I=I,JX=JX)
+
+ oc <- names(object$other)
+ if(!missing(i) || !missing(j)) for(a in oc) object$other[[a]] <- object$other[[a]][i,j,drop=FALSE]
+
object
})
assign("[.EListRaw", get("[.EList"))
assign("[.MArrayLM",
-function(object, i, j, ...)
+function(object, i, j)
# Subsetting for MArrayLM objects
# Gordon Smyth
-# 26 April 2005. Last modified 28 September 2013.
+# 26 April 2005. Last modified 11 Dec 2013.
{
- if(nargs() != 3) stop("Two subscripts required",call.=FALSE)
- if(!is.null(object$coefficients)) object$coefficients <- as.matrix(object$coefficients)
- if(!is.null(object$stdev.unscaled)) object$stdev.unscaled <- as.matrix(object$stdev.unscaled)
- if(!is.null(object$weights)) object$weights <- as.matrix(object$weights)
- if(!is.null(object$p.value)) object$p.value <- as.matrix(object$p.value)
- if(!is.null(object$lods)) object$lods <- as.matrix(object$lods)
- if(!is.null(object$targets)) object$targets <- as.data.frame(object$targets)
- if(!is.null(object$cov.coefficients)) object$cov.coefficients <- as.matrix(object$cov.coefficients)
- if(!is.null(object$contrasts)) object$contrasts <- as.matrix(object$contrasts)
- if(is.null(object$contrasts) && !is.null(object$coefficients)) {
+# Recognized components
+ IJ <- c("coefficients","stdev.unscaled","t","p.value","lods","weights")
+ IX <- "genes"
+ JX <- c("targets","design")
+ I <- c("Amean","sigma","df.residual","df.prior","df.total","s2.post","F","F.p.value")
+ JJ <- "cov.coefficients"
+ XJ <- "contrasts"
+ J <- "var.prior"
+
+# After subsetting by columns, a contrast component should always be present
+# so that output is equivalent to that from contrasts.fit()
+ if(!missing(j) && is.null(object$contrasts) && !is.null(object$coefficients)) {
object$contrasts <- diag(ncol(object$coefficients))
- rownames(object$contrasts) <- colnames(object$contrasts) <- colnames(object$coefficients)
+ cn <- colnames(object$contrasts)
+ dimnames(object$contrasts) <- list(Coefficient=cn,Contrast=cn)
}
- if(!is.null(object$genes)) object$genes <- as.data.frame(object$genes)
- if(missing(i)) {
- if(missing(j))
- return(object)
- else {
- object$coefficients <- object$coefficients[,j,drop=FALSE]
- object$stdev.unscaled <- object$stdev.unscaled[,j,drop=FALSE]
- object$t <- object$t[,j,drop=FALSE]
- object$weights <- object$weights[,j,drop=FALSE]
- object$p.value <- object$p.value[,j,drop=FALSE]
- object$lods <- object$lods[,j,drop=FALSE]
- object$cov.coefficients <- object$cov.coefficients[j,j,drop=FALSE]
- object$contrasts <- object$contrasts[,j,drop=FALSE]
- object$var.prior <- object$var.prior[j]
- }
- } else {
- if(is.character(i)) {
- i <- match(i,rownames(object))
- i <- i[!is.na(i)]
- }
- if(missing(j)) {
- object$coefficients <- object$coefficients[i,,drop=FALSE]
- object$stdev.unscaled <- object$stdev.unscaled[i,,drop=FALSE]
- object$t <- object$t[i,,drop=FALSE]
- object$weights <- object$weights[i,,drop=FALSE]
- object$p.value <- object$p.value[i,,drop=FALSE]
- object$lods <- object$lods[i,,drop=FALSE]
- object$genes <- object$genes[i,,drop=FALSE]
- } else {
- object$coefficients <- object$coefficients[i,j,drop=FALSE]
- object$stdev.unscaled <- object$stdev.unscaled[i,j,drop=FALSE]
- object$t <- object$t[i,j,drop=FALSE]
- object$weights <- object$weights[i,j,drop=FALSE]
- object$p.value <- object$p.value[i,j,drop=FALSE]
- object$lods <- object$lods[i,j,drop=FALSE]
- object$genes <- object$genes[i,,drop=FALSE]
- object$cov.coefficients <- object$cov.coefficients[j,j,drop=FALSE]
- object$contrasts <- object$contrasts[,j,drop=FALSE]
- object$var.prior <- object$var.prior[j]
- }
- object$df.residual <- object$df.residual[i]
- if(length(object$df.prior)>1) object$df.prior <- object$df.prior[i]
- object$df.total <- object$df.total[i]
- object$sigma <- object$sigma[i]
- object$s2.post <- object$s2.post[i]
- object$Amean <- object$Amean[i]
+
+# Ensure matrix or data.frame objects not dropped to vectors
+ for (a in c(IJ,JJ,XJ)) if(!is.null(object[[a]])) object[[a]] <- as.matrix(object[[a]])
+ for (a in c("targets","genes")) if(!is.null(object[[a]]) && is.null(dim(object[[a]]))) object[[a]] <- data.frame(object[[a]])
+
+ object <- subsetListOfArrays(object,i,j,IJ=IJ,IX=IX,I=I,JX=JX)
+
+# Special treatment for JJ,XJ,J
+ if(!missing(j)) {
+ object$cov.coefficients <- object$cov.coefficients[j,j,drop=FALSE]
+ object$contrasts <- object$contrasts[,j,drop=FALSE]
+ object$var.prior <- object$var.prior[j]
}
- if(!is.null(object$F))
- if(missing(j)) {
- object$F <- object$F[i]
- object$F.p.value <- object$F.p.value[i]
- } else {
- F.stat <- classifyTestsF(object,fstat.only=TRUE)
- object$F <- as.vector(F.stat)
- df1 <- attr(F.stat,"df1")
- df2 <- attr(F.stat,"df2")
- if (df2[1] > 1e+06)
- object$F.p.value <- pchisq(df1*object$F,df1,lower.tail=FALSE)
- else
- object$F.p.value <- pf(object$F,df1,df2,lower.tail=FALSE)
- }
+
+# If columns have been subsetted, need to re-generate F
+ if(!is.null(object[["F"]]) && !missing(j)) {
+ F.stat <- classifyTestsF(object,fstat.only=TRUE)
+ object$F <- as.vector(F.stat)
+ df1 <- attr(F.stat,"df1")
+ df2 <- attr(F.stat,"df2")
+ if (df2[1] > 1e6)
+ object$F.p.value <- pchisq(df1*object$F,df1,lower.tail=FALSE)
+ else
+ object$F.p.value <- pf(object$F,df1,df2,lower.tail=FALSE)
+ }
+
object
})
diff --git a/R/toptable.R b/R/toptable.R
index 0b57ce4..5998e3d 100644
--- a/R/toptable.R
+++ b/R/toptable.R
@@ -3,12 +3,12 @@
topTable <- function(fit,coef=NULL,number=10,genelist=fit$genes,adjust.method="BH",sort.by="B",resort.by=NULL,p.value=1,lfc=0,confint=FALSE)
# Summary table of top genes, object-orientated version
# Gordon Smyth
-# 4 August 2003. Last modified 7 April 2013.
+# 4 August 2003. Last modified 7 Dec 2013.
{
# Check fit
if(!is(fit,"MArrayLM")) stop("fit must be an MArrayLM object")
if(is.null(fit$coefficients)) stop("coefficients not found in fit object")
- if(is.null(fit$stdev.unscaled)) stop("stdev.unscaled not found in fit object")
+ if(confint && is.null(fit$stdev.unscaled)) stop("stdev.unscaled not found in fit object")
if(is.null(coef)) coef <- 1:ncol(fit)
if(length(coef)>1) {
@@ -120,7 +120,7 @@ topTableF <- function(fit,number=10,genelist=fit$genes,adjust.method="BH",sort.b
toptable <- function(fit,coef=1,number=10,genelist=NULL,A=NULL,eb=NULL,adjust.method="BH",sort.by="B",resort.by=NULL,p.value=1,lfc=0,confint=FALSE,...)
# Summary table of top genes
# Gordon Smyth
-# 21 Nov 2002. Last revised 16 April 2013.
+# 21 Nov 2002. Last revised 7 Dec 2013.
{
# Check fit
fit$coefficients <- as.matrix(fit$coefficients)
@@ -181,11 +181,20 @@ toptable <- function(fit,coef=1,number=10,genelist=NULL,A=NULL,eb=NULL,adjust.me
coef <- 1
}
+# Check for lods compoent
+ if(is.null(eb$lods)) {
+ if(sort.by=="B") stop("Trying to sort.by B, but B-statistic (lods) not found in MArrayLM object",.call=FALSE)
+ if(!is.null(resort.by)) if(resort.by=="B") stop("Trying to resort.by B, but B-statistic (lods) not found in MArrayLM object",.call=FALSE)
+ include.B <- FALSE
+ } else {
+ include.B <- TRUE
+ }
+
# Extract statistics for table
M <- fit$coefficients[,coef]
tstat <- as.matrix(eb$t)[,coef]
P.Value <- as.matrix(eb$p.value)[,coef]
- B <- as.matrix(eb$lods)[,coef]
+ if(include.B) B <- as.matrix(eb$lods)[,coef]
# Apply multiple testing adjustment
adj.P.Value <- p.adjust(P.Value,method=adjust.method)
@@ -201,7 +210,7 @@ toptable <- function(fit,coef=1,number=10,genelist=NULL,A=NULL,eb=NULL,adjust.me
tstat <- tstat[sig]
P.Value <- P.Value[sig]
adj.P.Value <- adj.P.Value[sig]
- B <- B[sig]
+ if(include.B) B <- B[sig]
rn <- rn[sig]
}
@@ -232,7 +241,8 @@ toptable <- function(fit,coef=1,number=10,genelist=NULL,A=NULL,eb=NULL,adjust.me
tab$CI.975 <- M[top]+margin.error
}
if(!is.null(A)) tab$AveExpr <- A[top]
- tab <- data.frame(tab,t=tstat[top],P.Value=P.Value[top],adj.P.Val=adj.P.Value[top],B=B[top])
+ tab <- data.frame(tab,t=tstat[top],P.Value=P.Value[top],adj.P.Val=adj.P.Value[top])
+ if(include.B) tab$B <- B[top]
rownames(tab) <- rn[top]
# Resort table
diff --git a/R/treat.R b/R/treat.R
index ebaf34e..8ffd95a 100644
--- a/R/treat.R
+++ b/R/treat.R
@@ -1,6 +1,6 @@
### treat.R
-treat <- function(fit, lfc=0, trend=FALSE)
+treat <- function(fit, lfc=0, trend=FALSE, robust=FALSE, winsor.tail.p=c(0.05,0.1))
# Moderated t-statistics with threshold
# Davis McCarthy, Gordon Smyth
# 25 July 2008. Last revised 7 April 2013.
@@ -27,7 +27,7 @@ treat <- function(fit, lfc=0, trend=FALSE)
} else {
covariate <- NULL
}
- sv <- squeezeVar(sigma^2, df.residual, covariate=covariate)
+ sv <- squeezeVar(sigma^2, df.residual, covariate=covariate, robust=robust, winsor.tail.p=winsor.tail.p)
fit$df.prior <- sv$df.prior
fit$s2.prior <- sv$var.prior
fit$s2.post <- sv$var.post
@@ -55,110 +55,16 @@ topTreat <- function(fit,coef=1,number=10,genelist=fit$genes,adjust.method="BH",
# Gordon Smyth
# 15 June 2009. Last modified 20 April 2013.
{
-# Check fit object
- M <- as.matrix(fit$coefficients)
- rn <- rownames(M)
- A <- fit$Amean
- if(is.null(A)) {
- if(sort.by=="A") stop("Cannot sort by A-values as these have not been given")
- } else {
- if(NCOL(A)>1) A <- rowMeans(A,na.rm=TRUE)
- }
-
# Check coef is length 1
if(length(coef)>1) {
coef <- coef[1]
warning("Treat is for single coefficients: only first value of coef being used")
}
-# Ensure genelist is a data.frame
- if(!is.null(genelist) && is.null(dim(genelist))) genelist <- data.frame(ID=genelist,stringsAsFactors=FALSE)
-
-# Check rownames
- if(is.null(rn))
- rn <- 1:nrow(M)
- else
- if(anyDuplicated(rn)) {
- if(is.null(genelist))
- genelist <- data.frame(ID=rn,stringsAsFactors=FALSE)
- else
- if("ID" %in% names(genelist))
- genelist$ID0 <- rn
- else
- genelist$ID <- rn
- rn <- 1:nrow(M)
- }
-
-# Check sort.by
- sort.by <- match.arg(sort.by,c("logFC","M","A","Amean","AveExpr","P","p","T","t","none"))
- if(sort.by=="M") sort.by="logFC"
- if(sort.by=="A" || sort.by=="Amean") sort.by <- "AveExpr"
- if(sort.by=="T") sort.by <- "t"
- if(sort.by=="p") sort.by <- "P"
-
-# Check resort.by
- if(!is.null(resort.by)) {
- resort.by <- match.arg(resort.by,c("logFC","M","A","Amean","AveExpr","P","p","T","t"))
- if(resort.by=="M") resort.by <- "logFC"
- if(resort.by=="A" || resort.by=="Amean") resort.by <- "AveExpr"
- if(resort.by=="p") resort.by <- "P"
- if(resort.by=="T") resort.by <- "t"
- }
-
-# Extract columns from fit
- M <- M[,coef]
- tstat <- as.matrix(fit$t)[,coef]
- P.Value <- as.matrix(fit$p.value)[,coef]
-
-# Apply multiple testing adjustment
- adj.P.Value <- p.adjust(P.Value,method=adjust.method)
-
-# Apply p.value threshold
- if(p.value < 1) {
- sig <- (adj.P.Value < p.value)
- if(any(is.na(sig))) sig[is.na(sig)] <- FALSE
- if(all(!sig)) return(data.frame())
- genelist <- genelist[sig,,drop=FALSE]
- M <- M[sig]
- A <- A[sig]
- tstat <- tstat[sig]
- P.Value <- P.Value[sig]
- adj.P.Value <- adj.P.Value[sig]
- }
- ord <- switch(sort.by,
- logFC=order(abs(M),decreasing=TRUE),
- AveExpr=order(A,decreasing=TRUE),
- P=order(P.Value,decreasing=FALSE),
- t=order(abs(tstat),decreasing=TRUE),
- none=1:length(M)
- )
-
-# Enough rows left?
- if(length(M) < number) number <- length(M)
- if(number < 1) return(data.frame())
-
-# Assemble data.frame of top genes
- top <- ord[1:number]
- if(is.null(genelist))
- tab <- data.frame(logFC=M[top])
- else {
- tab <- data.frame(genelist[top,,drop=FALSE],logFC=M[top],stringsAsFactors=FALSE)
- }
- if(!is.null(A)) tab <- data.frame(tab,AveExpr=A[top])
- tab <- data.frame(tab,t=tstat[top],P.Value=P.Value[top],adj.P.Val=adj.P.Value[top])
- rownames(tab) <- rn[top]
-
-# Resort
- if(!is.null(resort.by)) {
- ord <- switch(resort.by,
- logFC=order(tab$logFC,decreasing=TRUE),
- AveExpr=order(tab$AveExpr,decreasing=TRUE),
- P=order(tab$P.Value,decreasing=FALSE),
- t=order(tab$t,decreasing=TRUE)
- )
- tab <- tab[ord,]
- }
+# Check sort.by and resort.by
+ if(sort.by=="B") stop("Trying to sort.by B, but treat doesn't produce a B-statistic")
+ if(!is.null(resort.by)) if(resort.by=="B") stop("Trying to resort.by B, but treat doesn't produce a B-statistic")
- tab
+ topTable(fit=fit,coef=coef,number=number,genelist=genelist,adjust.method=adjust.method,sort.by=sort.by,resort.by=sort.by,p.value=p.value)
}
diff --git a/R/venn.R b/R/venn.R
index 5453ee8..18dd6a2 100755
--- a/R/venn.R
+++ b/R/venn.R
@@ -30,9 +30,9 @@ vennCounts <- function(x,include="both") {
vennDiagram <- function(object,include="both",names=NULL,mar=rep(1,4),cex=c(1.5,1,0.7),lwd=1,circle.col=NULL,counts.col=NULL,show.include=NULL,...)
# Plot Venn diagram
-# Gordon Smyth, James Wettenhall.
-# Capabilities for multiple counts and colors by Francois Pepin.
-# 4 July 2003. Last modified 10 March 2013.
+# Gordon Smyth, James Wettenhall, Yifang Hu.
+# Capabilities for multiple counts and colors uses code by Francois Pepin.
+# 4 July 2003. Last modified 31 January 2014.
{
# Check include
include <- as.character(include)
@@ -140,7 +140,7 @@ vennDiagram <- function(object,include="both",names=NULL,mar=rep(1,4),cex=c(1.5,
##############################################
# Open plot
- plot(c(-20, 420), c(-20, 420), type="n", axes=FALSE, ylab="", xlab="")
+ plot(c(-20, 420), c(-20, 420), type="n", axes=FALSE, ylab="", xlab="", ...)
# Function to turn and move ellipses
relocate_elp <- function(e, alpha, x, y)
diff --git a/inst/doc/changelog.txt b/inst/doc/changelog.txt
index cfa5cd9..b352173 100755
--- a/inst/doc/changelog.txt
+++ b/inst/doc/changelog.txt
@@ -1,3 +1,56 @@
+31 January 2014: limma 3.18.10
+
+- treat() now has new arguments robust and winsor.tail.p which are
+ passed through to robust empirical Bayes estimation.
+
+- vennDiagram() wasn't passing extra arguments (...) to plot() when
+ the number of sets was greater than 3. Now fixed.
+
+- read.ilmn() now sets the same probe IDs as rownames for both the
+ expression matrix E and the annotation data.frame genes.
+
+- romer() now uses propTrueNull(method="lfdr") instead of convest().
+ This makes romer() substantially faster when the number of genes
+ is large.
+
+11 Jan 2014: limma 3.18.9
+
+- update citation of Law et al (2014) reference for voom().
+
+11 Jan 2014: limma 3.18.8
+
+- roast() and mroast() now permit array weights and observation
+ weights to both be specified.
+
+- bug fix to mroast() which was not accepting array.weights.
+
+17 Dec 2013: limma 3.18.7
+
+- make sure that F is recomputed when the columns of an MArrayLM
+ are subsetted.
+
+12 Dec 2013: limma 3.18.6
+
+- linear model fit functions lm.series(), mrlm.series() and
+ gls.series() no longer drop the dimensions of the components of
+ the fitted object when there is just coefficient or just one gene.
+ Previously this was done inconsistently in some cases but not
+ others. Now the matrix components always keep dimensions.
+
+- New function subsetListOfArrays(), which is used to simply the
+ subsetting code for RGList, MAList, EList, EListRaw and MArrayLM
+ objects.
+
+ 7 Dec 2013: limma 3.18.5
+
+- Bug fix to topTreat(). Rownames were incorrectly ordered when p<1.
+
+- topTable() will now work on an MArrayLM fit object that is missing
+ the lods component, for example as produced by treat().
+
+- lmFit(), eBayes(), tmixture.vector() and so on now work even when
+ there is only one gene (one row of data).
+
2 Dec 2013: limma 3.18.4
- bug fix to genas(), which was not handling vector df.prior
diff --git a/inst/doc/intro.pdf b/inst/doc/intro.pdf
index b0e97a8..8a231a7 100644
Binary files a/inst/doc/intro.pdf and b/inst/doc/intro.pdf differ
diff --git a/man/camera.Rd b/man/camera.Rd
index f481dae..89d4507 100644
--- a/man/camera.Rd
+++ b/man/camera.Rd
@@ -15,7 +15,7 @@ interGeneCorrelation(y, design)
\arguments{
\item{y}{numeric matrix giving log-expression or log-ratio values for a series of microarrays, or any data 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.}
- \item{index}{an index vector or a list of index vectors. Can be any vector such that \code{y[index,]} selects the rows corresponding to the test set.}
+ \item{index}{an index vector or a list of index vectors. Can be any vector such that \code{y[index,]} selects the rows corresponding to the test set. The list can be made using \link{symbols2indices}.}
\item{design}{design matrix.}
\item{contrast}{contrast of the linear model coefficients for which the test is required. Can be an integer specifying a column of \code{design}, or else a numeric vector of same length as the number of columns of \code{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)}.}
@@ -71,7 +71,8 @@ Analyzing gene expression data in terms of gene sets: methodological issues.
\code{\link{rankSumTestWithCorrelation}},
\code{\link{geneSetTest}},
\code{\link{roast}},
-\code{\link{romer}}.
+\code{\link{romer}},
+\code{\link{symbols2indices}}.
}
\examples{
diff --git a/man/ebayes.Rd b/man/ebayes.Rd
index 4f340d7..28130d3 100755
--- a/man/ebayes.Rd
+++ b/man/ebayes.Rd
@@ -7,7 +7,7 @@
\usage{
ebayes(fit, proportion=0.01, stdev.coef.lim=c(0.1,4), trend=FALSE, robust=FALSE, winsor.tail.p=c(0.05,0.1))
eBayes(fit, proportion=0.01, stdev.coef.lim=c(0.1,4), trend=FALSE, robust=FALSE, winsor.tail.p=c(0.05,0.1))
-treat(fit, lfc=0, trend=FALSE)
+treat(fit, lfc=0, trend=FALSE, robust=FALSE, winsor.tail.p=c(0.05,0.1))
}
\arguments{
\item{fit}{an \code{MArrayLM} fitted model object produced by \code{lmFit} or \code{contrasts.fit}.
diff --git a/man/normalizebetweenarrays.Rd b/man/normalizebetweenarrays.Rd
index 08c4557..1479cf9 100755
--- a/man/normalizebetweenarrays.Rd
+++ b/man/normalizebetweenarrays.Rd
@@ -12,7 +12,7 @@ normalizeBetweenArrays(object, method=NULL, targets=NULL, cyclic.method="fast",
\arguments{
\item{object}{a numeric \code{matrix}, \code{\link[limma:EList]{EListRaw}}, \code{\link[limma:rglist]{RGList}} or \code{\link[limma:malist]{MAList}} object containing un-normalized expression data.
- \code{matrix} or \code{EListRaw} objects will be assumed to contain single-channel data, usually log-intensities, whereas \code{RGList} or \code{MAList} objects should contain two-color data.}
+ If a matrix, then it is assumed to contain log-transformed single-channel data.}
\item{method}{character string specifying the normalization method to be used.
Choices for single-channel data are \code{"none"}, \code{"scale"}, \code{"quantile"} or \code{"cyclicloess"}.
Choices for two-color data are those previously mentioned plus \code{"Aquantile"}, \code{"Gquantile"}, \code{"Rquantile"} or \code{"Tquantile"}.
diff --git a/man/roast.Rd b/man/roast.Rd
index cae3ca0..3143845 100644
--- a/man/roast.Rd
+++ b/man/roast.Rd
@@ -28,7 +28,7 @@ Rotation gene set testing for linear models.
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.}
+ For \code{mroast}, \code{index} is a list of index vectors. The list can be made using \link{symbols2indices}.}
\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{set.statistic}{summary set statistic. Possibilities are \code{"mean"},\code{"floormean"},\code{"mean50"} or \code{"msq"}.}
@@ -111,7 +111,7 @@ The default setting for the set statistic was changed in limma 3.5.9 (3 June 201
}
\seealso{
-\code{\link{camera}}, \code{\link{romer}}, \code{\link{geneSetTest}}.
+\code{\link{camera}}, \code{\link{romer}}, \code{\link{geneSetTest}}, \code{\link{symbols2indices}}.
An overview of tests in limma is given in \link{08.Tests}.
}
diff --git a/man/subsetting.Rd b/man/subsetting.Rd
index 4e72f67..129b07e 100755
--- a/man/subsetting.Rd
+++ b/man/subsetting.Rd
@@ -5,30 +5,43 @@
\alias{[.EListRaw}
\alias{[.EList}
\alias{[.MArrayLM}
-\title{Subset RGList, MAList, EList or MArrayLM Objects}
+\alias{subsetListOfArrays}
+\title{Subset RGList, MAList, EListRaw, EList or MArrayLM Objects}
\description{
-Extract a subset of an \code{RGList}, \code{MAList}, \code{EList} or \code{MArrayLM} object.
+Return an \code{RGList}, \code{MAList}, \code{EListRaw}, \code{EList} or \code{MArrayLM} object with only selected rows and columns of the original object.
}
\usage{
-\method{[}{RGList}(object, i, j, \ldots)
+\method{[}{RGList}(object, i, j)
+subsetListOfArrays(object, i, j, IJ, IX, I, JX)
}
\arguments{
- \item{object}{object of class \code{RGList}, \code{MAList}, \code{EList} or \code{MArrayLM}}
- \item{i,j}{elements to extract. \code{i} subsets the probes or spots while \code{j} subsets the arrays}
- \item{\ldots}{not used}
+ \item{object}{object of class \code{RGList}, \code{MAList}, \code{EListRaw}, \code{EList} or \code{MArrayLM}.}
+ \item{i,j}{elements to extract. \code{i} subsets the probes or spots while \code{j} subsets the arrays.}
+ \item{IJ}{character vector giving names of components that should be subsetted by \code{i} and \code{j}.}
+ \item{IX}{character vector giving names of 2-dimensional components that should be subsetted by \code{i} only.}
+ \item{I}{character vector giving names of vector components that should be subsetted by \code{i}.}
+ \item{JX}{character vector giving names of 2-dimensional components whose row dimension corresponds to \code{j}.}
}
\details{
\code{i,j} may take any values acceptable for the matrix components of \code{object}.
+Either or both can be missing.
See the \link[base]{Extract} help entry for more details on subsetting matrices.
+
+\code{object[]} will return the whole object unchanged.
+A single index \code{object[i]} will be taken to subset rows, so \code{object[i]} and \code{object[i,]} are equivalent.
+
+\code{subsetListOfArrays} is used internally as a utility function by the subsetting operations.
+It is not intended to be called directly by users.
+Values must be supplied for all arguments other than \code{i} and \code{j}.
}
\value{
-An object of the same class as \code{object} holding data from the specified subset of genes and arrays.
+An object the same as \code{object} but containing data from the specified subset of rows and columns only.
}
\author{Gordon Smyth}
\seealso{
\code{\link[base]{Extract}} in the base package.
- \link{03.ReadingData} gives an overview of data input and manipulation functions in LIMMA.
+ \link{02.Classes} for a summary of the different data classes.
}
\examples{
M <- A <- matrix(11:14,4,2)
@@ -36,7 +49,7 @@ rownames(M) <- rownames(A) <- c("a","b","c","d")
colnames(M) <- colnames(A) <- c("A","B")
MA <- new("MAList",list(M=M,A=A))
MA[1:2,]
+MA[c("a","b"),]
MA[1:2,2]
MA[,2]
}
-\keyword{manip}
diff --git a/man/voom.Rd b/man/voom.Rd
index f0d88ac..17796f8 100644
--- a/man/voom.Rd
+++ b/man/voom.Rd
@@ -13,7 +13,7 @@ voom(counts, design = NULL, lib.size = NULL, normalize.method = "none", plot = F
\item{counts}{a numeric \code{matrix} containing raw counts, or an \code{ExpressionSet} containing raw counts, or a \code{DGEList} object.}
\item{design}{design matrix with rows corresponding to samples and columns to coefficients to be estimated. Defaults to the unit vector meaning that samples are treated as replicates.}
\item{lib.size}{numeric vector containing total library sizes for each sample.
- If \code{NULL} and \code{counts} is a \code{DGEList} then, then normalized library sizes are taken from \code{counts}.
+ If \code{NULL} and \code{counts} is a \code{DGEList} then, the normalized library sizes are taken from \code{counts}.
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.}
@@ -53,9 +53,9 @@ Law, CW (2013).
PhD Thesis. University of Melbourne, Australia.
\url{http://repository.unimelb.edu.au/10187/17598}
-Law, CW, Chen, Y, Shi, W, Smyth, GK (2013).
-Voom! Precision weights unlock linear model analysis tools for RNA-seq read counts.
-Technical Report 1 May 2013, Bioinformatics Division, Walter and Eliza Hall Institute of Medical Reseach, Melbourne, Australia.
+Law, CW, Chen, Y, Shi, W, Smyth, GK (2014).
+Voom: precision weights unlock linear model analysis tools for RNA-seq read counts.
+\emph{Genome Biology} (Accepted 9 January 2014).
\url{http://www.statsci.org/smyth/pubs/VoomPreprint.pdf}
}
diff --git a/tests/limma-Tests.R b/tests/limma-Tests.R
index b188f63..46f8b2d 100755
--- a/tests/limma-Tests.R
+++ b/tests/limma-Tests.R
@@ -218,8 +218,12 @@ y[iset1,3:4] <- y[iset1,3:4]+3
iset2 <- 6:10
roast(y=y,iset1,design,contrast=2)
roast(y=y,iset1,design,contrast=2,array.weights=c(0.5,1,0.5,1))
-mroast(y=y,list(set1=iset1,set2=iset2),design,contrast=2)
+w <- matrix(runif(100*4),100,4)
+roast(y=y,iset1,design,contrast=2,weights=w)
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))
### camera
diff --git a/tests/limma-Tests.Rout.save b/tests/limma-Tests.Rout.save
index baffee5..ef6c8f3 100755
--- a/tests/limma-Tests.Rout.save
+++ b/tests/limma-Tests.Rout.save
@@ -1152,14 +1152,28 @@ Mixed 1 0.008
Down 0 0.999
Up 1 0.002
Mixed 1 0.003
-> mroast(y=y,list(set1=iset1,set2=iset2),design,contrast=2)
- NGenes PropDown PropUp Direction PValue FDR PValue.Mixed FDR.Mixed
-set1 5 0 1 Up 0.006 0.010 0.003 0.0050
-set2 5 0 0 Up 0.566 0.565 0.547 0.5465
+> w <- matrix(runif(100*4),100,4)
+> roast(y=y,iset1,design,contrast=2,weights=w)
+ Active.Prop P.Value
+Down 0 0.999
+Up 1 0.002
+Mixed 1 0.002
> 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.004 0.006 0.009 0.0170
-set2 5 0 0 Up 0.364 0.363 0.205 0.2045
+set1 5 0 1 Up 0.006 0.010 0.008 0.0150
+set2 5 0 0 Up 0.948 0.947 0.687 0.6865
+> 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.006 0.010 0.004 0.0070
+set2 5 0 0 Up 0.716 0.715 0.658 0.6575
+> 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.004 0.006 0.003 0.0050
+set2 5 0.2 0 Down 0.984 0.983 0.250 0.2495
+> 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.002 0.002 0.001 0.0010
+set2 5 0 0 Down 0.772 0.771 0.146 0.1455
>
> ### camera
>
@@ -1176,9 +1190,9 @@ set2 5 0.1719094 Down 0.9068364381 0.90683644
> y <- new("EList",list(E=y))
> roast(y=y,iset1,design,contrast=2)
Active.Prop P.Value
-Down 0 0.999
-Up 1 0.002
-Mixed 1 0.002
+Down 0 0.998
+Up 1 0.003
+Mixed 1 0.006
> camera(y=y,iset1,design,contrast=2)
NGenes Correlation Direction PValue
set1 5 -0.2481655 Up 0.0009047749
@@ -1271,21 +1285,21 @@ Loading required package: splines
[1] "E" "weights" "design" "targets"
> summary(v$E)
V1 V2 V3 V4
- Min. :12.07 Min. :12.55 Min. :12.24 Min. :12.22
- 1st Qu.:13.11 1st Qu.:13.14 1st Qu.:13.12 1st Qu.:13.11
- Median :13.34 Median :13.32 Median :13.35 Median :13.40
- Mean :13.29 Mean :13.29 Mean :13.29 Mean :13.29
- 3rd Qu.:13.47 3rd Qu.:13.52 3rd Qu.:13.50 3rd Qu.:13.53
- Max. :14.09 Max. :13.91 Max. :13.93 Max. :13.81
+ Min. :12.25 Min. :12.58 Min. :12.19 Min. :12.24
+ 1st Qu.:13.13 1st Qu.:13.07 1st Qu.:13.15 1st Qu.:13.03
+ Median :13.29 Median :13.30 Median :13.30 Median :13.27
+ Mean :13.28 Mean :13.29 Mean :13.29 Mean :13.28
+ 3rd Qu.:13.49 3rd Qu.:13.51 3rd Qu.:13.50 3rd Qu.:13.50
+ Max. :14.23 Max. :14.28 Max. :13.97 Max. :13.96
> summary(v$weights)
V1 V2 V3 V4
- Min. : 7.242 Min. : 7.236 Min. : 7.255 Min. : 7.230
- 1st Qu.: 9.212 1st Qu.: 8.952 1st Qu.: 9.068 1st Qu.: 9.172
- Median :14.353 Median :10.616 Median :13.268 Median :15.137
- Mean :18.936 Mean :16.836 Mean :19.262 Mean :19.886
- 3rd Qu.:31.683 3rd Qu.:29.349 3rd Qu.:31.948 3rd Qu.:31.975
- Max. :32.629 Max. :32.621 Max. :32.627 Max. :32.636
+ Min. : 5.935 Min. : 5.935 Min. : 5.935 Min. : 5.935
+ 1st Qu.: 6.788 1st Qu.: 7.049 1st Qu.: 7.207 1st Qu.: 6.825
+ Median :11.066 Median :10.443 Median :10.606 Median :10.414
+ Mean :10.421 Mean :10.485 Mean :10.571 Mean :10.532
+ 3rd Qu.:13.485 3rd Qu.:14.155 3rd Qu.:13.859 3rd Qu.:14.121
+ Max. :15.083 Max. :15.101 Max. :15.095 Max. :15.063
>
> proc.time()
user system elapsed
- 1.34 0.10 1.45
+ 1.45 0.09 1.52
--
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