Abstract:

This package has developed a tool for performing a novel pathway enrichment analysis based on Bayesian network (BNrich) to investigate the topology features of the pathways. This algorithm as a biologically intuitive, method, analyzes of most structural data acquired from signaling pathways such as causal relationships between genes using the property of Bayesian networks and also infer finalized networks more conveniently by simplifying networks in the early stages and using Least Absolute Shrinkage Selector Operator (LASSO). impacted pathways are ultimately prioritized the by Fisher’s Exact Test on significant parameters. Here, we provide an instance code that applies BNrich in all of the fields described above.

Introduction

This document offers an introductory overview of how to use the package. The BNrich tool uses Bayesian Network (BN) in a new topology-based pathway analysis (TPA) method. The BN has been demonstrated as a beneficial technique for integrating and modeling biological data into causal relationships (1–4). The proposed method utilizes BN to model variations in downstream components (children) as a consequence of the change in upstream components (parents). For this purpose, The method employs 187 KEGG human non-metabolic pathways (5–7) which their cycles were eliminated manually by a biological intuitive, as BN structures and gene expression data to estimate its parameters (8,9). The cycles of inferred networks were eliminated on the basis of biologically intuitive rules instead of using computing algorithms (10). The inferred networks are simplified in two steps; unifying genes and LASSO. Similarly, the originally continuous gene expression data is used to BN parameters learning, rather than discretized data (8). The algorithm estimates regression coefficients by continuous data based on the parameter learning techniques in the BN (11,12). The final impacted pathways are gained by Fisher’s exact test. This method can represent effective genes and biological relations in impacted pathways based on a significant level.

Quick Start

install.packages("BNrich_0.1.0.tar.gz", type="source", repos=NULL)
library("BNrich")

destfile = tempfile("files", fileext = ".rda")
files <- fetch_data_file()
load(destfile)

Data <- system.file("extdata", "Test_DATA.RData", package = "BNrich", mustWork = TRUE)
load(Data)
head(dataH)

head(dataD)

Unify data, the first step of simplification

Initially, we need to unify gene products based on 187 imported signaling pathways (mapkG list) in two states disease (dataD) and control (dataH). This is the first simplification step, unifying nodes in signaling pathways with genes those exist in gene expression data.

unify_results <- unify_path(dataH, dataD, mapkG, pathway.id)

The unify_path function performs the following processes: • Split datasets into KEGG pathways • Delete all gene expression data are not in pathways • Removes all gene products in pathways are not in dataset platforms • Remove any pathways with the number of edges is less than 5 This function returns a list contain data_h,data_d,mapkG1 and pathway.id1. data_h and data_d are lists contain data frames related to control and disease objects unified for any signaling pathways. The mapkG1 is a list contains unified signaling pathways and pathway.id1 is new pathway ID vector based on remained pathways. In the example dataset, the number of edges in the one pathway becomes less than 5 and are removed:

length(mapkG)

187

mapkG1 <- unify_results$mapkG1
length(mapkG1)

186

As well, the number of edges reduces in the remaining pathways. In first pathway hsa:01521 the number of edges from 230 reduces to 204:

pathway.id[1]

"hsa:01521"

mapkG[[1]]

A graphNEL graph with directed edges 
Number of Nodes = 79 
Number of Edges = 230

pathway.id1 <- unify_results$pathway.id1
pathway.id1[1]

"hsa:01521"

mapkG1[[1]]

A graphNEL graph with directed edges
Number of Nodes = 71 
Number of Edges = 204

BN: construct structures and estimate parameters

BN <- BN_struct(unify_results$mapkG1)

data_h <- unify_results$data_h
data_d <- unify_results$data_d
LASSO_results <- LASSO_BN(BN, data_h, data_d)

nrow(arcs(BN[[1]]))

204

nrow(arcs(LASSO_results$BN_H[[1]]))

116

nrow(arcs(LASSO_results$BN_D[[1]]))

116

BN_H <- LASSO_results$BN_H
BN_D <- LASSO_results$BN_D
esti_results <- esti_par(BN_H, BN_D, data_h, data_d)

esti_results$BNs_h[[1]]$` hsa:1978`

Parameters of node hsa:1978 (Gaussian distribution)
Conditional density: hsa:1978 | hsa:2475
Coefficients:
(Intercept)     hsa:2475  
  2.8841264    0.6958609  
Standard deviation of the residuals: 0.3489612 

esti_results$BNs_d[[1]]$`hsa:1978`

Parameters of node hsa:1978 (Gaussian distribution)
Conditional density: hsa:1978 | hsa:2475
Coefficients:
(Intercept)     hsa:2475  
  0.9046357   1.1870730  
Standard deviation of the residuals: 0.2713789

Testing the equality BNs parameters

BN_h <- esti_results$BNs_h
BN_d <- esti_results$BNs_d
coef_h <- esti_results$coef_h
coef_d <- esti_results$coef_d
var_mat_results<- var_mat (data_h, coef_h, BN_h, data_d, coef_d, BN_d)

(var_mat_results$var_mat_Bh[[1]])[5]

          [,1]      [,2]
[1,] 10.177073      -3.630152
[2,] -3.630152      1.296990

(var_mat_results$var_mat_Bd[[1]])[5]

          [,1]      [,2]
[1,]  3.549338      -1.0392040
[2,] -1.039204      0.3053785

