# Gene regulatory network inference

Fabricio Almeida-Silva1 and Thiago Motta Venancio1

1Universidade Estadual do Norte Fluminense Darcy Ribeiro, RJ, Brazil

#### 4 September 2022

# Contents

# 1 Installation

    if(!requireNamespace('BiocManager', quietly = TRUE))
      install.packages('BiocManager')

    BiocManager::install("BioNERO")

    # Load package after installation
    library(BioNERO)
    set.seed(123) # for reproducibility

# 2 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).

# 3 Data preprocessing

Before inferring the GRN, we will preprocess the expression data the same way we did in the previous vignette.

    # Load example data set
    data(zma.se)

    # Preprocess the expression data
    final_exp <- exp_preprocess(
        zma.se,
        min_exp = 10,
        variance_filter = TRUE,
        n = 2000
    )
    ## Number of removed samples: 1

# 4 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    TCP

## 4.1 Consensus 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
    ## 291 Zm00001d041474 Zm00001d027374
    ## 304 Zm00001d042263 Zm00001d045042
    ## 280 Zm00001d041474 Zm00001d006602
    ## 281 Zm00001d041474 Zm00001d006942
    ## 289 Zm00001d041474 Zm00001d018668
    ## 145 Zm00001d028432 Zm00001d014349

## 4.2 Algorithm-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
    ## 49836 Zm00001d039989 Zm00001d043942 0.4559930
    ## 37975 Zm00001d030617 Zm00001d034368 0.4293781
    ## 34029 Zm00001d038832 Zm00001d031640 0.4287395
    ## 8526  Zm00001d045323 Zm00001d008594 0.4084538
    ## 39247 Zm00001d045323 Zm00001d035094 0.4084167
    ## 34629 Zm00001d052229 Zm00001d031963 0.3899265
    dim(genie3)
    ## [1] 60140     3

ARACNE: 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   3

CLR: 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     3

Users 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 1 1 NOTE: Under the hood, `exp2grn()` uses `exp_combined()` followed by averaging ranks with `grn_average_rank()` and filtering with `grn_filter()`.

    grn_list <- grn_combined(final_exp, regulators = zma.tfs$Gene, nTrees = 10)
    head(grn_list$genie3)
    ##                Node1          Node2    Weight
    ## 12013 Zm00001d041474 Zm00001d011541 0.4629469
    ## 22057 Zm00001d041474 Zm00001d018986 0.4020641
    ## 22456 Zm00001d028432 Zm00001d019479 0.3941906
    ## 21133 Zm00001d038832 Zm00001d018403 0.3926708
    ## 56194 Zm00001d038832 Zm00001d049552 0.3923130
    ## 45622 Zm00001d036148 Zm00001d040163 0.3895359
    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.20075

# 5 Gene regulatory network analysis

After inferring the GRN, `BioNERO` allows users to perform some common downstream analyses.

## 5.1 Hub gene identification

GRN hubs are defined as the top 10% most highly connected regulators, but this percentile is flexible in `BioNERO`.2 2 2 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. They can be identified with `get_hubs_grn()`.

    hubs <- get_hubs_grn(grn)
    hubs
    ##              Gene Degree
    ## 1  Zm00001d038832     19
    ## 2  Zm00001d041474     13
    ## 3  Zm00001d025339     11
    ## 4  Zm00001d046937     11
    ## 5  Zm00001d052229     11
    ## 6  Zm00001d010227     11
    ## 7  Zm00001d030617     11
    ## 8  Zm00001d039989     10
    ## 9  Zm00001d038227     10
    ## 10 Zm00001d011139      9
    ## 11 Zm00001d044315      9
    ## 12 Zm00001d028974      9
    ## 13 Zm00001d013777      9
    ## 14 Zm00001d003822      8
    ## 15 Zm00001d046568      8
    ## 16 Zm00001d020020      7
    ## 17 Zm00001d014377      7
    ## 18 Zm00001d034751      7
    ## 19 Zm00001d054038      7
    ## 20 Zm00001d031655      6
    ## 21 Zm00001d042267      6
    ## 22 Zm00001d036148      6
    ## 23 Zm00001d018081      6
    ## 24 Zm00001d035440      6
    ## 25 Zm00001d027957      6
    ## 26 Zm00001d042263      5
    ## 27 Zm00001d002654      5

