ns <- 40 ## 40 samples
cla <- rep(c("Trt","Ctr"),each=ns/2)
ngene <- 10 ## 10 genes per group
npath <- 10 ## 10 groups
nreal <- 3 ## alter groups ##
nnull <- npath-nreal ## null groups ##
pname <- c(paste("RealP",1:nreal, sep=""), paste("NullP",1:nnull, sep=""))
## Three main inputs in the function ##
## [1] Simulate (gene) expression matrix (emat) ##
rmv <- function(mn, covm, nr, nc){
sigma <- diag(nr)
sigma[sigma==0] <- covm
x1 <- rmvnorm(nc/2, mean=mn, sigma=sigma)
x0 <- rmvnorm(nc/2, mean=rep(0,nr), sigma=sigma)
mat <- t(rbind(x1,x0))
return(mat)
}
covm <- 0.9 ##covariance
ct <- c(6,8,10) ##mean
library(mvtnorm)
emat <- c()
for (i in 1:nreal) emat <- rbind(emat, rmv(rep(ct[i],ngene),covm=covm, ngene, ns)) # for alt pathways
for (i in 1:(npath-nreal)) emat <- rbind(emat, rmv(mn=rep(0,ngene),covm=covm, nr=ngene, nc=ns))
dimnames(emat) <- list(paste("Gene", 1:(ngene*npath),sep=""), cla)
## [2] class label ##
cla
## [3] indicator matrix (row: genes and col: pathways)
imat <- kronecker(diag(npath),rep(1,ngene))
dimnames(imat) <- list(paste("Gene",1:(ngene*npath), sep=""), pname)
results.pcot2 <- pcot2(emat, cla, imat)
results.pcot2$res.sig
results.pcot2$res.all
Run the code above in your browser using DataLab