Gene regulatory network inference
Installation
Introduction and algorithm description
In the previous vignette, we explored all aspects of gene coexpression networks (GCNs), which are represented as undirected weighted graphs. It is undirected because, for a given link between gene A and gene B, we can only say that these genes are coexpressed, but we cannot know whether gene A controls gene B or otherwise. Further, weighted means that some coexpression relationships between gene pairs are stronger than others. In this vignette, we will demonstrate how to infer gene regulatory networks (GRNs) from expression data with BioNERO. GRNs display interactions between regulators (e.g., transcription factors or miRNAs) and their targets (e.g., genes). Hence, they are represented as directed unweighted graphs.
Numerous algorithms have been developed to infer GRNs from expression
data. However, the algorithm performances are highly dependent on the
benchmark data set. To solve this uncertainty, Marbach et al. (2012) proposed the application
of the “wisdom of the crowds” principle to GRN inference. This
approach consists in inferring GRNs with different algorithms, ranking
the interactions identified by each method, and calculating the average
rank for each interaction across all algorithms used. This way, we can
have consensus, high-confidence edges to be used in biological
interpretations. For that, BioNERO implements three popular
algorithms: GENIE3 (Huynh-Thu et al.
2010), ARACNE (Margolin et al.
2006) and CLR (Faith et al.
2007).
Data preprocessing
Before inferring the GRN, we will preprocess the expression data the same way we did in the previous vignette.
Gene regulatory network inference
BioNERO requires only 2 objects for GRN inference: the
expression data (SummarizedExperiment, matrix or data
frame) and a character vector of regulators
(transcription factors or miRNAs). The transcription factors used in
this vignette were downloaded from PlantTFDB 4.0 (Jin et al. 2017).
data(zma.tfs)
head(zma.tfs)
## Gene Family
## 6 Zm00001d022525 Dof
## 25 Zm00001d037605 GATA
## 28 Zm00001d049540 NAC
## 45 Zm00001d042287 MYB
## 46 Zm00001d042288 NAC
## 54 Zm00001d039371 TCPConsensus GRN inference
Inferring GRNs based on the wisdom of the crowds principle
can be done with a single function: exp2grn(). This
function will infer GRNs with GENIE3, ARACNE and CLR, calculate average
ranks for each interaction and filter the resulting network based on the
optimal scale-free topology (SFT) fit. In the filtering step, n
different networks are created by subsetting the top n
quantiles. For instance, if a network of 10,000 edges is given as input
with nsplit = 10, 10 different networks will be created:
the first with 1,000 edges, the second with 2,000 edges, and so on, with
the last network being the original input network. Then, for each
network, the function will calculate the SFT fit and select the best
fit.
# Using 10 trees for demonstration purposes. Use the default: 1000
grn <- exp2grn(
exp = final_exp,
regulators = zma.tfs$Gene,
nTrees = 10
)
## The top number of edges that best fits the scale-free topology is 247
head(grn)
## Regulator Target
## 290 Zm00001d041474 Zm00001d018986
## 280 Zm00001d041474 Zm00001d006602
## 281 Zm00001d041474 Zm00001d006942
## 325 Zm00001d044315 Zm00001d043497
## 65 Zm00001d013777 Zm00001d046996
## 252 Zm00001d038832 Zm00001d021147Algorithm-specific GRN inference
This section is directed to users who, for some reason (e.g., comparison, exploration), want to infer GRNs with particular algorithms. The available algorithms are:
GENIE3: a regression-tree based algorithm that decomposes the prediction of GRNs for n genes into n regression problems. For each regression problem, the expression profile of a target gene is predicted from the expression profiles of all other genes using random forests (default) or extra-trees.
# Using 10 trees for demonstration purposes. Use the default: 1000
genie3 <- grn_infer(
final_exp,
method = "genie3",
regulators = zma.tfs$Gene,
nTrees = 10)
head(genie3)
## Node1 Node2 Weight
## 20352 Zm00001d041474 Zm00001d017881 0.5439514
## 41340 Zm00001d034751 Zm00001d037111 0.5322394
## 13037 Zm00001d034751 Zm00001d012407 0.4348469
## 13207 Zm00001d045323 Zm00001d012513 0.4203583
## 55378 Zm00001d028432 Zm00001d048693 0.4071160
## 50200 Zm00001d013777 Zm00001d044212 0.3957483
dim(genie3)
## [1] 60136 3ARACNE: information-theoretic algorithm that aims to remove indirect interactions inferred by coexpression.
aracne <- grn_infer(final_exp, method = "aracne", regulators = zma.tfs$Gene)
head(aracne)
## Node1 Node2 Weight
## 23861 Zm00001d038832 Zm00001d021147 1.789818
## 1758 Zm00001d038832 Zm00001d000432 1.692232
## 11337 Zm00001d038832 Zm00001d011086 1.692232
## 27014 Zm00001d011139 Zm00001d024274 1.674840
## 51070 Zm00001d011139 Zm00001d045069 1.658043
## 28387 Zm00001d038832 Zm00001d025784 1.641802
dim(aracne)
## [1] 411 3CLR: extension of the relevance networks algorithm that uses mutual information to identify regulatory interactions.
clr <- grn_infer(final_exp, method = "clr", regulators = zma.tfs$Gene)
head(clr)
## Node1 Node2 Weight
## 26302 Zm00001d046937 Zm00001d023376 12.70216
## 11267 Zm00001d046937 Zm00001d011080 12.25336
## 12540 Zm00001d041474 Zm00001d012007 10.74023
## 51019 Zm00001d042263 Zm00001d045042 10.50925
## 17810 Zm00001d041474 Zm00001d015811 10.33216
## 29278 Zm00001d046937 Zm00001d026632 10.20075
dim(clr)
## [1] 26657 3Users can also infer GRNs with the 3 algorithms at once using the
function exp_combined(). The resulting edge lists are
stored in a list of 3 elements. 1
grn_list <- grn_combined(final_exp, regulators = zma.tfs$Gene, nTrees = 10)
head(grn_list$genie3)
## Node1 Node2 Weight
## 12013 Zm00001d041474 Zm00001d011541 0.4629469
## 30418 Zm00001d046568 Zm00001d027841 0.4289222
## 33403 Zm00001d041474 Zm00001d030748 0.4140894
## 6910 Zm00001d044315 Zm00001d006725 0.4103733
## 22057 Zm00001d041474 Zm00001d018986 0.4020641
## 45153 Zm00001d034751 Zm00001d039733 0.3935705
head(grn_list$aracne)
## Node1 Node2 Weight
## 23861 Zm00001d038832 Zm00001d021147 1.789818
## 1758 Zm00001d038832 Zm00001d000432 1.692232
## 11337 Zm00001d038832 Zm00001d011086 1.692232
## 27014 Zm00001d011139 Zm00001d024274 1.674840
## 51070 Zm00001d011139 Zm00001d045069 1.658043
## 28387 Zm00001d038832 Zm00001d025784 1.641802
head(grn_list$clr)
## Node1 Node2 Weight
## 26302 Zm00001d046937 Zm00001d023376 12.70216
## 11267 Zm00001d046937 Zm00001d011080 12.25336
## 12540 Zm00001d041474 Zm00001d012007 10.74023
## 51019 Zm00001d042263 Zm00001d045042 10.50925
## 17810 Zm00001d041474 Zm00001d015811 10.33216
## 29278 Zm00001d046937 Zm00001d026632 10.20075Gene regulatory network analysis
After inferring the GRN, BioNERO allows users to perform
some common downstream analyses.
Hub gene identification
GRN hubs are defined as the top 10% most highly connected regulators,
but this percentile is flexible in BioNERO.2 They can be identified
with get_hubs_grn().
hubs <- get_hubs_grn(grn)
hubs
## Gene Degree
## 1 Zm00001d038832 16
## 2 Zm00001d041474 13
## 3 Zm00001d046937 13
## 4 Zm00001d011139 12
## 5 Zm00001d052229 11
## 6 Zm00001d013777 10
## 7 Zm00001d039989 10
## 8 Zm00001d038227 10
## 9 Zm00001d030617 10
## 10 Zm00001d044315 9
## 11 Zm00001d003822 9
## 12 Zm00001d020020 9
## 13 Zm00001d046568 9
## 14 Zm00001d010227 9
## 15 Zm00001d025339 8
## 16 Zm00001d028974 8
## 17 Zm00001d042267 7
## 18 Zm00001d014377 7
## 19 Zm00001d054038 6
## 20 Zm00001d042263 6
## 21 Zm00001d035440 6
## 22 Zm00001d036148 6
## 23 Zm00001d031655 6
## 24 Zm00001d034751 6
## 25 Zm00001d018081 6
## 26 Zm00001d027957 5Network visualization
GRNs can also be visualized interactively for exploratory purposes.
Finally, BioNERO can also be used for visualization and
hub identification in protein-protein (PPI) interaction networks. The
functions get_hubs_ppi() and plot_ppi() work
the same way as their equivalents for GRNs (get_hubs_grn()
and plot_grn()).
Session information
This vignette was created under the following conditions:
sessionInfo()
## 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] BioNERO_1.15.0 BiocStyle_2.33.1
##
## loaded via a namespace (and not attached):
## [1] RColorBrewer_1.1-3 ggdendro_0.2.0
## [3] sys_3.4.3 rstudioapi_0.17.1
## [5] jsonlite_1.8.9 shape_1.4.6.1
## [7] NetRep_1.2.7 magrittr_2.0.3
## [9] farver_2.1.2 rmarkdown_2.28
## [11] GlobalOptions_0.1.2 zlibbioc_1.51.2
## [13] vctrs_0.6.5 memoise_2.0.1
## [15] base64enc_0.1-3 htmltools_0.5.8.1
## [17] S4Arrays_1.5.11 dynamicTreeCut_1.63-1
## [19] SparseArray_1.5.45 Formula_1.2-5
## [21] sass_0.4.9 bslib_0.8.0
## [23] htmlwidgets_1.6.4 plyr_1.8.9
## [25] impute_1.79.0 cachem_1.1.0
## [27] networkD3_0.4 buildtools_1.0.0
## [29] igraph_2.1.1 lifecycle_1.0.4
## [31] ggnetwork_0.5.13 iterators_1.0.14
## [33] pkgconfig_2.0.3 Matrix_1.7-1
## [35] R6_2.5.1 fastmap_1.2.0
## [37] GenomeInfoDbData_1.2.13 MatrixGenerics_1.17.1
## [39] clue_0.3-65 digest_0.6.37
## [41] colorspace_2.1-1 patchwork_1.3.0
## [43] AnnotationDbi_1.69.0 S4Vectors_0.43.2
## [45] GENIE3_1.27.0 Hmisc_5.2-0
## [47] GenomicRanges_1.57.2 RSQLite_2.3.7
## [49] labeling_0.4.3 fansi_1.0.6
## [51] mgcv_1.9-1 httr_1.4.7
## [53] abind_1.4-8 compiler_4.4.1
## [55] withr_3.0.2 bit64_4.5.2
## [57] doParallel_1.0.17 htmlTable_2.4.3
## [59] backports_1.5.0 BiocParallel_1.39.0
## [61] DBI_1.2.3 intergraph_2.0-4
## [63] highr_0.11 MASS_7.3-61
## [65] DelayedArray_0.31.14 rjson_0.2.23
## [67] tools_4.4.1 foreign_0.8-87
## [69] nnet_7.3-19 glue_1.8.0
## [71] nlme_3.1-166 grid_4.4.1
## [73] checkmate_2.3.2 cluster_2.1.6
## [75] reshape2_1.4.4 sva_3.53.0
## [77] generics_0.1.3 gtable_0.3.6
## [79] preprocessCore_1.67.1 data.table_1.16.2
## [81] WGCNA_1.73 utf8_1.2.4
## [83] XVector_0.45.0 BiocGenerics_0.51.3
## [85] ggrepel_0.9.6 foreach_1.5.2
## [87] pillar_1.9.0 stringr_1.5.1
## [89] limma_3.61.12 genefilter_1.87.0
## [91] circlize_0.4.16 splines_4.4.1
## [93] dplyr_1.1.4 lattice_0.22-6
## [95] survival_3.7-0 bit_4.5.0
## [97] annotate_1.85.0 tidyselect_1.2.1
## [99] locfit_1.5-9.10 GO.db_3.20.0
## [101] ComplexHeatmap_2.21.1 maketools_1.3.1
## [103] Biostrings_2.73.2 knitr_1.48
## [105] gridExtra_2.3 IRanges_2.39.2
## [107] edgeR_4.3.21 SummarizedExperiment_1.35.5
## [109] RhpcBLASctl_0.23-42 stats4_4.4.1
## [111] xfun_0.48 Biobase_2.65.1
## [113] statmod_1.5.0 matrixStats_1.4.1
## [115] stringi_1.8.4 UCSC.utils_1.1.0
## [117] statnet.common_4.10.0 yaml_2.3.10
## [119] minet_3.63.0 evaluate_1.0.1
## [121] codetools_0.2-20 tibble_3.2.1
## [123] BiocManager_1.30.25 cli_3.6.3
## [125] rpart_4.1.23 xtable_1.8-4
## [127] munsell_0.5.1 jquerylib_0.1.4
## [129] network_1.18.2 Rcpp_1.0.13
## [131] GenomeInfoDb_1.41.2 coda_0.19-4.1
## [133] png_0.1-8 XML_3.99-0.17
## [135] fastcluster_1.2.6 parallel_4.4.1
## [137] ggplot2_3.5.1 blob_1.2.4
## [139] scales_1.3.0 crayon_1.5.3
## [141] GetoptLong_1.0.5 rlang_1.1.4
## [143] KEGGREST_1.45.1References
NOTE: Under the hood,
exp2grn()usesexp_combined()followed by averaging ranks withgrn_average_rank()and filtering withgrn_filter().↩︎NOTE: Remember: GRNs are represented as directed graphs. This implies that only regulators are taken into account when identifying hubs. The goal here is to identify regulators (e.g., transcription factors) that control the expression of several genes.↩︎