## 5.2 Network visualization

    plot_grn(grn)

GRNs can also be visualized interactively for exploratory purposes.

    plot_grn(grn, interactive = TRUE, dim_interactive = c(500,500))

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.2.1 (2022-06-23 ucrt)
    ## Platform: x86_64-w64-mingw32/x64 (64-bit)
    ## Running under: Windows Server x64 (build 20348)
    ##
    ## Matrix products: default
    ##
    ## locale:
    ## [1] LC_COLLATE=C
    ## [2] LC_CTYPE=English_United States.utf8
    ## [3] LC_MONETARY=English_United States.utf8
    ## [4] LC_NUMERIC=C
    ## [5] LC_TIME=English_United States.utf8
    ##
    ## attached base packages:
    ## [1] stats     graphics  grDevices utils     datasets  methods   base
    ##
    ## other attached packages:
    ## [1] BioNERO_1.4.2    BiocStyle_2.24.0
    ##
    ## loaded via a namespace (and not attached):
    ##   [1] backports_1.4.1             circlize_0.4.15
    ##   [3] Hmisc_4.7-1                 plyr_1.8.7
    ##   [5] igraph_1.3.4                splines_4.2.1
    ##   [7] GENIE3_1.18.0               BiocParallel_1.30.3
    ##   [9] GenomeInfoDb_1.32.3         ggnetwork_0.5.10
    ##  [11] ggplot2_3.3.6               sva_3.44.0
    ##  [13] digest_0.6.29               foreach_1.5.2
    ##  [15] htmltools_0.5.3             magick_2.7.3
    ##  [17] GO.db_3.15.0                fansi_1.0.3
    ##  [19] magrittr_2.0.3              checkmate_2.1.0
    ##  [21] memoise_2.0.1               cluster_2.1.4
    ##  [23] doParallel_1.0.17           limma_3.52.2
    ##  [25] ComplexHeatmap_2.12.1       fastcluster_1.2.3
    ##  [27] Biostrings_2.64.1           annotate_1.74.0
    ##  [29] matrixStats_0.62.0          jpeg_0.1-9
    ##  [31] colorspace_2.0-3            blob_1.2.3
    ##  [33] ggrepel_0.9.1               xfun_0.32
    ##  [35] dplyr_1.0.10                crayon_1.5.1
    ##  [37] RCurl_1.98-1.8              jsonlite_1.8.0
    ##  [39] genefilter_1.78.0           impute_1.70.0
    ##  [41] survival_3.4-0              iterators_1.0.14
    ##  [43] glue_1.6.2                  gtable_0.3.1
    ##  [45] zlibbioc_1.42.0             XVector_0.36.0
    ##  [47] GetoptLong_1.0.5            DelayedArray_0.22.0
    ##  [49] shape_1.4.6                 BiocGenerics_0.42.0
    ##  [51] abind_1.4-5                 scales_1.2.1
    ##  [53] edgeR_3.38.4                DBI_1.1.3
    ##  [55] Rcpp_1.0.9                  xtable_1.8-4
    ##  [57] htmlTable_2.4.1             clue_0.3-61
    ##  [59] foreign_0.8-82              bit_4.0.4
    ##  [61] preprocessCore_1.58.0       Formula_1.2-4
    ##  [63] stats4_4.2.1                htmlwidgets_1.5.4
    ##  [65] httr_1.4.4                  RColorBrewer_1.1-3
    ##  [67] farver_2.1.1                pkgconfig_2.0.3
    ##  [69] XML_3.99-0.10               nnet_7.3-17
    ##  [71] sass_0.4.2                  deldir_1.0-6
    ##  [73] locfit_1.5-9.6              utf8_1.2.2
    ##  [75] dynamicTreeCut_1.63-1       labeling_0.4.2
    ##  [77] tidyselect_1.1.2            rlang_1.0.5
    ##  [79] reshape2_1.4.4              AnnotationDbi_1.58.0
    ##  [81] munsell_0.5.0               tools_4.2.1
    ##  [83] cachem_1.0.6                cli_3.3.0
    ##  [85] generics_0.1.3              RSQLite_2.2.16
    ##  [87] statnet.common_4.6.0        evaluate_0.16
    ##  [89] stringr_1.4.1               fastmap_1.1.0
    ##  [91] yaml_2.3.5                  RhpcBLASctl_0.21-247.1
    ##  [93] knitr_1.40                  bit64_4.0.5
    ##  [95] purrr_0.3.4                 KEGGREST_1.36.3
    ##  [97] nlme_3.1-159                compiler_4.2.1
    ##  [99] rstudioapi_0.14             png_0.1-7
    ## [101] minet_3.54.0                tibble_3.1.8
    ## [103] statmod_1.4.37              bslib_0.4.0
    ## [105] stringi_1.7.8               highr_0.9
    ## [107] lattice_0.20-45             Matrix_1.4-1
    ## [109] vctrs_0.4.1                 pillar_1.8.1
    ## [111] lifecycle_1.0.1             networkD3_0.4
    ## [113] BiocManager_1.30.18         jquerylib_0.1.4
    ## [115] GlobalOptions_0.1.2         data.table_1.14.2
    ## [117] bitops_1.0-7                patchwork_1.1.2
    ## [119] GenomicRanges_1.48.0        R6_2.5.1
    ## [121] latticeExtra_0.6-30         bookdown_0.28
    ## [123] network_1.17.2              gridExtra_2.3
    ## [125] IRanges_2.30.1              codetools_0.2-18
    ## [127] assertthat_0.2.1            SummarizedExperiment_1.26.1
    ## [129] rjson_0.2.21                S4Vectors_0.34.0
    ## [131] GenomeInfoDbData_1.2.8      intergraph_2.0-2
    ## [133] mgcv_1.8-40                 parallel_4.2.1
    ## [135] grid_4.2.1                  rpart_4.1.16
    ## [137] NetRep_1.2.4                coda_0.19-4
    ## [139] rmarkdown_2.16              MatrixGenerics_1.8.1
    ## [141] Cairo_1.6-0                 ggnewscale_0.4.7
    ## [143] Biobase_2.56.0              WGCNA_1.71
    ## [145] base64enc_0.1-3             interp_1.1-3

