[med-svn] [r-bioc-limma] 02/04: Imported Upstream version 3.20.7+dfsg
Andreas Tille
tille at debian.org
Fri Jun 27 17:52:13 UTC 2014
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-bioc-limma.
commit fc787934936850e8f501f2c5daf105a178efc191
Author: Andreas Tille <tille at debian.org>
Date: Fri Jun 27 19:01:14 2014 +0200
Imported Upstream version 3.20.7+dfsg
---
DESCRIPTION | 25 ++--
NAMESPACE | 12 +-
R/barcodeplot.R | 183 ++++++++++++++++++------
R/beadCountWeights.R | 28 ++--
R/diffSplice.R | 210 ++++++++++++++++++++++++++++
R/dups.R | 7 +-
R/fitFDistRobustly.R | 13 +-
R/geneset-camera.R | 53 ++++---
R/geneset-roast.R | 325 +++++++++++++++++++++++--------------------
R/geneset-romer.R | 24 ++--
R/lmfit.R | 6 +-
R/loessFit.R | 106 ++++++++++----
R/plot.R | 70 ++++++++++
R/plotMDS.R | 53 ++++---
R/plots-ma.R | 28 +++-
R/read-ilmn.R | 10 +-
R/read.idat.R | 41 ++++++
R/removeBatchEffect.R | 16 +--
R/subsetting.R | 6 +-
R/toptable.R | 31 +++--
R/treat.R | 16 ++-
R/tricubeMovingAverage.R | 18 +++
R/utility.R | 6 +-
R/weightedLowess.R | 49 +++++++
R/zscore.R | 25 +++-
build/vignette.rds | Bin 229 -> 230 bytes
inst/NEWS.Rd | 117 ++++++++++++++++
inst/doc/changelog.txt | 251 ++++++++++++++++++++++++++++++---
inst/doc/intro.pdf | Bin 46191 -> 46191 bytes
man/01Introduction.Rd | 30 +++-
man/02classes.Rd | 2 +-
man/03reading.Rd | 4 +-
man/04Background.Rd | 2 +-
man/05Normalization.Rd | 2 +-
man/06linearmodels.Rd | 2 +-
man/07SingleChannel.Rd | 2 +-
man/08Tests.Rd | 2 +-
man/09Diagnostics.Rd | 4 +-
man/10GeneSetTests.Rd | 36 +++++
man/10Other.Rd | 10 --
man/11RNAseq.Rd | 31 +++++
man/EList.Rd | 25 ++--
man/PrintLayout.Rd | 9 +-
man/alias2Symbol.Rd | 7 +-
man/arrayWeights.Rd | 11 +-
man/arrayWeightsQuick.Rd | 6 +-
man/asMatrixWeights.Rd | 1 -
man/avearrays.Rd | 3 +-
man/barcodeplot.Rd | 17 ++-
man/beadCountWeights.Rd | 8 +-
man/blockDiag.Rd | 12 +-
man/camera.Rd | 22 +--
man/changelog.Rd | 4 +
man/classifytests.Rd | 2 +
man/controlStatus.Rd | 12 +-
man/diffSplice.Rd | 41 ++++++
man/ebayes.Rd | 6 +-
man/genas.Rd | 2 +-
man/geneSetTest.Rd | 7 +-
man/heatdiagram.Rd | 8 +-
man/imageplot3by2.Rd | 3 +-
man/kooperberg.Rd | 6 +-
man/lmFit.Rd | 14 +-
man/loessfit.Rd | 37 +++--
man/nec.Rd | 10 +-
man/normalizeWithinArrays.Rd | 4 +-
man/normalizeprintorder.Rd | 9 +-
man/{plotma.Rd => plot.Rd} | 55 +++++---
man/plotMDS.Rd | 17 ++-
man/plotSA.Rd | 3 +-
man/plotSplice.Rd | 29 ++++
man/plotma.Rd | 13 +-
man/plotma3by2.Rd | 5 +-
man/poolvar.Rd | 6 -
man/printtipWeights.Rd | 34 +++--
man/qqt.Rd | 6 +-
man/read.columns.Rd | 4 +-
man/read.idat.Rd | 70 ++++++++++
man/read.ilmn.Rd | 27 ++--
man/read.ilmn.targets.Rd | 2 +-
man/read.maimages.Rd | 5 +-
man/readGPRHeader.Rd | 12 +-
man/removeBatchEffect.Rd | 19 ++-
man/roast.Rd | 24 ++--
man/romer.Rd | 8 +-
man/symbols2indices.Rd | 2 +
man/targetsA2C.Rd | 5 +-
man/topRomer.Rd | 7 +-
man/topSplice.Rd | 34 +++++
man/toptable.Rd | 13 +-
man/tricubeMovingAverage.Rd | 36 +++++
man/volcanoplot.Rd | 4 +-
man/voom.Rd | 7 +-
man/vooma.Rd | 3 +-
man/weightedLowess.Rd | 72 ++++++++++
man/writefit.Rd | 3 +-
man/zscore.Rd | 40 ++++--
src/weighted_lowess.c | 271 ++++++++++++++++++++++++++++++++++++
tests/limma-Tests.R | 17 ++-
tests/limma-Tests.Rout.save | 139 ++++++++++--------
100 files changed, 2434 insertions(+), 700 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index acc3ea4..557914d 100755
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,17 +1,22 @@
Package: limma
-Version: 3.18.10
-Date: 2014/01/31
+Version: 3.20.7
+Date: 2014/06/24
Title: Linear Models for Microarray Data
-Author: Gordon Smyth [cre,aut], Matthew Ritchie [ctb], Jeremy Silver [ctb], James Wettenhall [ctb], Natalie Thorne [ctb], Mette Langaas [ctb], Egil Ferkingstad [ctb], Marcus Davy [ctb], Francois Pepin [ctb], Dongseok Choi [ctb], Davis McCarthy [ctb], Di Wu [ctb], Alicia Oshlack [ctb], Carolyn de Graaf [ctb], Yifang Hu [ctb], Wei Shi [ctb], Belinda Phipson [ctb]
+Description: Data analysis, linear models and differential expression for microarray data.
+Author: Gordon Smyth [cre,aut], Matthew Ritchie [ctb], Jeremy Silver [ctb], James Wettenhall [ctb], Natalie Thorne [ctb], Davis McCarthy [ctb], Di Wu [ctb], Yifang Hu [ctb], Wei Shi [ctb], Belinda Phipson [ctb], Alicia Oshlack [ctb], Carolyn de Graaf [ctb], Mette Langaas [ctb], Egil Ferkingstad [ctb], Marcus Davy [ctb], Francois Pepin [ctb], Dongseok Choi [ctb]
Maintainer: Gordon Smyth <smyth at wehi.edu.au>
+License: GPL (>=2)
Depends: R (>= 2.3.0), methods
Suggests: statmod (>= 1.2.2), splines, locfit, MASS, ellipse, affy,
- vsn, AnnotationDbi, org.Hs.eg.db
+ vsn, AnnotationDbi, org.Hs.eg.db, illuminaio
LazyLoad: yes
-Description: Data analysis, linear models and differential expression for microarray data.
-License: GPL (>=2)
URL: http://bioinf.wehi.edu.au/limma
-biocViews: Microarray, OneChannel, TwoChannel, DataImport,
- QualityControl, Preprocessing, Bioinformatics,
- DifferentialExpression, MultipleComparisons, TimeCourse
-Packaged: 2014-02-01 04:15:50 UTC; biocbuild
+biocViews: ExonArray, GeneExpression, Transcription,
+ AlternativeSplicing, DifferentialExpression,
+ DifferentialSplicing, GeneSetEnrichment, DataImport, Genetics,
+ Bayesian, Clustering, Regression, TimeCourse, Microarray,
+ microRNAArray, mRNAMicroarray, OneChannel,
+ ProprietaryPlatforms, TwoChannel, RNASeq, BatchEffect,
+ MultipleComparison, Normalization, Preprocessing,
+ QualityControl
+Packaged: 2014-06-26 02:22:20 UTC; biocbuild
diff --git a/NAMESPACE b/NAMESPACE
index 475997c..357432b 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -4,6 +4,8 @@ exportPattern("^[^\\.]")
exportClasses("RGList","MAList","EListRaw","EList","MArrayLM","TestResults","LargeDataObject","Roast","MDS")
exportMethods("show")
+import(methods)
+
S3method("[",MAList)
S3method("[",MArrayLM)
S3method("[",RGList)
@@ -70,6 +72,10 @@ S3method(summary,EListRaw)
S3method(summary,TestResults)
S3method(normalizeVSN,RGList)
S3method(normalizeVSN,EListRaw)
+S3method(plot,RGList)
+S3method(plot,MAList)
+S3method(plot,EList)
+S3method(plot,MArrayLM)
S3method(plotMA,RGList)
S3method(plotMA,MAList)
S3method(plotMA,EList)
@@ -81,9 +87,3 @@ S3method(plotDensities,RGList)
S3method(plotDensities,MAList)
S3method(plotDensities,EListRaw)
S3method(plotDensities,EList)
-S3method(roast,MAList)
-S3method(roast,EList)
-S3method(mroast,MAList)
-S3method(mroast,EList)
-S3method(camera,MAList)
-S3method(camera,EList)
diff --git a/R/barcodeplot.R b/R/barcodeplot.R
index fcb498d..a222e43 100644
--- a/R/barcodeplot.R
+++ b/R/barcodeplot.R
@@ -1,59 +1,162 @@
## BARCODEPLOT.R
-barcodeplot <- function(statistics,index,index2=NULL,labels=c("Group 1 (highest statistic)","Group 2 (lowest statistic)"),quantiles=c(-1,1),col.bars=NULL,offset.bars=!is.null(index2), ...)
-# Barcode plot of one or two gene sets.
-# Gordon Smyth and Di Wu
-# 20 October 2008. Last revised 11 Feb 2012.
+barcodeplot <- function (statistics, index, index2 = NULL, labels = c("Up", "Down"), quantiles = c(-1, 1), col.bars = NULL, worm=TRUE, span.worm=0.45, ...)
+# Barcode plot of one or two gene sets.
+# Gordon Smyth, Di Wu and Yifang Hu
+# 20 October 2008. Last revised 4 Feb 2014.
{
+# Are there up and down sets?
TWO <- !is.null(index2)
- rankstat <- rank(statistics,na.last=TRUE)
- r <- rankstat[index]
- if(TWO) r2 <- rankstat[index2]
- isna <- is.na(statistics)
- n <- sum(!isna)
- if(any(isna)) {
- rankstat <- rankstat[1:n]
- r <- r[r<=n]
- if(TWO) r2 <- r2[r2<=n]
+
+# Convert indexes to logical
+ set1 <- set2 <- rep.int(FALSE,length(statistics))
+ names(set1) <- names(set2) <- names(statistics)
+ set1[index] <- TRUE
+ if(TWO) set2[index2] <- TRUE
+
+# Sort statistics and indexes
+ ostat <- order(statistics, na.last = TRUE, decreasing=TRUE)
+ statistics <- statistics[ostat]
+ set1 <- set1[ostat]
+ if(TWO) set2 <- set2[ostat]
+
+# Trim off missing values
+ n <- sum(!is.na(statistics))
+ if(n==0L) {
+ message("No valid statistics")
+ return(invisible())
}
- if(is.null(col.bars))
- if(is.null(index2))
- col.bars=c("black","black")
- else
- col.bars=c("red","blue")
+ if (n < length(statistics)) {
+ statistics <- statistics[1:n]
+ set1 <- set1[1:n]
+ if (TWO) set2 <- set2[1:n]
+ }
+
+# Convert indexes to integer
+ r <- which(set1)
+ if (TWO) {
+ r2 <- which(set2)
+ if(!length(r2)) TWO <- FALSE
+ }
+
+# Check there is something to plot
+ if(!length(r))
+ if(TWO) {
+ r <- r2
+ set1 <- set2
+ TWO <- FALSE
+ message("Using index2 as primary set")
+ } else {
+ message("No selected genes to plot")
+ return(invisible())
+ }
+
+# Check other arguments
quantiles <- sort(quantiles)
+ if (is.null(col.bars))
+ if (TWO)
+ col.bars = c("red", "blue")
+ else
+ col.bars = c("black", "black")
- ylim <- c(-1,1)
- if(offset.bars) {
- ylim[2] <- ylim[2]+0.5
- if(!is.null(index2)) ylim[1] <- ylim[1]-0.5
+ # worm plot setting
+ ylim.worm <- ylim <- c(-1, 1)
+ ylab.worm <- ""
+ xlab.worm <- "statistics"
+
+ if(!TWO) ylim.worm <- c(0, 1)
+
+ if(worm) {
+ ylim.worm <- c(-2.1, 2.1)
+ if(!TWO) ylim.worm <- c(0, 2.1)
}
- n <- length(statistics)
- plot(1:n,xlim=c(0,n),ylim=ylim,type="n",axes=FALSE,xlab="",ylab="",...)
+
+ ylim[2] <- ylim[2] + 0.5
+ if (TWO) ylim[1] <- ylim[1] - 0.5
+
+ plot(1:n, xlim = c(0, n), ylim = ylim.worm, type = "n", axes = FALSE, xlab = xlab.worm, ylab = ylab.worm, ...)
+
npos <- sum(statistics > quantiles[2])
nneg <- sum(statistics < quantiles[1])
- rect(npos+0.5,-0.5,n-nneg+0.5,0.5,col="lightgray",border=NA)
- if(npos) rect(0.5,-0.5,npos+0.5,0.5,col="pink",border=NA)
- if(nneg) rect(n-nneg+0.5,-0.5,n+0.5,0.5,col="lightblue",border=NA)
- r <- n+1-r
+ rect.yb <- -0.5
+ if(!TWO) rect.yb <- 0
+
+ rect(npos + 0.5, rect.yb, n - nneg + 0.5, 0.5, col = "lightgray", border = NA)
+ if (npos) rect(0.5, rect.yb, npos + 0.5, 0.5, col = "pink", border = NA)
+ if (nneg) rect(n - nneg + 0.5, rect.yb, n + 0.5, 0.5, col = "lightblue", border = NA)
+
lwd <- 50/length(r)
- lwd <- min(2,lwd)
- lwd <- max(0.1,lwd)
- barlim <- ylim[2]-c(1.5,0.5)
- segments(r,barlim[1],r,barlim[2],lwd=lwd,col=col.bars[1])
+ lwd <- min(2, lwd)
+ lwd <- max(0.1, lwd)
+
+ barlim <- ylim[2] - c(1.5, 0.5)
+ segments(r, barlim[1], r, barlim[2], lwd = lwd, col = col.bars[1])
if(TWO) {
- r2 <- n+1-r2
lwd2 <- 50/length(r2)
- lwd2 <- min(2,lwd2)
- lwd2 <- max(0.1,lwd2)
- barlim2 <- ylim[1]+c(0.5,1.5)
- segments(r2,barlim2[1],r2,barlim2[2],lwd=lwd2,col=col.bars[2])
+ lwd2 <- min(2, lwd2)
+ lwd2 <- max(0.1, lwd2)
+ barlim2 <- ylim[1] + c(0.5, 1.5)
+ segments(r2, barlim2[1], r2, barlim2[2], lwd = lwd2, col = col.bars[2])
+ }
+
+ lab.at <- 0
+ if(!TWO) lab.at <- 0.5
+ axis(side = 2, at = lab.at, padj = 3, cex.axis = 0.85, labels = labels[1], tick = FALSE)
+ axis(side = 4, at = lab.at, padj = -3, cex.axis = 0.85, labels = labels[2], tick = FALSE)
+
+ # label statistics on x axis
+ prob <- (10:0)/10
+ axis(at = seq(1,n,len=11), side = 1, cex.axis = 0.7, las = 2, labels = format(quantile(statistics, p = prob), digits = 1))
+
+ # create worm
+ if(worm) {
+ # rescale x to new range
+ rescale <- function(x, newrange, oldrange=range(x)) {
+ newrange[1] + (x-oldrange[1]) / (oldrange[2]-oldrange[1]) * (newrange[2] - newrange[1])
+ }
+
+ # calculate enrichment
+ ave.enrich1 <- length(r)/n
+ worm1 <- tricubeMovingAverage(set1,span=span.worm)/ave.enrich1
+
+ if(TWO) {
+ ave.enrich2 <- length(r2)/n
+ worm2 <- tricubeMovingAverage(-set2,span=span.worm)/ave.enrich2
+ }
+
+ # rescale worm
+ max.worm1 <- max(worm1)
+ r.worm1 <- c(0,max.worm1)
+ worm1.scale <- rescale(worm1, newrange=c(1.1,2.1), oldrange=r.worm1)
+
+ if(TWO) {
+ min.worm2 <- min(worm2)
+ r.worm2 <- c(min.worm2,0)
+ worm2.scale <- rescale(worm2, newrange=c(-2.1,-1.1), oldrange=r.worm2)
+ }
+
+ # plot worms
+
+ if(!TWO) {
+ lines(x = 1:n, y = worm1.scale, col = col.bars[1], lwd = 2)
+ abline(h = rescale(1,newrange=c(1.1,2.1),oldrange=r.worm1), lty=2)
+ axis(side = 2, at = c(1.1, 2.1), cex.axis = 0.8, labels = c(0, format(max.worm1, digits = 2)))
+ axis(side = 2, labels = "Enrichment", at = 1.6, padj = -0.6, tick = FALSE, cex.axis = 0.8)
+ }
+
+ if(TWO) {
+ lines(x = 1:n, y = worm1.scale, col = col.bars[1], lwd = 2)
+ abline(h = rescale(1,newrange=c(1.1,2.1),oldrange=r.worm1), lty=2)
+ lines(x = 1:n, y = worm2.scale, col = col.bars[2], lwd = 2)
+ abline(h = rescale(-1,newrange=c(-2.1,-1.1),oldrange=r.worm2), lty=2)
+ axis(side = 2, at = c(1.1, 2.1), cex.axis = 0.7, labels = c(0, format(max.worm1, digits = 2)))
+ axis(side = 2, at = c(-1.1, -2.1), cex.axis = 0.7, labels = c(0, format(-min.worm2, digits = 2)))
+ axis(side = 2, labels = "Enrichment", at = -1.6, tick = FALSE, padj = -0.6, cex.axis = 0.7)
+ axis(side = 2, labels = "Enrichment", at = 1.6, tick = FALSE, padj = -0.6, cex.axis = 0.7)
+ }
}
- mtext(labels[1],side=2,line=-0.5,col="black")
- mtext(labels[2],side=4,line=-1,col="black")
invisible()
}
-
diff --git a/R/beadCountWeights.R b/R/beadCountWeights.R
index 1a8a6b8..73d7de1 100644
--- a/R/beadCountWeights.R
+++ b/R/beadCountWeights.R
@@ -1,18 +1,26 @@
-beadCountWeights <- function(y, x, design=NULL, bead.stdev=NULL, nbeads=NULL, array.cv=TRUE, scale=FALSE)
+beadCountWeights <- function(y, x, design=NULL, bead.stdev=NULL, bead.stderr=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.
+# 4 August 2010. Last modified 19 Dec 2013.
{
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")
- }
+ if(is.null(nbeads)) {
+ nbeads <- y$other$Avg_NBEADS
+ if(is.null(nbeads)) stop("NBEADS not found in data object")
+ }
+ if(is.null(bead.stdev)) {
+ if (is.null(bead.stderr)) {
+ if (is.null(y$other$BEAD_STDEV)) {
+ y$other$BEAD_STDEV <- y$other$BEAD_STDERR*sqrt(nbeads)
+ } else {
+ bead.stdev <- y$other$BEAD_STDEV
+ }
+ } else {
+ bead.stdev <- bead.stderr*sqrt(nbeads)
+ }
+ if(is.null(bead.stdev)) stop("BEAD_STDEV and BEAD_STDERR are missing. At least one is required.")
+ }
P <- nrow(E)
A <- ncol(E)
if(nrow(E.raw) != P) stop("dimensions don't match")
diff --git a/R/diffSplice.R b/R/diffSplice.R
new file mode 100644
index 0000000..c2981a5
--- /dev/null
+++ b/R/diffSplice.R
@@ -0,0 +1,210 @@
+diffSplice <- function(fit,geneid,exonid=NULL,verbose=TRUE)
+# Test for splicing variants between conditions
+# using linear model fit of exon data.
+# Charity Law and Gordon Smyth
+# Created 13 Dec 2013. Last modified 18 Feb 2014.
+{
+ exon.genes <- fit$genes
+ if(is.null(exon.genes)) exon.genes <- data.frame(ExonID=1:nrow(fit))
+ if(length(geneid)==1) {
+ genecolname <- as.character(geneid)
+ geneid <- exon.genes[[genecolname]]
+ } else {
+ exon.genes$GeneID <- geneid
+ genecolname <- "GeneID"
+ }
+ if(!is.null(exonid))
+ if(length(exonid)==1) {
+ exoncolname <- as.character(exonid)
+ exonid <- exon.genes[[exoncolname]]
+ } else {
+ exon.genes$ExonID <- exonid
+ exoncolname <- "ExonID"
+ }
+ else
+ exoncolname <- NULL
+
+# Sort by geneid
+ if(is.null(exonid))
+ o <- order(geneid)
+ else
+ o <- order(geneid,exonid)
+ geneid <- geneid[o]
+ exon.genes <- exon.genes[o,,drop=FALSE]
+ exon.coefficients <- fit$coefficients[o,,drop=FALSE]
+ exon.stdev.unscaled <- fit$stdev.unscaled[o,,drop=FALSE]
+ exon.df.residual <- fit$df.residual[o]
+ exon.s2 <- fit$sigma[o]^2
+
+# Count exons and get genewise variances
+ exon.stat <- cbind(1,exon.df.residual,exon.s2)
+ gene.sum <- rowsum(exon.stat,geneid,reorder=FALSE)
+ gene.nexons <- gene.sum[,1]
+ gene.df.residual <- gene.sum[,2]
+ gene.s2 <- gene.sum[,3] / gene.sum[,1]
+ if(verbose) {
+ cat("Total number of exons: ", length(geneid), "\n")
+ cat("Total number of genes: ", length(gene.nexons), "\n")
+ cat("Number of genes with 1 exon: ", sum(gene.nexons==1), "\n")
+ cat("Mean number of exons in a gene: ", round(mean(gene.nexons),0), "\n")
+ cat("Max number of exons in a gene: ", max(gene.nexons), "\n")
+ }
+
+# Posterior genewise variances
+ squeeze <- squeezeVar(var=gene.s2, df=gene.df.residual)
+
+# Remove genes with only 1 exon
+ gene.keep <- gene.nexons>1
+ ngenes <- sum(gene.keep)
+ if(ngenes==0) stop("No genes with more than one exon")
+
+ exon.keep <- rep(gene.keep,gene.nexons)
+ geneid <- geneid[exon.keep]
+ exon.genes <- exon.genes[exon.keep,,drop=FALSE]
+ exon.coefficients <- exon.coefficients[exon.keep,,drop=FALSE]
+ exon.stdev.unscaled <- exon.stdev.unscaled[exon.keep,,drop=FALSE]
+ exon.df.residual <- exon.df.residual[exon.keep]
+
+ gene.nexons <- gene.nexons[gene.keep]
+ gene.df.test <- gene.nexons-1
+ gene.df.residual <- gene.df.residual[gene.keep]
+ gene.df.total <- gene.df.residual+squeeze$df.prior
+ gene.df.total <- pmin(gene.df.total,sum(gene.df.residual))
+ gene.s2.post <- squeeze$var.post[gene.keep]
+
+# Genewise betas
+ u2 <- 1/exon.stdev.unscaled^2
+ u2.rowsum <- rowsum(u2,geneid,reorder=FALSE)
+ gene.betabar <- rowsum(exon.coefficients*u2,geneid,reorder=FALSE) / u2.rowsum
+
+# T-statistics for exon-level tests
+ g <- rep(1:ngenes,times=gene.nexons)
+ exon.coefficients <- exon.coefficients-gene.betabar[g,,drop=FALSE]
+ exon.t <- exon.coefficients / exon.stdev.unscaled / sqrt(gene.s2.post[g])
+ gene.F <- rowsum(exon.t^2,geneid,reorder=FALSE) / gene.df.test
+ exon.1mleverage <- 1 - (u2 / u2.rowsum[g,,drop=FALSE])
+ exon.coefficients <- exon.coefficients / exon.1mleverage
+ exon.t <- exon.t / sqrt(exon.1mleverage)
+ exon.p.value <- 2 * pt(abs(exon.t), df=gene.df.total[g], lower.tail=FALSE)
+ gene.F.p.value <- pf(gene.F, df1=gene.df.test, df2=gene.df.total, lower.tail=FALSE)
+
+# Exon level output
+ out <- new("MArrayLM",list())
+ out$genes <- exon.genes
+ out$genecolname <- genecolname
+ out$exoncolname <- exoncolname
+ out$coefficients <- exon.coefficients
+ out$t <- exon.t
+ out$p.value <- exon.p.value
+
+# Gene level output
+ out$gene.df.prior <- squeeze$df.prior
+ out$gene.df.residual <- gene.df.residual
+ out$gene.df.total <- gene.df.total
+ out$gene.s2.post <- gene.s2.post
+ out$gene.F <- gene.F
+ out$gene.F.p.value <- gene.F.p.value
+
+# Which columns of exon.genes contain gene level annotation?
+ exon.lastexon <- cumsum(gene.nexons)
+ exon.firstexon <- exon.lastexon-gene.nexons+1
+ no <- logical(nrow(exon.genes))
+ isdup <- vapply(exon.genes,duplicated,no)[-exon.firstexon,,drop=FALSE]
+ isgenelevel <- apply(isdup,2,all)
+ out$gene.genes <- exon.genes[exon.lastexon,isgenelevel, drop=FALSE]
+ out$gene.genes$NExons <- gene.nexons
+
+ out
+}
+
+topSplice <- function(fit, coef=ncol(fit), level="exon", number=10, FDR=1)
+# Collate voomex results in data.frame, ordered from most significant at top
+# Gordon Smyth
+# Created 18 Dec 2013. Last modified 17 March 2014.
+{
+ coef <- coef[1]
+ exon <- match.arg(level,c("exon","gene"))
+ if(level=="exon") {
+ number <- min(number,nrow(fit$coefficients))
+ P <- fit$p.value[,coef]
+ BH <- p.adjust(P, method="BH")
+ if(FDR<1) number <- min(number,sum(BH<FDR))
+ o <- order(P)[1:number]
+ data.frame(fit$genes[o,,drop=FALSE],logFC=fit$coefficients[o,coef],t=fit$t[o,coef],P.Value=P[o],FDR=BH[o])
+ } else {
+ number <- min(number,nrow(fit$gene.F))
+ P <- fit$gene.F.p.value[,coef]
+ BH <- p.adjust(P, method="BH")
+ if(FDR<1) number <- min(number,sum(BH<FDR))
+ o <- order(P)[1:number]
+ data.frame(fit$gene.genes[o,,drop=FALSE],F=fit$gene.F[o,coef],P.Value=P[o],FDR=BH[o])
+ }
+}
+
+plotSplice <- function(fit, coef=ncol(fit), geneid=NULL, rank=1L, FDR = 0.05)
+# Plot exons of most differentially spliced gene
+# Gordon Smyth and Yifang Hu
+# Created 3 Jan 2014. Last modified 19 March 2014.
+{
+ # Gene labelling including gene symbol
+ genecolname <- fit$genecolname
+ genelab <- grep(paste0(genecolname,"|Symbol|symbol"), colnames(fit$gene.genes), value = T)
+
+ if(is.null(geneid)) {
+ if(rank==1L)
+ i <- which.min(fit$gene.F.p.value[,coef])
+ else
+ i <- order(fit$gene.F.p.value[,coef])[rank]
+
+ geneid <- paste(fit$gene.genes[i,genelab], collapse = ".")
+
+ } else {
+ i <- which(fit$gene.genes[,fit$genecolname]==geneid)
+
+ geneid <- paste(fit$gene.genes[i,genelab], collapse = ".")
+
+ if(!length(i)) stop(paste("geneid",geneid,"not found"))
+ }
+ exon.lastexon <- cumsum(fit$gene.genes$NExons[1:i])
+ j <- (exon.lastexon[i-1]+1):exon.lastexon[i]
+
+ exoncolname <- fit$exoncolname
+
+ if(is.null(exoncolname)){
+
+ plot(fit$coefficients[j,coef], xlab = "Exon", ylab = "logFC (this exon vs rest)", main = geneid, type = "b")
+
+ }
+
+ # Plot exons and mark exon ids on the axis
+ if(!is.null(exoncolname)) {
+
+ exon.id <- fit$genes[j, exoncolname]
+ xlab <- paste("Exon", exoncolname, sep = " ")
+
+ plot(fit$coefficients[j, coef], xlab = "", ylab = "logFC (this exon vs rest)", main = geneid,type = "b", xaxt = "n")
+ axis(1, at = 1:length(j), labels = exon.id, las = 2, cex.axis = 0.6)
+ mtext(xlab, side = 1, padj = 5.2)
+
+ # Mark the topSpliced exons
+ top <- topSplice(fit, coef = coef, number = Inf, level = "exon", FDR = FDR)
+ m <- which(top[,genecolname] %in% fit$gene.genes[i,genecolname])
+
+ if(length(m) > 0){
+
+ if(length(m) == 1) cex <- 1.5 else{
+
+ abs.fdr <- abs(log10(top$FDR[m]))
+ from <- range(abs.fdr)
+ to <- c(1,2)
+ cex <- (abs.fdr - from[1])/diff(from) * diff(to) + to[1]
+
+ }
+ mark <- match(top[m, exoncolname], exon.id)
+ points((1:length(j))[mark], fit$coefficients[j[mark], coef], col = "red", pch = 16, cex = cex)
+ }
+ }
+
+ abline(h=0,lty=2)
+
+}
\ No newline at end of file
diff --git a/R/dups.R b/R/dups.R
index 9657a29..1118a1d 100755
--- a/R/dups.R
+++ b/R/dups.R
@@ -99,6 +99,7 @@ duplicateCorrelation <- function(object,design=NULL,ndups=2,spacing=1,block=NULL
require( "statmod" ) # need mixedModel2Fit function
rho <- rep(NA,ngenes)
+ nafun <- function(e) NA
for (i in 1:ngenes) {
y <- drop(M[i,])
o <- is.finite(y)
@@ -111,10 +112,10 @@ duplicateCorrelation <- function(object,design=NULL,ndups=2,spacing=1,block=NULL
Z <- model.matrix(~0+A)
if(!is.null(weights)) {
w <- drop(weights[i,])[o]
- s <- mixedModel2Fit(y,X,Z,w,only.varcomp=TRUE,maxit=20)$varcomp
+ s <- tryCatch(mixedModel2Fit(y,X,Z,w,only.varcomp=TRUE,maxit=20)$varcomp,error=nafun)
} else
- s <- mixedModel2Fit(y,X,Z,only.varcomp=TRUE,maxit=20)$varcomp
- rho[i] <- s[2]/sum(s)
+ s <- tryCatch(mixedModel2Fit(y,X,Z,only.varcomp=TRUE,maxit=20)$varcomp,error=nafun)
+ if(!is.na(s[1])) rho[i] <- s[2]/sum(s)
}
}
arho <- atanh(pmax(-1,rho))
diff --git a/R/fitFDistRobustly.R b/R/fitFDistRobustly.R
index 3165b07..1280db5 100644
--- a/R/fitFDistRobustly.R
+++ b/R/fitFDistRobustly.R
@@ -3,7 +3,7 @@ fitFDistRobustly <- function(x,df1,covariate=NULL,winsor.tail.p=c(0.05,0.1),trac
# given the first degrees of freedom, using first and second
# moments of Winsorized z-values
# Gordon Smyth and Belinda Phipson
-# 8 Sept 2002. Last revised 25 April 2013.
+# 8 Sept 2002. Last revised 20 February 2014.
{
# Check x
n <- length(x)
@@ -29,7 +29,7 @@ fitFDistRobustly <- function(x,df1,covariate=NULL,winsor.tail.p=c(0.05,0.1),trac
ok <- !is.na(x) & is.finite(df1) & (df1 > 1e-6)
notallok <- !all(ok)
if(notallok) {
- scale <- df2.shrunk <- x
+ df2.shrunk <- x
x <- x[ok]
if(length(df1)>1) df1 <- df1[ok]
if(!is.null(covariate)) {
@@ -39,8 +39,13 @@ fitFDistRobustly <- function(x,df1,covariate=NULL,winsor.tail.p=c(0.05,0.1),trac
fit <- Recall(x=x,df1=df1,covariate=covariate,winsor.tail.p=winsor.tail.p,trace=trace)
df2.shrunk[ok] <- fit$df2.shrunk
df2.shrunk[!ok] <- fit$df2
- scale[ok] <- fit$scale
- scale[!ok] <- exp(approx(covariate,log(fit$scale),xout=covariate2,rule=2)$y)
+ if(is.null(covariate))
+ scale <- fit$scale
+ else {
+ scale <- x
+ scale[ok] <- fit$scale
+ scale[!ok] <- exp(approx(covariate,log(fit$scale),xout=covariate2,rule=2)$y)
+ }
return(list(scale=scale,df2=fit$df2,df2.shrunk=df2.shrunk))
}
diff --git a/R/geneset-camera.R b/R/geneset-camera.R
index 696ac34..5912bb5 100644
--- a/R/geneset-camera.R
+++ b/R/geneset-camera.R
@@ -19,49 +19,44 @@ interGeneCorrelation <- function(y, design)
list(vif=vif,correlation=correlation)
}
-camera <- function(y,index,design=NULL,contrast=ncol(design),weights=NULL,use.ranks=FALSE,allow.neg.cor=TRUE,trend.var=FALSE,sort=TRUE)
-UseMethod("camera")
+camera <- function(y,...) UseMethod("camera")
-camera.EList <- function(y,index,design=NULL,contrast=ncol(design),weights=NULL,use.ranks=FALSE,allow.neg.cor=TRUE,trend.var=FALSE,sort=TRUE)
-# Gordon Smyth
-# Created 4 Jan 2013. Last modified 22 Jan 2013
-{
- if(is.null(design)) design <- y$design
- if(is.null(weights)) weights <- y$weights
- y <- y$E
- camera(y=y,index=index,design=design,contrast=contrast,weights=weights,use.ranks=use.ranks,allow.neg.cor=allow.neg.cor,trend.var=trend.var,sort=sort)
-}
-
-camera.MAList <- function(y,index,design=NULL,contrast=ncol(design),weights=NULL,use.ranks=FALSE,allow.neg.cor=TRUE,trend.var=FALSE,sort=TRUE)
-# Gordon Smyth
-# Created 4 Jan 2013. Last modified 22 Jan 2013
-{
- if(is.null(design)) design <- y$design
- if(is.null(weights)) weights <- y$weights
- y <- y$M
- camera(y=y,index=index,design=design,contrast=contrast,weights=weights,use.ranks=use.ranks,allow.neg.cor=allow.neg.cor,trend.var=trend.var,sort=sort)
-}
-
-camera.default <- function(y,index,design=NULL,contrast=ncol(design),weights=NULL,use.ranks=FALSE,allow.neg.cor=TRUE,trend.var=FALSE,sort=TRUE)
+camera.default <- function(y,index,design=NULL,contrast=ncol(design),weights=NULL,use.ranks=FALSE,allow.neg.cor=TRUE,trend.var=FALSE,sort=TRUE,...)
# Competitive gene set test allowing for correlation between genes
# Gordon Smyth and Di Wu
-# Created 2007. Last modified 22 Jan 2013
+# Created 2007. Last modified 25 Feb 2014
{
-# Check y
- y <- as.matrix(y)
- G <- nrow(y)
- n <- ncol(y)
+# Extract components from y
+ y <- getEAWP(y)
+ G <- nrow(y$exprs)
+ n <- ncol(y$exprs)
# Check index
if(!is.list(index)) index <- list(set1=index)
# Check design
- if(is.null(design)) stop("no design matrix")
+ if(is.null(design)) design <- y$design
+ if(is.null(design))
+ design <- matrix(1,n,1)
+ else {
+ design <- as.matrix(design)
+ if(mode(design) != "numeric") stop("design must be a numeric matrix")
+ }
+ if(nrow(design) != n) stop("row dimension of design matrix must match column dimension of data")
+ ne <- nonEstimable(design)
+ if(!is.null(ne)) cat("Coefficients not estimable:",paste(ne,collapse=" "),"\n")
+
p <- ncol(design)
df.residual <- n-p
df.camera <- min(df.residual,G-2)
# Check weights
+ if(is.null(weights)) weights <- y$weights
+
+# Reduce to numeric expression matrix
+ y <- y$exprs
+
+# Check weights
if(!is.null(weights)) {
if(any(weights<=0)) stop("weights must be positive")
if(length(weights)==n) {
diff --git a/R/geneset-roast.R b/R/geneset-roast.R
index 85e27f0..5f1e78c 100644
--- a/R/geneset-roast.R
+++ b/R/geneset-roast.R
@@ -11,47 +11,36 @@ setMethod("show","Roast",
function(object) print(object$p.value)
)
-roast <- 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)
-UseMethod("roast")
+roast <- function(y,...) UseMethod("roast")
-roast.EList <- 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)
-# Gordon Smyth
-# Created 4 Jan 2013. Last modified 22 Jan 2013.
-{
- if(is.null(design)) design <- y$design
- if(is.null(weights) && is.null(array.weights)) weights <- y$weights
- y <- as.matrix(y)
- roast(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)
-}
-
-roast.MAList <- 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)
-# Gordon Smyth
-# Created 4 Jan 2013. Last modified 22 Jan 2013.
-{
- if(is.null(design)) design <- y$design
- if(is.null(weights) && is.null(array.weights)) weights <- y$weights
- y <- as.matrix(y)
- roast(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)
-}
-
-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)
+roast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.statistic="mean",gene.weights=NULL,array.weights=NULL,weights=NULL,block=NULL,correlation,var.prior=NULL,df.prior=NULL,trend.var=FALSE,nrot=999,approx.zscore=TRUE,...)
# Rotation gene set testing for linear models
# Gordon Smyth and Di Wu
-# Created 24 Apr 2008. Last modified 9 Jan 2014.
+# Created 24 Apr 2008. Last modified 18 June 2014.
{
-# Check y
- y <- as.matrix(y)
- ngenes <- nrow(y)
- 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))
+
+# Extract components from y
+ y <- getEAWP(y)
+ ngenes <- nrow(y$exprs)
+ n <- ncol(y$exprs)
+
+# Check index
if(is.null(index)) index <- rep.int(TRUE,ngenes)
# Check design
- if(is.null(design)) stop("no design matrix")
- design <- as.matrix(design)
+ if(is.null(design)) design <- y$design
+ if(is.null(design))
+ design <- matrix(1,n,1)
+ else {
+ design <- as.matrix(design)
+ if(mode(design) != "numeric") stop("design must be a numeric matrix")
+ }
if(nrow(design) != n) stop("row dimension of design matrix must match column dimension of data")
+ ne <- nonEstimable(design)
+ if(!is.null(ne)) cat("Coefficients not estimable:",paste(ne,collapse=" "),"\n")
+
p <- ncol(design)
p0 <- p-1L
d <- n-p
@@ -63,24 +52,28 @@ roast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.sta
if(!is.null(var.prior) && var.prior<0) stop("var.prior must be non-negative")
if(!is.null(df.prior) && df.prior<0) stop("df.prior must be non-negative")
-# Check observation weights
- if(!is.null(weights)) {
- weights <- as.matrix(weights)
- dimw <- dim(weights)
- if(dimw[1]!=ngenes || dimw[2]!=n) stop("weights must have same dimensions as y")
- if(any(weights<=0)) stop("weights must be positive")
- }
-
# Check array weights
if(!is.null(array.weights)) {
if(length(array.weights) != n) stop("Length of array.weights doesn't match number of array")
if(any(array.weights <= 0)) stop("array.weights must be positive")
- if(!is.null(weights)) {
- weights <- .matvec(weights,array.weights)
+ }
+
+# Check observational weights
+ if(is.null(weights) && is.null(array.weights)) weights <- y$weights
+ if(!is.null(weights)) {
+ weights <- as.matrix(weights)
+ dimw <- dim(weights)
+ if(dimw[1] != ngenes || dimw[2] != n) stop("weights must have same dimensions as y")
+ if(any(weights <= 0)) stop("weights must be positive")
+ if(!is.null(array.weights)) {
+ weights <- .matvec(weights, array.weights)
array.weights <- NULL
}
}
+# Reduce to numeric expression matrix
+ y <- y$exprs
+
# Divide out array weights
if(!is.null(array.weights)) {
sw <- sqrt(array.weights)
@@ -104,7 +97,7 @@ roast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.sta
}
# Check contrast
- if(all(contrast==0)) stop("contrast all zero")
+ if(all(contrast == 0)) stop("contrast all zero")
# Reform design matrix so that contrast is last coefficient
if(length(contrast) == 1) {
@@ -187,14 +180,14 @@ roast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.sta
B <- Y[1,]
modt <- signc*B/sd.post
- statobs <- p <- rep(0,3)
- names(statobs) <- names(p) <- c("down","up","mixed")
- statrot <- array(0,c(nrot,3),dimnames=list(NULL,names(p)))
+ statobs <- p <- rep(0,4)
+ names(statobs) <- names(p) <- c("down","up","upordown","mixed")
+ statrot <- array(0,c(nrot,4),dimnames=list(NULL,names(p)))
# Convert to z-scores
- modt <- zscoreT(modt,df=d0+d)
+ modt <- zscoreT(modt,df=d0+d,approx=approx.zscore)
-# Active proportions
+# Active proportions
if(!is.null(gene.weights)) {
lgw <- length(gene.weights)
if(lgw > nset && lgw==ngenes) {
@@ -220,48 +213,31 @@ roast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.sta
else
sdr.post <- sqrt(s02)
modtr <- signc*Br/sdr.post
- modtr <- zscoreT(modtr,df=d0+d)
-
- if(set.statistic=="msq") {
-# Observed statistics
- modt2 <- modt^2
- if(!is.null(gene.weights)) {
- modt2 <- abs(gene.weights)*modt2
- modt <- gene.weights*modt
- }
- statobs["mixed"] <- mean(modt2)
- statobs["up"] <- sum(modt2[modt > 0])/nset
- statobs["down"] <- sum(modt2[modt < 0])/nset
-# Simulated statistics
- if(!is.null(gene.weights)) {
- gene.weights <- sqrt(abs(gene.weights))
- modtr <- t(gene.weights*t(modtr))
- }
- statrot[,"mixed"] <- rowMeans(modtr^2)
- statrot[,"up"] <- rowMeans(pmax(modtr,0)^2)
- statrot[,"down"] <- rowMeans(pmax(-modtr,0)^2)
-# p-values
- p <- (rowSums(t(statrot) >= statobs)+1)/(nrot+1)
- }
+ modtr <- zscoreT(modtr,df=d0+d,approx=approx.zscore)
- if(set.statistic=="mean50") {
+ switch(set.statistic,
+ "mean" = {
# Observed statistics
if(!is.null(gene.weights)) modt <- gene.weights*modt
- half <- nset %/% 2L +1L
- meanTop <- function(x,n) mean(sort(x,partial=n)[n:length(x)])
- statobs["mixed"] <- meanTop(abs(modt),half)
- statobs["up"] <- meanTop(modt,half)
- statobs["down"] <- meanTop(-modt,half)
+ m <- mean(modt)
+ statobs["down"] <- -m
+ statobs["up"] <- m
+ statobs["mixed"] <- mean(abs(modt))
# Simulated statistics
if(!is.null(gene.weights)) modtr <- t(gene.weights*t(modtr))
- statrot[,"mixed"] <- apply(abs(modtr),1,meanTop,n=half)
- statrot[,"up"] <- apply(modtr,1,meanTop,n=half)
- statrot[,"down"] <- apply(-modtr,1,meanTop,n=half)
+ m <- rowMeans(modtr)
+ statrot[,"down"] <- -m
+ statrot[,"up"] <- m
+ statrot[,"mixed"] <- rowMeans(abs(modtr))
# p-values
- p <- (rowSums(t(statrot) >= statobs) + 1)/(nrot + 1)
- }
-
- if(set.statistic=="floormean") {
+ p["down"] <- sum(statrot[,c("down","up")] > statobs["down"])
+ p["up"] <- sum(statrot[,c("down","up")] > statobs["up"])
+ p["upordown"] <- min(p[c("down","up")])
+ p["mixed"] <- sum(statrot[,c("mixed")] > statobs["mixed"])
+ p <- (p+1) / (c(2,2,1,1)*nrot + 1)
+ },
+
+ "floormean" = {
# Observed statistics
chimed <- qchisq(0.5,df=1)
amodt <- pmax(abs(modt),chimed)
@@ -269,107 +245,154 @@ roast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.sta
amodt <- gene.weights*amodt
modt <- gene.weights*modt
}
- statobs["mixed"] <- mean(amodt)
- statobs["up"] <- mean(pmax(modt,0))
statobs["down"] <- mean(pmax(-modt,0))
+ statobs["up"] <- mean(pmax(modt,0))
+ statobs["upordown"] <- max(statobs[c("down","up")])
+ statobs["mixed"] <- mean(amodt)
# Simulated statistics
amodtr <- pmax(abs(modtr),chimed)
if(!is.null(gene.weights)) {
amodtr <- t(gene.weights*t(amodtr))
modtr <- t(gene.weights*t(modtr))
}
- statrot[,"mixed"] <- rowMeans(amodtr)
- statrot[,"up"] <- rowMeans(pmax(modtr,0))
statrot[,"down"] <- rowMeans(pmax(-modtr,0))
+ statrot[,"up"] <- rowMeans(pmax(modtr,0))
+ i <- statrot[,"up"] > statrot[,"down"]
+ statrot[i,"upordown"] <- statrot[i,"up"]
+ statrot[!i,"upordown"] <- statrot[!i,"down"]
+ statrot[,"mixed"] <- rowMeans(amodtr)
# p-values
- p <- (rowSums(t(statrot) >= statobs) + 1)/(nrot + 1)
- }
-
- if(set.statistic=="mean") {
+ p["down"] <- sum(statrot[,c("down","up")] > statobs["down"])
+ p["up"] <- sum(statrot[,c("down","up")] > statobs["up"])
+ p["upordown"] <- sum(statrot[,c("upordown")] > statobs["upordown"])
+ p["mixed"] <- sum(statrot[,c("mixed")] > statobs["mixed"])
+ p <- (p+1) / (c(2,2,1,1)*nrot + 1)
+ },
+
+ "mean50" = {
+ if(nset%%2L == 0L) {
+ half1 <- nset %/% 2L
+ half2 <- half1 + 1L
+ } else {
+ half1 <- half2 <- nset %/% 2L + 1L
+ }
# Observed statistics
if(!is.null(gene.weights)) modt <- gene.weights*modt
- m <- mean(modt)
- statobs["mixed"] <- mean(abs(modt))
- statobs["up"] <- m
- statobs["down"] <- -m
+ s <- sort(modt,partial=half2)
+ statobs["down"] <- -mean(s[1:half1])
+ statobs["up"] <- mean(s[half2:nset])
+ statobs["upordown"] <- max(statobs[c("down","up")])
+ s <- sort(abs(modt),partial=half2)
+ statobs["mixed"] <- mean(s[half2:nset])
# Simulated statistics
if(!is.null(gene.weights)) modtr <- t(gene.weights*t(modtr))
- m <- rowMeans(modtr)
- statrot[,"mixed"] <- rowMeans(abs(modtr))
- statrot[,"up"] <- m
- statrot[,"down"] <- -m
+ for (g in 1L:nrot) {
+ s <- sort(modtr[g,],partial=half2)
+ statrot[g,"down"] <- -mean(s[1:half1])
+ statrot[g,"up"] <- mean(s[half2:nset])
+ statrot[g,"upordown"] <- max(statrot[g,c("down","up")])
+ s <- sort(abs(modtr[g,]),partial=half2)
+ statrot[g,"mixed"] <- mean(s[half2:nset])
+ }
# p-values
- p <- (rowSums(t(statrot) >= statobs) + 1)/(nrot + 1)
- }
+ p["down"] <- sum(statrot[,c("down","up")] > statobs["down"])
+ p["up"] <- sum(statrot[,c("down","up")] > statobs["up"])
+ p["upordown"] <- sum(statrot[,c("upordown")] > statobs["upordown"])
+ p["mixed"] <- sum(statrot[,c("mixed")] > statobs["mixed"])
+ p <- (p+1) / (c(2,2,1,1)*nrot + 1)
+ },
+
+ "msq" = {
+# Observed statistics
+ modt2 <- modt^2
+ if(!is.null(gene.weights)) {
+ modt2 <- abs(gene.weights)*modt2
+ modt <- gene.weights*modt
+ }
+ statobs["down"] <- sum(modt2[modt < 0])/nset
+ statobs["up"] <- sum(modt2[modt > 0])/nset
+ statobs["upordown"] <- max(statobs[c("down","up")])
+ statobs["mixed"] <- mean(modt2)
+# Simulated statistics
+ if(!is.null(gene.weights)) {
+ gene.weights <- sqrt(abs(gene.weights))
+ modtr <- t(gene.weights*t(modtr))
+ }
+ statrot[,"down"] <- rowMeans(pmax(-modtr,0)^2)
+ statrot[,"up"] <- rowMeans(pmax(modtr,0)^2)
+ i <- statrot[,"up"] > statrot[,"down"]
+ statrot[i,"upordown"] <- statrot[i,"up"]
+ statrot[!i,"upordown"] <- statrot[!i,"down"]
+ statrot[,"mixed"] <- rowMeans(modtr^2)
+# p-values
+ p["down"] <- sum(statrot[,c("down","up")] > statobs["down"])
+ p["up"] <- sum(statrot[,c("down","up")] > statobs["up"])
+ p["upordown"] <- sum(statrot[,c("upordown")] > statobs["upordown"])
+ p["mixed"] <- sum(statrot[,c("mixed")] > statobs["mixed"])
+ p <- (p+1) / (c(2,2,1,1)*nrot + 1)
+ })
# Output
- out <- data.frame(c(r2,r1,r1+r2),p)
- dimnames(out) <- list(c("Down","Up","Mixed"),c("Active.Prop","P.Value"))
+ out <- data.frame(c(r2,r1,max(r1,r2),r1+r2),p)
+ dimnames(out) <- list(c("Down","Up","UpOrDown","Mixed"),c("Active.Prop","P.Value"))
new("Roast",list(p.value=out,var.prior=s02,df.prior=d0,ngenes.in.set=nset))
}
-mroast <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.statistic="mean",gene.weights=NULL,array.weights=NULL,weights=NULL,block=NULL,correlation,var.prior=NULL,df.prior=NULL,trend.var=FALSE,nrot=999,adjust.method="BH",midp=TRUE,sort="directional")
-UseMethod("mroast")
-
-mroast.EList <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.statistic="mean",gene.weights=NULL,array.weights=NULL,weights=NULL,block=NULL,correlation,var.prior=NULL,df.prior=NULL,trend.var=FALSE,nrot=999,adjust.method="BH",midp=TRUE,sort="directional")
-# Gordon Smyth
-# Created 8 Jan 2013.
-{
- if(is.null(design)) design <- y$design
- if(is.null(weights) && is.null(array.weights)) weights <- y$weights
- y <- as.matrix(y)
- 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)
-}
-
-mroast.MAList <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.statistic="mean",gene.weights=NULL,array.weights=NULL,weights=NULL,block=NULL,correlation,var.prior=NULL,df.prior=NULL,trend.var=FALSE,nrot=999,adjust.method="BH",midp=TRUE,sort="directional")
-# Gordon Smyth
-# Created 8 Jan 2013.
-{
- if(is.null(design)) design <- y$design
- if(is.null(weights) && is.null(array.weights)) weights <- y$weights
- y <- as.matrix(y)
- 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)
-}
+mroast <- function(y,...) UseMethod("mroast")
-mroast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.statistic="mean",gene.weights=NULL,array.weights=NULL,weights=NULL,block=NULL,correlation,var.prior=NULL,df.prior=NULL,trend.var=FALSE,nrot=999,adjust.method="BH",midp=TRUE,sort="directional")
+mroast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.statistic="mean",gene.weights=NULL,array.weights=NULL,weights=NULL,block=NULL,correlation,var.prior=NULL,df.prior=NULL,trend.var=FALSE,nrot=999,approx.zscore=TRUE,adjust.method="BH",midp=TRUE,sort="directional",...)
# Rotation gene set testing with multiple sets
# Gordon Smyth and Di Wu
-# Created 28 Jan 2010. Last revised 9 Jan 2014.
+# Created 28 Jan 2010. Last revised 18 Feb 2014.
{
-# Check y
- y <- as.matrix(y)
- ngenes <- nrow(y)
- n <- ncol(y)
+# Extract components from y
+ y <- getEAWP(y)
+ ngenes <- nrow(y$exprs)
+ n <- ncol(y$exprs)
# Check index
- if(is.null(index)) index <- rep(TRUE,nrow(y))
+ if(is.null(index)) index <- rep(TRUE,ngenes)
if(!is.list(index)) index <- list(set = index)
nsets <- length(index)
if(is.null(names(index))) names(index) <- paste("set",1:nsets,sep="")
-# Check design
- if(is.null(design)) stop("No design matrix")
- design <- as.matrix(design)
+# Check design matrix
+ if(is.null(design)) design <- y$design
+ if(is.null(design))
+ design <- matrix(1,n,1)
+ else {
+ design <- as.matrix(design)
+ if(mode(design) != "numeric") stop("design must be a numeric matrix")
+ }
+ if(nrow(design) != n) stop("row dimension of design matrix must match column dimension of data")
+ ne <- nonEstimable(design)
+ if(!is.null(ne)) cat("Coefficients not estimable:",paste(ne,collapse=" "),"\n")
# Check gene.weights
if(!is.null(gene.weights)) if(length(gene.weights) != ngenes) stop("gene.weights must have length equal to nrow(y)")
# Check array.weights
if(!is.null(array.weights)) {
- if(length(array.weights) != ncol(y)) stop("array.weights wrong length")
- if(any(array.weights<=0)) stop("array.weights must be positive")
+ if(length(array.weights) != n) stop("array.weights wrong length")
+ if(any(array.weights <= 0)) stop("array.weights must be positive")
}
# Check weights
+ if(is.null(weights) && is.null(array.weights)) weights <- y$weights
if(!is.null(weights)) {
weights <- as.matrix(weights)
- if(any(dim(weights) != dim(y))) stop("weights must have same dimensions as y")
+ dimw <- dim(weights)
+ if(dimw[1] != ngenes || dimw[2] != n) stop("weights must have same dimensions as y")
+ if(any(weights <= 0)) stop("weights must be positive")
if(!is.null(array.weights)) {
weights <- .matvec(weights,array.weights)
array.weights <- NULL
}
}
+# Reduce to numeric expression matrix
+ y <- y$exprs
+
# Divide out array.weights
if(!is.null(array.weights)) {
sw <- sqrt(array.weights)
@@ -407,20 +430,19 @@ mroast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.st
df.prior <- sv$df.prior
}
- pv <- adjpv <- active <- array(0,c(nsets,3),dimnames=list(names(index),c("Down","Up","Mixed")))
+ pv <- adjpv <- active <- array(0,c(nsets,4),dimnames=list(names(index),c("Down","Up","UpOrDown","Mixed")))
NGenes <- rep(0,nsets)
if(nsets<1) return(pv)
for(i in 1:nsets) {
- out <- roast(y=y,index=index[[i]],design=design,contrast=contrast,set.statistic=set.statistic,gene.weights=gene.weights,weights=weights,var.prior=var.prior,df.prior=df.prior,nrot=nrot)
+ out <- roast(y=y,index=index[[i]],design=design,contrast=contrast,set.statistic=set.statistic,gene.weights=gene.weights,weights=weights,var.prior=var.prior,df.prior=df.prior,nrot=nrot,approx.zscore=approx.zscore,...)
pv[i,] <- out$p.value$P.Value
active[i,] <- out$p.value$Active.Prop
NGenes[i] <- out$ngenes.in.set
}
# Use mid-p-values or ordinary p-values?
- pv2 <- pv
- if(midp) pv2 <- pv2-1/2/(nrot+1)
-
+# pv2 <- pv
+# if(midp) pv2 <- pv2-1/2/(nrot+1)
# adjpv[,"Down"] <- p.adjust(pv2[,"Down"], method=adjust.method)
# adjpv[,"Up"] <- p.adjust(pv2[,"Up"], method=adjust.method)
# adjpv[,"Mixed"] <- p.adjust(pv2[,"Mixed"], method=adjust.method)
@@ -429,17 +451,22 @@ mroast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.st
# New-style output
Up <- pv[,"Up"] < pv[,"Down"]
Direction <- rep.int("Down",nsets); Direction[Up] <- "Up"
- TwoSidedP <- pv[,"Down"]; TwoSidedP[Up] <- pv[Up,"Up"]; TwoSidedP <- 2*TwoSidedP
- TwoSidedP2 <- pv2[,"Down"]; TwoSidedP2[Up] <- pv2[Up,"Up"]; TwoSidedP2 <- 2*TwoSidedP2
+ TwoSidedP2 <- pv[,"UpOrDown"]
+ MixedP2 <- pv[,"Mixed"]
+ if(midp) {
+ TwoSidedP2 <- TwoSidedP2 - 1/2/(nrot+1)
+ MixedP2 <- MixedP2 - 1/2/(nrot+1)
+ }
+
tab <- data.frame(
NGenes=NGenes,
PropDown=active[,"Down"],
PropUp=active[,"Up"],
Direction=Direction,
- PValue=TwoSidedP,
+ PValue=pv[,"UpOrDown"],
FDR=p.adjust(TwoSidedP2,method="BH"),
PValue.Mixed=pv[,"Mixed"],
- FDR.Mixed=p.adjust(pv2[,"Mixed"],method="BH"),
+ FDR.Mixed=p.adjust(MixedP2,method="BH"),
row.names=names(index),
stringsAsFactors=FALSE
)
@@ -448,9 +475,9 @@ mroast.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),set.st
sort <- match.arg(sort,c("directional","mixed","none"))
if(sort=="none") return(tab)
if(sort=="directional")
- o <- order(tab$PValue)
+ o <- order(tab$PValue,tab$PValue.Mixed,-tab$NGenes)
else
- o <- order(tab$PValue.Mixed)
+ o <- order(tab$PValue.Mixed,tab$PValue,-tab$NGenes)
tab[o,,drop=FALSE]
}
diff --git a/R/geneset-romer.R b/R/geneset-romer.R
index dbc1ef0..19f178d 100644
--- a/R/geneset-romer.R
+++ b/R/geneset-romer.R
@@ -116,29 +116,29 @@ romer <- function(index,y,design,contrast=ncol(design),array.weights=NULL,block=
mh<-.meanHalf(s.r[index[[i]]],m[i])
s.rank.up[i] <-mh[2]
s.rank.down[i]<-mh[1]
- s.rank.mixed[i]<-.meanHalf(s.abs.r[index[[i]]],m[i])[2]
+ s.rank.mixed[i]<-.meanHalf(s.abs.r[index[[i]]],m[i])[2]
}
# Estimate hyper-parameters p0
p.obs <- 2*pt(abs(modt),df=d0+d,lower.tail=FALSE)
- p0.obs <- 1-convest(p.obs) # proportion of DE probes
-
+ p0.obs <- 1-propTrueNull(p.obs) # proportion of DE probes
+
# Estimate hyper-paremeter v0
covmat <- chol2inv(qr$qr, size = qr$rank)
stdev.unscaled <- rep(sqrt(covmat[qr$rank,qr$rank]),ngenes)
-
- proportion<-p0.obs
+
+ proportion <- p0.obs
stdev.coef.lim <- c(0.1, 4)
-
+
# get v0
df.total <- rep(d,ngenes) + sv$df.prior
var.prior.lim <- stdev.coef.lim^2/sv$var.prior
- var.prior <- tmixture.vector(modt, stdev.unscaled, df.total,proportion, var.prior.lim)
+ var.prior <- tmixture.vector(modt, stdev.unscaled, df.total, proportion, var.prior.lim)
if (any(is.na(var.prior))) {
var.prior[is.na(var.prior)] <- 1/sv$var.prior
warning("Estimation of var.prior failed - set to default value")
}
-
+
# Estimate posterior probability of DE
r <- rep(1, ngenes) %o% var.prior
r <- (stdev.unscaled^2 + r)/stdev.unscaled^2
@@ -180,7 +180,7 @@ romer <- function(index,y,design,contrast=ncol(design),array.weights=NULL,block=
s.rank.up.2 <-mh.2[2]
s.rank.down.2 <-mh.2[1]
- s.rank.mixed.2 <-.meanHalf(s.abs.r.2[index[[j]]],m[j])[2]
+ s.rank.mixed.2 <-.meanHalf(s.abs.r.2[index[[j]]],m[j])[2]
if(s.rank.mixed.2>=s.rank.mixed[j]) p.value[j,1]<-p.value[j,1]+1
if(s.rank.up.2>=s.rank.up[j]) p.value[j,2]<-p.value[j,2]+1
@@ -256,7 +256,7 @@ romer <- function(index,y,design,contrast=ncol(design),array.weights=NULL,block=
R <- chol(cormatrix)
y <- t(backsolve(R, t(y), transpose = TRUE))
design <- backsolve(R, design, transpose = TRUE)
- }
+ }
# Reform design matrix so that contrast of interest is last column
qr <- qr(contrast)
@@ -288,7 +288,7 @@ romer <- function(index,y,design,contrast=ncol(design),array.weights=NULL,block=
# Estimate hyper-parameter p0
p.obs <- 2*pt(abs(modt),df=d0+d,lower.tail=FALSE)
- p0.obs <- 1-convest(p.obs) # proportion of DE probes
+ p0.obs <- 1-propTrueNull(p.obs) # proportion of DE probes
# Estimate hyper-paremeter v0
covmat <- chol2inv(qr$qr, size = qr$rank)
@@ -375,7 +375,7 @@ romer <- function(index,y,design,contrast=ncol(design),array.weights=NULL,block=
cbind(NGenes=set.sizes,p.value)
}
-topRomer<-function(x,n=10,alternative="up")
+topRomer <- function(x,n=10,alternative="up")
# extracts a number of top gene sets results from the romer output.
# Gordon Smyth and Yifang Hu.
# 22 Mar 2010. Last modified 5 Sep 2011.
diff --git a/R/lmfit.R b/R/lmfit.R
index aecb915..4ba39ab 100644
--- a/R/lmfit.R
+++ b/R/lmfit.R
@@ -387,10 +387,10 @@ residuals.MArrayLM <- function(object,y,...)
}
getEAWP <- function(object)
-# Given any microarray data object, extract basic information needed for
-# linear modelling.
+# Given any microarray data object, extract basic information needed
+# for linear modelling.
# Gordon Smyth
-# 9 March 2008. Last modified 26 June 2013.
+# 9 March 2008. Last modified 26 June 2013.
{
y <- list()
diff --git a/R/loessFit.R b/R/loessFit.R
index 788c8be..cb50521 100644
--- a/R/loessFit.R
+++ b/R/loessFit.R
@@ -1,66 +1,112 @@
# LOESS FUNCTIONS
-loessFit <- function(y, x, weights=NULL, span=0.3, iterations=4, min.weight=1e-5, max.weight=1e5, equal.weights.as.null=TRUE)
+loessFit <- function(y, x, weights=NULL, span=0.3, iterations=4L, min.weight=1e-5, max.weight=1e5, equal.weights.as.null=TRUE, method="weightedLowess")
# Fast lowess fit for univariate x and y allowing for weights
-# Uses lowess() if weights=NULL and locfit.raw() otherwise
+# Uses lowess() if weights=NULL and weightedLowess(), locfit.raw() or loess() otherwise
# Gordon Smyth
-# 28 June 2003. Last revised 2 Oct 2013.
+# 28 June 2003. Last revised 17 Apr 2014.
{
+# Check x and y
n <- length(y)
+ if(length(x) != n) stop("y and x have different lengths")
out <- list(fitted=rep(NA,n),residuals=rep(NA,n))
+
obs <- is.finite(y) & is.finite(x)
xobs <- x[obs]
yobs <- y[obs]
nobs <- length(yobs)
+
+# If no good obs, exit straight away
if(nobs==0) return(out)
- if(is.null(weights)) {
- o <- order(xobs)
- 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 {
+
+# Check span
+ if(span < 1/nobs) {
+ out$fitted[obs] <- y[obs]
+ out$residuals[obs] <- 0
+ return(out)
+ }
+
+# Check min.weight
+ if(min.weight<0) min.weight <- 0
+
+# Check weights
+ if(!is.null(weights)) {
+ if(length(weights) != n) stop("y and weights have different lengths")
wobs <- weights[obs]
wobs[is.na(wobs)] <- 0
wobs <- pmax(wobs,min.weight)
wobs <- pmin(wobs,max.weight)
-# Test whether weights are equal
+# If weights all equal, treat as NULL
if(equal.weights.as.null) {
r <- range(wobs)
- if(r[2]-r[1] < 1e-15) return(Recall(y=y,x=x,span=span,iterations=iterations))
+ if(r[2]-r[1] < 1e-15) weights <- NULL
}
-# Count number of observations with positive weights
- if(min.weight>0)
- nwobs <- nobs
- else
- nwobs <- sum(wobs>0)
- if(nwobs < 4+1/span) {
- if(nwobs>1) {
- fit <- lm.wfit(cbind(1,xobs),yobs,wobs)
- } else {
- fit <- list()
- fit$fitted <- rep(sum(wobs*yobs)/sum(wobs),nobs)
- fit$residuals <- yobs-fit$fitted
- }
+ }
+
+# If no weights, so use classic lowess algorithm
+ if(is.null(weights)) {
+ o <- order(xobs)
+ lo <- lowess(x=xobs,y=yobs,f=span,iter=iterations-1L)
+ out$fitted[obs][o] <- lo$y
+ out$residuals[obs] <- yobs-out$fitted[obs]
+ return(out)
+ }
+
+# Count number of observations with positive weights (must always be positive)
+ if(min.weight>0)
+ nwobs <- nobs
+ else
+ nwobs <- sum(wobs>0)
+
+# Check whether too few obs to estimate lowess curve
+ if(nwobs < 4+1/span) {
+ if(nwobs==1L) {
+ out$fitted[obs] <- yobs[wobs>0]
+ out$residuals[obs] <- yobs-out$fitted[obs]
} else {
+ fit <- lm.wfit(cbind(1,xobs),yobs,wobs)
+ out$fitted[obs] <- fit$fitted
+ out$residuals[obs] <- fit$residuals
+ }
+ return(out)
+ }
+
+# Need to compute lowess with unequal weights
+ method <- match.arg(method, c("weightedLowess","locfit","loess"))
+ switch(method,
+ "weightedLowess" = {
+ fit <- weightedLowess(x=xobs,y=yobs,weights=wobs,span=span,iterations=iterations,npts=200)
+ out$fitted[obs] <- fit$fitted
+ out$residuals[obs] <- fit$residuals
+ },
+ "locfit" = {
# 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)
+ if(inherits(loadresult,"error")) stop("locfit package 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)
+ 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] <- fitted(fit)
+ out$residuals[obs] <- res
+ },
+ "loess" = {
+# Suppress warning "k-d tree limited by memory"
+ oldopt <- options(warn=-1)
+ on.exit(options(oldopt))
+ bin <- 0.01
+ fit <- loess(yobs~xobs,weights=wobs,span=span,degree=1,parametric=FALSE,normalize=FALSE,statistics="approximate",surface="interpolate",cell=bin/span,iterations=iterations,trace.hat="approximate")
+ out$fitted[obs] <- fit$fitted
+ out$residuals[obs] <- fit$residuals
}
- out$fitted[obs] <- fitted(fit)
- out$residuals[obs] <- res
- }
+ )
out
}
diff --git a/R/plot.R b/R/plot.R
new file mode 100644
index 0000000..252a43d
--- /dev/null
+++ b/R/plot.R
@@ -0,0 +1,70 @@
+## PLOT.R
+## plot() produces MA plots on expression objects
+
+plot.RGList <- function(x, y, array=1, xlab="A", ylab="M", main=colnames(x)[array], xlim=NULL, ylim=NULL, status, values, pch, col, cex, legend=TRUE, zero.weights=FALSE, ...)
+# MA-plot with color coding for controls
+# Gordon Smyth
+# Created 21 March 2014
+{
+ MA <- MA.RG(x[,array])
+ plot.MAList(MA,array=1,xlab=xlab,ylab=ylab,main=main,xlim=xlim,ylim=ylim,status=status,values=values,pch=pch,col=col,cex=cex,legend=legend,zero.weights=zero.weights,...)
+}
+
+plot.MAList <- function(x, y, array=1, xlab="A", ylab="M", main=colnames(x)[array], xlim=NULL, ylim=NULL, status, values, pch, col, cex, legend=TRUE, zero.weights=FALSE, ...)
+# MA-plot with color coding for controls
+# Gordon Smyth
+# Created 21 March 2014
+{
+ MA <- x
+ 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=xlab,ylab=ylab,main=main,xlim=xlim,ylim=ylim,status=status,values=values,pch=pch,col=col,cex=cex,legend=legend, ...)
+}
+
+plot.MArrayLM <- function(x, y, coef=ncol(x), xlab="AveExpr", ylab="logFC", main=colnames(x)[coef], xlim=NULL, ylim=NULL, status, values, pch, col, cex, legend=TRUE, zero.weights=FALSE, ...)
+# MA-plot with color coding for controls
+# Gordon Smyth
+# Created 21 March 2014
+{
+ fit <- x
+ if(is.null(fit$Amean)) stop("MA-plot not possible because Amean component is absent.")
+ x <- fit$Amean
+ y <- as.matrix(fit$coef)[,coef]
+ if(is.null(fit$weights)) w <- NULL else w <- as.matrix(fit$weights)[,coef]
+ if(missing(status)) status <- fit$genes$Status
+ if(!is.null(w) && !zero.weights) {
+ i <- is.na(w) | (w <= 0)
+ y[i] <- NA
+ }
+ .plotMAxy(x,y,xlab=xlab,ylab=ylab,main=main,xlim=xlim,ylim=ylim,status=status,values=values,pch=pch,col=col,cex=cex,legend=legend, ...)
+}
+
+plot.EList <- function(x, y, array=1, xlab="A", ylab="M", main=colnames(x)[array], xlim=NULL, ylim=NULL, status, values, pch, col, cex, legend=TRUE, zero.weights=FALSE, ...)
+# MA-plot with color coding for controls
+# Gordon Smyth
+# Created 21 March 2014
+{
+ E <- x
+ E$E <- as.matrix(E$E)
+ narrays <- ncol(E$E)
+ if(narrays < 2) stop("Need at least two arrays")
+ if(narrays > 5)
+ x <- apply(E$E,1,median,na.rm=TRUE)
+ else
+ x <- rowMeans(E$E,na.rm=TRUE)
+ y <- E$E[,array]-x
+ if(is.null(E$weights)) w <- NULL else w <- as.matrix(E$weights)[,array]
+ if(missing(status)) status <- E$genes$Status
+
+ if(!is.null(w) && !zero.weights) {
+ i <- is.na(w) | (w <= 0)
+ y[i] <- NA
+ }
+ .plotMAxy(x,y,xlab=xlab,ylab=ylab,main=main,xlim=xlim,ylim=ylim,status=status,values=values,pch=pch,col=col,cex=cex,legend=legend, ...)
+}
diff --git a/R/plotMDS.R b/R/plotMDS.R
index ccdfaee..0999a01 100644
--- a/R/plotMDS.R
+++ b/R/plotMDS.R
@@ -10,12 +10,18 @@ setMethod("show","MDS",function(object) {
plotMDS <- function(x,...) UseMethod("plotMDS")
-plotMDS.MDS <- function(x,labels=colnames(x$distance.matrix),col=NULL,cex=1,dim.plot=x$dim.plot,xlab=paste("Dimension",dim.plot[1]),ylab=paste("Dimension",dim.plot[2]),...)
+plotMDS.MDS <- function(x,labels=NULL,pch=NULL,col=NULL,cex=1,dim.plot=x$dim.plot,xlab=paste("Dimension",dim.plot[1]),ylab=paste("Dimension",dim.plot[2]),...)
# Method for MDS objects
# Create a new plot using MDS coordinates or distances previously created
-# Gordon Smyth
-# 21 May 2011. Last modified 6 Sep 2012.
+# Gordon Smyth and Yifang Hu
+# 21 May 2011. Last modified 2 June 2014
{
+# Check labels
+ if(is.null(labels) & is.null(pch)) {
+ labels <- colnames(x$distance.matrix)
+ if(is.null(labels)) labels <- 1:length(x$x)
+ }
+
# Are new dimensions requested?
if(!all(dim.plot==x$dim.plot)) {
ndim <- max(dim.plot)
@@ -24,24 +30,28 @@ plotMDS.MDS <- function(x,labels=colnames(x$distance.matrix),col=NULL,cex=1,dim.
x$y <- x$cmdscale.out[,dim.plot[2]]
}
-# Estimate width of labels in plot coordinates.
-# Estimate will be ok for default plot width, but maybe too small for smaller plots.
- if(is.null(labels)) labels <- 1:length(x$x)
- labels <- as.character(labels)
- StringRadius <- 0.01*cex*nchar(labels)
- left.x <- x$x-StringRadius
- right.x <- x$x+StringRadius
+# Make the plot
+ if(is.null(labels)){
+# Plot symbols instead of text
+ plot(x$x, x$y, pch = pch, xlab = xlab, ylab = ylab, col = col, cex = cex, ...)
+ } else {
+# Plot text. Need to estimate width of labels in plot coordinates.
+# Estimate will be ok for default plot width, but maybe too small for smaller plots.
+ labels <- as.character(labels)
+ StringRadius <- 0.01*cex*nchar(labels)
+ left.x <- x$x-StringRadius
+ right.x <- x$x+StringRadius
+ plot(c(left.x, right.x), c(x$y, x$y), type = "n", xlab = xlab, ylab = ylab, ...)
+ text(x$x, x$y, labels = labels, col = col, cex = cex)
+ }
-# Redo plot
- plot(c(left.x,right.x),c(x$y,x$y),type="n",xlab=xlab,ylab=ylab,...)
- text(x$x,x$y,labels=labels,col=col,cex=cex)
return(invisible(x))
}
-plotMDS.default <- function(x,top=500,labels=colnames(x),col=NULL,cex=1,dim.plot=c(1,2),ndim=max(dim.plot),gene.selection="pairwise",xlab=paste("Dimension",dim.plot[1]),ylab=paste("Dimension",dim.plot[2]),...)
+plotMDS.default <- function(x,top=500,labels=NULL,pch=NULL,col=NULL,cex=1,dim.plot=c(1,2),ndim=max(dim.plot),gene.selection="pairwise",xlab=paste("Dimension",dim.plot[1]),ylab=paste("Dimension",dim.plot[2]),...)
# Multi-dimensional scaling with top-distance
# Di Wu and Gordon Smyth
-# 19 March 2009. Last modified 29 May 2013.
+# 19 March 2009. Last modified 22 May 2014
{
# Check x
x <- as.matrix(x)
@@ -56,13 +66,16 @@ plotMDS.default <- function(x,top=500,labels=colnames(x),col=NULL,cex=1,dim.plot
# Check top
top <- min(top,nprobes)
-# Check labels
- if(is.null(labels)) labels <- 1:nsamples
- labels <- as.character(labels)
+# Check labels and pch
+ if(is.null(pch) & is.null(labels)) {
+ labels <- colnames(x)
+ if(is.null(labels)) labels <- 1:nsamples
+ }
+ if(!is.null(labels)) labels <- as.character(labels)
# Check dim
if(ndim < 2) stop("Need at least two dim.plot")
- if(nsamples < ndim) stop("Two few samples")
+ if(nsamples < ndim) stop("Too few samples")
if(nprobes < ndim) stop("Too few rows")
# Check gene.selection
@@ -92,5 +105,5 @@ plotMDS.default <- function(x,top=500,labels=colnames(x),col=NULL,cex=1,dim.plot
mds <- new("MDS",list(dim.plot=dim.plot,distance.matrix=dd,cmdscale.out=a1,top=top,gene.selection=gene.selection))
mds$x <- a1[,dim.plot[1]]
mds$y <- a1[,dim.plot[2]]
- plotMDS(mds,labels=labels,col=col,cex=cex,xlab=xlab,ylab=ylab,...)
+ plotMDS(mds,labels=labels,pch=pch,col=col,cex=cex,xlab=xlab,ylab=ylab,...)
}
diff --git a/R/plots-ma.R b/R/plots-ma.R
index 0ce4aec..954456c 100755
--- a/R/plots-ma.R
+++ b/R/plots-ma.R
@@ -28,15 +28,15 @@ plotMA.MAList <- function(MA, array=1, xlab="A", ylab="M", main=colnames(MA)[arr
.plotMAxy(x,y,xlab=xlab,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, ...)
+plotMA.MArrayLM <- function(MA, coef=ncol(MA), xlab="AveExpr", ylab="logFC", main=colnames(MA)[coef], 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.
+# Last modified 21 March 2014.
{
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]
+ y <- as.matrix(MA$coef)[,coef]
+ if(is.null(MA$weights)) w <- NULL else w <- as.matrix(MA$weights)[,coef]
if(missing(status)) status <- MA$genes$Status
if(!is.null(w) && !zero.weights) {
i <- is.na(w) | (w <= 0)
@@ -179,9 +179,25 @@ plotMA.default <- function(MA, array=1, xlab="A", ylab="M", main=colnames(MA)[ar
plotMA3by2 <- function(MA, prefix="MA", path=NULL, main=colnames(MA), zero.weights=FALSE, common.lim=TRUE, device="png", ...)
# Make files of MA-plots, six to a page
-# Gordon Smyth 27 May 2004. Last modified 9 June 2007.
+# Gordon Smyth 27 May 2004. Last modified 12 Feb 2014.
{
if(is(MA,"RGList")) MA <- MA.RG(MA)
+ if(is(MA,"EListRaw")) {
+ A <- rowMeans(log2(MA$E))
+ MA$A <- array(A,dim(MA))
+ MA$M <- log2(MA$E)-MA$A
+ MA <- new("MAList",unclass(MA))
+ }
+ if(is(MA,"EList")) {
+ A <- rowMeans(MA$E)
+ MA$A <- array(A,dim(MA))
+ MA$M <- MA$E-MA$A
+ MA <- new("MAList",unclass(MA))
+ }
+ if(is.matrix(MA)) {
+ A <- rowMeans(MA)
+ MA <- new("MAList",list(M=MA-A, A=array(A,dim(MA))))
+ }
if(is.null(path)) path <- "."
prefix <- file.path(path,prefix)
narrays <- ncol(MA)
@@ -236,7 +252,7 @@ plotPrintTipLoess <- function(object,layout,array=1,span=0.4,...) {
}
mdplot <- function(x,...)
-# Mean-different plot
+# Mean-difference plot
# Gordon Smyth
# 16 March 2005
{
diff --git a/R/read-ilmn.R b/R/read-ilmn.R
index 952f317..8530aba 100644
--- a/R/read-ilmn.R
+++ b/R/read-ilmn.R
@@ -3,7 +3,7 @@
read.ilmn <- function(files=NULL, ctrlfiles=NULL, path=NULL, ctrlpath=NULL, probeid="Probe", annotation=c("TargetID", "SYMBOL"), expr="AVG_Signal", other.columns="Detection",sep="\t", quote="\"", verbose=TRUE, ...)
# Read one or more files of Illumina BeadStudio output
# Wei Shi and Gordon Smyth.
-# Created 15 July 2009. Last modified 24 January 2014.
+# Created 15 July 2009. Last modified 27 November 2013.
{
if(!is.null(files)){
f <- unique(files)
@@ -63,10 +63,10 @@ read.ilmn.targets <- function(targets, ...)
.read.oneilmnfile <- function(fname, probeid, annotation, expr, other.columns, sep, quote, verbose, ...)
# Read a single file of Illumina BeadStudio output
-# Wei Shi and Gordon Smyth
-# Created 15 July 2009. Last modified 5 Jan 2011.
+# Wei Shi and Gordon Smyth
+# Created 15 July 2009. Last modified 16 June 2014.
{
- h <- readGenericHeader(fname,columns=expr)
+ h <- readGenericHeader(fname,columns=expr,sep=sep)
skip <- h$NHeaderRecords
header <- h$ColumnNames
@@ -101,7 +101,7 @@ read.ilmn.targets <- function(targets, ...)
# Add probe annotation
if(length(anncol)) {
elist$genes <- x[,anncol,drop=FALSE]
- row.names(elist$genes) <- pids
+ if(!any(duplicated(pids))) row.names(elist$genes) <- pids
}
# elist$targets <- data.frame(SampleNames=snames, stringsAsFactors=FALSE)
diff --git a/R/read.idat.R b/R/read.idat.R
new file mode 100644
index 0000000..9145b21
--- /dev/null
+++ b/R/read.idat.R
@@ -0,0 +1,41 @@
+read.idat <- function(idatfiles, bgxfile, dateinfo=FALSE)
+# Function to read idat data from gene expression BeadArrays
+# Matt Ritchie, 30 September 2013
+{
+ require(illuminaio)
+ cat("Reading manifest file", bgxfile, "... ")
+ bgx <- readBGX(bgxfile)
+ cat("Done\n")
+ nregprobes <-nrow(bgx$probes)
+ nctrlprobes <-nrow(bgx$control)
+ nprobes <- nregprobes+nctrlprobes
+ nsamples <- length(idatfiles)
+ elist <- new("EListRaw")
+ elist$source <- "illumina"
+ elist$genes <- rbind(bgx$probes[,c("Probe_Id","Array_Address_Id")], bgx$controls[,c("Probe_Id","Array_Address_Id")])
+ elist$genes$Status <- "regular"
+ elist$genes$Status[(nregprobes+1):nprobes] <- bgx$controls[,"Reporter_Group_Name"]
+ elist$targets <- data.frame("IDATfile"=idatfiles, "DecodeInfo"=rep(NA, nsamples), "ScanInfo"=rep(NA, nsamples))
+ tmp <- matrix(NA, nprobes, nsamples)
+ colnames(tmp) <- idatfiles
+ rownames(tmp) <- elist$genes[,"Array_Address_Id"]
+ elist$E <- elist$other$STDEV <- elist$other$NumBeads <- tmp
+ rm(tmp)
+ for(i in 1:nsamples) {
+ cat("\t", idatfiles[i], "... ")
+ tmp <- readIDAT(idatfiles[i])
+ cat("Done\n")
+ ind <- match(elist$genes[,"Array_Address_Id"], tmp$Quants[,'IllumicodeBinData'])
+ if(sum(is.na(ind))>0)
+ stop("Can't match ids in manifest with those in idat file", idatfiles[i], "- please check that you have the right files\n")
+ elist$E[,i] <- round(tmp$Quants[ind,'MeanBinData'],1) # intensity data
+ elist$other$STDEV[,i] <- tmp$Quants[ind,'DevBinData'] # Bead STDEV
+ elist$other$NumBeads[,i] <- tmp$Quants[ind,'NumGoodBeadsBinData'] # NumBeads
+ if(dateinfo) {
+ elist$targets[i,"DecodeInfo"] = paste(tmp$RunInfo[1,], sep="", collapse=" ")
+ elist$targets[i,"ScanInfo"] = paste(tmp$RunInfo[2,], sep="", collapse=" ")
+ }
+ }
+ cat("Finished reading data.\n")
+ return(elist)
+}
diff --git a/R/removeBatchEffect.R b/R/removeBatchEffect.R
index 9ed1b37..6ed2466 100644
--- a/R/removeBatchEffect.R
+++ b/R/removeBatchEffect.R
@@ -3,13 +3,12 @@
# A refinement would be to empirical Bayes shrink
# the batch effects before subtracting them.
-removeBatchEffect <- function(x,batch=NULL,batch2=NULL,covariates=NULL,design=matrix(1,ncol(x),1))
+removeBatchEffect <- function(x,batch=NULL,batch2=NULL,covariates=NULL,design=matrix(1,ncol(x),1),...)
# Remove batch effects from matrix of expression data
# Gordon Smyth and Carolyn de Graaf
-# Created 1 Aug 2008. Last revised 4 Feb 2012.
+# Created 1 Aug 2008. Last revised 1 June 2014.
{
- x <- as.matrix(x)
- if(is.null(batch) && is.null(batch2) && is.null(covariates)) return(x)
+ if(is.null(batch) && is.null(batch2) && is.null(covariates)) return(as.matrix(x))
if(!is.null(batch)) {
batch <- as.factor(batch)
contrasts(batch) <- contr.sum(levels(batch))
@@ -21,8 +20,9 @@ removeBatchEffect <- function(x,batch=NULL,batch2=NULL,covariates=NULL,design=ma
batch2 <- model.matrix(~batch2)[,-1,drop=FALSE]
}
if(!is.null(covariates)) covariates <- as.matrix(covariates)
- X <- cbind(batch,batch2,covariates)
- X <- qr.resid(qr(design),X)
- t(qr.resid(qr(X),t(x)))
+ X.batch <- cbind(batch,batch2,covariates)
+ fit <- lmFit(x,cbind(design,X.batch),...)
+ beta <- fit$coefficients[,-(1:ncol(design)),drop=FALSE]
+ beta[is.na(beta)] <- 0
+ as.matrix(x) - beta %*% t(X.batch)
}
-
diff --git a/R/subsetting.R b/R/subsetting.R
index 7a6bb3c..8495612 100755
--- a/R/subsetting.R
+++ b/R/subsetting.R
@@ -102,12 +102,12 @@ assign("[.MArrayLM",
function(object, i, j)
# Subsetting for MArrayLM objects
# Gordon Smyth
-# 26 April 2005. Last modified 11 Dec 2013.
+# 26 April 2005. Last modified 14 May 2014.
{
# Recognized components
IJ <- c("coefficients","stdev.unscaled","t","p.value","lods","weights")
IX <- "genes"
- JX <- c("targets","design")
+ JX <- character(0)
I <- c("Amean","sigma","df.residual","df.prior","df.total","s2.post","F","F.p.value")
JJ <- "cov.coefficients"
XJ <- "contrasts"
@@ -117,7 +117,7 @@ function(object, i, j)
# so that output is equivalent to that from contrasts.fit()
if(!missing(j) && is.null(object$contrasts) && !is.null(object$coefficients)) {
object$contrasts <- diag(ncol(object$coefficients))
- cn <- colnames(object$contrasts)
+ cn <- colnames(object$coefficients)
dimnames(object$contrasts) <- list(Coefficient=cn,Contrast=cn)
}
diff --git a/R/toptable.R b/R/toptable.R
index 5998e3d..b87e637 100644
--- a/R/toptable.R
+++ b/R/toptable.R
@@ -3,17 +3,31 @@
topTable <- function(fit,coef=NULL,number=10,genelist=fit$genes,adjust.method="BH",sort.by="B",resort.by=NULL,p.value=1,lfc=0,confint=FALSE)
# Summary table of top genes, object-orientated version
# Gordon Smyth
-# 4 August 2003. Last modified 7 Dec 2013.
+# 4 August 2003. Last modified 27 February 2014.
{
# Check fit
if(!is(fit,"MArrayLM")) stop("fit must be an MArrayLM object")
if(is.null(fit$coefficients)) stop("coefficients not found in fit object")
if(confint && is.null(fit$stdev.unscaled)) stop("stdev.unscaled not found in fit object")
- if(is.null(coef)) coef <- 1:ncol(fit)
+ if(is.null(coef)) {
+ if(is.null(fit$treat.lfc)) {
+ coef <- 1:ncol(fit)
+ cn <- colnames(fit)
+ if(!is.null(cn)) {
+ i <- which(cn == "(Intercept)")
+ if(length(i)) {
+ coef <- coef[-i]
+ message("Removing intercept from test coefficients")
+ }
+ }
+ } else
+ coef <- ncol(fit)
+ }
if(length(coef)>1) {
+ if(!is.null(fit$treat.lfc)) stop("Treat p-values can only be displayed for single coefficients")
coef <- unique(coef)
- if(length(fit$coef[1,coef]) < ncol(fit)) fit <- eBayes(fit[,coef])
+ if(length(fit$coef[1,coef]) < ncol(fit)) fit <- fit[,coef]
if(sort.by=="B") sort.by <- "F"
return(topTableF(fit,number=number,genelist=genelist,adjust.method=adjust.method,sort.by=sort.by,p.value=p.value,lfc=lfc))
}
@@ -37,7 +51,7 @@ 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 4 November 2013.
+# 27 August 2006. Last modified 24 June 2014.
{
# Check fit
if(is.null(fit$coefficients)) stop("Coefficients not found in fit")
@@ -111,7 +125,7 @@ topTableF <- function(fit,number=10,genelist=fit$genes,adjust.method="BH",sort.b
tab <- data.frame(M[o,,drop=FALSE])
else
tab <- data.frame(genelist[o,,drop=FALSE],M[o,,drop=FALSE])
- tab$AveExpr=fit$Amean[o]
+ tab$AveExpr <- 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
@@ -236,9 +250,10 @@ toptable <- function(fit,coef=1,number=10,genelist=NULL,A=NULL,eb=NULL,adjust.me
tab <- data.frame(genelist[top,,drop=FALSE],logFC=M[top],stringsAsFactors=FALSE)
}
if(confint) {
- margin.error <- sqrt(eb$s2.post[top])*fit$stdev.unscaled[top]*qnorm(0.975)
- tab$CI.025 <- M[top]-margin.error
- tab$CI.975 <- M[top]+margin.error
+ if(is.numeric(confint)) alpha <- (1+confint[1])/2 else alpha <- 0.975
+ margin.error <- sqrt(eb$s2.post[top])*fit$stdev.unscaled[top,coef]*qt(alpha,df=eb$df.total[top])
+ tab$CI.L <- M[top]-margin.error
+ tab$CI.R <- M[top]+margin.error
}
if(!is.null(A)) tab$AveExpr <- A[top]
tab <- data.frame(tab,t=tstat[top],P.Value=P.Value[top],adj.P.Val=adj.P.Value[top])
diff --git a/R/treat.R b/R/treat.R
index 8ffd95a..610cbbc 100644
--- a/R/treat.R
+++ b/R/treat.R
@@ -1,14 +1,15 @@
### treat.R
treat <- function(fit, lfc=0, trend=FALSE, robust=FALSE, winsor.tail.p=c(0.05,0.1))
-# Moderated t-statistics with threshold
-# Davis McCarthy, Gordon Smyth
-# 25 July 2008. Last revised 7 April 2013.
+# Moderated t-statistics with threshold
+# Davis McCarthy, Gordon Smyth
+# 25 July 2008. Last revised 27 February 2014.
{
# Check fit
if(!is(fit,"MArrayLM")) stop("fit must be an MArrayLM object")
if(is.null(fit$coefficients)) stop("coefficients not found in fit object")
if(is.null(fit$stdev.unscaled)) stop("stdev.unscaled not found in fit object")
+ fit$lods <- NULL
coefficients <- as.matrix(fit$coefficients)
stdev.unscaled <- as.matrix(fit$stdev.unscaled)
@@ -47,13 +48,16 @@ treat <- function(fit, lfc=0, trend=FALSE, robust=FALSE, winsor.tail.p=c(0.05,0.
fc.down <- (coefficients < -lfc)
fit$t[fc.up] <- tstat.right[fc.up]
fit$t[fc.down] <- -tstat.right[fc.down]
+ fit$treat.lfc <- lfc
fit
}
-topTreat <- function(fit,coef=1,number=10,genelist=fit$genes,adjust.method="BH",sort.by="p",resort.by=NULL,p.value=1)
+topTreat <- function(fit,coef=1,sort.by="p",resort.by=NULL,...)
# Summary table of top genes by treat
# Gordon Smyth
-# 15 June 2009. Last modified 20 April 2013.
+# 15 June 2009. Last modified 27 February 2014.
+
+# NOTE: This function may become deprecated soon as topTable() takes over
{
# Check coef is length 1
if(length(coef)>1) {
@@ -65,6 +69,6 @@ topTreat <- function(fit,coef=1,number=10,genelist=fit$genes,adjust.method="BH",
if(sort.by=="B") stop("Trying to sort.by B, but treat doesn't produce a B-statistic")
if(!is.null(resort.by)) if(resort.by=="B") stop("Trying to resort.by B, but treat doesn't produce a B-statistic")
- topTable(fit=fit,coef=coef,number=number,genelist=genelist,adjust.method=adjust.method,sort.by=sort.by,resort.by=sort.by,p.value=p.value)
+ topTable(fit=fit,coef=coef,sort.by=sort.by,...)
}
diff --git a/R/tricubeMovingAverage.R b/R/tricubeMovingAverage.R
new file mode 100644
index 0000000..0ab341b
--- /dev/null
+++ b/R/tricubeMovingAverage.R
@@ -0,0 +1,18 @@
+tricubeMovingAverage <- function(x,span=0.5,full.length=TRUE) {
+ n <- length(x)
+ width <- span*n
+# Round width of smoothing window to nearest odd number
+ hwidth <- as.integer(width %/% 2L)
+ if(hwidth <= 0L) return(x)
+ width <- 2L * hwidth + 1L
+ u <- seq(from=-1,to=1,length=width) * width / (width+1)
+ tricube.weights <- (1-abs(u)^3)^3
+ tricube.weights <- tricube.weights/sum(tricube.weights)
+ if(!full.length) return(as.vector(filter(x,tricube.weights),mode="numeric")[(hwidth+1):(n-hwidth)])
+ z <- numeric(hwidth)
+ x <- as.vector(filter(c(z,x,z),tricube.weights),mode="numeric")[(hwidth+1):(n+hwidth)]
+ cw <- cumsum(tricube.weights)
+ x[1:hwidth] <- x[1:hwidth] / cw[(width-hwidth):(width-1)]
+ x[(n-hwidth+1):n] <- x[(n-hwidth+1):n] / cw[(width-1):(width-hwidth)]
+ x
+}
diff --git a/R/utility.R b/R/utility.R
index 012d1a7..2006942 100755
--- a/R/utility.R
+++ b/R/utility.R
@@ -58,7 +58,7 @@ blockDiag <- function(...)
helpMethods <- function(genericFunction) {
# Prompt user for help topics on methods for generic function
# Gordon Smyth
-# 21 April 2003. Last revised 28 Oct 2003.
+# 21 April 2003. Last revised 24 June 2014.
objectclass <- class(genericFunction)
if(objectclass != "standardGeneric") {
@@ -70,7 +70,7 @@ helpMethods <- function(genericFunction) {
}
}
functionname <- genericFunction at generic
- methodnames <- names(getMethods(genericFunction)@methods)
+ methodnames <- names(findMethods(genericFunction))
nmethods <- length(methodnames)
if(nmethods == 0) {
cat("No available methods\n")
@@ -80,7 +80,7 @@ helpMethods <- function(genericFunction) {
for (i in 1:nmethods) cat(i,": ",aliasnames[i],"\n",sep="")
cat("Type number to choose help topic: ")
n <- as.integer(readline())
- if(n > 0 && n <= nmethods)
+ if(!is.na(n) && n > 0 && n <= nmethods)
eval(parse(text=paste("help(\"",aliasnames[n],"\")",sep="")))
else {
cat("No topic chosen\n")
diff --git a/R/weightedLowess.R b/R/weightedLowess.R
new file mode 100644
index 0000000..dfd7ab1
--- /dev/null
+++ b/R/weightedLowess.R
@@ -0,0 +1,49 @@
+weightedLowess <- function(x, y, weights=rep(1, length(y)), delta=NULL, npts=200, span=0.3, iterations=4)
+# This function clusters points by average linkage and fits a lowess curve of degree 1 to the cluster
+# midpoints. Fitted values are computed by linear interpolation of the fitted coefficients (i.e.
+# quadratic interpolation between points. Several iterations of robustification are also performed
+# using the fitted residuals.
+#
+# written by Aaron Lun
+{
+ len<-length(y)
+ x<-as.double(x)
+ y<-as.double(y)
+ if (length(x)!=len || length(weights)!=len) { stop("vectors have unequal lengths"); }
+ weights<-as.double(weights)
+ iterations<-as.integer(iterations+0.5);
+
+ # Choosing an appropriate 'delta' for approximation. We assume that the covariates
+ # have some cluster structure, where clusters are defined by partitioning on the
+ # 'numclusters' largest distances between points. We want 'npts' evenly spaced points
+ # across the cluster range (i.e. the total range minus the partitioned distances).
+ # Each cluster must also have at least one point; so, we compute the spacing (and
+ # thus delta) as the ratio of the cluster range to the remaining number of points
+ # (after each additional cluster beyond the required first one has eaten 1 point).
+ o<-order(x)
+ x<-x[o]
+ if (is.null(delta)) {
+ npts<-as.integer(npts+0.5)
+ if (npts < 1L) { stop("number of points should be a positive integer") }
+ if (npts >= length(x)) {
+ delta<-0
+ } else {
+ dx<-sort(diff(x))
+ cumrange<-cumsum(dx)
+ numclusters<-1:npts-1L
+ delta<-min(cumrange[length(dx)-numclusters]/(npts-numclusters))
+ }
+ }
+ delta<-as.double(delta)
+
+ # Running the smoothing procedure with specified values.
+ out<-.Call("weighted_lowess", x, y[o], weights[o], span, iterations, delta, PACKAGE="limma")
+ names(out)<-c("fitted", "weights")
+ out$fitted[o]<-out$fitted
+ out$residuals<-y-out$fitted
+ out$weights[o]<-out$weights
+ out$delta<-delta
+
+ return(out);
+}
+
diff --git a/R/zscore.R b/R/zscore.R
index 328bc01..c38e28f 100755
--- a/R/zscore.R
+++ b/R/zscore.R
@@ -31,16 +31,27 @@ zscoreGamma <- function(q, shape, rate = 1, scale = 1/rate)
z
}
-zscoreT <- function(x, df)
+zscoreT <- function(x, df, approx=FALSE)
# Z-score equivalents for t distribution deviates
# Gordon Smyth
-# 24 August 2003
+# 24 August 2003. Last modified 3 June 2014.
{
- z <- x
- df <- rep(df,length.out=length(x))
- pos <- x>0
- if(any(pos)) z[pos] <- qnorm(pt(x[pos],df=df[pos],lower.tail=FALSE,log.p=TRUE),lower.tail=FALSE,log.p=TRUE)
- if(any(!pos)) z[!pos] <- qnorm(pt(x[!pos],df=df[!pos],lower.tail=TRUE,log.p=TRUE),lower.tail=TRUE,log.p=TRUE)
+ if(approx) {
+ if(all(df > 10000)) return(x)
+
+# Approximation from Hill (1970)
+ A <- df-0.5
+ B <- 48*A*A
+ W <- A*log1p(x*x/df)
+ z <- (((((-0.4*W-3.3)*W-24)*W-85.5)/(0.8*W*W+100+B)+W+3)/B+1)*sqrt(W)
+ z[x<0] <- -z[x<0]
+ } else {
+ z <- x
+ df <- rep(df,length.out=length(x))
+ pos <- x>0
+ if(any(pos)) z[pos] <- qnorm(pt(x[pos],df=df[pos],lower.tail=FALSE,log.p=TRUE),lower.tail=FALSE,log.p=TRUE)
+ if(any(!pos)) z[!pos] <- qnorm(pt(x[!pos],df=df[!pos],lower.tail=TRUE,log.p=TRUE),lower.tail=TRUE,log.p=TRUE)
+ }
z
}
diff --git a/build/vignette.rds b/build/vignette.rds
index 4a915c2..253b32f 100644
Binary files a/build/vignette.rds and b/build/vignette.rds differ
diff --git a/inst/NEWS.Rd b/inst/NEWS.Rd
index 45eabec..db768f7 100644
--- a/inst/NEWS.Rd
+++ b/inst/NEWS.Rd
@@ -2,6 +2,123 @@
\title{limma News}
\encoding{UTF-8}
+\section{Version 3.20.0}{\itemize{
+
+\item
+New functions diffSplice(), topSplice() and plotSplice() provide functionality to analyse differential splicing using exon-level expression data from either microarrays or RNA-seq.
+
+\item
+New Pasilla case study added to User's Guide, demonstrating differential splicing analysis of RNA-seq data.
+
+\item
+new function weightedLowess() which fits a lowess curve with prior weights.
+Unlike previous implementations of lowess or loess, the weights are used in calculating which neighbouring points to include in each local regression as well as in the local regression itself.
+
+\item
+weightedLoess() now becomes the default method used by loessFit() to fit the loess curve when there are weights.
+The previous locfit and loess() methods are offered as options.
+
+\item
+linear model fit functions lm.series(), mrlm.series() and gls.series() no longer drop the dimensions of the components of the fitted object when there is just coefficient or just one gene.
+Previously this was done inconsistently in some cases but not others.
+Now the matrix components always keep dimensions.
+
+\item
+The functions lmFit(), eBayes() and tmixture.vector() now work even when there is just one gene (one row of data).
+
+\item
+New function subsetListOfArrays(), which is used to simplify the subsetting code for RGList, MAList, EList, EListRaw and MArrayLM objects.
+
+\item
+new function tricubeMovingAverage() for smoothing a time series.
+
+\item
+barcodeplot() has a new option to add enrichment worms to the plot, making use of tricubeMovingAverage().
+
+\item
+New plot() methods for RGList, MAList, EList and MArrayLM class objects. In each case, this produces a similar result to plotMA().
+When using plot() or plotMA() on an MArrayLM object, the column is now specified by the 'coef' argument instead of by 'array'.
+
+\item
+plotMA3by2() now works on single channel data objects as well as on MAList objects.
+
+\item
+New function read.idat() to read files from Illumina expression beadarrays in IDAT format.
+
+\item
+The ctrlpath argument of read.ilmn() now defaults to the same as path for regular probes.
+This means that only one path setting is required if the regular and control probe profiles are in the same directory.
+
+\item
+read.ilmn() now sets the same probe IDs as rownames for both the expression matrix E and the annotation data.frame genes, providing that the probe IDs are unique.
+
+\item
+beadCountWeights() can now work with either probe-wise standard errors or probe-wise standard deviations.
+
+\item
+treat() has new arguments robust and winsor.tail.p which are passed through to robust empirical Bayes estimation.
+
+\item
+topTreat() now includes ... argument which is passed to topTable().
+
+\item
+topTable() with confint=TRUE now produces confidence intervals based on the t-distribution instead of on the normal distribution.
+It also now accepts a numeric value for the confint argument to specify a confidence level other the default of 0.95.
+
+\item
+topTable() will now work on an MArrayLM fit object that is missing the lods component, for example as produced by treat().
+
+\item
+roast() and mroast() now permit array weights and observation weights to both be specified.
+
+\item
+camera(), roast() and mroast() now use getEAWP() to interpret the data object. This means that they now work on any class of data object that lmFit() will.
+
+\item
+romer() now uses propTrueNull(method="lfdr") instead of convest().
+This makes it substantially faster when the number of genes is large.
+
+\item
+genas() now uses fit$df.total from the MArrayLM object.
+This prevents df.total from exceeding the total pooled residual df for the dataset.
+The genas() results will change slightly for datasets for which df.prior was very lage.
+
+\item
+plotDensities() is now an S3 generic function with methods for RGList, MAList, EListRaw and EList objects.
+
+\item
+plotFB is now an S3 generic function with methods for RGList and EList data objects.
+
+\item
+New topic help pages 10GeneSetTests.Rd and 11RNAseq.Rd.
+The page 10Other.Rd is deleted.
+All topic help pages are now listed under 'See also' in the package introduction page accessible by ?limma.
+
+\item
+avereps() was never intended to be applied to RGList or EListRaw objects.
+It now gives an error when applied to these objects instead of returning a matrix of questionable value.
+
+\item
+Bug fix: fitFDistRobustly() was failing when there were missing values or zero df values and covariate was NULL.
+
+\item
+Bug fix: vennDiagram() wasn't passing extra arguments (...) to plot() when the number of sets was greater than 3.
+
+\item
+Bug fix to topTreat(). Rownames were incorrectly ordered when p<1.
+
+\item
+bug fix to genas(), which was not handling vector df.prior correctly when the fit object was generated using robust=TRUE.
+
+\item
+bug fix to squeezeVar().
+Previously there was an error when robust=TRUE and trend=FALSE and some of the estimated df.prior were infinite.
+
+\item
+bug fix to topTable() and topTableF() when sorting by F-statistic combined with p-value or lfc cutoffs.
+}}
+
+
\section{Version 3.18.0}{\itemize{
\item
diff --git a/inst/doc/changelog.txt b/inst/doc/changelog.txt
index b352173..f075d59 100755
--- a/inst/doc/changelog.txt
+++ b/inst/doc/changelog.txt
@@ -1,35 +1,218 @@
-31 January 2014: limma 3.18.10
+24 June 2014: limma 3.20.7
+
+- update to helpMethods().
+
+- bug fix to topTableF() and topTable() ordering of Amean values when
+ ranking is by F-statistic and a lfc or p.value filter has been set.
+
+18 June 2014: limma 3.20.6
+
+- improvements to roast() and mroast(). The directional (up and
+ down) tests done by roast() now use both the original rotations
+ and their opposite signs, effectively doubling the number of
+ effective rotations. The two-sided tests are now done explicitly
+ by rotation instead of doubling the smallest one-sided p-value.
+ The two-sided p-value is now called "UpOrDown" in the roast()
+ output.
+
+16 June 2014: limma 3.20.5
+
+- bug fix to read.ilmn() when sep=",".
+
+- biocViews keywords updated.
+
+- roast() and mroast() have a new argument 'approx.zscore'. By
+ default they now use approx=TRUE when calling zscoreT().
+
+- zscoreT() can now optionally use a fast approximation instead of
+ the slower exact calculation.
+
+- plotMDS() can now optionally plot samples using symbols instead of
+ text labels.
+
+- many Rd files have been reformated to avoid over long lines in the
+ pdf manual.
+
+- Further change to removeBatchEffect(), which now calls lmFit() to
+ estimate the batch effect. This means that removeBatchEffect() can
+ now take into account weights and other arguments that will affect
+ the linear model fit. This also corrects a problem, introduced by
+ the previous change, when the batch effect was confounded with the
+ design matrix.
+
+20 May 2014: limma 3.20.4
+
+- Change to removeBatchEffect() when there is a non-trivial design
+ matrix. Previously batch adjustments were made only within the
+ treatment levels defined by the design matrix. The behavior of
+ removeBatchEffect() is now more consistent that of linear modelling
+ with lmFit() with batches as additive effects.
+
+19 May 2014: limma 3.20.3
+
+- Remove unnecessary call to eBayes() from topTable when ranking by
+ F-statistic for a subset of the coefficients.
+
+13 May 2014: limma 3.20.2
+
+- Subsetting columns of a MArrayLM object should not subset the
+ design matrix. Now fixed.
+
+12 April 2014: limma 3.20.1
+
+- Update to Pasilla case study in User's Guide.
+
+12 April 2014: limma 3.20.0 (Bioconductor 2.14 Release Branch)
+
+ 7 April 2014: limma 3.19.30
+
+- Update to Pasilla case study in User's Guide.
+
+30 March 2014: limma 3.19.29
+
+- duplicateCorrelation() now traps errors arising from mixedModel2()
+ or La.svd(). Rows that produce errors are assigned an NA
+ correlation which does not contribute to the consensus value.
+
+21 March 2014: limma 3.19.28
+
+- updating NEWS.Rd ready for Bioconductor 2.14 release.
+
+21 March 2014: limma 3.19.27
+
+- new function read.idat() to read files from Illumina expression
+ beadarrays in IDAT format.
+
+- new plot() methods for RGList, MAList, EList and MArrayLM class
+ objects. In each case, produces similar result to plotMA().
+
+- the array argument of plotMA.MArrayLM has been renamed to coef, and
+ the default is changed from 1 to ncol(MA).
+
+21 March 2014: limma 3.19.26
+
+- Updates to help pages. The topic page 10Other.Rd is deleted.
+ The topic pages GeneSetTesting.Rd and RNAseq.Rd renamed to
+ 10GeneSetTests.Rd and 11RNAseq.Rd and their content is expanded.
+ The example for weightedLowess.Rd is revised to take less time.
+
+20 March 2014: limma 3.19.25
+
+- improvements to plotSplice(), which now highlights signficant exons
+ and plots exon position annotation.
+
+- New Pasilla case study added to User's Guide, demonstrating
+ differential splicing analysis of RNA-seq data.
+
+- New topic help pages for Gene Set Testing and RNA-seq analysis.
+ The topic pages are now linked from the package introduction page
+ accessible by ?limma.
+
+15 March 2014: limma 3.19.24
+
+- new functions diffSplice(), topSplice() and plotSplice() to detect
+ differentially spliced genes from exon-level expression data.
+
+5 March 2014: limma 3.19.23
+
+- biocViews updated in DESCRIPTION file.
+
+28 February 2014: limma 3.19.22
+
+- topTreat() now includes ... argument which is passed to topTable().
+
+- camera(), roast() and mroast() restored as S3 generic functions
+ with a default method only.
+
+27 February 2014: limma 3.19.21
+
+- camera(), roast() and mroast() now use getEAWP() to interpret the
+ data object. This means that they now work on any class of data
+ object that lmFit() will. They are no longer generic functions.
+
+- topTable() with confint=TRUE now produces confidence intervals
+ based on the t-distribution instead of on the normal distribution.
+ It also now accepts a numeric value for the confint argument to
+ specify a confidence level other the default of 0.95.
+
+21 February 2014: limma 3.19.20
+
+- Bug fix: fitFDistRobustly() was failing when there were missing
+ values or zero df values and covariate was NULL.
+
+18 February 2014: limma 3.19.19
+
+- new function tricubeMovingAverage() for smoothing a time series.
+
+- The enrichment worms shown by barcodeplot() now use
+ tricubeMovingAverage(). This makes the worms much smoother than
+ previously.
+
+- plotMA3by2() now accepts single channel data objects as well as
+ two color.
+
+- Bug fix: the mroast methods for EList and MAList objects were
+ ignoring the midp, p.adjust and sorting arguments.
+
+11 February 2014: limma 3.19.18
+
+- The new subsetting code of 13 Dec 2013 introduced a bug whereby
+ MArrayLM objects couldn't be subsetted by column name. Now fixed.
+
+ 5 February 2014: limma 3.19.17
+
+- The change to read.ilmn() in 3.19.16 caused an error when the
+ Illumina probe IDs were not all unique. The probe IDs are now used
+ as row.names for the annotation data.frame only when they are
+ unique.
+
+31 January 2014: limma 3.19.16
- treat() now has new arguments robust and winsor.tail.p which are
passed through to robust empirical Bayes estimation.
+- barcodeplot() has a new option to add enrichment worms to the plot.
+
- vennDiagram() wasn't passing extra arguments (...) to plot() when
the number of sets was greater than 3. Now fixed.
- read.ilmn() now sets the same probe IDs as rownames for both the
expression matrix E and the annotation data.frame genes.
+- beadCountWeights() can now work with either probe-wise standard
+ errors or probe-wise standard deviations.
+
+- new function weightedLoess() which fits a lowess curve with
+ prior weights. Unlike previous implementations of lowess or loess,
+ the weights are used in calculating which neighbouring points to
+ include in each local regression as well as in the local
+ regression itself.
+
+- weightedLoess() now becomes the default method used by loessFit()
+ to fit the loess curve when there are weights. The previous locfit
+ and loess() methods are offered as options.
+
+14 Jan 2014: limma 3.19.15
+
- romer() now uses propTrueNull(method="lfdr") instead of convest().
This makes romer() substantially faster when the number of genes
is large.
-11 Jan 2014: limma 3.18.9
+11 Jan 2014: limma 3.19.14
- update citation of Law et al (2014) reference for voom().
-11 Jan 2014: limma 3.18.8
-
- roast() and mroast() now permit array weights and observation
weights to both be specified.
- bug fix to mroast() which was not accepting array.weights.
-17 Dec 2013: limma 3.18.7
+17 Dec 2013: limma 3.19.13
- make sure that F is recomputed when the columns of an MArrayLM
are subsetted.
-12 Dec 2013: limma 3.18.6
+12 Dec 2013: limma 3.19.12
- linear model fit functions lm.series(), mrlm.series() and
gls.series() no longer drop the dimensions of the components of
@@ -37,21 +220,30 @@
Previously this was done inconsistently in some cases but not
others. Now the matrix components always keep dimensions.
-- New function subsetListOfArrays(), which is used to simply the
+11 Dec 2013: limma 3.19.11
+
+- New function subsetListOfArrays(), which is used to simplify the
subsetting code for RGList, MAList, EList, EListRaw and MArrayLM
objects.
- 7 Dec 2013: limma 3.18.5
+ 7 Dec 2013: limma 3.19.10
- Bug fix to topTreat(). Rownames were incorrectly ordered when p<1.
- topTable() will now work on an MArrayLM fit object that is missing
the lods component, for example as produced by treat().
-- lmFit(), eBayes(), tmixture.vector() and so on now work even when
- there is only one gene (one row of data).
+ 6 Dec 2013: limma 3.19.9
+
+- Fix to tmixture.vector() to allow it work even with just one gene.
+ This allows the eBayes() to work on fit objects with just one row.
+
+ 3 Dec 2013: limma 3.19.8
- 2 Dec 2013: limma 3.18.4
+- lmFit() now works even when there is only one gene (one row of
+ data).
+
+ 2 Dec 2013: limma 3.19.7
- bug fix to genas(), which was not handling vector df.prior
correctly when the fit object was generated using robust=TRUE.
@@ -74,33 +266,45 @@
objects. It now gives an error when applied to these objects
instead of returning a matrix of questionable value.
+18 Nov 2013: limma 3.19.6
+
- plotDensities() is now an S3 generic function with methods for
RGList, MAList, EListRaw and EList objects.
-13 Nov 2013: limma 3.18.3
+13 Nov 2013: limma 3.19.5
+
+- NAMESPACE now exports plotFB methods.
+
+- bug fix to x-axis label for some plotMA methods.
+
+13 Nov 2013: limma 3.19.4
- plotFB now an S3 generic function with methods for RGList and EList
data objects.
-- edits to help pages for read.ilmn, neqc, beadCountWeights, voom and
- vooma. Update cross link to affy::normalize-method in
- normalizeBetweenArrays.Rd. Remove cross reference to ROC package
- in auROC.Rd.
+10 Nov 2013: limma 3.19.3
-- Vignette limma.Rnw renamed to intro.Rnw and moved to vignette
- directory.
-
-- bug fix to x-axis label for some plotMA methods.
+- edits to help pages for read.ilmn, neqc, beadCountWeights, voom and
+ vooma.
- 4 Nov 2013: limma 3.18.2
+ 4 Nov 2013: limma 3.19.2
- bug fix to topTable() and topTableF() when sorting by F-statistic
combined with p-value or lfc cutoffs.
-23 Oct 2013: limma 3.18.1
+23 Oct 2013: limma 3.19.1
+
+- Update cross reference to affy package in
+ normalizeBetweenArrays.Rd.
+
+- Revise auROC.Rd help page. Remove cross reference to ROC package.
+
+- Vignette limma.Rnw renamed to intro.Rnw and moved to vignette
+ directory.
- NEWS.Rd updated to reflect changes in Bioconductor 2.13 Release.
+15 Oct 2013: limma 3.19.0 (Bioconductor 2.14 Developmental Branch)
15 Oct 2013: limma 3.18.0 (Bioconductor 2.13 Release Branch)
11 Oct 2013: limma 3.17.28
@@ -365,6 +569,9 @@
are included in the test set. Previously this argument was called
'iset' for roast() and romer() and 'indices' for camera().
+- When sorting results, mroast() now tries to break ties more
+ sensibly rather than keeping original order.
+
- section on time course experiments with many time points added to
User's Guide.
diff --git a/inst/doc/intro.pdf b/inst/doc/intro.pdf
index 8a231a7..78b9b67 100644
Binary files a/inst/doc/intro.pdf and b/inst/doc/intro.pdf differ
diff --git a/man/01Introduction.Rd b/man/01Introduction.Rd
index 0afd381..d008fc5 100644
--- a/man/01Introduction.Rd
+++ b/man/01Introduction.Rd
@@ -14,17 +14,17 @@ The linear model and differential expression functions apply to all microarrays
\details{
There are three types of documentation available:
-\tabular{ll}{
-(1) \tab
+\enumerate{
+\item
The \emph{LIMMA User's Guide} can be reached through the "User
Guides and Package Vignettes" links at the top of the LIMMA
contents page. The function \code{\link{limmaUsersGuide}} gives
the file location of the User's Guide.\cr
-(2) \tab
+\item
An overview of limma functions grouped by purpose is contained
-in the numbered chapters at the top of the LIMMA contents page,
+in the numbered chapters at the foot of the LIMMA package index page,
of which this page is the first.\cr
-(3) \tab
+\item
The LIMMA contents page gives an
alphabetical index of detailed help topics.\cr
}
@@ -32,9 +32,11 @@ alphabetical index of detailed help topics.\cr
The function \code{\link{changeLog}} displays the record of changes to the package.
}
-\author{Gordon Smyth}
+\author{Gordon Smyth, with contributions from many colleagues}
+
\references{
-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., 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.
@@ -44,4 +46,18 @@ Smyth, G. K. (2005). Limma: linear models for microarray data.
In: \emph{Bioinformatics and Computational Biology Solutions using R and Bioconductor}.
R. Gentleman, V. Carey, S. Dudoit, R. Irizarry, W. Huber (eds), Springer, New York, 2005.
}
+
+\seealso{
+\link{02.Classes},
+\link{03.ReadingData},
+\link{04.Background},
+\link{05.Normalization},
+\link{06.LinearModels},
+\link{07.SingleChannel},
+\link{08.Tests},
+\link{09.Diagnostics},
+\link{10.GeneSetTests},
+\link{11.RNAseq}
+}
+
\keyword{package}
diff --git a/man/02classes.Rd b/man/02classes.Rd
index 8ef7839..5c308f2 100644
--- a/man/02classes.Rd
+++ b/man/02classes.Rd
@@ -1,6 +1,6 @@
\name{02.Classes}
\alias{02.Classes}
-\title{Classes Defined by this Package}
+\title{Topic: Classes Defined by this Package}
\description{
diff --git a/man/03reading.Rd b/man/03reading.Rd
index b152c5d..46c9961 100644
--- a/man/03reading.Rd
+++ b/man/03reading.Rd
@@ -1,6 +1,6 @@
\name{03.ReadingData}
\alias{03.ReadingData}
-\title{Reading Microarray Data from Files}
+\title{Topic: Reading Microarray Data from Files}
\description{
This help page gives an overview of LIMMA functions used to read data from files.
@@ -25,6 +25,8 @@ There are also a series of utility functions which read the header information f
\code{\link{read.ilmn}} reads probe or gene summary profile files from Illumina BeadChips,
and produces an \code{ElistRaw} object.
+\code{\link{read.idat}} reads Illumina files in IDAT format, and produces an \code{EListRaw} object.
+
The function \link{as.MAList} can be used to convert a \code{marrayNorm} object to an \code{MAList} object if the data was read and normalized using the marray and marrayNorm packages.
}
diff --git a/man/04Background.Rd b/man/04Background.Rd
index cb4b8e7..3979a48 100644
--- a/man/04Background.Rd
+++ b/man/04Background.Rd
@@ -1,6 +1,6 @@
\name{04.Background}
\alias{04.Background}
-\title{Background Correction}
+\title{Topic: Background Correction}
\description{
This page deals with background correction methods provided by the \code{\link{backgroundCorrect}}, \code{\link{kooperberg}} or \code{\link{neqc}} functions.
diff --git a/man/05Normalization.Rd b/man/05Normalization.Rd
index 5b85bc6..b715086 100644
--- a/man/05Normalization.Rd
+++ b/man/05Normalization.Rd
@@ -1,6 +1,6 @@
\name{05.Normalization}
\alias{05.Normalization}
-\title{Normalization of Microarray Data}
+\title{Topic: Normalization of Microarray Data}
\description{
This page gives an overview of the LIMMA functions available to normalize data from single-channel or two-colour microarrays.
diff --git a/man/06linearmodels.Rd b/man/06linearmodels.Rd
index 1d46547..8992d56 100644
--- a/man/06linearmodels.Rd
+++ b/man/06linearmodels.Rd
@@ -1,6 +1,6 @@
\name{06.LinearModels}
\alias{06.LinearModels}
-\title{Linear Models for Microarrays}
+\title{Topic: Linear Models for Microarrays}
\description{
This page gives an overview of the LIMMA functions available to fit linear models and to interpret the results.
diff --git a/man/07SingleChannel.Rd b/man/07SingleChannel.Rd
index 520cab4..0a0da0e 100644
--- a/man/07SingleChannel.Rd
+++ b/man/07SingleChannel.Rd
@@ -1,6 +1,6 @@
\name{07.SingleChannel}
\alias{07.SingleChannel}
-\title{Individual Channel Analysis of Two-Color Microarrays}
+\title{Topic: Individual Channel Analysis of Two-Color Microarrays}
\description{
This page gives an overview of the LIMMA functions fit linear models to two-color microarray data in terms of the log-intensities rather than log-ratios.
diff --git a/man/08Tests.Rd b/man/08Tests.Rd
index 4fbf79c..fcde541 100644
--- a/man/08Tests.Rd
+++ b/man/08Tests.Rd
@@ -1,6 +1,6 @@
\name{08.Tests}
\alias{08.Tests}
-\title{Hypothesis Testing for Linear Models}
+\title{Topic: Hypothesis Testing for Linear Models}
\description{
LIMMA provides a number of functions for multiple testing across both contrasts and genes.
diff --git a/man/09Diagnostics.Rd b/man/09Diagnostics.Rd
index a34abea..97754b9 100644
--- a/man/09Diagnostics.Rd
+++ b/man/09Diagnostics.Rd
@@ -1,6 +1,6 @@
\name{09.Diagnostics}
\alias{09.Diagnostics}
-\title{Diagnostics and Quality Assessment}
+\title{Topic: Diagnostics and Quality Assessment}
\description{
This page gives an overview of the LIMMA functions available for microarray quality assessment and diagnostic plots.
@@ -22,7 +22,7 @@ Produces a spatial picture of any spot-specific measure from an array image. If
\item{ \code{\link{plotFB}} }{
Plots foreground versus background log-intensies.}
-\item{ \code{\link{plotMA}} }{
+\item{ \code{\link{plotMA}} or \code{\link{plot.MAList}} }{
MA-plots.
One of the most useful plots of a two-color array.
\code{\link{plotMA3by2}} will write MA-plots to files, six plots to a page.
diff --git a/man/10GeneSetTests.Rd b/man/10GeneSetTests.Rd
new file mode 100644
index 0000000..f991624
--- /dev/null
+++ b/man/10GeneSetTests.Rd
@@ -0,0 +1,36 @@
+\name{10.GeneSetTests}
+\alias{10.GeneSetTests}
+\title{Topic: Gene Set Tests}
+
+\description{
+This page gives an overview of the LIMMA functions for gene set testing and pathway analysis.
+
+\describe{
+\item{ \code{\link{roast}} }{
+ Self-contained gene set testing for one set.}
+
+\item{ \code{\link{mroast}} }{
+ Self-contained gene set testing for many sets.}
+
+\item{ \code{\link{camera}} }{
+ Competitive gene set testing.}
+
+\item{ \code{\link{romer}} }{
+ Gene set enrichment analysis.}
+
+\item{ \code{\link{symbols2indices}} }{
+ Convert gene sets consisting of vectors of symbols into a list of indices suitable for use in the above functions.}
+
+\item{ \code{\link{alias2Symbol}} and \code{\link{alias2SymbolTable}} }{
+ Convert gene symbols or aliases to current official symbols.}
+
+\item{ \code{\link{topRomer}} }{
+ Display results from romer tests as a top-table.}
+
+\item{ \code{\link{geneSetTest}} or \code{\link{wilcoxGST}} }{
+ Simple gene set testing based on gene or probe permutation.}
+
+\item{ \code{\link{barcodeplot}} }{
+ Enrichment plot of a gene set.}
+}
+}
diff --git a/man/10Other.Rd b/man/10Other.Rd
deleted file mode 100644
index 3b008a5..0000000
--- a/man/10Other.Rd
+++ /dev/null
@@ -1,10 +0,0 @@
-\name{10.Other}
-\alias{10.Other}
-\title{Other Functions}
-
-\description{
-This page describes some functions not covered in the previous numbered pages, so far only \code{\link{blockDiag}} and \code{\link{poolVar}} which are not used in the package yet but are part of the development of methods to handle technical and biological replicates.
-}
-
-\author{Gordon Smyth}
-\keyword{documentation}
diff --git a/man/11RNAseq.Rd b/man/11RNAseq.Rd
new file mode 100644
index 0000000..d66e456
--- /dev/null
+++ b/man/11RNAseq.Rd
@@ -0,0 +1,31 @@
+\name{11.RNAseq}
+\alias{11.RNAseq}
+\title{Topic: Analysis of RNA-seq Data}
+
+\description{
+This page gives an overview of the LIMMA functions available to analyze RNA-seq data.
+
+\describe{
+\item{ \code{\link{voom}} }{
+ Transform RNA-seq or ChIP-seq counts to log counts per million (log-cpm) with associated precision weights.
+ After this tranformation, RNA-seq or ChIP-seq data can be analyzed using the same functions as would be used for microarray data.}
+
+\item{ \code{\link{diffSplice}} }{
+ Test for differential splicing of exons between experimental conditions.}
+
+\item{ \code{\link{topSplice}} }{
+ Show a data.frame of top results from \code{diffSplice}.}
+
+\item{ \code{\link{plotSplice}} }{
+ Plot results from \code{diffSplice}.}
+}
+}
+
+\references{
+Law, CW, Chen, Y, Shi, W, Smyth, GK (2014).
+Voom: precision weights unlock linear model analysis tools for RNA-seq read counts.
+\emph{Genome Biology} 15, R29.
+\url{http://genomebiology.com/2014/15/2/R29}
+}
+
+\keyword{documentation}
diff --git a/man/EList.Rd b/man/EList.Rd
index 63f31b8..1f05e52 100644
--- a/man/EList.Rd
+++ b/man/EList.Rd
@@ -13,20 +13,21 @@ A list-based S4 classes for storing expression values (E-values) for a set of on
}
\section{Components}{
-\code{EList} objects can be created by \code{new("EList",E)} where \code{E} is a list.
-These classes contains no slots (other than \code{.Data}), but objects should contain a list component \code{E} as follows:
-\tabular{ll}{
- \code{E} \tab numeric matrix containing expression values. Rows correspond to probes and columns to arrays.
- In an \code{EListRaw} object, these are unlogged values.
- In an \code{EList} object, these are log2 values.
-}
+\code{EList} objects can be created by \code{new("EList",x)} where \code{x} is a list.
+These classes contains no slots (other than \code{.Data}), but objects should contain a list component \code{E},
+which is a numeric matrix containing expression values.
+Rows correspond to probes and columns to arrays.
+In an \code{EListRaw} object, \code{E} contains unlogged values.
+In an \code{EList} object, \code{E} contains log2 values.
+
Optional components include:
-\tabular{ll}{
- \code{weights} \tab numeric matrix of same dimensions as \code{E} containing relative spot quality weights. Elements should be non-negative.\cr
- \code{other} \tab list containing other matrices, all of the same dimensions as \code{E}.\cr
- \code{genes} \tab data.frame containing probe information. Should have one row for each probe. May have any number of columns.\cr
- \code{targets} \tab data.frame containing information on the target RNA samples. Rows correspond to arrays. May have any number of columns.
+\describe{
+ \item{\code{weights}}{numeric matrix of same dimensions as \code{E} containing relative spot quality weights. Elements should be non-negative.}
+ \item{\code{other}}{list containing other matrices, all of the same dimensions as \code{E}.}
+ \item{\code{genes}}{data.frame containing probe information. Should have one row for each probe. May have any number of columns.}
+ \item{\code{targets}}{data.frame containing information on the target RNA samples. Rows correspond to arrays. May have any number of columns.}
}
+
Valid \code{EList} or \code{EListRaw} objects may contain other optional components, but all probe or array information should be contained in the above components.
}
diff --git a/man/PrintLayout.Rd b/man/PrintLayout.Rd
index 0a7fa73..78c5c2e 100755
--- a/man/PrintLayout.Rd
+++ b/man/PrintLayout.Rd
@@ -35,15 +35,18 @@ Objects of this class contains no slots but should contain the following list co
\examples{
# Settings for Swirl and ApoAI example data sets in User's Guide
-printer <- list(ngrid.r=4, ngrid.c=4, nspot.r=22, nspot.c=24, ndups=1, spacing=1, npins=16, start="topleft")
+printer <- list(ngrid.r=4, ngrid.c=4, nspot.r=22, nspot.c=24,
+ ndups=1, spacing=1, npins=16, start="topleft")
# Typical settings at the Australian Genome Research Facility
# Full pin set, duplicates side-by-side on same row
-printer <- list(ngrid.r=12, ngrid.c=4, nspot.r=20, nspot.c=20, ndups=2, spacing=1, npins=48, start="topright")
+printer <- list(ngrid.r=12, ngrid.c=4, nspot.r=20, nspot.c=20,
+ ndups=2, spacing=1, npins=48, start="topright")
# Half pin set, duplicates in top and lower half of slide
-printer <- list(ngrid.r=12, ngrid.c=4, nspot.r=20, nspot.c=20, ndups=2, spacing=9600, npins=24, start="topright")
+printer <- list(ngrid.r=12, ngrid.c=4, nspot.r=20, nspot.c=20,
+ ndups=2, spacing=9600, npins=24, start="topright")
}
\keyword{classes}
diff --git a/man/alias2Symbol.Rd b/man/alias2Symbol.Rd
index 6e9b708..c1e1510 100644
--- a/man/alias2Symbol.Rd
+++ b/man/alias2Symbol.Rd
@@ -1,7 +1,7 @@
\name{alias2Symbol}
\alias{alias2Symbol}
\alias{alias2SymbolTable}
-\title{Convert Gene Alias to Official Gene Symbols}
+\title{Convert Gene Aliases to Official Gene Symbols}
\description{
Maps gene alias names to official gene symbols.
}
@@ -23,7 +23,7 @@ The user needs to have the appropriate Bioconductor organism package installed.
\code{alias2Symbol} maps a set of aliases to a set of symbols, without necessarily preserving order.
The output vector may be longer or shorter than the original vector, because some aliases might not be found and some aliases may map to more than one symbol.
\code{alias2SymbolTable} maps each alias to a gene symbol and returns a table with one row for each alias.
-If an alias maps to more than one symbol, then the first one found will be returned.
+If an alias maps to more than one symbol, then the first one found is returned.
}
\value{
Character vector of gene symbols.
@@ -34,9 +34,8 @@ Character vector of gene symbols.
\author{Gordon Smyth and Yifang Hu}
\seealso{
This function is often used to assist gene set testing, see
-\link{08.Tests}.
+\link{10.GeneSetTests}.
}
\examples{
if(require("org.Hs.eg.db")) alias2Symbol(c("PUMA","NOXA","BIM"))
}
-\keyword{character}
diff --git a/man/arrayWeights.Rd b/man/arrayWeights.Rd
index 3cffdf3..18fc209 100644
--- a/man/arrayWeights.Rd
+++ b/man/arrayWeights.Rd
@@ -7,8 +7,10 @@ Estimates relative quality weights for each array in a multi-array
experiment.
}
\usage{
-arrayWeights(object, design = NULL, weights = NULL, method = "genebygene", maxiter = 50, tol = 1e-10, trace=FALSE)
-arrayWeightsSimple(object, design = NULL, maxiter = 100, tol = 1e-6, maxratio = 100, trace=FALSE)
+arrayWeights(object, design = NULL, weights = NULL, method = "genebygene", maxiter = 50,
+ tol = 1e-10, trace=FALSE)
+arrayWeightsSimple(object, design = NULL,
+ maxiter = 100, tol = 1e-6, maxratio = 100, trace=FALSE)
}
\arguments{
\item{object}{object of class \code{numeric}, \code{matrix}, \code{MAList}, \code{marrayNorm},
@@ -52,7 +54,10 @@ the slots or components in the data \code{object}.
A vector of array weights.
}
\references{
-Ritchie, M. E., Diyagama, D., Neilson, van Laar, R., J., Dobrovic, A., Holloway, A., and Smyth, G. K. (2006). Empirical array quality weights in the analysis of microarray data. BMC Bioinformatics 7, 261. \url{http://www.biomedcentral.com/1471-2105/7/261/abstract}
+Ritchie, M. E., Diyagama, D., Neilson, van Laar, R., J., Dobrovic, A., Holloway, A., and Smyth, G. K. (2006).
+Empirical array quality weights in the analysis of microarray data.
+BMC Bioinformatics 7, 261.
+\url{http://www.biomedcentral.com/1471-2105/7/261/abstract}
}
\seealso{
An overview of linear model functions in limma is given by \link{06.LinearModels}.
diff --git a/man/arrayWeightsQuick.Rd b/man/arrayWeightsQuick.Rd
index 7eb4eb0..4328877 100644
--- a/man/arrayWeightsQuick.Rd
+++ b/man/arrayWeightsQuick.Rd
@@ -21,7 +21,10 @@ This is a quick and dirty version of \code{\link{arrayWeights}}.
Numeric vector of weights of length \code{ncol(fit)}.
}
\references{
-Ritchie, M. E., Diyagama, D., Neilson, van Laar, R., J., Dobrovic, A., Holloway, A., and Smyth, G. K. (2006). Empirical array quality weights for microarray data. BMC Bioinformatics. (Accepted 11 April 2006)
+Ritchie, M. E., Diyagama, D., Neilson, van Laar, R., J., Dobrovic, A., Holloway, A., and Smyth, G. K. (2006).
+Empirical array quality weights in the analysis of microarray data.
+BMC Bioinformatics 7, 261.
+\url{http://www.biomedcentral.com/1471-2105/7/261/abstract}
}
\author{Gordon Smyth}
\seealso{
@@ -34,4 +37,3 @@ fit <- lmFit(y, design)
arrayWeightsQuick(y, fit)
}
}
-\keyword{regression}
diff --git a/man/asMatrixWeights.Rd b/man/asMatrixWeights.Rd
index 5429a26..d1b2371 100644
--- a/man/asMatrixWeights.Rd
+++ b/man/asMatrixWeights.Rd
@@ -36,4 +36,3 @@ asMatrixWeights(1:4,c(4,3))
An overview of functions in LIMMA used for fitting linear models is given in \link{06.LinearModels}.
}
-\keyword{hplot}
diff --git a/man/avearrays.Rd b/man/avearrays.Rd
index 784e860..f1ddcd7 100644
--- a/man/avearrays.Rd
+++ b/man/avearrays.Rd
@@ -27,7 +27,7 @@ For an \code{MAList} object, the components \code{M} and \code{A} are both avera
If \code{x} is of mode \code{"character"}, then the replicate values are assumed to be equal and the first is taken as the average.
}
\value{
-A data object of the same class as \code{x} with a row for each unique value of \code{ID}.
+A data object of the same class as \code{x} with a column for each unique value of \code{ID}.
}
\author{Gordon Smyth}
\seealso{
@@ -41,4 +41,3 @@ x <- matrix(rnorm(8*3),8,3)
colnames(x) <- c("a","a","b")
avearrays(x)
}
-\keyword{array}
diff --git a/man/barcodeplot.Rd b/man/barcodeplot.Rd
index 56037de..216f91b 100644
--- a/man/barcodeplot.Rd
+++ b/man/barcodeplot.Rd
@@ -5,17 +5,18 @@
Plot positions of one or two gene sets in a ranked list of statistics.
}
\usage{
-barcodeplot(statistics,index,index2=NULL,labels=c("Group 1 (highest statistic)","Group 2 (lowest statistic)"),
- quantiles=c(-1,1),col.bars=NULL,offset.bars=!is.null(index2), ...)
+barcodeplot(statistics, index, index2=NULL, labels=c("Up","Down"), quantiles=c(-1,1),
+ col.bars=NULL, worm=TRUE, span.worm=0.45, ...)
}
\arguments{
\item{index}{index vector for the gene set. This can be a vector of indices, or a logical vector of the same length as \code{statistics} or, in general, any vector such that \code{statistic[index]} gives the statistic values for the gene set to be tested.}
\item{index2}{index vector for a second gene set. Usually used to specify down-regulated genes when \code{index} is used for up-regulated genes.vector specifying the elements of \code{statistic} in the test group.}
\item{statistics}{numeric vector giving the values of statistics to rank genes by.}
- \item{labels}{character vector of labels for the two groups of RNA samples that are being compared by the statistics. First labe is associated with high statistics and is displayed at the left end of the plot. Second label is associated with low or negative statistics and is displayed at the right end of the plot.}
+ \item{labels}{character vector of labels for the two groups of RNA samples that are being compared by the statistics. First label is associated with high statistics and is displayed at the left end of the plot. Second label is associated with low or negative statistics and is displayed at the right end of the plot.}
\item{quantiles}{numeric vector of length 2, giving cutoff values for \code{statistics} considered small or large respectively. Used to color the rectangle of the barcodeplot.}
\item{col.bars}{character vector giving colors for the bars on the barcodeplot. Defaults to \code{"black"} for one set or \code{c("red","blue")} for two sets.}
- \item{offset.bars}{logical. When there are two sets, bars for the first set are normally offset up and and bars for the second set are offset down from the rectangle of the barcodeplot.}
+ \item{worm}{logical, should enrichment worms be plotted?}
+ \item{span.worm}{loess span for enrichment worms.}
\item{\ldots}{other arguments are passed to \code{plot}.}
}
@@ -25,16 +26,18 @@ No value is returned but a plot is produced as a side effect.
\details{
Rank statistics left to right from largest to smallest, and show positions of one or two specified subsets.
-This plot is typically used in conjunction with gene set tests.
+This plot is intended to be used in conjunction with gene set tests, and shows the enrichment of gene sets amongst high or low ranked genes.
It first appeared in the literature in Lim et al (2009).
It was inspired by the set location plot of Subramanian et al (2005).
+
+The enrichment worm is calculated by tricube-weighted moving average.
}
\seealso{
-\code{\link{camera}}, \code{\link{wilcox.test}}
+\code{\link{tricubeMovingAverage}}, \code{\link{camera}}, \code{\link{wilcox.test}}
}
-\author{Gordon Smyth and Di Wu}
+\author{Gordon Smyth, Di Wu and Yifang Hu}
\references{
Lim E, Vaillant F, Wu D, Forrest NC, Pal B, Hart AH, Asselin-Labat ML, Gyorki DE, Ward T, Partanen A, Feleppa F, Huschtscha LI, Thorne HJ; kConFab; Fox SB, Yan M, French JD, Brown MA, Smyth GK, Visvader JE, Lindeman GJ (2009).
diff --git a/man/beadCountWeights.Rd b/man/beadCountWeights.Rd
index ff26b5f..2bbbc7c 100644
--- a/man/beadCountWeights.Rd
+++ b/man/beadCountWeights.Rd
@@ -8,17 +8,19 @@ Estimates weights which account for biological variation and technical variation
}
\usage{
-beadCountWeights(y, x, design = NULL, bead.stdev = NULL, nbeads = NULL, array.cv = TRUE, scale = FALSE)
+beadCountWeights(y, x, design = NULL, bead.stdev = NULL, bead.stderr = NULL,
+ nbeads = NULL, array.cv = TRUE, scale = FALSE)
}
\arguments{
\item{y}{normalized log2-expression values.}
- \item{x}{raw expression values.}
+ \item{x}{raw expression values, with the same dimensions as \code{y}.}
\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{bead.stderr}{numeric matrix containing bead-level standard errors.}
\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.}
@@ -42,7 +44,7 @@ Coefficients of variation are calculated using raw expression values.
The biological variance for each probe across the arrays are estimated using a Newton iteration, with the assumption that the total residual deviance for each probe from \code{lmFit} is inversely proportional to the sum of the technical variance and biological variance.
-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}.
+If any of the arguments \code{design}, \code{bead.stdev}, \code{bead.stderr} 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}, \code{bead.stderr} and \code{nbeads} are not set explicitly in the call, it is necessary that they are available for extraction [...]
}
\value{
diff --git a/man/blockDiag.Rd b/man/blockDiag.Rd
index 3e3c2e9..dc011f4 100755
--- a/man/blockDiag.Rd
+++ b/man/blockDiag.Rd
@@ -5,11 +5,15 @@
\description{Form a block diagonal matrix from the given blocks.}
\usage{
-blockDiag(...)
+blockDiag(\dots)
}
\arguments{
-\item{...}{numeric matrices}
+\item{\dots}{numeric matrices}
+}
+
+\details{
+This function is sometimes useful for constructing a design matrix for a disconnected two-color microarray experiment in conjunction with \code{modelMatrix}.
}
\value{A block diagonal matrix with dimensions equal to the sum of the input dimensions}
@@ -23,7 +27,5 @@ blockDiag(a,b)
\author{Gordon Smyth}
\seealso{
-\link{10.Other}
+\code{\link{modelMatrix}}
}
-
-\keyword{array}
diff --git a/man/camera.Rd b/man/camera.Rd
index 89d4507..204eaec 100644
--- a/man/camera.Rd
+++ b/man/camera.Rd
@@ -1,7 +1,5 @@
\name{camera}
\alias{camera}
-\alias{camera.EList}
-\alias{camera.MAList}
\alias{camera.default}
\alias{interGeneCorrelation}
\title{Competitive Gene Set Test Accounting for Inter-gene Correlation}
@@ -9,20 +7,24 @@
Test whether a set of genes is highly ranked relative to other genes in terms of differential expression, accounting for inter-gene correlation.
}
\usage{
-\method{camera}{EList}(y, index, design, contrast=ncol(design), weights=NULL, use.ranks=FALSE, allow.neg.cor=TRUE, trend.var=FALSE, sort=TRUE)
+\S3method{camera}{default}(y, index, design, contrast=ncol(design), weights=NULL,
+ use.ranks=FALSE, allow.neg.cor=TRUE, trend.var=FALSE, sort=TRUE, \dots)
interGeneCorrelation(y, design)
}
\arguments{
- \item{y}{numeric matrix giving log-expression or log-ratio values for a series of microarrays, or any data object that can coerced to a matrix including \code{ExpressionSet}, \code{MAList}, \code{EList} or \code{PLMSet} objects.
- Rows correspond to probes and columns to samples.}
- \item{index}{an index vector or a list of index vectors. Can be any vector such that \code{y[index,]} selects the rows corresponding to the test set. The list can be made using \link{symbols2indices}.}
+ \item{y}{a numeric matrix of log-expression values or log-ratios of expression values, or any data object containing such a matrix.
+ Rows correspond to probes and columns to samples.
+ Any type of object that can be processed by \code{\link{getEAWP}} is acceptable.
+ }
+ \item{index}{an index vector or a list of index vectors. Can be any vector such that \code{y[index,]} selects the rows corresponding to the test set. The list can be made using \code{\link{symbols2indices}}.}
\item{design}{design matrix.}
\item{contrast}{contrast of the linear model coefficients for which the test is required. Can be an integer specifying a column of \code{design}, or else a numeric vector of same length as the number of columns of \code{design}.}
- \item{weights}{can be a numeric matrix of individual weights, of same size as \code{y}, or a numeric vector of array weights with length equal to \code{ncol(y)}.}
- \item{use.ranks}{do a rank-based test (\code{TRUE}) or a parametric test (\code{FALSE}?}
+ \item{weights}{can be a numeric matrix of individual weights, of same size as \code{y}, or a numeric vector of array weights with length equal to \code{ncol(y)}, or a numeric vector of gene weights with length equal to \code{nrow(y)}.}
+ \item{use.ranks}{do a rank-based test (\code{TRUE}) or a parametric test (\code{FALSE})?}
\item{allow.neg.cor}{should reduced variance inflation factors be allowed for negative correlations?}
\item{trend.var}{logical, should an empirical Bayes trend be estimated? See \code{\link{eBayes}} for details.}
\item{sort}{logical, should the results be sorted by p-value?}
+ \item{\dots}{other arguments are not currently used}
}
\details{
@@ -68,11 +70,15 @@ Analyzing gene expression data in terms of gene sets: methodological issues.
}
\seealso{
+\code{\link{getEAWP}}
+
\code{\link{rankSumTestWithCorrelation}},
\code{\link{geneSetTest}},
\code{\link{roast}},
\code{\link{romer}},
\code{\link{symbols2indices}}.
+
+There is a topic page on \link{10.GeneSetTests}.
}
\examples{
diff --git a/man/changelog.Rd b/man/changelog.Rd
index c33a69f..913d3e4 100755
--- a/man/changelog.Rd
+++ b/man/changelog.Rd
@@ -16,6 +16,10 @@ changeLog(n=20)
\author{Gordon Smyth}
+\examples{
+changeLog()
+}
+
\seealso{
\link{01.Introduction}
}
diff --git a/man/classifytests.Rd b/man/classifytests.Rd
index 400c16b..409f31d 100755
--- a/man/classifytests.Rd
+++ b/man/classifytests.Rd
@@ -57,7 +57,9 @@ If \code{tstat} is an \code{MArrayLM} object, then all arguments except for \cod
\code{cor.matrix} is the same as the correlation matrix of the coefficients from which the t-statistics are calculated.
If \code{cor.matrix} is not specified, then it is calculated from \code{design} and \code{contrasts} if at least \code{design} is specified or else defaults to the identity matrix.
In terms of \code{design} and \code{contrasts}, \code{cor.matrix} is obtained by standardizing the matrix
+
\code{ t(contrasts) \%*\% solve(t(design) \%*\% design) \%*\% contrasts }
+
to a correlation matrix.
}
\seealso{
diff --git a/man/controlStatus.Rd b/man/controlStatus.Rd
index ef959a9..c699703 100755
--- a/man/controlStatus.Rd
+++ b/man/controlStatus.Rd
@@ -40,9 +40,15 @@ Attributes contain plotting parameters associated with each spot type.
An overview of LIMMA functions for reading data is given in \link{03.ReadingData}.
}
\examples{
-genes <- data.frame(ID=c("Control","Control","Control","Control","AA1","AA2","AA3","AA4"),
-Name=c("Ratio 1","Ratio 2","House keeping 1","House keeping 2","Gene 1","Gene 2","Gene 3","Gene 4"))
-types <- data.frame(SpotType=c("Gene","Ratio","Housekeeping"),ID=c("*","Control","Control"),Name=c("*","Ratio*","House keeping*"),col=c("black","red","blue"))
+genes <- data.frame(
+ ID=c("Control","Control","Control","Control","AA1","AA2","AA3","AA4"),
+ Name=c("Ratio 1","Ratio 2","House keeping 1","House keeping 2",
+ "Gene 1","Gene 2","Gene 3","Gene 4"))
+types <- data.frame(
+ SpotType=c("Gene","Ratio","Housekeeping"),
+ ID=c("*","Control","Control"),
+ Name=c("*","Ratio*","House keeping*"),
+ col=c("black","red","blue"))
status <- controlStatus(types,genes)
}
\keyword{IO}
diff --git a/man/diffSplice.Rd b/man/diffSplice.Rd
new file mode 100644
index 0000000..a94e4fa
--- /dev/null
+++ b/man/diffSplice.Rd
@@ -0,0 +1,41 @@
+\name{diffSplice}
+\alias{diffSplice}
+\title{Test for Differential Splicing}
+\description{Given a linear model fit at the exon level, test for differences in exon retention between experimental conditions.}
+\usage{
+diffSplice(fit, geneid, exonid=NULL, verbose=TRUE)
+}
+\arguments{
+ \item{fit}{an \code{MArrayLM} fitted model object produced by \code{lmFit} or \code{contrasts.fit}. Rows should correspond to exons.}
+ \item{geneid}{gene identifiers. Either a vector of length \code{nrow(fit)} or the name of the column of \code{fit$genes} containing the gene identifiers. Rows with the same ID are assumed to belong to the same gene.}
+ \item{exonid}{exon identifiers. Either a vector of length \code{nrow(fit)} or the name of the column of \code{fit$genes} containing the exon identifiers.}
+ \item{verbose}{logical, if \code{TRUE} some diagnostic information about the number of genes and exons is output.}
+}
+\value{
+An object of class \code{MArrayLM} containing both exon level and gene level tests.
+Results are sorted by geneid and by exonid within gene.
+ \item{coefficients}{numeric matrix of coefficients of same dimensions as \code{fit}. Each coefficient is the difference between the log-fold-change for that exon versus the average log-fold-change for all other exons for the same gene.}
+ \item{t}{numeric matrix of moderated t-statistics}
+ \item{p.value}{numeric vector of p-values corresponding to the t-statistics}
+ \item{genes}{data.frame of exon annotation}
+ \item{genecolname}{character string giving the name of the column of \code{genes} containing gene IDs}
+ \item{gene.F}{numeric vector of moderated F-statistics, one for each gene.}
+ \item{gene.F.p.value}{numeric vector of p-values corresponding to \code{gene.F}}
+}
+\details{
+This function is often used in conjunction with voom.
+}
+\seealso{
+\code{\link{voom}}
+}
+\author{Gordon Smyth and Charity Law}
+
+\examples{
+\dontrun{
+v <- voom(dge,design)
+fit <- lmFit(v,design)
+ex <- diffSplice(fit,geneid="EntrezID")
+topSplice(ex)
+plotSplice(ex)
+}
+}
diff --git a/man/ebayes.Rd b/man/ebayes.Rd
index 28130d3..17a3e16 100755
--- a/man/ebayes.Rd
+++ b/man/ebayes.Rd
@@ -5,8 +5,10 @@
\title{Empirical Bayes Statistics for Differential Expression}
\description{Given a microarray linear model fit, compute moderated t-statistics, moderated F-statistic, and log-odds of differential expression by empirical Bayes moderation of the standard errors towards a common value.}
\usage{
-ebayes(fit, proportion=0.01, stdev.coef.lim=c(0.1,4), trend=FALSE, robust=FALSE, winsor.tail.p=c(0.05,0.1))
-eBayes(fit, proportion=0.01, stdev.coef.lim=c(0.1,4), trend=FALSE, robust=FALSE, winsor.tail.p=c(0.05,0.1))
+ebayes(fit, proportion=0.01, stdev.coef.lim=c(0.1,4),
+ trend=FALSE, robust=FALSE, winsor.tail.p=c(0.05,0.1))
+eBayes(fit, proportion=0.01, stdev.coef.lim=c(0.1,4),
+ trend=FALSE, robust=FALSE, winsor.tail.p=c(0.05,0.1))
treat(fit, lfc=0, trend=FALSE, robust=FALSE, winsor.tail.p=c(0.05,0.1))
}
\arguments{
diff --git a/man/genas.Rd b/man/genas.Rd
index 5ceed38..c1efe81 100644
--- a/man/genas.Rd
+++ b/man/genas.Rd
@@ -43,7 +43,7 @@ The biological and technical correlations are overlaid on the scatterplot using
}
\value{
- \code{genas} produces a list with the following components.
+ \code{genas} produces a list with the following components:
\item{technical.correlation}{estimate of the technical correlation}
\item{biological.correlation}{estimate of the biological correlation}
\item{covariance.matrix}{estimate of the covariance matrix from which the biological correlation is obtained}
diff --git a/man/geneSetTest.Rd b/man/geneSetTest.Rd
index f743e5e..25b5ad5 100755
--- a/man/geneSetTest.Rd
+++ b/man/geneSetTest.Rd
@@ -7,7 +7,8 @@ Test whether a set of genes is highly ranked relative to other genes in terms of
Genes are assumed to be independent.
}
\usage{
-geneSetTest(index, statistics, alternative="mixed", type="auto", ranks.only=TRUE, nsim=9999)
+geneSetTest(index, statistics, alternative = "mixed", type= "auto",
+ ranks.only = TRUE, nsim=9999)
wilcoxGST(index, statistics, ...)
}
@@ -78,7 +79,7 @@ For this reason, the alternative \code{camera} is now recommended over \code{gen
\seealso{
\code{\link{camera}}, \code{\link{roast}}, \code{\link{romer}}, \code{\link{wilcox.test}}, \code{\link{barcodeplot}}
-An overview of tests in limma is given in \link{08.Tests}.
+There is a topic page on \link{10.GeneSetTests}.
}
\author{Gordon Smyth and Di Wu}
@@ -103,5 +104,3 @@ stat <- rnorm(100)
sel <- 1:10; stat[sel] <- stat[sel]+1
wilcoxGST(sel,stat)
}
-
-\keyword{htest}
diff --git a/man/heatdiagram.Rd b/man/heatdiagram.Rd
index 2588018..0f39551 100755
--- a/man/heatdiagram.Rd
+++ b/man/heatdiagram.Rd
@@ -6,8 +6,12 @@
Creates a heat diagram showing the co-regulation of genes under one condition with a range of other conditions.
}
\usage{
-heatDiagram(results,coef,primary=1,names=NULL,treatments=colnames(coef),limit=NULL,orientation="landscape",low="green",high="red",cex=1,mar=NULL,ncolors=123,...)
-heatdiagram(stat,coef,primary=1,names=NULL,treatments=colnames(stat),critical.primary=4,critical.other=3,limit=NULL,orientation="landscape",low="green",high="red",cex=1,mar=NULL,ncolors=123,...)
+heatDiagram(results, coef, primary=1, names=NULL, treatments=colnames(coef), limit=NULL,
+ orientation="landscape", low="green", high="red", cex=1, mar=NULL,
+ ncolors=123, ...)
+heatdiagram(stat, coef, primary=1, names=NULL, treatments=colnames(stat),
+ critical.primary=4, critical.other=3, limit=NULL, orientation="landscape",
+ low="green", high="red", cex=1, mar=NULL, ncolors=123, ...)
}
\arguments{
\item{results}{\code{TestResults} matrix, containing elements -1, 0 or 1, from \code{\link{decideTests}}}
diff --git a/man/imageplot3by2.Rd b/man/imageplot3by2.Rd
index 70987e1..9902511 100755
--- a/man/imageplot3by2.Rd
+++ b/man/imageplot3by2.Rd
@@ -5,7 +5,8 @@
Write imageplots to files in PNG format, six plots to a file in a 3 by 2 grid arrangement.
}
\usage{
-imageplot3by2(RG, z="Gb", prefix=paste("image",z,sep="-"), path=NULL, zlim=NULL, common.lim=TRUE, ...)
+imageplot3by2(RG, z="Gb", prefix=paste("image",z,sep="-"), path=NULL,
+ zlim=NULL, common.lim=TRUE, ...)
}
\arguments{
\item{RG}{an \code{RGList} or \code{MAList} object, or any list with component named by \code{z}}
diff --git a/man/kooperberg.Rd b/man/kooperberg.Rd
index 83658da..3b09b20 100755
--- a/man/kooperberg.Rd
+++ b/man/kooperberg.Rd
@@ -56,8 +56,10 @@ A comparison of background correction methods for two-colour microarrays.
# given GenePix Results (gpr) files in the working directory (data not
# provided).
\dontrun{
-genepixFiles <- dir(pattern="*\\\\.gpr$") # get the names of the GenePix image analysis output files in the current directory
-RG <- read.maimages(genepixFiles, source="genepix", other.columns=c("F635 SD","B635 SD","F532 SD","B532 SD","B532 Mean","B635 Mean","F Pixels","B Pixels"))
+# get the names of the GenePix image analysis output files in the current directory
+genepixFiles <- dir(pattern="*\\\\.gpr$")
+RG <- read.maimages(genepixFiles, source="genepix", other.columns=c("F635 SD","B635 SD",
+ "F532 SD","B532 SD","B532 Mean","B635 Mean","F Pixels","B Pixels"))
RGmodel <- kooperberg(RG)
MA <- normalizeWithinArrays(RGmodel)
}
diff --git a/man/lmFit.Rd b/man/lmFit.Rd
index 8a181d2..d8fb620 100755
--- a/man/lmFit.Rd
+++ b/man/lmFit.Rd
@@ -3,17 +3,18 @@
\title{Linear Model for Series of Arrays}
\description{Fit linear model for each gene given a series of arrays}
\usage{
-lmFit(object,design=NULL,ndups=1,spacing=1,block=NULL,correlation,weights=NULL,method="ls",...)
+lmFit(object, design=NULL, ndups=1, spacing=1, block=NULL, correlation, weights=NULL,
+ method="ls", ...)
}
\arguments{
\item{object}{object of class \code{numeric}, \code{matrix}, \code{MAList}, \code{EList}, \code{marrayNorm}, \code{ExpressionSet} or \code{PLMset} containing log-ratios or log-values of expression for a series of microarrays}
\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{ndups}{positive integer giving the number of times each gene is printed on an array}
- \item{spacing}{positive integer giving the spacing between duplicate spots, \code{spacing=1} for consecutive spots}
+ \item{ndups}{positive integer giving the number of times each distinct probe is printed on each array.}
+ \item{spacing}{positive integer giving the spacing between duplicate occurrences of the same probe, \code{spacing=1} for consecutive rows.}
\item{block}{vector or factor specifying a blocking variable on the arrays. Has length equal to the number of arrays. Must be \code{NULL} if \code{ndups>2}.}
\item{correlation}{the inter-duplicate or inter-technical replicate correlation}
- \item{weights}{optional numeric matrix containing weights for each spot}
- \item{method}{character string, \code{"ls"} for least squares or \code{"robust"} for robust regression}
+ \item{weights}{non-negative observation weights. Can be a numeric matrix of individual weights, of same size as the object expression matrix, or a numeric vector of array weights with length equal to \code{ncol} of the expression matrix, or a numeric vector of gene weights with length equal to \code{nrow} of the expression matrix.}
+ \item{method}{fitting method; \code{"ls"} for least squares or \code{"robust"} for robust regression}
\item{...}{other optional arguments to be passed to \code{lm.series}, \code{gls.series} or \code{mrlm}}
}
@@ -98,6 +99,9 @@ topTreat(fit2,coef=2)
# Volcano plot
volcanoplot(fit,coef=2,highlight=2)
+# MA plot
+plot(fit,coef=2)
+
# Q-Q plot of moderated t-statistics
qqt(fit$t[,2],df=fit$df.residual+fit$df.prior)
abline(0,1)
diff --git a/man/loessfit.Rd b/man/loessfit.Rd
index 4301a47..65e0dbb 100755
--- a/man/loessfit.Rd
+++ b/man/loessfit.Rd
@@ -8,38 +8,44 @@ Returns fitted values and residuals.
}
\usage{
-loessFit(y, x, weights=NULL, span=0.3, iterations=4, min.weight=1e-5, max.weight=1e5, equal.weights.as.null=TRUE)
+loessFit(y, x, weights=NULL, span=0.3, iterations=4L, min.weight=1e-5, max.weight=1e5,
+ equal.weights.as.null=TRUE, method="weightedLowess")
}
\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{weights}{numeric vector of non-negative prior 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{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}{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.}
+ \item{method}{method used for weighted lowess. Possibilities are \code{"weightedLowess"}, \code{"loess"} or \code{"locfit"}.}
}
\details{
-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.
+This function is essentially a wrapper function for \code{lowess} and \code{weightedLowess} with added error checking.
+The idea is to provide the classic univariate lowess algorithm 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.
+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 prior 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).
+By default, \code{weights} that are all equal (even all zero) are treated as if they were \code{NULL}, so \code{lowess} is called in this case also.
-By default, \code{weights} that are all equal (even all zero) are treated as if they were \code{NULL}.
+When unequal \code{weights} are provided, this function calls \code{weightedLowess} by default, although two other possibilities are also provided.
+\code{weightedLowess} implements a similar algorithm to \code{lowess} except that it uses the prior weights both in the local regressions and in determining which other observations to include in the local neighbourhood of each observation.
-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.
+Two alternative algorithms for weighted lowess curve fitting are provided as options.
+If \code{method="loess"}, then a call is made to \code{loess(y~x,weights=weights,span=span,degree=1,family="symmetric",...)}.
+This method differs from \code{weightedLowess} in that the prior weights are ignored when determining the neighbourhood of each observation.
-The arguments \code{span} and \code{iterations} here have the same meaning as for \code{loess}.
+If \code{method="locfit"}, then repeated calls are made to \code{locfit:::locfit.raw} with \code{deg=1}.
+In principle, this is similar to \code{"loess"}, but \code{"locfit"} makes some approximations and is very much faster and uses much less memory than \code{"loess"} for long data vectors.
+
+The arguments \code{span} and \code{iterations} here have the same meaning as for \code{weightedLowess} and \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.
@@ -47,6 +53,12 @@ When there are insufficient observations to estimate the loess curve, \code{loes
This mimics the behavior of \code{lowess} but not that of \code{loess} or \code{locfit.raw}.
}
+\note{
+With unequal weights, \code{"loess"} was the default method prior to limma version 3.17.25.
+The default was changed to \code{"locfit"} in limma 3.17.25, and then to \code{"weightedLowess"} in limma 3.19.16.
+\code{"weightedLowess"} will potentially give somewhat different results to the older algorithms because the local neighbourhood of each observation is determined differently (more carefully).
+}
+
\value{
A list with components
\item{fitted}{numeric vector of same length as \code{y} giving the loess fit}
@@ -62,7 +74,8 @@ Robust locally weighted regression and smoothing scatterplots.
}
\seealso{
-This function uses \code{\link{lowess}} and \code{\link[locfit]{locfit.raw}}.
+If \code{weights=NULL}, this function calls \code{\link{lowess}}.
+Otherwise it calls \code{\link{weightedLowess}}, \code{\link[locfit]{locfit.raw}} or \code{\link{loess}}.
See the help pages of those functions for references and credits.
Compare with \code{\link{loess}} in the stats package.
diff --git a/man/nec.Rd b/man/nec.Rd
index 57cf5a1..d746d71 100644
--- a/man/nec.Rd
+++ b/man/nec.Rd
@@ -16,21 +16,21 @@ neqc(x, status=NULL, negctrl="negative", regular="regular", offset=16,
\item{regular}{character string identifier for regular probes, i.e., all probes other than control probes.}
\item{offset}{numeric value added to the intensities after background correction.}
\item{robust}{logical. Should robust estimators be used for the background mean and standard deviation?}
- \item{detection.p}{a character string giving the name of the component which contains detection p value information in \code{x} or a numeric matrix giving detection p values, \code{Detection} by default}
+ \item{detection.p}{dection p-values. Only used when no negative control probes can be found in the data. Can be a numeric matrix or a character string giving the name of the component of \code{x$other} containing the matrix.}
\item{...}{any other arguments are passed to \code{normalizeBetweenArrays.}}
}
\details{
-\code{neqc} performs background correction followed by quantile normalization, using negative control probes for background correction and both negative and positive controls for normalization.
+\code{neqc} performs background correction followed by quantile normalization, using negative control probes for background correction and both negative and positive controls for normalization (Shi et al, 2010).
\code{nec} is similar but performs background correction only.
When control data are available, these function call \code{\link{normexp.fit.control}} to estimate the parameters required by normal+exponential(normexp) convolution model with the help of negative control probes, followed by \code{\link{normexp.signal}} to perform the background correction.
If \code{x} contains background intensities \code{x$Eb}, then these are first subtracted from the foreground intensities, prior to normexp background correction.
After background correction, an \code{offset} is added to the data.
-When control data are not available, these functions call \code{\link{normexp.fit.detection.p}} to estimate the normexp parameters.
-\code{\link{normexp.fit.detection.p}} infers negative control probe intensities from regular probes by taking advantage of their detection p value information.
+When expression values for negative controls are not available, the \code{detection.p} argument is used instead,
+In that case, these functions call \code{\link{normexp.fit.detection.p}}, which infers the negative control probe intensities from the detection p-values associated with the regular probes.
-For more descriptions to parameters \code{x}, \code{status}, \code{negctrl}, \code{regular} and \code{detection.p}, please refer to functions \code{\link{normexp.fit.control}}, \code{\link{normexp.fit.detection.p}} and \code{\link{read.ilmn}}.
+For more detailed descriptions of the arguments \code{x}, \code{status}, \code{negctrl}, \code{regular} and \code{detection.p}, please refer to functions \code{\link{normexp.fit.control}}, \code{\link{normexp.fit.detection.p}} and \code{\link{read.ilmn}}.
Both \code{nec} and \code{neqc} perform the above steps.
\code{neqc} continues on to quantile normalize the background-corrected intensities, including control probes.
diff --git a/man/normalizeWithinArrays.Rd b/man/normalizeWithinArrays.Rd
index 64b89a9..7bb879e 100755
--- a/man/normalizeWithinArrays.Rd
+++ b/man/normalizeWithinArrays.Rd
@@ -7,7 +7,9 @@
Normalize the expression log-ratios for one or more two-colour spotted microarray experiments so that the log-ratios average to zero within each array or sub-array.
}
\usage{
-normalizeWithinArrays(object, layout, method="printtiploess", weights=object$weights, span=0.3, iterations=4, controlspots=NULL, df=5, robust="M", bc.method="subtract", offset=0)
+normalizeWithinArrays(object, layout, method="printtiploess", weights=object$weights,
+ span=0.3, iterations=4, controlspots=NULL, df=5, robust="M",
+ bc.method="subtract", offset=0)
MA.RG(object, bc.method="subtract", offset=0)
RG.MA(object)
}
diff --git a/man/normalizeprintorder.Rd b/man/normalizeprintorder.Rd
index dbd290a..a8b5bce 100755
--- a/man/normalizeprintorder.Rd
+++ b/man/normalizeprintorder.Rd
@@ -8,9 +8,12 @@
Normalize intensity values on one or more spotted microarrays to adjust for print-order effects.
}
\usage{
-normalizeForPrintorder(object, layout, start="topleft", method = "loess", separate.channels = FALSE, span = 0.1, plate.size = 32)
-normalizeForPrintorder.rg(R, G, printorder, method = "loess", separate.channels = FALSE, span = 0.1, plate.size = 32, plot = FALSE)
-plotPrintorder(object, layout, start="topleft", slide = 1, method = "loess", separate.channels = FALSE, span = 0.1, plate.size = 32)
+normalizeForPrintorder(object, layout, start="topleft", method = "loess",
+ separate.channels = FALSE, span = 0.1, plate.size = 32)
+normalizeForPrintorder.rg(R, G, printorder, method = "loess", separate.channels = FALSE,
+ span = 0.1, plate.size = 32, plot = FALSE)
+plotPrintorder(object, layout, start="topleft", slide = 1, method = "loess",
+ separate.channels = FALSE, span = 0.1, plate.size = 32)
}
\arguments{
\item{object}{an \code{RGList} or \code{list} object containing components \code{R} and \code{G} which are matrices containing the red and green channel intensities for a series of arrays}
diff --git a/man/plotma.Rd b/man/plot.Rd
old mode 100755
new mode 100644
similarity index 60%
copy from man/plotma.Rd
copy to man/plot.Rd
index 3ab86b2..3e7b953
--- a/man/plotma.Rd
+++ b/man/plot.Rd
@@ -1,21 +1,26 @@
-\title{MA-Plot}
-\name{plotMA}
-\alias{plotMA}
-\alias{plotMA.RGList}
-\alias{plotMA.MAList}
-\alias{plotMA.EList}
-\alias{plotMA.MArrayLM}
-\alias{plotMA.default}
+\title{plot.MArrayLM}
+\name{plot.MArrayLM}
+\alias{plot.RGList}
+\alias{plot.MAList}
+\alias{plot.EList}
+\alias{plot.MArrayLM}
\description{
-Creates an MA-plot with color coding for control spots.
+Plot methods for data or model fit objects.
+These represent the object by creating an MA-plot, with optional color coding for control spots.
}
\usage{
-\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, ...)
+\method{plot}{MAList}(x, y, array = 1, xlab= "A", ylab= "M", main = colnames(x)[array],
+ xlim = NULL, ylim = NULL, status, values, pch, col, cex, legend = TRUE,
+ zero.weights = FALSE, \dots)
+\method{plot}{MArrayLM}(x, y, coef = ncol(x), xlab="AveExpr", ylab="logFC",
+ main = colnames(x)[coef], xlim=NULL, ylim=NULL, status, values, pch, col, cex,
+ legend = TRUE, zero.weights = FALSE, \dots)
}
\arguments{
- \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{x}{an \code{RGList}, \code{MAList}, \code{EList} or \code{MArrayLM} object.}
+ \item{y}{not used.}
+ \item{array}{integer giving the array to be plotted (if \code{x} is an \code{RGList}, \code{MAList} or \code{EList} object).}
+ \item{coef}{integer giving the linear model coefficient to be plotted.}
\item{xlab}{character string giving label for x-axis}
\item{ylab}{character string giving label for y-axis}
\item{main}{character string giving title for plot}
@@ -34,15 +39,15 @@ Creates an MA-plot with color coding for control spots.
Ignored if there is no \code{status} vector.}
\item{legend}{logical, should a legend of plotting symbols and colors be included. Can also be a character string giving position to place legend. Ignored if there is no \code{status} vector.}
\item{zero.weights}{logical, should spots with zero or negative weights be plotted?}
- \item{...}{any other arguments are passed to \code{plot}}
+ \item{\dots}{any other arguments are passed to \code{plot.default}}
}
\details{
An MA-plot is a plot of log-intensity ratios (M-values) versus log-intensity averages (A-values).
-If \code{MA} is an \code{RGList} or \code{MAList} then this function produces an ordinary within-array MA-plot.
-If \code{MA} is an \code{MArrayLM} object, then the plot is an fitted model MA-plot in which the estimated coefficient is on the y-axis and the average A-value is on the x-axis.
+If \code{x} is an \code{RGList} or \code{MAList} then this function produces an ordinary within-array MA-plot.
+If \code{x} is an \code{MArrayLM} object, then the plot is an fitted model MA-plot in which the estimated coefficient is on the y-axis and the average A-value is on the x-axis.
-If \code{MA} is a \code{matrix} or \code{ExpressionSet} object, then this function produces a between-array MA-plot.
+If \code{x} is a \code{EList} object, then this function produces a between-array MA-plot.
In this case the A-values in the plot are the average log-intensities across the arrays and the M-values are the deviations of the log-intensities for the specified array from the average.
If there are more than five arays, then the average is computed robustly using medians.
With five or fewer arrays, it is computed by means.
@@ -51,7 +56,7 @@ The \code{status} vector is intended to specify the control status of each spot,
The vector is usually computed using the function \code{\link{controlStatus}} and a spot-types file.
However the function may be used to highlight any subset of spots.
-The \code{status} can be included as the component \code{MA$genes$Status} instead of being passed as an argument to \code{plotMA}.
+The \code{status} can be included as the component \code{x$genes$Status} instead of being passed as an argument to \code{plot}.
The arguments \code{values}, \code{pch}, \code{col} and \code{cex} can be included as attributes to \code{status} instead of being passed as arguments to \code{plotMA}.
See \code{\link[graphics]{points}} for possible values for \code{pch}, \code{col} and \code{cex}.
@@ -61,6 +66,8 @@ See \code{\link[graphics]{points}} for possible values for \code{pch}, \code{col
\references{See \url{http://www.statsci.org/micrarra/refs/maplots.html}}
\author{Gordon Smyth}
\examples{
+# See also lmFit example
+
MA <- new("MAList")
MA$A <- runif(300,4,16)
MA$M <- rt(300,df=3)
@@ -71,18 +78,20 @@ status[4:6] <- "M=3"
MA$M[4:6] <- 3
status[7:9] <- "M=-3"
MA$M[7:9] <- -3
-plotMA(MA,main="MA-Plot with Simulated Data",status=status,values=c("M=0","M=3","M=-3"),col=c("blue","red","green"))
+plot(MA,main="MA-Plot with Simulated Data",status=status,
+ values=c("M=0","M=3","M=-3"),col=c("blue","red","green"))
# Same as above
attr(status,"values") <- c("M=0","M=3","M=-3")
attr(status,"col") <- c("blue","red","green")
-plotMA(MA,main="MA-Plot with Simulated Data",status=status)
+plot(MA,main="MA-Plot with Simulated Data",status=status)
# Same as above
MA$genes$Status <- status
-plotMA(MA,main="MA-Plot with Simulated Data")
+plot(MA,main="MA-Plot with Simulated Data")
}
\seealso{
-An overview of diagnostic functions available in LIMMA is given in \link{09.Diagnostics}.
+\code{\link[limma:plotma]{plotMA}}, \code{\link{plotFB}}, \code{\link{plotMDS}}, \code{\link{plotSA}}
+
+An overview of diagnostic plots available in LIMMA is given in \link{09.Diagnostics}.
}
-\keyword{hplot}
diff --git a/man/plotMDS.Rd b/man/plotMDS.Rd
index 258fbf8..d6a7062 100644
--- a/man/plotMDS.Rd
+++ b/man/plotMDS.Rd
@@ -5,20 +5,25 @@
\alias{plotMDS.default}
\alias{MDS-class}
\alias{show,MDS-method}
+
\description{
-Plot the sample relations based on MDS.
+Plot the sample relations based on multi-dimensional scaling (MDS).
Distances on the plot can be interpreted in terms of \emph{leading log2-fold-change}.
}
-\usage{
-\method{plotMDS}{default}(x, top=500, labels=colnames(x), col=NULL, cex=1, dim.plot=c(1,2), ndim=max(dim.plot), gene.selection="pairwise",
- xlab=paste("Dimension",dim.plot[1]), ylab=paste("Dimension",dim.plot[2]), ...)
-\method{plotMDS}{MDS}(x, labels=colnames(x$distance.matrix), col=NULL, cex=1, dim.plot=x$dim.plot, xlab=paste("Dimension",dim.plot[1]), ylab=paste("Dimension",dim.plot[2]),...)
+\usage{
+\method{plotMDS}{default}(x, top=500, labels=NULL, pch=NULL, col=NULL, cex=1,
+ dim.plot=c(1,2), ndim=max(dim.plot), gene.selection="pairwise",
+ xlab=paste("Dimension",dim.plot[1]), ylab=paste("Dimension",dim.plot[2]), ...)
+\method{plotMDS}{MDS}(x, labels=NULL, pch=NULL, col=NULL, cex=1, dim.plot=x$dim.plot,
+ xlab=paste("Dimension",dim.plot[1]), ylab=paste("Dimension",dim.plot[2]),...)
}
+
\arguments{
\item{x}{any data object which can be coerced to a matrix, such as \code{ExpressionSet} or \code{EList}.}
\item{top}{number of top genes used to calculate pairwise distances.}
\item{labels}{character vector of sample names or labels. If \code{x} has no column names, then defaults the index of the samples.}
+ \item{pch}{plotting symbol or symbols. See \code{\link{points}} for possible values. Ignored if \code{labels} is non-\code{NULL}.}
\item{col}{numeric or character vector of colors for the plotting characters.}
\item{cex}{numeric vector of plot symbol expansions.}
\item{dim.plot}{which two dimensions should be plotted, numeric vector of length two.}
@@ -76,5 +81,3 @@ mds <- plotMDS(x, col=c(rep("black",3), rep("red",3)) )
# or labels can be provided, here group indicators:
plotMDS(mds, col=c(rep("black",3), rep("red",3)), labels= c(rep("Grp1",3), rep("Grp2",3)))
}
-
-\keyword{hplot}
diff --git a/man/plotSA.Rd b/man/plotSA.Rd
index cd855df..f75a3be 100755
--- a/man/plotSA.Rd
+++ b/man/plotSA.Rd
@@ -5,7 +5,8 @@
Plot log residual standard deviation versus average log expression for a fitted microarray linear model.
}
\usage{
-plotSA(fit, xlab="Average log-expression", ylab="log2(sigma)", zero.weights=FALSE, pch=16, cex=0.2, ...)
+plotSA(fit, xlab="Average log-expression", ylab="log2(sigma)",
+ zero.weights=FALSE, pch=16, cex=0.2, ...)
}
\arguments{
\item{fit}{an \code{MArrayLM} object.}
diff --git a/man/plotSplice.Rd b/man/plotSplice.Rd
new file mode 100644
index 0000000..7be592e
--- /dev/null
+++ b/man/plotSplice.Rd
@@ -0,0 +1,29 @@
+\title{Plot exons on differentially spliced gene}
+\name{plotSplice}
+\alias{plotSplice}
+\description{
+Plot exons of differentially spliced gene.
+}
+\usage{
+plotSplice(fit, coef=ncol(fit), geneid=NULL, rank=1L, FDR = 0.05)
+}
+\arguments{
+ \item{fit}{\code{MArrayLM} fit object produced by \code{diffSplice}.}
+ \item{coef}{the coefficient (column) of fit for which differentially splicing is assessed.}
+ \item{geneid}{character string, ID of the gene to plot.}
+ \item{rank}{integer, if \code{geneid=NULL} then this ranked gene will be plotted.}
+ \item{FDR}{numeric, mark exons with false discovery rate less than this cutoff.}
+}
+
+\details{
+Plots interaction log2-fold-change by exon for the specified gene.
+}
+
+\value{A plot is created on the current graphics device.}
+\author{Gordon Smyth and Yifang Hu}
+\seealso{
+\code{\link{diffSplice}}
+
+An overview of diagnostic functions available in LIMMA is given in \link{09.Diagnostics}.
+}
+\examples{# See diffSplice}
diff --git a/man/plotma.Rd b/man/plotma.Rd
index 3ab86b2..f8924d2 100755
--- a/man/plotma.Rd
+++ b/man/plotma.Rd
@@ -10,12 +10,18 @@
Creates an MA-plot with color coding for control spots.
}
\usage{
-\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, ...)
+\method{plotMA}{default}(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}{MArrayLM}(MA, coef = ncol(MA), xlab = "AveExpr", ylab = "logFC",
+ main = colnames(MA)[coef], xlim = NULL, ylim = NULL, status, values, pch, col,
+ cex, legend = TRUE, zero.weights = FALSE, ...)
}
\arguments{
\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{array}{integer giving the array to be plotted.}
+ \item{coef}{integer giving the linear model coefficient to be plotted.}
\item{xlab}{character string giving label for x-axis}
\item{ylab}{character string giving label for y-axis}
\item{main}{character string giving title for plot}
@@ -71,7 +77,8 @@ status[4:6] <- "M=3"
MA$M[4:6] <- 3
status[7:9] <- "M=-3"
MA$M[7:9] <- -3
-plotMA(MA,main="MA-Plot with Simulated Data",status=status,values=c("M=0","M=3","M=-3"),col=c("blue","red","green"))
+plotMA(MA,main="MA-Plot with Simulated Data",
+ status=status,values=c("M=0","M=3","M=-3"),col=c("blue","red","green"))
# Same as above
attr(status,"values") <- c("M=0","M=3","M=-3")
diff --git a/man/plotma3by2.Rd b/man/plotma3by2.Rd
index 51d22ac..20d0a2d 100755
--- a/man/plotma3by2.Rd
+++ b/man/plotma3by2.Rd
@@ -5,10 +5,11 @@
Write MA-plots to files in PNG format, six plots to a file in a 3 by 2 grid arrangement.
}
\usage{
-plotMA3by2(MA, prefix="MA", path=NULL, main=colnames(MA), zero.weights=FALSE, common.lim=TRUE, device="png", ...)
+plotMA3by2(MA, prefix="MA", path=NULL, main=colnames(MA),
+ zero.weights=FALSE, common.lim=TRUE, device="png", ...)
}
\arguments{
- \item{MA}{an \code{MAList} or \code{RGList} object, or any list with components \code{M} containing log-ratios and \code{A} containing average intensities}
+ \item{MA}{an \code{MAList}, \code{RGList}, \code{EListRaw} or \code{EList} object, or a matrix containing log-intensities.}
\item{prefix}{character string giving prefix to attach to file names}
\item{path}{character string specifying directory for output files}
\item{main}{character vector giving titles for plots}
diff --git a/man/poolvar.Rd b/man/poolvar.Rd
index 8f33a59..7804c89 100755
--- a/man/poolvar.Rd
+++ b/man/poolvar.Rd
@@ -48,10 +48,6 @@ Welch, B. L. (1947). The generalization of 'Student's' problem when several diff
Welch, B. L. (1949). Further note on Mrs. Aspin's tables and on certain approximations to the tabled function. \emph{Biometrika} \bold{36}, 293-296.
}
-\seealso{
-\link{10.Other}
-}
-
\examples{
# Welch's t-test with unequal variances
x <- rnorm(10,mean=1,sd=2)
@@ -63,5 +59,3 @@ tstat <- (mean(x)-mean(y)) / sqrt(out$var*out$multiplier)
pvalue <- 2*pt(-abs(tstat),df=out$df)
# Equivalent to t.test(x,y)
}
-
-\keyword{htest}
diff --git a/man/printtipWeights.Rd b/man/printtipWeights.Rd
index ea099a3..b013f69 100755
--- a/man/printtipWeights.Rd
+++ b/man/printtipWeights.Rd
@@ -2,14 +2,16 @@
\alias{printtipWeights}
\alias{printtipWeightsSimple}
\title{Sub-array Quality Weights}
+
\description{
-Estimates relative quality weights for each sub-array in a multi-array
-experiment.
+Estimates relative quality weights for each sub-array in a multi-array experiment.
}
+
\usage{
-printtipWeights(object, design = NULL, weights = NULL, method = "genebygene", layout, maxiter = 50, tol = 1e-10, trace=FALSE)
-
+printtipWeights(object, design = NULL, weights = NULL, method = "genebygene", layout,
+ maxiter = 50, tol = 1e-10, trace=FALSE)
}
+
\arguments{
\item{object}{object of class \code{numeric}, \code{matrix}, \code{MAList}, \code{marrayNorm},
or \code{ExpressionSet} containing log-ratios or log-values of
@@ -24,17 +26,17 @@ printtipWeights(object, design = NULL, weights = NULL, method = "genebygene", la
\item{layout}{list specifying the dimensions of the spot matrix and the grid matrix. For details see \code{\link[limma:PrintLayout]{PrintLayout-class}}.}
\item{maxiter}{maximum number of iterations allowed.}
\item{tol}{convergence tolerance.}
- \item{trace}{logical variable. If true then output diagnostic information
- at each iteration of '"reml"' algorithm.}
- }
+ \item{trace}{logical variable. If true then output diagnostic information at each iteration of \code{"reml"} algorithm.}
+}
+
\details{
The relative reliability of each sub-array (print-tip group) is estimated by measuring how
well the expression values for that sub-array follow the linear model.
The method described in Ritchie et al (2006) and implemented in
the \code{arrayWeights} function is adapted for this purpose.
A heteroscedastic model is fitted to the expression values for
-each gene by calling the function \code{lm.wfit}. The dispersion model
-is fitted to the squared residuals from the mean fit, and is set up to
+each gene by calling the function \code{lm.wfit}.
+The dispersion model is fitted to the squared residuals from the mean fit, and is set up to
have sub-array specific coefficients, which are updated in either full REML
scoring iterations, or using an efficient gene-by-gene update algorithm.
The final estimates of the sub-array variances are converted to weights.
@@ -44,17 +46,21 @@ In particular, the arguments \code{design}, \code{weights} and \code{layout} wil
be extracted from the data \code{object} if available and do not normally need to
be set explicitly in the call; if any of these are set in the call then they
will over-ride the slots or components in the data \code{object}.
-
}
-\value{
- A matrix of sub-array weights which can be passed to \code{lmFit}.
- }
+
+\value{A matrix of sub-array weights.}
+
\references{
-Ritchie, M. E., Diyagama, D., Neilson, van Laar, R., J., Dobrovic, A., Holloway, A., and Smyth, G. K. (2006). Empirical array quality weights in the analysis of microarray data. BMC Bioinformatics 7, 261. \url{http://www.biomedcentral.com/1471-2105/7/261/abstract}
+Ritchie, M. E., Diyagama, D., Neilson, van Laar, R., J., Dobrovic, A., Holloway, A., and Smyth, G. K. (2006).
+Empirical array quality weights in the analysis of microarray data.
+\emph{BMC Bioinformatics} 7, 261.
+\url{http://www.biomedcentral.com/1471-2105/7/261/abstract}
}
+
\seealso{
An overview of linear model functions in limma is given by \link{06.LinearModels}.
}
+
\examples{
\dontrun{
# This example is designed for work on a subset of the data
diff --git a/man/qqt.Rd b/man/qqt.Rd
index 6457700..cde8fd5 100755
--- a/man/qqt.Rd
+++ b/man/qqt.Rd
@@ -4,9 +4,9 @@
\title{Student's t Quantile-Quantile Plot}
\description{Plots the quantiles of a data sample against the theoretical quantiles of a Student's t distribution.}
\usage{qqt(y, df = Inf, ylim = range(y), main = "Student's t Q-Q Plot",
- xlab = "Theoretical Quantiles", ylab = "Sample Quantiles", plot.it = TRUE, ...)
+ xlab = "Theoretical Quantiles", ylab = "Sample Quantiles", plot.it = TRUE, \dots)
qqf(y, df1, df2, ylim=range(y), main= "F Distribution Q-Q Plot",
- xlab = "Theoretical Quantiles", ylab = "Sample Quantiles", plot.it = TRUE,...)
+ xlab = "Theoretical Quantiles", ylab = "Sample Quantiles", plot.it = TRUE, \dots)
}
\arguments{
\item{y}{a numeric vector or array containing the data sample}
@@ -18,7 +18,7 @@ qqf(y, df1, df2, ylim=range(y), main= "F Distribution Q-Q Plot",
\item{xlab}{x-axis title for the plot}
\item{ylab}{y-axis title for the plot}
\item{plot.it}{whether or not to produce a plot}
-\item{...}{other arguments to be passed to \code{plot}}
+\item{\dots}{other arguments to be passed to \code{plot}}
}
\value{A list is invisibly returned containing the values plotted in the QQ-plot:
\item{x}{theoretical quantiles of the t-distribution or F-distribution}
diff --git a/man/read.columns.Rd b/man/read.columns.Rd
index 22b3e68..16b4b08 100644
--- a/man/read.columns.Rd
+++ b/man/read.columns.Rd
@@ -5,7 +5,8 @@
Reads specified columns from a file in table format and creates a data frame from it, with cases corresponding to lines and variables to fields in the file.
}
\usage{
-read.columns(file,required.col=NULL,text.to.search="",sep="\t",quote="\"",skip=0,fill=TRUE,blank.lines.skip=TRUE,comment.char="",allowEscapes=FALSE,...)
+read.columns(file, required.col=NULL, text.to.search="", sep="\t", quote="\"", skip=0,
+ fill=TRUE, blank.lines.skip=TRUE, comment.char="", allowEscapes=FALSE, ...)
}
\arguments{
\item{file}{the name of the file which the data are to be read from.}
@@ -46,4 +47,5 @@ A data frame (data.frame) containing a representation of the data in the file.
An overview of LIMMA functions for reading data is given in \link{03.ReadingData}.
}
+\keyword{IO}
\keyword{file}
diff --git a/man/read.idat.Rd b/man/read.idat.Rd
new file mode 100644
index 0000000..7a6f716
--- /dev/null
+++ b/man/read.idat.Rd
@@ -0,0 +1,70 @@
+\name{read.idat}
+\alias{read.idat}
+
+\title{Read Illumina expression data directly from IDAT files}
+
+\description{Read Illumina BeadArray data from IDAT and manifest (.bgx) files for gene expression platforms.}
+
+\usage{
+read.idat(idatfiles, bgxfile, dateinfo=FALSE)
+}
+
+\arguments{
+ \item{idatfiles}{character vector specifying idat files to be read in.}
+ \item{bgxfile}{character string specifying bead manifest file (.bgx) to be read in.}
+ \item{dateinfo}{logical. Should date and software version info be read in?}
+}
+
+\details{
+ Illumina's BeadScan/iScan software ouputs probe intensities in IDAT
+ format (encrypted XML files) and probe info in a platform specific manifest file (.bgx).
+ These files can be processed using the low-level functions \code{readIDAT} and \code{readBGX}
+ from the \code{illuminaio} package (Smith et al. 2013).
+
+ The \code{read.idat} function provides a convenient way to read these files.
+ into R and store them in an \code{EListRaw-class} object (similar to \code{read.ilmn},
+ which imports data output by Illumina's GenomeStudio software) that can be used
+ by downstream processing functions in \code{limma}.
+
+ Probe types are indicated in the \code{Status} column of the \code{genes}
+ component of the \code{EListRaw-class} object.
+}
+
+\value{
+ An \code{EListRaw-class} object with the following components:
+ \item{E}{ numeric matrix of raw intensities.}
+ \item{other}{ list containing matrices of \code{NumBeads} and \code{STDEV} for each probe.}
+ \item{genes}{ data.frame of probe annotation.}
+ \item{targets}{ data.frame of sample information.}
+}
+
+\references{
+ Smith ML, Baggerly KA, Bengtsson H, Ritchie ME, Hansen KD (2013).
+ illuminaio: An open source IDAT parsing tool. \emph{F1000 Research}, 2:264.
+ \url{http://f1000research.com/articles/2-264/v1}
+}
+
+\author{Matt Ritchie}
+
+\seealso{
+ \code{read.ilmn} imports gene expression data output by GenomeStudio.
+
+ \code{neqc} performs normexp by control background correction, log
+ transformation and quantile between-array normalization for
+ Illumina expression data.
+
+ \code{propexpr} estimates the proportion of expressed probes in a microarray.
+}
+
+\examples{
+\dontrun{
+idatfiles = dir(pattern="idat")
+bgxfile = dir(pattern="bgx")
+data = read.idat(idatfiles, bgxfile)
+propexpr(data)
+datanorm = neqc(data)
+}
+}
+\keyword{file}
+\concept{illumina microarrays}
+\concept{microarray data file}
diff --git a/man/read.ilmn.Rd b/man/read.ilmn.Rd
index 4254a58..56b10a5 100644
--- a/man/read.ilmn.Rd
+++ b/man/read.ilmn.Rd
@@ -3,19 +3,19 @@
\title{Read Illumina Expression Data}
\description{Read Illumina summary probe profile files and summary control probe profile files}
\usage{
-read.ilmn(files=NULL, ctrlfiles=NULL, path=NULL, ctrlpath=NULL,
-probeid="Probe", annotation=c("TargetID", "SYMBOL"), expr="AVG_Signal", other.columns="Detection",
-sep="\t", quote="\"", verbose=TRUE, \dots)
+read.ilmn(files=NULL, ctrlfiles=NULL, path=NULL, ctrlpath=NULL, probeid="Probe",
+ annotation=c("TargetID", "SYMBOL"), expr="AVG_Signal",
+ other.columns="Detection", sep="\t", quote="\"", verbose=TRUE, \dots)
}
\arguments{
\item{files}{character vector giving the names of the summary probe profile files.}
\item{ctrlfiles}{character vector giving the names of the summary control probe profile files.}
- \item{path}{character string giving the directory containing the summary probe profile files. The default is the current working directory.}
- \item{ctrlpath}{character string giving the directory containing the summary control probe profile files. Defaults to the same directory as for the probe profile files.}
+ \item{path}{character string giving the directory containing the summary probe profile files. Default is the current working directory.}
+ \item{ctrlpath}{character string giving the directory containing the summary control probe profile files. Default is the same directory as for the probe profile files.}
\item{probeid}{character string giving the name of the probe identifier column.}
- \item{annotation}{character vector giving possible names of the annotation column. It could be called "TargetID" or "SYMBOL" depending on which version of BeadStudio is used.}
- \item{expr}{character string giving the keyword in the names of the expression intensity columns.}
- \item{other.columns}{character vector giving the keywords in the names of extra columns required, such as "Detection", "Avg_NBEADS", "BEAD_STDEV" etc. Each keyword corresponds to one type of columns. The detection p value columns will be read in by default.}
+ \item{annotation}{character vector giving possible column names for probe annotation.}
+ \item{expr}{character string giving a keyword identifying the expression intensity columns. Any input column with column name containing this key will be read as containing intensity values.}
+ \item{other.columns}{character vector giving keywords sufficient to identify any extra data columns that should be read in, such as "Detection", "Avg_NBEADS", "BEAD_STDEV" etc. The default of \code{Detection} is usually sufficient to identify the columns containing detection p-values.}
\item{sep}{the field separator character.}
\item{quote}{character string of characters to be treated as quote marks.}
\item{verbose}{logical, \code{TRUE} to report names of profile files being read.}
@@ -38,10 +38,9 @@ Examples of keywords are "Detection", "Avg_NBEADS", "BEAD_STDEV" etc.
\value{
An \code{\link{EListRaw-class}} object with the following components:
-\item{E}{ numeric matrix of raw intensities.}
-\item{genes}{ data.frame of probe annotation.}
-\item{targets}{ data.frame of sample information.}
-\item{other}{ list of other column data.}
+\item{E}{numeric matrix of intensities.}
+\item{genes}{data.frame of probe annotation. Contains any columns specified by \code{annotation} that are found in the input files.}
+\item{other}{a list of matrices corresponding to any \code{other.columns} found in the input files.}
}
\author{Wei Shi and Gordon K Smyth}
@@ -65,4 +64,6 @@ x <- read.ilmn(files="sample probe profile.txt",
# See neqc and beadCountWeights for other examples using read.ilmn
}
-\keyword{IO}
+\keyword{file}
+\concept{illumina microarrays}
+\concept{microarray data file}
diff --git a/man/read.ilmn.targets.Rd b/man/read.ilmn.targets.Rd
index 404dff4..e547f61 100644
--- a/man/read.ilmn.targets.Rd
+++ b/man/read.ilmn.targets.Rd
@@ -27,4 +27,4 @@ An \code{\link{EListRaw-class}} object. See return value of the function \code{\
\code{\link{read.ilmn}}
}
-\keyword{IO}
+\concept{illumina microarrays}
diff --git a/man/read.maimages.Rd b/man/read.maimages.Rd
index 24d601e..90a820e 100755
--- a/man/read.maimages.Rd
+++ b/man/read.maimages.Rd
@@ -15,6 +15,7 @@ read.imagene(files, path=NULL, ext=NULL, names=NULL, columns=NULL, other.columns
}
\arguments{
\item{files}{character vector giving the names of the files containing image analysis output or, for Imagene data, a character matrix of names of files.
+ Alternatively, it can be a data.frame containing a column called \code{FileName}.
If omitted, then all files with extension \code{ext} in the specified directory will be read in alphabetical order.}
\item{source}{character string specifying the image analysis program which produced the output files. Choices are \code{"generic"}, \code{"agilent"}, \code{"agilent.median"}, \code{"agilent.mean"}, \code{"arrayvision"}, \code{"arrayvision.ARM"}, \code{"arrayvision.MTM"}, \code{"bluefuse"}, \code{"genepix"}, \code{"genepix.custom"}, \code{"genepix.median"}, \code{"imagene"}, \code{"imagene9"}, \code{"quantarray"}, \code{"scanarrayexpress"}, \code{"smd.old"}, \code{"smd"}, \code{"spot"} [...]
\item{path}{character string giving the directory containing the files.
@@ -112,7 +113,7 @@ For two-color data, an \code{\link[limma:rglist]{RGList}} object containing the
\item{weights}{spot quality weights, if \code{wt.fun} is given}
\item{other}{list containing matrices corresponding to \code{other.columns} if given}
\item{genes}{data frame containing annotation information about the probes, for example gene names and IDs and spatial positions on the array, currently set only if \code{source} is \code{"agilent"}, \code{"genepix"} or \code{source="imagene"} or if the \code{annotation} argument is set}
- \item{targets}{data frame with column \code{FileName} giving the names of the files read}
+ \item{targets}{data frame with column \code{FileName} giving the names of the files read. If \code{files} was a data.frame on input, then the whole data.frame is stored here on output.}
\item{source}{character string giving the image analysis program name}
\item{printer}{list of class \code{\link[=PrintLayout-class]{PrintLayout}}, currently set only if \code{source="imagene"}}
}
@@ -141,4 +142,6 @@ RG <- read.maimages(files,"genepix",wt.fun=wtflags(0.1))}
\dontrun{files <- dir(pattern="*\\\\.spot$")
RG <- read.maimages(files,"spot",wt.fun=wtarea(150))}
}
+
\keyword{file}
+\concept{microarray data file}
diff --git a/man/readGPRHeader.Rd b/man/readGPRHeader.Rd
index 947188a..f404fe7 100755
--- a/man/readGPRHeader.Rd
+++ b/man/readGPRHeader.Rd
@@ -2,9 +2,9 @@
\alias{readGenericHeader}
\alias{readGPRHeader}
\alias{readSMDHeader}
-\title{Read Header Information from Image Analysis Raw Data File}
+\title{Read Header Information from Microarray Raw Data File}
\description{
-Read the header information from a GenePix Results (GPR) file or from an SMD raw data file.
+Read the header information from a microarray raw data file, as output from an image analysis software program such as GenePix.
These functions are used internally by \code{read.maimages} and are not usually called directly by users.
}
\usage{
@@ -29,11 +29,9 @@ A key component is \code{NHeaderRecords} which gives the number of lines in the
All other components are character vectors.
}
\references{
-See \url{http://www.axon.com/gn_GenePix_File_Formats.html} for GenePix formats.
+See \url{http://www.moleculardevices.com/Products/Software/GenePix-Pro.html} for GenePix formats.
-See \url{http://www.bluegnome.co.uk} for information on BlueFuse.
-
-See \url{http://genome-www.stanford.edu/Microarray} for the SMD.
+See \url{http://http://smd.princeton.edu} for the SMD.
}
\author{Gordon Smyth}
\seealso{\code{\link{read.maimages}}
@@ -41,3 +39,5 @@ See \url{http://genome-www.stanford.edu/Microarray} for the SMD.
An overview of LIMMA functions to read data is given in \link{03.ReadingData}.
}
\keyword{file}
+\concept{microarray data file}
+
diff --git a/man/removeBatchEffect.Rd b/man/removeBatchEffect.Rd
index aa8e25d..1899149 100644
--- a/man/removeBatchEffect.Rd
+++ b/man/removeBatchEffect.Rd
@@ -5,15 +5,17 @@
Remove batch effects from expression data.
}
\usage{
-removeBatchEffect(x,batch=NULL,batch2=NULL,covariates=NULL,design=matrix(1,ncol(x),1))
+removeBatchEffect(x, batch=NULL, batch2=NULL, covariates=NULL,
+ design=matrix(1,ncol(x),1), \dots)
}
\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. All values must be non-missing.}
+ \item{x}{numeric matrix, or any object that can be coerced to a matrix by \code{as.matrix(x)}, containing log-expression values for a series of samples.
+ Rows correspond to probes and columns to samples.}
\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.}
\item{design}{optional design matrix relating to treatment conditions to be preserved}
+ \item{\dots}{other arguments are passed to \code{\link{lmFit}}.}
}
\value{
A numeric matrix of log-expression values with batch and covariate effects removed.
@@ -26,6 +28,9 @@ For linear modelling, it is better to include the batch factors in the linear mo
The design matrix is used to describe comparisons between the samples, for example treatment effects, which should not be removed.
The function (in effect) fits a linear model to the data, including both batches and regular treatments, then removes the component due to the batch effects.
+
+The data object \code{x} can be of any class for which \code{lmFit} and \code{as.matrix} work.
+If \code{x} contains weights, then these will be used in estimating the batch effects.
}
\seealso{
@@ -33,4 +38,10 @@ The function (in effect) fits a linear model to the data, including both batches
}
\author{Gordon Smyth and Carolyn de Graaf}
-\keyword{regression}
+\examples{
+y <- matrix(rnorm(10*6),10,6)
+colnames(y) <- c("A1","A2","A3","B1","B2","B3")
+y[,1:3] <- y[,1:3] + 10
+y
+removeBatchEffect(y,batch=c("A","A","A","B","B","B"))
+}
diff --git a/man/roast.Rd b/man/roast.Rd
index 3143845..1e97d67 100644
--- a/man/roast.Rd
+++ b/man/roast.Rd
@@ -1,11 +1,7 @@
\name{roast}
\alias{roast}
-\alias{roast.EList}
-\alias{roast.MAList}
\alias{roast.default}
\alias{mroast}
-\alias{mroast.EList}
-\alias{mroast.MAList}
\alias{mroast.default}
\alias{Roast-class}
\alias{show,Roast-method}
@@ -15,12 +11,13 @@ Rotation gene set testing for linear models.
}
\usage{
-\method{roast}{EList}(y, index=NULL, design=NULL, contrast=ncol(design), set.statistic="mean",
+\S3method{roast}{default}(y, index=NULL, design=NULL, contrast=ncol(design), set.statistic="mean",
gene.weights=NULL, array.weights=NULL, weights=NULL, block=NULL, correlation,
- var.prior=NULL, df.prior=NULL, trend.var=FALSE, nrot=999)
-\method{mroast}{EList}(y, index=NULL, design=NULL, contrast=ncol(design), set.statistic="mean",
+ var.prior=NULL, df.prior=NULL, trend.var=FALSE, nrot=999, approx.zscore=TRUE, \dots)
+\S3method{mroast}{default}(y, index=NULL, design=NULL, contrast=ncol(design), set.statistic="mean",
gene.weights=NULL, array.weights=NULL, weights=NULL, block=NULL, correlation,
- var.prior=NULL, df.prior=NULL, trend.var=FALSE, nrot=999, adjust.method="BH", midp=TRUE, sort="directional")
+ var.prior=NULL, df.prior=NULL, trend.var=FALSE, nrot=999, approx.zscore=TRUE,
+ adjust.method="BH", midp=TRUE, sort="directional", \dots)
}
\arguments{
@@ -43,14 +40,16 @@ Rotation gene set testing for linear models.
\item{nrot}{number of rotations used to estimate the p-values.}
\item{adjust.method}{method used to adjust the p-values for multiple testing. See \code{\link{p.adjust}} for possible values.}
\item{midp}{logical, should mid-p-values be used in instead of ordinary p-values when adjusting for multiple testing?}
- \item{sort}{character, whether to sort output table by directional p-values (\code{"directional"}), non-directional p-value (\code{"mixed"}), or not at all (\code{"none"}).}
+ \item{sort}{character, whether to sort output table by directional p-value (\code{"directional"}), non-directional p-value (\code{"mixed"}), or not at all (\code{"none"}).}
+ \item{approx.zscore}{logical, if \code{TRUE} then a fast approximation is used to convert t-statistics into z-scores prior to computing set statistics. If \code{FALSE}, z-scores will be exact.}
+ \item{\dots}{other arguments not currently used.}
}
\value{
\code{roast} produces an object of class \code{"Roast"}.
This consists of a list with the following components:
- \item{p.value}{data.frame with columns \code{Active.Prop} and \code{P.Value}, giving the proportion of genes in the set contributing meaningfully to significance and estimated p-values, respectively.
-Rows correspond to the alternative hypotheses mixed, up or down.}
+ \item{p.value}{data.frame with columns \code{Active.Prop} and \code{P.Value}, giving the proportion of genes in the set contributing materially to significance and estimated p-values, respectively.
+Rows correspond to the alternative hypotheses Down, Up, UpOrDown (two-sided) and Mixed.}
\item{var.prior}{prior value for residual variances.}
\item{df.prior}{prior degrees of freedom for residual variances.}
@@ -113,7 +112,7 @@ The default setting for the set statistic was changed in limma 3.5.9 (3 June 201
\seealso{
\code{\link{camera}}, \code{\link{romer}}, \code{\link{geneSetTest}}, \code{\link{symbols2indices}}.
-An overview of tests in limma is given in \link{08.Tests}.
+There is a topic page on \link{10.GeneSetTests}.
}
\author{Gordon Smyth and Di Wu}
@@ -155,4 +154,3 @@ roast(y,index1,design,contrast=2)
mroast(y,index1,design,contrast=2)
mroast(y,list(set1=index1,set2=index2),design,contrast=2)
}
-\keyword{htest}
diff --git a/man/romer.Rd b/man/romer.Rd
index dcb05ee..fc03f6a 100644
--- a/man/romer.Rd
+++ b/man/romer.Rd
@@ -1,12 +1,12 @@
\name{romer}
\alias{romer}
-\alias{romer2}
\title{Rotation Gene Set Enrichment Analysis}
\description{
Gene set enrichment analysis for linear models using rotation tests (ROtation testing using MEan Ranks).
}
\usage{
-romer(index,y,design,contrast=ncol(design),array.weights=NULL,block=NULL,correlation,set.statistic="mean",nrot=9999)
+romer(index, y, design, contrast=ncol(design), array.weights=NULL, block=NULL,
+ correlation, set.statistic="mean", nrot=9999)
}
\arguments{
\item{index}{list of indices specifying the rows of \code{y} in the gene sets. The list can be made using \link{symbols2indices}.}
@@ -64,7 +64,7 @@ This statistic performs well in practice but is slightly slower to compute.
\code{\link{roast}},
\code{\link{wilcoxGST}}
-An overview of tests in limma is given in \link{08.Tests}.
+There is a topic page on \link{10.GeneSetTests}.
}
\author{Yifang Hu and Gordon Smyth}
@@ -101,4 +101,4 @@ r
topRomer(r,alt="up")
topRomer(r,alt="down")
}
-\keyword{htest}
+\concept{gene set test}
diff --git a/man/symbols2indices.Rd b/man/symbols2indices.Rd
index 77fcedc..6a2d3b4 100644
--- a/man/symbols2indices.Rd
+++ b/man/symbols2indices.Rd
@@ -22,5 +22,7 @@ Typically, \code{symbols} is the vector of symbols of genes on a microarray, and
\seealso{
\code{\link{romer}}, \code{\link{mroast}}, \code{\link{camera}}
+
+There is a topic page on \link{10.GeneSetTests}.
}
\author{Gordon Smyth and Yifang Hu}
\ No newline at end of file
diff --git a/man/targetsA2C.Rd b/man/targetsA2C.Rd
index 1b26f51..e14a682 100755
--- a/man/targetsA2C.Rd
+++ b/man/targetsA2C.Rd
@@ -6,7 +6,8 @@
Convert a two-color targets dataframe with one row per array to one with one row per channel.
}
\usage{
-targetsA2C(targets, channel.codes=c(1,2), channel.columns=list(Target=c("Cy3","Cy5")), grep=FALSE)
+targetsA2C(targets, channel.codes = c(1,2), channel.columns = list(Target=c("Cy3","Cy5")),
+ grep = FALSE)
}
\arguments{
\item{targets}{data.frame with one row per array giving information about target samples associated covariates.}
@@ -38,4 +39,4 @@ targets <- data.frame(FileName=c("file1.gpr","file2.gpr"),Cy3=c("WT","KO"),Cy5=c
targetsA2C(targets)
}
\author{Gordon Smyth}
-\keyword{htest}
+\concept{separate channel analysis}
diff --git a/man/topRomer.Rd b/man/topRomer.Rd
index dacd47e..62f9f51 100644
--- a/man/topRomer.Rd
+++ b/man/topRomer.Rd
@@ -24,5 +24,10 @@ This function takes the results from romer and returns a number of top gene set
# See romer for examples
}
-\seealso{\code{\link{romer}}}
+\seealso{
+\code{\link{romer}}
+
+There is a topic page on \link{10.GeneSetTests}.
+}
+
\author{Gordon Smyth and Yifang Hu}
\ No newline at end of file
diff --git a/man/topSplice.Rd b/man/topSplice.Rd
new file mode 100644
index 0000000..3362bc8
--- /dev/null
+++ b/man/topSplice.Rd
@@ -0,0 +1,34 @@
+\title{Top table of differentially spliced genes or exons}
+\name{topSplice}
+\alias{topSplice}
+\description{
+Top table ranking the most differentially spliced genes or exons.
+}
+\usage{
+topSplice(fit, coef=ncol(fit), level="exon", number=10, FDR=1)
+}
+\arguments{
+ \item{fit}{\code{MArrayLM} fit object produced by \code{diffSplice}.}
+ \item{coef}{the coefficient (column) of fit for which differentially splicing is assessed.}
+ \item{level}{character string, should the table be by \code{"exon"} or by \code{"gene"}.}
+ \item{number}{integer, maximum number of rows to output.}
+ \item{FDR}{numeric, only show exons or genes with false discovery rate less than this cutoff.}
+}
+
+\details{
+Ranks genes by the Plots interaction log2-fold-change by exon for the specified gene.
+}
+
+\value{A data.frame with any annotation columns found in \code{fit} plus the following columns
+ \item{logFC}{log2-fold change of exon vs other exons for the same gene (if \code{level="exon"})}
+ \item{t}{moderated t-statistic (if \code{level="exon"})}
+ \item{F}{moderated F-statistic (if \code{level="gene"})}
+ \item{P.Value}{p-value}
+ \item{FDR}{false discovery rate}
+}
+
+\author{Gordon Smyth}
+\seealso{
+An overview of diagnostic functions available in LIMMA is given in \link{09.Diagnostics}.
+}
+\examples{# See diffSplice}
diff --git a/man/toptable.Rd b/man/toptable.Rd
index da794cb..5c03114 100755
--- a/man/toptable.Rd
+++ b/man/toptable.Rd
@@ -11,11 +11,10 @@ Extract a table of the top-ranked genes from a linear model fit.
topTable(fit, coef=NULL, number=10, genelist=fit$genes, adjust.method="BH",
sort.by="B", resort.by=NULL, p.value=1, lfc=0, confint=FALSE)
toptable(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, ...)
+ sort.by="B", resort.by=NULL, p.value=1, lfc=0, confint=FALSE, \dots)
topTableF(fit, number=10, genelist=fit$genes, adjust.method="BH",
sort.by="F", p.value=1, lfc=0)
-topTreat(fit, coef=1, number=10, genelist=fit$genes, adjust.method="BH",
- sort.by="p", resort.by=NULL, p.value=1)
+topTreat(fit, coef=1, sort.by="p", resort.by=NULL, \dots)
}
\arguments{
\item{fit}{list containing a linear model fit produced by \code{lmFit}, \code{lm.series}, \code{gls.series} or \code{mrlm}.
@@ -38,15 +37,15 @@ topTreat(fit, coef=1, number=10, genelist=fit$genes, adjust.method="BH",
\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}.}
- \item{confint}{logical, should 95\% confidence intervals be output for \code{logFC}?}
- \item{...}{any other arguments are passed to \code{ebayes} if \code{eb} is \code{NULL}}
+ \item{confint}{logical, should confidence 95\% intervals be output for \code{logFC}? Alternatively, can take a numeric value between zero and one specifying the confidence level required.}
+ \item{\dots}{For \code{toptable}, other arguments are passed to \code{ebayes} (if \code{eb=NULL}). For \code{topTreat}, other arguments are passed to \code{topTable}.}
}
\value{
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})}
- \item{CI.975}{right limit of confidence interval for \code{logFC} (if \code{confint=TRUE})}
+ \item{CI.L}{left limit of confidence interval for \code{logFC} (if \code{confint=TRUE} or \code{confint} is numeric)}
+ \item{CI.R}{right limit of confidence interval for \code{logFC} (if \code{confint=TRUE} or \code{confint} is numeric)}
\item{AveExpr}{average log2-expression for the probe over all arrays and channels, same as \code{Amean} in the \code{MarrayLM} object}
\item{t}{moderated t-statistic (omitted for \code{topTableF})}
\item{F}{moderated F-statistic (omitted for \code{topTable} unless more than one coef is specified)}
diff --git a/man/tricubeMovingAverage.Rd b/man/tricubeMovingAverage.Rd
new file mode 100644
index 0000000..25b8661
--- /dev/null
+++ b/man/tricubeMovingAverage.Rd
@@ -0,0 +1,36 @@
+\name{tricubeMovingAverage}
+\alias{tricubeMovingAverage}
+\alias{tricubeMovingAverage}
+\title{Moving Average Smoother With Tricube Weights}
+\description{
+Apply a moving average smoother to the columns of a matrix.
+}
+\usage{
+tricubeMovingAverage(x,span=0.5,full.length=TRUE)
+}
+\arguments{
+ \item{x}{numeric vector}
+ \item{span}{proportion of points included in the local window}
+ \item{full.length}{logical value, should output have same number of length as input?}
+}
+\details{
+This function smooths a vector (considered as a time series) using moving average with tricube weights.
+It is similar to a loess curve of degree zero, with a couple of differences:
+a continuity correction is applied when computing the neighbouring points and, when \code{full.length=TRUE}, the span halves at the end points.
+
+The \code{filter} function in the stats package is called to do the low-level calculations.
+}
+\value{
+Numeric vector of smoothed values.
+If \code{full.length=TRUE}, of same length as \code{x}.
+If \code{full.length=FALSE}, has \code{width-1} fewer rows than \code{x}.
+}
+\examples{
+x <- rbinom(100,size=1,prob=0.5)
+plot(1:100,tricubeMovingAverage(x))
+}
+\seealso{
+\code{\link{filter}}
+}
+\author{Gordon Smyth}
+\keyword{smooth}
diff --git a/man/volcanoplot.Rd b/man/volcanoplot.Rd
index 953388d..046e765 100755
--- a/man/volcanoplot.Rd
+++ b/man/volcanoplot.Rd
@@ -5,7 +5,8 @@
Creates a volcano plot of log-fold changes versus log-odds of differential expression.
}
\usage{
-volcanoplot(fit, coef=1, highlight=0, names=fit$genes$ID, xlab="Log Fold Change", ylab="Log Odds", pch=16, cex=0.35, ...)
+volcanoplot(fit, coef=1, highlight=0, names=fit$genes$ID,
+ xlab="Log Fold Change", ylab="Log Odds", pch=16, cex=0.35, ...)
}
\arguments{
\item{fit}{an \code{MArrayLM} fitted linear model object}
@@ -31,4 +32,3 @@ An overview of presentation plots following the fitting of a linear model in LIM
\examples{
# See lmFit examples
}
-\keyword{hplot}
diff --git a/man/voom.Rd b/man/voom.Rd
index 17796f8..b8cf5c6 100644
--- a/man/voom.Rd
+++ b/man/voom.Rd
@@ -7,7 +7,8 @@ The data are then ready for linear modelling.
}
\usage{
-voom(counts, design = NULL, lib.size = NULL, normalize.method = "none", plot = FALSE, span=0.5, ...)
+voom(counts, design = NULL, lib.size = NULL, normalize.method = "none",
+ plot = FALSE, span=0.5, ...)
}
\arguments{
\item{counts}{a numeric \code{matrix} containing raw counts, or an \code{ExpressionSet} containing raw counts, or a \code{DGEList} object.}
@@ -55,8 +56,8 @@ PhD Thesis. University of Melbourne, Australia.
Law, CW, Chen, Y, Shi, W, Smyth, GK (2014).
Voom: precision weights unlock linear model analysis tools for RNA-seq read counts.
-\emph{Genome Biology} (Accepted 9 January 2014).
-\url{http://www.statsci.org/smyth/pubs/VoomPreprint.pdf}
+\emph{Genome Biology} 15, R29.
+\url{http://genomebiology.com/2014/15/2/R29}
}
\seealso{
diff --git a/man/vooma.Rd b/man/vooma.Rd
index 1b6157c..f09e3e8 100644
--- a/man/vooma.Rd
+++ b/man/vooma.Rd
@@ -10,7 +10,8 @@ Estimate the mean-variance relationship and use this to compute appropriate obse
\usage{
vooma(y, design=NULL, correlation, block=NULL, plot=FALSE, span=NULL)
-voomaByGroup(y, group, design=NULL, correlation, block=NULL, plot=FALSE, span=NULL, col=NULL)
+voomaByGroup(y, group, design=NULL, correlation, block=NULL,
+ plot=FALSE, span=NULL, col=NULL)
}
\arguments{
diff --git a/man/weightedLowess.Rd b/man/weightedLowess.Rd
new file mode 100644
index 0000000..091697c
--- /dev/null
+++ b/man/weightedLowess.Rd
@@ -0,0 +1,72 @@
+\name{weightedLoess}
+\alias{weightedLowess}
+
+\title{Lowess fit with weighting}
+
+\description{Fit robust lowess curves of degree 1 to weighted covariates and responses.}
+
+\usage{
+weightedLowess(x, y, weights = rep(1, length(y)),
+ delta=NULL, npts = 200, span = 0.3, iterations = 4)
+}
+
+\arguments{
+\item{x}{a numeric vector of covariates}
+\item{y}{a numeric vector of response values}
+\item{weights}{a numeric vector containing frequency weights for each covariate}
+\item{delta}{a numeric scalar specifying the maximum distance between adjacent points}
+\item{npts}{an integer scalar specifying the approximate number of points to use when computing \code{delta}}
+\item{span}{a numeric scalar specifying the width of the smoothing window as a proportion of the total weight}
+\item{iterations}{an integer scalar specifying the number of robustifying iterations}
+}
+
+\value{
+A list of numeric vectors for the fitted responses, the residuals, the robustifying weights and the chosen delta.
+}
+
+\details{
+This function extends the lowess algorithm to handle non-negative prior weights. These weights are
+used during span calculations such that the span distance for each point must include the
+specified proportion of all weights. They are also applied during weighted linear regression to
+compute the fitted value (in addition to the tricube weights determined by \code{span}). For integer
+weights, the prior weights are equivalent to using \code{rep(..., w)} on \code{x} and \code{y} prior to fitting.
+
+For large vectors, running time is reduced by only performing locally weighted regression for
+several points. Fitted values for all points adjacent to the chosen points are computed by
+linear interpolation between the chosen points. For this purpose, the first and last points are always
+chosen. Note that the regression itself uses all (neighbouring) points.
+
+Points are defined as adjacent to a chosen point if the distance to the latter is positive and less
+than \code{delta}. The first chosen point is that corresponding to the smallest covariate; the
+next chosen point is then the next non-adjacent point, and so on. By default, the smallest \code{delta}
+is chosen to obtain a number of chosen points approximately equal to the specified \code{npts}.
+Increasing \code{npts} or supplying a small \code{delta} will improve the accuracy of the fit (i.e.
+closer to the full lowess procedure) at the cost of running time.
+
+Robustification is performed using the magnitude of the residuals. Residuals greater than 6 times the
+median residual are assigned weights of zero. Otherwise, Tukey's biweight function is applied.
+Weights are then used for weighted linear regression. Greater values of \code{iterations} will
+provide greater robustness.
+}
+
+\author{Aaron Lun}
+
+\seealso{
+\code{\link{lowess}}
+}
+
+\examples{
+y <- rt(100,df=4)
+x <- runif(100)
+w <- runif(100)
+out <- weightedLowess(x, y, w, span=0.7)
+plot(x,y,cex=w)
+o <- order(x)
+lines(x[o],out$fitted[o],col="red")
+}
+
+\references{
+Cleveland, W.S. (1979).
+Robust Locally Weighted Regression and Smoothing Scatterplots.
+\emph{Journal of the American Statistical Association} 74, 829-836.
+}
diff --git a/man/writefit.Rd b/man/writefit.Rd
index ae9ddf0..5bf6d16 100755
--- a/man/writefit.Rd
+++ b/man/writefit.Rd
@@ -8,7 +8,8 @@ Write a microarray linear model fit to a file.
}
\usage{
-write.fit(fit, results=NULL, file, digits=3, adjust="none", method="separate", F.adjust="none", sep="\t", ...)
+write.fit(fit, results=NULL, file, digits=3, adjust="none", method="separate",
+ F.adjust="none", sep="\t", ...)
}
\arguments{
diff --git a/man/zscore.Rd b/man/zscore.Rd
index 9862452..edb8368 100755
--- a/man/zscore.Rd
+++ b/man/zscore.Rd
@@ -14,22 +14,23 @@ Compute z-score equivalents of non-normal random deviates.
\usage{
zscore(q, distribution, ...)
zscoreGamma(q, shape, rate = 1, scale = 1/rate)
-zscoreT(x, df)
+zscoreT(x, df, approx=FALSE)
tZscore(x, df)
zscoreHyper(q, m, n, k)
}
\arguments{
-\item{q, x}{numeric vector or matrix giving deviates of a random variable}
-\item{distribution}{character name of probabability distribution for which a cumulative distribution function exists}
-\item{\ldots}{other arguments specify distributional parameters and are passed to the cumulative distribution}
-\item{shape}{gamma shape parameter (>0)}
-\item{rate}{gamma rate parameter (>0)}
-\item{scale}{gamma scale parameter (>0)}
-\item{df}{degrees of freedom (>0 for \code{zscoreT} or >=1 for \code{tZscore})}
-\item{m}{as for \code{\link{qhyper}}}
-\item{n}{as for \code{\link{qhyper}}}
-\item{k}{as for \code{\link{qhyper}}}
+ \item{q, x}{numeric vector or matrix giving deviates of a random variable}
+ \item{distribution}{character name of probabability distribution for which a cumulative distribution function exists}
+ \item{\ldots}{other arguments specify distributional parameters and are passed to the cumulative distribution function}
+ \item{shape}{gamma shape parameter (>0)}
+ \item{rate}{gamma rate parameter (>0)}
+ \item{scale}{gamma scale parameter (>0)}
+ \item{df}{degrees of freedom (>0 for \code{zscoreT} or >=1 for \code{tZscore})}
+ \item{approx}{logical, if \code{TRUE} then a fast approximation is used to convert t-statistics into z-scores. If \code{FALSE}, z-scores will be exact.}
+ \item{m}{as for \code{\link{qhyper}}}
+ \item{n}{as for \code{\link{qhyper}}}
+ \item{k}{as for \code{\link{qhyper}}}
}
\value{
@@ -48,11 +49,24 @@ The argument \code{distribution} is the name of the cumulative distribution func
\code{tZscore} is the inverse of \code{zscoreT}, and computes t-distribution equivalents for standard normal deviates.
-Care is taken to do the computations accurately in both tails of the distributions.
+The transformation to z-scores is done by converting to log tail probabilities, and then using \code{qnorm}.
+For numerical accuracy, the left or right tail is used, depending on which is likely to be smaller.
+
+If \code{approx=TRUE}, then the approximation from Hill (1970) is used to convert t-statistics to z-scores directly without computing tail probabilities.
+Brophy (1987) showed this to be most accurate of a variety of possible closed-form transformations.
}
\author{Gordon Smyth}
+\references{
+Hill, GW (1970). Algorithm 395: Student's t-distribution.
+\emph{Communications of the ACM} 13, 617-620.
+
+Brophy, AL (1987).
+Efficient estimation of probabilities in the t distribution.
+\emph{Behavior Research Methods} 19, 462--466.
+}
+
\seealso{
\code{\link{qnorm}}, \code{\link{pgamma}}, \code{\link{pt}} in the stats package.
}
@@ -66,5 +80,3 @@ zscoreGamma(c(1,2.5), shape=0.5, scale=2)
zscoreT(2, df=3)
tZscore(2, df=3)
}
-
-\keyword{distribution}
diff --git a/src/weighted_lowess.c b/src/weighted_lowess.c
new file mode 100644
index 0000000..da02c0c
--- /dev/null
+++ b/src/weighted_lowess.c
@@ -0,0 +1,271 @@
+#include "R.h"
+#include "Rdefines.h"
+#include <math.h>
+
+#define THRESHOLD 0.0000001
+
+/* This function determines the number of points to be used in the analysis,
+ * based on the chosen delta. It returns that number as well as the index
+ * values of those points in a pointer reference. Note that the first
+ * and last point (i.e. smallest and largest) are always considered.
+ */
+
+void find_seeds (int ** indices, int * number, const double* xptr, const int npts, const double delta) {
+ int pt, last_pt=0;
+ int total=2;
+ for (pt=1; pt<npts-1; ++pt) {
+ if (xptr[pt] - xptr[last_pt] > delta) {
+ ++total;
+ last_pt=pt;
+ }
+ }
+ (*number)=total;
+
+ /* Second pass to actually record the indices at which these events happen. */
+ int* idptr=(int*)R_alloc(*number, sizeof(int));
+ idptr[0]=0;
+ total=1;
+ last_pt=0;
+ for (pt=1; pt<npts-1; ++pt) {
+ if (xptr[pt] - xptr[last_pt] > delta) {
+ idptr[total]=pt;
+ last_pt=pt;
+ ++total;
+ }
+ }
+ idptr[total]=npts-1;
+ (*indices)=idptr;
+ return;
+}
+
+/* This function identifies the start and end index in the span for each chosen sampling
+ * point. It returns two arrays via reference containing said indices. It also returns
+ * an array containing the maximum distance between points at each span.
+ *
+ * We don't use the update-based algorithm in Cleveland's paper, as it ceases to be
+ * numerically stable once you throw in double-precision weights. It's not particularly
+ * amenable to updating through cycles of addition and subtraction. At any rate, the
+ * algorithm as a whole remains quadratic (as weights must be recomputed) so there's no
+ * damage to scalability.
+ */
+
+void find_limits (const int* indices, const int num, const double* xptr, const double* wptr,
+ const int npts, const double spanweight, int** start, int** end, double** dist) {
+ int* spbegin=(int*)R_alloc(num, sizeof(int));
+ int* spend=(int*)R_alloc(num, sizeof(int));
+ double* spdist=(double*)R_alloc(num, sizeof(double));
+
+ int curx;
+ for (curx=0; curx<num; ++curx) {
+ const int curpt=indices[curx];
+ int left=curpt, right=curpt;
+ double curw=wptr[curpt];
+ int ende=(curpt==npts-1), ends=(curpt==0);
+ double mdist=0, ldist, rdist;
+
+ while (curw < spanweight && (!ende || !ends)) {
+ if (ende) {
+ /* Can only extend backwards. */
+ --left;
+ curw+=wptr[left];
+ if (left==0) { ends=1; }
+ ldist=xptr[curpt]-xptr[left];
+ if (mdist < ldist) { mdist=ldist; }
+ } else if (ends) {
+ /* Can only extend forwards. */
+ ++right;
+ curw+=wptr[right];
+ if (right==npts-1) { ende=1; }
+ rdist=xptr[right]-xptr[curpt];
+ if (mdist < rdist) { mdist=rdist; }
+ } else {
+ /* Can do either; extending by the one that minimizes the curpt mdist. */
+ ldist=xptr[curpt]-xptr[left-1];
+ rdist=xptr[right+1]-xptr[curpt];
+ if (ldist < rdist) {
+ --left;
+ curw+=wptr[left];
+ if (left==0) { ends=1; }
+ if (mdist < ldist) { mdist=ldist; }
+ } else {
+ ++right;
+ curw+=wptr[right];
+ if (right==npts-1) { ende=1; }
+ if (mdist < rdist) { mdist=rdist; }
+ }
+ }
+ }
+
+ /* Extending to ties. */
+ while (left>0 && xptr[left]==xptr[left-1]) { --left; }
+ while (right<npts-1 && xptr[right]==xptr[right+1]) { ++right; }
+
+ /* Recording */
+ spbegin[curx]=left;
+ spend[curx]=right;
+ spdist[curx]=mdist;
+ }
+
+ (*start)=spbegin;
+ (*end)=spend;
+ (*dist)=spdist;
+ return;
+}
+
+/* Computes the lowess fit at a given point using linear regression with a combination of tricube,
+ * prior and robustness weighting. Some additional effort is put in to avoid numerical instability
+ * and undefined values when divisors are near zero.
+ */
+
+double lowess_fit (const double* xptr, const double* yptr, const double* wptr, const double* rwptr,
+ const int npts, const int curpt, const int left, const int right, const double dist, double* work) {
+ double ymean=0, allweight=0;
+ int pt;
+ if (dist < THRESHOLD) {
+ for (pt=left; pt<=right; ++pt) {
+ work[pt]=wptr[pt]*rwptr[pt];
+ ymean+=yptr[pt]*work[pt];
+ allweight+=work[pt];
+ }
+ ymean/=allweight;
+ return ymean;
+ }
+ double xmean=0;
+ for (pt=left; pt<=right; ++pt) {
+ work[pt]=pow(1-pow(fabs(xptr[curpt]-xptr[pt])/dist, 3.0), 3.0)*wptr[pt]*rwptr[pt];
+ xmean+=work[pt]*xptr[pt];
+ ymean+=work[pt]*yptr[pt];
+ allweight+=work[pt];
+ }
+ xmean/=allweight;
+ ymean/=allweight;
+
+ double var=0, covar=0, temp;
+ for (pt=left; pt<=right; ++pt) {
+ temp=xptr[pt]-xmean;
+ var+=temp*temp*work[pt];
+ covar+=temp*(yptr[pt]-ymean)*work[pt];
+ }
+ if (var < THRESHOLD) { return ymean; }
+
+ const double slope=covar/var;
+ const double intercept=ymean-slope*xmean;
+ return slope*xptr[curpt]+intercept;
+}
+
+/* This is a C version of the local weighted regression (lowess) trend fitting algorithm,
+ * based on the Fortran code in lowess.f from http://www.netlib.org/go written by Cleveland.
+ * Consideration of non-equal prior weights is added to the span calculations and linear
+ * regression. These weights are intended to have the equivalent effect of frequency weights
+ * (at least, in the integer case; extended by analogy to all non-negative values).
+ */
+
+SEXP weighted_lowess(SEXP covariate, SEXP response, SEXP weight, SEXP span, SEXP iter, SEXP delta) {
+ if (!IS_NUMERIC(covariate)) { error("covariates must be double precision"); }
+ if (!IS_NUMERIC(response)) { error("responses must be double precision"); }
+ if (!IS_NUMERIC(weight)) { error("weights must be double precision"); }
+
+ const int npts=LENGTH(covariate);
+ if (npts!=LENGTH(response) || npts!=LENGTH(weight)) { error("weight, covariate and response vectors have unequal lengths"); }
+ if (npts<2) { error("need at least two points"); }
+ const double* covptr=NUMERIC_POINTER(covariate);
+ const double* resptr=NUMERIC_POINTER(response);
+ const double* weiptr=NUMERIC_POINTER(weight);
+
+ if (!IS_NUMERIC(span) || LENGTH(span)!=1) { error("span should be a double-precision scalar"); }
+ const double spv=NUMERIC_VALUE(span);
+ if (!IS_INTEGER(iter) || LENGTH(iter)!=1) { error("number of robustness iterations should be an integer scalar"); }
+ const int niter=INTEGER_VALUE(iter);
+ if (niter<=0) { error("number of robustness iterations should be positive"); }
+ if (!IS_NUMERIC(delta) || LENGTH(delta)!=1) { error("delta should be a double-precision scalar"); }
+ const double dv=NUMERIC_VALUE(delta);
+
+ /*** NO MORE ERRORS AT THIS POINT, MEMORY ASSIGNMENTS ARE ACTIVE. ***/
+
+ /* Computing the span weight that each span must achieve. */
+ double totalweight=0;
+ int pt;
+ for (pt=0; pt<npts; ++pt) { totalweight+=weiptr[pt]; }
+ double spanweight=totalweight*spv;
+
+ /* Setting up the indices of points for sampling; the frame start and end for those indices, and the max dist. */
+ int *seed_index;
+ int nseeds;
+ find_seeds(&seed_index, &nseeds, covptr, npts, dv);
+ int *frame_start, *frame_end;
+ double* max_dist;
+ find_limits (seed_index, nseeds, covptr, weiptr, npts, spanweight, &frame_start, &frame_end, &max_dist);
+
+ /* Setting up arrays to hold the fitted values, residuals and robustness weights. */
+ SEXP output=PROTECT(NEW_LIST(2));
+ SET_VECTOR_ELT(output, 0, NEW_NUMERIC(npts));
+ double* fitptr=NUMERIC_POINTER(VECTOR_ELT(output, 0));
+ double* rsdptr=(double*)R_alloc(npts, sizeof(double));
+ SET_VECTOR_ELT(output, 1, NEW_NUMERIC(npts));
+ double* robptr=NUMERIC_POINTER(VECTOR_ELT(output, 1));
+ int* rorptr=(int*)R_alloc(npts, sizeof(int));
+ for (pt=0; pt<npts; ++pt) { robptr[pt]=1; }
+
+ /* Robustness iterations. */
+ int it=0;
+ for (it=0; it<niter; ++it) {
+ int cur_seed, last_pt=0, subpt;
+ double current;
+
+ /* Computing fitted values for seed points, and interpolating to the intervening points. */
+ fitptr[0]=lowess_fit(covptr, resptr, weiptr, robptr, npts, 0, frame_start[0], frame_end[0], max_dist[0], rsdptr);
+ for (cur_seed=1; cur_seed<nseeds; ++cur_seed) {
+ pt=seed_index[cur_seed];
+ fitptr[pt]=lowess_fit(covptr, resptr, weiptr, robptr, npts, pt, frame_start[cur_seed],
+ frame_end[cur_seed], max_dist[cur_seed], rsdptr); /* using rsdptr as a holding cell. */
+
+ if (pt-last_pt > 1) {
+ /* Some protection is provided against infinite slopes. This shouldn't be
+ * a problem for non-zero delta; the only concern is at the final point
+ * where the covariate distance may be zero. Besides, if delta is not
+ * positive, pt-last_pt could never be 1 so we'd never reach this point.
+ */
+ current = covptr[pt]-covptr[last_pt];
+ if (current > THRESHOLD) {
+ const double slope=(fitptr[pt]-fitptr[last_pt])/current;
+ const double intercept=fitptr[pt] - slope*covptr[pt];
+ for (subpt=last_pt+1; subpt<pt; ++subpt) { fitptr[subpt]=slope*covptr[subpt]+intercept; }
+ } else {
+ const double endave=0.5*(fitptr[pt]+fitptr[last_pt]);
+ for (subpt=last_pt+1; subpt<pt; ++subpt) { fitptr[subpt]=endave; }
+ }
+ }
+ last_pt=pt;
+ }
+
+ /* Computing the weighted MAD of the absolute values of the residuals. */
+ for (pt=0; pt<npts; ++pt) {
+ rsdptr[pt]=fabs(resptr[pt]-fitptr[pt]);
+ rorptr[pt]=pt;
+ }
+ rsort_with_index(rsdptr, rorptr, npts);
+ current=0;
+ double cmad=THRESHOLD;
+ const double halfweight=totalweight/2;
+ for (pt=0; pt<npts; ++pt) {
+ current+=weiptr[rorptr[pt]];
+ if (current==halfweight) { /* In the unlikely event of an exact match. */
+ cmad=3*(rsdptr[pt]+rsdptr[pt+1]);
+ break;
+ } else if (current>halfweight) {
+ cmad=6*rsdptr[pt];
+ break;
+ }
+ }
+
+ /* Computing the robustness weights. */
+ for (pt=0; pt<npts; ++pt) {
+ if (rsdptr[pt]<cmad) {
+ robptr[rorptr[pt]]=pow(1-pow(rsdptr[pt]/cmad, 2.0), 2.0);
+ } else { robptr[rorptr[pt]]=0; }
+ }
+ }
+
+ UNPROTECT(1);
+ return output;
+}
diff --git a/tests/limma-Tests.R b/tests/limma-Tests.R
index 46f8b2d..7dc83d3 100755
--- a/tests/limma-Tests.R
+++ b/tests/limma-Tests.R
@@ -45,15 +45,22 @@ f1 <- quantile(out$fitted)
r1 <- quantile(out$residual)
w <- rep(1,100)
w[1:50] <- 0.5
-out <- loessFit(y,x,weights=w)
+out <- loessFit(y,x,weights=w,method="weightedLowess")
f2 <- quantile(out$fitted)
r2 <- quantile(out$residual)
-w <- rep(1,100)
-w[2*(1:50)] <- 0
-out <- loessFit(y,x,weights=w)
+out <- loessFit(y,x,weights=w,method="locfit")
f3 <- quantile(out$fitted)
r3 <- quantile(out$residual)
-data.frame(f1,f2,f3,r1,r2,r3)
+out <- loessFit(y,x,weights=w,method="loess")
+f4 <- quantile(out$fitted)
+r4 <- quantile(out$residual)
+w <- rep(1,100)
+w[2*(1:50)] <- 0
+out <- loessFit(y,x,weights=w,method="weightedLowess")
+f5 <- quantile(out$fitted)
+r5 <- quantile(out$residual)
+data.frame(f1,f2,f3,f4,f5)
+data.frame(r1,r2,r3,r4,r5)
### normalizeWithinArrays
diff --git a/tests/limma-Tests.Rout.save b/tests/limma-Tests.Rout.save
index ef6c8f3..06f2f23 100755
--- a/tests/limma-Tests.Rout.save
+++ b/tests/limma-Tests.Rout.save
@@ -1,7 +1,7 @@
-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 version 3.1.0 (2014-04-10) -- "Spring Dance"
+Copyright (C) 2014 The R Foundation for Statistical Computing
+Platform: x86_64-w64-mingw32/x64 (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
@@ -591,21 +591,34 @@ $G
> r1 <- quantile(out$residual)
> w <- rep(1,100)
> w[1:50] <- 0.5
-> out <- loessFit(y,x,weights=w)
+> out <- loessFit(y,x,weights=w,method="weightedLowess")
> f2 <- quantile(out$fitted)
> r2 <- quantile(out$residual)
-> w <- rep(1,100)
-> w[2*(1:50)] <- 0
-> out <- loessFit(y,x,weights=w)
+> out <- loessFit(y,x,weights=w,method="locfit")
> 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.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
+> out <- loessFit(y,x,weights=w,method="loess")
+> f4 <- quantile(out$fitted)
+> r4 <- quantile(out$residual)
+> w <- rep(1,100)
+> w[2*(1:50)] <- 0
+> out <- loessFit(y,x,weights=w,method="weightedLowess")
+> f5 <- quantile(out$fitted)
+> r5 <- quantile(out$residual)
+> data.frame(f1,f2,f3,f4,f5)
+ f1 f2 f3 f4 f5
+0% -0.78835384 -0.687432210 -0.78957137 -0.76756060 -0.63778292
+25% -0.18340154 -0.179683572 -0.18979269 -0.16773223 -0.38064318
+50% -0.11492924 -0.114796040 -0.12087983 -0.07185314 -0.15971879
+75% 0.01507921 -0.008145125 -0.01857508 0.04030634 0.07839396
+100% 0.21653837 0.145106033 0.19214597 0.21417361 0.51836274
+> data.frame(r1,r2,r3,r4,r5)
+ r1 r2 r3 r4 r5
+0% -2.04434053 -2.05132680 -2.02404318 -2.101242874 -2.22280633
+25% -0.59321065 -0.57200209 -0.58975649 -0.577887481 -0.71037756
+50% 0.05874864 0.04514326 0.08335198 -0.001769806 0.06785517
+75% 0.56010750 0.55124530 0.57618740 0.561454370 0.65383830
+100% 2.57936026 2.64549799 2.57549257 2.402324533 2.28648835
>
> ### normalizeWithinArrays
>
@@ -643,12 +656,12 @@ Array 2 corrected
> MA <- normalizeWithinArrays(RGb,method="loess")
> summary(MA$M)
V1 V2
- Min. :-5.85163 Min. :-5.69877
- 1st Qu.:-1.18482 1st Qu.:-1.55421
- Median :-0.21631 Median : 0.06267
- Mean : 0.03613 Mean :-0.05369
- 3rd Qu.: 1.49673 3rd Qu.: 1.41900
- Max. : 7.07528 Max. : 6.28902
+ Min. :-5.88044 Min. :-5.66985
+ 1st Qu.:-1.18483 1st Qu.:-1.57014
+ Median :-0.21632 Median : 0.04823
+ Mean : 0.03487 Mean :-0.05481
+ 3rd Qu.: 1.49669 3rd Qu.: 1.45113
+ Max. : 7.07324 Max. : 6.19744
> #MA <- normalizeWithinArrays(RG[,1:2], mouse.setup, method="robustspline")
> #MA$M[1:5,]
> #MA <- normalizeWithinArrays(mouse.data, mouse.setup)
@@ -659,11 +672,11 @@ Array 2 corrected
> MA2 <- normalizeBetweenArrays(MA,method="scale")
> MA$M[1:5,]
[,1] [,2]
-[1,] -1.1623390 4.5343276
-[2,] 0.8971391 0.3495635
-[3,] 2.8247455 1.4459533
-[4,] -1.8533288 0.4894799
-[5,] 1.9158416 -5.5363732
+[1,] -1.1689588 4.5558123
+[2,] 0.8971363 0.3296544
+[3,] 2.8247439 1.4249960
+[4,] -1.8533240 0.4804851
+[5,] 1.9158459 -5.5087631
> MA$A[1:5,]
[,1] [,2]
[1,] -2.48465011 -2.4041550
@@ -674,11 +687,11 @@ Array 2 corrected
> MA2 <- normalizeBetweenArrays(MA,method="quantile")
> MA$M[1:5,]
[,1] [,2]
-[1,] -1.1623390 4.5343276
-[2,] 0.8971391 0.3495635
-[3,] 2.8247455 1.4459533
-[4,] -1.8533288 0.4894799
-[5,] 1.9158416 -5.5363732
+[1,] -1.1689588 4.5558123
+[2,] 0.8971363 0.3296544
+[3,] 2.8247439 1.4249960
+[4,] -1.8533240 0.4804851
+[5,] 1.9158459 -5.5087631
> MA$A[1:5,]
[,1] [,2]
[1,] -2.48465011 -2.4041550
@@ -1008,8 +1021,8 @@ $df.residual
$cov.coefficients
alpha-beta mu+alpha mu+beta
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
+mu+alpha 0.3333333 3.333333e-01 5.551115e-17
+mu+beta -0.3333333 5.551115e-17 3.333333e-01
$Amean
[1] 0.2034961 0.1954604 -0.2863347 0.1188659 0.1784593
@@ -1143,37 +1156,40 @@ c 0.3482980 -0.4820695 -0.3841313
> y[iset1,3:4] <- y[iset1,3:4]+3
> iset2 <- 6:10
> roast(y=y,iset1,design,contrast=2)
- Active.Prop P.Value
-Down 0 0.996
-Up 1 0.005
-Mixed 1 0.008
+ Active.Prop P.Value
+Down 0 0.996498249
+Up 1 0.004002001
+UpOrDown 1 0.008000000
+Mixed 1 0.008000000
> roast(y=y,iset1,design,contrast=2,array.weights=c(0.5,1,0.5,1))
- Active.Prop P.Value
-Down 0 0.999
-Up 1 0.002
-Mixed 1 0.003
+ Active.Prop P.Value
+Down 0 0.99899950
+Up 1 0.00150075
+UpOrDown 1 0.00300000
+Mixed 1 0.00300000
> w <- matrix(runif(100*4),100,4)
> roast(y=y,iset1,design,contrast=2,weights=w)
- Active.Prop P.Value
-Down 0 0.999
-Up 1 0.002
-Mixed 1 0.002
+ Active.Prop P.Value
+Down 0 0.9994997
+Up 1 0.0010005
+UpOrDown 1 0.0020000
+Mixed 1 0.0020000
> mroast(y=y,list(set1=iset1,set2=iset2),design,contrast=2,gene.weights=runif(100))
- NGenes PropDown PropUp Direction PValue FDR PValue.Mixed FDR.Mixed
-set1 5 0 1 Up 0.006 0.010 0.008 0.0150
-set2 5 0 0 Up 0.948 0.947 0.687 0.6865
+ NGenes PropDown PropUp Direction PValue FDR PValue.Mixed FDR.Mixed
+set1 5 0 1 Up 0.008 0.0150 0.008 0.0150
+set2 5 0 0 Up 0.959 0.9585 0.687 0.6865
> mroast(y=y,list(set1=iset1,set2=iset2),design,contrast=2,array.weights=c(0.5,1,0.5,1))
- NGenes PropDown PropUp Direction PValue FDR PValue.Mixed FDR.Mixed
-set1 5 0 1 Up 0.006 0.010 0.004 0.0070
-set2 5 0 0 Up 0.716 0.715 0.658 0.6575
+ NGenes PropDown PropUp Direction PValue FDR PValue.Mixed FDR.Mixed
+set1 5 0 1 Up 0.004 0.0070 0.004 0.0070
+set2 5 0 0 Up 0.679 0.6785 0.658 0.6575
> mroast(y=y,list(set1=iset1,set2=iset2),design,contrast=2,weights=w)
- NGenes PropDown PropUp Direction PValue FDR PValue.Mixed FDR.Mixed
-set1 5 0.0 1 Up 0.004 0.006 0.003 0.0050
-set2 5 0.2 0 Down 0.984 0.983 0.250 0.2495
+ NGenes PropDown PropUp Direction PValue FDR PValue.Mixed FDR.Mixed
+set1 5 0.0 1 Up 0.003 0.0050 0.003 0.0050
+set2 5 0.2 0 Down 0.950 0.9495 0.250 0.2495
> mroast(y=y,list(set1=iset1,set2=iset2),design,contrast=2,weights=w,array.weights=c(0.5,1,0.5,1))
- NGenes PropDown PropUp Direction PValue FDR PValue.Mixed FDR.Mixed
-set1 5 0 1 Up 0.002 0.002 0.001 0.0010
-set2 5 0 0 Down 0.772 0.771 0.146 0.1455
+ NGenes PropDown PropUp Direction PValue FDR PValue.Mixed FDR.Mixed
+set1 5 0 1 Up 0.001 0.0010 0.001 0.0010
+set2 5 0 0 Down 0.791 0.7905 0.146 0.1455
>
> ### camera
>
@@ -1189,10 +1205,11 @@ set2 5 0.1719094 Down 0.9068364381 0.90683644
>
> y <- new("EList",list(E=y))
> roast(y=y,iset1,design,contrast=2)
- Active.Prop P.Value
-Down 0 0.998
-Up 1 0.003
-Mixed 1 0.006
+ Active.Prop P.Value
+Down 0 0.997498749
+Up 1 0.003001501
+UpOrDown 1 0.006000000
+Mixed 1 0.006000000
> camera(y=y,iset1,design,contrast=2)
NGenes Correlation Direction PValue
set1 5 -0.2481655 Up 0.0009047749
@@ -1302,4 +1319,4 @@ Loading required package: splines
>
> proc.time()
user system elapsed
- 1.45 0.09 1.52
+ 1.23 0.12 1.34
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/r-bioc-limma.git
More information about the debian-med-commit
mailing list