Browse code

8 Jan 2016: limma 3.31.8

- decideTests() is now an S3 generic function with a default method
and a method for MArrayLM objects.


git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/limma@125765 bc3139a8-67e5-0310-9ffc-ced21a209358

Gordon Smyth authored on 08/01/2017 05:49:49
Showing 5 changed files

... ...
@@ -1,6 +1,6 @@
1 1
 Package: limma
2
-Version: 3.31.7
3
-Date: 2016-12-14
2
+Version: 3.31.8
3
+Date: 2017-01-08
4 4
 Title: Linear Models for Microarray Data
5 5
 Description: Data analysis, linear models and differential expression for microarray data.
6 6
 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]
... ...
@@ -59,6 +59,8 @@ S3method(cbind,MAList)
59 59
 S3method(cbind,RGList)
60 60
 S3method(cbind,EList)
61 61
 S3method(cbind,EListRaw)
62
+S3method(decideTests,default)
63
+S3method(decideTests,MArrayLM)
62 64
 S3method(detectionPValues,default)
63 65
 S3method(detectionPValues,EListRaw)
64 66
 S3method(dim,MAList)
... ...
@@ -4,13 +4,14 @@ setClass("TestResults",representation("matrix"))
4 4
 
5 5
 summary.TestResults <- function(object,...)
6 6
 #	Gordon Smyth
7
-#	26 Feb 2004.  Last modified 31 Jan 2005.
7
+#	Created 26 Feb 2004.  Last modified 6 Jan 2017.
8 8
 {
9
-#	apply(object,2,table)
10
-	tab <- array(0,c(3,ncol(object)),dimnames=list(c("-1","0","1"),colnames(object)))
11
-	tab[1,] <- colSums(object== -1,na.rm=TRUE)
12
-	tab[2,] <- colSums(object== 0,na.rm=TRUE)
13
-	tab[3,] <- colSums(object== 1,na.rm=TRUE)
9
+	Levels <- attr(object,"levels")
10
+	if(is.null(Levels)) Levels <- c(-1L,0L,1L)
11
+	nlevels <- length(Levels)
12
+	tab <- matrix(0L,nlevels,ncol(object))
13
+	dimnames(tab) <- list(as.character(Levels),colnames(object))
14
+	for (i in 1:nlevels) tab[i,] <- colSums(object==Levels[i],na.rm=TRUE)
14 15
 	class(tab) <- "table"
15 16
 	tab
16 17
 }
... ...
@@ -20,12 +21,84 @@ setMethod("show","TestResults",function(object) {
20 21
 	printHead([email protected])
21 22
 })
22 23
 
23
-decideTests <- function(object,method="separate",adjust.method="BH",p.value=0.05,lfc=0)
24
+decideTests <- function(object,...) UseMethod("decideTests")
25
+
26
+decideTests.default <- function(object,method="separate",adjust.method="BH",p.value=0.05,lfc=0,coefficients=NULL,cor.matrix=NULL,tstat=NULL,df=Inf,genewise.p.value=NULL,...)
27
+#	Accept or reject hypothesis tests across genes and contrasts
28
+#	Gordon Smyth
29
+#	17 Aug 2004. Last modified 8 Jan 2017.
30
+{
31
+	method <- match.arg(method,c("separate","global","hierarchical","nestedF"))
32
+	if(method=="nestedF") stop("nestedF adjust method requires an MArrayLM object",call.=FALSE)
33
+
34
+	adjust.method <- match.arg(adjust.method,c("none","bonferroni","holm","BH","fdr","BY"))
35
+	if(adjust.method=="fdr") adjust.method <- "BH"
36
+
37
+	p <- as.matrix(object)
38
+	if(any(p>1) || any(p<0)) stop("object doesn't appear to be a matrix of p-values")
39
+
40
+	switch(method,
41
+
42
+	  separate={
43
+		for (i in 1:ncol(p)) p[,i] <- p.adjust(p[,i],method=adjust.method)
44
+
45
+	},global={
46
+		p[] <- p.adjust(p[],method=adjust.method)
47
+
48
+	},hierarchical={
49
+		if(is.null(genewise.p.value)) {
50
+#			Apply Simes' method by rows to get genewise p-values
51
+			genewise.p.value <- rep_len(1,nrow(p))
52
+			ngenes <- nrow(p)
53
+			ncontrasts <- ncol(p)
54
+			Simes.multiplier <- ncontrasts/(1:ncontrasts)
55
+			for (g in 1:ngenes) {
56
+				op <- sort(p[g,],na.last=TRUE)
57
+				genewise.p.value[g] <- min(op*Simes.multiplier,na.rm=TRUE)
58
+			}
59
+		}
60
+#		Adjust genewise p-values
61
+		DEgene <- p.adjust(genewise.p.value,method=adjust.method) <= p.value
62
+#		Adjust row-wise p-values
63
+		p[!DEgene,] <- 1
64
+		gDE <- which(DEgene)
65
+		for (g in gDE) p[g,] <- p.adjust(p[g,],method=adjust.method)
66
+#		Adjust p-value cutoff for number of DE genes
67
+		nDE <- length(gDE)
68
+		a <- switch(adjust.method,
69
+			none=1,
70
+			bonferroni=1/ngenes,
71
+			holm=1/(ngenes-nDE+1),
72
+			BH=nDE/ngenes,
73
+			BY=nDE/ngenes/sum(1/(1:ngenes))
74
+		)
75
+		p.value <- a*p.value
76
+	},nestedF={
77
+		stop("nestedF adjust method requires an MArrayLM object",call.=FALSE)
78
+	})
79
+
80
+	isDE <- array(0L,dim(p),dimnames=dimnames(p))
81
+	isDE[p <= p.value] <- 1L
82
+	if(is.null(coefficients)) coefficients <- tstat
83
+	if(is.null(coefficients)) {
84
+		attr(isDE,"levels") <- c(0L,1L)
85
+	} else {
86
+		attr(isDE,"levels") <- c(-1L,0L,1L)
87
+		coefficients <- as.matrix(coefficients)
88
+		if( !all(dim(coefficients)==dim(p)) ) stop("dim(object) disagrees with dim(coefficients)")
89
+		i <- coefficients<0
90
+		isDE[i] <- -isDE[i]
91
+		if(lfc>0) isDE[ abs(coefficients)<lfc ] <- 0L
92
+	}
93
+
94
+	new("TestResults",isDE)
95
+}
96
+
97
+decideTests.MArrayLM <- function(object,method="separate",adjust.method="BH",p.value=0.05,lfc=0,...)
24 98
 #	Accept or reject hypothesis tests across genes and contrasts
