Nested Effects Models-based perturbation inference:
Inference of unobserved perturbations from gene expression profiles.
Inference of unobserved perturbations from gene expression profiles.
2022-11-01
Contents
1 Introduction
If many genes are perturbed in a population of cells, this can lead to diseases like cancer. The perturbations can happen in different ways, e.g. via mutations, copy number abberations or methylation. However, not all perturbations are observed in all samples.
Nested Effects Model-based perturbation inference (NEM\(\pi\)) uses observed perturbation profiles and gene expression data to infer unobserved perturbations and augment observed ones. The causal network of the perturbed genes (P-genes) is modelled as an adjacency matrix \(\phi\) and the genes with observed gene expression (E-genes) are modelled with the attachment \(\theta\) with \(\theta_{ij}=1\), if E-gene \(j\) is attached to S-gene \(i\). If E-gene \(j\) is attached to P-gene \(i\), \(j\) shows an effect for a perturbation of P-gene \(i\). Hence, \(\phi\theta\) predicts gene expression profiles, which can be compared to the real data. NEM\(\pi\) iteratively infers a network \(\phi\) based on gene expression profiles and a perturbation profile, and the perturbation profile based on a network \(\phi\).
2 Installation and loading
Use devtools to install the latest version from github or use the BiocManahger to install the package from bioconductor.
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install("nempi")Load the package with the library function.
library(nempi)3 Small example
We look at a small example for which we first simulate data and then infer the unobserved perturbations. We draw a random network for the perturbed genes (P-genes). Then we simulate gene expression for each sample given the subset of P-genes that has been perturbed in the sample. E.g., if P-gene A is upstream of P-gene B and P-gene A has been perturbed, all E-genes attached to B also show an effect.
3.1 Data simulation
library(mnem)
seed <- 8675309
Pgenes <- 10
Egenes <- 5
samples <- 100
edgeprob <- 0.5
uninform <- floor((Pgenes * Egenes) * 0.1)
Nems <- mw <- 1
noise <- 1
multi <- c(0.2, 0.1)
set.seed(seed)
simmini <- simData(Sgenes = Pgenes, Egenes = Egenes, Nems = Nems, mw = mw, nCells = samples,
uninform = uninform, multi = multi, badCells = floor(samples * 0.1), edgeprob = edgeprob)
data <- simmini$data
ones <- which(data == 1)
zeros <- which(data == 0)
data[ones] <- rnorm(length(ones), 1, noise)
data[zeros] <- rnorm(length(zeros), -1, noise)
epiNEM::HeatmapOP(data, col = "RdBu", cexRow = 0.75, cexCol = 0.75, bordercol = "transparent",
xrot = 0, dendrogram = "both")
Figure 1: Heatmap of the simulated log odds
Effects are blue and no effects
are red. Rows denote the observed E-genes and columns the samples annoted by
P-genes. Each P-gene
has been perturbed in many cells. The E-genes are annotated as how they are
attached in the ground truth. E.g. E-genes named ‘1’ are attached to S-gene
‘1’ in the ground truth.
The typical data input for NEM\(\pi\) consists of a data matrix with samples as columns and E-genes as rows. The columns are either labeled by their perturbed gene(s) or unlabeled (default: ""). After the data simulation all samples are labeled. We unlabel \(50\%\) of the sample and pretend we do not know, which P-gene has been perturbed.
lost <- sample(1:ncol(data), floor(ncol(data) * 0.5))
colnames(data)[lost] <- ""3.2 Perturbation inference
We use NEM\(\pi\) and other methods to infer the perturbations. We use the area under the precision-recall curve as the two measures of accuracy. We also plot the NEM\(\pi\) result.
res <- nempi(data)
fit <- pifit(res, simmini, data)
print(fit$auc)## [1] 0.7788576ressvm <- classpi(data)
fit <- pifit(ressvm, simmini, data, propagate = FALSE)
print(fit$auc)## [1] 0.7090377resnn <- classpi(data, method = "nnet")
fit <- pifit(resnn, simmini, data, propagate = FALSE)
print(fit$auc)## [1] 0.6273228resrf <- classpi(data, method = "randomForest")
fit <- pifit(resrf, simmini, data, propagate = FALSE)
print(fit$auc)## [1] 0.6987766col <- rgb(seq(0, 1, length.out = 10), seq(1, 0, length.out = 10), seq(1, 0, length.out = 10))
plot(res, heatlist = list(col = "RdBu"), barlist = list(col = col))
Compared to support vector machines (svm), neural nets (nnet) and random forest (rf) class prediction, NEM\(\pi\) achieves a higher accuracy.
Note that NEM\(\pi\) is in general more powerful, if the P-genes are connected in a denser network. The other methods perform equally well for sparse or even disconnected network \(\phi\). However, they usually profit from combinatorial perturbations.
3.3 Prior matrix
Alternatively, we can also provide NEM\(\pi\) with a probabilistic perturbation matrix \(\Gamma\) as a prior. In the rows are the potentially perturbed genes and in the columns are the samples. The entries are between \(0\) and \(1\), with the sum of all entries of a sample summing to \(1\).
Gamma <- matrix(0, Pgenes, ncol(data))
rownames(Gamma) <- seq_len(Pgenes)
colnames(Gamma) <- colnames(data)
for (i in seq_len(Pgenes)) {
Gamma[i, grep(paste0("^", i, "_|_", i, "$|_", i, "_|^", i, "$"), colnames(data))] <- 1
}
Gamma <- apply(Gamma, 2, function(x) return(x/sum(x)))
Gamma[is.na(Gamma)] <- 0
epiNEM::HeatmapOP(Gamma, col = "RdBu", cexRow = 0.75, cexCol = 0.75, bordercol = "transparent",
xrot = 0, dendrogram = "both")
Figure 2: Heatmap of the probabilsitic perturbation matrix
colnames(data) <- sample(seq_len(Pgenes), ncol(data), replace = TRUE)
res <- nempi(data, Gamma = Gamma)
fit <- pifit(res, simmini, data)
print(fit$auc)## [1] 0.77885763.4 Final perturbation matrix
The final perturbation matrix \(\Omega\) over all samples is slightly different from the matrix \(\Gamma\) we used as input. \(\Gamma\) only denotes the most upstream P-gene perturbed. E.g., if A is upstream of B, than \(\Gamma\) only denotes a perturbation of A, even though B is also perturbed in every samples in which A is perturbed. We can compute it by the matrix multiplication \(\Omega = \phi^T \times \Gamma\).
Omega <- t(mnem::transitive.closure(res$res$adj)) %*% res$Gamma
epiNEM::HeatmapOP(Omega, col = "RdBu", cexRow = 0.75, cexCol = 0.75, bordercol = "transparent",
xrot = 0, dendrogram = "both")
4 Session information
sessionInfo()## R version 4.2.1 (2022-06-23)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.5 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.16-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.16-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] nempi_1.6.0 mnem_1.14.0 BiocStyle_2.26.0
##
## loaded via a namespace (and not attached):
## [1] amap_0.8-19 Rtsne_0.16 colorspace_2.0-3
## [4] RcppEigen_0.3.3.9.2 deldir_1.0-6 class_7.3-20
## [7] modeltools_0.2-23 mclust_6.0.0 Linnorm_2.22.0
## [10] corpcor_1.6.10 RcppArmadillo_0.11.4.0.1 ggdendro_0.1.23
## [13] clue_0.3-62 proxy_0.4-27 flexmix_2.3-18
## [16] fansi_1.0.3 splines_4.2.1 cachem_1.0.6
## [19] robustbase_0.95-0 knitr_1.40 wesanderson_0.3.6
## [22] jsonlite_1.8.3 apcluster_1.4.10 cluster_2.1.4
## [25] kernlab_0.9-31 sfsmisc_1.1-13 png_0.1-7
## [28] latex2exp_0.9.5 snowfall_1.84-6.2 graph_1.76.0
## [31] BiocManager_1.30.19 compiler_4.2.1 assertthat_0.2.1
## [34] Matrix_1.5-1 fastmap_1.1.0 limma_3.54.0
## [37] cli_3.4.1 formatR_1.12 htmltools_0.5.3
## [40] tools_4.2.1 BoolNet_2.1.7 igraph_1.3.5
## [43] gtable_0.3.1 glue_1.6.2 epiNEM_1.22.0
## [46] dplyr_1.0.10 gmodels_2.18.1.1 Rcpp_1.0.9
## [49] jquerylib_0.1.4 vctrs_0.5.0 gdata_2.18.0.1
## [52] nlme_3.1-160 fpc_2.2-9 xfun_0.34
## [55] fastcluster_1.2.3 stringr_1.4.1 ggm_2.5
## [58] lifecycle_1.0.3 gtools_3.9.3 XML_3.99-0.12
## [61] statmod_1.4.37 DEoptimR_1.0-11 MASS_7.3-58.1
## [64] zoo_1.8-11 scales_1.2.1 minet_3.56.0
## [67] parallel_4.2.1 RBGL_1.74.0 RColorBrewer_1.1-3
## [70] yaml_2.3.6 ggplot2_3.3.6 sass_0.4.2
## [73] bdsmatrix_1.3-6 naturalsort_0.1.3 fastICA_1.2-3
## [76] latticeExtra_0.6-30 stringi_1.7.8 highr_0.9
## [79] randomForest_4.7-1.1 e1071_1.7-12 permute_0.9-7
## [82] BiocGenerics_0.44.0 pcalg_2.7-7 rlang_1.0.6
## [85] pkgconfig_2.0.3 prabclus_2.3-2 matrixStats_0.62.0
## [88] evaluate_0.17 lattice_0.20-45 tidyselect_1.2.0
## [91] magrittr_2.0.3 bookdown_0.29 R6_2.5.1
## [94] magick_2.7.3 snow_0.4-4 generics_0.1.3
## [97] DBI_1.1.3 pillar_1.8.1 mgcv_1.8-41
## [100] abind_1.4-5 nnet_7.3-18 tibble_3.1.8
## [103] tsne_0.1-3.1 interp_1.1-3 utf8_1.2.2
## [106] ellipse_0.4.3 rmarkdown_2.17 jpeg_0.1-9
## [109] grid_4.2.1 flexclust_1.4-1 data.table_1.14.4
## [112] vegan_2.6-4 Rgraphviz_2.42.0 digest_0.6.30
## [115] diptest_0.76-0 stats4_4.2.1 munsell_0.5.0
## [118] bslib_0.4.05 References:
Markowetz, F., Bloch, J., and Spang, R. (2005). Non-transcriptional pathway features reconstructed from secondary effects of rna interference. Bioinformatics, 21(21), 4026–4032.
Markowetz, F., Kostka, D., Troyanskaya, O. G., and Spang, R. (2007). Nested effects models for high-dimensional phenotyping screens. Bioinformatics, 23(13), i305–i312.
Pirkl, M., Beerenwinkel, N.; Single cell network analysis with a mixture of Nested Effects Models, Bioinformatics, Volume 34, Issue 17, 1 September 2018, Pages i964–i971, https://doi.org/10.1093/bioinformatics/bty602.
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK (2015). “limma powers differential expression analyses for RNA-sequencing and microarray studies.” Nucleic Acids Research, 43(7), e47.