[med-svn] [r-bioc-limma] 01/03: Imported Upstream version 3.18.0~dfsg
Charles Plessy
plessy at alioth.debian.org
Wed Oct 16 06:44:44 UTC 2013
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 5c7fc4d452ad9e4fa48e2fc5503808976838aa5c
Author: Charles Plessy <plessy at debian.org>
Date: Wed Oct 16 14:29:25 2013 +0900
Imported Upstream version 3.18.0~dfsg
---
DESCRIPTION | 12 +--
NAMESPACE | 4 +
R/alias2Symbol.R | 6 +-
R/beadCountWeights.R | 99 +++++++++++++++++
R/contrastAsCoef.R | 26 +++++
R/genas.R | 228 +++++++++++++++++++++------------------
R/geneset-roast.R | 27 ++---
R/lmfit.R | 98 ++++++++++++-----
R/loessFit.R | 50 ++++-----
R/norm.R | 4 +-
R/plots-ma.R | 251 ++++++++++++++++++++++++++-----------------
R/predFCm.R | 2 +-
R/read-ilmn.R | 36 +++++--
R/subsetting.R | 3 +-
R/toptable.R | 184 +++++++++++++++++++++----------
R/treat.R | 63 ++++++++---
build/vignette.rds | Bin 0 -> 219 bytes
inst/doc/changelog.txt | 180 +++++++++++++++++++++++++++----
inst/doc/limma.R | 2 -
inst/doc/limma.pdf | Bin 45597 -> 45597 bytes
man/01Introduction.Rd | 3 +-
man/06linearmodels.Rd | 3 +-
man/asmalist.Rd | 7 ++
man/beadCountWeights.Rd | 61 +++++++++++
man/contrastAsCoef.Rd | 48 +++++++++
man/contrasts.fit.Rd | 2 +-
man/ebayes.Rd | 15 ++-
man/fitfdist.Rd | 4 +-
man/genas.Rd | 88 ++++++++++-----
man/getEAWP.Rd | 8 +-
man/lmFit.Rd | 3 +
man/loessfit.Rd | 86 +++++++++------
man/plotma.Rd | 9 +-
man/predFCm.Rd | 1 +
man/propTrueNull.Rd | 1 +
man/removeBatchEffect.Rd | 4 +-
man/roast.Rd | 1 +
man/squeezeVar.Rd | 23 +++-
man/targetsA2C.Rd | 2 +-
man/toptable.Rd | 15 ++-
src/normexp.c | 4 +-
tests/limma-Tests.R | 13 ++-
tests/limma-Tests.Rout.save | 229 +++++++++++++++++++++++----------------
43 files changed, 1347 insertions(+), 558 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 9dd2969..6bfd4f1 100755
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,12 +1,12 @@
Package: limma
-Version: 3.16.8
-Date: 2013/09/19
+Version: 3.18.0
+Date: 2013/10/11
Title: Linear Models for Microarray Data
-Author: Gordon Smyth with contributions from Matthew Ritchie, Jeremy Silver, James Wettenhall, Natalie Thorne, Mette Langaas, Egil Ferkingstad, Marcus Davy, Francois Pepin, Dongseok Choi, Davis McCarthy, Di Wu, Alicia Oshlack, Carolyn de Graaf, Yifang Hu, Wei Shi and Belinda Phipson.
+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>
Depends: R (>= 2.3.0), methods
-Suggests: affy, MASS, org.Hs.eg.db, splines, statmod (>= 1.2.2), vsn,
- ellipse
+Suggests: statmod (>= 1.2.2), splines, locfit, MASS, ellipse, affy,
+ vsn, AnnotationDbi, org.Hs.eg.db
LazyLoad: yes
Description: Data analysis, linear models and differential expression for microarray data.
License: GPL (>=2)
@@ -14,4 +14,4 @@ URL: http://bioinf.wehi.edu.au/limma
biocViews: Microarray, OneChannel, TwoChannel, DataImport,
QualityControl, Preprocessing, Bioinformatics,
DifferentialExpression, MultipleComparisons, TimeCourse
-Packaged: 2013-09-19 04:07:58 UTC; biocbuild
+Packaged: 2013-10-15 03:16:09 UTC; biocbuild
diff --git a/NAMESPACE b/NAMESPACE
index dd7d911..ea6a1d5 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -68,6 +68,10 @@ S3method(summary,EListRaw)
S3method(summary,TestResults)
S3method(normalizeVSN,RGList)
S3method(normalizeVSN,EListRaw)
+S3method(plotMA,RGList)
+S3method(plotMA,MAList)
+S3method(plotMA,EList)
+S3method(plotMA,MArrayLM)
S3method(plotMDS,MDS)
S3method(roast,MAList)
S3method(roast,EList)
diff --git a/R/alias2Symbol.R b/R/alias2Symbol.R
index 3ad50ed..ecac981 100644
--- a/R/alias2Symbol.R
+++ b/R/alias2Symbol.R
@@ -4,13 +4,14 @@ alias2Symbol <- function(alias,species="Hs",expand.symbols=FALSE)
# Convert a set of alias names to official gene symbols
# via Entrez Gene identifiers
# Di Wu, Gordon Smyth and Yifang Hu
-# 4 Sep 2008. Last revised 15 Jan 2009
+# 4 Sep 2008. Last revised 29 September 2013
{
alias <- as.character(alias)
species <- match.arg(species,c("Dm","Hs","Mm","Rn"))
DB <- paste("org",species,"eg","db",sep=".")
ALIAS2EG <- paste("org",species,"egALIAS2EG",sep=".")
SYMBOL <- paste("org",species,"egSYMBOL",sep=".")
+ suppressPackageStartupMessages(require("AnnotationDbi",character.only=TRUE))
suppressPackageStartupMessages(require(DB,character.only=TRUE))
if(expand.symbols)
{
@@ -32,13 +33,14 @@ alias2SymbolTable <- function(alias,species="Hs")
# Convert a vector of alias names to the vector of corresponding official gene symbols
# via Entrez Gene identifiers
# Di Wu, Gordon Smyth and Yifang Hu
-# Created 3 Sep 2009. Last modified 23 Dec 2012.
+# Created 3 Sep 2009. Last modified 29 September Dec 2013.
{
alias <- as.character(alias)
species <- match.arg(species,c("Dm","Hs","Mm","Rn"))
DB <- paste("org",species,"eg","db",sep=".")
ALIAS2EG <- paste("org",species,"egALIAS2EG",sep=".")
SYMBOL <- paste("org",species,"egSYMBOL",sep=".")
+ suppressPackageStartupMessages(require("AnnotationDbi",character.only=TRUE))
suppressPackageStartupMessages(require(DB,character.only=TRUE))
isSymbol <- alias %in% Rkeys(get(SYMBOL))
diff --git a/R/beadCountWeights.R b/R/beadCountWeights.R
new file mode 100644
index 0000000..1a8a6b8
--- /dev/null
+++ b/R/beadCountWeights.R
@@ -0,0 +1,99 @@
+beadCountWeights <- function(y, x, design=NULL, bead.stdev=NULL, nbeads=NULL, array.cv=TRUE, scale=FALSE)
+# Compute weights for BeadChips based on bead-level counts and standard deviations
+# Charity Law and Gordon Smyth
+# 4 August 2010. Last modified 17 July 2012.
+{
+ E <- as.matrix(y)
+ E.raw <- as.matrix(x)
+ if(is.null(bead.stdev)) {
+ bead.stdev <- y$other$BEAD_STDEV
+ if(is.null(bead.stdev)) stop("BEAD_STDEV not found in data object")
+ }
+ if(is.null(nbeads)) {
+ nbeads <- y$other$Avg_NBEADS
+ if(is.null(nbeads)) stop("NBEADS not found in data object")
+ }
+ P <- nrow(E)
+ A <- ncol(E)
+ if(nrow(E.raw) != P) stop("dimensions don't match")
+ if(ncol(E.raw) != A) stop("dimensions don't match")
+ if(nrow(bead.stdev) != P) stop("dimensions don't match")
+ if(ncol(bead.stdev) != A) stop("dimensions don't match")
+ if(nrow(nbeads) != P) stop("dimensions don't match")
+ if(ncol(nbeads) != A) stop("dimensions don't match")
+
+# Coefficient of variation of bead-level observations
+# Array-specific or constant.
+ cv <- bead.stdev/E.raw
+ cv.array <- apply(sqrt(cv), 2, function(k) mean(k, trim=0.125, na.rm=TRUE))^2
+ cv.constant <- mean(sqrt(cv), trim=0.125, na.rm=TRUE)^2
+ if(array.cv) {cv <- rep(cv.array, each=nrow(y))}
+ else cv <- cv.constant
+
+# Predicted variance of normalized probe-level values
+ tv <- log(cv^2/nbeads+1)/log(2)^2
+
+# Squared-residuals to calculate biological variance
+ qrX <- qr(design)
+ res <- qr.resid(qrX, t(E))
+ r2 <- (t(res))^2
+
+# Leverages
+ Q <- qr.Q(qrX)
+ h <- rowSums(Q^2)
+ h <- matrix(rep(h, each=P), nrow=P, ncol=A)
+
+ bv <- .ilmn.biological.variance(var.from.beads=tv, sq.residuals=r2, leverage=h)
+
+ out <- list()
+ out$var.biological <- bv
+ out$var.technical <- tv
+ out$cv.array <- cv.array
+ out$cv.constant <- cv.constant
+ out$weights <- 1/(tv + bv)
+ out$weights <- out$weights/rowMeans(out$weights)
+ if(scale) out$weights <- out$weights*rowMeans(1/tv)
+ out
+}
+
+
+.ilmn.biological.variance <- function(var.from.beads, sq.residuals, leverage)
+# Find the component of the between-array variance
+# not explained by bead-level variability
+# Gordon Smyth and Charity Law
+# 28 July 2010. Modified 16 July 2012.
+{
+ if(any(var.from.beads < 0)) stop("negative variances not allowed")
+ if(any(sq.residuals < 0)) stop("negative variances not allowed")
+ tv <- as.matrix(var.from.beads)
+ r2 <- as.matrix(sq.residuals)
+ h <- as.matrix(leverage)
+ P <- nrow(tv)
+ A <- ncol(tv)
+ if(nrow(r2) != P) stop("dimensions don't match")
+ if(ncol(r2) != A) stop("dimensions don't match")
+ if(nrow(h) != P) stop("dimensions don't match")
+ if(ncol(h) != A) stop("dimensions don't match")
+
+ # Newton's iteration
+ F <- rowMeans(r2/(2*tv^2)-(1-h)/(2*tv))
+ bv <- rep(0,length=P)
+ i <- (F > 0)
+ iter <- 0
+ while(any(i)) {
+ iter <- iter+1
+ if(iter > 200) {
+ warning("More than 200 iterations required")
+ return(bv)
+ }
+ Fdash <- - rowMeans(r2[i,,drop=FALSE]/(bv[i]+tv[i,,drop=FALSE])^3 - (1-h[i,,drop=FALSE])/(2*(bv[i]+tv[i,,drop=FALSE])^2))
+ step <- - F[i]/Fdash
+ bv[i] <- bv[i] + step
+ F[i] <- rowMeans(r2[i,,drop=FALSE]/(2*(bv[i]+tv[i,,drop=FALSE])^2) - (1-h[i,,drop=FALSE])/(2*(bv[i]+tv[i,,drop=FALSE])))
+ i[i] <- (step > 1e-5)
+ print(sum(i))
+ }
+ bv
+}
+
+
diff --git a/R/contrastAsCoef.R b/R/contrastAsCoef.R
new file mode 100644
index 0000000..efa3dbc
--- /dev/null
+++ b/R/contrastAsCoef.R
@@ -0,0 +1,26 @@
+contrastAsCoef <- function(design, contrast=NULL, first=TRUE)
+# Reform a design matrix so that one or more contrasts become simple coefficients
+# Gordon Smyth
+# 31 August 2013
+{
+ design <- as.matrix(design)
+ if(is.null(contrast)) return(design)
+ contrast <- as.matrix(contrast)
+ if(ncol(design) != nrow(contrast)) stop("Length of contrast doesn't match ncol(design)")
+ qrc <- qr(contrast)
+ ncontrasts <- qrc$rank
+ if(ncontrasts==0) stop("contrast is all zero")
+ coef <- 1:ncontrasts
+ Dvec <- diag(qrc$qr)[coef]
+ design <- t(qr.qty(qrc,t(design)))
+ colnames(design) <- paste("Q",1:ncol(design),sep="")
+ cn <- colnames(contrast)
+ if(is.null(cn)) cn <- paste("C",qrc$pivot[coef],sep="")
+ colnames(design)[coef] <- cn
+ if(!first) {
+ design <- cbind(design[,-coef,drop=FALSE],design[,coef,drop=FALSE])
+ coef <- rev( ncol(design)-coef+1 )
+ }
+ design[,coef] <- t(t(design[,coef])/Dvec)
+ list(design=design,coef=coef,qr=qrc)
+}
diff --git a/R/genas.R b/R/genas.R
index 4829f53..787baab 100644
--- a/R/genas.R
+++ b/R/genas.R
@@ -1,28 +1,58 @@
## GENAS.R
-genas <- function(fit,coef=c(1,2),chooseMethod=NULL,plot=FALSE,alpha=0.4)
+genas <- function(fit,coef=c(1,2),subset="all",plot=FALSE,alpha=0.4)
# Genuine association of gene expression profiles
# Belinda Phipson and Gordon Smyth
-# 21 September 2009. Last modified 2 December 2012.
+# 21 September 2009. Last modified 26 July 2013.
{
- if(ncol(fit)>2) fit<-fit[,coef]
- if(length(fit$s2.prior)==1) trend<-FALSE else trend<-TRUE
- fit <- eBayes(fit,trend=trend)
- if(is.null(chooseMethod)) chooseMethod <- "n"
-
- x1<-fitGammaIntercept(fit$coeff[,1]^2/fit$s2.post,offset=fit$cov.coeff[1,1])
- x2<-fitGammaIntercept(fit$coeff[,2]^2/fit$s2.post,offset=fit$cov.coeff[2,2])
- if(x1 > 0 & x2 > 0){
- v0null<-matrix(c(x1,0,0,x2),2,2)
- C<-chol(v0null)
- x<-log(diag(C))
+ out <- list(
+ technical.correlation=NA,
+ covariance.matrix=matrix(NA,2,2),
+ biological.correlation=NA,
+ deviance=0,
+ p.value=1,
+ n=0
+ )
+
+# Check fit
+ if(nrow(fit)<1) {
+ message("fit object has zero rows")
+ return(out)
}
- else x<-c(0,0)
+ if(is.null(fit$s2.post)) fit <- eBayes(fit)
+ trend <- (length(fit$s2.prior) > 1)
+ robust <- (length(fit$df.prior) > 1)
+
+# Check coef
+ if(length(coef) != 2) stop("coef should contain two column numbers")
+
+# Check subset
+ if(subset=="n") subset <- "all" # for backward compatibility
+ subset <- match.arg(subset, c("all","Fpval","p.union","p.int","logFC","predFC"))
- m<-2
+# Keep only the two fit contrasts to be correlated
+ if(ncol(fit)>2) fit <- fit[,coef]
fit.plot <- fit
- fit <- .whichGenes(fit,chooseMethod)
- fit <- eBayes(fit,trend=trend)
+
+ x1 <- fitGammaIntercept(fit$coeff[,1]^2/fit$s2.post,offset=fit$cov.coeff[1,1])
+ x2 <- fitGammaIntercept(fit$coeff[,2]^2/fit$s2.post,offset=fit$cov.coeff[2,2])
+ if(x1 > 0 && x2 > 0) {
+ v0null <- matrix(c(x1,0,0,x2),2,2)
+ C <- chol(v0null)
+ x <- log(diag(C))
+ } else
+ x <- c(0,0)
+ m <- 2
+
+# Subset to genes that show some differential expression
+ if(subset != "all") {
+ fit <- .whichGenes(fit,subset)
+ if(nrow(fit)<1) {
+ message("No genes met criterion for inclusion in analysis")
+ return(out)
+ }
+ fit <- eBayes(fit,trend=trend,robust=robust)
+ }
Q2 <- optim(x, .multTLogLikNull, fit = fit, m = m)
Q1 <- optim(c(Q2$par[1], Q2$par[2], 0),.multTLogLik,fit=fit,m=m)
@@ -35,33 +65,32 @@ genas <- function(fit,coef=c(1,2),chooseMethod=NULL,plot=FALSE,alpha=0.4)
V<-fit$cov.coefficients
rhotech<-V[2,1]/sqrt(V[1,1]*V[2,2])
- if(plot){
- require(ellipse)
- lim<-mean(c(sd(fit.plot$coeff[,1]),sd(fit.plot$coeff[,2])))
- if(nrow(fit)<500) lim <- 1.5*lim else lim <- 2*lim
- max<-max(c(fit.plot$coeff[,1],fit.plot$coeff[,2]))
- min<-min(c(fit.plot$coeff[,1],fit.plot$coeff[,2]))
- max<-sign(max)*max(abs(min),abs(max))
- min<-sign(min)*max(abs(min),abs(max))
- if(abs(rhobiol)>abs(rhotech)){
- plot(fit.plot$coeff[,1],fit.plot$coeff[,2],pch=16,cex=0.4,ylim=c(min,max),xlim=c(min,max),xlab=colnames(fit.plot$coeff)[1],ylab=colnames(fit.plot$coeff)[2])
- polygon(ellipse(rhotech,scale=c(lim,lim)),col=rgb(0,0,1,alpha=alpha),border=rgb(0,0,1,alpha=alpha))
- polygon(ellipse(rhobiol,scale=c(lim,lim)),col=rgb(0,1,0,alpha=alpha),border=rgb(0,1,0,alpha=alpha))
- points(fit.plot$coeff[,1],fit.plot$coeff[,2],pch=16,cex=0.4)
- abline(h=0,v=0)
- legend(min,max,legend=bquote(rho[biol]==.(round(rhobiol,3))),col=rgb(0,1,0,alpha=alpha),pch=16,bty="n",cex=0.8)
- legend(min,max-lim/2,legend=bquote(rho[tech]==.(round(rhotech,3))),col=rgb(0,0,1,alpha=alpha),pch=16,bty="n",cex=0.8)
- }
- else{
- plot(fit.plot$coeff[,1],fit.plot$coeff[,2],pch=16,cex=0.4,ylim=c(min,max),xlim=c(min,max),xlab=colnames(fit.plot$coeff)[1],ylab=colnames(fit.plot$coeff)[2])
- polygon(ellipse(rhobiol,scale=c(lim,lim)),col=rgb(0,1,0,alpha=alpha),border=rgb(0,1,0,alpha=alpha))
- polygon(ellipse(rhotech,scale=c(lim,lim)),col=rgb(0,0,1,alpha=alpha),border=rgb(0,0,1,alpha=alpha))
- points(fit.plot$coeff[,1],fit.plot$coeff[,2],pch=16,cex=0.4)
- abline(h=0,v=0)
- legend(min,max,legend=bquote(rho[biol]==.(round(rhobiol,3))),col=rgb(0,1,0,alpha=alpha),pch=16,bty="n",cex=0.8)
- legend(min,max-lim/2,legend=bquote(rho[tech]==.(round(rhotech,3))),col=rgb(0,0,1,alpha=alpha),pch=16,bty="n",cex=0.8)
- }
+ if(plot) {
+ require(ellipse)
+ lim <- mean(c(sd(fit.plot$coeff[,1]),sd(fit.plot$coeff[,2])))
+ if(nrow(fit)<500) lim <- 1.5*lim else lim <- 2*lim
+ max <- max(c(fit.plot$coeff[,1],fit.plot$coeff[,2]))
+ min <- min(c(fit.plot$coeff[,1],fit.plot$coeff[,2]))
+ max <- sign(max)*max(abs(min),abs(max))
+ min <- sign(min)*max(abs(min),abs(max))
+ if(abs(rhobiol)>abs(rhotech)){
+ plot(fit.plot$coeff[,1],fit.plot$coeff[,2],pch=16,cex=0.4,ylim=c(min,max),xlim=c(min,max),xlab=colnames(fit.plot$coeff)[1],ylab=colnames(fit.plot$coeff)[2])
+ polygon(ellipse(rhotech,scale=c(lim,lim)),col=rgb(0,0,1,alpha=alpha),border=rgb(0,0,1,alpha=alpha))
+ polygon(ellipse(rhobiol,scale=c(lim,lim)),col=rgb(0,1,0,alpha=alpha),border=rgb(0,1,0,alpha=alpha))
+ points(fit.plot$coeff[,1],fit.plot$coeff[,2],pch=16,cex=0.4)
+ abline(h=0,v=0)
+ legend(min,max,legend=bquote(rho[biol]==.(round(rhobiol,3))),col=rgb(0,1,0,alpha=alpha),pch=16,bty="n",cex=0.8)
+ legend(min,max-lim/2,legend=bquote(rho[tech]==.(round(rhotech,3))),col=rgb(0,0,1,alpha=alpha),pch=16,bty="n",cex=0.8)
+ } else {
+ plot(fit.plot$coeff[,1],fit.plot$coeff[,2],pch=16,cex=0.4,ylim=c(min,max),xlim=c(min,max),xlab=colnames(fit.plot$coeff)[1],ylab=colnames(fit.plot$coeff)[2])
+ polygon(ellipse(rhobiol,scale=c(lim,lim)),col=rgb(0,1,0,alpha=alpha),border=rgb(0,1,0,alpha=alpha))
+ polygon(ellipse(rhotech,scale=c(lim,lim)),col=rgb(0,0,1,alpha=alpha),border=rgb(0,0,1,alpha=alpha))
+ points(fit.plot$coeff[,1],fit.plot$coeff[,2],pch=16,cex=0.4)
+ abline(h=0,v=0)
+ legend(min,max,legend=bquote(rho[biol]==.(round(rhobiol,3))),col=rgb(0,1,0,alpha=alpha),pch=16,bty="n",cex=0.8)
+ legend(min,max-lim/2,legend=bquote(rho[tech]==.(round(rhotech,3))),col=rgb(0,0,1,alpha=alpha),pch=16,bty="n",cex=0.8)
}
+ }
D<-abs(2*(Q2$value-Q1$value))
p.val<-pchisq(D,df=1,lower.tail=FALSE)
@@ -99,9 +128,9 @@ genas <- function(fit,coef=c(1,2),chooseMethod=NULL,plot=FALSE,alpha=0.4)
.multTLogLik <- function(x,fit,m)
-# Calculate the log-likelihood with biological correlation
-# Belinda Phipson and Gordon Smyth
-# 21 September 2009. Last modified 21 September 2009.
+# Calculate the log-likelihood with biological correlation
+# Belinda Phipson and Gordon Smyth
+# 21 September 2009. Last modified 21 September 2009.
{
d0<-fit$df.prior
d<-fit$df.residual
@@ -130,72 +159,61 @@ genas <- function(fit,coef=c(1,2),chooseMethod=NULL,plot=FALSE,alpha=0.4)
}
-.whichGenes <- function(fit,chooseMethod){
- if(chooseMethod=="Fpval"){
- p <- 1-propTrueNull(fit$F.p.value)
- R <- rank(fit$F.p.value)
- cut <- p*nrow(fit)
- genes <- R <= cut
- if(sum(genes)==0){
- # message("There is no evidence of differential expression. LogFC cut-off used instead.")
- chooseMethod<-"logFC"
- }
- else {
- if(sum(genes)<500) message("Less than 500 DE genes. Correlation estimate may be inaccurate.")
- }
+.whichGenes <- function(fit,subset)
+{
+ if(subset=="all") return(genes)
+
+ if(subset=="Fpval") {
+ p <- 1-propTrueNull(fit$F.p.value)
+ R <- rank(fit$F.p.value)
+ cut <- p*nrow(fit)
+ genes <- (R <= cut)
}
- if(chooseMethod=="p.union"){
- p1 <- 1-propTrueNull(fit$p.value[,1])
- p2 <- 1-propTrueNull(fit$p.value[,2])
- cut1 <- p1*nrow(fit)
- cut2 <- p2*nrow(fit)
- if(p1==0 & p2==0){
- # message("There is no evidence of differential expression. LogFC cut-off used instead.")
- chooseMethod<-"logFC"
+
+ if(subset=="p.union"){
+ p1 <- 1-propTrueNull(fit$p.value[,1])
+ p2 <- 1-propTrueNull(fit$p.value[,2])
+ cut1 <- p1*nrow(fit)
+ cut2 <- p2*nrow(fit)
+ if(p1==0 & p2==0){
+ genes <- FALSE
+ } else {
+ R1 <- rank(fit$p.value[,1])
+ R2 <- rank(fit$p.value[,2])
+ genes <- R1 <= cut1 | R2 <= cut2
}
- else {
- if(cut1+cut2 < 500) message("Less than 500 DE genes. Correlation estimate may be inaccurate.")
+ }
+
+ if(subset=="p.int"){
+ p1 <- 1-propTrueNull(fit$p.value[,1])
+ p2 <- 1-propTrueNull(fit$p.value[,2])
R1 <- rank(fit$p.value[,1])
R2 <- rank(fit$p.value[,2])
- genes <- R1 <= cut1 | R2 <= cut2
- }
+ cut1 <- p1*nrow(fit)
+ cut2 <- p2*nrow(fit)
+ genes <- R1 <= cut1 & R2 <= cut2
}
- if(chooseMethod=="p.int"){
- p1 <- 1-propTrueNull(fit$p.value[,1])
- p2 <- 1-propTrueNull(fit$p.value[,2])
- R1 <- rank(fit$p.value[,1])
- R2 <- rank(fit$p.value[,2])
- cut1 <- p1*nrow(fit)
- cut2 <- p2*nrow(fit)
- genes <- R1 <= cut1 & R2 <= cut2
- if(sum(genes)==0){
- # message("There is no evidence of differential expression. LogFC cut-off used instead.")
- chooseMethod<-"logFC"
- }
- else {
- if(sum(genes)<500) message("Less than 500 DE genes. Correlation estimate may be inaccurate.")
- }
+
+ if(subset=="logFC") {
+ q1 <- quantile(abs(fit$coeff[,1]),probs=0.9)
+ q2 <- quantile(abs(fit$coeff[,2]),probs=0.9)
+ genes <- abs(fit$coeff[,1]) > q1 | abs(fit$coeff[,2]) > q2
+ fit$coeff[,1] <- sign(fit$coeff[,1]) * (abs(fit$coeff[,1])-q1)
+ fit$coeff[,2] <- sign(fit$coeff[,2]) * (abs(fit$coeff[,2])-q2)
}
- if(chooseMethod=="logFC"){
- q1 <- quantile(abs(fit$coeff[,1]),probs=0.9)
- q2 <- quantile(abs(fit$coeff[,2]),probs=0.9)
- genes <- abs(fit$coeff[,1]) > q1 | abs(fit$coeff[,2]) > q2
- fit$coeff[,1] <- sign(fit$coeff[,1]) * (abs(fit$coeff[,1])-q1)
- fit$coeff[,2] <- sign(fit$coeff[,2]) * (abs(fit$coeff[,2])-q2)
+
+ if(subset=="predFC") {
+ pfc1 <- predFCm(fit, coef=1)
+ pfc2 <- predFCm(fit, coef=2)
+ q1 <- quantile(abs(pfc1),probs=0.9)
+ q2 <- quantile(abs(pfc2),probs=0.9)
+ genes <- abs(pfc1) > q1 | abs(pfc2) > q2
+ fit$coeff[,1] <- pfc1
+ fit$coeff[,2] <- pfc2
+ fit$coeff[,1] <- sign(fit$coeff[,1]) * (abs(fit$coeff[,1])-q1)
+ fit$coeff[,2] <- sign(fit$coeff[,2]) * (abs(fit$coeff[,2])-q2)
}
- if(chooseMethod=="predFC"){
- pfc1 <- predFCm(fit, coef=1)
- pfc2 <- predFCm(fit, coef=2)
- q1 <- quantile(abs(pfc1),probs=0.9)
- q2 <- quantile(abs(pfc2),probs=0.9)
- genes <- abs(pfc1) > q1 | abs(pfc2) > q2
- fit$coeff[,1] <- pfc1
- fit$coeff[,2] <- pfc2
- fit$coeff[,1] <- sign(fit$coeff[,1]) * (abs(fit$coeff[,1])-q1)
- fit$coeff[,2] <- sign(fit$coeff[,2]) * (abs(fit$coeff[,2])-q2)
- }
- if(chooseMethod=="n") genes <- c(rep(TRUE,nrow(fit)))
- if(length(fit$df.prior)!=1) fit$df.prior <- fit$df.prior[genes]
+
fit[genes,]
}
diff --git a/R/geneset-roast.R b/R/geneset-roast.R
index 0445ba2..80fcca5 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 22 Jan 2013.
+# Created 24 Apr 2008. Last modified 19 Sep 2013.
{
# Check y
y <- as.matrix(y)
@@ -45,6 +45,7 @@ roast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.sta
n <- ncol(y)
# Check index
+ if(is.list(index)) return(mroast(y=y,index=index,design=design,contrast=contrast,set.statistic=set.statistic,gene.weights=gene.weights,array.weights=array.weights,weights=weights,block=block,correlation=correlation,var.prior=var.prior,df.prior=df.prior,trend.var=trend.var,nrot=nrot))
if(is.null(index)) index <- rep.int(TRUE,ngenes)
# Check design
@@ -104,13 +105,14 @@ roast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.sta
}
# 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
+# 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)
- sign2 <- sign(qr$qr[p,p])
+ signc <- sign(qr$qr[p,p])
if(is.null(var.prior) || is.null(df.prior)) {
# Fit model to all genes
@@ -119,15 +121,15 @@ roast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.sta
} else {
ws <- sqrt(weights)
effects <- matrix(0,n,ngenes)
- sign2 <- rep.int(0,ngenes)
+ signc <- rep.int(0,ngenes)
for (g in 1:ngenes) {
wX <- X*ws[g,]
wy <- y[g,]*ws[g,]
qrX <- qr(wX)
- sign2[g] <- sign(qrX$qr[p,p])
+ signc[g] <- sign(qrX$qr[p,p])
effects[,g] <- qr.qty(qrX,wy)
}
- sign2 <- sign2[index]
+ signc <- signc[index]
}
# Estimate global parameters s0 and d0
s2 <- colMeans(effects[-(1:p),,drop=FALSE]^2)
@@ -152,12 +154,12 @@ roast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.sta
ws <- sqrt(weights[index,,drop=FALSE])
nset <- nrow(y)
effects <- matrix(0,n,nset)
- sign2 <- rep.int(0,nset)
+ signc <- rep.int(0,nset)
for (g in 1:nset) {
wX <- X*ws[g,]
wy <- y[g,]*ws[g,]
qrX <- qr(wX)
- sign2[g] <- sign(qrX$qr[p,p])
+ signc[g] <- sign(qrX$qr[p,p])
effects[,g] <- qr.qty(qrX,wy)
}
}
@@ -176,7 +178,6 @@ roast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.sta
Y <- effects
YY <- colSums(Y^2)
B <- Y[1,]
- signc <- sign1*sign2
modt <- signc*B/sd.post
statobs <- p <- rep(0,3)
diff --git a/R/lmfit.R b/R/lmfit.R
index 902c349..1ea3104 100644
--- a/R/lmfit.R
+++ b/R/lmfit.R
@@ -1,25 +1,34 @@
# LINEAR MODELS
lmFit <- function(object,design=NULL,ndups=1,spacing=1,block=NULL,correlation,weights=NULL,method="ls",...)
-# Fit linear model
+# Fit genewise linear models
# Gordon Smyth
-# 30 June 2003. Last modified 6 Feb 2009.
+# 30 June 2003. Last modified 28 July 2013.
{
-# Check arguments
+# Extract components from y
y <- getEAWP(object)
+
+# Check design matrix
+ if(is.null(design)) design <- y$design
if(is.null(design))
- if(is.null(y$design))
- design <- matrix(1,ncol(y$exprs),1)
- else
- design <- as.matrix(object$design)
- else
+ design <- matrix(1,ncol(y$exprs),1)
+ else {
design <- as.matrix(design)
- if(mode(design) != "numeric") stop("design must be a numeric matrix")
+ if(mode(design) != "numeric") stop("design must be a numeric matrix")
+ }
ne <- nonEstimable(design)
if(!is.null(ne)) cat("Coefficients not estimable:",paste(ne,collapse=" "),"\n")
+
+# Weights and spacing arguments can be specified in call or stored in y
+# Precedence for these arguments is
+# 1. Specified in function call
+# 2. Stored in object
+# 3. Default values
if(missing(ndups) && !is.null(y$printer$ndups)) ndups <- y$printer$ndups
if(missing(spacing) && !is.null(y$printer$spacing)) spacing <- y$printer$spacing
if(missing(weights) && !is.null(y$weights)) weights <- y$weights
+
+# Check method
method <- match.arg(method,c("ls","robust"))
# If duplicates are present, reduce probe-annotation and Amean to correct length
@@ -28,7 +37,7 @@ lmFit <- function(object,design=NULL,ndups=1,spacing=1,block=NULL,correlation,we
if(!is.null(y$Amean)) y$Amean <- rowMeans(unwrapdups(as.matrix(y$Amean),ndups=ndups,spacing=spacing),na.rm=TRUE)
}
-# Dispath fitting algorithms
+# Dispatch fitting algorithms
if(method=="robust")
fit <- mrlm(y$exprs,design=design,ndups=ndups,spacing=spacing,weights=weights,...)
else
@@ -56,31 +65,43 @@ lmFit <- function(object,design=NULL,ndups=1,spacing=1,block=NULL,correlation,we
}
lm.series <- function(M,design=NULL,ndups=1,spacing=1,weights=NULL)
-{
# Fit linear model for each gene to a series of arrays
# Gordon Smyth
-# 18 Apr 2002. Revised 10 July 2008.
-
+# 18 Apr 2002. Revised 16 Apr 2013.
+{
+# Check expression matrix
M <- as.matrix(M)
narrays <- ncol(M)
- if(is.null(design)) design <- matrix(1,narrays,1)
- design <- as.matrix(design)
+
+# Check design
+ if(is.null(design))
+ design <- matrix(1,narrays,1)
+ else
+ design <- as.matrix(design)
nbeta <- ncol(design)
coef.names <- colnames(design)
+ if(is.null(coef.names)) coef.names <- paste("x",1:nbeta,sep="")
+
+# Check weights
if(!is.null(weights)) {
weights <- asMatrixWeights(weights,dim(M))
weights[weights <= 0] <- NA
M[!is.finite(weights)] <- NA
}
+
+# Reform duplicated rows into columns
if(ndups>1) {
M <- unwrapdups(M,ndups=ndups,spacing=spacing)
design <- design %x% rep(1,ndups)
if(!is.null(weights)) weights <- unwrapdups(weights,ndups=ndups,spacing=spacing)
}
+
+# Initialize standard errors
ngenes <- nrow(M)
stdev.unscaled <- beta <- matrix(NA,ngenes,nbeta,dimnames=list(rownames(M),coef.names))
- sigma <- rep(NA,ngenes)
- df.residual <- rep(0,ngenes)
+
+# Check whether QR-decomposition is constant for all genes
+# If so, fit all genes in one sweep
NoProbeWts <- !any(is.na(M)) && (is.null(weights) || !is.null(attr(weights,"arrayweights")))
if(NoProbeWts) {
if(is.null(weights))
@@ -105,6 +126,11 @@ lm.series <- function(M,design=NULL,ndups=1,spacing=1,weights=NULL)
fit$pivot <- fit$qr$pivot
return(fit)
}
+
+# Genewise QR-decompositions are required, so iteration through genes
+ beta <- stdev.unscaled
+ sigma <- rep(NA,ngenes)
+ df.residual <- rep(0,ngenes)
for (i in 1:ngenes) {
y <- as.vector(M[i,])
obs <- is.finite(y)
@@ -128,10 +154,13 @@ lm.series <- function(M,design=NULL,ndups=1,spacing=1,weights=NULL)
sigma[i] <- sqrt(sum(w*out$residuals^2)/out$df.residual)
}
}
+
+# Correlation matrix of coefficients
QR <- qr(design)
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)
}
@@ -359,7 +388,7 @@ getEAWP <- function(object)
# Given any microarray data object, extract basic information needed for
# linear modelling.
# Gordon Smyth
-# 9 March 2008. Last modified 9 Feb 2012.
+# 9 March 2008. Last modified 26 June 2013.
{
y <- list()
@@ -375,22 +404,17 @@ getEAWP <- function(object)
}
y$weights <- object$weights
y$probes <- object$genes
- if(is.null(y$probes) && !is.null(rownames(y$exprs))) y$probes <- data.frame(ID=rownames(y$exprs),stringsAsFactors=FALSE)
y$design <- object$design
} else {
if(is(object,"ExpressionSet")) {
y$exprs <- exprs(object)
- if(length(object at featureData@data))
- y$probes <- object at featureData@data
- else
- y$probes <- data.frame(ID=rownames(y$exprs),stringsAsFactors=FALSE)
+ if(length(object at featureData@data)) y$probes <- object at featureData@data
y$Amean <- rowMeans(y$exprs,na.rm=TRUE)
} else {
if(is(object,"PLMset")) {
y$exprs <- object at chip.coefs
if(length(y$exprs)==0) stop("chip.coefs has length zero")
if(length(object at se.chip.coefs)) y$weights <- 1/pmax(object at se.chip.coefs,1e-5)^2
- if(!is.null(rownames(y$exprs))) y$probes <- data.frame(ID=rownames(y$exprs),stringsAsFactors=FALSE)
y$Amean <- rowMeans(y$exprs,na.rm=TRUE)
} else {
if(is(object,"marrayNorm")) {
@@ -404,11 +428,33 @@ getEAWP <- function(object)
} else {
# Default method for matrices, data.frames, vsn objects etc.
y$exprs <- as.matrix(object)
- if(!is.null(rownames(y$exprs))) y$probes <- data.frame(ID=rownames(y$exprs),stringsAsFactors=FALSE)
# If exprs are positive, assume they are log-intensities rather than log-ratios
- if(all(y$exprs>=0,na.rm=TRUE)) y$Amean <- rowMeans(y$exprs,na.rm=TRUE)
+# if(all(y$exprs>=0,na.rm=TRUE)) y$Amean <- rowMeans(y$exprs,na.rm=TRUE)
+ 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")
+
+# Get rownames from probes?
+ if(is.null(rownames(y$exprs)) && !is.null(row.names(y$probes))) rownames(y$exprs) <- row.names(y$probes)
+
+# Check rownames are unique
+# rn <- rownames(y$exprs)
+# if(is.null(rn))
+# rownames(y$exprs) <- 1:nrow(y$exprs)
+# else
+# if(anyDuplicated(rn)>0) {
+# rownames(y$exprs) <- 1:nrow(y$exprs)
+# if(is.null(y$probes))
+# y$probes <- data.frame(ID=rn,stringsAsFactors=FALSE)
+# else
+# if("ID" %in% names(y$probes))
+# y$probes$ID0 <- rn
+# else
+# y$probes$ID <- rn
+# }
+
y
}
diff --git a/R/loessFit.R b/R/loessFit.R
index 84b6ec3..788c8be 100644
--- a/R/loessFit.R
+++ b/R/loessFit.R
@@ -1,11 +1,10 @@
# LOESS FUNCTIONS
-loessFit <- function(y, x, weights=NULL, span=0.3, bin=NULL, iterations=4, min.weight=1e-5, max.weight=1e5, equal.weights.as.null=TRUE)
-# Fast loess fit for simple x and y
-# This function uses stats:::lowess if no weights and stats:::loess otherwise.
-# It is intended to give a streamlined common interface to the two functions.
+loessFit <- function(y, x, weights=NULL, span=0.3, iterations=4, min.weight=1e-5, max.weight=1e5, equal.weights.as.null=TRUE)
+# Fast lowess fit for univariate x and y allowing for weights
+# Uses lowess() if weights=NULL and locfit.raw() otherwise
# Gordon Smyth
-# 28 June 2003. Last revised 8 Oct 2012.
+# 28 June 2003. Last revised 2 Oct 2013.
{
n <- length(y)
out <- list(fitted=rep(NA,n),residuals=rep(NA,n))
@@ -15,22 +14,11 @@ loessFit <- function(y, x, weights=NULL, span=0.3, bin=NULL, iterations=4, min.w
nobs <- length(yobs)
if(nobs==0) return(out)
if(is.null(weights)) {
- iter <- iterations-1
- if(is.null(bin)) bin <- 0.01
- delta = bin * diff(range(xobs))
o <- order(xobs)
-# The .C("lowess" call is copied from stats:::lowess
-# lo <- .C("lowess", x = as.double(xobs[o]), as.double(yobs[o]),
-# nobs, as.double(span), as.integer(iter), as.double(delta),
-# y = double(nobs), double(nobs), double(nobs), PACKAGE = "stats")
-# out$fitted[obs][o] <- lo$y
-# out$residuals[obs][o] <- lo[[9]]
-# For R 2.16.X
- lo <- lowess(x=xobs,y=yobs,f=span,iter=iter,delta=delta)
+ lo <- lowess(x=xobs,y=yobs,f=span,iter=iterations-1)
out$fitted[obs][o] <- lo$y
out$residuals[obs] <- yobs-out$fitted[obs]
} else {
- if(is.null(bin)) bin <- 0.005
wobs <- weights[obs]
wobs[is.na(wobs)] <- 0
wobs <- pmax(wobs,min.weight)
@@ -38,7 +26,7 @@ loessFit <- function(y, x, weights=NULL, span=0.3, bin=NULL, iterations=4, min.w
# Test whether weights are equal
if(equal.weights.as.null) {
r <- range(wobs)
- if(r[2]-r[1] < 1e-15) return(Recall(y,x,span=span,bin=bin,iterations=iterations))
+ if(r[2]-r[1] < 1e-15) return(Recall(y=y,x=x,span=span,iterations=iterations))
}
# Count number of observations with positive weights
if(min.weight>0)
@@ -54,16 +42,24 @@ loessFit <- function(y, x, weights=NULL, span=0.3, bin=NULL, iterations=4, min.w
fit$residuals <- yobs-fit$fitted
}
} else {
-# Suppress warning "k-d tree limited by memory"
-# oldopt <- options(warning.expression=expression())
- oldopt <- options(warn=-1)
- on.exit(options(oldopt))
-# fit <- .vsimpleLoess(y=yobs, x=xobs, weights=wobs, span=span, degree=1, cell=bin/span, iterations=iterations)
-# For R 2.16.X
- fit <- stats:::simpleLoess(y=yobs,x=xobs,weights=wobs,span=span,degree=1,parametric=FALSE,normalize=FALSE,statistics="none",surface="interpolate",cell=bin/span,iterations=iterations,trace.hat="approximate")
+# Check for locfit package
+ loaded <- ( "package:locfit" %in% search() )
+ if(!loaded) {
+ loadresult <- tryCatch(suppressPackageStartupMessages(library("locfit",character.only=TRUE,quietly=TRUE)),error=function(e) e)
+ if(inherits(loadresult,"error")) stop("limma:::loessFit with weights requires locfit package, which is not available",call.=FALSE)
+ }
+
+# Weighted lowess with robustifying iterations
+ biweights <- rep(1,nobs)
+ for (i in 1:iterations) {
+ fit <- locfit.raw(x=xobs,y=yobs,weights=wobs*biweights,alpha=span,deg=1)
+ res <- residuals(fit,type="raw")
+ s <- median(abs(res))
+ biweights <- pmax(1-(res/(6*s))^2,0)^2
+ }
}
- out$fitted[obs] <- fit$fitted
- out$residuals[obs] <- fit$residuals
+ out$fitted[obs] <- fitted(fit)
+ out$residuals[obs] <- res
}
out
}
diff --git a/R/norm.R b/R/norm.R
index 5f6a9f8..1e81cd6 100755
--- a/R/norm.R
+++ b/R/norm.R
@@ -332,7 +332,7 @@ plotPrintorder <- function(object,layout,start="topleft",slide=1,method="loess",
normalizeBetweenArrays <- function(object, method=NULL, targets=NULL, cyclic.method="fast", ...) {
# Normalize between arrays
# Gordon Smyth
-# 12 Apri 2003. Last revised 2 Feb 2012.
+# 12 Apri 2003. Last revised 21 July 2013.
# Default method
if(is.null(method)) {
@@ -350,7 +350,7 @@ normalizeBetweenArrays <- function(object, method=NULL, targets=NULL, cyclic.met
# Methods for matrices
if(is(object,"matrix")) {
- if(!(method %in% c("none","scale","quantile","cyclicloess"))) stop("method not applicable to matrix objects")
+ if(!(method %in% c("none","scale","quantile","cyclicloess"))) stop("method not applicable to single-channel data")
return(switch(method,
none = object,
scale = normalizeMedianValues(object),
diff --git a/R/plots-ma.R b/R/plots-ma.R
index 9e2ed63..76fcdea 100755
--- a/R/plots-ma.R
+++ b/R/plots-ma.R
@@ -1,20 +1,102 @@
-# M-A PLOTS
+## PLOTS-MA.R
+## M-A PLOTS
-plotMA <- function(MA, array=1, xlab="A", ylab="M", main=colnames(MA)[array], xlim=NULL, ylim=NULL, status, values, pch, col, cex, legend=TRUE, zero.weights=FALSE, ...)
+plotMA <- function(MA,...) UseMethod("plotMA")
+
+plotMA.RGList <- function(MA, array=1, xlab="A", ylab="M", main=colnames(MA)[array], xlim=NULL, ylim=NULL, status, values, pch, col, cex, legend=TRUE, zero.weights=FALSE, ...)
# MA-plot with color coding for controls
# Gordon Smyth 7 April 2003, James Wettenhall 27 June 2003.
# Last modified 23 April 2013.
{
-# Convert to MAList if possible
- if(class(MA)=="list") MA <- new("MAList",MA)
- if(is(MA,"RGList")) {
- MA <- MA.RG(MA[,array])
- array <- 1
+ MA <- MA.RG(MA[,array])
+ plotMA(MA,array=1,xlab=ylab,ylab=ylab,main=main,xlim=xlim,ylim=ylim,status=status,values=values,pch=pch,col=col,cex=cex,legend=legend,zero.weights=zero.weights,...)
+}
+
+plotMA.MAList <- function(MA, array=1, xlab="A", ylab="M", main=colnames(MA)[array], xlim=NULL, ylim=NULL, status, values, pch, col, cex, legend=TRUE, zero.weights=FALSE, ...)
+# MA-plot with color coding for controls
+# Gordon Smyth 7 April 2003, James Wettenhall 27 June 2003.
+# Last modified 23 April 2013.
+{
+ x <- as.matrix(MA$A)[,array]
+ y <- as.matrix(MA$M)[,array]
+ if(is.null(MA$weights)) w <- NULL else w <- as.matrix(MA$weights)[,array]
+ if(missing(status)) status <- MA$genes$Status
+ if(!is.null(w) && !zero.weights) {
+ i <- is.na(w) | (w <= 0)
+ y[i] <- NA
+ }
+ .plotMAxy(x,y,xlab=ylab,ylab=ylab,main=main,xlim=xlim,ylim=ylim,status=status,values=values,pch=pch,col=col,cex=cex,legend=legend, ...)
+}
+
+plotMA.MArrayLM <- function(MA, array=1, xlab="A", ylab="M", main=colnames(MA)[array], xlim=NULL, ylim=NULL, status, values, pch, col, cex, legend=TRUE, zero.weights=FALSE, ...)
+# MA-plot with color coding for controls
+# Gordon Smyth 7 April 2003, James Wettenhall 27 June 2003.
+# Last modified 23 April 2013.
+{
+ if(is.null(MA$Amean)) stop("MA-plot not possible because Amean component is absent.")
+ x <- MA$Amean
+ y <- as.matrix(MA$coef)[,array]
+ if(is.null(MA$weights)) w <- NULL else w <- as.matrix(MA$weights)[,array]
+ if(missing(status)) status <- MA$genes$Status
+ if(!is.null(w) && !zero.weights) {
+ i <- is.na(w) | (w <= 0)
+ y[i] <- NA
}
- if(is(MA,"EListRaw")) {
- MA <- normalizeBetweenArrays(MA,method="none")
+ .plotMAxy(x,y,xlab=ylab,ylab=ylab,main=main,xlim=xlim,ylim=ylim,status=status,values=values,pch=pch,col=col,cex=cex,legend=legend, ...)
+}
+
+plotMA.EList <- function(MA, array=1, xlab="A", ylab="M", main=colnames(MA)[array], xlim=NULL, ylim=NULL, status, values, pch, col, cex, legend=TRUE, zero.weights=FALSE, ...)
+# MA-plot with color coding for controls
+# Gordon Smyth 7 April 2003, James Wettenhall 27 June 2003.
+# Last modified 23 April 2013.
+{
+ MA$E <- as.matrix(MA$E)
+ narrays <- ncol(MA$E)
+ if(narrays < 2) stop("Need at least two arrays")
+ if(narrays > 5)
+ x <- apply(MA$E,1,median,na.rm=TRUE)
+ else
+ x <- rowMeans(MA$E,na.rm=TRUE)
+ y <- MA$E[,array]-x
+ if(is.null(MA$weights)) w <- NULL else w <- as.matrix(MA$weights)[,array]
+ if(missing(status)) status <- MA$genes$Status
+
+ if(!is.null(w) && !zero.weights) {
+ i <- is.na(w) | (w <= 0)
+ y[i] <- NA
+ }
+ .plotMAxy(x,y,xlab=ylab,ylab=ylab,main=main,xlim=xlim,ylim=ylim,status=status,values=values,pch=pch,col=col,cex=cex,legend=legend, ...)
+}
+
+plotMA.default <- function(MA, array=1, xlab="A", ylab="M", main=colnames(MA)[array], xlim=NULL, ylim=NULL, status, values, pch, col, cex, legend=TRUE, zero.weights=FALSE, ...)
+# MA-plot with color coding for controls
+# Gordon Smyth 7 April 2003, James Wettenhall 27 June 2003.
+# Last modified 23 April 2013.
+{
+# Data is assumed to be single-channel
+ MA <- as.matrix(MA)
+ narrays <- ncol(MA)
+ if(narrays < 2) stop("Need at least two arrays")
+ if(narrays > 5)
+ x <- apply(MA,1,median,na.rm=TRUE)
+ else
+ x <- rowMeans(MA,na.rm=TRUE)
+ y <- MA[,array]-x
+ w <- NULL
+ if(missing(status)) status <- NULL
+
+ if(!is.null(w) && !zero.weights) {
+ i <- is.na(w) | (w <= 0)
+ y[i] <- NA
}
+ .plotMAxy(x,y,xlab=ylab,ylab=ylab,main=main,xlim=xlim,ylim=ylim,status=status,values=values,pch=pch,col=col,cex=cex,legend=legend, ...)
+}
+.plotMAxy <- function(x, y, xlab="A", ylab="M", main=NULL, xlim=NULL, ylim=NULL, status, values, pch, col, cex, legend=TRUE, ...)
+# MA-plot with color coding for controls
+# Gordon Smyth 7 April 2003, James Wettenhall 27 June 2003.
+# Last modified 23 April 2013.
+{
# Check legend
legend.position <- "topleft"
if(!is.logical(legend)) {
@@ -23,99 +105,74 @@ plotMA <- function(MA, array=1, xlab="A", ylab="M", main=colnames(MA)[array], xl
}
legend.position <- match.arg(legend.position,c("bottomright","bottom","bottomleft","left","topleft","top","topright","right","center"))
- if(is(MA,"MAList")) {
-# Data is two-color
- x <- as.matrix(MA$A)[,array]
- y <- as.matrix(MA$M)[,array]
- if(is.null(MA$weights)) w <- NULL else w <- as.matrix(MA$weights)[,array]
- if(missing(status)) status <- MA$genes$Status
- } else if(is(MA,"MArrayLM")) {
- if(is.null(MA$Amean)) stop("MA-plot not possible because Amean component is absent.")
- x <- MA$Amean
- y <- as.matrix(MA$coef)[,array]
- if(is.null(MA$weights)) w <- NULL else w <- as.matrix(MA$weights)[,array]
- if(missing(status)) status <- MA$genes$Status
- } else if(is(MA,"EList")) {
- MA$E <- as.matrix(MA$E)
- narrays <- ncol(MA$E)
- if(narrays < 2) stop("Need at least two arrays")
- if(narrays > 5)
- x <- apply(MA$E,1,median,na.rm=TRUE)
- else
- x <- rowMeans(MA$E,na.rm=TRUE)
- y <- MA$E[,array]-x
- if(is.null(MA$weights)) w <- NULL else w <- as.matrix(MA$weights)[,array]
- if(missing(status)) status <- MA$genes$Status
- } else {
-# Data is assumed to be single-channel
- MA <- as.matrix(MA)
- narrays <- ncol(MA)
- if(narrays < 2) stop("Need at least two arrays")
- if(narrays > 5)
- x <- apply(MA,1,median,na.rm=TRUE)
- else
- x <- rowMeans(MA,na.rm=TRUE)
- y <- MA[,array]-x
- w <- NULL
- if(missing(status)) status <- NULL
- }
-
- if(!is.null(w) && !zero.weights) {
- i <- is.na(w) | (w <= 0)
- y[i] <- NA
- }
+# Check xlim and ylim
if(is.null(xlim)) xlim <- range(x,na.rm=TRUE)
if(is.null(ylim)) ylim <- range(y,na.rm=TRUE)
+
+# Setup plot axes
plot(x,y,xlab=xlab,ylab=ylab,main=main,xlim=xlim,ylim=ylim,type="n",...)
+
+# If no status information, just plot points normally
if(is.null(status) || all(is.na(status))) {
- if(missing(pch)) pch=16
- if(missing(cex)) cex=0.3
+ if(missing(pch)) pch <- 16
+ if(missing(cex)) cex <- 0.3
points(x,y,pch=pch[[1]],cex=cex[1])
- } else {
- if(missing(values)) {
- if(is.null(attr(status,"values")))
- values <- names(sort(table(status),decreasing=TRUE))
- else
- values <- attr(status,"values")
- }
-# Non-highlighted points
- sel <- !(status %in% values)
- nonhi <- any(sel)
- if(nonhi) points(x[sel],y[sel],pch=16,cex=0.3)
-
- nvalues <- length(values)
- if(missing(pch)) {
- if(is.null(attr(status,"pch")))
- pch <- rep(16,nvalues)
- else
- pch <- attr(status,"pch")
- }
- if(missing(cex)) {
- if(is.null(attr(status,"cex"))) {
- cex <- rep(1,nvalues)
- if(!nonhi) cex[1] <- 0.3
- } else
- cex <- attr(status,"cex")
- }
- if(missing(col)) {
- if(is.null(attr(status,"col"))) {
- col <- nonhi + 1:nvalues
- } else
- col <- attr(status,"col")
- }
- pch <- rep(pch,length=nvalues)
- col <- rep(col,length=nvalues)
- cex <- rep(cex,length=nvalues)
- for (i in 1:nvalues) {
- sel <- status==values[i]
- points(x[sel],y[sel],pch=pch[[i]],cex=cex[i],col=col[i])
- }
- if(legend) {
- if(is.list(pch))
- legend(legend.position,legend=values,fill=col,col=col,cex=0.9)
- else
- legend(legend.position,legend=values,pch=pch,,col=col,cex=0.9)
- }
+ return(invisible())
+ }
+
+# From here, status is not NULL and not all missing
+
+# Check values
+ if(missing(values)) {
+ if(is.null(attr(status,"values")))
+ values <- names(sort(table(status),decreasing=TRUE))
+ else
+ values <- attr(status,"values")
+ }
+ nvalues <- length(values)
+
+# Plot non-highlighted points
+ sel <- !(status %in% values)
+ nonhi <- any(sel)
+ if(nonhi) points(x[sel],y[sel],pch=16,cex=0.3)
+
+ if(missing(pch)) {
+ if(is.null(attr(status,"pch")))
+ pch <- rep(16,nvalues)
+ else
+ pch <- attr(status,"pch")
+ }
+
+ if(missing(cex)) {
+ if(is.null(attr(status,"cex"))) {
+ cex <- rep(1,nvalues)
+ if(!nonhi) cex[1] <- 0.3
+ } else
+ cex <- attr(status,"cex")
+ }
+
+ if(missing(col)) {
+ if(is.null(attr(status,"col"))) {
+ col <- nonhi + 1:nvalues
+ } else
+ col <- attr(status,"col")
+ }
+
+ pch <- rep(pch,length=nvalues)
+ col <- rep(col,length=nvalues)
+ cex <- rep(cex,length=nvalues)
+
+# Plot highlighted classes of points
+ for (i in 1:nvalues) {
+ sel <- status==values[i]
+ points(x[sel],y[sel],pch=pch[[i]],cex=cex[i],col=col[i])
+ }
+
+ if(legend) {
+ if(is.list(pch))
+ legend(legend.position,legend=values,fill=col,col=col,cex=0.9)
+ else
+ legend(legend.position,legend=values,pch=pch,,col=col,cex=0.9)
}
invisible()
}
diff --git a/R/predFCm.R b/R/predFCm.R
index acaaadd..e2fe0c1 100644
--- a/R/predFCm.R
+++ b/R/predFCm.R
@@ -1,7 +1,7 @@
predFCm <- function(fit,coef=2,prob=TRUE,VarRel=NULL)
# Belinda Phipson 29 May 2012. Updated 8 January 2013.
{
- p <- 1-propTrueNull(fit$p.value[,coef])
+ p <- 1-propTrueNull(fit$p.value[,coef],method="lfdr")
if(p==0) p<-1e-8
if(length(fit$s2.prior)==1) trend<-FALSE else trend<-TRUE
fit <- eBayes(fit,proportion = p,trend=trend)
diff --git a/R/read-ilmn.R b/R/read-ilmn.R
index 204fee4..5cfdf59 100644
--- a/R/read-ilmn.R
+++ b/R/read-ilmn.R
@@ -68,28 +68,41 @@ read.ilmn.targets <- function(targets, ...)
h <- readGenericHeader(fname,columns=expr)
skip <- h$NHeaderRecords
header <- h$ColumnNames
-
+
elist <- new("EListRaw")
elist$source <- "illumina"
reqcol <- header[grep(tolower(paste(c(probeid, annotation, expr, other.columns), collapse="|")), tolower(header))]
reqcol <- trimWhiteSpace(reqcol)
-
+
x <- read.columns(file=fname, required.col=reqcol, skip=skip, sep=sep, quote=quote, stringsAsFactors=FALSE, ...)
nprobes <- nrow(x)
-
- pids <- x[, grep(tolower(probeid), tolower(colnames(x)))]
- snames <- colnames(x)[grep(tolower(expr), tolower(colnames(x)))]
+
+# Match column names to find column numbers
+ cn <- tolower(colnames(x))
+ idcol <- grep(tolower(probeid), cn)
+ anncol <- grep(tolower(paste(annotation,collapse="|")), cn)
+ exprcol <- grep(tolower(expr), cn)
+
+# Probe IDs
+ pids <- x[,idcol]
+
+# Sample names
+ snames <- colnames(x)[exprcol]
snames <- unlist(strsplit(snames, paste("[.]*", expr, "-*", sep="")))
snames <- snames[snames != ""]
nsamples <- length(snames)
-
- elist$E <- data.matrix(x[, grep(tolower(expr), tolower(colnames(x)))])
+
+# Expression matrix
+ elist$E <- data.matrix(x[,exprcol])
colnames(elist$E) <- snames
rownames(elist$E) <- pids
-
- elist$genes <- x[, c(grep(tolower(probeid), tolower(colnames(x))), grep(tolower(paste(annotation,collapse="|")), tolower(colnames(x)))), drop=FALSE]
- elist$targets <- data.frame(SampleNames=snames, stringsAsFactors=FALSE)
-
+
+# Add probe annotation
+ if(length(anncol)) elist$genes <- x[,anncol,drop=FALSE]
+
+# elist$targets <- data.frame(SampleNames=snames, stringsAsFactors=FALSE)
+
+# Add other columns if required
if(!is.null(other.columns)){
elist$other <- vector("list", length(other.columns))
for(i in 1:length(other.columns)){
@@ -102,5 +115,6 @@ read.ilmn.targets <- function(targets, ...)
}
names(elist$other) <- other.columns
}
+
elist
}
diff --git a/R/subsetting.R b/R/subsetting.R
index 835d5d5..4cb97a8 100755
--- a/R/subsetting.R
+++ b/R/subsetting.R
@@ -151,7 +151,7 @@ assign("[.MArrayLM",
function(object, i, j, ...)
# Subsetting for MArrayLM objects
# Gordon Smyth
-# 26 April 2005. Last modified 13 January 2010.
+# 26 April 2005. Last modified 28 September 2013.
{
if(nargs() != 3) stop("Two subscripts required",call.=FALSE)
if(!is.null(object$coefficients)) object$coefficients <- as.matrix(object$coefficients)
@@ -208,6 +208,7 @@ function(object, i, 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]
diff --git a/R/toptable.R b/R/toptable.R
index 8efae1c..1bc60f4 100644
--- a/R/toptable.R
+++ b/R/toptable.R
@@ -37,48 +37,82 @@ topTable <- function(fit,coef=NULL,number=10,genelist=fit$genes,adjust.method="B
topTableF <- function(fit,number=10,genelist=fit$genes,adjust.method="BH",sort.by="F",p.value=1,lfc=0)
# Summary table of top genes by F-statistic
# Gordon Smyth
-# 27 August 2006. Last modified 19 Jan 2012.
+# 27 August 2006. Last modified 22 April 2013.
{
-# Check input
- if(is.null(fit$F)) stop("F-statistics not available. Try topTable for individual coef instead.")
+# Check fit
+ if(is.null(fit$coefficients)) stop("Coefficients not found in fit")
+ M <- as.matrix(fit$coefficients)
+ rn <- rownames(M)
+ if(is.null(colnames(M))) colnames(M) <- paste("Coef",1:ncol(M),sep="")
+ Amean <- fit$Amean
+ Fstat <- fit$F
+ Fp <- fit$F.p.value
+ if(is.null(Fstat)) stop("F-statistics not found in fit")
+
+# Ensure genelist is a data.frame
if(!is.null(genelist) && is.null(dim(genelist))) genelist <- data.frame(ProbeID=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("F","none"))
- M <- as.matrix(fit$coefficients)
- if(is.null(colnames(M))) colnames(M) <- paste("Coef",1:ncol(M),sep="")
- rn <- rownames(M)
- if(is.null(rn) || length(rn) > length(unique(rn))) rn <- 1:nrow(M)
- adj.P.Value <- p.adjust(fit$F.p.value,method=adjust.method)
-
- if(lfc > 0) {
- big <- apply(abs(M)>lfc,1,any)
- big[is.na(big)] <- FALSE
- if(any(!big)) {
- fit <- fit[big,]
- if(!is.null(genelist)) genelist <- genelist[big,,drop=FALSE]
- M <- M[big,]
- rn <- rn[big]
- adj.P.Value <- adj.P.Value[big]
+# Apply multiple testing adjustment
+ adj.P.Value <- p.adjust(Fp,method=adjust.method)
+
+# Thin out fit by lfc and p.value thresholds
+ if(lfc > 0 || p.value < 1) {
+ if(lfc>0)
+ big <- rowSums(abs(M)>lfc,na.rm=TRUE)>0
+ else
+ big <- TRUE
+ if(p.value<1) {
+ sig <- adj.P.Value <= p.value
+ sig[is.na(sig)] <- FALSE
+ } else
+ sig <- TRUE
+ keep <- big & sig
+ if(!all(keep)) {
+ M <- M[keep,,drop=FALSE]
+ rn <- rn[keep]
+ Amean <- Amean[keep]
+ Fstat <- Fstat[keep]
+ Fp <- p.value[keep]
+ genelist <- genelist[keep,,drop=FALSE]
+ adj.P.Value <- adj.P.Value[keep]
}
}
+
+# Enough rows left?
if(nrow(fit) < number) number <- nrow(fit)
+ if(number < 1) return(data.frame())
+# Find rows of top genes
if(sort.by=="F")
- o <- order(fit$F.p.value,decreasing=FALSE)[1:number]
+ o <- order(Fp,decreasing=FALSE)[1:number]
else
o <- 1:number
- if(p.value < 1) {
- nsig <- sum(adj.P.Value[o] <= p.value,na.rm=TRUE)
- if(nsig < number) o <- o[1:nsig]
- }
-
+# Assemble data.frame
if(is.null(genelist))
tab <- data.frame(M[o,,drop=FALSE])
else
tab <- data.frame(genelist[o,,drop=FALSE],M[o,,drop=FALSE])
- if(!is.null(fit$Amean)) tab <- data.frame(tab,AveExpr=fit$Amean[o])
- tab <- data.frame(tab,F=fit$F[o],P.Value=fit$F.p.value[o],adj.P.Val=adj.P.Value[o])
+ tab$AveExpr=fit$Amean[o]
+ tab <- data.frame(tab,F=Fstat[o],P.Value=Fp[o],adj.P.Val=adj.P.Value[o])
rownames(tab) <- rn[o]
tab
}
@@ -86,37 +120,77 @@ 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 30 Jan 2011.
+# 21 Nov 2002. Last revised 16 April 2013.
{
-# Check input
- if(length(coef)>1) coef <- coef[1]
+# Check fit
+ fit$coefficients <- as.matrix(fit$coefficients)
+ rn <- rownames(fit$coefficients)
+
+# 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)
- if(is.null(eb)) {
- fit$coefficients <- as.matrix(fit$coefficients)[,coef]
- fit$stdev.unscaled <- as.matrix(fit$stdev.unscaled)[,coef]
- eb <- ebayes(fit,...)
- coef <- 1
+
+# Check rownames
+ if(is.null(rn))
+ rn <- 1:nrow(fit$coefficients)
+ 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(fit$coefficients)
+ }
+
+# Check sort.by
+ sort.by <- match.arg(sort.by,c("logFC","M","A","Amean","AveExpr","P","p","T","t","B","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","B"))
+ 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"
}
- M <- as.matrix(fit$coefficients)[,coef]
+
+# Check A
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)
}
+
+# Compute eb if not given, compute just the one column required
+ if(is.null(eb)) {
+ fit$coefficients <- fit$coefficients[,coef,drop=FALSE]
+ fit$stdev.unscaled <- as.matrix(fit$stdev.unscaled)[,coef,drop=FALSE]
+ eb <- ebayes(fit,...)
+ coef <- 1
+ }
+
+# 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]
- rownum <- 1:length(M)
- sort.by <- match.arg(sort.by,c("logFC","M","A","Amean","AveExpr","P","p","T","t","B","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"
# Apply multiple testing adjustment
adj.P.Value <- p.adjust(P.Value,method=adjust.method)
-# Apply p.value or lfc thresholds
+# Thin out fit by p.value and lfc thresholds
if(p.value < 1 | lfc > 0) {
sig <- (adj.P.Value < p.value) & (abs(M) > lfc)
if(any(is.na(sig))) sig[is.na(sig)] <- FALSE
@@ -128,8 +202,14 @@ toptable <- function(fit,coef=1,number=10,genelist=NULL,A=NULL,eb=NULL,adjust.me
P.Value <- P.Value[sig]
adj.P.Value <- adj.P.Value[sig]
B <- B[sig]
- rownum <- rownum[sig]
+ rn <- rn[sig]
}
+
+# Are enough rows left?
+ if(length(M) < number) number <- length(M)
+ if(number < 1) return(data.frame())
+
+# Select top rows
ord <- switch(sort.by,
logFC=order(abs(M),decreasing=TRUE),
AveExpr=order(A,decreasing=TRUE),
@@ -138,9 +218,9 @@ toptable <- function(fit,coef=1,number=10,genelist=NULL,A=NULL,eb=NULL,adjust.me
B=order(B,decreasing=TRUE),
none=1:length(M)
)
- if(length(M) < number) number <- length(M)
- if(number < 1) return(data.frame())
top <- ord[1:number]
+
+# Assemble output data.frame
if(is.null(genelist))
tab <- data.frame(logFC=M[top])
else {
@@ -151,15 +231,12 @@ toptable <- function(fit,coef=1,number=10,genelist=NULL,A=NULL,eb=NULL,adjust.me
tab$CI.025 <- M[top]-margin.error
tab$CI.975 <- M[top]+margin.error
}
- if(!is.null(A)) tab <- data.frame(tab,AveExpr=A[top])
+ 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])
- rownames(tab) <- as.character(rownum)[top]
+ rownames(tab) <- rn[top]
+
+# Resort table
if(!is.null(resort.by)) {
- resort.by <- match.arg(resort.by,c("logFC","M","A","Amean","AveExpr","P","p","T","t","B"))
- 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"
ord <- switch(resort.by,
logFC=order(tab$logFC,decreasing=TRUE),
AveExpr=order(tab$AveExpr,decreasing=TRUE),
@@ -169,8 +246,7 @@ toptable <- function(fit,coef=1,number=10,genelist=NULL,A=NULL,eb=NULL,adjust.me
)
tab <- tab[ord,]
}
-# attr(tab,"coef") <- coef
-# attr(tab,"adjust.method") <- adjust.method
+
tab
}
diff --git a/R/treat.R b/R/treat.R
index 733a4a3..ebaf34e 100644
--- a/R/treat.R
+++ b/R/treat.R
@@ -53,26 +53,63 @@ treat <- function(fit, lfc=0, trend=FALSE)
topTreat <- function(fit,coef=1,number=10,genelist=fit$genes,adjust.method="BH",sort.by="p",resort.by=NULL,p.value=1)
# Summary table of top genes by treat
# Gordon Smyth
-# 15 June 2009. Last modified 17 March 2010.
+# 15 June 2009. Last modified 20 April 2013.
{
-# Check input
- if(length(coef)>1) coef <- coef[1]
- M <- as.matrix(fit$coefficients)[,coef]
+# 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)
}
- tstat <- as.matrix(fit$t)[,coef]
- P.Value <- as.matrix(fit$p.value)[,coef]
+
+# 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)
@@ -95,8 +132,12 @@ topTreat <- function(fit,coef=1,number=10,genelist=fit$genes,adjust.method="BH",
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])
@@ -105,13 +146,10 @@ topTreat <- function(fit,coef=1,number=10,genelist=fit$genes,adjust.method="BH",
}
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) <- as.character(1:length(M))[top]
+ rownames(tab) <- rn[top]
+
+# Resort
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"
ord <- switch(resort.by,
logFC=order(tab$logFC,decreasing=TRUE),
AveExpr=order(tab$AveExpr,decreasing=TRUE),
@@ -120,6 +158,7 @@ topTreat <- function(fit,coef=1,number=10,genelist=fit$genes,adjust.method="BH",
)
tab <- tab[ord,]
}
+
tab
}
diff --git a/build/vignette.rds b/build/vignette.rds
new file mode 100644
index 0000000..e3a034b
Binary files /dev/null and b/build/vignette.rds differ
diff --git a/inst/doc/changelog.txt b/inst/doc/changelog.txt
index a5f23de..5fbadff 100755
--- a/inst/doc/changelog.txt
+++ b/inst/doc/changelog.txt
@@ -1,22 +1,96 @@
-19 September 2013: limma 3.16.8
+11 Oct 2013: limma 3.17.28
-- bug fix to mroast() which was not using weights correctly.
+- Add comment to RNA-seq case study in User's Guide about conflict
+ with plotMA() function in BiocGenerics.
-29 July 2013: limma 3.16.7
+11 Oct 2013: limma 3.17.27
+
+- plotMA() is now an S3 generic function.
+
+2 Oct 2013: limma 3.17.26
+
+- Improved error message from loessFit() when locfit package is
+ needed but not available.
+
+29 Sept 2013: limma 3.17.25
+
+- loessFit() now uses the locfit.raw in the locfit package when
+ weights are provided instead of loess in the stats package. The
+ function now runs very efficiently even on very long data vectors.
+ The output results will change slightly when weights are provided.
+
+- Roles of contributors now specified in author field of DESCRIPTION
+ file using standard codes.
+
+- Bug fix to subsetting for MArrayLM objects: the df.total component
+ was not being subsetted.
+
+- Removing defunct Berkeley Press links to published Smyth (2004)
+ article in several Rd files. Replacing with link to Preprint.
+
+- Adding link to Phipson (2013) thesis in two Rd files.
+
+- Minor cleaning up of normexp.c source code. Remove unused variable
+ and keep constant representing infinity within floating point
+ range.
+
+19 Sept 2013: limma 3.17.24
+
+- roast() now calls mroast() if the index vector is a list.
+
+- Bug fix to mroast(), which had been ignoring the weights.
+
+1 Sept 2013: limma 3.17.23
+
+- New function contrastAsCoef(), which reforms a design matrix so
+ that one or more specified contrasts become coefficients. This
+ function is called by roast().
+
+22 August 2013: limma 3.17.22
+
+- Add Majewski et al reference to genas.Rd.
+
+- Add Phipson et al and Sartor et al references to squeezeVar.Rd.
+
+- Add Phipson et al reference to eBayes.Rd.
+
+3 August 2013: limma 3.17.21
+
+- fix typo in User's Guide.
+
+29 July 2013: limma 3.17.20
- duplicateCorrelation() now uses the weights matrix when block is
set. Previously the weights were ignored when block was used.
+26 July 2013: limma 3.17.19
+
+- genas() now returns NA results and a message when no genes satisfy
+ the criterion for inclusion in the analysis.
+
+26 July 2013: limma 3.17.18
+
+- argument chooseMethod for genas() renamed to subset. Option "n"
+ renamed to "all". Some editing of help page and streamlining of
+ code for genas().
+
+21 July 2013: limma 3.17.17
+
- edits to normalizeBetweenArrays help page (i) to further clarify
which normalization methods are available for single-channel data
and which are available for two-color data and (ii) to give a cross
citation to the neqc() function for Illumina BeadChips.
-14 July 2013: limma 3.16.6
+26 June 2013: limma 3.17.16
- bug fix to eBayes(robust=TRUE) when some of the df.prior values are
infinite.
+ 5 June 2013: limma 3.17.15
+
+- new function beadCountWeights() to estimate quantitative weights
+ from the bead counts for each probe for Illumina BeadArrays.
+
- voom() now outputs lib.size as a column of targets instead of as a
separate component.
@@ -26,39 +100,105 @@
- plotMDS() now checks explicitly that there are at least 3 samples
to plot.
-- The legend argument to plotMA() can now take a character value
- giving the position to place the legend.
+27 May 2013: limma 3.17.14
-- New merge methods for EList and EListRaw objects.
+- normexp.fit.detection.p() now tolerates some non-monotonicity in
+ the detection p-pvalues as a function of expression.
-- Fix minor typo in roast.Rd.
+27 May 2013: limma 3.17.13
-27 May 2013: limma 3.16.5
+- Update lmscFit and voom references.
-- normexp.fit.detection.p() now tolerates some non-monotonicity in
- the detection p-pvalues as a function of expression.
+- Update voomByGroup documentation.
+
+12 May 2013: limma 3.17.12
-- A number of help pages have been updated, including a new reference
- for lmscFit.
+- New references added for separate channel analysis and voom.
-17 May 2013: limma 3.16.4
+ 9 May 2013: limma 3.17.11
-- Update references on help page for voom().
+- New merge methods for EList and EListRaw objects.
- 1 May 2013: limma 3.16.3
+ 5 May 2013: limma 3.17.10
+
+- read.ilmn() no longer adds the probe IDs to the gene annotation
+ data.frame, leaving them instead as rownames of the expression
+ matrix. It longer creates a targets file since the sample names
+ are already preserved as column names of the expression matrix.
+
+- Update mammmary stem cell case study in User's Guide. As well
+ as reflecting changes to read.ilmn() and topTable(), this now
+ demonstrates how to find signature genes for particular cell type.
+
+ 1 May 2013: limma 3.17.9
- Bug fix to ebayes(), which was not passing the 'robust' argument
correctly on to squeezeVar().
-25 April 2013: limma 3.16.2
+- Update of mammary stem cell case study in User's Guide.
+
+25 April 2013: limma 3.17.8
+
+- Bug fix to fitFDistRobustly(), which affected the estimated scale
+ when df2 is estimated to be Inf.
+
+23 April 2013: limma 3.17.7
+
+- Update Estrogen case study in User's Guide to reflect change in
+ rownames of topTable output.
+
+- Update Nigerian RNA-seq case study in User's Guide, adding new
+ MA-plot.
+
+- The legend argument to plotMA() can now take a character value
+ giving the position to place the legend.
+
+22 April 2013: limma 3.17.6
+
+- Further fixes to toptable(), topTableF() and topTreat() to
+ implement the changes begun in 3.17.5.
+
+- documentation about rownames and column names and the use of
+ rownames(fit) and colnames(fit) added to lmFit.Rd.
+
+20 April 2013: limma 3.17.5
+
+- Reverting most of change to getEAWP() in version 3.17.3. getEAWP()
+ now simply preserves the rownames found in object, specifically
+ those found in the expression matrix. Null rownames are left null,
+ unless the probe annotation data.frame has row.names, in which case
+ these are used.
+
+- toptable(), topTable() and topTreat() now preserve the rownames of
+ the fit object, unless the fit object has duplicated names, in
+ which case the rownames are copied to the ID column. Empty
+ rownames are replaced with 1:nrow(fit).
+
+15 April 2013: limma 3.17.4
+
+- toptable() and topTable() now preserve rownames of fit object in
+ the rownames of the toptable.
+
+15 April 2013: limma 3.17.3
+
+- getEAWP() now forces the output exprs matrix to have unique row
+ names.
+
+- fitFDistRobustly() now uses a smoother for the smallest df.prior
+ values. This may result in smaller tail values than before when
+ a group of input x values appear to be outliers but the largest
+ value is not individually a stand-out value.
+
+10 April 2013: limma 3.17.2
-- bug fix and update to fitFDistRobustly().
+- improvements to help pages for data classes.
- 7 April 2013: limma 3.16.1
+ 7 April 2013: limma 3.17.1
- topTable() and treat() now give more informative error messages
when the argument fit is not a valid MArrayLM object.
+ 4 April 2013: limma 3.17.0 (Bioconductor 2.13 Developmental Branch)
4 April 2013: limma 3.16.0 (Bioconductor 2.12 Release Branch)
31 March 2013: limma 3.15.21
@@ -74,7 +214,7 @@
29 March 2013: limma 3.15.19
-- revise RNA-seq case study in User's Guide.
+- revise Nigerian RNA-seq case study in User's Guide.
- bug fix to lmscFit() when the residual df = 1.
diff --git a/inst/doc/limma.R b/inst/doc/limma.R
deleted file mode 100644
index ba4b3f9..0000000
--- a/inst/doc/limma.R
+++ /dev/null
@@ -1,2 +0,0 @@
-### R code from vignette source 'limma.Rnw'
-
diff --git a/inst/doc/limma.pdf b/inst/doc/limma.pdf
index a53a2ad..5d8cf4a 100644
Binary files a/inst/doc/limma.pdf and b/inst/doc/limma.pdf differ
diff --git a/man/01Introduction.Rd b/man/01Introduction.Rd
index c434bb6..0afd381 100644
--- a/man/01Introduction.Rd
+++ b/man/01Introduction.Rd
@@ -37,7 +37,8 @@ The function \code{\link{changeLog}} displays the record of changes to the packa
Smyth, G. K., Yang, Y.-H., Speed, T. P. (2003). Statistical issues in microarray data analysis. \emph{Methods in Molecular Biology} 224, 111-136.
Smyth, G. K. (2004). Linear models and empirical Bayes methods for assessing differential expression in microarray experiments.
-\emph{Statistical Applications in Genetics and Molecular Biology}, Volume \bold{3}, Article 3. \url{http://www.bepress.com/sagmb/vol3/iss1/art3}
+\emph{Statistical Applications in Genetics and Molecular Biology}, Volume \bold{3}, Article 3.
+\url{http://www.statsci.org/smyth/pubs/ebayes.pdf}
Smyth, G. K. (2005). Limma: linear models for microarray data.
In: \emph{Bioinformatics and Computational Biology Solutions using R and Bioconductor}.
diff --git a/man/06linearmodels.Rd b/man/06linearmodels.Rd
index 20ec6b9..1d46547 100644
--- a/man/06linearmodels.Rd
+++ b/man/06linearmodels.Rd
@@ -96,7 +96,8 @@ For multiple testing functions which operate on linear model fits, see \link{08.
\author{Gordon Smyth}
\references{
Smyth, G. K. (2004). Linear models and empirical Bayes methods for assessing differential expression in microarray experiments.
-\emph{Statistical Applications in Genetics and Molecular Biology}, \bold{3}, No. 1, Article 3. \url{http://www.bepress.com/sagmb/vol3/iss1/art3}
+\emph{Statistical Applications in Genetics and Molecular Biology}, \bold{3}, No. 1, Article 3.
+\url{http://www.statsci.org/smyth/pubs/ebayes.pdf}
Smyth, G. K., Michaud, J., and Scott, H. (2005). The use of within-array replicate spots for assessing differential expression in microarray experiments. Bioinformatics 21(9), 2067-2075.
}
diff --git a/man/asmalist.Rd b/man/asmalist.Rd
index 8897800..284621f 100755
--- a/man/asmalist.Rd
+++ b/man/asmalist.Rd
@@ -18,6 +18,13 @@ as.MAList(object)
Object of class \code{\link[=MAList-class]{MAList}}
}
+\details{
+The \code{marrayNorm} class is defined in the \code{marray} package.
+This function converts a normalized two color microarray data object created by the \code{marray} package into the corresponding limma data object.
+
+Note that such conversion is not necessary to access the limma linear modelling functions, because \code{lmFit} will operate on a \code{marrayNorm} data object directly.
+}
+
\author{Gordon Smyth}
\seealso{
diff --git a/man/beadCountWeights.Rd b/man/beadCountWeights.Rd
new file mode 100644
index 0000000..226a215
--- /dev/null
+++ b/man/beadCountWeights.Rd
@@ -0,0 +1,61 @@
+\name{beadCountWeights}
+\alias{beadCountWeights}
+
+\title{Bead Count Weights for Illumina BeadChips}
+
+\description{
+Estimates weights which account for biological variation and technical variation resulting from varying bead numbers.
+}
+
+\usage{
+beadCountWeights(y, x, design = NULL, bead.stdev = NULL, nbeads = NULL, array.cv = TRUE, scale = FALSE)
+}
+
+\arguments{
+ \item{y}{any matrix-like object containing log-expression values that can be coerced to a matrix.}
+ \item{x}{any matrix-like object containing raw expression values that can be coerced to a matrix.}
+ \item{design}{the design matrix of the microarray experiment, with rows
+ corresponding to arrays and columns to coefficients to be
+ estimated. Defaults to the unit vector meaning that the
+ arrays are treated as replicates.}
+ \item{bead.stdev}{numeric matrix containing bead-level standard deviations.}
+ \item{nbeads}{numeric matrix containing number of beads.}
+ \item{array.cv}{logical, should technical variation for each observation be calculated from a constant or array-specific coefficient of variation? The default is to use array-specific coefficients of variation.}
+ \item{scale}{logical, should weights be scaled so that the average weight size is the mean of the inverse technical variance along a probe? By default, weights are scaled so that the average weight size along a probe is 1.}
+}
+
+\details{
+The relative reliability of each gene on each array is measured by estimating its technical and biological variability.
+
+The technical variance for each gene on each array is inversely proportional to the number of beads and is estimated using array-specific bead-level coefficients of variation.
+
+Coefficients of variation are calculated using raw expression values.
+
+The biological variance for each gene across the arrays are estimated using Newton's iterations, with the assumption that the total residual deviance for each gene from \code{lmFit} is inversely proportional to the sum of the technical variance and biological variance.
+
+The bead number weights are an inverse of the sum of estimates for technical variances and biological variances.
+
+If any of the arguments \code{design}, \code{bead.stdev} or \code{nbeads} are set explicitly in the call they will over-ride the slots or components in the data \code{object}. The argument \code{design} does not normally need to be set in the call but will be extracted from the data \code{object} if available. If arguments \code{bead.stdev} and \code{nbeads} are not set explicitly in the call, it is necessary that they are available for extraction from the data \code{object}.
+}
+
+\value{
+A list object with the following components:
+
+\item{weights}{numeric matrix of bead number weights}
+\item{cv.constant}{numeric value of constant bead-level coefficient of variation}
+\item{cv.array}{numeric vector of array-specific bead-level coefficient of variation}
+\item{var.technical}{numeric matrix of technical variance}
+\item{var.biological}{numeric vector of biological variance}
+}
+
+\references{
+Law, CW (2013).
+\emph{Precision weights for gene expression analysis}.
+PhD Thesis. University of Melbourne, Australia.
+}
+
+\author{Charity Law and Gordon Smyth}
+
+\seealso{
+An overview of linear model functions in limma is given by \link{06.LinearModels}.
+}
diff --git a/man/contrastAsCoef.Rd b/man/contrastAsCoef.Rd
new file mode 100644
index 0000000..35f8841
--- /dev/null
+++ b/man/contrastAsCoef.Rd
@@ -0,0 +1,48 @@
+\name{contrastAsCoef}
+\alias{contrastAsCoef}
+\title{Reform a Design Matrix to that Contrasts Become Coefficients}
+\description{
+Reform a design matrix so that one or more coefficients from the new matrix correspond to specified contrasts of coefficients from the old matrix.
+}
+\usage{
+contrastAsCoef(design, contrast=NULL, first=TRUE)
+}
+\arguments{
+ \item{design}{numeric design matrix.}
+ \item{contrast}{numeric matrix with rows corresponding to columns of the design matrix (coefficients) and columns containing contrasts. May be a vector if there is only one contrast.}
+ \item{first}{logical, should coefficients corresponding to contrasts be the first columns (\code{TRUE}) or last columns (\code{FALSE}) of the output design matrix.}
+}
+
+\details{
+If \code{contrast} doesn't have full column rank, then superfluous columns are dropped.
+}
+
+\value{
+A list with components
+\item{design}{reformed design matrix}
+\item{coef}{columns of design matrix which hold the meaningful coefficients}
+\item{qr}{QR-decomposition of contrast matrix}
+}
+
+\seealso{
+\code{\link[stats]{model.matrix}} in the stats package.
+
+An overview of linear model functions in limma is given by \link{06.LinearModels}.
+}
+
+\author{Gordon Smyth}
+
+\examples{
+design <- cbind(1,c(0,0,1,1,0,0),c(0,0,0,0,1,1))
+cont <- c(0,-1,1)
+design2 <- contrastAsCoef(design, cont)$design
+
+# Original coef[3]-coef[2] becomes coef[1]
+y <- rnorm(6)
+fit1 <- lm(y~0+design)
+fit2 <- lm(y~0+design2)
+coef(fit1)
+coef(fit2)
+}
+
+\keyword{regression}
diff --git a/man/contrasts.fit.Rd b/man/contrasts.fit.Rd
index 84aef1b..f535e6c 100755
--- a/man/contrasts.fit.Rd
+++ b/man/contrasts.fit.Rd
@@ -9,7 +9,7 @@ contrasts.fit(fit, contrasts=NULL, coefficients=NULL)
}
\arguments{
\item{fit}{an \code{\link[limma:marraylm]{MArrayLM}} object or a list object produced by the function \code{lm.series} or equivalent. Must contain components \code{coefficients} and \code{stdev.unscaled}.}
- \item{contrasts}{numeric matrix with row corresponding to coefficients in \code{fit} and columns containing contrasts. May be a vector if there is only one contrast.}
+ \item{contrasts}{numeric matrix with rows corresponding to coefficients in \code{fit} and columns containing contrasts. May be a vector if there is only one contrast.}
\item{coefficients}{vector indicating which coefficients are to be kept in the revised fit object. An alternative way to specify the \code{contrasts}.}
}
\value{
diff --git a/man/ebayes.Rd b/man/ebayes.Rd
index 41f0e6e..4f340d7 100755
--- a/man/ebayes.Rd
+++ b/man/ebayes.Rd
@@ -64,6 +64,13 @@ Use \code{\link{topTreat}} to summarize output from \code{treat}.
Instead of testing for genes which have log-fold-changes different from zero, it tests whether the log2-fold-change is greater than \code{lfc} in absolute value (McCarthy and Smyth, 2009).
\code{treat} is concerned with p-values rather than posterior odds, so it does not compute the B-statistic \code{lods}.
The idea of thresholding doesn't apply to F-statistics in a straightforward way, so moderated F-statistics are also not computed.
+
+If \code{trend=TRUE} then an intensity-dependent trend is fitted to the prior variances \code{s2.prior}.
+Specifically, \code{squeezeVar} is called with the \code{covariate} equal to \code{Amean}, the average log2-intensity for each gene.
+See \code{\link{squeezeVar}} for more details.
+
+If \code{robust=TRUE} then the robust empirical Bayes procedure of Phipson et al (2013) is used.
+See \code{\link{squeezeVar}} for more details.
}
\seealso{
\code{\link{squeezeVar}}, \code{\link{fitFDist}}, \code{\link{tmixture.matrix}}.
@@ -77,8 +84,14 @@ McCarthy, D. J., and Smyth, G. K. (2009). Testing significance relative to a fol
Loennstedt, I., and Speed, T. P. (2002). Replicated microarray data. \emph{Statistica Sinica} \bold{12}, 31-46.
+Phipson, B, Lee, S, Majewski, IJ, Alexander, WS, and Smyth, GK (2013).
+Empirical Bayes in the presence of exceptional cases, with application to microarray data.
+Bioinformatics Division, Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia.
+\url{http://www.statsci.org/smyth/pubs/RobustEBayesPreprint.pdf}
+
Smyth, G. K. (2004). Linear models and empirical Bayes methods for assessing differential expression in microarray experiments.
-\emph{Statistical Applications in Genetics and Molecular Biology}, Volume \bold{3}, Article 3. \url{http://www.bepress.com/sagmb/vol3/iss1/art3}
+\emph{Statistical Applications in Genetics and Molecular Biology}, Volume \bold{3}, Article 3.
+\url{http://www.statsci.org/smyth/pubs/ebayes.pdf}
}
\examples{
# See also lmFit examples
diff --git a/man/fitfdist.Rd b/man/fitfdist.Rd
index dfe9e28..d7a14f5 100755
--- a/man/fitfdist.Rd
+++ b/man/fitfdist.Rd
@@ -37,11 +37,13 @@ A list containing the components
Smyth, G. K. (2004). Linear models and empirical Bayes methods for
assessing differential expression in microarray experiments.
\emph{Statistical Applications in Genetics and Molecular Biology}, \bold{3},
- No. 1, Article 3. \url{http://www.bepress.com/sagmb/vol3/iss1/art3}
+ No. 1, Article 3.
+ \url{http://www.statsci.org/smyth/pubs/ebayes.pdf}
Phipson, B. (2013).
\emph{Empirical Bayes modelling of expression profiles and their associations}.
PhD Thesis. University of Melbourne, Australia.
+\url{http://repository.unimelb.edu.au/10187/17614}
Phipson, B., and Smyth, G. K. (2013).
\emph{Robust empirical Bayes estimation protetcts against hyper-variable genes and improves power to detect differential expression in RNA-seq data}.
diff --git a/man/genas.Rd b/man/genas.Rd
index 8d5abc9..5ceed38 100644
--- a/man/genas.Rd
+++ b/man/genas.Rd
@@ -6,29 +6,47 @@
Calculates biological correlation between two gene expression profiles.
}
\usage{
-genas(fit, coef=c(1,2), chooseMethod=NULL,plot=FALSE,alpha=0.4)
+genas(fit, coef=c(1,2), subset="all", plot=FALSE, alpha=0.4)
}
+
\arguments{
- \item{fit}{an \code{MArrayLM} fitted model object produced by \code{lmFit} or \code{contrasts.fit} and followed by \code{eBayes}}
- \item{coef}{numeric vector of length 2 to indicate which contrasts/columns in the fit object are to be used}
- \item{chooseMethod}{character string, "n" for none, "Fpval" to subset using the F p-value, "p.union" to subset based on the union of the two contrasts' significant moderated t p-value, "p.int" to subset based on the intersection of the two contrasts' significant moderated t p-value, "logFC" to subset based on genes that have absolute logFC greater than the 90th quantile, "predFC" to subset based on genes that have absolute predictive logFC greater than the 90th quantile}
- \item{plot}{logical, if true a logFC versus logFC plot is outputted with biological and technical correlation represented by ellipses}
- \item{alpha}{plot option, a numeric value between 0 and 1 which determines the transparency of the ellipses}
- }
-\details{
-The biological correlation between the log fold changes of pairs of comparisons is computed. This method is to be applied when multiple groups (such as treatment groups, mutants or knock-outs) are compared back to the same control group.
+ \item{fit}{an \code{MArrayLM} fitted model object produced by \code{lmFit} or \code{contrasts.fit} and followed by \code{eBayes}.}
+ \item{coef}{numeric vector of length 2 indicating which columns in the fit object are to be correlated.}
+ \item{subset}{character string indicating which subset of genes to include in the correlation analysis.
+ Choices are \code{"all"}, \code{"Fpval"}, \code{"p.union"}, \code{"p.int"}, \code{"logFC"} or \code{"predFC"}.}
+ \item{plot}{logical, should a scatterplot be produced summarizing the correlation analysis?}
+ \item{alpha}{numeric value between 0 and 1 determining the transparency of the technical and biological ellipses if a plot is produced.
+ \code{alpha=0} indicates fully transparent and \code{alpha=1} indicates fully opague.}
+}
-This method is an extension of the empirical Bayes method of \code{limma}. It aims to separate the technical correlation, which comes from comparing multiple treatment/mutant/knock-out groups to the same control group, from biological correlation, which is the true correlation of the gene expression profiles between two treatment/mutant/knock-out groups.
+\details{
+The function estimates the biological correlation between two different contrasts in a linear model.
+By biological correlation, we mean the correlation that would exist between the log2-fold changes (logFC) for the two contrasts, if measurement error could be eliminated and the true log-fold-changes were known.
+This function is motivated by the fact that different contrasts for a linear model are often strongly correlated in a technical sense.
+For example, the estimated logFC for multiple treatment conditions compared back to the same control group will be positively correlated even in the absence of any biological effect.
+This function aims to separate the biological from the technical components of the correlation.
+The method is explained briefly in Majewski et al (2010) and in full detail in Phipson (2013).
-The \code{chooseMethod} argument specifies whether and how the fit object should be subsetted. The default is "n", which uses all genes in the fit object to estimate the biological correlation. Only genes that display evidence of differential expression can be used to estimate the biological correlation. The option "Fpval" chooses genes based on how many F p-values are estimated to be truly significant using the method \code{propNotDE}. This should capture genes that display any evidence [...]
+The \code{subset} argument specifies whether and how the fit object should be subsetted.
+Ideally, only genes that are truly differentially expressed for one or both of the contrasts should be used estimate the biological correlation.
+The default is \code{"all"}, which uses all genes in the fit object to estimate the biological correlation.
+The option \code{"Fpval"} chooses genes based on how many F-test p-values are estimated to be truly significant using the method \code{propNotDE}.
+This should capture genes that display any evidence of differential expression in either of the two contrasts.
+The options \code{"p.union"} and \code{"p.int"} are based on the moderated t p-values from both contrasts.
+From the \code{propNotDE} method an estimate of the number of p-values truly significant in either of the two contrasts can be obtained.
+"p.union" takes the union of these genes and \code{"p.int"} takes the intersection of these genes.
+The other options, \code{"logFC"} and \code{"predFC"} subsets on genes that attain a logFC or predFC at least as large as the 90th percentile of the log fold changes or predictive log fold changes on the absolute scale.
-The \code{plot} option is a logical argument that specifies whether or not to plot the log fold changes of the two contrasts. The biological and technical correlations are overlaid on the log fold change versus log fold change scatterplot using transparent ellipses. \code{library(ellipse)} is required to enable the plotting of ellipses. The \code{alpha} argument takes values between 0 and 1 and controls how transparent the ellipses are.
+The \code{plot} option is a logical argument that specifies whether or not to plot a scatter plot of log-fold-changes for the two contrasts.
+The biological and technical correlations are overlaid on the scatterplot using semi-transparent ellipses.
+\code{library(ellipse)} is required to enable the plotting of ellipses.
}
+
\value{
\code{genas} produces a list with the following components.
\item{technical.correlation}{estimate of the technical correlation}
- \item{covariance.matrix}{estimate of the covariance matrix from which the biological correlation is obtained}
\item{biological.correlation}{estimate of the biological correlation}
+ \item{covariance.matrix}{estimate of the covariance matrix from which the biological correlation is obtained}
\item{deviance}{the likelihood ratio test statistic used to test whether the biological correlation is equal to 0}
\item{p.value}{the p.value associated with \code{deviance}}
\item{n}{the number of genes used to estimate the biological correlation}
@@ -40,29 +58,43 @@ The \code{plot} option is a logical argument that specifies whether or not to pl
\author{Belinda Phipson and Gordon Smyth}
\references{
+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}
+
Phipson, B. (2013).
\emph{Empirical Bayes modelling of expression profiles and their associations}.
PhD Thesis. University of Melbourne, Australia.
+\url{http://repository.unimelb.edu.au/10187/17614}
}
\examples{
-library(limma)
-# Simulate gene expression data,
-# 6 microarrays with 1000 genes on each array
-set.seed(2004)
-y<-matrix(rnorm(6000),ncol=6)
+# Simulate gene expression data
+
+# Three conditions (Control, A and B) and 1000 genes
+ngene <- 1000
+mu.A <- mu.B <- mu.ctrl <- rep(5,ngene)
+
+# 200 genes are differentially expressed.
+# All are up in condition A and down in B
+# so the biological correlation is negative.
+mu.A[1:200] <- mu.ctrl[1:200]+2
+mu.B[1:200] <- mu.ctrl[1:200]-2
+
+# Two microarrays for each condition
+mu <- cbind(mu.ctrl,mu.ctrl,mu.A,mu.A,mu.B,mu.B)
+y <- matrix(rnorm(6000,mean=mu,sd=1),ngene,6)
# two experimental groups and one control group with two replicates each
-group<-factor(c("A","A","B","B","control","control"))
-design<-model.matrix(~0+group)
-colnames(design)<-c("A","B","control")
+group <- factor(c("Ctrl","Ctrl","A","A","B","B"), levels=c("Ctrl","A","B"))
+design <- model.matrix(~group)
# fit a linear model
-fit<-lmFit(y,design)
-contrasts<-makeContrasts(A-control,B-control,levels=design)
-fit2<-contrasts.fit(fit,contrasts)
-fit2<-eBayes(fit2)
+fit <- lmFit(y,design)
+fit <- eBayes(fit)
-# calculate biological correlation between the gene expression profiles of (A vs control) and (B vs control)
-genas(fit2)
-}
\ No newline at end of file
+# Estimate biological correlation between the logFC profiles
+# for A-vs-Ctrl and B-vs-Ctrl
+genas(fit, coef=c(2,3), plot=TRUE, subset="F")
+}
diff --git a/man/getEAWP.Rd b/man/getEAWP.Rd
index 0882b03..396c022 100644
--- a/man/getEAWP.Rd
+++ b/man/getEAWP.Rd
@@ -18,7 +18,11 @@ For single-channel objects, \code{Amean} is computed from the matrix of expressi
\code{PLMset}, \code{vsn} and \code{ExpressionSet} are assumed to be single-channel for this purpose.
If \code{object} is a matrix, it is assumed to contain log-intensities if the values are all positive and log-ratios otherwise.
-\code{Amean} is computed in the former case but not the latter.
+\code{Amean} is computed in the former case but not the latter.
+
+From April 2013, the output \code{exprs} matrix is ensured to have unique row names.
+If \code{object} has no row names, then the output row names of \code{exprs} are 1 to the number of rows.
+If \code{object} has row names but with duplicated names, then row names of \code{exprs} are set to 1 up to the number of rows and the original row names are preserved in the \code{ID} column of \code{probes}.
}
\value{
A list with components
@@ -26,6 +30,8 @@ A list with components
\item{weights}{numeric matrix of weights}
\item{probes}{data.frame of probe-annotation}
\item{Amean}{numeric vector of average log-expression for each probe}
+\code{exprs} is the only required component.
+The other components will be \code{NULL} if not found in the input object.
}
\author{Gordon Smyth}
\seealso{
diff --git a/man/lmFit.Rd b/man/lmFit.Rd
index b4269c8..8a181d2 100755
--- a/man/lmFit.Rd
+++ b/man/lmFit.Rd
@@ -19,6 +19,9 @@ lmFit(object,design=NULL,ndups=1,spacing=1,block=NULL,correlation,weights=NULL,m
\value{
An \code{\link[limma:marraylm]{MArrayLM}} object containing the result of the fits.
+
+The rownames of \code{object} are preserved in the fit object and can be retrieved by \code{rownames(fit)} where \code{fit} is output from \code{lmFit}.
+The column names of \code{design} are preserved as column names and can be retrieved by \code{colnames(fit)}.
}
\details{
diff --git a/man/loessfit.Rd b/man/loessfit.Rd
index ca12eff..4301a47 100755
--- a/man/loessfit.Rd
+++ b/man/loessfit.Rd
@@ -1,46 +1,52 @@
\name{loessFit}
\alias{loessFit}
-\title{Fast Simple Loess}
+\title{Univariate Lowess With Prior Weights}
+
\description{
-Locally weighted linear regression when there is only one x-variable and only the fitted values and residuals are required.
+Univariate locally weighted linear regression allowing for prior weights.
+Returns fitted values and residuals.
}
+
\usage{
-loessFit(y, x, weights=NULL, span=0.3, bin=NULL, iterations=4, min.weight=1e-5, max.weight=1e5, equal.weights.as.null=TRUE)
+loessFit(y, x, weights=NULL, span=0.3, iterations=4, min.weight=1e-5, max.weight=1e5, equal.weights.as.null=TRUE)
}
+
\arguments{
\item{y}{numeric vector of response values. Missing values are allowed.}
\item{x}{numeric vector of predictor values Missing values are allowed.}
\item{weights}{numeric vector of non-negative weights. Missing values are treated as zero.}
\item{span}{positive numeric value between 0 and 1 specifying proportion of data to be used in the local regression moving window.
Larger numbers give smoother fits.}
- \item{bin}{numeric value between 0 and 1 giving the proportion of the data that can be grouped in a single bin when doing local regression fit.
- \code{bin=0} forces an exact local regression fit with no interpolation. Default is 0.005 with weights and 0.01 without weights.}
- \item{iterations}{number of iterations of loess fit. Values greater than 1 produce robust fits.}
+ \item{iterations}{number of local regression fits. Values greater than 1 produce robust fits.}
\item{min.weight}{minimum weight. Any lower weights will be reset.}
- \item{max.weight}{maxium weight. Any higher weights will be reset.}
- \item{equal.weights.as.null}{should equal weights be treated as if weights were \code{NULL}, so that \code{lowess} is used instead of \code{loess}? Applies even all weights are all zero.}
+ \item{max.weight}{maximum weight. Any higher weights will be reset.}
+ \item{equal.weights.as.null}{should equal weights be treated as if weights were \code{NULL}, so that \code{lowess} is called? Applies even if all weights are all zero.}
}
\details{
-This function calls the \code{lowess} function when \code{weights=NULL} and \code{stats:::simpleLoess} otherwise.
-Its purpose is to give a unified and streamlined interface to \code{lowess} and \code{loess} for use in \code{\link{normalizeWithinArrays}}.
-When \code{weights} is null, this function call \code{lowess} in the stats package with appropropriate choice of tuning parameters.
-When \code{weights} is non-null, it is in effect a call to \code{loess} with \code{degree=1}.
-See the help pages for those functions for references and credits.
+This function is essentially a wrapper function for \code{lowess} and \code{locfit.raw} with added error checking.
+The idea is to provide the classic univariate lowess algoritm of Cleveland (1979) but allowing for prior weights and missing values.
+
+The venerable \code{lowess} code is fast, uses little memory and has an accurate interpolation scheme, so it is an advantage to use it when weights are not needed.
+This functions calls \code{lowess} when \code{weights=NULL}, but returns values in original rather than sorted order and allows missing values.
+The treatment of missing values is analogous to \code{na.exclude}.
+
+When \code{weights} are provided, this function makes repeated calls to \code{locfit.raw} with \code{deg=1}.
+Each iteration applies robustifying weights computed from the previous fit, as described by Cleveland (1981).
-Note that \code{lowess} is faster, needs less memory and has a more accurate interpolation scheme than \code{loess}, so it is desirable to use \code{lowess} whenever \code{loess} is not needed to handle quantitative weights.
+By default, \code{weights} that are all equal (even all zero) are treated as if they were \code{NULL}.
-The arguments \code{span}, \code{cell} and \code{iterations} here have the same meaning as in \code{loess}.
-\code{span} is equivalent to the argument \code{f} of \code{lowess} and \code{iterations} is equivalent to \code{iter+1}.
+In principle this function gives analogous results to \code{loess(y~x,weights=weights,span=span,degree=1,family="symmetric")}.
+However \code{lowess} and \code{locfit.raw} scale up much more efficiently for long data vectors.
-The parameter \code{bin} is intended to give a uniform interface to the \code{delta} argument of \code{lowess} and the \code{cell} argument of \code{loess}.
-\code{bin} translates to \code{delta=bin*diff(range(x))} in a call to \code{lowess} or to \code{cell=bin/span} in a call to \code{loess}.
-This is an attempt to put the \code{delta} and \code{cell} arguments on comparable scales.
+The arguments \code{span} and \code{iterations} here have the same meaning as for \code{loess}.
+\code{span} is equivalent to the argument \code{f} of \code{lowess} while \code{iterations} is equivalent to \code{iter+1} for \code{lowess}.
+It gives the total number of fits rather than the number of robustifying fits.
-Unlike \code{lowess}, \code{loessFit} returns values in original rather than sorted order.
-Also unlike \code{lowess}, \code{loessFit} allows missing values, the treatment being analogous to \code{na.exclude}.
-Unlike \code{loess}, \code{loessFit} returns a linear regression fit if there are insufficient observations to estimate the loess curve.
+When there are insufficient observations to estimate the loess curve, \code{loessFit} returns a linear regression fit.
+This mimics the behavior of \code{lowess} but not that of \code{loess} or \code{locfit.raw}.
}
+
\value{
A list with components
\item{fitted}{numeric vector of same length as \code{y} giving the loess fit}
@@ -48,24 +54,38 @@ A list with components
}
\author{Gordon Smyth}
+
+\references{
+Cleveland, W. S. (1979).
+Robust locally weighted regression and smoothing scatterplots.
+\emph{Journal of the American Statistical Association} 74, 829-836.
+}
\seealso{
-\code{\link[stats]{lowess}} and \code{\link[stats]{loess}} in the stats package.
+This function uses \code{\link{lowess}} and \code{\link[locfit]{locfit.raw}}.
+See the help pages of those functions for references and credits.
+
+Compare with \code{\link{loess}} in the stats package.
See \link{05.Normalization} for an outline of the limma package normalization functions.
}
\examples{
-y <- rnorm(1000)
-x <- rnorm(1000)
-w <- rep(1,1000)
-# The following are equivalent apart from execution time
-# and interpolation inaccuracies
-fit1 <- loessFit(y,x)$fitted
-fit2 <- loessFit(y,x,w)$fitted
-fit3 <- fitted(loess(y~x,weights=w,span=0.3,family="symmetric",iterations=4))
-# The same but with sorted x-values
-fit4 <- lowess(x,y,f=0.3)$y
+x <- (1:100)/101
+y <- sin(2*pi*x)+rnorm(100,sd=0.4)
+out <- loessFit(y,x)
+plot(x,y)
+lines(x,out$fitted,col="red")
+
+# Example using weights
+
+y <- x-0.5
+w <- rep(c(0,1),50)
+y[w==0] <- rnorm(50,sd=0.1)
+pch <- ifelse(w>0,16,1)
+plot(x,y,pch=pch)
+out <- loessFit(y,x,weights=w)
+lines(x,out$fitted,col="red")
}
\keyword{models}
diff --git a/man/plotma.Rd b/man/plotma.Rd
index e825a55..3ab86b2 100755
--- a/man/plotma.Rd
+++ b/man/plotma.Rd
@@ -1,14 +1,19 @@
\title{MA-Plot}
\name{plotMA}
\alias{plotMA}
+\alias{plotMA.RGList}
+\alias{plotMA.MAList}
+\alias{plotMA.EList}
+\alias{plotMA.MArrayLM}
+\alias{plotMA.default}
\description{
Creates an MA-plot with color coding for control spots.
}
\usage{
-plotMA(MA, array=1, xlab="A", ylab="M", main=colnames(MA)[array], xlim=NULL, ylim=NULL, status, values, pch, col, cex, legend=TRUE, zero.weights=FALSE, ...)
+\method{plotMA}{MAList}(MA, array=1, xlab="A", ylab="M", main=colnames(MA)[array], xlim=NULL, ylim=NULL, status, values, pch, col, cex, legend=TRUE, zero.weights=FALSE, ...)
}
\arguments{
- \item{MA}{an \code{RGList}, \code{MAList} or \code{MArrayLM} object, or any list with components \code{M} containing log-ratios and \code{A} containing average intensities.
+ \item{MA}{an \code{RGList}, \code{MAList}, \code{EList} or \code{MArrayLM} object.
Alternatively a \code{matrix} or \code{ExpressionSet} object.}
\item{array}{integer giving the array to be plotted. Corresponds to columns of \code{M} and \code{A}.}
\item{xlab}{character string giving label for x-axis}
diff --git a/man/predFCm.Rd b/man/predFCm.Rd
index 6f3e3bc..614b6e4 100644
--- a/man/predFCm.Rd
+++ b/man/predFCm.Rd
@@ -30,6 +30,7 @@ The predictive log fold changes are calculated as the posterior log fold changes
Phipson, B. (2013).
\emph{Empirical Bayes modelling of expression profiles and their associations}.
PhD Thesis. University of Melbourne, Australia.
+\url{http://repository.unimelb.edu.au/10187/17614}
}
\examples{
diff --git a/man/propTrueNull.Rd b/man/propTrueNull.Rd
index 2a88120..f02d7f5 100644
--- a/man/propTrueNull.Rd
+++ b/man/propTrueNull.Rd
@@ -60,6 +60,7 @@ Estimating the number of true null hypotheses from a histogram of p values.
Phipson, B (2013).
Empirical Bayes Modelling of Expression Profiles and Their Associations.
PhD Thesis, University of Melbourne, Australia.
+\url{http://repository.unimelb.edu.au/10187/17614}
}
\author{Belinda Phipson and Gordon Smyth for \code{propTrueNull}; Egil Ferkingstad, Mette Langaas and Marcus Davy for \code{convest}}
diff --git a/man/removeBatchEffect.Rd b/man/removeBatchEffect.Rd
index ea02c0f..aa8e25d 100644
--- a/man/removeBatchEffect.Rd
+++ b/man/removeBatchEffect.Rd
@@ -9,7 +9,7 @@ removeBatchEffect(x,batch=NULL,batch2=NULL,covariates=NULL,design=matrix(1,ncol(
}
\arguments{
\item{x}{numeric matrix, or any object that can be coerced to a matrix by \code{as.matrix(x)}, containing log-expression intensities for a series of microarrays.
- Rows correspond to probes and columns to arrays.}
+ Rows correspond to probes and columns to arrays. All values must be non-missing.}
\item{batch}{factor or vector indicating batches.}
\item{batch2}{factor or vector indicating batches.}
\item{covariates}{matrix or vector of covariates to be adjusted for.}
@@ -19,7 +19,7 @@ removeBatchEffect(x,batch=NULL,batch2=NULL,covariates=NULL,design=matrix(1,ncol(
A numeric matrix of log-expression values with batch and covariate effects removed.
}
\details{
-This function is useful for removing batch effects, associated with hybridization time or other technical variables, prior to clustering or unsupervised analysis such as PCA or heatmaps.
+This function is useful for removing batch effects, associated with hybridization time or other technical variables, prior to clustering or unsupervised analysis such as PCA, MDS or heatmaps.
It is not intended to use with linear modelling.
For linear modelling, it is better to include the batch factors in the linear model.
diff --git a/man/roast.Rd b/man/roast.Rd
index 1c705df..cae3ca0 100644
--- a/man/roast.Rd
+++ b/man/roast.Rd
@@ -129,6 +129,7 @@ Rotation tests.
Phipson B, and Smyth GK (2010).
Permutation P-values should never be zero: calculating exact P-values when permutations are randomly drawn.
\emph{Statistical Applications in Genetics and Molecular Biology}, Volume 9, Article 39.
+\url{http://www.statsci.org/smyth/pubs/PermPValuesPreprint.pdf}
Routledge, RD (1994).
Practicing safe statistics with the mid-p.
diff --git a/man/squeezeVar.Rd b/man/squeezeVar.Rd
index 263a20c..6a2c7b5 100755
--- a/man/squeezeVar.Rd
+++ b/man/squeezeVar.Rd
@@ -29,8 +29,13 @@ Specifically, the sample variances \code{var} are assumed to follow scaled chi-s
and an scaled inverse chi-squared prior is assumed for the true variances.
The scale and degrees of freedom of this prior distribution are estimated from the values of \code{var}.
-The effect of this function is to squeeze the variances towards a common value.
+The effect of this function is to squeeze the variances towards a common value, or to a global trend if a \code{covariate} is provided.
The squeezed variances have a smaller expected mean square error to the true variances than do the sample variances themselves.
+
+If \code{covariate} is non-null, then the scale parameter of the prior distribution is assumed to depend on the covariate.
+If the covariate is average log-expression, then the effect is an intensity-dependent trend similar to that in Sartor et al (2006).
+
+\code{robust=TRUE} implements the robust empirical Bayes procedure of Phipson et al (2013) which allows some of the \code{var} values to be outliers.
}
\note{
@@ -48,10 +53,20 @@ A list with components
\author{Gordon Smyth}
\references{
+Phipson, B, Lee, S, Majewski, IJ, Alexander, WS, and Smyth, GK (2013).
+Empirical Bayes in the presence of exceptional cases, with application to microarray data.
+Bioinformatics Division, Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia.
+\url{http://www.statsci.org/smyth/pubs/RobustEBayesPreprint.pdf}
+
+Sartor MA, Tomlinson CR, Wesselkamper SC, Sivaganesan S, Leikauf GD, Medvedovic M (2006).
+Intensity-based hierarchical Bayes method improves testing for differentially expressed genes in microarray experiments.
+BMC bioinformatics 7, 538.
+
Smyth, G. K. (2004). Linear models and empirical Bayes methods for
- assessing differential expression in microarray experiments.
- \emph{Statistical Applications in Genetics and Molecular Biology}, \bold{3},
- No. 1, Article 3. \url{http://www.bepress.com/sagmb/vol3/iss1/art3}
+assessing differential expression in microarray experiments.
+\emph{Statistical Applications in Genetics and Molecular Biology}, \bold{3},
+No. 1, Article 3.
+\url{http://www.statsci.org/smyth/pubs/ebayes.pdf}
}
\seealso{
diff --git a/man/targetsA2C.Rd b/man/targetsA2C.Rd
index 8dceafa..1b26f51 100755
--- a/man/targetsA2C.Rd
+++ b/man/targetsA2C.Rd
@@ -29,7 +29,7 @@ data.frame with twice as many rows as \code{targets}.
Any pair of columns named by \code{channel.columns} will now be one column.
}
\seealso{
-\code{targetsA2C} is used by the \code{\link[convert:coerce]{coerce}} method from \code{RGList} for \code{ExpressionSet} in the convert package.
+\code{targetsA2C} is used by the \code{coerce} method from \code{RGList} to \code{ExpressionSet} in the convert package.
An overview of methods for single channel analysis in limma is given by \link{07.SingleChannel}.
}
diff --git a/man/toptable.Rd b/man/toptable.Rd
index bc78d36..da794cb 100755
--- a/man/toptable.Rd
+++ b/man/toptable.Rd
@@ -30,8 +30,11 @@ topTreat(fit, coef=1, number=10, genelist=fit$genes, adjust.method="BH",
If \code{NULL}, this will be automatically generated.}
\item{adjust.method}{method used to adjust the p-values for multiple testing. Options, in increasing conservatism, include \code{"none"}, \code{"BH"}, \code{"BY"} and \code{"holm"}.
See \code{\link{p.adjust}} for the complete list of options. A \code{NULL} value will result in the default adjustment method, which is \code{"BH"}.}
- \item{sort.by}{character string specifying statistic to rank genes by. Possibilities for \code{topTable} and \code{toptable} are \code{"logFC"}, \code{"AveExpr"}, \code{"t"}, \code{"P"}, \code{"p"}, \code{"B"} or \code{"none"}. \code{"M"} is allowed as a synonym for \code{"logFC"} for backward compatibility. Other permitted synonyms are \code{"A"} or \code{"Amean"} for \code{"AveExpr"}, \code{"T"} for \code{"t"} and \code{"p"} for \code{"P"}. Possibilities for \code{topTableF} are \ [...]
- Possibilities for \code{topTreat} are as for \code{topTable} minus \code{"B"}.}
+ \item{sort.by}{character string specifying statistic to rank genes by.
+ Possible values for \code{topTable} and \code{toptable} are \code{"logFC"}, \code{"AveExpr"}, \code{"t"}, \code{"P"}, \code{"p"}, \code{"B"} or \code{"none"}.
+ (Permitted synonyms are \code{"M"} for \code{"logFC"}, \code{"A"} or \code{"Amean"} for \code{"AveExpr"}, \code{"T"} for \code{"t"} and \code{"p"} for \code{"P"}.)
+ Possibilities for \code{topTableF} are \code{"F"} or \code{"none"}.
+ Possibilities for \code{topTreat} are as for \code{topTable} except for \code{"B"}.}
\item{resort.by}{character string specifying statistic to sort the selected genes by in the output data.frame. Possibilities are the same as for \code{sort.by}.}
\item{p.value}{cutoff value for adjusted p-values. Only genes with lower p-values are listed.}
\item{lfc}{minimum absolute log2-fold-change required. \code{topTable} and \code{topTableF} include only genes with (at least one) absolute log-fold-changes greater than \code{lfc}. \code{topTreat} does not remove genes but ranks genes by evidence that their log-fold-change exceeds \code{lfc}.}
@@ -39,7 +42,7 @@ topTreat(fit, coef=1, number=10, genelist=fit$genes, adjust.method="BH",
\item{...}{any other arguments are passed to \code{ebayes} if \code{eb} is \code{NULL}}
}
\value{
- Produces a dataframe with a row for the \code{number} top genes and the following columns:
+ A dataframe with a row for the \code{number} top genes and the following columns:
\item{genelist}{one or more columns of probe annotation, if genelist was included as input}
\item{logFC}{estimate of the log2-fold-change corresponding to the effect or contrast (for \code{topTableF} there may be several columns of log-fold-changes)}
\item{CI.025}{left limit of confidence interval for \code{logFC} (if \code{confint=TRUE})}
@@ -50,11 +53,15 @@ topTreat(fit, coef=1, number=10, genelist=fit$genes, adjust.method="BH",
\item{P.Value}{raw p-value}
\item{adj.P.Value}{adjusted p-value or q-value}
\item{B}{log-odds that the gene is differentially expressed (omitted for \code{topTreat})}
+
+If \code{fit} had unique rownames, then the row.names of the above data.frame are the same in sorted order.
+Otherwise, the row.names of the data.frame indicate the row number in \code{fit}.
+If \code{fit} had duplicated row names, then these are preserved in the \code{ID} column of the data.frame, or in \code{ID0} if \code{genelist} already contained an \code{ID} column.
}
\details{
\code{toptable} is an earlier interface and is retained only for backward compatibility.
-This function summarizes a linear model fit object produced by \code{lmFit}, \code{lm.series}, \code{gls.series} or \code{mrlm} by selecting the top-ranked genes for any given contrast.
+These functions summarize the linear model fit object produced by \code{lmFit}, \code{lm.series}, \code{gls.series} or \code{mrlm} by selecting the top-ranked genes for any given contrast.
\code{topTable} and \code{topTableF} assume that the linear model fit has already been processed by \code{\link{eBayes}}.
\code{topTreat} assumes that the fit has been processed by \code{\link{treat}}.
diff --git a/src/normexp.c b/src/normexp.c
index 54ce074..faccacf 100644
--- a/src/normexp.c
+++ b/src/normexp.c
@@ -46,7 +46,6 @@ double normexp_m2loglik_saddle(int m, double *par, void *ex){
double thetaQuadratic;
int keepRepeating = 1;
int j,i;
- double maxDeviation;
int *hasConverged; // vector of 0/1 indicators,
// hasConverged[i] = 0 if theta[i] has not yet converged,
// hasConverged[i] = 1 if theta[i] has converged,
@@ -84,7 +83,6 @@ double normexp_m2loglik_saddle(int m, double *par, void *ex){
j = 0;
while(keepRepeating == 1){
j++;
- maxDeviation = 0.0;
for(i = 0; i < *n; i++){
// Only loop over entries of theta[.] that haven't yet converged
if(hasConverged[i] == 0){
@@ -150,7 +148,7 @@ void fit_saddle_nelder_mead(double *par, double *X, int *N, int *fail, int *fnco
parsOut[0] = par[0];
parsOut[1] = par[1];
parsOut[2] = par[2];
- double abstol = -1e500; // infinity
+ double abstol = -1e308; // infinity
double intol = 1.490116e-08; // square root of machine precision
void ex();
double alpha = 1.0;
diff --git a/tests/limma-Tests.R b/tests/limma-Tests.R
index 43363c4..b188f63 100755
--- a/tests/limma-Tests.R
+++ b/tests/limma-Tests.R
@@ -48,7 +48,8 @@ w[1:50] <- 0.5
out <- loessFit(y,x,weights=w)
f2 <- quantile(out$fitted)
r2 <- quantile(out$residual)
-w[1:80] <- 0
+w <- rep(1,100)
+w[2*(1:50)] <- 0
out <- loessFit(y,x,weights=w)
f3 <- quantile(out$fitted)
r3 <- quantile(out$residual)
@@ -117,6 +118,7 @@ toptable(fit)
### topTable
+rownames(M) <- LETTERS[1:10]
fit <- lmFit(M,design)
fit2 <- eBayes(contrasts.fit(fit,contrasts=contrast.matrix))
topTable(fit2)
@@ -248,3 +250,12 @@ fit$df.residual[1]
fit$df.prior
fit$s2.prior
summary(fit$s2.post)
+
+### voom
+
+y <- matrix(rpois(100*4,lambda=20),100,4)
+design <- cbind(Int=1,x=c(0,0,1,1))
+v <- voom(y,design)
+names(v)
+summary(v$E)
+summary(v$weights)
diff --git a/tests/limma-Tests.Rout.save b/tests/limma-Tests.Rout.save
index 896c075..baffee5 100755
--- a/tests/limma-Tests.Rout.save
+++ b/tests/limma-Tests.Rout.save
@@ -1,7 +1,6 @@
-R version 2.15.2 (2012-10-26) -- "Trick or Treat"
-Copyright (C) 2012 The R Foundation for Statistical Computing
-ISBN 3-900051-07-0
+R version 3.0.2 (2013-09-25) -- "Frisbee Sailing"
+Copyright (C) 2013 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
@@ -595,17 +594,18 @@ $G
> out <- loessFit(y,x,weights=w)
> f2 <- quantile(out$fitted)
> r2 <- quantile(out$residual)
-> w[1:80] <- 0
+> w <- rep(1,100)
+> w[2*(1:50)] <- 0
> out <- loessFit(y,x,weights=w)
> f3 <- quantile(out$fitted)
> r3 <- quantile(out$residual)
> data.frame(f1,f2,f3,r1,r2,r3)
f1 f2 f3 r1 r2 r3
-0% -0.78835384 -0.78097895 -0.78367166 -2.04434053 -2.04315267 -11.5610680
-25% -0.18340154 -0.18907787 -0.15525721 -0.59321065 -0.59309327 -0.8249477
-50% -0.11492924 -0.12136183 -0.03316003 0.05874864 0.08898459 -0.2466309
-75% 0.01507921 -0.01000344 0.13229151 0.56010750 0.56606786 0.4502908
-100% 0.21653837 0.21604173 11.69912073 2.57936026 2.56259812 2.5149556
+0% -0.78835384 -0.78957137 -0.64724211 -2.04434053 -2.02404318 -2.17625016
+25% -0.18340154 -0.18979269 -0.38092221 -0.59321065 -0.58975649 -0.68935397
+50% -0.11492924 -0.12087983 -0.16793705 0.05874864 0.08335198 0.07929677
+75% 0.01507921 -0.01857508 0.05351707 0.56010750 0.57618740 0.64831434
+100% 0.21653837 0.19214597 0.54678196 2.57936026 2.57549257 2.30629358
>
> ### normalizeWithinArrays
>
@@ -791,80 +791,92 @@ Array 2 corrected
>
> ### topTable
>
+> rownames(M) <- LETTERS[1:10]
> fit <- lmFit(M,design)
> fit2 <- eBayes(contrasts.fit(fit,contrasts=contrast.matrix))
> topTable(fit2)
- First3 Last3 Last3.First3 F P.Value adj.P.Val
-1 1.77602021 0.06025114 -1.71576906 42.67587166 2.924856e-19 2.924856e-18
-4 -0.05454069 0.39127869 0.44581938 2.10919528 1.213356e-01 5.350369e-01
-6 -0.16249607 -0.33009728 -0.16760121 1.82939237 1.605111e-01 5.350369e-01
-7 0.30852468 -0.06873462 -0.37725930 1.35021957 2.591833e-01 6.479584e-01
-8 -0.16942269 0.20578118 0.37520387 0.96017584 3.828256e-01 7.656511e-01
-10 0.21417623 0.07074940 -0.14342683 0.68755586 5.028035e-01 8.039105e-01
-3 -0.12236781 0.15095948 0.27332729 0.51032807 6.002986e-01 8.039105e-01
-2 -0.11982833 0.13529287 0.25512120 0.44141085 6.431284e-01 8.039105e-01
-5 0.01897934 0.10434934 0.08536999 0.15202008 8.589710e-01 9.496107e-01
-9 -0.04720963 0.03996397 0.08717360 0.05170315 9.496107e-01 9.496107e-01
+ First3 Last3 Last3.First3 AveExpr F P.Value
+A 1.77602021 0.06025114 -1.71576906 0.918135675 42.67587166 2.924856e-19
+D -0.05454069 0.39127869 0.44581938 0.168369004 2.10919528 1.213356e-01
+F -0.16249607 -0.33009728 -0.16760121 -0.246296671 1.82939237 1.605111e-01
+G 0.30852468 -0.06873462 -0.37725930 0.119895035 1.35021957 2.591833e-01
+H -0.16942269 0.20578118 0.37520387 0.018179245 0.96017584 3.828256e-01
+J 0.21417623 0.07074940 -0.14342683 0.142462814 0.68755586 5.028035e-01
+C -0.12236781 0.15095948 0.27332729 0.014295836 0.51032807 6.002986e-01
+B -0.11982833 0.13529287 0.25512120 0.007732271 0.44141085 6.431284e-01
+E 0.01897934 0.10434934 0.08536999 0.061664340 0.15202008 8.589710e-01
+I -0.04720963 0.03996397 0.08717360 -0.003622829 0.05170315 9.496107e-01
+ adj.P.Val
+A 2.924856e-18
+D 5.350369e-01
+F 5.350369e-01
+G 6.479584e-01
+H 7.656511e-01
+J 8.039105e-01
+C 8.039105e-01
+B 8.039105e-01
+E 9.496107e-01
+I 9.496107e-01
> topTable(fit2,coef=3,resort.by="logFC")
- logFC t P.Value adj.P.Val B
-4 0.44581938 1.6389001 1.090771e-01 4.386591e-01 -5.494162
-8 0.37520387 1.3793067 1.754637e-01 4.386591e-01 -5.881353
-3 0.27332729 1.0047928 3.210366e-01 5.899046e-01 -6.322593
-2 0.25512120 0.9378645 3.539428e-01 5.899046e-01 -6.386845
-9 0.08717360 0.3204634 7.502852e-01 7.552787e-01 -6.770768
-5 0.08536999 0.3138331 7.552787e-01 7.552787e-01 -6.772846
-10 -0.14342683 -0.5272589 6.009262e-01 7.511577e-01 -6.684136
-6 -0.16760121 -0.6161276 5.413009e-01 7.511577e-01 -6.633922
-7 -0.37725930 -1.3868627 1.731642e-01 4.386591e-01 -5.871024
-1 -1.71576906 -6.3074289 1.747263e-07 1.747263e-06 12.838652
+ logFC AveExpr t P.Value adj.P.Val B
+D 0.44581938 0.168369004 1.6389001 1.090771e-01 4.386591e-01 -5.494162
+H 0.37520387 0.018179245 1.3793067 1.754637e-01 4.386591e-01 -5.881353
+C 0.27332729 0.014295836 1.0047928 3.210366e-01 5.899046e-01 -6.322593
+B 0.25512120 0.007732271 0.9378645 3.539428e-01 5.899046e-01 -6.386845
+I 0.08717360 -0.003622829 0.3204634 7.502852e-01 7.552787e-01 -6.770768
+E 0.08536999 0.061664340 0.3138331 7.552787e-01 7.552787e-01 -6.772846
+J -0.14342683 0.142462814 -0.5272589 6.009262e-01 7.511577e-01 -6.684136
+F -0.16760121 -0.246296671 -0.6161276 5.413009e-01 7.511577e-01 -6.633922
+G -0.37725930 0.119895035 -1.3868627 1.731642e-01 4.386591e-01 -5.871024
+A -1.71576906 0.918135675 -6.3074289 1.747263e-07 1.747263e-06 12.838652
> topTable(fit2,coef=3,resort.by="p")
- logFC t P.Value adj.P.Val B
-1 -1.71576906 -6.3074289 1.747263e-07 1.747263e-06 12.838652
-4 0.44581938 1.6389001 1.090771e-01 4.386591e-01 -5.494162
-7 -0.37725930 -1.3868627 1.731642e-01 4.386591e-01 -5.871024
-8 0.37520387 1.3793067 1.754637e-01 4.386591e-01 -5.881353
-3 0.27332729 1.0047928 3.210366e-01 5.899046e-01 -6.322593
-2 0.25512120 0.9378645 3.539428e-01 5.899046e-01 -6.386845
-6 -0.16760121 -0.6161276 5.413009e-01 7.511577e-01 -6.633922
-10 -0.14342683 -0.5272589 6.009262e-01 7.511577e-01 -6.684136
-9 0.08717360 0.3204634 7.502852e-01 7.552787e-01 -6.770768
-5 0.08536999 0.3138331 7.552787e-01 7.552787e-01 -6.772846
+ logFC AveExpr t P.Value adj.P.Val B
+A -1.71576906 0.918135675 -6.3074289 1.747263e-07 1.747263e-06 12.838652
+D 0.44581938 0.168369004 1.6389001 1.090771e-01 4.386591e-01 -5.494162
+G -0.37725930 0.119895035 -1.3868627 1.731642e-01 4.386591e-01 -5.871024
+H 0.37520387 0.018179245 1.3793067 1.754637e-01 4.386591e-01 -5.881353
+C 0.27332729 0.014295836 1.0047928 3.210366e-01 5.899046e-01 -6.322593
+B 0.25512120 0.007732271 0.9378645 3.539428e-01 5.899046e-01 -6.386845
+F -0.16760121 -0.246296671 -0.6161276 5.413009e-01 7.511577e-01 -6.633922
+J -0.14342683 0.142462814 -0.5272589 6.009262e-01 7.511577e-01 -6.684136
+I 0.08717360 -0.003622829 0.3204634 7.502852e-01 7.552787e-01 -6.770768
+E 0.08536999 0.061664340 0.3138331 7.552787e-01 7.552787e-01 -6.772846
> topTable(fit2,coef=3,sort="logFC",resort.by="t")
- logFC t P.Value adj.P.Val B
-4 0.44581938 1.6389001 1.090771e-01 4.386591e-01 -5.494162
-8 0.37520387 1.3793067 1.754637e-01 4.386591e-01 -5.881353
-3 0.27332729 1.0047928 3.210366e-01 5.899046e-01 -6.322593
-2 0.25512120 0.9378645 3.539428e-01 5.899046e-01 -6.386845
-9 0.08717360 0.3204634 7.502852e-01 7.552787e-01 -6.770768
-5 0.08536999 0.3138331 7.552787e-01 7.552787e-01 -6.772846
-10 -0.14342683 -0.5272589 6.009262e-01 7.511577e-01 -6.684136
-6 -0.16760121 -0.6161276 5.413009e-01 7.511577e-01 -6.633922
-7 -0.37725930 -1.3868627 1.731642e-01 4.386591e-01 -5.871024
-1 -1.71576906 -6.3074289 1.747263e-07 1.747263e-06 12.838652
+ logFC AveExpr t P.Value adj.P.Val B
+D 0.44581938 0.168369004 1.6389001 1.090771e-01 4.386591e-01 -5.494162
+H 0.37520387 0.018179245 1.3793067 1.754637e-01 4.386591e-01 -5.881353
+C 0.27332729 0.014295836 1.0047928 3.210366e-01 5.899046e-01 -6.322593
+B 0.25512120 0.007732271 0.9378645 3.539428e-01 5.899046e-01 -6.386845
+I 0.08717360 -0.003622829 0.3204634 7.502852e-01 7.552787e-01 -6.770768
+E 0.08536999 0.061664340 0.3138331 7.552787e-01 7.552787e-01 -6.772846
+J -0.14342683 0.142462814 -0.5272589 6.009262e-01 7.511577e-01 -6.684136
+F -0.16760121 -0.246296671 -0.6161276 5.413009e-01 7.511577e-01 -6.633922
+G -0.37725930 0.119895035 -1.3868627 1.731642e-01 4.386591e-01 -5.871024
+A -1.71576906 0.918135675 -6.3074289 1.747263e-07 1.747263e-06 12.838652
> topTable(fit2,coef=3,resort.by="B")
- logFC t P.Value adj.P.Val B
-1 -1.71576906 -6.3074289 1.747263e-07 1.747263e-06 12.838652
-4 0.44581938 1.6389001 1.090771e-01 4.386591e-01 -5.494162
-7 -0.37725930 -1.3868627 1.731642e-01 4.386591e-01 -5.871024
-8 0.37520387 1.3793067 1.754637e-01 4.386591e-01 -5.881353
-3 0.27332729 1.0047928 3.210366e-01 5.899046e-01 -6.322593
-2 0.25512120 0.9378645 3.539428e-01 5.899046e-01 -6.386845
-6 -0.16760121 -0.6161276 5.413009e-01 7.511577e-01 -6.633922
-10 -0.14342683 -0.5272589 6.009262e-01 7.511577e-01 -6.684136
-9 0.08717360 0.3204634 7.502852e-01 7.552787e-01 -6.770768
-5 0.08536999 0.3138331 7.552787e-01 7.552787e-01 -6.772846
+ logFC AveExpr t P.Value adj.P.Val B
+A -1.71576906 0.918135675 -6.3074289 1.747263e-07 1.747263e-06 12.838652
+D 0.44581938 0.168369004 1.6389001 1.090771e-01 4.386591e-01 -5.494162
+G -0.37725930 0.119895035 -1.3868627 1.731642e-01 4.386591e-01 -5.871024
+H 0.37520387 0.018179245 1.3793067 1.754637e-01 4.386591e-01 -5.881353
+C 0.27332729 0.014295836 1.0047928 3.210366e-01 5.899046e-01 -6.322593
+B 0.25512120 0.007732271 0.9378645 3.539428e-01 5.899046e-01 -6.386845
+F -0.16760121 -0.246296671 -0.6161276 5.413009e-01 7.511577e-01 -6.633922
+J -0.14342683 0.142462814 -0.5272589 6.009262e-01 7.511577e-01 -6.684136
+I 0.08717360 -0.003622829 0.3204634 7.502852e-01 7.552787e-01 -6.770768
+E 0.08536999 0.061664340 0.3138331 7.552787e-01 7.552787e-01 -6.772846
> topTable(fit2,coef=3,lfc=1)
- logFC t P.Value adj.P.Val B
-1 -1.715769 -6.307429 1.747263e-07 1.747263e-06 12.83865
+ logFC AveExpr t P.Value adj.P.Val B
+A -1.715769 0.9181357 -6.307429 1.747263e-07 1.747263e-06 12.83865
> topTable(fit2,coef=3,p=0.2)
- logFC t P.Value adj.P.Val B
-1 -1.715769 -6.307429 1.747263e-07 1.747263e-06 12.83865
+ logFC AveExpr t P.Value adj.P.Val B
+A -1.715769 0.9181357 -6.307429 1.747263e-07 1.747263e-06 12.83865
> topTable(fit2,coef=3,p=0.2,lfc=0.5)
- logFC t P.Value adj.P.Val B
-1 -1.715769 -6.307429 1.747263e-07 1.747263e-06 12.83865
+ logFC AveExpr t P.Value adj.P.Val B
+A -1.715769 0.9181357 -6.307429 1.747263e-07 1.747263e-06 12.83865
> topTable(fit2,coef=3,p=0.2,lfc=0.5,sort="none")
- logFC t P.Value adj.P.Val B
-1 -1.715769 -6.307429 1.747263e-07 1.747263e-06 12.83865
+ logFC AveExpr t P.Value adj.P.Val B
+A -1.715769 0.9181357 -6.307429 1.747263e-07 1.747263e-06 12.83865
>
> designlist <- list(Null=matrix(1,6,1),Two=design,Three=cbind(1,c(0,0,1,1,0,0),c(0,0,0,0,1,1)))
> out <- selectModel(M,designlist)
@@ -923,29 +935,29 @@ Warning message:
In rlm.default(x = X, y = y, weights = w, ...) :
'rlm' failed to converge in 20 steps
> fit$coef
- First3Arrays Last3Arrays
- [1,] 1.75138894 0.06025114
- [2,] -0.11982833 0.10322039
- [3,] -0.09302502 0.15095948
- [4,] -0.05454069 0.33700045
- [5,] 0.07927938 0.10434934
- [6,] -0.16249607 -0.34010852
- [7,] 0.30852468 -0.06873462
- [8,] -0.16942269 0.24392984
- [9,] -0.04720963 0.03996397
-[10,] 0.21417623 -0.05679272
+ First3Arrays Last3Arrays
+A 1.75138894 0.06025114
+B -0.11982833 0.10322039
+C -0.09302502 0.15095948
+D -0.05454069 0.33700045
+E 0.07927938 0.10434934
+F -0.16249607 -0.34010852
+G 0.30852468 -0.06873462
+H -0.16942269 0.24392984
+I -0.04720963 0.03996397
+J 0.21417623 -0.05679272
> fit$stdev.unscaled
- First3Arrays Last3Arrays
- [1,] 0.5933418 0.5773503
- [2,] 0.5773503 0.6096497
- [3,] 0.6017444 0.5773503
- [4,] 0.5773503 0.6266021
- [5,] 0.6307703 0.5773503
- [6,] 0.5773503 0.5846707
- [7,] 0.5773503 0.5773503
- [8,] 0.5773503 0.6544564
- [9,] 0.5773503 0.5773503
-[10,] 0.5773503 0.6689776
+ First3Arrays Last3Arrays
+A 0.5933418 0.5773503
+B 0.5773503 0.6096497
+C 0.6017444 0.5773503
+D 0.5773503 0.6266021
+E 0.6307703 0.5773503
+F 0.5773503 0.5846707
+G 0.5773503 0.5773503
+H 0.5773503 0.6544564
+I 0.5773503 0.5773503
+J 0.5773503 0.6689776
> fit$sigma
[1] 0.2894294 0.2679396 0.2090236 0.1461395 0.2309018 0.2827476 0.2285945
[8] 0.2267556 0.3537469 0.2172409
@@ -999,6 +1011,9 @@ alpha-beta 0.6666667 3.333333e-01 -3.333333e-01
mu+alpha 0.3333333 3.333333e-01 9.280771e-17
mu+beta -0.3333333 9.280771e-17 3.333333e-01
+$Amean
+[1] 0.2034961 0.1954604 -0.2863347 0.1188659 0.1784593
+
$method
[1] "ls"
@@ -1247,6 +1262,30 @@ Loading required package: splines
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.2495 0.2749 0.3091 0.3512 0.3649 0.9233
>
+> ### voom
+>
+> y <- matrix(rpois(100*4,lambda=20),100,4)
+> design <- cbind(Int=1,x=c(0,0,1,1))
+> v <- voom(y,design)
+> names(v)
+[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
+> 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
+>
> proc.time()
user system elapsed
- 1.02 0.07 1.09
+ 1.34 0.10 1.45
--
Alioth's /git/debian-med/git-commit-notice on /srv/git.debian.org/git/debian-med/r-bioc-limma.git
More information about the debian-med-commit
mailing list