EBcoexpressVignette.R
In EBcoexpress: EBcoexpress for Differential Co-Expression Analysis

### R code from vignette source 'EBcoexpressVignette.Rnw'

###################################################
### code chunk number 1: foo
###################################################
options(keep.source = TRUE, width = 60)
foo <- packageDescription("EBcoexpress")


###################################################
### code chunk number 2: packageLoad
###################################################
library(EBcoexpress)
data(fiftyGenes) # A 50-by-125 expression matrix
tinyCond <- c(rep(1,100),rep(2,25))
tinyPat <- ebPatterns(c("1,1","1,2"))


###################################################
### code chunk number 3: Dcreation
###################################################
D <- makeMyD(fiftyGenes, tinyCond, useBWMC=TRUE)


###################################################
### code chunk number 4: initialization
###################################################
set.seed(3)
initHP <- initializeHP(D, tinyCond)


###################################################
### code chunk number 5: lookAtInitHP
###################################################
print(initHP)


###################################################
### code chunk number 6: EMruns
###################################################
zout <- ebCoexpressZeroStep(D, tinyCond, tinyPat, initHP)
oout <- ebCoexpressOneStep(D, tinyCond, tinyPat, initHP)


###################################################
### code chunk number 7: EMrunsNotRun (eval = FALSE)
###################################################
## fout <- ebCoexpressFullTCAECM(D, tinyCond, tinyPat, initHP)


###################################################
### code chunk number 8: pulloff
###################################################
result0 <- zout$POSTPROBS
result1 <- oout$POSTPROBS


###################################################
### code chunk number 9: pulloffNR (eval = FALSE)
###################################################
## resultF <- fout$POSTPROBS


###################################################
### code chunk number 10: f1a
###################################################
priorDiagnostic(D, tinyCond, zout, 1)


###################################################
### code chunk number 11: f1b
###################################################
priorDiagnostic(D, tinyCond, zout, 2)


###################################################
### code chunk number 12: f1c
###################################################
priorDiagnostic(D, tinyCond, oout, 1)


###################################################
### code chunk number 13: f1d
###################################################
priorDiagnostic(D, tinyCond, oout, 2)


###################################################
### code chunk number 14: grabbingPairs
###################################################
ppbDC1 <- 1-result1[,1]
crit_s <- crit.fun(result1[,1], 0.05)
kept_s <- ppbDC1[ppbDC1 >= crit_s]
kept_h <- ppbDC1[ppbDC1 >= 0.95]
klabs_s <- names(kept_s)
  # DC pair names, under soft thresholding
klabs_h <- names(kept_h)
  # DC pair names, under hard thresholding


###################################################
### code chunk number 15: TPgeneration
###################################################
ii <- c(rep(1:24, times=24:1), rep(26:49, times=24:1))
base1 <- 2:25
for(j in 3:25) base1 <- c(base1, j:25)
base2 <- 27:50
for(j in 28:50) base2 <- c(base2, j:50)
jj <- c(base1, base2)
TP <- mapply("X",ii,"~","X",jj,FUN=paste,sep="")
names(TP) <- NULL
numDC <- length(TP)


###################################################
### code chunk number 16: moreGrabbing
###################################################
nk_s <- length(kept_s) # No. taken as DC (soft)
nk_h <- length(kept_h) # No. taken as DC (hard)
nY_s <- sum(klabs_s %in% TP)
  # Number of TP taken as DC under soft thresholding
nY_h <- sum(klabs_h %in% TP)
  # Number of TP taken as DC under soft thresholding

(nk_s - nY_s)/nk_s           # Soft threshold Obs. FDR
(nk_h - nY_h)/nk_h           # Hard threshold Obs. FDR
nY_s/numDC                   # Soft threshold Obs. Power
nY_h/numDC                   # Hard threshold Obs. Power


###################################################
### code chunk number 17: rankThem
###################################################
hubs <- rankMyGenes(oout)
print(hubs)


###################################################
### code chunk number 18: newInit
###################################################
newInitHP <- initHP
newInitHP$G <- 3
newInitHP$MUS <- c(-0.5, 0, 0.5)
newInitHP$TAUS <- c(0.1, 0.2, 0.2)
newInitHP$WEIGHTS <- c(0.25, 0.5, 0.25)


###################################################
### code chunk number 19: dummy1 (eval = FALSE)
###################################################
## ebCoexpressMeta(DList, conditionsList, pattern, hpEstsList)


###################################################
### code chunk number 20: dummy2
###################################################
D1 <- D
D2 <- D
DList <- list(D1, D2)
cond1 <- tinyCond
cond2 <- tinyCond
conditionsList <- list(cond1, cond2)
pattern <- ebPatterns(c("1,1","1,2"))
initHP1 <- initHP
initHP2 <- initHP
out1 <- ebCoexpressZeroStep(D1, cond1, pattern, initHP1)
out2 <- ebCoexpressZeroStep(D2, cond2, pattern, initHP2)
hpEstsList <- list(out1$MODEL$HPS, out2$MODEL$HPS)

metaResults <- ebCoexpressMeta(
    DList, conditionsList, pattern, hpEstsList)


###################################################
### code chunk number 21: grab20
###################################################
twentyGeneNames <- dimnames(fiftyGenes)[[1]][c(1:10,26:35)]


###################################################
### code chunk number 22: f2a
###################################################
showNetwork(twentyGeneNames, D, condFocus = 1, gsep = "~",
  layout = "kamada.kawai", seed = 5, vertex.shape="circle",
  vertex.label.cex=1, vertex.color="white", edge.width=2,
  vertex.frame.color="black", vertex.size=20,
  vertex.label.color="black", vertex.label.family="sans")


###################################################
### code chunk number 23: f2b
###################################################
showNetwork(twentyGeneNames, D, condFocus = 2, gsep = "~",
  layout = "kamada.kawai", seed = 5, vertex.shape="circle",
  vertex.label.cex=1, vertex.color="white", edge.width=2,
  vertex.frame.color="black", vertex.size=20,
  vertex.label.color="black", vertex.label.family="sans")


###################################################
### code chunk number 24: f3a
###################################################
showNetwork(twentyGeneNames, D, condFocus = 1, gsep = "~",
  layout = "kamada.kawai", seed = 5, vertex.shape="circle",
  vertex.label.cex=1, vertex.color="white", edge.width=2,
  vertex.frame.color="black", vertex.size=20,
  vertex.label.color="black", vertex.label.family="sans",
  hidingThreshold=0.3)


###################################################
### code chunk number 25: f3b
###################################################
showNetwork(twentyGeneNames, D, condFocus = 2, gsep = "~",
  layout = "kamada.kawai", seed = 5, vertex.shape="circle",
  vertex.label.cex=1, vertex.color="white", edge.width=2,
  vertex.frame.color="black", vertex.size=20,
  vertex.label.color="black", vertex.label.family="sans",
  hidingThreshold=0.3)


###################################################
### code chunk number 26: fig4
###################################################
showPair("X1~X2", fiftyGenes, tinyCond, pch=20,
  xlim=c(-4,4), ylim=c(-4,4))