[med-svn] [r-bioc-limma] 01/02: Imported Upstream version 3.28.4+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 147341f83eaa19960dab51545d6448e31144f56c
Author: Andreas Tille <tille at debian.org>
Date: Wed Jun 22 17:29:22 2016 +0200
Imported Upstream version 3.28.4+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, 143 insertions(+), 572 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index cb99d7f..8b4864d 100755
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
Package: limma
-Version: 3.28.10
-Date: 2016-06-15
+Version: 3.28.4
+Date: 2016-05-11
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-06-16 01:27:03 UTC; biocbuild
+Packaged: 2016-05-11 01:26:45 UTC; biocbuild
diff --git a/NAMESPACE b/NAMESPACE
index ea342dc..b63b6c6 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -59,8 +59,6 @@ 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 46ef938..09d719e 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 installed (or can't be loaded)")
+ if(!OK) stop("AnnotationDbi package required but not available")
# 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 installed (or can't be loaded)")
+ if(!OK) stop("org.Hs.eg.db package required but not available")
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 installed (or can't be loaded)")
+ if(!OK) stop("org.Mm.eg.db package required but not available")
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 installed (or can't be loaded)")
+ if(!OK) stop("org.Rn.eg.db package required but not available")
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 installed (or can't be loaded)")
+ if(!OK) stop("org.Dm.eg.db package required but not available")
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 installed (or can't be loaded)")
+ if(!OK) stop("AnnotationDbi package required but not available")
# 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 installed (or can't be loaded)")
+ if(!OK) stop("org.Hs.eg.db package required but not available")
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 installed (or can't be loaded)")
+ if(!OK) stop("org.Mm.eg.db package required but not available")
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 installed (or can't be loaded)")
+ if(!OK) stop("org.Rn.eg.db package required but not available")
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 installed (or can't be loaded)")
+ if(!OK) stop("org.Dm.eg.db package required but not available")
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
deleted file mode 100644
index a0f8a8f..0000000
--- a/R/detectionPValues.R
+++ /dev/null
@@ -1,41 +0,0 @@
-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 d6e8802..90430bf 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 installed")
+ if(!requireNamespace("statmod",quietly=TRUE)) stop("statmod package required but is not available")
rho <- rep(NA,ngenes)
nafun <- function(e) NA
for (i in 1:ngenes) {
diff --git a/R/fitFDistRobustly.R b/R/fitFDistRobustly.R
index 0d82ad0..9c2b0c6 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 installed")
+ if(!requireNamespace("statmod",quietly=TRUE)) stop("statmod package required but is not available")
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
deleted file mode 100644
index 0d3ba10..0000000
--- a/R/fitmixture.R
+++ /dev/null
@@ -1,48 +0,0 @@
-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 8c9a3a5..04449b2 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=FALSE,inter.gene.cor=0.01,trend.var=FALSE,sort=TRUE,...)
+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,...)
# Competitive gene set test allowing for correlation between genes
# Gordon Smyth and Di Wu
-# Created 2007. Last modified 4 June 2016.
+# Created 2007. Last modified 22 Dec 2015
{
# 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.na(inter.gene.cor) || is.null(inter.gene.cor))
+ fixed.cor <- !is.null(inter.gene.cor)
# Set df for camera tests
if(fixed.cor) {
@@ -197,9 +197,6 @@ 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 18f5989..f5260cd 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 11 May 2016
+# Created 30 January 2015. Last modified 9 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 installed")
+ if(!OK) stop("statmod package required but not available")
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=="p2") {
+ if(standardize=="posterior.sd") {
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 but apply squeezing to robust variances
- if(standardize=="posterior.sd") {
+# Empirical Bayes moderation method II: estimate hyperparameters from residual variances
+ if(standardize=="p2") {
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 f0fd1e3..5319f41 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 installed (or can't be loaded)")
+ if(!OK) stop("GO.db package required but not available")
# Get access to required annotation functions
suppressPackageStartupMessages(OK <- requireNamespace("AnnotationDbi",quietly=TRUE))
- if(!OK) stop("AnnotationDbi package required but not installed (or can't be loaded)")
+ if(!OK) stop("AnnotationDbi package required but not available")
# 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 installed (or can't be loaded)")
+ if(!OK) stop("org.Hs.eg.db package required but not available")
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 installed (or can't be loaded)")
+ if(!OK) stop("org.Mm.eg.db package required but not available")
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 installed (or can't be loaded)")
+ if(!OK) stop("org.Rn.eg.db package required but not available")
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 installed (or can't be loaded)")
+ if(!OK) stop("org.Dm.eg.db package required but not available")
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 installed (or can't be loaded)")
+ if(!OK) stop("org.Pt.eg.db package required but not available")
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 installed (or can't be loaded)")
+ if(!requireNamespace("BiasedUrn",quietly=TRUE)) stop("BiasedUrn package required but is not available")
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 4691704..3d18152 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 installed (or can't be loaded)")
+ if(!requireNamespace("locfit",quietly=TRUE)) stop("locfit required but is not available")
# 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 dbdf393..5c2f660 100755
--- a/R/norm.R
+++ b/R/norm.R
@@ -333,12 +333,7 @@ plotPrintorder <- function(object,layout,start="topleft",slide=1,method="loess",
normalizeBetweenArrays <- function(object, method=NULL, targets=NULL, cyclic.method="fast", ...) {
# Normalize between arrays
# Gordon Smyth
-# 12 Apri 2003. Last revised 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")
- }
+# 12 Apri 2003. Last revised 21 July 2013.
# Default method
if(is.null(method)) {
diff --git a/R/read.idat.R b/R/read.idat.R
index 4c5f607..8e040b5 100644
--- a/R/read.idat.R
+++ b/R/read.idat.R
@@ -1,88 +1,43 @@
-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.
+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.
{
-# 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)
+ 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)
}
diff --git a/R/sepchannel.R b/R/sepchannel.R
index cd0305c..1991562 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 installed")
+ if(!requireNamespace("statmod",quietly=TRUE)) stop("statmod package required but is not available")
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 66c0037..500b25b 100755
--- a/inst/doc/changelog.txt
+++ b/inst/doc/changelog.txt
@@ -1,42 +1,3 @@
-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.
@@ -61,7 +22,7 @@
4 May 2016: limma 3.28.1
-- unnecessary argument 'design' removed from fitted.MArrayLM.
+- unnecesary 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 c3ef58e..ab7d1b2 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 e55cacd..1eb4283 100644
--- a/man/03reading.Rd
+++ b/man/03reading.Rd
@@ -26,7 +26,6 @@ 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 116dfab..6670298 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 approximation to \code{mroast}, especially useful when heteroscedasticity of genes can be ignored.}
+ Fast version of \code{mroast} 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 ce35342..93def95 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{
-\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)
+\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)
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 \code{NA} or \code{NULL}, then an inter-gene correlation will be estimated for each tested set.}
+ \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{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,30 +46,20 @@ 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}.
-
-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.
+This produces a less conservative test that performs well in many practical situations.
}
-\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 (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.}
+\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}
\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}
@@ -90,7 +80,6 @@ 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}}.
@@ -111,7 +100,7 @@ index2 <- 21:40
camera(y, index1, design)
camera(y, index2, design)
-camera(y, list(set1=index1,set2=index2), design, inter.gene.cor=NA)
+camera(y, list(set1=index1,set2=index2), design)
camera(y, list(set1=index1,set2=index2), design, inter.gene.cor=0.01)
}
diff --git a/man/detectionPValue.Rd b/man/detectionPValue.Rd
deleted file mode 100644
index aca5353..0000000
--- a/man/detectionPValue.Rd
+++ /dev/null
@@ -1,58 +0,0 @@
-\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 ee170ad..6054cc9 100755
--- a/man/dupcor.Rd
+++ b/man/dupcor.Rd
@@ -52,40 +52,16 @@ 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://bioinformatics.oxfordjournals.org/content/21/9/2067}]
-[Preprint with corrections: \url{http://www.statsci.org/smyth/pubs/dupcor.pdf}]
+\url{http://www.statsci.org/smyth/pubs/dupcor.pdf}
}
-
\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)
-
-# 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)
+# Also see lmFit examples
-# 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)
+\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)
+}
}
-
\keyword{multivariate}
diff --git a/man/fitmixture.Rd b/man/fitmixture.Rd
deleted file mode 100644
index 5a7db5b..0000000
--- a/man/fitmixture.Rd
+++ /dev/null
@@ -1,66 +0,0 @@
-\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
deleted file mode 100644
index d64cfaa..0000000
--- a/man/logcosh.Rd
+++ /dev/null
@@ -1,26 +0,0 @@
-\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 f3ce8b6..6cb3c70 100644
--- a/man/propexpr.Rd
+++ b/man/propexpr.Rd
@@ -8,25 +8,21 @@ 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 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}.}
+ \item{status}{character vector giving probe types.}
+ \item{labels}{character vector giving probe type identifiers.}
}
-
\details{
-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.
-
+This function estimates the proportion of expressed in a microarray by utilizing the negative control probes.
Illumina BeadChip arrays contain 750~1600 negative control 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.
+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.
}
-
\value{
Numeric vector giving the proportions of expressed probes in each array.
}
@@ -46,16 +42,7 @@ Description to the control probes in Illumina BeadChips can be found in \code{\l
\examples{
\dontrun{
-# 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")
+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 39f66a7..a2e247d 100644
--- a/man/read.idat.Rd
+++ b/man/read.idat.Rd
@@ -6,93 +6,64 @@
\description{Read Illumina BeadArray data from IDAT and manifest (.bgx) files for gene expression platforms.}
\usage{
-read.idat(idatfiles, bgxfile, dateinfo = FALSE, annotation = "Symbol",
- tolerance = 0L, verbose = TRUE)
+read.idat(idatfiles, bgxfile, dateinfo=FALSE, tolerance=0)
}
\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 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?}
+ \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.}
}
\details{
- 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).
+ 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).
+ 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} 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.
+ 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}.
- If more than \code{tolerance} probes in the manifest cannot be found in an IDAT file then the function will return an error.
+ Probe types are indicated in the \code{Status} column of the \code{genes}
+ component of the \code{EListRaw-class} object.
}
\value{
- An \code{EListRaw} object with the following components:
+ An \code{EListRaw-class} object with the following components:
\item{E}{ numeric matrix of raw intensities.}
- \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}.}
+ \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/}
+\url{http://f1000research.com/articles/2-264/v1}
}
\author{Matt Ritchie}
\seealso{
- \code{\link{read.ilmn}} imports gene expression data output by GenomeStudio.
+ \code{read.ilmn} imports gene expression data output by GenomeStudio.
- \code{\link{neqc}} performs normexp by control background correction, log
+ \code{neqc} performs normexp by control background correction, log
transformation and quantile between-array normalization for
Illumina expression data.
- \code{\link{propexpr}} estimates the proportion of expressed probes in a microarray.
-
- \code{\link{detectionPValues}} computes detection p-values from the negative controls.
+ \code{propexpr} estimates the proportion of expressed probes in a microarray.
}
\examples{
\dontrun{
-idatfiles <- dir(pattern="idat")
-bgxfile <- dir(pattern="bgx")
-x <- read.idat(idatfiles, bgxfile)
-x$other$Detection <- detectionPValues(x)
+idatfiles = dir(pattern="idat")
+bgxfile = dir(pattern="bgx")
+data = read.idat(idatfiles, bgxfile)
propexpr(data)
-y <- neqc(data)
+datanorm = neqc(data)
}
}
diff --git a/man/roast.Rd b/man/roast.Rd
index fd336dc..51e50c9 100644
--- a/man/roast.Rd
+++ b/man/roast.Rd
@@ -13,14 +13,14 @@ Rotation gene set testing for linear models.
}
\usage{
-\method{roast}{default}(y, index = NULL, design = NULL, contrast = ncol(design), geneid = NULL,
+\S3method{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)
-\method{mroast}{default}(y, index = NULL, design = NULL, contrast = ncol(design), geneid = NULL,
+\S3method{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)
-\method{fry}{default}(y, index = NULL, design = NULL, contrast = ncol(design), geneid = NULL,
+\S3method{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 000ad2f..12cb45f 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{
-\method{romer}{default}(y, index, design = NULL, contrast = ncol(design),
+\S3method{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 974fff2..1bcf9d1 100755
--- a/tests/limma-Tests.R
+++ b/tests/limma-Tests.R
@@ -231,12 +231,9 @@ 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)
@@ -244,7 +241,6 @@ 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 f085afd..dbff537 100755
--- a/tests/limma-Tests.Rout.save
+++ b/tests/limma-Tests.Rout.save
@@ -1,5 +1,5 @@
-R version 3.3.0 (2016-05-03) -- "Supposedly Educational"
+R Under development (unstable) (2016-03-14 r70331) -- "Unsuffered Consequences"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)
@@ -1191,27 +1191,16 @@ 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),allow.neg.cor=TRUE,inter.gene.cor=NA)
+> camera(y=y,iset1,design,contrast=2,weights=c(0.5,1,0.5,1))
NGenes Correlation Direction PValue
set1 5 -0.2481655 Up 0.001050253
-> camera(y=y,list(set1=iset1,set2=iset2),design,contrast=2,allow.neg.cor=TRUE,inter.gene.cor=NA)
+> camera(y=y,list(set1=iset1,set2=iset2),design,contrast=2)
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
>
@@ -1222,12 +1211,9 @@ Down 0 0.997498749
Up 1 0.003001501
UpOrDown 1 0.006000000
Mixed 1 0.006000000
-> camera(y=y,iset1,design,contrast=2,allow.neg.cor=TRUE,inter.gene.cor=NA)
+> camera(y=y,iset1,design,contrast=2)
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
>
@@ -1382,4 +1368,4 @@ GO:0035295 tube development BP 3 0 2 1.0000000 0.0122943723
>
> proc.time()
user system elapsed
- 2.73 0.20 2.94
+ 2.76 0.18 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