# References

Faith, Jeremiah J., Boris Hayete, Joshua T. Thaden, Ilaria Mogno, Jamey Wierzbowski, Guillaume Cottarel, Simon Kasif, James J. Collins, and Timothy S. Gardner. 2007. “Large-scale mapping and validation of Escherichia coli transcriptional regulation from a compendium of expression profiles.” PLoS Biology 5 (1): 0054–0066. https://doi.org/10.1371/journal.pbio.0050008.

Huynh-Thu, Vân Anh, Alexandre Irrthum, Louis Wehenkel, and Pierre Geurts. 2010. “Inferring regulatory networks from expression data using tree-based methods.” PLoS ONE 5 (9): 1–10. https://doi.org/10.1371/journal.pone.0012776.

Jin, J, F Tian, D C Yang, Y Q Meng, L Kong, J Luo, and G Gao. 2017. “PlantTFDB 4.0: toward a central hub for transcription factors and regulatory interactions in plants.” Nucleic Acids Res 45 (D1): D1040–D1045. https://doi.org/10.1093/nar/gkw982.

Marbach, Daniel, James C. Costello, Robert Küffner, Nicole M. Vega, Robert J. Prill, Diogo M. Camacho, Kyle R. Allison, et al. 2012. “Wisdom of crowds for robust gene network inference.” Nature Methods 9 (8): 796–804. https://doi.org/10.1038/nmeth.2016.

Margolin, Adam A., Ilya Nemenman, Katia Basso, Chris Wiggins, Gustavo Stolovitzky, Riccardo Dalla Favera, and Andrea Califano. 2006. “ARACNE: An algorithm for the reconstruction of gene regulatory networks in a mammalian cellular context.” BMC Bioinformatics 7 (SUPPL.1): 1–15. https://doi.org/10.1186/1471-2105-7-S1-S7.