@@ -1,6 +1,6 @@

 Package: limma

-Version: 3.23.12

-Date: 2015/04/04

+Version: 3.23.13

+Date: 2015/04/08

 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]

@@ -12,6 +12,8 @@ S3method("[",RGList)

 S3method("[",EList)

 S3method("[",EListRaw)

 S3method(anova,MAList)

+S3method(as.data.frame,EList)

+S3method(as.data.frame,MAList)

 S3method(as.data.frame,MArrayLM)

 S3method(as.matrix,MAList)

 S3method(as.matrix,EList)

@@ -54,6 +56,7 @@ S3method("dimnames<-",RGList)

 S3method("dimnames<-",EList)

 S3method("dimnames<-",EListRaw)

 S3method(fitted,MArrayLM)

+S3method(fry,default)

 S3method(goana,default)

 S3method(goana,MArrayLM)

 S3method(length,MAList)

@@ -65,6 +68,7 @@ S3method(merge,MAList)

 S3method(merge,RGList)

 S3method(merge,EList)

 S3method(merge,EListRaw)

+S3method(mroast,default)

 S3method(roast,default)

 S3method(romer,default)

 S3method(rbind,MAList)

new file mode 100644

@@ -0,0 +1,136 @@

+fry <- function(y,...) UseMethod("fry")

+

+fry.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),weights=NULL,sort=TRUE,...)

+#	Quick version of roast gene set test assuming equal variances between genes

+#	The up and down p-values are equivalent to those from roast with nrot=Inf

+#	in the special case of prior.df=Inf.

+#	Gordon Smyth and Goknur Giner

+#	Created 30 January 2015.  Last modified 8 April 2015

+{

+#	Issue warning if extra arguments found

+	dots <- names(list(...))

+	if(length(dots)) warning("Extra arguments disregarded: ",sQuote(dots))

+

+#	Extract components from y

+	y <- getEAWP(y)

+	G <- nrow(y$exprs)

+	n <- ncol(y$exprs)

+

+#	Check index

+	if(is.null(index)) index <- list(set1=1L:G)

+	if(!is.list(index)) index <- list(set1=index)

+

+#	Check design

+	if(is.null(design)) design <- y$design

+	if(is.null(design))

+		design <- matrix(1,n,1)

+	else {

+		design <- as.matrix(design)

+		if(mode(design) != "numeric") stop("design must be a numeric matrix")

+	}

+	if(nrow(design) != n) stop("row dimension of design matrix must match column dimension of data")

+	ne <- nonEstimable(design)

+	if(!is.null(ne)) cat("Coefficients not estimable:",paste(ne,collapse=" "),"\n")

+

+	p <- ncol(design)

+	df.residual <- n-p

+	df.camera <- min(df.residual,G-2)

+

+#	Check weights

+	if(is.null(weights)) weights <- y$weights

+

+#	Reduce to numeric expression matrix

+	y <- y$exprs

+

+#	Check weights

+	if(!is.null(weights)) {

+		if(any(weights<=0)) stop("weights must be positive")

+		if(length(weights)==n) {

+			sw <- sqrt(weights)

+			y <- t(t(y)*sw)

+			design <- design*sw

+			weights <- NULL

+		}

+	}

+	if(!is.null(weights)) {

+		if(length(weights)==G) weights <- matrix(weights,G,n)

+		weights <- as.matrix(weights)

+		if(any( dim(weights) != dim(y) )) stop("weights not conformal with y")

+	}

+

+#	Reform design matrix so that contrast of interest is last column

+	if(is.character(contrast)) {

+		contrast <- which(contrast==colnames(design))

+		if(length(contrast)==0) stop("coef ",contrast," not found")

+	}

+	if(length(contrast)==1) {

+		j <- c((1:p)[-contrast], contrast)

+		if(contrast<p) design <- design[,j]

+	} else {

+		QR <- qr(contrast)

+		design <- t(qr.qty(QR,t(design)))

+		if(sign(QR$qr[1,1]<0)) design[,1] <- -design[,1]

+		design <- design[,c(2:p,1)]

+	}

+

+#	Compute effects matrix

+	if(is.null(weights)) {

+		QR <- qr(design)

+		if(QR$rank<p) stop("design matrix is not of full rank")

+		effects <- qr.qty(QR,t(y))

+		unscaledt <- effects[p,]

+		if(QR$qr[p,p]<0) unscaledt <- -unscaledt

+	} else {

+		effects <- matrix(0,n,G)

+		unscaledt <- rep(0,n)

+		sw <- sqrt(weights)

+		yw <- y*sw

+		for (g in 1:G) {

+			xw <- design*sw[g,]

+			QR <- qr(xw)

+			if(QR$rank<p) stop("weighted design matrix not of full rank for gene ",g)

+			effects[,g] <- qr.qty(QR,yw[g,])

+			unscaledt[g] <- effects[p,g]

+			if(QR$qr[p,p]<0) unscaledt[g] <- -unscaledt[g]

+		}

+	}

+

+#	Standardized residuals

+	U <- effects[-(1:(p-1)),,drop=FALSE]

+

+#	Global statistics

+	nsets <- length(index)

+	NGenes <- rep.int(0L,nsets)

+	Direction <- rep.int("",nsets)

+	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)

