data(fruitfly)
colnames(fruitfly) ## check if arrays are arranged in the default order
gnames <- rownames(fruitfly)
assay <- rep(c("A", "B", "C"), each = 12)
time.grp <- rep(c(1:12), 3)
size <- rep(3, nrow(fruitfly))
out1 <- mb.long(fruitfly, times=12, reps=size, rep.grp = assay, time.grp = time.grp)
summary(out1)
plotProfile(out1, type="b", gnames=gnames, legloc=c(2,15), pch=c("A","B","C"), xlab="Hour")
## Simulate gene expression data
## Note: this simulation is for demonstration purpose only,
## and does not necessarily reflect the real
## features of longitudinal time course data
## one biological condition, 5 time points, 3 replicates
## 500 genes, 10 genes change over time
SS <- matrix(c( 0.01, -0.0008, -0.003, 0.007, 0.002,
-0.0008, 0.02, 0.002, -0.0004, -0.001,
-0.003, 0.002, 0.03, -0.0054, -0.009,
0.007, -0.0004, -0.00538, 0.02, 0.0008,
0.002, -0.001, -0.009, 0.0008, 0.07), ncol=5)
sim.Sigma <- function()
{
S <- matrix(rep(0,25),ncol=5)
x <- mvrnorm(n=10, mu=rep(0,5), Sigma=10*SS)
for(i in 1:10)
S <- S+crossprod(t(x[i,]))
solve(S)
}
sim.data1 <- function(x, indx=1)
{
mu <- rep(runif(1,8,x[1]),5)
if(indx==1) res <- as.numeric(t(mvrnorm(n=3, mu=mu+rnorm(5,sd=4), Sigma=sim.Sigma())))
if(indx==0) res <- as.numeric(t(mvrnorm(n=3, mu=mu, Sigma=sim.Sigma())))
res
}
M1 <- matrix(rep(14,500*15), ncol=15)
M1[1:10,] <- t(apply(M1[1:10,],1,sim.data1))
M1[11:500,] <- t(apply(M1[11:500,],1,sim.data1, 0))
## Which genes are nonconstant?
MB.1D1 <- mb.long(M1, times=5, reps=rep(3, 500))
MB.1D1$percent # check the percent of moderation
plotProfile(MB.1D1,type="b") # plots the no. 1 gene
plotProfile(MB.1D1,type="b",ranking=10) # plots the no. 10 gene
genenames <- as.character(1:500)
plotProfile(MB.1D1, type="b", gid="8", gnames=genenames) #plots the gene with ID "8"
##
MB.1D1.r <- mb.long(M1, type="r", times=5, reps=rep(3, 500))
plotProfile(MB.1D1.r,type="b",gnames=genenames)
plotProfile(MB.1D1.r,type="b", gid="1", gnames=genenames) #plots the gene with ID "1"
## assign the following labellings to columns of M1
## which is actually the same as the default
## Not Run
trt <- rep("wildtype", 15)
assay <- rep(c("A","B","C"), rep(5,3))
time.grp <- rep(c(0, 1, 3, 4, 6), 3)
## MB.1D2 should give the same results as MB.1D1
#MB.1D2 <- mb.long(M1, times=5, reps=rep(3, 500), condition.grp = trt, rep.grp = assay,
#time.grp=time.grp)
## suppose now the replicates are in this order instead
assay <- rep(c("A","C","B"), rep(5,3))
## then
MB.1D3 <- mb.long(M1, times=5, reps=rep(3, 500), condition.grp = trt, rep.grp = assay, time.grp=time.grp)
MB.1D3$rep.group #check the replicate and time group
MB.1D3$time.group
## Now let's simulate another dataset with two biological conditions
## 500 genes also, 10 of them have different expected time course profiles
## between these two biological conditions
## 3 replicates, 5 time points for each condition
sim.data2 <- function(x, indx=1)
{
mu <- rep(runif(1,8,x[1]),5)
if(indx==1)
res <- c(as.numeric(t(mvrnorm(n=3, mu=mu+rnorm(5,sd=5), Sigma=sim.Sigma()))),
as.numeric(t(mvrnorm(n=3, mu=mu+rnorm(5,sd=3.2), Sigma=sim.Sigma()))))
if(indx==0) res <- as.numeric(t(mvrnorm(n=6, mu=mu+rnorm(5,sd=3), Sigma=sim.Sigma())))
res
}
M2 <- matrix(rep(14,500*30), ncol=30)
M2[1:10,] <- t(apply(M2[1:10,],1,sim.data2))
M2[11:500,] <- t(apply(M2[11:500,],1,sim.data2, 0))
## assume it is a paired two-sample problem
trt <- rep(c("wt","mt"),each=15)
assay <- rep(rep(c("1.2.04","2.4.04","3.5.04"),each=5),2)
size <- matrix(3, nrow=500, ncol=2)
MB.paired <- mb.long(M2, method="paired", times=5, reps=size, condition.grp=trt, rep.grp=assay)
MB.paired$con.group # check the condition, replicate and time groups
MB.paired$rep.group
MB.paired$time.group
plotProfile(MB.paired, type="b")
genenames <- as.character(1:500)
plotProfile(MB.paired, gid="12", type="b", gnames=genenames) #plots the gene with ID "12"
### assume it is a unpaired two-sample problem
assay <- rep(c("1.2.04","2.4.04","3.5.04","5.21.04","7.17.04","8.4.04"),each=5)
MB.2D <- mb.long(M2, method="2", times=5, reps=size, condition.grp=trt, rep.grp=assay)
MB.2D$con.group # check the condition, replicate and time groups
MB.2D$rep.group
MB.2D$time.group
plotProfile(MB.2D,type="b", gnames=genenames) # plot the no. 1 gene
## Now let's simulate another dataset with two biological conditions
## 500 genes also, 10 of them have different expected time course profiles
## between these two biological conditions
## the first condition has 3 replicates, while the second condition has 4 replicates,
## 5 time points for each condition
sim.data3 <- function(x, indx=1)
{
mu <- rep(runif(1,8,x[1]),5)
if(indx==1)
res <- c(as.numeric(t(mvrnorm(n=3, mu=mu+rnorm(5,sd=5), Sigma=sim.Sigma()))),
as.numeric(t(mvrnorm(n=4, mu=mu+rnorm(5,sd=3.2), Sigma=sim.Sigma()))))
if(indx==0) res <- as.numeric(t(mvrnorm(n=7, mu=mu+rnorm(5,sd=3), Sigma=sim.Sigma())))
res
}
M3 <- matrix(rep(14,500*35), ncol=35)
M3[1:10,] <- t(apply(M3[1:10,],1,sim.data3))
M3[11:500,] <- t(apply(M3[11:500,],1,sim.data3, 0))
assay <- rep(c("1.2.04","2.4.04","3.5.04","5.21.04","7.17.04","9.10.04","12.1.04"),each=5)
trt <- c(rep(c("wildtype","mutant"),each=15),rep("mutant",5))
## Note that "mutant" < "wildtype", the sample sizes are (4, 3)
size <- matrix(c(4,3), nrow=500, ncol=2, byrow=TRUE)
MB.2D.2 <- mb.long(M3, method="2", times=5, reps=size, rep.grp=assay, condition.grp=trt)
MB.2D.2$con.group # check the condition, replicate and time groups
MB.2D.2$rep.group
MB.2D.2$time.group
plotProfile(MB.2D.2, type="b") # plot the no. 1 gene
Run the code above in your browser using DataLab