vignette_lpNet.R
In lpNet: Linear Programming Model for Network Inference

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

###################################################
### code chunk number 1: setup1
###################################################
library("lpNet")


###################################################
### code chunk number 2: data1
###################################################
n <- 5 # number of genes
K <- 7 # number of perturbation experiments
obs <- matrix(rnorm(n*K), nrow=n, ncol=K)
delta <- apply(obs, 1, mean)
delta_type <- "perGene"
# perturbation vector (0 if gene inactivated and 1 otherwise)
b <- c(0,1,1,1,1,    # perturbation exp1
       1,0,1,1,1,     # perturbation exp2
       1,1,0,1,1,     # perturbation exp3...
       1,1,1,0,1,
       1,1,1,1,0,
       1,0,0,1,1,
       1,1,1,1,1)


###################################################
### code chunk number 3: lp
###################################################
res1 <- doILP(obs, delta, lambda=1, b, n, K, T_=NULL,
              annot=getEdgeAnnot(n), delta_type)
adja1 <- getAdja(res1, n)


###################################################
### code chunk number 4: lp_pos
###################################################
res2 <- doILP(obs, delta, lambda=1, b, n, K, T_=NULL, 
              annot=getEdgeAnnot(n,allpos=TRUE), delta_type, 
              all.pos=TRUE)
adja2 <- getAdja(res2,n)


###################################################
### code chunk number 5: data_ts
###################################################
T_ <- 4 # number of time points
obs_ts <- array(rnorm(n*K*T_), c(n,K,T_))


###################################################
### code chunk number 6: lp_ts
###################################################
res1 <- doILP(obs_ts, delta, lambda=1, b, n, K, T_, 
              annot=getEdgeAnnot(n), delta_type, 
              flag_time_series=TRUE)
adja1 <- getAdja(res1, n)


###################################################
### code chunk number 7: gaussian
###################################################
active_mu <- 1.1
inactive_mu <- 0.9
active_sd <- inactive_sd <- 0.01
mu_type <- "single"


###################################################
### code chunk number 8: loocv
###################################################
times <- 5 # number of times the removed observation is sampled
annot_node <- seq(1,n) # annotation of the nodes
loocv_res <- loocv(kfold=NULL, times, obs, delta, lambda=1, 
                   b, n, K, T_=NULL, annot=getEdgeAnnot(n), 
                   annot_node, active_mu, active_sd, 
                   inactive_mu, inactive_sd, mu_type,  
                   delta_type)

loocv_res$MSE


###################################################
### code chunk number 9: loocv
###################################################
adja_loocv <- getSampleAdjaMAD(loocv_res$edges_all, n, 
                               annot_node)


###################################################
### code chunk number 10: rangeLambda
###################################################
lambda <- calcRangeLambda(obs, delta, delta_type)
MSE <- Inf
for (lamd in lambda) {
  loocv_res <- loocv(kfold=NULL, times, obs, delta, lambda=lamd, 
                     b, n, K, T_=NULL, annot=getEdgeAnnot(n), 
                     annot_node, active_mu, active_sd, 
                     inactive_mu, inactive_sd, mu_type, 
                     delta_type)
  if (loocv_res$MSE < MSE) {
		MSE <- loocv_res$MSE
		edges_all <- loocv_res$edges_all
		bestLambda <- lamd
  }
}
adja_bestLambda <- getSampleAdjaMAD(edges_all, n, annot_node)


###################################################
### code chunk number 11: foldcv
###################################################
kfold <- 5
MSE <- Inf
for (lamd in lambda) {
  kcv_res <- kfoldCV(kfold, times, obs, delta, lambda=lamd, 
                     b, n, K, T_=NULL, annot=getEdgeAnnot(n),
                     annot_node, active_mu, active_sd, 
                     inactive_mu, inactive_sd, mu_type, 
                     delta_type)
  if (kcv_res$MSE < MSE) {
		MSE <- kcv_res$MSE
		edges_all <- kcv_res$edges_all
		bestLambda <- lamd
  }
}
adja_bestLambda <- getSampleAdjaMAD(edges_all, n, annot_node)


