A statistical method for measuring activation of gene regulatory networks

2.1 𝒢 statistic distribution

This section addresses the deduction of the 𝒢 statistic, given by Eq. (1), distribution. Let P_s be a random variable that gives the p value from the hypothesis test for the edge s in the graph. Assuming the null hypothesis that the gene regulatory network represented by graph G is inactive, P_s may assume any value within the [0, 1] interval independently of the other edges in the graph. In other words, we can assume P_s ∼ U(0, 1), meaning that P_s follows a uniform distribution over [0, 1] for all s = 1, 2, … , S.

Taking Y = −log(P_s), we find that

indicating that f_Y(y) = e^−y, implying that Y = −log(P_s) follows an exponential distribution for the parameter 1 (one) for all s = 1, …, S. Then, defining a new random variable Z=∑s=1S−log⁡(Ps) and using the fact that the moment generating function (m.g.f.) of an exponential distribution with the parameter α is given by

we can find the distribution of Z by taking its m.g.f.

However, g_Z(t) = (1/(1 − t))^S is the m.g.f. of a gamma (Γ) distribution with parameters S and 1 (one). Thus, ∑s=1S−log⁡(Ps)∼Γ(S,1) and has a density probability function given by

Finally, we can see that 𝒢 = Z/S, yielding

and

which is the probability density function of a gamma distribution with parameters S and 1/S. Thus, we have shown that 𝒢 ∼ Γ(S, 1/S). Figure 1 illustrates this density function for certain values of S.

Figure 1:

Probability density functions from Γ(S,1/S) for S=1,2,7,15.

Note that a gamma distribution with parameters α and β has a mean and variance given by αβ and αβ2, respectively. We can see that 𝒢 has a mean of one and a variance of 1/S for any value of S. Note that the variance decreases as the number of edges increases, which can be observed in Figure 1. Even with the fact that 𝒢 has a mean of one for all values of S, the values obtained for this statistic may not be comparable between different networks due to differences in the variance relative to the number of edges.

To validate this theoretical probability distribution, we used a gene co-expression network generated from a real dataset previously analyzed by Gomes et al. (2005) and random permutations as described in the Results section. We next present the construction of hypothesis tests to verify the profile of activation and inactivation in any network G.

One problem associated with the distribution of the 𝒢 statistic is the fact that it is deduced with the assumption of independence between the edges on the graph G. While this could be a reasonable assumption for edges that do not have a common node, this is not true for edges in this situation. In Section 2 from Supplementary Material, we present a simulation study to evaluate these assumption. In summary, we show that breaking this assumption do not significantly affects the distribution of 𝒢, but can result in a little increase in the probability of false positives occurrence (Sup. Figures 3 and 4).

2.2 Hypothesis test construction

From the 𝒢 statistic, given by Eq. (1), as well as its probability distribution under the null hypothesis of network inactivation, it is possible to perform two types of statistical hypothesis tests aimed at evaluating the status of the network. The more important of these is testing whether the network G is active under the proposed biological model (H_A) versus its inactivation (H₀). In other words, we can test the following hypothesis:

Supposing that H₀ is true, 𝒢 follows a gamma distribution as described above, and its value should remain relatively close to one. If, however, H_A is true, the 𝒢 statistic should return large positive values [greater than one, as expected from the Γ(S, 1/S) distribution], providing evidence for H_A and indicating statistical significance that the network appears do be “turned on” or active at that particular moment in time. This is therefore a right-tailed test.

For the second case, we can test whether the network is in a state of activation (in H_A), but in an opposite manner relative to the biological knowledge represented by the network. It is therefore possible to test the two hypothesis:

In this case, if H_A is true, 𝒢 results in small positive values (near zero); thus, this test is left-tailed. However, in both situations, the null distribution of the 𝒢 statistic is given by a gamma distribution parameterized by S and 1/S, where S is the total number of edges in the biological network under consideration. In this manner, the p values for both tests can be easily computed from this distribution.

The first test presented above is paramount, in the sense that it evaluates network activation according to the proposed biological model. However, it is nevertheless possible to test a kind of opposite activation of the network in cases were the resulting 𝒢 values are very close to zero. This can be done using the second test presented above.

For both tests, if several gene regulatory networks are to be tested simultaneously it is important to do the p values adjustment for multiple tests. There are several methods for this procedure, and here we used the false discovery rate (FDR) to this end (Benjamini & Hochberg, 1995).

