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 BRCA1 MRE11 0.40 0.4375000
FALSE 2 BRCA1 RFC2 0.70 0.9285714
FALSE 3 BRCA1 RAD51 0.20 1.0000000
FALSE 4 BRCA1 PALB2 0.30 0.3333333
FALSE 5 BRCA1 RAD51AP1 1.00 0.8500000
FALSE 6 BRCA1 RFC5 0.15 0.8333333
The resulting network can be converted to igraph
object using bnlearn::as.igraph()
.
g <- bnlearn::as.igraph(ret$av)
igraph::evcent(g)$vector
## BRCA1 MRE11 RFC2 RAD51 PALB2 RAD51AP1
## 3.504669e-01 2.695253e-01 1.118264e-01 4.413717e-01 1.337147e-01 7.598065e-01
## RFC5 XRCC3 RFC3 BRIP1 DNA2 BRCA2
## 4.190653e-01 1.796466e-01 3.547601e-01 6.475161e-01 7.431875e-02 3.601397e-01
## CHEK1 TOPBP1 RFC4 RHNO1 EXO1 ATR
## 1.000000e+00 2.045001e-01 8.896507e-02 3.832701e-01 1.439518e-01 1.971832e-01
## RMI2 RMI1 XRCC2 BLM
## 1.564404e-17 3.190784e-01 2.267392e-01 5.471876e-01
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.0 RC (2023-04-13 r84269)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.2 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.17-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.68.0 org.Hs.eg.db_3.17.0
## [3] AnnotationDbi_1.62.0 DESeq2_1.40.0
## [5] SummarizedExperiment_1.30.0 Biobase_2.60.0
## [7] MatrixGenerics_1.12.0 matrixStats_0.63.0
## [9] GenomicRanges_1.52.0 GenomeInfoDb_1.36.0
## [11] IRanges_2.34.0 S4Vectors_0.38.0
## [13] BiocGenerics_0.46.0 bnlearn_4.8.1
## [15] CBNplot_1.0.0 BiocStyle_2.28.0
##
## loaded via a namespace (and not attached):
## [1] splines_4.3.0 later_1.3.0
## [3] bitops_1.0-7 ggplotify_0.1.0
## [5] filelock_1.0.2 tibble_3.2.1
## [7] oaqc_1.0 polyclip_1.10-4
## [9] graph_1.78.0 lifecycle_1.0.3
## [11] lattice_0.21-8 MASS_7.3-59
## [13] ggdist_3.2.1 magrittr_2.0.3
## [15] limma_3.56.0 sass_0.4.5
## [17] rmarkdown_2.21 jquerylib_0.1.4
## [19] yaml_2.3.7 httpuv_1.6.9
## [21] cowplot_1.1.1 DBI_1.1.3
## [23] RColorBrewer_1.1-3 zlibbioc_1.46.0
## [25] purrr_1.0.1 ggraph_2.1.0
## [27] RCurl_1.98-1.12 yulab.utils_0.0.6
## [29] tweenr_2.0.2 rappdirs_0.3.3
## [31] pvclust_2.2-0 GenomeInfoDbData_1.2.10
## [33] enrichplot_1.20.0 ggrepel_0.9.3
## [35] tidytree_0.4.2 reactome.db_1.84.0
## [37] codetools_0.2-19 DelayedArray_0.26.0
## [39] xml2_1.3.3 DOSE_3.26.0
## [41] ggforce_0.4.1 tidyselect_1.2.0
## [43] aplot_0.1.10 farver_2.1.1
## [45] gmp_0.7-1 viridis_0.6.2
## [47] BiocFileCache_2.8.0 jsonlite_1.8.4
## [49] ellipsis_0.3.2 tidygraph_1.2.3
## [51] tools_4.3.0 treeio_1.24.0
## [53] Rcpp_1.0.10 glue_1.6.2
## [55] gridExtra_2.3 xfun_0.39
## [57] qvalue_2.32.0 distributional_0.3.2
## [59] dplyr_1.1.2 withr_2.5.0
## [61] BiocManager_1.30.20 fastmap_1.1.1
## [63] fansi_1.0.4 digest_0.6.31
## [65] R6_2.5.1 mime_0.12
## [67] gridGraphics_0.5-1 colorspace_2.1-0
## [69] GO.db_3.17.0 RSQLite_2.3.1
## [71] utf8_1.2.3 tidyr_1.3.0
## [73] generics_0.1.3 data.table_1.14.8
## [75] graphlayouts_0.8.4 httr_1.4.5
## [77] scatterpie_0.1.9 graphite_1.46.0
## [79] pkgconfig_2.0.3 gtable_0.3.3
## [81] Rmpfr_0.9-2 blob_1.2.4
## [83] XVector_0.40.0 clusterProfiler_4.8.0
## [85] shadowtext_0.1.2 htmltools_0.5.5
## [87] bookdown_0.33 fgsea_1.26.0
## [89] scales_1.2.1 png_0.1-8
## [91] ggfun_0.0.9 knitr_1.42
## [93] tzdb_0.3.0 reshape2_1.4.4
## [95] nlme_3.1-162 curl_5.0.0
## [97] cachem_1.0.7 stringr_1.5.0
## [99] BiocVersion_3.17.1 parallel_4.3.0
## [101] HDO.db_0.99.1 ReactomePA_1.44.0
## [103] pillar_1.9.0 grid_4.3.0
## [105] vctrs_0.6.2 promises_1.2.0.1
## [107] dbplyr_2.3.2 xtable_1.8-4
## [109] evaluate_0.20 magick_2.7.4
## [111] readr_2.1.4 cli_3.6.1
## [113] locfit_1.5-9.7 compiler_4.3.0
## [115] rlang_1.1.0 crayon_1.5.2
## [117] labeling_0.4.2 plyr_1.8.8
## [119] stringi_1.7.12 viridisLite_0.4.1
## [121] BiocParallel_1.34.0 munsell_0.5.0
## [123] Biostrings_2.68.0 lazyeval_0.2.2
## [125] GOSemSim_2.26.0 Matrix_1.5-4
## [127] ExperimentHub_2.8.0 hms_1.1.3
## [129] patchwork_1.1.2 bit64_4.0.5
## [131] ggplot2_3.4.2 KEGGREST_1.40.0
## [133] shiny_1.7.4 highr_0.10
## [135] interactiveDisplayBase_1.38.0 AnnotationHub_3.8.0
## [137] igraph_1.4.2 memoise_2.0.1
## [139] bslib_0.4.2 ggtree_3.8.0
## [141] fastmatch_1.1-3 bit_4.0.5
## [143] downloader_0.4 ape_5.7-1
## [145] gson_0.1.0