###################################################
### code chunk number 12: prior1
###################################################
res3 <- doILP(obs, delta, lambda=1, b, n, K, T_=NULL, 
              annot=getEdgeAnnot(n), delta_type, 
              sourceNode=1, sinkNode=5)

adja3 <- getAdja(res3,n)


###################################################
### code chunk number 13: prior2
###################################################
prior <- list(c("w+_1_2", 1, ">", 1))
res4 <- doILP(obs,delta, lambda=1, b, n, K, T_=NULL, 
              annot=getEdgeAnnot(n), delta_type, 
              prior=prior)
adja4 <- getAdja(res4, n)


###################################################
### code chunk number 14: prior3
###################################################
prior <- list(c("w+_1_2", 1/0.9, ">", 1))
res5 <- doILP(obs, delta, lambda=1, b, n, K, T_=NULL, 
              annot=getEdgeAnnot(n), delta_type, 
              prior=prior)
              
adja5 <- getAdja(res5, n)


###################################################
### code chunk number 15: prior4
###################################################
library("KEGGgraph")

toyKGML <- system.file("extdata/kgml-ed-toy.xml", package="KEGGgraph")
toyGraph <- parseKGML2Graph(toyKGML, genesOnly=FALSE)
adja <- as(toyGraph,"matrix")
entries <- which(adja!=0, arr.ind=TRUE)

### use apply to set the prior from a given adjacency matrix
myFun <- function(el, sign, confidence, rhs) {
  prior <- c(sprintf("w+_%s_%s", el[[1]][1], el[[1]][2], 
             adja[el[[1]][1],el[[1]][2]]), confidence, sign, rhs) 
}
prior <- lapply(apply(entries,1,list), myFun, ">", 1, 1)

res5 <- doILP(obs, delta, lambda=1, b, n, K, T_=NULL, 
              annot=getEdgeAnnot(n), delta_type, 
              prior=prior)
              
adja5 <- getAdja(res5, n)


###################################################
### code chunk number 16: sahinData
###################################################
data("SahinRNAi2008")
dataStim <- dat.normalized[dat.normalized[ ,17] == 1,-17]
dataUnstim <- dat.normalized[dat.normalized[ ,17] == 0,-17]
# summarize replicates; transpose: rows=genes, cols=experiments
dataSt <- t(summarizeRepl(dataStim, type=mean))
dataUnst <- t(summarizeRepl(dataUnstim, type=mean))


###################################################
### code chunk number 17: sahinData_parameter
###################################################
n <- 16 # number of genes
K <- 16 # number of experiments
annot <- getEdgeAnnot(n)
# perturbation vector;       kd of:
b <- c(0,rep(1,15),          # erbb1
       0,0,rep(1,14),           # erbb1 & 2
       0,rep(1,14),0,	   # erbb1 & 3
       1,0,rep(1,13),0,	 # erbb2 & 3
       rep(1,10),0,1,1,1,1,1,   # IGFR
       rep(1,11),0,1,1,1,1,     # ERalpha
       rep(1,12),0,1,1,1,       # MYC
       rep(1,7),0,rep(1,8),     # AKT1
       rep(1,8),0,rep(1,7),     # MEK1
       rep(1,5),0,rep(1,10),    # CDK2
       rep(1,6),0,rep(1,9),     # CDK4
       rep(1,13),0,1,1,         # CDK6
       1,1,0,rep(1,13),         # p21
       1,1,1,0,rep(1,12),       # p27
       rep(1,4),0,rep(1,11),    # Cyclin D1
       rep(1,14),0,1)           # Cyclin E1


###################################################
### code chunk number 18: sahinData_delta
###################################################
delta <- as.double(dataUnst[ ,1])
delta[11:16] <- mean(dataUnst[ ,1], na.rm=T)


###################################################
### code chunk number 19: sahinData_lp
###################################################
resERBB <- doILP(dataSt[ ,-1], delta, lambda=1.83, 
                 b, n, K, T_=NULL, annot, 
                 delta_type, all.pos=FALSE, 
                 sourceNode=c(1,2,16), sinkNode=10)
adjaERBB <- getAdja(resERBB,n)