This statistical method has been implemented in the activeNet function of the maigesPack R package, which integrates other gene expression data analysis methods and is part of the BioConductor project (available at http://www.bioconductor.org).

2.3 Microarray dataset

Gomes et al. (2005) used microarrays to analyze the gene expression profile of gastric-esophageal cells. The expression of 4800 cDNA clones representing 4565 full-length human genes were observed in the RNA samples of 71 patients: 39 esophagus and gastric-esophageal junction (GEJ) samples [9 normal esophageal mucosa, 6 esophagitis mucosa, 10 Barrett’s mucosa (4 of the long type and 6 of the short type), 5 adenocarcinomas of the esophagus, and 9 adenocarcinomas of the GEJ] and 32 stomach samples (6 normal body and antrum mucosa, 5 normal cardiac mucosa, 9 intestinal metaplasia mucosa, 7 samples of intestinal-type adenocarcinoma, and 5 samples of diffuse-type carcinoma). More details can be obtained in the original work.

Normalized data, as described in the original article, were used to validate the 𝒢 statistic distribution and to exemplify the hypothesis tests using a small gene network and two KEGG pathways, as described below.

The KEGG pathways investigated here were cytokine-cytokine receptor interaction (id 04060) and glycerolipid metabolism (id 00561), which were obtained from KEGG database on 08/17/2015 using retrieveKGML function and parsed to graphs using the parseKGML2Graph function, both from the KEGGgraph package available for the R software for statistical computing (Ihaka & Gentleman, 1996). These two pathways were used because they appeared to be activated by Segal’s method in the original work (Gomes et al., 2005).

The cytokine-cytokine receptor interaction network originally had 265 nodes and 254 edges. However, after matching the network with genes available in the dataset, we identified 45 nodes and 11 edges for this network, yielding coverages of approximately 17% and 4% for nodes and edges, respectively. For glycerolipid metabolism, we began with 59 nodes and 763 edges. However, a comparison with the gene expression data resulted in 17 nodes and 57 edges, approximately 29% and 7% coverage, respectively.

Thus, we used these two subnetworks obtained from the KEGG database mainly for the purpose of illustrating the operation of the method proposed here. Biologically, it is not reliable to discuss in depth the results obtained due to the poor coverage of both networks for this particular database. We also used a public database obtained through GEO repository to replicate the analysis of these two KEGG networks using an alternative gene expression data, as presented in the Section 1 from Supplementary Material.

3.1 Validation of 𝒢 distribution

One way to validate the theoretical probabilistic distribution of 𝒢 is data permutation, where we shuffled gene observations and recalculate the statistic for a network known to be active. In the process, we break any associations present into the network. Thus, for each data permutation, the 𝒢 statistic given by Eq. (1) is recalculated, and the process is repeated a number of times to estimate the empirical distribution for the 𝒢 statistic and to compare it with the theoretical distribution given in the Methods section.

For this, we require an activated network, so we constructed a co-expression network using 30 random genes with pairwise absolute correlation values above 0.75 from microarray dataset cited above (Gomes et al., 2005). This number of 30 random genes were defined empirically, as it is near the mean nodes for the two KEGG subnetworks used to illustrate the method. This generated a network with 30 nodes and 16 edges, with ten nodes without any in or out edges, as seen in Figure 2 (nodes without edges are not represented in the figure). Because the network was obtained from the dataset, all correlation values for each edge are real. 𝒢 statistic calculation resulted in an approximate value of 4.4, which is a rather unlikely value from a Γ(16, 1/16) distribution, leading to a p value very close to zero (on the order of 10⁻¹⁵) in the activation test for this co-expression network, as would be expected since it was constructed from the gene expression dataset.

Figure 2:

Gene co-expression network constructed for 30 random genes. Ten nodes without edges are not represented in the figure. Edge labels represent estimated Pearson correlation values.

Then, gene expression values for this network were randomly permuted 2000 times, recalculating the 𝒢 value each time. Figure 3 presents a histogram of these 2000 values, which gives the empirical distribution for 𝒢 superimposed with the density probability function of Γ (16, 1/16) in dashed line. As we can see in the figure, there is a good fit between the histogram and the theoretical density, indicating that the 𝒢 statistic really seems to follow the Γ (16, 1/16) distribution.

Figure 3:

Histogram from the G statistic applied over a co-expression network with 30 nodes and 16 edges after 2 thousand permutations from gene expression values. Indicated by the dashed line is the Γ(16,1/16) theoretical distribution.

Another way to verify the theoretical model is to use qqplot. This type of graphic allows a comparison between quantiles from two groups of observations to verify whether they follow the same pattern. Thus, we can compare the sample quantile from the permuted 𝒢 statistic with that of Γ (16, 1/16). Figure 4 shows the good fit of the qqplot.

The main focus of this section was to show the good fit between the theoretical probability distribution deduced previously and the empirical one estimated by permutation methods and none of the results presented here are meant to be biologically meaningful or relevant.

3.2 Activation measuring of a small gene network

This subsection presents the results of the classification of a small and simple hypothetical gene regulatory network. According to the previous results of Gomes et al. (2005), it is possible to suppose that an increase in the expression of the pro-inflammatory genes CCL20, CCL18 and IFNAR2 in normal tissues must be accompanied by an increase in expression of genes ADH1B, AKR1B10, ALDH3A2 and IL1R2 in order to promote the detoxification of products generated during the inflammatory reaction. When these interactions are broken, the accumulation of toxic products synthesized via the inflammatory process – in particular, the production of aldehyde – can favor the accumulation of mutations in DNA and an eventual malignant transformation.

Figure 4:

Qqplot comparing the sample quantiles from the permuted G statistic with the Γ(16,1/16) theoretical ones.

From the idea behind the above explanation, it is possible to suppose that the expression values observed for the genes CCL20, CCL18 and IFNAR2 must be positively correlated with that ones for the genes ADH1B, AKR1B10, ALDH3A2 and IL1R2 in normal tissues and not associated for tumor ones; which can be represented by the graph in Figure 5.

Figure 5:

Graph representing a simple gene regulatory network, where the expression levels for genes CCL20, CCL18 and IFNAR2 must be positively associated with the the ones for genes ADH1B, AKR1B10, ALDH3A2 and IL1R2 in normal tissues.

The graph at Figure 5 represents a gene regulatory network that could be active in normal tissues, both for stomach and esophagus and inactive for tumors. Obviously, some additional biological experimentation is needed to validate this hypothesis, but even so, this graph can be used to test the functionality of the statistical methodology proposed here.

In this way, we applied the 𝒢 statistic and hypothesis test proposed above to verify the activation profile of this network using the original dataset, comparing tissue types classified as Normal (normal tissues from esophagus, GEJ and stomach), Metaplasia (Barrett’s esophagus and intestinal metaplasia from the stomach) and Adenocarcinoma (all adenocarcinomas from the esophagus and stomach, both intestinal and diffuse types). We found that for a p value less than 0.001 the 𝒢 statistic of approximately 3.4 in normal tissues is significant, and it is possible to say that there is a strong statistical evidence for network activation for this biological condition, while it is lacking for metaplasia and adenocarcinomas (results in Table 1), as hypothesized above.

From these results, we could conclude that our findings seems to be in accordance with the biological explanation underlying the construction of the proposed network and demonstrate the great potential of the statistical method proposed in this work. Biologically, as mentioned above, some additional experimentation is needed to validate these results.

3.3 Activation measuring of two KEGG networks

Thus far, we have presented a simulation-based result for the validation of a probabilistic model of the 𝒢 statistic using a co-expression network and the resulting classification of a simple hypothetical network based on the results of original work by Gomes et al. (2005), whose behavior could be considered known among normal and tumor tissues.

Now we present the results for the activation measuring of two networks obtained from KEGG (cytokine-cytokine receptor interaction and glycerolipid metabolism) that appeared to be significantly altered in the original work using the gene group classification method proposed originally by Segal et al. (2004). Cytokine-cytokine receptor interaction was considered inactive in normal esophagus samples, whereas glycerolipid metabolism was inactive in adenocarcinomas and active in intestinal metaplasia of the stomach.

Using these two KEGG pathways presents a great opportunity for testing the statistical methodology proposed here under a more realistic scenario by calculating the 𝒢 statistic and applying the respective activation test to both networks. To this end, we used the two subnetworks whose genes were present in the database. Figure 6 and Figure 7 show graphs representing these two subnetworks. In Figure 6, some nodes without edges are not represented in the figure.

Figure 6:

Graph representing the subnetwork for cytokine-cytokine receptor interaction, after the match with genes from the dataset. Some nodes without edges are not represented in the figure.

For the cytokine-cytokine receptor interaction pathway, Table 2 shows that, at 5% significance level, the 𝒢 statistic presented significant values only for normal tissue type (𝒢 approximately equal to 1.76, 1.27 and 1.10, for normal, metaplasia and adenocarcinomas, respectively), indicating that the subnetwork for this pathway seems to be activated in normal tissue and inactivated for the other two tissue types. In a similar manner, the glycerolipid metabolism subnetwork also presented similar results for the three tissue type (𝒢 approximately equal to 1.96, 0.92 and 0.97, for normal, metaplasia and adenocarcinomas, respectively). These results can be found in Table 3.

The method proposed by Segal et al. (2004) and used in the work of Gomes et al. (2005) classifies each sample type as active or inactive based on the number of genes over- or under-expressed, respectively, using a hypergeometric distribution to calculate the significance of the result. Thus, it may be classified as a kind of GSEA method, which does not consider the interrelation between genes that constitute the group. So the results obtained by Gomes et al. (2005) are not comparable with those ones presented in this section.

Figure 7:

Graph representing the subnetwork for glycerolipid metabolism, after the match with genes from the dataset.

It is also important to note that any interpretations regarding the biology of these two networks based on the results obtained here should be made with caution, because only a small portion of the original edges (and thus, number of genes) are represented by the genes available in this database. Furthermore our results need biological validation to be discussed in more depth.

The same type of analysis for both KEGG networks were done using a public database, available through GEO (access code GSE33651). The subnetworks identified in this case are shown in Sup. Figures 1 and 2, and the results obtained comparing normal and tumor gastric tissues are seen at Sup. Tables 1 and 2 (Supplementary Material).