README.md
In cbg-ethz/pcNEM: Probabilistic combinatorial nested effects model

pcNEM

Probabilistic combinatorial nested effects model (pc-NEM), is a variant of NEM, which probabilistically models the perturbations for network reconstruction from combinatorial gene knockdown data. The model and the inference algorithm are implemented as part of the R/Bioconductor package nem.

install.packages("devtools") 
library(devtools) 
install_github("cbg-ethz/pcNEM")

A small toy example with 6 S-genes, 6 knockdown experiments, and 90 E-genes. First, we sample a random network and then generate a perturbation map with on-target probabilities set to 0.8 and off-target probabilities set to 0.2. Next, we simulate data using the network and the perturbation map with FPR = 0.05 and FNR = 0.01. Once we have simulated the data, we use pc-NEM to learn the network and noise parameters.

library(pcnem)
library(Rgraphviz)

set.seed(42)

# Sample a toy network with N S-genes and M E-genes
N      <- 6
M      <- 90
Sgenes <- paste0("S", 1:N)
Phi    <- pcnem:::sampleRndNetwork(Sgenes = Sgenes, trans.close = TRUE)

# Generate knockout map with one experiment for each Sgene
KOmap           <- matrix(0, nrow = N, ncol = N )
rownames(KOmap) <- paste0("K", 1:N)
colnames(KOmap) <- Sgenes
# On target probabilities set to 0.8 in this example
diag(KOmap)     <- 0.8
# Off-target probabilities set to 0.2 in this example
KOmap[sample(1:N^2, 4)] <- 0.2

# Generating data
alpha <- 0.05
beta  <- 0.01
D     <- pcnem:::sampleData.pcnem(Phi = Phi, M = M, map = KOmap, typeI.err = alpha, typeII.err = beta)$D

# Setting all the control parameters and runnin pc-nem
control            <- set.default.parameters(unique(colnames(D)), type = "mLL", pcombi = TRUE, trans.close = FALSE)
control$map        <- as.matrix(KOmap)
control$iterations <- 10000
pcnem_mle          <- nem(D, inference = "AdaSimAnneal", control = control, verbose = FALSE)

# True network likelihood
control$para <- c(alpha,beta)
true_mle     <- nem(D, inference = "search", control = control, verbose = FALSE, models = list(Phi))

# Plotting true and inferred network
par(mfrow = c(1,2))
plot(as(Phi, "graphNEL"), main = "True network")
plot(pcnem_mle$graph, main = "MLE graph")

# Estimated noise
cat("Estimated type I error:", pcnem_mle$typeIEst, "\n")
cat("Estimated type II error:", pcnem_mle$typeIIEst)

The FNR is slightly underestimated due to finite sample size from small number of experiments and effects in this example.

You can tune several hyperparameters for your case using the set.default.parameters() function.

pcombi : Logical parameter set to TRUE for pc-NEM. Default set to FALSE. temper : Binary parameter to choose between two variant implementations of adaptive simulated annealing. FALSE corresponds to adaptation of temperatures at varying intervals but cooled at a fixed rate. TRUE corresponds to the scheme described in the paper. Both are very similar in performance. Default is set to FALSE. iterations : Number of iterations for adaptive simulated annealing (ASA). Deafult is 2e4. stepsave : The length of intervals to adapt tempertature and noise. Default is 1e2. revallowed : Binary parameter for including reversal moves. Default allows reversals and is set to 1. AcceptRate : The ideal acceptance rate for ASA. Temp : Initial temperature. Default is 50. AdaptRate : Rate of adaptation of temperature. Default is 0.3. noiseEst : Binary parameter to include estimation of noise parameters. Default is TRUE. moveprobs : Probability of moving between DAG space and noise space if noiseEst = TRUE. Default = c(0.6,0.4). moveprobsNoise : Probability of moving between alpa and beta space for noise estimation. Default = c(0.5,0.5). sigma : Initial covariance matrix for noise estimation.