[med-svn] [r-bioc-limma] 02/02: Imported Upstream version 3.28.10+dfsg
Andreas Tille
tille at debian.org
Wed Jun 22 15:36:27 UTC 2016
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to annotated tag upstream/3.28.10+dfsg
in repository r-bioc-limma.
commit f1fbb8beedeea8e3d28e24f1da52fb514ec2cd04
Author: Andreas Tille <tille at debian.org>
Date: Wed Jun 22 17:29:35 2016 +0200
Imported Upstream version 3.28.10+dfsg
---
DESCRIPTION | 6 +--
NAMESPACE | 2 +
R/alias2Symbol.R | 20 +++----
R/detectionPValues.R | 41 ++++++++++++++
R/dups.R | 2 +-
R/fitFDistRobustly.R | 2 +-
R/fitmixture.R | 48 +++++++++++++++++
R/geneset-camera.R | 9 ++--
R/geneset-fry.R | 10 ++--
R/goana.R | 16 +++---
R/loessFit.R | 2 +-
R/norm.R | 7 ++-
R/read.idat.R | 127 ++++++++++++++++++++++++++++++--------------
R/sepchannel.R | 2 +-
inst/doc/changelog.txt | 41 +++++++++++++-
inst/doc/intro.pdf | Bin 46191 -> 46191 bytes
man/03reading.Rd | 1 +
man/10GeneSetTests.Rd | 2 +-
man/camera.Rd | 37 ++++++++-----
man/detectionPValue.Rd | 58 ++++++++++++++++++++
man/dupcor.Rd | 40 +++++++++++---
man/fitmixture.Rd | 66 +++++++++++++++++++++++
man/logcosh.Rd | 26 +++++++++
man/propexpr.Rd | 37 ++++++++-----
man/read.idat.Rd | 77 ++++++++++++++++++---------
man/roast.Rd | 6 +--
man/romer.Rd | 2 +-
tests/limma-Tests.R | 4 ++
tests/limma-Tests.Rout.save | 24 +++++++--
29 files changed, 572 insertions(+), 143 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 8b4864d..cb99d7f 100755
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
Package: limma
-Version: 3.28.4
-Date: 2016-05-11
+Version: 3.28.10
+Date: 2016-06-15
Title: Linear Models for Microarray Data
Description: Data analysis, linear models and differential expression for microarray data.
Author: Gordon Smyth [cre,aut], Yifang Hu [ctb], Matthew Ritchie [ctb], Jeremy Silver [ctb], James Wettenhall [ctb], Davis McCarthy [ctb], Di Wu [ctb], Wei Shi [ctb], Belinda Phipson [ctb], Aaron Lun [ctb], Natalie Thorne [ctb], Alicia Oshlack [ctb], Carolyn de Graaf [ctb], Yunshun Chen [ctb], Mette Langaas [ctb], Egil Ferkingstad [ctb], Marcus Davy [ctb], Francois Pepin [ctb], Dongseok Choi [ctb]
@@ -22,4 +22,4 @@ biocViews: ExonArray, GeneExpression, Transcription,
MultipleComparison, Normalization, Preprocessing,
QualityControl
NeedsCompilation: yes
-Packaged: 2016-05-11 01:26:45 UTC; biocbuild
+Packaged: 2016-06-16 01:27:03 UTC; biocbuild
diff --git a/NAMESPACE b/NAMESPACE
index b63b6c6..ea342dc 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -59,6 +59,8 @@ S3method(cbind,MAList)
S3method(cbind,RGList)
S3method(cbind,EList)
S3method(cbind,EListRaw)
+S3method(detectionPValues,default)
+S3method(detectionPValues,EListRaw)
S3method(dim,MAList)
S3method(dim,MArrayLM)
S3method(dim,RGList)
diff --git a/R/alias2Symbol.R b/R/alias2Symbol.R
index 09d719e..46ef938 100644
--- a/R/alias2Symbol.R
+++ b/R/alias2Symbol.R
@@ -11,28 +11,28 @@ alias2Symbol <- function(alias,species="Hs",expand.symbols=FALSE)
# Get access to required annotation functions
suppressPackageStartupMessages(OK <- requireNamespace("AnnotationDbi",quietly=TRUE))
- if(!OK) stop("AnnotationDbi package required but not available")
+ if(!OK) stop("AnnotationDbi package required but not installed (or can't be loaded)")
# Get alias to symbol mappings
switch(species,
Hs = {
suppressPackageStartupMessages(OK <- requireNamespace("org.Hs.eg.db",quietly=TRUE))
- if(!OK) stop("org.Hs.eg.db package required but not available")
+ if(!OK) stop("org.Hs.eg.db package required but not installed (or can't be loaded)")
ALIAS2EG <- org.Hs.eg.db::org.Hs.egALIAS2EG
SYMBOL <- org.Hs.eg.db::org.Hs.egSYMBOL
}, Mm = {
suppressPackageStartupMessages(OK <- requireNamespace("org.Mm.eg.db",quietly=TRUE))
- if(!OK) stop("org.Mm.eg.db package required but not available")
+ if(!OK) stop("org.Mm.eg.db package required but not installed (or can't be loaded)")
ALIAS2EG <- org.Mm.eg.db::org.Mm.egALIAS2EG
SYMBOL <- org.Mm.eg.db::org.Mm.egSYMBOL
}, Rn = {
suppressPackageStartupMessages(OK <- requireNamespace("org.Rn.eg.db",quietly=TRUE))
- if(!OK) stop("org.Rn.eg.db package required but not available")
+ if(!OK) stop("org.Rn.eg.db package required but not installed (or can't be loaded)")
ALIAS2EG <- org.Rn.eg.db::org.Rn.egALIAS2EG
SYMBOL <- org.Rn.eg.db::org.Rn.egSYMBOL
}, Dm = {
suppressPackageStartupMessages(OK <- requireNamespace("org.Dm.eg.db",quietly=TRUE))
- if(!OK) stop("org.Dm.eg.db package required but not available")
+ if(!OK) stop("org.Dm.eg.db package required but not installed (or can't be loaded)")
ALIAS2EG <- org.Dm.eg.db::org.Dm.egALIAS2EG
SYMBOL <- org.Dm.eg.db::org.Dm.egSYMBOL
}
@@ -61,28 +61,28 @@ alias2SymbolTable <- function(alias,species="Hs")
# Get access to required annotation functions
suppressPackageStartupMessages(OK <- requireNamespace("AnnotationDbi",quietly=TRUE))
- if(!OK) stop("AnnotationDbi package required but not available")
+ if(!OK) stop("AnnotationDbi package required but not installed (or can't be loaded)")
# Get alias to symbol mappings
switch(species,
Hs = {
suppressPackageStartupMessages(OK <- requireNamespace("org.Hs.eg.db",quietly=TRUE))
- if(!OK) stop("org.Hs.eg.db package required but not available")
+ if(!OK) stop("org.Hs.eg.db package required but not installed (or can't be loaded)")
ALIAS2EG <- org.Hs.eg.db::org.Hs.egALIAS2EG
SYMBOL <- org.Hs.eg.db::org.Hs.egSYMBOL
}, Mm = {
suppressPackageStartupMessages(OK <- requireNamespace("org.Mm.eg.db",quietly=TRUE))
- if(!OK) stop("org.Mm.eg.db package required but not available")
+ if(!OK) stop("org.Mm.eg.db package required but not installed (or can't be loaded)")
ALIAS2EG <- org.Mm.eg.db::org.Mm.egALIAS2EG
SYMBOL <- org.Mm.eg.db::org.Mm.egSYMBOL
}, Rn = {
suppressPackageStartupMessages(OK <- requireNamespace("org.Rn.eg.db",quietly=TRUE))
- if(!OK) stop("org.Rn.eg.db package required but not available")
+ if(!OK) stop("org.Rn.eg.db package required but not installed (or can't be loaded)")
ALIAS2EG <- org.Rn.eg.db::org.Rn.egALIAS2EG
SYMBOL <- org.Rn.eg.db::org.Rn.egSYMBOL
}, Dm = {
suppressPackageStartupMessages(OK <- requireNamespace("org.Dm.eg.db",quietly=TRUE))
- if(!OK) stop("org.Dm.eg.db package required but not available")
+ if(!OK) stop("org.Dm.eg.db package required but not installed (or can't be loaded)")
ALIAS2EG <- org.Dm.eg.db::org.Dm.egALIAS2EG
SYMBOL <- org.Dm.eg.db::org.Dm.egSYMBOL
}
diff --git a/R/detectionPValues.R b/R/detectionPValues.R
new file mode 100644
index 0000000..a0f8a8f
--- /dev/null
+++ b/R/detectionPValues.R
@@ -0,0 +1,41 @@
+detectionPValues <- function(x,...) UseMethod("detectionPValues")
+
+detectionPValues.EListRaw <- function(x,status=NULL,...)
+# Detection p-values from negative controls for EListRaw
+# Gordon Smyth
+# Created 13 June 2016
+{
+ if(is.null(status)) status <- x$genes$Status
+ detectionPValues(x=x$E, status=status, ...)
+}
+
+detectionPValues.default <- function(x, status, negctrl="negative",...)
+# Detection p-values from negative controls
+# Gordon Smyth
+# Created 13 June 2016. Last modified 14 June 2016.
+{
+ x <- as.matrix(x)
+
+# Identify negative control rows
+ if(missing(status) || is.null(status)) stop("need to specify status vector")
+ isneg <- as.integer(status==negctrl)
+ notneg <- 1L-isneg
+ nneg <- sum(isneg)
+
+# Count the number of negative controls greater than each value
+# cs1 counts number of negative controls >= the observed value
+# cs2 counts number of negative controls > the observed value
+# By taking the average of cs1 and cs2, we count ties as 1/2
+ DetectionPValue <- x
+ cs1 <- cs2 <- rep_len(0,length.out=nrow(x))
+ for (j in 1:ncol(x)) {
+ o1 <- order(x[,j],isneg,decreasing=TRUE)
+ o2 <- order(x[,j],notneg,decreasing=TRUE)
+ cs1[o1] <- cumsum(isneg[o1])
+ cs2[o2] <- cumsum(isneg[o2])
+ DetectionPValue[,j] <- cs1+cs2
+ }
+
+# Convert to p-values
+ DetectionPValue/(2*nneg)
+}
diff --git a/R/dups.R b/R/dups.R
index 90430bf..d6e8802 100755
--- a/R/dups.R
+++ b/R/dups.R
@@ -97,7 +97,7 @@ duplicateCorrelation <- function(object,design=NULL,ndups=2,spacing=1,block=NULL
design <- design %x% rep(1,ndups)
}
- if(!requireNamespace("statmod",quietly=TRUE)) stop("statmod package required but is not available")
+ if(!requireNamespace("statmod",quietly=TRUE)) stop("statmod package required but is not installed")
rho <- rep(NA,ngenes)
nafun <- function(e) NA
for (i in 1:ngenes) {
diff --git a/R/fitFDistRobustly.R b/R/fitFDistRobustly.R
index 9c2b0c6..0d82ad0 100644
--- a/R/fitFDistRobustly.R
+++ b/R/fitFDistRobustly.R
@@ -107,7 +107,7 @@ fitFDistRobustly <- function(x,df1,covariate=NULL,winsor.tail.p=c(0.05,0.1),trac
if(trace) cat("Variance of Winsorized Fisher-z",zwvar,"\n")
# Theoretical Winsorized moments
- if(!requireNamespace("statmod",quietly=TRUE)) stop("statmod package required but is not available")
+ if(!requireNamespace("statmod",quietly=TRUE)) stop("statmod package required but is not installed")
g <- statmod::gauss.quad.prob(128,dist="uniform")
linkfun <- function(x) x/(1+x)
linkinv <- function(x) x/(1-x)
diff --git a/R/fitmixture.R b/R/fitmixture.R
new file mode 100644
index 0000000..0d3ba10
--- /dev/null
+++ b/R/fitmixture.R
@@ -0,0 +1,48 @@
+logcosh <- function(x)
+# Compute log(cosh(x)) without floating over or underflow
+# Gordon Smyth 8 April 2007
+{
+ y <- abs(x)-log(2)
+ i <- abs(x) < 1e-4
+ y[i] <- 0.5*x[i]^2
+ i <- !i & (abs(x) < 17)
+ y[i] <- log(cosh(x[i]))
+ y
+}
+
+fitmixture <- function(log2e,mixprop,niter=4,trace=FALSE)
+# Fit mixture model by non-linear least squares
+# Gordon Smyth 9 April 2007
+{
+ log2e <- as.matrix(log2e)
+ mixprop <- as.numeric(mixprop)
+ narrays <- ncol(log2e)
+ nprobes <- nrow(log2e)
+ y <- 2^log2e
+ start <- lm.fit(cbind(mixprop,1-mixprop),t(y))$coef
+ start <- pmax(start,1)
+ b <- as.vector(c(1,-1) %*% log(start))
+ z <- log(y)
+ pm <- matrix(1,nprobes,1) %*% matrix(2*mixprop-1,1,narrays)
+ for (i in 1:niter) {
+ mub <- logcosh(b/2)+log(1+tanh(b/2)*pm)
+ a <- rowMeans(z-mub)
+ mu <- a+mub
+ if(trace) {
+ s <- sqrt(narrays/(narrays-2)*rowMeans((z-mu)^2))
+ if(i>1) cat("stdev changes",summary(sold-s),"\n")
+ sold <- s
+ }
+ tb <- tanh(b/2)
+ dmu <- (tb+pm)/(1+tb*pm)/2
+ b <- b + rowMeans(dmu*(z-mu)) / (1e-6+rowMeans((dmu-rowMeans(dmu))^2))
+ if(trace) cat("betas",summary(b),"\n")
+ }
+ mub <- logcosh(b/2)+log(1+tanh(b/2)*pm)
+ a <- rowMeans(z-mub)
+ mu <- a+mub
+ s <- sqrt(narrays/(narrays-2)*rowMeans((z-mu)^2))
+ l2 <- log(2)
+ list(A=a/l2,M=b/l2,stdev=s/l2)
+}
+
diff --git a/R/geneset-camera.R b/R/geneset-camera.R
index 04449b2..8c9a3a5 100644
--- a/R/geneset-camera.R
+++ b/R/geneset-camera.R
@@ -21,10 +21,10 @@ interGeneCorrelation <- function(y, design)
camera <- function(y,...) UseMethod("camera")
-camera.default <- function(y,index,design=NULL,contrast=ncol(design),weights=NULL,use.ranks=FALSE,allow.neg.cor=TRUE,inter.gene.cor=NULL,trend.var=FALSE,sort=TRUE,...)
+camera.default <- function(y,index,design=NULL,contrast=ncol(design),weights=NULL,use.ranks=FALSE,allow.neg.cor=FALSE,inter.gene.cor=0.01,trend.var=FALSE,sort=TRUE,...)
# Competitive gene set test allowing for correlation between genes
# Gordon Smyth and Di Wu
-# Created 2007. Last modified 22 Dec 2015
+# Created 2007. Last modified 4 June 2016.
{
# Issue warning if extra arguments found
dots <- names(list(...))
@@ -59,7 +59,7 @@ camera.default <- function(y,index,design=NULL,contrast=ncol(design),weights=NUL
if(is.null(weights)) weights <- y$weights
# Check inter.gene.cor
- fixed.cor <- !is.null(inter.gene.cor)
+ fixed.cor <- !(is.na(inter.gene.cor) || is.null(inter.gene.cor))
# Set df for camera tests
if(fixed.cor) {
@@ -197,6 +197,9 @@ camera.default <- function(y,index,design=NULL,contrast=ncol(design),weights=NUL
tab$PValue <- tab$TwoSided
tab$Down <- tab$Up <- tab$TwoSided <- NULL
+# Remove correlation column if it was not estimated
+ if(fixed.cor) tab$Correlation <- NULL
+
# Add FDR
if(nsets>1) tab$FDR <- p.adjust(tab$PValue,method="BH")
diff --git a/R/geneset-fry.R b/R/geneset-fry.R
index f5260cd..18f5989 100644
--- a/R/geneset-fry.R
+++ b/R/geneset-fry.R
@@ -5,7 +5,7 @@ fry.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),geneid=NU
# The up and down p-values are equivalent to those from roast with nrot=Inf
# in the special case of prior.df=Inf.
# Gordon Smyth and Goknur Giner
-# Created 30 January 2015. Last modified 9 May 2016
+# Created 30 January 2015. Last modified 11 May 2016
{
# Partial matching of extra arguments
Dots <- list(...)
@@ -38,7 +38,7 @@ fry.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),geneid=NU
# Estimate genewise sds robustly
OK <- requireNamespace("statmod",quietly=TRUE)
- if(!OK) stop("statmod package required but not available")
+ if(!OK) stop("statmod package required but not installed")
gq <- statmod::gauss.quad.prob(128,"uniform")
df.residual <- ncol(Effects)-1
Eu2max <- sum( (df.residual+1)*gq$nodes^df.residual*qchisq(gq$nodes,df=1)*gq$weights )
@@ -46,7 +46,7 @@ fry.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),geneid=NU
s2.robust <- (rowSums(Effects^2)-u2max) / (df.residual+1-Eu2max)
# Empirical Bayes moderation method I: estimate hyperparameters from robust variances
- if(standardize=="posterior.sd") {
+ if(standardize=="p2") {
sv <- squeezeVar(s2.robust, df=0.92*df.residual,
covariate=covariate,
robust=Dots$robust,
@@ -54,8 +54,8 @@ fry.default <- function(y,index=NULL,design=NULL,contrast=ncol(design),geneid=NU
s2.robust <- sv$var.post
}
-# Empirical Bayes moderation method II: estimate hyperparameters from residual variances
- if(standardize=="p2") {
+# Empirical Bayes moderation method II: estimate hyperparameters from residual variances but apply squeezing to robust variances
+ if(standardize=="posterior.sd") {
s2 <- rowMeans(Effects[,-1,drop=FALSE]^2)
if(Dots$robust) {
fit <- fitFDistRobustly(s2, df1=df.residual, covariate=covariate, winsor.tail.p=Dots$winsor.tail.p)
diff --git a/R/goana.R b/R/goana.R
index 5319f41..f0fd1e3 100644
--- a/R/goana.R
+++ b/R/goana.R
@@ -114,33 +114,33 @@ goana.default <- function(de, universe = NULL, species = "Hs", prior.prob = NULL
# Get access to package of GO terms
suppressPackageStartupMessages(OK <- requireNamespace("GO.db",quietly=TRUE))
- if(!OK) stop("GO.db package required but not available")
+ if(!OK) stop("GO.db package required but not installed (or can't be loaded)")
# Get access to required annotation functions
suppressPackageStartupMessages(OK <- requireNamespace("AnnotationDbi",quietly=TRUE))
- if(!OK) stop("AnnotationDbi package required but not available")
+ if(!OK) stop("AnnotationDbi package required but not installed (or can't be loaded)")
# Get gene-GOterm mappings
switch(species,
Hs = {
suppressPackageStartupMessages(OK <- requireNamespace("org.Hs.eg.db",quietly=TRUE))
- if(!OK) stop("org.Hs.eg.db package required but not available")
+ if(!OK) stop("org.Hs.eg.db package required but not installed (or can't be loaded)")
GO2ALLEGS <- org.Hs.eg.db::org.Hs.egGO2ALLEGS
}, Mm = {
suppressPackageStartupMessages(OK <- requireNamespace("org.Mm.eg.db",quietly=TRUE))
- if(!OK) stop("org.Mm.eg.db package required but not available")
+ if(!OK) stop("org.Mm.eg.db package required but not installed (or can't be loaded)")
GO2ALLEGS <- org.Mm.eg.db::org.Mm.egGO2ALLEGS
}, Rn = {
suppressPackageStartupMessages(OK <- requireNamespace("org.Rn.eg.db",quietly=TRUE))
- if(!OK) stop("org.Rn.eg.db package required but not available")
+ if(!OK) stop("org.Rn.eg.db package required but not installed (or can't be loaded)")
GO2ALLEGS <- org.Rn.eg.db::org.Rn.egGO2ALLEGS
}, Dm = {
suppressPackageStartupMessages(OK <- requireNamespace("org.Dm.eg.db",quietly=TRUE))
- if(!OK) stop("org.Dm.eg.db package required but not available")
+ if(!OK) stop("org.Dm.eg.db package required but not installed (or can't be loaded)")
GO2ALLEGS <- org.Dm.eg.db::org.Dm.egGO2ALLEGS
}, Pt = {
suppressPackageStartupMessages(OK <- requireNamespace("org.Pt.eg.db",quietly=TRUE))
- if(!OK) stop("org.Pt.eg.db package required but not available")
+ if(!OK) stop("org.Pt.eg.db package required but not installed (or can't be loaded)")
GO2ALLEGS <- org.Pt.eg.db::org.Pt.egGO2ALLEGS
}
)
@@ -209,7 +209,7 @@ goana.default <- function(de, universe = NULL, species = "Hs", prior.prob = NULL
W <- AVE.PW*(Total-S[,"N"])/(PW.ALL-S[,"N"]*AVE.PW)
# Wallenius' noncentral hypergeometric test
- if(!requireNamespace("BiasedUrn",quietly=TRUE)) stop("BiasedUrn package required but is not available")
+ if(!requireNamespace("BiasedUrn",quietly=TRUE)) stop("BiasedUrn package required but is not installed (or can't be loaded)")
for(j in 1:nsets) for(i in 1:nrow(S))
P[i,j] <- BiasedUrn::pWNCHypergeo(S[i,1+j], S[i,"N"], Total-S[i,"N"], TotalDE[[j]], W[i],lower.tail=FALSE) + BiasedUrn::dWNCHypergeo(S[i,1+j], S[i,"N"], Total-S[i,"N"], TotalDE[[j]], W[i])
S <- S[,-ncol(S)]
diff --git a/R/loessFit.R b/R/loessFit.R
index 3d18152..4691704 100644
--- a/R/loessFit.R
+++ b/R/loessFit.R
@@ -81,7 +81,7 @@ loessFit <- function(y, x, weights=NULL, span=0.3, iterations=4L, min.weight=1e-
},
"locfit" = {
# Check for locfit package
- if(!requireNamespace("locfit",quietly=TRUE)) stop("locfit required but is not available")
+ if(!requireNamespace("locfit",quietly=TRUE)) stop("locfit required but is not installed (or can't be loaded)")
# Weighted lowess with robustifying iterations
biweights <- rep(1,nobs)
for (i in 1:iterations) {
diff --git a/R/norm.R b/R/norm.R
index 5c2f660..dbdf393 100755
--- a/R/norm.R
+++ b/R/norm.R
@@ -333,7 +333,12 @@ plotPrintorder <- function(object,layout,start="topleft",slide=1,method="loess",
normalizeBetweenArrays <- function(object, method=NULL, targets=NULL, cyclic.method="fast", ...) {
# Normalize between arrays
# Gordon Smyth
-# 12 Apri 2003. Last revised 21 July 2013.
+# 12 Apri 2003. Last revised 13 June 2016.
+
+ if(is.data.frame(object)) {
+ object <- as.matrix(object)
+ if(mode(object) != "numeric") stop("'object' is a data.frame and not all columns are numeric")
+ }
# Default method
if(is.null(method)) {
diff --git a/R/read.idat.R b/R/read.idat.R
index 8e040b5..4c5f607 100644
--- a/R/read.idat.R
+++ b/R/read.idat.R
@@ -1,43 +1,88 @@
-read.idat <- function(idatfiles, bgxfile, dateinfo=FALSE, tolerance=0)
-# Read idat data from gene expression BeadChips
-# Matt Ritchie
-# Created 30 September 2013. Modified on 14 January 2015 and 26 June 2015.
+read.idat <- function(idatfiles, bgxfile, dateinfo=FALSE, annotation="Symbol", tolerance=0L, verbose=TRUE)
+# Read GenomeStudio IDAT files for Illumina gene expression BeadChips
+# Matt Ritchie and Gordon Smyth
+# Created 30 September 2013. Last modified 14 June 2016.
{
- if(!requireNamespace("illuminaio",quietly=TRUE)) stop("illuminaio package required but is not available")
- cat("Reading manifest file", bgxfile, "... ")
- bgx <- illuminaio::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 <- illuminaio::readIDAT(idatfiles[i])
- cat("Done\n")
- ind <- match(elist$genes[,"Array_Address_Id"], tmp$Quants[,'IllumicodeBinData'])
- if(sum(is.na(ind))>tolerance)
- stop("Can't match all ids in manifest with those in idat file ", idatfiles[i], "\n", sum(is.na(ind)),
- " missing - please check that you have the right files, or consider setting \'tolerance\'=", sum(is.na(ind)))
- elist$E[!is.na(ind),i] <- round(tmp$Quants[ind[!is.na(ind)],'MeanBinData'],1) # intensity data
- elist$other$STDEV[!is.na(ind),i] <- tmp$Quants[ind[!is.na(ind)],'DevBinData'] # Bead STDEV
- elist$other$NumBeads[!is.na(ind),i] <- tmp$Quants[ind[!is.na(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)
+# Need illuminaio package
+ OK <- requireNamespace("illuminaio",quietly=TRUE)
+ if(!OK) stop("illuminaio package required but is not installed (or can't be loaded)")
+
+# Check idatfiles
+ idatafiles <- as.character(idatfiles)
+ nsamples <- length(idatfiles)
+
+# Initialize EListRaw object
+ elist <- new("EListRaw")
+ elist$source <- "illumina"
+ elist$targets <- data.frame("IDATfile"=idatfiles,stringsAsFactors=FALSE)
+ if(dateinfo) elist$targets$DecodeInfo <- elist$targets$ScanInfo <- rep_len("", nsamples)
+
+# Read bead manifest file
+ if(verbose) cat("Reading manifest file", bgxfile, "... ")
+ bgx <- illuminaio::readBGX(bgxfile)
+ if(verbose) cat("Done\n")
+ nregprobes <- nrow(bgx$probes)
+ nctrlprobes <- nrow(bgx$control)
+ nprobes <- nregprobes+nctrlprobes
+
+# Assemble gene annotation
+ elist$genes <- rbind(bgx$probes[,c("Probe_Id","Array_Address_Id")], bgx$controls[,c("Probe_Id","Array_Address_Id")])
+
+# Set probe control status
+ elist$genes$Status <- "regular"
+ elist$genes$Status[(nregprobes+1):nprobes] <- bgx$controls[,"Reporter_Group_Name"]
+
+# Add optional annotation columns
+ if(!is.null(annotation) && !is.na(annotation)) {
+ annotation <- as.character(annotation)
+ annotation <- intersect(names(bgx$probes),annotation)
+ }
+ if(length(annotation)) {
+ ac <- annotation %in% names(bgx$controls)
+ for (i in 1:length(annotation)) {
+ elist$genes[[annotation[i]]] <- NA_character_
+ elist$genes[[annotation[i]]][1:nregprobes] <- bgx$probes[[annotation[i]]]
+ if(ac[i]) elist$genes[[annotation[i]]][(nregprobes+1L):nprobes] <- bgx$controls[[annotation[i]]]
+ }
+ }
+
+# Initalize expression matrices
+ elist$E <- matrix(0, nprobes, nsamples)
+ elist$E[] <- NA
+ colnames(elist$E) <- removeExt(idatfiles)
+ rownames(elist$E) <- elist$genes[,"Array_Address_Id"]
+ elist$other$STDEV <- elist$other$NumBeads <- elist$E
+
+# Read IDAT files
+ for(j in 1:nsamples) {
+ if(verbose) cat("\t", idatfiles[j], "... ")
+ tmp <- illuminaio::readIDAT(idatfiles[j])
+ if(verbose) cat("Done\n")
+ ind <- match(elist$genes$Array_Address_Id, tmp$Quants$IllumicodeBinData)
+
+# Check for whether values are available for all probes
+ if(anyNA(ind)) {
+ nna <- sum(is.na(ind))
+ if(nna > tolerance)
+ stop("Can't match all ids in manifest with those in idat file ", idatfiles[i], "\n", sum(is.na(ind)),
+ " missing - please check that you have the right files, or consider setting \'tolerance\'=", sum(is.na(ind)))
+ i <- which(!is.na(ind))
+ ind <- ind[i]
+ elist$E[i,j] <- tmp$Quants$MeanBinData[ind]
+ elist$other$STDEV[i,j] <- tmp$Quants$DevBinData[ind]
+ elist$other$NumBeads[i,j] <- tmp$Quants$NumGoodBeadsBinData[ind]
+ } else {
+ elist$E[,j] <- tmp$Quants$MeanBinData[ind]
+ elist$other$STDEV[,j] <- tmp$Quants$DevBinData[ind]
+ elist$other$NumBeads[,j] <- tmp$Quants$NumGoodBeadsBinData[ind]
+ }
+
+ if(dateinfo) {
+ elist$targets$DecodeInfo[j] = paste(tmp$RunInfo[1,], collapse=" ")
+ elist$targets$ScanInfo[j] = paste(tmp$RunInfo[2,], collapse=" ")
+ }
+ }
+
+ if(verbose) cat("Finished reading data.\n")
+ return(elist)
}
diff --git a/R/sepchannel.R b/R/sepchannel.R
index 1991562..cd0305c 100755
--- a/R/sepchannel.R
+++ b/R/sepchannel.R
@@ -87,7 +87,7 @@ intraspotCorrelation <- function(object,design,trim=0.15)
designA <- (Ident %x% matrix(c(0.5,0.5),1,2)) %*% design
X <- rbind(designM, designA)
Z <- diag(2) %x% rep(1,narrays)
- if(!requireNamespace("statmod",quietly=TRUE)) stop("statmod package required but is not available")
+ if(!requireNamespace("statmod",quietly=TRUE)) stop("statmod package required but is not installed")
arho <- rep(NA,ngenes)
degfre <- matrix(0,ngenes,2,dimnames=list(rownames(M),c("df.M","df.A")))
for (i in 1:ngenes) {
diff --git a/inst/doc/changelog.txt b/inst/doc/changelog.txt
index 500b25b..66c0037 100755
--- a/inst/doc/changelog.txt
+++ b/inst/doc/changelog.txt
@@ -1,3 +1,42 @@
+15 Jun 2016: limma 3.28.10
+
+- new argument 'annotation' for read.idat(), to allow users to read
+ any required columns from the manifest file.
+
+15 Jun 2016: limma 3.28.9
+
+- normalizeBetweenArrays() now checks whether the input 'object' is a
+ data.frame, and converts to a matrix if possible.
+
+13 Jun 2016: limma 3.28.8
+
+- new function detectionPValues() to compute detection p-values from
+ negative control probes.
+
+12 Jun 2016: limma 3.28.7
+
+- new argument 'verbose' for read.idat().
+
+- Example code added to duplicateCorrelation.Rd.
+
+ 8 Jun 2016: limma 3.28.6
+
+- The default settings for the 'inter.gene.cor' and 'use.neg.cor'
+ arguments of camera() have been changed. camera() now uses by
+ default a preset value for the inter-gene correlation. This has the
+ effect that it tends to rank highly co-regulated,
+ biologically-interpretable more highly than before.
+
+22 May 2016: limma 3.28.5
+
+- New function fitmixture() to estimate genewise fold changes and
+ expression values from mixture experiments. Previously this
+ funtion was only available as part of the illumina package available
+ from http://bioinf.wehi.edu.au/illumina
+
+- New function logcosh() to compute log(cosh(x)) accurately without
+ floating underflow or overflow.
+
11 May 2016: limma 3.28.4
- Slight change to fry() standardize="posterior.sd" method.
@@ -22,7 +61,7 @@
4 May 2016: limma 3.28.1
-- unnecesary argument 'design' removed from fitted.MArrayLM.
+- unnecessary argument 'design' removed from fitted.MArrayLM.
- bug fix to fry() with robust=TRUE.
diff --git a/inst/doc/intro.pdf b/inst/doc/intro.pdf
index ab7d1b2..c3ef58e 100644
Binary files a/inst/doc/intro.pdf and b/inst/doc/intro.pdf differ
diff --git a/man/03reading.Rd b/man/03reading.Rd
index 1eb4283..e55cacd 100644
--- a/man/03reading.Rd
+++ b/man/03reading.Rd
@@ -26,6 +26,7 @@ There are also a series of utility functions which read the header information f
and produces an \code{ElistRaw} object.
\code{\link{read.idat}} reads Illumina files in IDAT format, and produces an \code{EListRaw} object.
+\code{\link{detectionPValues}} can be used to add detection p-values.
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/10GeneSetTests.Rd b/man/10GeneSetTests.Rd
index 6670298..116dfab 100644
--- a/man/10GeneSetTests.Rd
+++ b/man/10GeneSetTests.Rd
@@ -13,7 +13,7 @@ This page gives an overview of the LIMMA functions for gene set testing and path
Self-contained gene set testing for many sets.}
\item{ \code{\link{fry}} }{
- Fast version of \code{mroast} when heteroscedasticity of genes can be ignored.}
+ Fast approximation to \code{mroast}, especially useful when heteroscedasticity of genes can be ignored.}
\item{ \code{\link{camera}} }{
Competitive gene set testing.}
diff --git a/man/camera.Rd b/man/camera.Rd
index 93def95..ce35342 100644
--- a/man/camera.Rd
+++ b/man/camera.Rd
@@ -7,9 +7,9 @@
Test whether a set of genes is highly ranked relative to other genes in terms of differential expression, accounting for inter-gene correlation.
}
\usage{
-\S3method{camera}{default}(y, index, design, contrast=ncol(design), weights=NULL,
- use.ranks=FALSE, allow.neg.cor=TRUE, inter.gene.cor=NULL, trend.var=FALSE,
- sort=TRUE, \dots)
+\method{camera}{default}(y, index, design, contrast = ncol(design), weights = NULL,
+ use.ranks = FALSE, allow.neg.cor=FALSE, inter.gene.cor=0.01, trend.var = FALSE,
+ sort = TRUE, \dots)
interGeneCorrelation(y, design)
}
\arguments{
@@ -23,7 +23,7 @@ interGeneCorrelation(y, 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)}, 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{inter.gene.cor}{numeric, optional preset value for the inter gene correlation within tested sets. If not \code{NULL}, then this value will over-ride estimation of the inter gene correlations.}
+ \item{inter.gene.cor}{numeric, optional preset value for the inter-gene correlation within tested sets. If \code{NA} or \code{NULL}, then an inter-gene correlation will be estimated for each tested set.}
\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}
@@ -46,20 +46,30 @@ The inflation factor depends on estimated genewise correlation and the number of
By default, \code{camera} uses \code{interGeneCorrelation} to estimate the mean pair-wise correlation within each set of genes.
\code{camera} can alternatively be used with a preset correlation specified by \code{inter.gene.cor} that is shared by all sets.
This usually works best with a small value, say \code{inter.gene.cor=0.01}.
-This produces a less conservative test that performs well in many practical situations.
+
+If \code{interGeneCorrelation=NA}, then \code{camera} will estimate the inter-gene correlation for each set.
+In this mode, \code{camera} gives rigorous error rate control for all sample sizes and all gene sets.
+However, in this mode, highly co-regulated gene sets that are biological interpretable may not always be ranked at the top of the list.
+
+With \code{interGeneCorrelation=0.01}, \code{camera} will rank biologically interpetable sets more highly.
+This gives a useful compromise between strict error rate control and interpretable gene set rankings.
}
+\note{The default settings for \code{inter.gene.cor} and \code{allow.neg.cor} were changed to the current values in limma 3.29.6.
+Previously, the default was to estimate an inter-gene correlation for each set.
+To reproduce the earlier default, use \code{allow.neg.cor=TRUE} and \code{inter.gene.cor=NA}.}
+
\value{
\code{camera} returns a data.frame with a row for each set and the following columns:
-\item{NGenes}{number of genes in set}
-\item{Correlation}{inter-gene correlation}
-\item{Direction}{direction of change (\code{"Up"} or \code{"Down"})}
-\item{PValue}{two-tailed p-value}
-\item{FDR}{Benjamini and Hochberg FDR adjusted p-value}
+\item{NGenes}{number of genes in set.}
+\item{Correlation}{inter-gene correlation (only included if the \code{inter.gene.cor} was not preset).}
+\item{Direction}{direction of change (\code{"Up"} or \code{"Down"}).}
+\item{PValue}{two-tailed p-value.}
+\item{FDR}{Benjamini and Hochberg FDR adjusted p-value.}
\code{interGeneCorrelation} returns a list with components:
-\item{vif}{variance inflation factor}
-\item{correlation}{inter-gene correlation}
+\item{vif}{variance inflation factor.}
+\item{correlation}{inter-gene correlation.}
}
\author{Di Wu and Gordon Smyth}
@@ -80,6 +90,7 @@ Analyzing gene expression data in terms of gene sets: methodological issues.
\code{\link{rankSumTestWithCorrelation}},
\code{\link{geneSetTest}},
\code{\link{roast}},
+\code{\link{fry}},
\code{\link{romer}},
\code{\link{ids2indices}}.
@@ -100,7 +111,7 @@ index2 <- 21:40
camera(y, index1, design)
camera(y, index2, design)
-camera(y, list(set1=index1,set2=index2), design)
+camera(y, list(set1=index1,set2=index2), design, inter.gene.cor=NA)
camera(y, list(set1=index1,set2=index2), design, inter.gene.cor=0.01)
}
diff --git a/man/detectionPValue.Rd b/man/detectionPValue.Rd
new file mode 100644
index 0000000..aca5353
--- /dev/null
+++ b/man/detectionPValue.Rd
@@ -0,0 +1,58 @@
+\name{detectionPValues}
+\alias{detectionPValues.default}
+\alias{detectionPValues.EListRaw}
+\alias{detectionPValues}
+
+\title{Detection P-Values from Negative Controls}
+
+\description{Compute the proportion of negative controls greater than each observed expression value.
+Particularly useful for Illumina BeadChips.}
+
+\usage{
+\method{detectionPValues}{EListRaw}(x, status = NULL, \dots)
+\method{detectionPValues}{default}(x, status, negctrl = "negative", \dots)
+}
+
+\arguments{
+ \item{x}{object of class \code{EListRaw} or a numeric \code{matrix} containing raw intensities for regular and control probes from a series of microarrays.}
+ \item{status}{character vector giving probe types. Defaults to \code{x$genes$Status} if \code{x} is an \code{EListRaw} object.}
+ \item{negctrl}{character string identifier for negative control probes.}
+ \item{\dots}{other arguments are not currently used.}
+ }
+
+\details{
+The rows of \code{x} for which \code{status == negctrl} are assumed to correspond to negative control probes.
+
+For each column of \code{x}, the detection p-values are defined as \code{(N.eq/2 + N.gt) / N.neg}, where \code{N.gt} is the number of negative controls with expression greater than the observed value, \code{N.eq} is the number of negative controls with expression equal to the observed value, and \code{N.neg} is the total number of negative controls.
+}
+
+\value{numeric matrix of same dimensions as \code{x} containing detection p-values.}
+
+\references{
+Shi, W, de Graaf, C, Kinkel, S, Achtman, A, Baldwin, T, Schofield, L, Scott, H, Hilton, D, Smyth, GK (2010).
+Estimating the proportion of microarray probes expressed in an RNA sample.
+\emph{Nucleic Acids Research} 38, 2168-2176.
+\url{http://nar.oxfordjournals.org/content/38/7/2168}
+}
+
+\author{Gordon Smyth}
+
+\seealso{
+ An overview of LIMMA functions to read expression data is given in \link{03.ReadingData}.
+
+ \code{\link{read.idat}} reads Illumina BeadChip expression data from binary IDAT files.
+
+ \code{\link{neqc}} performs normexp background correction and quantile normalization aided by control probes.
+}
+
+\examples{
+\dontrun{
+# Read Illumina binary IDAT files
+x <- read.idat(idat, bgx)
+x$genes$DectionPValue <- detectionPValues(x)
+y <- neqc(x)
+}
+}
+
+\keyword{background correction}
+\keyword{illumina beadchips}
diff --git a/man/dupcor.Rd b/man/dupcor.Rd
index 6054cc9..ee170ad 100755
--- a/man/dupcor.Rd
+++ b/man/dupcor.Rd
@@ -52,16 +52,40 @@ An overview of linear model functions in limma is given by \link{06.LinearModels
\author{Gordon Smyth}
\references{
Smyth, G. K., Michaud, J., and Scott, H. (2005). The use of within-array replicate spots for assessing differential expression in microarray experiments. \emph{Bioinformatics} 21(9), 2067-2075.
-\url{http://www.statsci.org/smyth/pubs/dupcor.pdf}
+[\url{http://bioinformatics.oxfordjournals.org/content/21/9/2067}]
+[Preprint with corrections: \url{http://www.statsci.org/smyth/pubs/dupcor.pdf}]
}
+
\examples{
-# Also see lmFit examples
+# Simulate gene expression data for 100 probes and 6 microarrays
+# Microarray are in two groups
+# First two probes are more highly expressed in second group
+# Std deviations vary between genes with prior df=4
+sd <- 0.3*sqrt(4/rchisq(100,df=4))
+y <- matrix(rnorm(100*6,sd=sd),100,6)
+rownames(y) <- paste("Gene",1:100)
+y[1:2,4:6] <- y[1:2,4:6] + 2
+design <- cbind(Grp1=1,Grp2vs1=c(0,0,0,1,1,1))
+options(digits=3)
-\dontrun{
-corfit <- duplicateCorrelation(MA, ndups=2, design)
-all.correlations <- tanh(corfit$atanh.correlations)
-boxplot(all.correlations)
-fit <- lmFit(MA, design, ndups=2, correlation=corfit$consensus)
-}
+# Fit with correlated arrays
+# Suppose each pair of arrays is a block
+block <- c(1,1,2,2,3,3)
+dupcor <- duplicateCorrelation(y,design,block=block)
+dupcor$consensus.correlation
+fit1 <- lmFit(y,design,block=block,correlation=dupcor$consensus)
+fit1 <- eBayes(fit1)
+topTable(fit1,coef=2)
+
+# Fit with duplicate probes
+# Suppose two side-by-side duplicates of each gene
+rownames(y) <- paste("Gene",rep(1:50,each=2))
+dupcor <- duplicateCorrelation(y,design,ndups=2)
+dupcor$consensus.correlation
+fit2 <- lmFit(y,design,ndups=2,correlation=dupcor$consensus)
+dim(fit2)
+fit2 <- eBayes(fit2)
+topTable(fit2,coef=2)
}
+
\keyword{multivariate}
diff --git a/man/fitmixture.Rd b/man/fitmixture.Rd
new file mode 100644
index 0000000..5a7db5b
--- /dev/null
+++ b/man/fitmixture.Rd
@@ -0,0 +1,66 @@
+\name{fitmixture}
+\alias{fitmixture}
+\title{Fit Mixture Model by Non-Linear Least Squares}
+\description{Fit Mixture Model by Non-Linear Least Squares}
+
+\usage{
+fitmixture(log2e, mixprop, niter = 4, trace = FALSE)
+}
+
+\arguments{
+ \item{log2e}{a numeric matrix containing log2 expression values. Rows correspond to probes for genes and columns to RNA samples.}
+ \item{mixprop}{a vector of length \code{ncol(log2e)} giving the mixing proportion (between 0 and 1) for each sample.}
+ \item{niter}{integer number of iterations.}
+ \item{trace}{logical. If \code{TRUE}, summary working estimates are output from each iteration.}
+}
+
+\details{
+A mixture experiment is one in which two reference RNA sources are mixed in different proportions to create experimental samples.
+Mixture experiments have been used to evaluate genomic technologies and analysis methods (Holloway et al, 2006).
+This function uses all the data for each gene to estimate the expression level of the gene in each of two pure samples.
+
+The function fits a nonlinear mixture model to the log2 expression values for each gene.
+The expected values of \code{log2e} for each gene are assumed to be of the form
+\code{log2( mixprop*Y1 + (1-mixprop)*Y2 )}
+where \code{Y1} and \code{Y2} are the expression levels of the gene in the two reference samples being mixed.
+The \code{mixprop} values are the same for each gene but \code{Y1} and \code{Y2} are specific to the gene.
+The function returns the estimated values \code{A=0.5*log2(Y1*Y2)} and \code{M=log2(Y2/Y1)} for each gene.
+
+The nonlinear estimation algorithm implemented in \code{fitmixture} uses a nested Gauss-Newton iteration (Smyth, 1996).
+It is fully vectorized so that the estimation is done for all genes simultaneously.
+}
+
+\value{List with three components:
+\item{A}{numeric vector giving the estimated average log2 expression of the two reference samples for each gene}
+\item{M}{numeric vector giving estimated log-ratio of expression between the two reference samples for each gene}
+\item{stdev}{standard deviation of the residual term in the mixture model for each gene}
+}
+
+\author{Gordon K Smyth}
+
+\references{
+Holloway, A. J., Oshlack, A., Diyagama, D. S., Bowtell, D. D. L., and Smyth, G. K. (2006).
+Statistical analysis of an RNA titration series evaluates microarray precision and sensitivity on a whole-array basis.
+\emph{BMC Bioinformatics} 7, Article 511.
+\url{http://www.biomedcentral.com/1471-2105/7/511}
+
+Smyth, G. K. (1996).
+Partitioned algorithms for maximum likelihood and other nonlinear estimation.
+\emph{Statistics and Computing}, 6, 201-216.
+\url{http://www.statsci.org/smyth/pubs/partitio.pdf}
+}
+
+\examples{
+ngenes <- 100
+TrueY1 <- rexp(ngenes)
+TrueY2 <- rexp(ngenes)
+mixprop <- matrix(c(0,0.25,0.75,1),1,4)
+TrueExpr <- TrueY1 %*% mixprop + TrueY2 %*% (1-mixprop)
+
+log2e <- log2(TrueExpr) + matrix(rnorm(ngenes*4),ngenes,4)*0.1
+out <- fitmixture(log2e,mixprop)
+
+# Plot true vs estimated log-ratios
+plot(log2(TrueY1/TrueY2), out$M)
+}
+
diff --git a/man/logcosh.Rd b/man/logcosh.Rd
new file mode 100644
index 0000000..d64cfaa
--- /dev/null
+++ b/man/logcosh.Rd
@@ -0,0 +1,26 @@
+\name{logcosh}
+\alias{logcosh}
+\title{Logarithm of cosh}
+\description{Compute \code{log(cosh(x))} without floating overflow or underflow}
+\usage{
+logcosh(x)
+}
+\arguments{
+ \item{x}{a numeric vector or matrix.}
+}
+
+\details{
+The computation uses asymptotic expressions for very large or very small arguments.
+For intermediate arguments, \code{log(cosh(x))} is returned.
+}
+
+\value{Numeric vector or matrix of same dimensions as \code{x}.}
+
+\author{Gordon K Smyth}
+
+\examples{
+x <- c(1e-8,1e-7,1e-6,1e-5,1e-4,1,3,50,800)
+logcosh(x)
+log(cosh(x))
+}
+
diff --git a/man/propexpr.Rd b/man/propexpr.Rd
index 6cb3c70..f3ce8b6 100644
--- a/man/propexpr.Rd
+++ b/man/propexpr.Rd
@@ -8,21 +8,25 @@ propexpr(x, neg.x=NULL, status=x$genes$Status, labels=c("negative","regular"))
\arguments{
\item{x}{matrix or similar object containing raw intensities for a set of arrays.}
\item{neg.x}{matrix or similar object containing raw intensities for negative control probes for the same arrays. If \code{NULL}, then negative controls must be provided in \code{x}.}
- \item{status}{character vector giving probe types.}
- \item{labels}{character vector giving probe type identifiers.}
+ \item{status}{character vector specifying control type of each probe. Only used if \code{neg.x} is \code{NULL}.}
+ \item{labels}{character vector giving the \code{status} values for negative control probes and regular (non-control) probes respectively. If of length 1, then all probes other than the negative controls are assumed to be regular. Only used if \code{neg.x} is \code{NULL}.}
}
+
\details{
-This function estimates the proportion of expressed in a microarray by utilizing the negative control probes.
+This function estimates the overall proportion of probes on each microarray that are correspond to expressed genes using the method of Shi et al (2010).
+The function is especially useful for Illumina BeadChips arrays, although it can in principle be applied to any platform with good quality negative controls.
+
+The negative controls can be supplied either as rows of \code{x} or as a separate matrix.
+If supplied as rows of \code{x}, then the negative controls are identified by the \code{status} vector.
+\code{x} might also include other types of control probes, but these will be ignored in the calculation.
+
Illumina BeadChip arrays contain 750~1600 negative control probes.
-The expression profile of these control probes can be saved to a separate file by the Illumina BeadStudio software when using it to output the expression profile for regular probes.
-The control probe profile could be re-generated if it was not generated when the regular probe profile was created by BeadStudio.
-Other microarray platforms can also use this function to estimate the proportion of expressed probes in each array, provided that they have a set of negative control probes.
-
-\code{labels} can include one or two probe type identifiers.
-Its first element should be the identifier for negative control probes (\code{negative} by default).
-If \code{labels} only contains one identifier, then it will be assumed to contain the identifier for negative control probes.
-By default, \code{regular} is the identifier for regular probes.
+If \code{read.idat} is used to read Illumina expression IDAT files, then the control probes will be populated as rows of the output \code{EListRaw} object, and the vector \code{x$genes$Status} will be set to identify control probes.
+
+Alternatively, expression values can be exported from Illumina's GenomeStudio software as tab-delimited text files.
+In this case, the control probes are usually written to a separate file from the regular probes.
}
+
\value{
Numeric vector giving the proportions of expressed probes in each array.
}
@@ -42,7 +46,16 @@ Description to the control probes in Illumina BeadChips can be found in \code{\l
\examples{
\dontrun{
-x <- read.ilmn(files="sample probe profile.txt",ctrlfiles="control probe profile.txt")
+# Read Illumina binary IDAT files
+x <- read.idat(idat, bgx)
+propexpr(x)
+
+# Read text files exported from GenomeStudio
+x <- read.ilmn(files = "sample probe profile.txt",
+ ctrlfiles = "control probe profile.txt")
propexpr(x)
}
}
+
+\keyword{background correction}
+\keyword{illumina beadchips}
diff --git a/man/read.idat.Rd b/man/read.idat.Rd
index a2e247d..39f66a7 100644
--- a/man/read.idat.Rd
+++ b/man/read.idat.Rd
@@ -6,64 +6,93 @@
\description{Read Illumina BeadArray data from IDAT and manifest (.bgx) files for gene expression platforms.}
\usage{
-read.idat(idatfiles, bgxfile, dateinfo=FALSE, tolerance=0)
+read.idat(idatfiles, bgxfile, dateinfo = FALSE, annotation = "Symbol",
+ tolerance = 0L, verbose = TRUE)
}
\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?}
- \item{tolerance}{numeric. The number of probe ID discrepancies allowed between the manifest and a given idat file.}
+ \item{dateinfo}{logical. Should date and software version information be read in?}
+ \item{annotation}{character vector of annotation columns to be read from the manifest file.}
+ \item{tolerance}{integer. The number of probe ID discrepancies allowed between the manifest and any of the IDAT files.}
+ \item{verbose}{logical. Should progress messages are sent to standard output?}
}
\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).
+ Illumina's BeadScan/iScan software outputs probe intensities in IDAT
+ format (encrypted XML files) and uses probe information stored 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}.
+ 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 to store them in an \code{EListRaw-class} object.
+ The function serves a similar purpose to \code{\link{read.ilmn}},
+ which reads text files exported by Illumina's GenomeStudio software,
+ but it reads the IDAT files directly without any need to convert them first to text.
+
+ The function reads information on control probes as well for regular probes.
Probe types are indicated in the \code{Status} column of the \code{genes}
- component of the \code{EListRaw-class} object.
+ component of the \code{EListRaw} object.
+
+ The \code{annotation} argument specifies probe annotation columns to be extracted from the manifest file.
+ The manifest typically contains the following columns:
+ \code{"Species"}, \code{"Source"}, \code{"Search_Key"}, \code{"Transcript"},
+ \code{"ILMN_Gene"}, \code{"Source_Reference_ID"}, \code{"RefSeq_ID"},
+ \code{"Unigene_ID"}, \code{"Entrez_Gene_ID"}, \code{"GI"},
+ \code{"Accession"}, \code{"Symbol"}, \code{"Protein_Product"},
+ \code{"Probe_Id"}, \code{"Array_Address_Id"}, \code{"Probe_Type"},
+ \code{"Probe_Start"}, \code{"Probe_Sequence"}, \code{"Chromosome"},
+ \code{"Probe_Chr_Orientation"}, \code{"Probe_Coordinates"}, \code{"Cytoband"},
+ \code{"Definition"}, \code{"Ontology_Component"}, \code{"Ontology_Process"},
+ \code{"Ontology_Function"}, \code{"Synonyms"}, \code{"Obsolete_Probe_Id"}.
+ Note that \code{"Probe_Id"} and \code{"Array_Address_Id"} are always extracted and
+ do not need to included in the \code{annotation} argument.
+
+ If more than \code{tolerance} probes in the manifest cannot be found in an IDAT file then the function will return an error.
}
\value{
- An \code{EListRaw-class} object with the following components:
+ An \code{EListRaw} 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.}
+ \item{other$NumBeads}{ numeric matrix of same dimensions as \code{E} giving number of beads used for each intensity value.}
+ \item{other$STDEV}{ numeric matrix of same dimensions as \code{E} giving bead-level standard deviation or standard error for each intensity value.}
+ \item{genes}{ data.frame of probe annotation.
+ This includes the \code{Probe_Id} and \code{Array_Address_Id} columns extracted from the manifest file,
+ plus a \code{Status} column identifying control probes,
+ plus any other columns specified by \code{annotation}.}
+ \item{targets}{ data.frame of sample information.
+ This includes the IDAT file names plus other columns if \code{dateinfo=TRUE}.}
}
\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}
+\url{http://f1000research.com/articles/2-264/}
}
\author{Matt Ritchie}
\seealso{
- \code{read.ilmn} imports gene expression data output by GenomeStudio.
+ \code{\link{read.ilmn}} imports gene expression data output by GenomeStudio.
- \code{neqc} performs normexp by control background correction, log
+ \code{\link{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.
+ \code{\link{propexpr}} estimates the proportion of expressed probes in a microarray.
+
+ \code{\link{detectionPValues}} computes detection p-values from the negative controls.
}
\examples{
\dontrun{
-idatfiles = dir(pattern="idat")
-bgxfile = dir(pattern="bgx")
-data = read.idat(idatfiles, bgxfile)
+idatfiles <- dir(pattern="idat")
+bgxfile <- dir(pattern="bgx")
+x <- read.idat(idatfiles, bgxfile)
+x$other$Detection <- detectionPValues(x)
propexpr(data)
-datanorm = neqc(data)
+y <- neqc(data)
}
}
diff --git a/man/roast.Rd b/man/roast.Rd
index 51e50c9..fd336dc 100644
--- a/man/roast.Rd
+++ b/man/roast.Rd
@@ -13,14 +13,14 @@ Rotation gene set testing for linear models.
}
\usage{
-\S3method{roast}{default}(y, index = NULL, design = NULL, contrast = ncol(design), geneid = NULL,
+\method{roast}{default}(y, index = NULL, design = NULL, contrast = ncol(design), geneid = NULL,
set.statistic = "mean", gene.weights = NULL, var.prior = NULL, df.prior = NULL,
nrot = 999, approx.zscore = TRUE, \dots)
-\S3method{mroast}{default}(y, index = NULL, design = NULL, contrast = ncol(design), geneid = NULL,
+\method{mroast}{default}(y, index = NULL, design = NULL, contrast = ncol(design), geneid = NULL,
set.statistic = "mean", gene.weights = NULL, var.prior = NULL, df.prior = NULL,
nrot = 999, approx.zscore = TRUE, adjust.method = "BH",
midp = TRUE, sort = "directional", \dots)
-\S3method{fry}{default}(y, index = NULL, design = NULL, contrast = ncol(design), geneid = NULL,
+\method{fry}{default}(y, index = NULL, design = NULL, contrast = ncol(design), geneid = NULL,
standardize = "posterior.sd", sort = "directional", \dots)
}
diff --git a/man/romer.Rd b/man/romer.Rd
index 12cb45f..000ad2f 100644
--- a/man/romer.Rd
+++ b/man/romer.Rd
@@ -7,7 +7,7 @@
Gene set enrichment analysis for linear models using rotation tests (ROtation testing using MEan Ranks).
}
\usage{
-\S3method{romer}{default}(y, index, design = NULL, contrast = ncol(design),
+\method{romer}{default}(y, index, design = NULL, contrast = ncol(design),
array.weights = NULL, block = NULL, correlation,
set.statistic = "mean", nrot = 9999, shrink.resid = TRUE, \dots)
}
diff --git a/tests/limma-Tests.R b/tests/limma-Tests.R
index 1bcf9d1..974fff2 100755
--- a/tests/limma-Tests.R
+++ b/tests/limma-Tests.R
@@ -231,9 +231,12 @@ mroast(y=y,list(set1=iset1,set2=iset2),design,contrast=2,gene.weights=runif(100)
mroast(y=y,list(set1=iset1,set2=iset2),design,contrast=2,array.weights=c(0.5,1,0.5,1))
mroast(y=y,list(set1=iset1,set2=iset2),design,contrast=2,weights=w)
mroast(y=y,list(set1=iset1,set2=iset2),design,contrast=2,weights=w,array.weights=c(0.5,1,0.5,1))
+fry(y=y,list(set1=iset1,set2=iset2),design,contrast=2,weights=w,array.weights=c(0.5,1,0.5,1))
### camera
+camera(y=y,iset1,design,contrast=2,weights=c(0.5,1,0.5,1),allow.neg.cor=TRUE,inter.gene.cor=NA)
+camera(y=y,list(set1=iset1,set2=iset2),design,contrast=2,allow.neg.cor=TRUE,inter.gene.cor=NA)
camera(y=y,iset1,design,contrast=2,weights=c(0.5,1,0.5,1))
camera(y=y,list(set1=iset1,set2=iset2),design,contrast=2)
@@ -241,6 +244,7 @@ camera(y=y,list(set1=iset1,set2=iset2),design,contrast=2)
y <- new("EList",list(E=y))
roast(y=y,iset1,design,contrast=2)
+camera(y=y,iset1,design,contrast=2,allow.neg.cor=TRUE,inter.gene.cor=NA)
camera(y=y,iset1,design,contrast=2)
### eBayes with trend
diff --git a/tests/limma-Tests.Rout.save b/tests/limma-Tests.Rout.save
index dbff537..f085afd 100755
--- a/tests/limma-Tests.Rout.save
+++ b/tests/limma-Tests.Rout.save
@@ -1,5 +1,5 @@
-R Under development (unstable) (2016-03-14 r70331) -- "Unsuffered Consequences"
+R version 3.3.0 (2016-05-03) -- "Supposedly Educational"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)
@@ -1191,16 +1191,27 @@ set2 5 0.2 0 Down 0.950 0.950 0.250 0.250
NGenes PropDown PropUp Direction PValue FDR PValue.Mixed FDR.Mixed
set1 5 0 1 Up 0.001 0.001 0.001 0.001
set2 5 0 0 Down 0.791 0.791 0.146 0.146
+> fry(y=y,list(set1=iset1,set2=iset2),design,contrast=2,weights=w,array.weights=c(0.5,1,0.5,1))
+ NGenes Direction PValue FDR PValue.Mixed FDR.Mixed
+set1 5 Up 0.0007432594 0.001486519 1.820548e-05 3.641096e-05
+set2 5 Down 0.8208140511 0.820814051 2.211837e-01 2.211837e-01
>
> ### camera
>
-> camera(y=y,iset1,design,contrast=2,weights=c(0.5,1,0.5,1))
+> camera(y=y,iset1,design,contrast=2,weights=c(0.5,1,0.5,1),allow.neg.cor=TRUE,inter.gene.cor=NA)
NGenes Correlation Direction PValue
set1 5 -0.2481655 Up 0.001050253
-> camera(y=y,list(set1=iset1,set2=iset2),design,contrast=2)
+> camera(y=y,list(set1=iset1,set2=iset2),design,contrast=2,allow.neg.cor=TRUE,inter.gene.cor=NA)
NGenes Correlation Direction PValue FDR
set1 5 -0.2481655 Up 0.0009047749 0.00180955
set2 5 0.1719094 Down 0.9068364378 0.90683644
+> camera(y=y,iset1,design,contrast=2,weights=c(0.5,1,0.5,1))
+ NGenes Direction PValue
+set1 5 Up 1.105329e-10
+> camera(y=y,list(set1=iset1,set2=iset2),design,contrast=2)
+ NGenes Direction PValue FDR
+set1 5 Up 7.334400e-12 1.466880e-11
+set2 5 Down 8.677115e-01 8.677115e-01
>
> ### with EList arg
>
@@ -1211,9 +1222,12 @@ Down 0 0.997498749
Up 1 0.003001501
UpOrDown 1 0.006000000
Mixed 1 0.006000000
-> camera(y=y,iset1,design,contrast=2)
+> camera(y=y,iset1,design,contrast=2,allow.neg.cor=TRUE,inter.gene.cor=NA)
NGenes Correlation Direction PValue
set1 5 -0.2481655 Up 0.0009047749
+> camera(y=y,iset1,design,contrast=2)
+ NGenes Direction PValue
+set1 5 Up 7.3344e-12
>
> ### eBayes with trend
>
@@ -1368,4 +1382,4 @@ GO:0035295 tube development BP 3 0 2 1.0000000 0.0122943723
>
> proc.time()
user system elapsed
- 2.76 0.18 2.94
+ 2.73 0.20 2.94
--
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