########################################################################################################################
### R-code used to compute the step-down P-values for the empirical Bayes adjusted t-tests for the Colon cancer data ###
###                                                                                                                  ###
### data.file is an asci file containing 10 columns                                                                  ###
### col1 = gene (probe) ORF                                                                                          ###
### col2 = individual 1, Normal tissue expression level                                                              ###
### col3 = individual 2, Normal tissue expression level                                                              ###
### col4 = individual 3, Normal tissue expression level                                                              ###
### col5 = individual 1, Adenoma tissue expression level                                                             ###
### col6 = individual 2, Adenoma tissue expression level                                                             ###
### col7 = individual 3, Adenoma tissue expression level                                                             ###
### col8 = individual 1, Carcinoma tissue expression level                                                           ###
### col9 = individual 2, Carcinoma tissue expression level                                                           ###
### col0 = individual 3, Carcinoma tissue expression level                                                           ###
########################################################################################################################




dat <- scan("data.file",list(gene=' ',cnnor=0,d5nor=0,ernor=0,cnade=0,d5ade=0,erade=0,cncar=0,d5car=0,ercar=0))
intensity<- c(dat$cnnor,dat$cnade,dat$cncar,dat$d5nor,dat$d5ade,dat$d5car,dat$ernor,dat$erade,dat$ercar)
intensity[intensity<20] <- 20
intensity <- log(intensity)
ngene <- length(dat$gene)
dim(intensity) <- c(ngene,3,3)
temp1 <- sweep(intensity,c(1,2),apply(intensity,c(1,2),mean))
stat1 <- apply(intensity,c(1,3),mean)
stat2 <- apply(intensity,1,mean)
stat3 <- apply(intensity,2,mean)
temp2 <- sweep(temp1,c(1,3),stat1)
res <- sweep(temp2,1,stat2,"+")
sigsq <- sum(res^2)/(2*2*ngene)
se <- sqrt(2*sigsq*(ngene-1)/(3*ngene))
ijmean <- apply(intensity,c(1,2),mean)
gmean <- mean(as.vector(intensity))
ijmean.row <- sweep(ijmean,1,apply(ijmean,1,mean))
ijmean.col <- sweep(ijmean,2,apply(ijmean,2,mean))
inter <- ijmean.row+ijmean.col-ijmean
inter <- inter+gmean
del12 <- inter[,2]-inter[,1]
del23 <- inter[,3]-inter[,2]
del12 <-  (del12)/se
del23 <-  (del23)/se

sigsq23 <- mean(del23^4)*mean(del23^2)/(2*mean(del23^4)-3*(mean(del23^2))^2)
nu23 <- (4*mean(del23^4)-6*(mean(del23^2))^2)/(mean(del23^4)-3*(mean(del23^2)^2))
del23 <- del23*(1-(nu23+1)/(del23^2+nu23*sigsq23))

sigsq12 <- mean(del12^4)*mean(del12^2)/(2*mean(del12^4)-3*(mean(del12^2))^2)
nu12 <- (4*mean(del12^4)-6*(mean(del12^2))^2)/(mean(del12^4)-3*(mean(del12^2)^2))
del12 <- del12*(1-(nu12+1)/(del12^2+nu12*sigsq12))

del12 <- abs(del12)
del23 <- abs(del23)
rk12 <- order(-del12)
rk23 <- order(-del23)
u12 <- del12[rk12]
u23 <- del23[rk23]


c12 <- c23 <- rep(0,ngene)

B <- 500
for (i in 1:B)
{
resid.n <- sample(res,replace=T)
dim(resid.n) <- c(ngene,3,3) 
intensity.n <- sweep(resid.n, MAR=c(1,3), stat1,FUN="+")
intensity.n <- sweep(intensity.n,2,stat3,"+")
intensity.n <- intensity.n -gmean
temp1.n <- sweep(intensity.n,c(1,2),apply(intensity.n,c(1,2),mean))
stat1.n <- apply(intensity.n,c(1,3),mean)
stat2.n <- apply(intensity.n,1,mean)
temp2.n <- sweep(temp1.n,c(1,3),stat1.n)
res.n <- sweep(temp2.n,1,stat2.n,"+")
sigsq.n <- sum(res.n^2)/(2*2*ngene)
se.n <- sqrt(2*sigsq.n*(ngene-1)/(3*ngene))
ijmean.n <- apply(intensity.n,c(1,2),mean)
gmean.n <- mean(as.vector(intensity.n))
ijmean.row.n <- sweep(ijmean.n,1,apply(ijmean.n,1,mean))
ijmean.col.n <- sweep(ijmean.n,2,apply(ijmean.n,2,mean))
inter.n <- ijmean.row.n+ijmean.col.n-ijmean.n
inter.n <- inter.n+gmean.n
del12.n <- inter.n[,2]-inter.n[,1]
del23.n <- inter.n[,3]-inter.n[,2]
del12.n <-  (del12.n/se.n)
del23.n <-  (del23.n/se.n)
sigsq23.n <- mean(del23.n^4)*mean(del23.n^2)/(2*mean(del23.n^4)-3*(mean(del23.n^2))^2)
nu23.n <- (4*mean(del23.n^4)-6*(mean(del23.n^2))^2)/(mean(del23.n^4)-3*(mean(del23.n^2)^2))
del23.n <- del23.n*(1-(nu23.n+1)/(del23.n^2+nu23.n*sigsq23.n))

sigsq12.n <- mean(del12.n^4)*mean(del12.n^2)/(2*mean(del12.n^4)-3*(mean(del12.n^2))^2)
nu12.n <- (4*mean(del12.n^4)-6*(mean(del12.n^2))^2)/(mean(del12.n^4)-3*(mean(del12.n^2)^2))
del12.n <- del12.n*(1-(nu12.n+1)/(del12.n^2+nu12.n*sigsq12.n))

del23.n <- abs(del23.n)
del12.n <- abs(del12.n)

u12.n <- del12.n[rk12]
u23.n <- del23.n[rk23]
for (ii in (ngene-1):1) {
   u12.n[ii] <- max(u12.n[ii+1],u12.n[ii])
   u23.n[ii] <- max(u23.n[ii+1],u23.n[ii])
                        }
for (ii in 1:ngene) {
if (u12.n[ii] >= u12[ii]) c12[ii] <- c12[ii]+1 
if (u23.n[ii] >= u23[ii]) c23[ii] <- c23[ii]+1 
}

cat(i,"\n")
}


p12.n <- c12/B
p23.n <- c23/B


for (ii in 2:ngene) p12.n[ii] <- max(p12.n[ii],p12.n[ii-1])
for (ii in 2:ngene) p23.n[ii] <- max(p23.n[ii],p23.n[ii-1])

genes12 <- dat$gene[rk12]
genes23 <- dat$gene[rk23]


NormalvsAdenoma <- data.frame(genes=genes12, "P-values" = p12.n)
AdenomavsCarcinoma <- data.frame(genes=genes23, "P-values" = p23.n)