25 99
 #	Gordon Smyth
26
-#	17 Aug 2004. Last modified 13 August 2010.
100
+#	17 Aug 2004. Last modified 8 Jan 2017.
27 101
 {
28
-	if(!is(object,"MArrayLM")) stop("Need MArrayLM object")
29 102
 	if(is.null(object$p.value)) object <- eBayes(object)
30 103
 	method <- match.arg(method,c("separate","global","hierarchical","nestedF"))
31 104
 	adjust.method <- match.arg(adjust.method,c("none","bonferroni","holm","BH","fdr","BY"))
... ...
@@ -1,3 +1,8 @@
1
+ 8 Jan 2016: limma 3.31.8
2
+
3
+- decideTests() is now an S3 generic function with a default method
4
+  and a method for MArrayLM objects.
5
+
1 6
 14 Dec 2016: limma 3.31.7
2 7
 
3 8
 - Bug fix for contrastAsCoef() when there is more than one contrast.
... ...
@@ -1,42 +1,55 @@
1 1
 \name{decideTests}
2 2
 \alias{decideTests}
3
+\alias{decideTests.default}
4
+\alias{decideTests.MArrayLM}
3 5
 \title{Multiple Testing Across Genes and Contrasts}
4 6
 \description{
5
-Classify a series of related t-statistics as up, down or not significant.
6
-A number of different multiple testing schemes are offered which adjust for multiple testing down the genes as well as across contrasts for each gene.
7
+Identify which genes are significantly differentially expressed for each contrast from a fit object containing p-values and test statistics.
8
+A number of different multiple testing strategies are offered that adjust for multiple testing down the genes as well as across contrasts for each gene.
7 9
 }
8 10
 \usage{
9
-decideTests(object,method="separate",adjust.method="BH",p.value=0.05,lfc=0)
11
+\method{decideTests}{MArrayLM}(object, method = "separate", adjust.method = "BH", p.value = 0.05, lfc = 0, \dots)
12
+\method{decideTests}{default}(object, method = "separate", adjust.method = "BH", p.value = 0.05, lfc = 0,
13
+           coefficients = NULL, cor.matrix = NULL, tstat = NULL, df = Inf, genewise.p.value = NULL, \dots)
10 14
 }
11 15
 \arguments{
12
-  \item{object}{\code{MArrayLM} object output from \code{eBayes} or \code{treat} from which the t-statistics may be extracted.}
13
-  \item{method}{character string specify how probes and contrasts are to be combined in the multiple testing strategy.  Choices are \code{"separate"}, \code{"global"}, \code{"hierarchical"}, \code{"nestedF"} or any partial string.}
16
+  \item{object}{a numeric matrix of p-values or an \code{MArrayLM} object from which p-values and t-statistics can be extracted.}
17
+  \item{method}{character string specifying how genes and contrasts are to be combined in the multiple testing scheme.  Choices are \code{"separate"}, \code{"global"}, \code{"hierarchical"} or \code{"nestedF"}.}
14 18
   \item{adjust.method}{character string specifying p-value adjustment method.  Possible values are \code{"none"}, \code{"BH"}, \code{"fdr"} (equivalent to \code{"BH"}), \code{"BY"} and \code{"holm"}. See \code{\link[stats]{p.adjust}} for details.}
15
-  \item{p.value}{numeric value between 0 and 1 giving the desired size of the test}
16
-  \item{lfc}{minimum log2-fold-change required}
19
+  \item{p.value}{numeric value between 0 and 1 giving the required family-wise error rate or false discovery rate.}
20
+  \item{lfc}{numeric, minimum absolute log2-fold-change required.}
21
+  \item{coefficients}{numeric matrix of coefficients or log2-fold-changes. Of same dimensions as \code{object}.}
22
+  \item{cor.matrix}{correlation matrix of coefficients. Square matrix of dimension \code{ncol(object)}.}
23
+  \item{tstat}{numeric matrix of t-statistics. Of same dimensions as \code{object}.}
24
+  \item{df}{numeric vector of length \code{nrow(object)} giving degrees of freedom for the t-statistics.}
25
+  \item{genewise.p.value}{numeric vector of length \code{nrow(object)} containing summary gene-level p-values for use with \code{method="hierarchical"}.}
26
+  \item{\dots}{other arguments are not used.}
17 27
 }
18 28
 \value{
19 29
 An object of class \code{\link[=TestResults-class]{TestResults}}.
20
-This is essentially a numeric matrix with elements \code{-1}, \code{0} or \code{1} depending on whether each t-statistic is classified as significantly negative, not significant or significantly positive respectively.
30
+This is essentially a numeric matrix with elements \code{-1}, \code{0} or \code{1} depending on whether each t-statistic is classified as significantly negative, not significant or significantly positive.
21 31
 
22 32
 If \code{lfc>0} then contrasts are judged significant only when the log2-fold change is at least this large in absolute value.
23 33
 For example, one might choose \code{lfc=log2(1.5)} to restrict to 50\% changes or \code{lfc=1} for 2-fold changes.
24 34
 In this case, contrasts must satisfy both the p-value and the fold-change cutoff to be judged significant.
25 35
 }
26 36
 \details{
27
-These functions implement multiple testing procedures for determining whether each statistic in a matrix of t-statistics should be considered significantly different from zero.
28
-Rows of \code{tstat} correspond to genes and columns to coefficients or contrasts.
37
+This function can be applied to a matrix of p-values but is more often applied to an \code{MArrayLM} fit object produced by \code{eBayes} or \code{treat}.
38
+In either case, rows of \code{object} correspond to genes and columns to coefficients or contrasts.
29 39
 
30
-The setting \code{method="separate"} is equivalent to using \code{topTable} separately for each coefficient in the linear model fit, and will give the same lists of probes if \code{adjust.method} is the same.
40
+This function applies a multiple testing procedure and a significance level cutoff to the statistics contained in \code{object}.
41
+It implements a number of multiple testing procedures for determining whether each statistic should be considered significantly different from zero.
42
+
43
+The setting \code{method="separate"} is equivalent to using \code{topTable} separately for each coefficient in the linear model fit, and will identify the same probes as significantly differentially expressed if \code{adjust.method} is the same.
31 44
 \code{method="global"} will treat the entire matrix of t-statistics as a single vector of unrelated tests.
32 45
 \code{method="hierarchical"} adjusts down genes and then across contrasts.
33 46
 \code{method="nestedF"} adjusts down genes and then uses \code{classifyTestsF} to classify contrasts as significant or not for the selected genes.
34 47
 Please see the limma User's Guide for a discussion of the statistical properties of these methods.
35 48
 }
36 49
 \note{
37
-Although this function enables users to set p-value and lfc cutoffs simultaneously, this is not generally recommended.
38
-If the fold changes and p-values are not highly correlated, then the use of a fold change cutoff can increase the false discovery rate above the nominal level.
39
-Users wanting to use fold change thresholding are recommended to use \code{treat} instead of \code{eBayes}, and to leave \code{lfc} at the default value when using \code{decideTests}.
50
+Although this function enables users to set p-value and lfc cutoffs simultaneously, this combination criterion not usually recommended.
51
+Unless the fold changes and p-values are very highly correlated, the addition of a fold change cutoff can increase the family-wise error rate or false discovery rate above the nominal level.
52
+Users wanting to use fold change thresholding are recommended to use \code{treat} instead of \code{eBayes} and to leave \code{lfc} at the default value when using \code{decideTests}.
40 53
 }
41 54
 \seealso{
42 55
 An overview of multiple testing functions is given in \link{08.Tests}.