The R package to infer and plot Bayesian networks. The network are inferred from expression data based on clusterProfiler or ReactomePA results. It makes use of libraries including clusterProfiler, ReactomePA, bnlearn, graphite and depmap. In this vignette, the description of functions and several use cases are depicted using GSE133624, which contains RNA-Seq data of bladder cancer. The more detail can be found on the book (https://noriakis.github.io/CBNplot/).
BiocManager::install("CBNplot")
library(CBNplot)
library(bnlearn)
library(DESeq2)
library(org.Hs.eg.db)
library(GEOquery)
## Load dataset and make metadata
filePaths <- getGEOSuppFiles("GSE133624")
counts = read.table(rownames(filePaths)[1], header=1, row.names=1)
meta = sapply(colnames(counts), function (x) substring(x,1,1))
meta = data.frame(meta)
colnames(meta) = c("Condition")
dds <- DESeqDataSetFromMatrix(countData = counts,
colData = meta,
design= ~ Condition)
## Prefiltering
filt <- rowSums(counts(dds) < 10) > dim(meta)[1]*0.9
dds <- dds[!filt,]
## Perform DESeq2()
dds = DESeq(dds)
res = results(dds, pAdjustMethod = "bonferroni")
## apply variance stabilizing transformation
v = vst(dds, blind=FALSE)
vsted = assay(v)
## Define the input genes, and use clusterProfiler::bitr to convert the ID.
sig = subset(res, padj<0.05)
cand.entrez = clusterProfiler::bitr(rownames(sig),
fromType="ENSEMBL", toType="ENTREZID", OrgDb=org.Hs.eg.db)$ENTREZID
## Perform enrichment analysis
pway = ReactomePA::enrichPathway(gene = cand.entrez)
pway = clusterProfiler::setReadable(pway, org.Hs.eg.db)
## Define including samples
incSample = rownames(subset(meta, Condition=="T"))
Then use CBNplot. Basically, you need to supply the enrichment analysis result, normalized expression value and samples to be included. For bngeneplot
, the pathway number in the result
slot of enrichment analysis results must be given.
bngeneplot(results = pway,exp = vsted,
expSample = incSample, pathNum = 15)
Data frame of raw values used in the inference, data frame containing strength and direction, averaged network, and plot can be obtained by specifying returnNet=TRUE
ret <- bngeneplot(results = pway,exp = vsted,
expSample = incSample, pathNum = 15, returnNet=TRUE)
ret$str |> head()
FALSE from to strength direction
FALSE 1 CENPQ SPDL1 0.40 0.7500000
FALSE 2 CENPQ NUP133 0.30 0.5000000
FALSE 3 CENPQ NUP37 0.45 0.6666667
FALSE 4 CENPQ NDC80 0.25 0.4000000
FALSE 5 CENPQ XPO1 0.40 0.5000000
FALSE 6 CENPQ BIRC5 0.30 0.5000000
The resulting network can be converted to igraph
object using bnlearn::as.igraph()
.
g <- bnlearn::as.igraph(ret$av)
igraph::evcent(g)$vector
## CENPQ SPDL1 NUP133 NUP37 NDC80 XPO1 BIRC5 PSMB5
## 0.4043407 0.2454753 0.5489709 0.2138337 0.3534869 0.5340994 0.3925472 0.1279739
## PSMA6 PSMA7 CENPI PSMA2 UBE2S NUP107 CENPA PSMD14
## 0.4173347 0.3032090 0.2188213 0.3115492 0.3782319 0.2363872 0.7303752 0.3345827
## CDC20 CENPF CENPL KIF18A PDS5A ZWINT CENPK NUP85
## 0.5345431 0.4975080 0.1678548 0.6729121 0.5689742 0.3450721 0.3997586 0.4273428
## SGO1 CDCA8 ESPL1 KNL1 CENPO CENPE KIF2C NUF2
## 0.4666727 0.8794683 0.6096135 0.2543809 0.6283373 0.4181057 0.6744517 0.3627395
## CDCA5 INCENP DSN1 CENPU SPC25 CENPH ANAPC1 AHCTF1
## 0.8284538 0.6608254 0.3850940 0.9349226 0.5191835 0.3064180 0.1884831 0.1738264
## BUB3 SKA1 BUB1B PSMD4 SPC24 SGO2 MAD2L1 PTTG1
## 0.3595813 0.8919014 0.5674659 0.2278502 0.5604701 0.2294218 1.0000000 0.5845318
## RAD21 PLK1 TUBA1C BUB1 ZWILCH UBE2C CKAP5 AURKB
## 0.2761026 0.1164351 0.3741177 0.7620561 0.3926130 0.6734771 0.4457428 0.5683693
## RCC2 KNTC1 ERCC6L
## 0.4743328 0.5470337 0.5785199
The relationship between pathways can be drawn by bnpathplot
. The number to be included in the inference can be specified by nCategory
.
bnpathplot(results = pway,exp = vsted,
expSample = incSample, nCategory=10, shadowText = TRUE)
bngeneplotCustom
and bnpathplotCustom
can be used to customize visualization with more flexibility, like highlighting the nodes and edges of interest by glowEdgeNum
and hub
.
bnpathplotCustom(results = pway, exp = vsted, expSample = incSample,
fontFamily="serif", glowEdgeNum=3, hub=3)
bngeneplotCustom(results = pway, exp = vsted, expSample = incSample,
pathNum=15, fontFamily="sans", glowEdgeNum=NULL, hub=3)
The detailed usage for the package, like including covariates to the plot and probabilistic reasoning is available in the package documentation (https://noriakis.github.io/CBNplot/).
sessionInfo()
## R version 4.3.2 Patched (2023-11-13 r85521)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.18-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
##
## 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
##
## time zone: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] GEOquery_2.70.0 org.Hs.eg.db_3.18.0
## [3] AnnotationDbi_1.64.1 DESeq2_1.42.0
## [5] SummarizedExperiment_1.32.0 Biobase_2.62.0
## [7] MatrixGenerics_1.14.0 matrixStats_1.2.0
## [9] GenomicRanges_1.54.1 GenomeInfoDb_1.38.3
## [11] IRanges_2.36.0 S4Vectors_0.40.2
## [13] BiocGenerics_0.48.1 bnlearn_4.9.1
## [15] CBNplot_1.2.1 BiocStyle_2.30.0
##
## loaded via a namespace (and not attached):
## [1] splines_4.3.2 later_1.3.2
## [3] bitops_1.0-7 ggplotify_0.1.2
## [5] filelock_1.0.3 tibble_3.2.1
## [7] oaqc_1.0 polyclip_1.10-6
## [9] graph_1.80.0 lifecycle_1.0.4
## [11] lattice_0.22-5 MASS_7.3-60
## [13] ggdist_3.3.1 magrittr_2.0.3
## [15] limma_3.58.1 sass_0.4.8
## [17] rmarkdown_2.25 jquerylib_0.1.4
## [19] yaml_2.3.8 httpuv_1.6.13
## [21] cowplot_1.1.2 DBI_1.2.0
## [23] RColorBrewer_1.1-3 abind_1.4-5
## [25] zlibbioc_1.48.0 purrr_1.0.2
## [27] ggraph_2.1.0 RCurl_1.98-1.13
## [29] yulab.utils_0.1.1 tweenr_2.0.2
## [31] rappdirs_0.3.3 pvclust_2.2-0
## [33] GenomeInfoDbData_1.2.11 enrichplot_1.22.0
## [35] ggrepel_0.9.4 tidytree_0.4.6
## [37] reactome.db_1.86.2 codetools_0.2-19
## [39] DelayedArray_0.28.0 xml2_1.3.6
## [41] DOSE_3.28.2 ggforce_0.4.1
## [43] tidyselect_1.2.0 aplot_0.2.2
## [45] farver_2.1.1 gmp_0.7-3
## [47] viridis_0.6.4 BiocFileCache_2.10.1
## [49] jsonlite_1.8.8 ellipsis_0.3.2
## [51] tidygraph_1.3.0 tools_4.3.2
## [53] treeio_1.26.0 Rcpp_1.0.11
## [55] glue_1.6.2 gridExtra_2.3
## [57] SparseArray_1.2.2 xfun_0.41
## [59] qvalue_2.34.0 distributional_0.3.2
## [61] dplyr_1.1.4 withr_2.5.2
## [63] BiocManager_1.30.22 fastmap_1.1.1
## [65] fansi_1.0.6 digest_0.6.33
## [67] R6_2.5.1 mime_0.12
## [69] gridGraphics_0.5-1 colorspace_2.1-0
## [71] GO.db_3.18.0 RSQLite_2.3.4
## [73] utf8_1.2.4 tidyr_1.3.0
## [75] generics_0.1.3 data.table_1.14.10
## [77] graphlayouts_1.0.2 httr_1.4.7
## [79] S4Arrays_1.2.0 scatterpie_0.2.1
## [81] graphite_1.48.0 pkgconfig_2.0.3
## [83] gtable_0.3.4 Rmpfr_0.9-4
## [85] blob_1.2.4 XVector_0.42.0
## [87] clusterProfiler_4.10.0 shadowtext_0.1.2
## [89] htmltools_0.5.7 bookdown_0.37
## [91] fgsea_1.28.0 scales_1.3.0
## [93] png_0.1-8 ggfun_0.1.3
## [95] knitr_1.45 tzdb_0.4.0
## [97] reshape2_1.4.4 nlme_3.1-164
## [99] curl_5.2.0 cachem_1.0.8
## [101] stringr_1.5.1 BiocVersion_3.18.1
## [103] parallel_4.3.2 HDO.db_0.99.1
## [105] ReactomePA_1.46.0 pillar_1.9.0
## [107] grid_4.3.2 vctrs_0.6.5
## [109] promises_1.2.1 dbplyr_2.4.0
## [111] xtable_1.8-4 evaluate_0.23
## [113] magick_2.8.2 readr_2.1.4
## [115] cli_3.6.2 locfit_1.5-9.8
## [117] compiler_4.3.2 rlang_1.1.2
## [119] crayon_1.5.2 labeling_0.4.3
## [121] plyr_1.8.9 fs_1.6.3
## [123] stringi_1.8.3 viridisLite_0.4.2
## [125] BiocParallel_1.36.0 munsell_0.5.0
## [127] Biostrings_2.70.1 lazyeval_0.2.2
## [129] GOSemSim_2.28.0 Matrix_1.6-4
## [131] ExperimentHub_2.10.0 hms_1.1.3
## [133] patchwork_1.1.3 bit64_4.0.5
## [135] ggplot2_3.4.4 statmod_1.5.0
## [137] KEGGREST_1.42.0 shiny_1.8.0
## [139] highr_0.10 interactiveDisplayBase_1.40.0
## [141] AnnotationHub_3.10.0 igraph_1.6.0
## [143] memoise_2.0.1 bslib_0.6.1
## [145] ggtree_3.10.0 fastmatch_1.1-4
## [147] bit_4.0.5 ape_5.7-1
## [149] gson_0.1.0