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
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
|
# Test suite
```{r functionDefs}
medianPercentShrinkage <- function(res,...) {
idx <- res$baseMean > 10
baseMean <- res$baseMean[idx]
qs <- quantile(baseMean, 0:10/10)
nms <- unname( round( (qs[-1] + qs[-length(qs)])/2 ) )
f <- cut(baseMean, qs)
delta <- 100 * ( (res$lfcMLE - res$log2FoldChange) / res$lfcMLE )[ idx ]
barplot(tapply(delta, f, median, na.rm=TRUE), las=2,
names=nms, xlab="mean expression", ylab="median percent LFC shrinkage",
...)
}
summarizeDESeqRun <- function(x,Wald=TRUE) {
name <- x$name
time <- x$time
res <- x$res
dds <- x$dds
cat(name,"\n")
cat(paste(as.character(design(dds)),collapse=""),"\n")
cat(paste(nrow(dds),"genes",ncol(dds),"samples \n"))
cat(paste(round(unname(time[3])),"seconds \n"))
summary(res)
if (Wald) {
par(mfrow=c(2,2))
yext <- max(abs(res$log2FoldChange),na.rm=TRUE)
plotMA(res,ylim=c(-yext,yext),main=name)
medianPercentShrinkage(res,main=name)
} else {
par(mfrow=c(1,2))
}
plotDispEsts(dds,main=name)
hist(res$pvalue[res$baseMean > 10],col="grey",main="p-values | base-mean > 10",xlab="")
}
library("GenomicRanges")
recount2SE <- function(name) {
filename <- paste0(name,"_eset.RData")
if (!file.exists(filename)) download.file(paste0(
"http://bowtie-bio.sourceforge.net/recount/ExpressionSets/",
filename),filename)
load(filename)
e <- get(paste0(name,".eset"))
se <- SummarizedExperiment(SimpleList(counts=exprs(e)),
colData=DataFrame(pData(e)))
se
}
```
```{r runAirway, cache=TRUE}
library("airway")
data(airway)
dds <- DESeqDataSet(airway, ~ cell + dex)
time <- system.time({ dds <- DESeq(dds) })
res <- results(dds,addMLE=TRUE)
airwayRes <- list(name="airway", time=time, res=res, dds=dds)
rm(time, res, dds)
```
```{r runPasilla, cache=TRUE}
library("pasilla")
library("Biobase")
data("pasillaGenes")
countData <- counts(pasillaGenes)
colData <- pData(pasillaGenes)[,c("condition","type")]
dds <- DESeqDataSetFromMatrix(countData = countData,
colData = colData,
design = ~ type + condition)
dds$condition <- relevel(dds$condition, "untreated","treated")
time <- system.time({ dds <- DESeq(dds) })
res <- results(dds,addMLE=TRUE)
pasillaRes <- list(name="pasilla", time=time, res=res, dds=dds)
rm(time, res, dds)
```
```{r runHammer, cache=TRUE}
se <- recount2SE("hammer")
se$Time[4] <- "2 months"
se$Time <- droplevels(se$Time)
dds <- DESeqDataSet(se, ~ Time + protocol)
dds$protocol <- relevel(dds$protocol, "control")
time <- system.time({ dds <- DESeq(dds) })
res <- results(dds,addMLE=TRUE)
hammerRes <- list(name="hammer", time=time, res=res, dds=dds)
rm(time, res, dds)
```
```{r runBottomly, cache=TRUE}
se <- recount2SE("bottomly")
dds <- DESeqDataSet(se, ~ strain)
time <- system.time({ dds <- DESeq(dds) })
res <- results(dds,addMLE=TRUE)
bottomlyRes <- list(name="bottomly", time=time, res=res, dds=dds)
rm(time, res, dds)
```
```{r runParathyroid, cache=TRUE}
library("DESeq2")
library("parathyroidSE")
data(parathyroidGenesSE)
se <- parathyroidGenesSE
dds0 <- DESeqDataSet(se, ~ patient + treatment)
dds0 <- dds0[,dds0$treatment != "OHT" & dds0$time == "48h"]
dds <- collapseReplicates(dds0, groupby = dds0$sample, run = dds0$run)
dds$treatment <- factor(dds$treatment, levels=c("Control","DPN"))
time <- system.time({ dds <- DESeq(dds) })
res <- results(dds,addMLE=TRUE)
parathyroidRes <- list(name="parathyroid", time=time, res=res, dds=dds)
rm(time, res, dds)
```
```{r runFission, cache=TRUE}
library("fission")
data(fission)
dds <- DESeqDataSet(fission, ~ strain + minute + strain:minute)
time <- system.time({ dds <- DESeq(dds, test="LRT", reduced= ~ strain + minute) })
res <- results(dds)
fissionRes <- list(name="fission", time=time, res=res, dds=dds)
rm(time, res, dds)
```
```{r plotAirway, fig.width=9, fig.height=9}
summarizeDESeqRun(airwayRes)
```
```{r plotPasilla, fig.width=9, fig.height=9}
summarizeDESeqRun(pasillaRes)
```
```{r plotHammer, fig.width=9, fig.height=9}
summarizeDESeqRun(hammerRes)
```
```{r plotBottomly, fig.width=9, fig.height=9}
summarizeDESeqRun(bottomlyRes)
```
```{r plotParathryoid, fig.width=9, fig.height=9}
summarizeDESeqRun(parathyroidRes)
```
```{r plotFission, fig.width=9, fig.height=4}
summarizeDESeqRun(fissionRes, Wald=FALSE)
```
```{r, fig.width=5, fig.height=5}
gene <- rownames(fissionRes$res)[which.min(fissionRes$res$pvalue)]
data <- plotCounts(fissionRes$dds, gene, intgroup=c("minute","strain"),
returnData=TRUE, transform=TRUE)
library("ggplot2")
ggplot(data, aes(minute, count, color=strain, group=strain)) +
ylab("log2 count") + geom_point() + geom_smooth(se=FALSE,method="loess") +
ggtitle(gene)
```
```{r}
sapply(list(airway=airwayRes, pasilla=pasillaRes, hammer=hammerRes,
bottomly=bottomlyRes, parathyroid=parathyroidRes, fission=fissionRes),
function(z) unname(z$time[3]))
```
```{r}
sessionInfo()
```
|
