1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
|
estimateGLMTrendedDisp <- function(y, ...)
UseMethod("estimateGLMTrendedDisp")
estimateGLMTrendedDisp.DGEList <- function(y, design=NULL, method="auto", ...)
{
if(is.null(y$AveLogCPM)) y$AveLogCPM <- aveLogCPM(y)
# Check method
method <- match.arg(method, c("auto", "bin.spline", "bin.loess", "power", "spline"))
ntags <- nrow(y$counts)
if(method=="auto"){
if(ntags < 200) {
method <- "power"
} else {
method <- "bin.spline"
}
}
y$trend.method <- method
d <- estimateGLMTrendedDisp(y=y$counts, design=design, offset=getOffset(y), AveLogCPM=y$AveLogCPM, method=method, weights=y$weights, ...)
y$trended.dispersion <- d
y
}
estimateGLMTrendedDisp.default <- function(y, design=NULL, offset=NULL, AveLogCPM=NULL, method="auto", weights=NULL, ...)
{
# Check y
y <- as.matrix(y)
ntags <- nrow(y)
if(ntags==0) return(numeric(0))
nlibs <- ncol(y)
# Check design
if(is.null(design)) {
design <- matrix(1,ncol(y),1)
rownames(design) <- colnames(y)
colnames(design) <- "Intercept"
} else {
design <- as.matrix(design)
}
if(ncol(design) >= ncol(y)) {
warning("No residual df: cannot estimate dispersion")
return(NA,ntags)
}
# Check offset
if(is.null(offset)) offset <- log(colSums(y))
# Check AveLogCPM
if(is.null(AveLogCPM)) AveLogCPM <- aveLogCPM(y,offset=offset,weights=weights)
# Check method
method <- match.arg(method,c("auto","bin.spline","bin.loess","power","spline"))
if(method=="auto"){
if(ntags < 200) {
method <- "power"
} else {
method <- "bin.spline"
}
}
# Call lower-level function
trend <- switch(method,
bin.spline=dispBinTrend(y, design, offset=offset, method.trend="spline", AveLogCPM=AveLogCPM, weights=weights, ...),
bin.loess=dispBinTrend(y, design, offset=offset, method.trend="loess", AveLogCPM=AveLogCPM, weights=weights, ...),
power=dispCoxReidPowerTrend(y, design, offset=offset, AveLogCPM=AveLogCPM, ...),
spline=dispCoxReidSplineTrend(y, design, offset=offset, AveLogCPM=AveLogCPM, ...)
)
trend$dispersion
}
|
