Nested Effects
Models-based perturbation inference:
Inference of unobserved perturbations from gene expression
profiles.
Inference of unobserved perturbations from gene expression profiles.
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π) 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 ϕ and the genes with observed gene expression (E-genes) are modelled with the attachment θ with θ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, ϕθ predicts gene expression profiles, which can be compared to the real data. NEMπ iteratively infers a network ϕ based on gene expression profiles and a perturbation profile, and the perturbation profile based on a network ϕ.
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.
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.
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")
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π 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.
Perturbation inference
We use NEMπ 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π result.
## [1] 0.7788576## [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π achieves a higher accuracy.
Note that NEMπ 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 ϕ. However, they usually profit from combinatorial perturbations.
Prior matrix
Alternatively, we can also provide NEMπ with a probabilistic perturbation matrix Γ 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")
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.7788576Final perturbation matrix
The final perturbation matrix Ω over all samples is slightly different from the matrix Γ we used as input. Γ only denotes the most upstream P-gene perturbed. E.g., if A is upstream of B, than Γ 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 Ω = ϕT × Γ.
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")
Session information
## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 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
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] nempi_1.15.0 mnem_1.21.0 BiocStyle_2.35.0
##
## loaded via a namespace (and not attached):
## [1] snowfall_1.84-6.3 RBGL_1.81.0
## [3] deldir_2.0-4 gridExtra_2.3
## [5] permute_0.9-7 rlang_1.1.4
## [7] magrittr_2.0.3 clue_0.3-65
## [9] flexclust_1.4-2 minet_3.63.0
## [11] matrixStats_1.4.1 e1071_1.7-16
## [13] compiler_4.4.1 mgcv_1.9-1
## [15] flexmix_2.3-19 gdata_3.0.1
## [17] png_0.1-8 vctrs_0.6.5
## [19] stringr_1.5.1 pkgconfig_2.0.3
## [21] fastmap_1.2.0 BoutrosLab.plotting.general_7.1.2
## [23] utf8_1.2.4 rmarkdown_2.28
## [25] graph_1.83.0 xfun_0.48
## [27] modeltools_0.2-23 randomForest_4.7-1.2
## [29] cachem_1.1.0 epiNEM_1.29.0
## [31] jsonlite_1.8.9 highr_0.11
## [33] fpc_2.2-13 jpeg_0.1-10
## [35] gmodels_2.19.1 parallel_4.4.1
## [37] prabclus_2.3-4 cluster_2.1.6
## [39] R6_2.5.1 stringi_1.8.4
## [41] bslib_0.8.0 RColorBrewer_1.1-3
## [43] limma_3.61.12 jquerylib_0.1.4
## [45] diptest_0.77-1 Rcpp_1.0.13
## [47] knitr_1.48 zoo_1.8-12
## [49] snow_0.4-4 ggm_2.5.1
## [51] BoolNet_2.1.9 Matrix_1.7-1
## [53] splines_4.4.1 nnet_7.3-19
## [55] igraph_2.1.1 abind_1.4-8
## [57] yaml_2.3.10 vegan_2.6-8
## [59] lattice_0.22-6 tibble_3.2.1
## [61] evaluate_1.0.1 Rtsne_0.17
## [63] tsne_0.1-3.1 proxy_0.4-27
## [65] RcppEigen_0.3.4.0.2 mclust_6.1.1
## [67] kernlab_0.9-33 pillar_1.9.0
## [69] BiocManager_1.30.25 stats4_4.4.1
## [71] ellipse_0.5.0 ggplot2_3.5.1
## [73] fastICA_1.2-5.1 munsell_0.5.1
## [75] scales_1.3.0 RcppArmadillo_14.0.2-1
## [77] gtools_3.9.5 class_7.3-22
## [79] apcluster_1.4.13 glue_1.8.0
## [81] maketools_1.3.1 Linnorm_2.29.0
## [83] tools_4.4.1 interp_1.1-6
## [85] hexbin_1.28.4 robustbase_0.99-4-1
## [87] sys_3.4.3 data.table_1.16.2
## [89] XML_3.99-0.17 buildtools_1.0.0
## [91] fastcluster_1.2.6 grid_4.4.1
## [93] amap_0.8-20 bdsmatrix_1.3-7
## [95] wesanderson_0.3.7 latticeExtra_0.6-30
## [97] naturalsort_0.1.3 colorspace_2.1-1
## [99] sfsmisc_1.1-19 nlme_3.1-166
## [101] pcalg_2.7-12 latex2exp_0.9.6
## [103] cli_3.6.3 fansi_1.0.6
## [105] ggdendro_0.2.0 corpcor_1.6.10
## [107] Rgraphviz_2.49.1 gtable_0.3.6
## [109] DEoptimR_1.1-3 sass_0.4.9
## [111] digest_0.6.37 BiocGenerics_0.53.0
## [113] htmltools_0.5.8.1 lifecycle_1.0.4
## [115] statmod_1.5.0 MASS_7.3-61References:
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.