Intrachromosomal regulation decay in breast cancer

One relevant problem in contemporary computational biology is the probabilistic inference of the best (i.e., the maximum-likelihood or maximum-entropy) set of regulatory interactions between genes starting from a large –but still partial – data corpus (, as given for instance, in RNA sequencing experiments over whole-genome transcriptomes. We will call this problem the gene regulatory program deconvolution, GRPD. Solving the GRPD problem involves large-scale probabilistic inference in highly noisy data sets, and it thus remains a challenge to common probabilistic modelling approaches.

In order to circumvent these limitations, a number of algorithms – some of them with information theoretical foundations – have been developed. These approaches include mutual information maximization, Markov random fields, use of the data processing inequality, minimum description length, and Kullback–Leibler divergence. Also relevant to the GRPD problem are machine learning techniques, as well as Monte Carlo methods, variational methods, and hidden Markov models or stochastic linear dynamical systems. Common to these latter models is the fact that they are parametrically conditioned on the hidden state vector; the past, present, and future observations are statistically independent [8].

Information theoretical-founded approaches are useful for the tasks of feature selection and network inference. Feature selection methods applied to transcriptomics aim to improve molecular diagnosis and prognosis in complex diseases (such as cancer) by identifying a (minimum) set (a molecular signature) of features that characterize the underlying biological phenomenon. Network inference, in contrast, usually tries to present the full set of statistical dependencies between genes by means of a probabilistic graphical model of a gene regulatory network.

A first step towards solving either the feature selection or network inference subproblems of a GRPD is to have a detailed knowledge of the joint and marginal gene expression probability distributions. In what follows, we present an information theoretical approach to the GRPD problem via mutual information distributions.

Let i ( {1, 2,..., N} and j ( {1, 2,..., N} be two identical sets of N random variables, representing the expression levels of N different genes in a transcriptome. For each duplex Di _{j (} (i, j) (representing a pair of genes), it is possible to define the mutual information function I(i, j) as follows [3]:

Here, I and J are the complete sampling spaces associated with the gene expression levels i and j, respectively – i.e., the sets of all possible values of i and j, within a given (large) experimental data corpus (, associated with a general probability triple ((, F, P). P(i, j) is the joint probability distribution of i and j in (, whereas P(i) and P(j) are the marginal probability distributions of i and j, respectively. As is widely known, the mutual information function I(i, j) quantifies the statistical dependence between two given random variables i and j [3].

We can also define the off-diagonal mutual information, I^†, as follows:

Here (ijis Kronecker’s delta. I^† (i, j) equals the mutual information I(i, j) everywhere, except in the set i ( j, i, j. The purpose of I^†( i, j) is to eliminate self-information from our calculations. From now on, we will drop the † superscript and we will always refer to the off-diagonal mutual information in all of our further calculations.

A GRP is the solution of a GRPD problem, i.e., the full set of interactions among genes that give rise to a transcriptional phenotype. Within the context of the theoretical and experimental settings we have just described, let us define what the solution of a GRPD problem is.

Definition 1. Here we define a gene regulatory program as the functional G[I(i, j)] of all the mutual information distribution functions for a given empirical transcriptomics sampling space Ω.

G [I(i, j)] contains all the information about the statistical dependencies among genes at a transcriptional level. G[I(i, j)] is thus a form of a Markov random field [10, 14], since it considers mutual information distributions at the pairwise sufficiency (hence, Markov) level [12].

Biological phenotypes are the result of a large set of regulatory interactions that are controlled by diverse genomic and epigenomic elements. Perturbation of these elements is involved in the origin and maintenance of the pathological phenotype observed in cancer. One of these elements is the spatial configuration of the genome. As we have previously mentioned, in a recent work [4], we have shown the existence of a major difference in the relationship between gene interactions and physical distance between breast cancer and healthy breast phenotypes.

Starting from the solution G [I(i, j)] of the breast cancer GRPD program, we are in a position to study in a more formal and less constrained way the phenomenon of loss of long-range regulation that our group has previously reported to exist in breast cancer [4].

The functional G [I(i, j)] is composed of the complete gene–gene mutual information distributions. Two distinctive issues distinguish the present approach from the studies we have previously made on gene regulatory networks [1, 4, 6]. First of all, no mutual information threshold has been established, so that even interactions with very small statistical dependencies are considered. Second, no ‘binning’ has been applied to the data. By removing these common constraints – at the cost of a highly increased computational burden in already complex calculations– in our analysis, we are in a position to establish (as we will see later on) that the loss of long-range regulation in cancer that we have been reporting is not the product of thresholding nor of binning. We have already shown in our previous work [4] that this phenomenon is not a sampling-size or network inference algorithm artefact either. Since the concept of physical chromosomal distance is only reasonable in the context of genes within the same chromosome, we will, from now on, consider chromosome-wise GRPs G^k[I(i, j)]; here, k ( {1, ²,..., 22, x, y} is an ind ex working as the chromosome label.

By analyzing G^k[I(i, j)], it is possible to develop a deeper understanding of the way correlation structure relates to functional features. For instance, in a previous work [6], we have observed the phenomenon of decay of long-range correlations. The visual inspection of the mutual information distribution dependence in the tumour plots in that work suggests a power-law decay of correlations.

Visual inspection, however, may be misleading us to attribute power-law behaviour to other long-tailed data distributions [2], and even common regression techniques may prove deceptive to establish functionality in heterogeneous variance settings [7, 25]. For these reasons, we decided to implement a comprehensive approach to model discrimination analysis.

To do this, we modelled the chromosomal gene–gene distance d (i, j) dependence of G^k[I(i, j)] by using the R implementation of Generalized Additive Models for Location Scale and Shape (GAMLSSs). GAMLSSs [16] are univariate distributional regression models. Under GAMLSSs, all the parameters of the assumed distribution for the response can be modelled as additive functions of the explanatory variables.

In the present context, the use of GAMLSSs allows us to perform differential goodness-of-fit analyses of the chromosome-wise gene–gene mutual information distributions. The data were adjusted to the following models: power law (PL), log-normal (LN), Weibull (WB), and linear exponential (LE), using polynomials up to order = 4. This analysis generated 4·4·23 = 368 models for the three high-likelihood distributions, under the constant, linear, logarithmic, and polynomial approximations for the 23 chromosomes – for the sake of this study, we have not considered chromosome y since the samples come from chromosomally female populations (xx), as the data come from breast cancer patients.

The model discrimination analysis we performed indicates that the best goodness of fit (by resorting to a combination of extensive coefficient of determination [R2] calculations and Bayesian information content predictors [16] – see Supplementary materials 2) corresponds to either the power law or the Weibull distributions (see Figure 2). The indicators for these distributions come so close that it is not unreasonable to consider that the G^k[I(i, j)] data in basal breast cancer indeed represents an intermediate regime with a mixture of power law and Weibull distribution features. Indeed, it can be proved, by means of the Glivenko–Cantelli theorem [23] that the two distributions belong to the same Glivenko–Cantelli classes of functions [24] and then asymptotically converge as you get more data. As already discussed, the healthy breast distributions are almost distance independent, as can be seen in Supplementary materials 2.

Power-law decay of spatial correlations is a well-established phenomenon in the physical and biological settings [17,22]. Power-law decay may be resulting from the action of a generalized central limit theorem [5] for multiplicative growth processes [13,15,18,20], which may [26,27] or may not [21] result from self-organization processes.

The Weibull distribution has been extensively used in recent times to characterize a variety of skew-distributed phenomena in the physical, biological, and economical sciences [9, 11, 19]. As in the case of the power law, the Weibull distribution may arise as a consequence of generalized central limit theorems, but in this case, for branching – instead of multiplicative growth – processes, as has been proved by Jo et al. [9] from the asymptotics of the Galton–Watson branching process of simple replicative systems. The authors have also shown that the branching process can be mapped into a process of aggregation of clusters. Weibull distributions may also arise in the context of percolation theory. Sornette [20] has shown how the existence of intermittency in continuum percolation changes the distribution from extreme exponential to a smoother Weibull-like form.

By considering the constructive processes associated with fat-tailed probability distributions – in particular, power law and Weibull-like decay – we can hypothesize that the behaviour of gene–gene correlations in breast cancer tumours may arise due to a combination of multiplicative growth, branching, and clustered aggregation processes. These hints are indeed relevant to arrive at the design of experimental strategies to probe what are the actual biological and physical processes behind the dramatic changes in the gene regulation patterns in cancer.