+		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)

+	}

+

+#	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)

+	} else {

+		tab <- data.frame(NGenes=NGenes,Direction=Direction,PValue=PValue)

+	}

+	rownames(tab) <- names(index)

+

+#	Sort by p-value

+	if(sort && nsets>1) {

+		o <- order(tab$PValue,-tab$NGenes)

+		tab <- tab[o,]

+	}

+

+	tab

+}

@@ -16,7 +16,7 @@ roast <- function(y,...) UseMethod("roast")

 roast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.statistic="mean",gene.weights=NULL,array.weights=NULL,weights=NULL,block=NULL,correlation,var.prior=NULL,df.prior=NULL,trend.var=FALSE,nrot=999,approx.zscore=TRUE,...)

 # Rotation gene set testing for linear models

 # Gordon Smyth and Di Wu

-# Created 24 Apr 2008.  Last modified 31 March 2015.

+# Created 24 Apr 2008.  Last modified 8 April 2015.

 {

 #	Issue warning if extra arguments found

 	dots <- names(list(...))

@@ -483,6 +483,11 @@ mroast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.st

 		stringsAsFactors=FALSE

 	)

+	if(midp) {

+		tab$FDR <- pmax(tab$FDR, pv[,"UpOrDown"])

+		tab$FDR.Mixed <- pmax(tab$FDR.Mixed, pv[,"Mixed"])

+	}

+

 #	Sort by p-value

 	sort <- match.arg(sort,c("directional","mixed","none"))

 	if(sort=="none") return(tab)

@@ -1,3 +1,11 @@

+ 8 April 2015: limma 3.23.13

+

+- new function fry(), which provides a fast version of mroast()

+  when there is little heteroscedasticity between genes.

+ 

+- minor change to mroast() output so that FDR values are never

+  smaller than p-values when midp=TRUE.

+ 

  4 April 2015: limma 3.23.12

 - camera(), roast(), mroast() and romer() now give a warning if

@@ -5,19 +5,25 @@

 \alias{mroast.default}

 \alias{Roast-class}

 \alias{show,Roast-method}

+\alias{fry}