var_mat_Bh <- var_mat_results $var_mat_Bh
var_mat_Bd <- var_mat_results $var_mat_Bd
Ttest_results <- parm_Ttest(data_h, coef_h, BN_h, data_d, coef_d, BN_d, var_mat_Bh, var_mat_Bd, pathway.id1)
head(Ttest_results)

Identification of enriched pathways

In the last step we can determine enriched pathways by own threshold on p-value or fdr. Hence we run the Fisher’s exact test for any final pathways. As stated above, the Ttest_results is a data frame contains T-test results for all parameters in final BNs achieved by parm_Ttest function and fdr.value A numeric threshold to determine significant parameters (default is 0.05).

BNrich_results <- BNrich(Ttest_results, pathway.id1, PathName_final, fdr.value = 0.05)
head(BNrich_results)

Session Info

The following package and versions were used in the production of this vignette.

 R version 3.6.1 (2019-07-05)
 Platform: x86_64-w64-mingw32/x64 (64-bit)
 Running under: Windows 7 x64 (build 7601) Service Pack 1

 Matrix products: default

 locale:
 [1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United Kingdom.1252   
 [3] LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C                           
 [5] LC_TIME=English_United Kingdom.1252    

 attached base packages:
 [1] stats     graphics  grDevices utils     datasets  methods   base     

 other attached packages:
 [1] BNrich_0.1.0

 loaded via a namespace (and not attached):
  [1] Rcpp_1.0.2          codetools_0.2-16    lattice_0.20-38     corpcor_1.6.9      
  [5] foreach_1.4.7       glmnet_2.0-18       digest_0.6.20       grid_3.6.1         
  [9] stats4_3.6.1        evaluate_0.14       graph_1.63.0        Matrix_1.2-17      
 [13] rmarkdown_1.14      bnlearn_4.5         iterators_1.0.12    tools_3.6.1        
 [17] parallel_3.6.1      xfun_0.8            yaml_2.2.0          rsconnect_0.8.15   
 [21] compiler_3.6.1      BiocGenerics_0.31.5 htmltools_0.3.6     knitr_1.24

References:

1. Yu J, Smith VA, Wang PP, Hartemink AJ, Jarvis ED. Advances to Bayesian network inference for generating causal networks from observational biological data. Bioinformatics . 2004 Dec 12;20(18):3594–603.
2. Gendelman R, Xing H, Mirzoeva OK, Sarde P, Curtis C, Feiler HS, et al. Bayesian network inference modeling identifies TRIB1 as a novel regulator of cell-cycle progression and survival in cancer cells. Cancer Res. 2017;77(7):1575–85.
3. Luo Y, El Naqa I, McShan DL, Ray D, Lohse I, Matuszak MM, et al. Unraveling biophysical interactions of radiation pneumonitis in non-small-cell lung cancer via Bayesian network analysis. Radiother Oncol . 2017 Apr 1 ;123(1):85–92.
4. Agrahari R, Foroushani A, Docking TR, Chang L, Duns G, Hudoba M, et al. Applications of Bayesian network models in predicting types of hematological malignancies. Sci Rep . 2018 Dec 3;8(1):6951.
5. Zhi-wei J, Zhen-lei Y, Cai-xiu Z, Li-ying W, Jun L, Hong-li W, et al. Comparison of the Network Structural Characteristics of Calcium Signaling Pathway in Cerebral Ischemia after Intervention by Different Components of Chinese Medicine. J Tradit Chinese Med. 2011;31(3):251–5.
6. Lou S, Ren L, Xiao J, Ding Q, Zhang W. Expression profiling based graph-clustering approach to determine renal carcinoma related pathway in response to kidney cancer. Eur Rev Med Pharmacol Sci. 2012;16(6):775–80.
7. Fu C, Deng S, Jin G, Wang X, Yu ZG. Bayesian network model for identification of pathways by integrating protein interaction with genetic interaction data. BMC Syst Biol. 2017;11.
8. Isci S, Ozturk C, Jones J, Otu HH. Pathway analysis of high-throughput biological data within a Bayesian network framework. 2011;27(12):1667–74.
9. Korucuoglu M, Isci S, Ozgur A, Otu HH. Bayesian pathway analysis of cancer microarray data. PLoS One. 2014;9(7):1–8.
10. Spirtes P, Richardson T. Directed Cyclic Graphical Representations of Feedback Models. Proc Elev Conf Uncertain Artif Intell. 1995;1–37.
11. Neapolitan RE. Learning Bayesian networks. first. Chicago: Pearson Prentice Hall; 2004. 291–425 p. 
12. Scutari M. Learning Bayesian Networks with the bnlearn R Package. J Stat Softw. 2010;35(3):1–22.
13. Hamm A, Prenen H, Van Delm W, Di Matteo M, Wenes M, Delamarre E, et al. Tumour-educated circulating monocytes are powerful candidate biomarkers for diagnosis and disease follow-up of colorectal cancer. Gut. 2016;65(6):990–1000.
14. Nagarajan R, Scutari M, Lèbre S. Bayesian Networks in R . New York, NY: Springer New York; 2013 [cited 2018 Apr 17].
15. Tibshirani R. The lasso method for variable selection in the cox model. Stat Med. 1997;16(4):385–95.
16. Buhlmann P, Geer S van de. Statistics for high-dimensional data: Methods, Theory and Applications . Springer Series in Statistics. 2011. 7–34 p.