-
Notifications
You must be signed in to change notification settings - Fork 8
/
diffExpr.R
63 lines (51 loc) · 2.01 KB
/
diffExpr.R
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
library("Biobase")
library("SummarizedExperiment")
# http://www-huber.embl.de/DESeq2paper/data/bottomly_sumexp.RData
load("bottomly_sumexp.RData")
bottomly <- updateObject(bottomly)
# http://www-huber.embl.de/DESeq2paper/script//bottomly/random_subsets.txt
randomSubsets <- read.table("random_subsets.txt",strings=FALSE)
se <- bottomly[,match(randomSubsets[1,],colnames(bottomly))]
colData(se)$run <- colnames(se)
eset <- ExpressionSet(assay(se),
AnnotatedDataFrame(as.data.frame(colData(se))))
###################################################
if (FALSE) {
# load Harold's count table
counts <- readRDS("counts_table_bottomly.rds")
# Harold's count table as eset:
pdata <- as.data.frame(colData(se)[colnames(counts),])
all(rownames(pdata) == colnames(counts))
eset <- ExpressionSet(counts, AnnotatedDataFrame(pdata))
}
###################################################
pData(eset)$condition <- pData(eset)$strain
levels(pData(eset)$condition) <- c("A","B")
pData(eset)$batch <- factor(pData(eset)$experiment.number)
levels(pData(eset)$batch) <- 1:3
library("DESeq2")
library("edgeR")
library("limma")
source("runScripts.R")
algos <- list("DESeq2"=runDESeq2,"edgeR"=runEdgeR,"edgeRQL"=runEdgeRQL,"limma.voom"=runVoom)
namesAlgos <- names(algos)
names(namesAlgos) <- namesAlgos
resTest <- list()
resHeldout <- list()
nreps <- 30
library(pbapply) # for progress bar
set.seed(1) # not needed AFAIK but for safety (e.g. SAMseq)
res <- pblapply(1:nreps, function(i) {
testSet <- as.character(randomSubsets[i,1:6])
heldOutSet <- as.character(randomSubsets[i,-(1:6)])
eTest <- eset[,testSet]
eHeldout <- eset[,heldOutSet]
resTest <- as.data.frame(lapply(namesAlgos,
function(n) algos[[n]](eTest)))
resHeldout <- as.data.frame(lapply(namesAlgos,
function(n) algos[[n]](eHeldout)))
list(resTest=resTest,resHeldout=resHeldout)
})
resTest <- lapply(res, "[[", "resTest")
resHeldout <- lapply(res, "[[", "resHeldout")
save(resTest,resHeldout,namesAlgos,file="sensFDR.rda")