## Simulate a dataset with 1000 genes and 20 arrays in a 2-sample design.
## The top 100 genes will be differentially expressed at varying levels
g.alt <- 100
g.null <- 900
n <- 20
data<-matrix(rnorm(n*(g.alt+g.null)),g.alt+g.null,n)
data[1:g.alt,1:(n/2)] <- data[1:g.alt,1:(n/2)] +
seq(2,2/g.alt,length=g.alt)
dimnames(data) <- list(c(paste("Alt",1:g.alt),
paste("Null",1:g.null)),
paste("Array",1:n))
## A treatment vector
trt <- rep(c("Trt","Ctr"),each=n/2)
## 2 alt. categories and 18 null categories of size 50
C.matrix <- kronecker(diag(20),rep(1,50))
dimnames(C.matrix) <- list(dimnames(data)[[1]],
c(paste("TrueCat",1:2),paste("NullCat",1:18)))
dim(C.matrix)
results <- safe(data,trt,C.mat = C.matrix,Pi.mat = 100)
results
## SAFE-plot made for the first category
if (interactive()) {
safeplot(results,"TrueCat 1")
}
Run the code above in your browser using DataLab