+\alias{fry.default}

 \title{Rotation Gene Set Tests}

 \description{

 Rotation gene set testing for linear models.

 }

 \usage{

-\S3method{roast}{default}(y, index=NULL, design=NULL, contrast=ncol(design), set.statistic="mean",

-     gene.weights=NULL, array.weights=NULL, weights=NULL, block=NULL, correlation,

-     var.prior=NULL, df.prior=NULL, trend.var=FALSE, nrot=999, approx.zscore=TRUE, \dots)

-\S3method{mroast}{default}(y, index=NULL, design=NULL, contrast=ncol(design), set.statistic="mean",

-     gene.weights=NULL, array.weights=NULL, weights=NULL, block=NULL, correlation,

-     var.prior=NULL, df.prior=NULL, trend.var=FALSE, nrot=999, approx.zscore=TRUE,

-     adjust.method="BH", midp=TRUE, sort="directional", \dots)

+\S3method{fry}{default}(y, index = NULL, design = NULL, contrast = ncol(design),

+      weights = NULL, sort = TRUE, \dots)

+\S3method{roast}{default}(y, index = NULL, design = NULL, contrast = ncol(design),

+      set.statistic = "mean", gene.weights = NULL, array.weights = NULL, weights = NULL,

+      block = NULL, correlation, var.prior = NULL, df.prior = NULL, trend.var = FALSE,

+      nrot = 999, approx.zscore = TRUE, \dots)

+\S3method{mroast}{default}(y, index = NULL, design = NULL, contrast = ncol(design),

+       set.statistic = "mean", gene.weights = NULL, array.weights = NULL, weights = NULL,

+       block = NULL, correlation, var.prior = NULL, df.prior = NULL, trend.var = FALSE,

+       nrot = 999, approx.zscore = TRUE, adjust.method = "BH", midp = TRUE,

+       sort = "directional", \dots)

 }

 \arguments{

@@ -40,7 +46,7 @@ Rotation gene set testing for linear models.

   \item{nrot}{number of rotations used to estimate the p-values.}

   \item{adjust.method}{method used to adjust the p-values for multiple testing. See \code{\link{p.adjust}} for possible values.}

   \item{midp}{logical, should mid-p-values be used in instead of ordinary p-values when adjusting for multiple testing?}

-  \item{sort}{character, whether to sort output table by directional p-value (\code{"directional"}), non-directional p-value (\code{"mixed"}), or not at all (\code{"none"}).}

+  \item{sort}{character, whether to sort output table by directional p-value (\code{"directional"}), non-directional p-value (\code{"mixed"}), or not at all (\code{"none"}).  For \code{fry} it is logical, whether to sort or not.}

   \item{approx.zscore}{logical, if \code{TRUE} then a fast approximation is used to convert t-statistics into z-scores prior to computing set statistics. If \code{FALSE}, z-scores will be exact.}

   \item{\dots}{other arguments not currently used.}

 }

@@ -103,6 +109,9 @@ This means that the smallest possible p-value is \code{1/(nrot+1)}.

 \code{mroast} does roast tests for multiple sets, including adjustment for multiple testing.

 By default, \code{mroast} reports ordinary p-values but uses mid-p-values (Routledge, 1994) at the multiple testing stage.

 Mid-p-values are probably a good choice when using false discovery rates (\code{adjust.method="BH"}) but not when controlling the family-wise type I error rate (\code{adjust.method="holm"}).

+

+\code{fry} is a fast version of \code{mroast} in the special case that \code{df.prior} is large and \code{set.statistic="mean"}.

+In this situation, it gives the same result as \code{mroast} with an infinite number of rotations.

 }

 \note{

@@ -1,5 +1,5 @@

-R version 3.2.0 alpha (2015-03-28 r68110)

+R version 3.2.0 beta (2015-04-02 r68137) -- "Full of Ingredients"

 Copyright (C) 2015 The R Foundation for Statistical Computing

 Platform: x86_64-w64-mingw32/x64 (64-bit)

@@ -1173,21 +1173,21 @@ 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.0150        0.008    0.0150

-set2      5        0      0        Up  0.959 0.9585        0.687    0.6865

+     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.0070        0.004    0.0070

-set2      5        0      0        Up  0.679 0.6785        0.658    0.6575

+     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.0050        0.003    0.0050

-set2      5      0.2      0      Down  0.950 0.9495        0.250    0.2495

+     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.0010        0.001    0.0010

-set2      5        0      0      Down  0.791 0.7905        0.146    0.1455

+     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

 > 

@@ -1376,4 +1376,4 @@ GO:0070062 0.055161144

 > 

 > proc.time()

    user  system elapsed 

-   2.94    0.23    3.16 

+   3.13    0.24    4.36