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
## 290 Zm00001d041474 Zm00001d018986
## 280 Zm00001d041474 Zm00001d006602
## 281 Zm00001d041474 Zm00001d006942
## 325 Zm00001d044315 Zm00001d043497
## 65  Zm00001d013777 Zm00001d046996
## 252 Zm00001d038832 Zm00001d021147

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
## 23436 Zm00001d030617 Zm00001d020827 0.4221923
## 2993  Zm00001d034751 Zm00001d002543 0.4121333
## 15973 Zm00001d003822 Zm00001d014519 0.4074259
## 29683 Zm00001d041474 Zm00001d027374 0.3993255
## 11749 Zm00001d030617 Zm00001d011401 0.3971345
## 40456 Zm00001d045323 Zm00001d036250 0.3932967
dim(genie3)
## [1] 60111     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 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
## 39975 Zm00001d041474 Zm00001d035696 0.5020457
## 3543  Zm00001d027957 Zm00001d003346 0.4897761
## 60968 Zm00001d038832 Zm00001d054089 0.4738519
## 34907 Zm00001d045323 Zm00001d032181 0.4284764
## 21728 Zm00001d044315 Zm00001d018758 0.4268109
## 18551 Zm00001d011139 Zm00001d016549 0.4189310
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 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     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      5

5.2 Network visualization

plot_grn(grn)