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
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
|
library(edgeR)
set.seed(0); u <- runif(100)
# generate raw counts from NB, create list object
y <- matrix(rnbinom(80,size=5,mu=10),nrow=20)
y <- rbind(0,c(0,0,2,2),y)
rownames(y) <- paste("Tag",1:nrow(y),sep=".")
d <- DGEList(counts=y,group=rep(1:2,each=2),lib.size=1001:1004)
# estimate common dispersion and find differences in expression
d <- estimateCommonDisp(d)
d$common.dispersion
de <- exactTest(d)
summary(de$table)
topTags(de)
d2 <- estimateTagwiseDisp(d,trend="none",prior.df=20)
summary(d2$tagwise.dispersion)
de <- exactTest(d2,dispersion="common")
topTags(de)
de <- exactTest(d2)
topTags(de)
d2 <- estimateTagwiseDisp(d,trend="movingave",span=0.4,prior.df=20)
summary(d2$tagwise.dispersion)
de <- exactTest(d2)
topTags(de)
summary(exactTest(d2,rejection="smallp")$table$PValue)
summary(exactTest(d2,rejection="deviance")$table$PValue)
d2 <- estimateTagwiseDisp(d,trend="loess",span=0.8,prior.df=20)
summary(d2$tagwise.dispersion)
de <- exactTest(d2)
topTags(de)
d2 <- estimateTagwiseDisp(d,trend="tricube",span=0.8,prior.df=20)
summary(d2$tagwise.dispersion)
de <- exactTest(d2)
topTags(de)
# mglmOneWay
design <- model.matrix(~group,data=d$samples)
mglmOneWay(d[1:10,],design,dispersion=0.2)
mglmOneWay(d[1:10,],design,dispersion=0)
fit <- glmFit(d,design,dispersion=d$common.dispersion,prior.count=0.5/4)
lrt <- glmLRT(fit,coef=2)
topTags(lrt)
fit <- glmFit(d,design,dispersion=d$common.dispersion,prior.count=0.5)
summary(fit$coef)
fit <- glmFit(d,design,prior.count=0.5/4)
lrt <- glmLRT(fit,coef=2)
topTags(lrt)
dglm <- estimateGLMCommonDisp(d,design)
dglm$common.dispersion
dglm <- estimateGLMTagwiseDisp(dglm,design,prior.df=20)
summary(dglm$tagwise.dispersion)
fit <- glmFit(dglm,design,prior.count=0.5/4)
lrt <- glmLRT(fit,coef=2)
topTags(lrt)
dglm <- estimateGLMTrendedDisp(dglm,design)
summary(dglm$trended.dispersion)
dglm <- estimateGLMTrendedDisp(dglm,design,method="power")
summary(dglm$trended.dispersion)
dglm <- estimateGLMTrendedDisp(dglm,design,method="spline")
summary(dglm$trended.dispersion)
dglm <- estimateGLMTrendedDisp(dglm,design,method="bin.spline")
summary(dglm$trended.dispersion)
dglm <- estimateGLMTagwiseDisp(dglm,design,prior.df=20)
summary(dglm$tagwise.dispersion)
# Continuous trend
nlibs <- 3
ntags <- 1000
dispersion.true <- 0.1
# Make first transcript respond to covariate x
x <- 0:2
design <- model.matrix(~x)
beta.true <- cbind(Beta1=2,Beta2=c(2,rep(0,ntags-1)))
mu.true <- 2^(beta.true %*% t(design))
# Generate count data
y <- rnbinom(ntags*nlibs,mu=mu.true,size=1/dispersion.true)
y <- matrix(y,ntags,nlibs)
colnames(y) <- c("x0","x1","x2")
rownames(y) <- paste("Gene",1:ntags,sep="")
d <- DGEList(y)
d <- calcNormFactors(d)
fit <- glmFit(d, design, dispersion=dispersion.true, prior.count=0.5/3)
results <- glmLRT(fit, coef=2)
topTags(results)
d1 <- estimateGLMCommonDisp(d, design, verbose=TRUE)
glmFit(d,design,dispersion=dispersion.true, prior.count=0.5/3)
# Exact tests
y <- matrix(rnbinom(20,mu=10,size=3/2),5,4)
group <- factor(c(1,1,2,2))
ys <- splitIntoGroupsPseudo(y,group,pair=c(1,2))
exactTestDoubleTail(ys$y1,ys$y2,dispersion=2/3)
y <- matrix(rnbinom(5*7,mu=10,size=3/2),5,7)
group <- factor(c(1,1,2,2,3,3,3))
ys <- splitIntoGroupsPseudo(y,group,pair=c(1,3))
exactTestDoubleTail(ys$y1,ys$y2,dispersion=2/3)
exactTestBetaApprox(ys$y1,ys$y2,dispersion=2/3)
y[1,3:4] <- 0
design <- model.matrix(~group)
fit <- glmFit(y,design,dispersion=2/3,prior.count=0.5/7)
summary(fit$coef)
lrt <- glmLRT(fit,contrast=cbind(c(0,1,0),c(0,0,1)))
topTags(lrt)
design <- model.matrix(~0+group)
fit <- glmFit(y,design,dispersion=2/3,prior.count=0.5/7)
lrt <- glmLRT(fit,contrast=cbind(c(-1,1,0),c(0,-1,1),c(-1,0,1)))
topTags(lrt)
# simple Good-Turing algorithm runs.
test1<-1:9
freq1<-c(2018046, 449721, 188933, 105668, 68379, 48190, 35709, 37710, 22280)
goodTuring(rep(test1, freq1))
test2<-c(312, 14491, 16401, 65124, 129797, 323321, 366051, 368599, 405261, 604962)
goodTuring(test2)
|
