Incorporating biological prior knowledge for Bayesian learning via maximal knowledge-driven information priors

Abstract

Background

Phenotypic classification is problematic because small samples are ubiquitous; and, for these, use of prior knowledge is critical. If knowledge concerning the feature-label distribution – for instance, genetic pathways – is available, then it can be used in learning. Optimal Bayesian classification provides optimal classification under model uncertainty. It differs from classical Bayesian methods in which a classification model is assumed and prior distributions are placed on model parameters. With optimal Bayesian classification, uncertainty is treated directly on the feature-label distribution, which assures full utilization of prior knowledge and is guaranteed to outperform classical methods.

Results

The salient problem confronting optimal Bayesian classification is prior construction. In this paper, we propose a new prior construction methodology based on a general framework of constraints in the form of conditional probability statements. We call this prior the maximal knowledge-driven information prior (MKDIP). The new constraint framework is more flexible than our previous methods as it naturally handles the potential inconsistency in archived regulatory relationships and conditioning can be augmented by other knowledge, such as population statistics. We also extend the application of prior construction to a multinomial mixture model when labels are unknown, which often occurs in practice. The performance of the proposed methods is examined on two important pathway families, the mammalian cell-cycle and a set of p53-related pathways, and also on a publicly available gene expression dataset of non-small cell lung cancer when combined with the existing prior knowledge on relevant signaling pathways.

Conclusion

The new proposed general prior construction framework extends the prior construction methodology to a more flexible framework that results in better inference when proper prior knowledge exists. Moreover, the extension of optimal Bayesian classification to multinomial mixtures where data sets are both small and unlabeled, enables superior classifier design using small, unstructured data sets. We have demonstrated the effectiveness of our approach using pathway information and available knowledge of gene regulating functions; however, the underlying theory can be applied to a wide variety of knowledge types, and other applications when there are small samples.

Background

Small samples are commonplace in phenotypic classification and, for these, prior knowledge is critical [1, 2]. If knowledge concerning the feature-label distribution is available, say, genetic pathways, then it can be used to design an optimal Bayesian classifier (OBC) for which uncertainty is treated directly on the feature-label distribution. As typical with Bayesian methods, the salient obstacle confronting OBC is prior construction. In this paper, we propose a new prior construction framework to incorporate gene regulatory knowledge via general types of constraints in the form of probability statements quantifying the probabilities of gene up- and down-regulation conditioned on the regulatory status of other genes. We extend the application of prior construction to a multinomial mixture model when labels are unknown, a key issue confronting the use of data arising from unplanned experiments in practice.

Regarding prior construction, E. T. Jaynes has remarked [3], “…there must exist a general formal theory of determination of priors by logical analysis of prior information – and that to develop it is today the top priority research problem of Bayesian theory”. It is precisely this kind of formal structure that is presented in this paper. The formal structure involves a constrained optimization in which the constraints incorporate existing scientific knowledge augmented by slackness variables. The constraints tighten the prior distribution in accordance with prior knowledge, while at the same time avoiding inadvertent over restriction of the prior, an important consideration with small samples.

Subsequent to the introduction of Jeffreys’ non-informative prior [4], there was a series of information-theoretic and statistical methods: Maximal data information priors (MDIP) [5], non-informative priors for integers [6], entropic priors [7], reference (non-informative) priors obtained through maximization of the missing information [8], and least-informative priors [9] (see also [10–12] and the references therein). The principle of maximum entropy can be seen as a method of constructing least-informative priors [13, 14], though it was first introduced in statistical mechanics for assigning probabilities. Except in the Jeffreys’ prior, almost all the methods are based on optimization: max- or min-imizing an objective function, usually an information theoretic one. The least-informative prior in [9] is found among a restricted set of distributions, where the feasible region is a set of convex combinations of certain types of distributions. In [15], several non-informative and informative priors for different problems are found. All of these methods emphasize the separation of prior knowledge and observed sample data.

Although the methods above are appropriate tools for generating prior probabilities, they are quite general methodologies without targeting any specific type of prior information. In that regard, the problem of prior selection, in any Bayesian paradigm, is usually treated conventionally (even “subjectively”) and independent of the real available prior knowledge and sample data.

Figure 1 shows a schematic view of the proposed mechanism for Bayesian operator design.

Fig. 1

A schematic illustration of the proposed Bayesian prior construction approach for a binary-classification problem. Information contained in the biological signaling pathways and their corresponding regulating functions is transformed to prior probabilities by MKDIP. Previously observed sample points (labeled or unlabeled) are used along with the constructed priors to design a Bayesian classifier to classify a new sample point (patient)

The a priori knowledge in the form of graphical models (e.g., Markov random fields) has been widely utilized in covariance matrix estimation in Gaussian graphical models. In these studies, using a given graphical model illustrating the interactions between variables, different problems have been addressed: e.g., constraints on the matrix structure [16, 17] or known independencies between variables [18, 19]. Nonetheless, these studies rely on a fundamental assumption: the given prior knowledge is complete and hence provides one single solution. However, in many applications including genomics, the given prior knowledge is uncertain, incomplete, and may be inconsistent. Therefore, instead of interpreting the prior knowledge as a single solution, e.g., a single deterministic covariance matrix, we aim at constructing a prior distribution on an uncertainty class.

In a different approach to prior knowledge, gene-gene relationships (pathway-based or protein-protein interaction (PPI) networks) are used to improve classification accuracy [20–26], consistency of biomarker discovery [27, 28], accuracy of identifying differentially expressed genes and regulatory target genes of a transcription factor [29–31], and targeted therapeutic strategies [32, 33]. The majority of these studies utilize gene expressions corresponding to sub-networks in PPI networks, for instance: mean or median of gene expression values in gene ontology network modules [20], probabilistic inference of pathway activity [24], and producing candidate sub-networks via a Markov clustering algorithm applied to high quality PPI networks [26, 34]. None of these methods incorporate the regulating mechanisms (activating or suppressing) into classification or feature-selection to the best of our knowledge.

The fundamental difference of the work presented in this paper is that we develop machinery to transform knowledge contained in biological signaling pathways to prior probabilities. We propose a general framework capable of incorporating any source of prior information by extending our previous prior construction methods [35–37]. We call the final prior distribution constructed via this framework, a maximal knowledge-driven information prior (MKDIP). The new MKDIP construction constitutes two steps: (1) Pairwise and functional information quantification: information in the biological pathways is quantified by an information theoretic formulation. (2) Objective-based Prior Selection: combining sample data and prior knowledge, we build an objective function, in which the expected mean log-likelihood is regularized by the quantified information in step 1. As a special case, where we do not have any sample data, or there is only one data point available for constructing the prior probability, the proposed framework is reduced to a regularized extension of the maximum entropy principle (MaxEnt) [38].

Owing to population heterogeneity we often face a mixture model, for example, when considering tumor sample heterogeneity where the assignment of a sample to any subtype or stage is not necessarily given. Thus, we derive the MKDIP construction and OBC for a mixture model. In this paper, we assume that data are categorical, e.g. binary or ternary gene-expression representations. Such categorical representations have many potential applications, including those wherein we only have access to a coarse set of measurements, e.g. epifluorescent imaging [39], rather than fine-resolution measurements such as microarray or RNA-Seq data. Finally, we emphasize that, in our framework, no single model is selected; instead, we consider all possible models as the uncertainty class that can be representative of the available prior information and assign probabilities to each model via the constructed prior.

Methods

Notation

Boldface lower case letters represent column vectors. Occasionally, concatenation of several vectors is also shown by boldface lower case letters. For a vector a, a ₀ represents the summation of all the elements and a _i denotes its i-th element. Probability sample spaces are shown by calligraphic uppercase letters. Uppercase letters are for sets and random variables (vectors). Probability measure over the random variable (vector) X is denoted by P(X), whether it be a probability density function or a probability mass function. E _X [f(X)] represents the expectation of f(X) with respect to X. P(x|y) denotes the conditional probability P(X=x|Y=y). θ represents generic parameters of a probability measure, for instance P(X|Y;θ) (or P _θ (X|Y)) is the conditional probability parameterized by θ. γ represents generic hyperparameter vectors. π(θ;γ) is the probability measure over the parameters θ governed by hyperparameters γ, the parameters themselves governing another probability measure over some random variables. Throughout the paper, the terms “pathway” and “network” are used interchangeably. Also, the terms “feature”’ and “variable” are used interchangeably. $\mathcal {M}ult(\boldsymbol {p};n)$ and $\mathcal {D}(\boldsymbol {\alpha })$ represent a multinomial distribution with vector parameter p and n samples, and a Dirichlet distribution with vector α, respectively.

Review of optimal Bayesian classification

Binary classification involves a feature vector X=(X ₁,X ₂,…,X _d )^T∈ℜ^d composed of random variables (features), a binary random variable (label) Y and a classifier ψ(X) to predict Y. The error is ε[ψ]=P(ψ(X)≠Y). An optimal classifier, ψ _bay, called a Bayes classifier, has minimal error, called the Bayes error, among all possible classifiers. The underlying probability model for classification is the joint feature-label distribution. It determines the class prior probabilities c ₀=c=P(Y=0) and c ₁=1−c=P(Y=1), and the class-conditional densities f ₀(x)=P(x|Y=0) and f ₁(x)=P(x|Y=1). A Bayes classifier is given by

$$ \psi_{\text{bay}}(\mathbf{x})\,=\, \left\{\begin{array}{ll} 1, & c_{1}f_{1\!}(\mathbf{x})\geq c_{0}f_{0}(\mathbf{x}), \\ 0, & \text{otherwise.} \end{array}\right. {.} $$

(1)

If the feature-label distribution is unknown but belongs to an uncertainty class of feature-label distributions parameterized by the vector θ∈Θ, then, given a random sample S _n , an optimal Bayeisan classifier (OBC) minimizes the expected error over Θ:

$$ \psi_{\text{OBC}}=\arg \min_{\psi \in \mathcal{C}}E_{\pi^{\ast }(\theta)}[\varepsilon_{\boldsymbol{\theta}}[\psi ]], $$

(2)

where the expectation is relative to the posterior distribution π ^∗(θ) over Θ, which is derived from the prior distribution π(θ) using Bayes’ rule [40, 41]. If we let θ ₀ and θ ₁ denote the class 0 and class 1 parameters, then we can write θ as θ=[c,θ ₀ ,θ ₁ ]. If we assume that c,θ ₀,θ ₁ are independent prior to observing the data, i.e. π(θ)=π(c)π(θ ₀)π(θ ₁), then the independence is preserved in the posterior distribution π ^∗(θ)=π ^∗(c)π ^∗(θ ₀)π ^∗(θ ₁) and the posteriors are given by $\pi ^{\ast }(\boldsymbol {\theta }_{y})\propto \pi (\boldsymbol {\theta } _{y})\prod _{i=1}^{n_{y}}f_{\boldsymbol {\theta }_{y}}(\mathbf {x}_{i}^{y}|y)$ for y=0,1, where $f_{\boldsymbol {\theta }_{y}}(\mathbf {x}_{i}^{y}|y)$ and n _y are the class-conditional density and number of sample points for class y, respectively [42].

Given a classifier ψ _n designed from random sample S _n , from the perspective of mean-square error, the best error estimate minimizes the MSE between its true error (a function of θ and ψ _n ) and an error estimate (a function of S _n and ψ _n ). This Bayesian minimum-mean-square-error (MMSE) estimate is given by the expected true error, $\widehat {\varepsilon }(\psi _{n},S_{n})=\mathrm {E}_{\boldsymbol {\theta } }[\varepsilon (\psi _{n},\boldsymbol {\theta })|S_{n}]$, where ε(ψ _n ,θ) is the error of ψ _n on the feature-label distribution parameterized by θ and the expectation is taken relative to the prior distribution π(θ) [42]. The expectation given the sample is over the posterior probability. Thus, $\widehat {\varepsilon }(\psi _{n},S_{n})=\mathrm {E}_{\pi ^{\ast }}[\varepsilon ]$.

The effective class-conditional density for class y is defined by

$$ f_{\boldsymbol{\Theta}}\left(\mathbf{x}|y\right) =\int_{\boldsymbol{\Theta}_{y}}f_{\boldsymbol{\theta} _{y}}\left(\mathbf{x}|y\right) \pi^{\ast }\left(\boldsymbol{\theta}_{y}\right) d\boldsymbol{\theta}_{y}, $$

(3)

Θ _y being the space for θ _y , and an OBC is given pointwise by [40]

$$ \begin{aligned} &\psi_{\text{OBC}}\left(\mathbf{x}\right) = \\ &\left\{ \begin{array}{ll} 0 & \text{if }\ \mathrm{E}_{\pi^{\ast }}[c]f_{\boldsymbol{\Theta}}\left(\mathbf{x} |0\right) \geq (1-\mathrm{E}_{\pi^{\ast }}[c])f_{\boldsymbol{\Theta}}\left(\mathbf{x} |1\right), \\ 1 & \text{otherwise}. \end{array} \right.. \end{aligned} $$

(4)

For discrete classification there is no loss in generality in assuming a single feature X taking values in the set {1,…,b} of “bins”. Classification is determined by the class 0 prior probability c and the class-conditional probability mass functions p _i =P(X=i|Y=0) and q _i =P(X=i|Y=1), for i=1,…,b. With uncertainty, we assume beta class priors and define the parameters θ ₀={p ₁,p ₂,…,p _b−1} and θ ₁={q ₁,q ₂,…,q _b−1}. The bin probabilities must be valid. Thus, {p ₁,p ₂,…,p _b−1}∈Θ ₀ if and only if 0≤p _i ≤1 for i=1,…,b−1 and $\sum _{i=1}^{b-1}p_{i} \leq 1$, in which case, $p_{b}=1-\sum _{i=1}^{b-1}p_{i}$. We use the Dirichlet priors

$$ \pi (\boldsymbol{\theta}_{0})\propto \prod_{i=1}^{b}p_{i}^{\alpha_{i}^{0}-1}\mathrm{\ and\ }\pi (\theta_{1})\propto \prod_{i=1}^{b}q_{i}^{\alpha_{i}^{1}-1}, \mathrm{\ } $$

(5)

where $\alpha _{i}^{y}>0$. These are conjugate priors, leading to the posteriors of the same form. The effective class-conditional densities are

$$ f_{\boldsymbol{\Theta}}\left(j|y\right) =\frac{U_{j}^{y}+\alpha_{j}^{y}}{n_{y}+\sum_{i=1}^{b} \alpha_{i}^{y}}, $$

(6)

for y=0,1, and the OBC is given by

$$ {\begin{aligned} \psi_{\text{OBC}}(\,j)=\left\{ \begin{array}{ll} 0,\quad\quad~\text{if }\ E_{\pi^{\ast }}[c]f_{\boldsymbol{\Theta}}\left(j|0\right)\geq (1-E_{\pi^{\ast }}[c])f_{\boldsymbol{\Theta}}\left(j|1\right); \\ 1,\quad\quad~\text{otherwise.} \end{array}\right. \end{aligned}} $$

(7)

where $U_{j}^{y}$ denotes the observed count for class y in bin j [40]. Hereafter, $\sum _{i=1}^{b}\alpha _{i}^{y}$ is represented by $\alpha _{0}^{y}$, i.e. $\alpha _{0}^{y} = \sum _{i=1}^{b}\alpha _{i}^{y}$, and is called the precision factor. In the sequel, the sub(super)-script relating to dependency on class y may be dropped; nonetheless, availability of prior knowledge for both classes is assumed.

Multinomial mixture model

In practice, data may not be labeled, due to potential tumor-tissue sample or stage heterogeneity, but still we want to classify a new sample point. A mixture model is a natural model for this scenario, assuming each sample point x _i arises from a mixture of multinomial distributions:

$$ P_{\boldsymbol{\theta}}(\mathbf{x}_{i}) = \sum_{j=0}^{M-1}c_{j} P_{\boldsymbol{\theta}_{j}}(\mathbf{x}_{i}), $$

(8)

where M is the number of components. When there exists two components, similar to binary classification, M=2. The conjugate prior distribution family for component probabilities (if unknown) is the Dirichlet distribution. In the mixture model, no closed-form analytical posterior distribution for the parameters exists, but Markov chain Monte Carlo (MCMC) methods [43] can be employed to numerically calculate the posterior distributions. Since the conditional distributions can be calculated analytically in the multinomial mixture model, Gibbs sampling [44, 45] can be employed for the Bayesian inference. If the prior probability distribution over the component probability vector (c=[c ₀,c ₁,…,c _M ]) is a Dirichlet distribution $\mathcal {D}(\boldsymbol {\phi })$ with parameter vector ϕ, the component-conditional probabilities are $\boldsymbol {\theta }_{j}=[ p^{j}_{1},p^{j}_{2},\ldots,p^{j}_{b}]$, and the prior probability distribution over them is Dirichlet $\mathcal {D}(\boldsymbol {\alpha ^{j}})$ with parameter vector α ^j (as in the classification problem), for j=1,…,M, the Gibbs updates are

$$\begin{aligned} & y_{i}^{(t)} \sim P\left(y_{i}=j|\boldsymbol{c}^{(t-1)},\boldsymbol{\theta}^{(t-1)},\mathbf{x}_{i}\right) \propto c_{j}^{(t-1)} {p_{\mathbf{x}_{i}}^{j,(t-1)}}\\ & \boldsymbol{c}^{(t)} \sim P\left(\boldsymbol{c}|\boldsymbol{\phi},\boldsymbol{y}^{(t)}\right) = \mathcal{D}\left(\boldsymbol{\phi}+ {\sum\nolimits}_{i=1}^{n}\left[I_{y_{i}^{(t)}=1},\ldots,I_{y_{i}^{(t)}=M}\right]\right)\\ & \boldsymbol{\theta_{j}}^{(t)} \sim P\left(\boldsymbol{\theta}_{j}|\mathbf{x},\boldsymbol{y}^{(t)},\boldsymbol{\alpha}_{j}\right)\\ &\quad= \mathcal{D}\left(\boldsymbol{\alpha}_{j}+ {\sum\nolimits}_{i=1:y_{i}^{(t)}=j}^{n}\left[I_{\mathbf{x}_{i}=1},\ldots,I_{\mathbf{x}_{i}=b}\right]\right), \end{aligned} $$

where the super-script in parentheses denotes the chain iteration number, I _w is one if w is true, and otherwise I _w is zero. In this framework, if the inference chain runs for Is iterations, then the numerical approximation of the OBC classification rule is

$$ \psi_{\text{OBC}}(k)\approx \arg \!\max_{y\in { \{1,\dots,M\}}}\sum_{t=1}^{Is}c_{y}^{(t)}p^{y,(t)}_{k}. $$

(9)

Without loss of generality the summation above can be over the iterations of the chain considering burn-in and thinning.

Prior construction: general framework

In this section, we propose a general framework for prior construction. We begin with introducing a knowledge-driven prior probability:

Definition 1

(Maximal Knowledge-driven Information Prior) If Π is a family of proper priors, then a maximal knowledge-driven information prior (MKDIP) is a solution to the following optimization problem:

$$ \arg\min_{\pi \in \Pi}E_{\pi}[C_{\boldsymbol{\theta}}(\xi,D)], $$

(10)

where C _θ (ξ,D)is a cost function that depends on (1) θ: the random vector parameterizing the underlying probability distribution, (2) ξ: state of (prior) knowledge, and (3) D: partial observation (part of the sample data).

Alternatively, by parameterizing the prior probability as π(θ;γ), with γ∈Γ denoting the hyperparameters, an MKDIP can be found by solving

$$ \arg\min_{\boldsymbol{\gamma}\in \Gamma}E_{\pi(\boldsymbol{\theta};\boldsymbol{\gamma})}[C_{\boldsymbol{\theta}}(\xi,D,\boldsymbol{\gamma})]. $$

(11)

In contrast to non-informative priors, the MKDIP incorporates available prior knowledge and even part of the data to construct an informative prior.

The MKDIP definition is very general because we want a general framework for prior construction. The next definition specializes it to cost functions of a specific form in a constrained optimization.

Definition 2

(MKDIP: Constrained Optimization with Additive Costs) As a special case in which C _θ can be decomposed into additive terms, the cost function is of the form:

$$C_{\boldsymbol{\theta}}(\xi,D,\boldsymbol{\gamma}) =(1-\beta)g^{(1)}_{\boldsymbol{\theta}}(\xi,\boldsymbol{\gamma}) + \beta g^{(2)}_{\boldsymbol{\theta}}(\xi,D), $$

where β is a non-negative regularization parameter. In this case, the MKDIP construction with additive costs and constraints involves solving the following optimization problem:

$$ \begin{aligned} &\arg\min_{\boldsymbol{\gamma}\in\Gamma}E_{\pi(\boldsymbol{\theta};\boldsymbol{\gamma})} \Big[(1-\beta)g^{(1)}_{\boldsymbol{\theta}}(\xi,\boldsymbol{\gamma}) + \beta g^{(2)}_{\boldsymbol{\theta}}(\xi,D)\Big]\\ & \text{Subject to:} \quad E_{\pi(\boldsymbol{\theta};\boldsymbol{\gamma})}[g^{(3)}_{\boldsymbol{\theta},i}(\xi)]=0;~i\in\{1,\ldots,n_{c}\}, \end{aligned} $$

(12)

where $g^{(3)}_{\boldsymbol {\theta },i}$, ∀i∈{1,…,n _c }, are constraints resulting from the state of knowledge ξvia a mapping:

$$\mathcal{T}: \xi \rightarrow E_{\pi\left(\boldsymbol{\theta};\boldsymbol{\gamma}\right)}\left[g^{(3)}_{\boldsymbol{\theta},i}(\xi)\right], \forall i \in \{1,\dots,n_{c}\}. $$

In the sequel, we will refer to g ⁽¹⁾(·) and g ⁽²⁾(·) as the cost functions, and $g_{i}^{(3)}(\cdot)$’s as the knowledge-driven constraints. We begin with introducing information-theoretic cost functions, and then we propose a general set of mapping rules, denoted by $\mathcal {T}$ in Definition 2, to convert biological pathway knowledge into mathematical forms. We then consider special cases with information-theoretic cost functions.

Information-theoretic cost functions

Instead of having least squares (or mean-squared error) as the standard cost functions in classical statistical inference problems, there is no universal cost function in the prior construction literature. That being said, in this paper, we utilize several widely used cost functions in the field:

• 1. (Maximum Entropy) The principle of maximum-entropy (MaxEnt) for probability construction [38] leads to the least informative prior given the constraints in order to prevent adding spurious information. Under our general framework MaxEnt can be formulated by setting:

  $$\beta=0,~g_{\boldsymbol{\theta}}^{(1)} = -H[\boldsymbol{\theta}], $$

  where H[.] denotes the Shannon entropy.

• 2. (Maximal Data Information) The maximal data information prior (MDIP) introduced by Zellner [46] as a choice of the objective function is a criterion for the constructed probability distribution to remain maximally committed to the data [47]. To achieve MDIP, we can set our general framework with:

  $$\begin{aligned} \beta=0,~g_{\boldsymbol{\theta}}^{(1)} &= \ln \pi(\boldsymbol{\theta};\boldsymbol{\gamma}) + H[P(x|\boldsymbol{\theta})]\\ &= \ln \pi(\boldsymbol{\theta};\boldsymbol{\gamma})-E_{x|\boldsymbol{\theta}}[\ln P(x|\boldsymbol{\theta})]. \end{aligned} $$

• 3. (Expected Mean Log-likelihood) The cost function introduced in [35] is the first one that utilizes part of the observed data for prior construction. In that, we have

  $$\beta =1,~g^{(2)}_{\boldsymbol{\theta}}=-\ell(\boldsymbol{\theta};D), $$

  where $\ell (\boldsymbol {\theta };D)=\frac {1}{n_{D}}\sum _{i=1}^{n_{D}}\log f(\boldsymbol {x}_{i}|\boldsymbol {\theta })$ is the mean log-likelihood function of the sample points used for prior construction (D), and n _D denotes the number of sample points in D. In [35], it is shown that this cost function is equivalent to the average Kullback-Leibler distance between the unknown distribution (empirically estimated by some part of the samples) and the uncertainty class of distributions.

As originally proposed, the preceding approaches did not involve expectation over the uncertainty class. They were extended to the general prior construction form in Definition 1, including the expectation, in [36] to produce the regularized maximum entropy prior (RMEP), the regularized maximal data information prior (RMDIP), and the regularized expected mean log-likelihood prior (REMLP). In all cases, optimization was subject to specialized constraints.

For MKDIP, we employ the same information-theoretic cost functions in the prior construction optimization framework. MKDIP-E, MKDIP-D, and MKDIP-R correspond to using the same cost functions as REMP, RMDIP, and REMLP, respectively, but with the new general types of constraints. To wit, we employ functional information from the signaling pathways and show that by adding these new constraints that can be readily derived from prior knowledge, we can improve both supervised (classification problem with labelled data) and unsupervised (mixture problem without labels) learning of Bayesian operators.

From prior knowledge to mathematical constraints

In this part, we present a general formulation for mapping the existing knowledge into a set of constraints. In most scientific problems, the prior knowledge is in the form of conditional probabilities. In the following, we consider a hypothetical gene network and show how each component in a given network can be converted into the corresponding inequalities as general constraints in MKDIP optimization.

Before proceeding we would like to say something about contextual effects on regulation. Because a regulatory model is not independent of cellular activity outside the model, complete control relations such as A→B in the model, meaning that gene B is up-regulated if and only if gene A is up-regulated (after some time delay), do not necessarily translate into conditional probability statements of the form P(X _B =1|X _A =1)=1, where X _A and X _B represent the binary gene values corresponding to genes A and B, respectively. Rather, what may be observed is P(X _B =1|X _A =1)=1−δ, where δ>0. The pathway A→B need not imply P(X _B =1|X _A =1)=1 because A→B is conditioned on the context of the cell, where by context we mean the overall state of the cell, not simply the activity being modeled. δ is called a conditioning parameter. In an analogous fashion, rather than P(X _B =1|X _A =0)=0, what may be observed is P(X _B =1|X _A =0)=η, where η>0, because there may be regulatory relations outside the model that up-regulate B. Such activity is referred to as cross-talk and η is called a crosstalk parameter. Conditioning and cross-talk effects can involve multiple genes and can be characterized analytically via context-dependent conditional probabilities [48].

Consider binary gene values X ₁,X ₂,…,X _m corresponding to genes g ₁,g ₂,…,g _m . There are m2^m−1 conditional probabilities of the form

$$\begin{array}{@{}rcl@{}} &&P(X_{i}=k_{i}| X_{1}=k_{1},\dots,X_{i-1}=k_{i-1},X_{i+1}= \\ &&\quad k_{i+1},\dots,X_{m}=k_{m})\\ &&=a^{k_{i}}_{i}(k_{1},\dots,k_{i-1},k_{i+1},\dots,k_{m}) \end{array} $$

(13)

to serve as constraints, the chosen constraints to be the conditional probabilities whose values are known (approximately). For instance, if g ₂ and g ₃ regulate g ₁, with X ₁=1 when X ₂=1 and X ₃=0, then, ignoring context effects,

$$a^{1}_{1}(1,0,k_{4},\dots,k_{m})=1 $$

for all k ₄,…,k _m . If, however, we take context conditioning into effect, then

$$a_{1}^{1}(1,0,k_{4},\dots,k_{m}) = 1-\delta_{1}(1,0,k_{4},\dots,k_{m}), $$

where δ ₁(1,0,k ₄,…,k _m ) is a conditioning parameter.

Moreover, ignoring context effects,

$$\begin{array}{@{}rcl@{}} a^{1}_{1}(1,1,k_{4},\dots,k_{m})&=& a^{1}_{1}(0,0,k_{4},\dots,k_{m})\\ &=&a^{1}_{1}(0,1,k_{4},\dots,k_{m}) = 0 \end{array} $$

for all k ₄,…,k _m . If, however, we take crosstalk into effect, then

$$\begin{array}{@{}rcl@{}} a^{1}_{1}(1, 1, k_{4},\dots, k_{m})& =& \eta_{1}(1, 1, k_{4},\dots, k_{m})\\ a^{1}_{1}(0, 0, k_{4},\dots, k_{m}) &=& \eta_{1}(0, 0, k_{4},\dots, k_{m})\\ a^{1}_{1}(0, 1, k_{4},\dots, k_{m}) &=& \eta_{1}(0, 1, k_{4},\dots, k_{m}), \end{array} $$

where η ₁(1,1,k ₄,…,k _m ), η ₁(0,0,k ₄,…,k _m ), and η ₁(0,0,k ₄,…,k _m ) are crosstalk parameters. In practice it is unlikely that we would know the conditioning and crosstalk parameters for all combinations of k ₄,…,k _m ; rather, we might just know the average, in which case, δ ₁(1,0,k ₄,…,k _m ) reduces to δ ₁(1,0), η ₁(1,1,k ₄,…,k _m ) reduces to η ₁(1,1), etc.

In this paradigm, the constraints resulting from our state of knowledge are of the following form:

$$ \begin{aligned} &g^{(3)}_{\boldsymbol{\theta},i}(\xi)=\\ &P(X_{i}=k_{i}| X_{1}=k_{1},\dots,X_{i-1}=k_{i-1},X_{i+1}=k_{i+1}, \\ & \quad \dots,X_{m}=k_{m}) -a^{k_{i}}_{i}(k_{1},\dots,k_{i-1},k_{i+1},\dots,k_{m}). \end{aligned} $$

(14)

The basic setting is very general and the conditional probabilities are what they are, whether or not they can be expressed in the regulatory form of conditioning or crosstalk parameters. The general scheme includes previous constraints and approaches proposed in [35] and [36] for the Gaussian and discrete setups. Moreover, in those we can drop the regulatory-set entropy because it is replaced by the set of conditional probabilities based on the regulatory set, whether forward (master predicting slaves) or backwards (slaves predicting masters) [48].

In this paradigm, the optimization constraints take the form

$$\begin{array}{@{}rcl@{}} && a^{k_{i}}_{i}(k_{1},\dots, k_{i-1}, k_{i+1},\dots, k_{m})- \\ && \quad \varepsilon_{i}(k_{1},\dots, k_{i-1}, k_{i+1},\dots, k_{m}) \\ && \leq E_{\pi(\boldsymbol{\theta};\boldsymbol{\gamma})}[P(X_{i} = k_{i}| X_{1} = k_{1},\dots, X_{i-1} = k_{i-1}, \\ && \quad \quad \quad \quad X_{i+1} = k_{i+1},\dots, X_{m} = k_{m})] \\ && \leq a^{k_{i}}_{i}(k_{1},\dots, k_{i-1}, k_{i+1},\dots, k_{m}) + \\ &&\quad \quad \varepsilon_{i}(k_{1},\dots, k_{i-1}, k_{i+1},\dots, k_{m}), \end{array} $$

(15)

where the expectation is with respect to the uncertainty in the model parameters, that is, the distribution of the model parameter θ, and ε _i is a slackness variable. Not all will be used, depending on our prior knowledge. In fact, the general conditional probabilities will not likely be used because they will likely not be known when there are too many conditioning variables. For instance, we may not know the probability in Eq. (13), but may know the conditioning on part of the variables which can be extracted from some interaction network (e.g. biological pathways). A slackness variable can be considered for each constraint to make the constraint framework more flexible, thereby allowing potential error or uncertainty in prior knowledge (allowing potential inconsistencies in prior knowledge). When using slackness variables, these variables also become optimization parameters, and a linear function (summation of all slackness variables) times a regulatory coefficient is added to the cost function of the optimization in Eq. (12). In other words, when having slackness variables, the optimization in Eq. (12) can be written as

$$ \begin{aligned} & \arg\min_{\boldsymbol{\gamma}\in\Gamma, \boldsymbol{\varepsilon}\in\mathcal{E}}E_{\pi(\boldsymbol{\theta};\boldsymbol{\gamma})}\Big[\lambda_{1}[(1-\beta)g^{(1)}_{\boldsymbol{\theta}}(\xi,\boldsymbol{\gamma}) + \beta g^{(2)}_{\boldsymbol{\theta}}(\xi,D)] \\ & \quad \quad \quad \quad \quad \quad \quad \quad \quad + \lambda_{2}\sum_{i=1}^{n_{c}}\varepsilon_{i}\Big]\\ & \text{Subject to:} -\varepsilon_{i}\leq E_{\pi(\boldsymbol{\theta};\boldsymbol{\gamma})}[g^{(3)}_{\boldsymbol{\theta},i}(\xi)]\leq\varepsilon_{i};~i\in\{1,\ldots,n_{c}\}, \end{aligned} $$

(16)

where λ ₁ and λ ₂ are non-negative regularization parameters, and ε and $\mathcal {E}$ represent the vector of all slackness variables and the feasible region for slackness variables, respectively. For each slackness variable, a possible range can be defined (note that all slackness variables are non-negative). The higher the uncertainty is about a constraint stemming from prior knowledge, the greater the possible range for the corresponding slackness variable can be (more on this in the “Results and discussion” section).

The new general type of constraints discussed here introduces a formal procedure for incorporating prior knowledge. It allows the incorporation of knowledge of the functional regulations in the signaling pathways, any constraints on the conditional probabilities, and also knowledge of the cross-talk and conditioning parameters (if present), unlike the previous work in [36] where only partial information contained in the edges of the pathways is used in an ad hoc way.

An illustrative example and connection with conditional entropy

Now, consider a hypothetical network depicted in Fig. 2. For instance, suppose we know that the expression of gene g ₁ is regulated by g ₂, g ₃, and g ₅. Then we have

$$P(X_{1} = 1| X_{2} = k_{2}, X_{3} = k_{3}, X_{5} = k_{5}) = a_{1}^{1}(k_{2}, k_{3}, k_{5}).$$

Fig. 2

An illustrative example showing the components directly connected to gene 1. In the Boolean function {AND, OR, NOT}={∧,∨,−}. Based on the regulating function of gene 1, it is up-regulated if gene 5 is up-regulated and genes 2 and 3 are down-regulated

As an example,

$$P(X_{1} = 1| X_{2} = 1, X_{3} = 1, X_{5} = 0) = a_{1}^{1}(1_{2}, 1_{3}, 0_{5}), $$

where the notation 1₂ denotes 1 for the second gene. Further, we might not know a ₁(k ₂,k ₃,k ₅) for all combinations of k ₂, k ₃, k ₅. Then we use the ones that we know. In the case of conditioning with g ₂, g ₃, and g ₅ regulating g ₁, with g ₁ on if the others are on,

$$a^{1}_{1}(1_{2}, 1_{3}, 1_{5}) = 1 -\delta_{1}(1_{2}, 1_{3}, 1_{5}). $$

If limiting to 3-gene predictors, g ₃, and g ₅ regulate g ₁, with g ₁ on if the other two are on, then

$$a^{1}_{1}(k_{2}, 1_{3}, 1_{5}) = 1 - \delta_{1}(k_{2}, 1_{3}, 1_{5}), $$

meaning that the conditioning parameter depends on whether X ₂=0 or 1.

Now, considering the conditional entropy, assuming that $\delta _{1} = \max _{(k_{2},k_{3},k_{5})} \delta _{1}(k_{2},k_{3},k_{5})$ and δ ₁<0.5, we may write

$$ \begin{aligned} &H[X_{1}|X_{2},X_{3},X_{5}] =\\ &-\left[\sum_{\mathcal{X}_{2},\mathcal{X}_{3},\mathcal{X}_{5}}\left[ P\left(X_{1}=0|X_{2}=x_{2},X_{3}=x_{3},X_{5}=x_{5}\right)\right.\right.\\ &\quad \quad\times P\left(X_{2}=x_{2},X_{3}=x_{3},X_{5}=x_{5}\right)\\ &\quad \quad\log\left[P\left(X_{1}=0|X_{2}=x_{2},X_{3}=x_{3},X_{5}=x_{5}\right)\right]\\ &\quad \quad +P\left(X_{1}=1|X_{2}=x_{2},X_{3}=x_{3},X_{5}=x_{5}\right)\\ &\quad \quad\times P\left(X_{2}=x_{2},X_{3}=x_{3},X_{5}=x_{5}\right)\\ &\quad \quad\left.\left.\log\left[P\left(X_{1}=1|X_{2}=x_{2}, X_{3}=x_{3},X_{5}=x_{5}\right)\right]\right]{\vphantom{\sum_{\mathcal{X}_{2}}}}\right]\\ &\quad\quad\quad\leq h(\delta_{1}), \end{aligned} $$

where h(δ)=−[δ log(δ)+(1−δ) log(1−δ)]. Hence, bounding the conditional probabilities, the conditional entropy is in turn bounded by h(δ ₁):

$${\lim}_{\delta_{1}\rightarrow 0^{+}} H\left[X_{1}|X_{2},X_{3},X_{5}\right] = 0. $$

It should be noted that constraining H[X ₁|X ₂,X ₃,X ₅] would not necessarily constrain the conditional probabilities, and may be considered as a more relaxed type of constraints. But, for example, in cases where there is no knowledge about the status of a gene given its regulator genes, constraining entropy is the only possible approach.

In our illustrative example, if we assume that the Boolean regulating function of X ₁ is known as shown in Fig. 2 and context effects exist, then the following knowledge constraints can be extracted from the pathway and regulating function:

$$\begin{aligned} &a^{0}_{1}\left(k_{2}, k_{3}, 0_{5}\right) = 1 - \delta_{1}\left(k_{2}, k_{3}, 0_{5}\right)\\ &a^{0}_{1}\left(k_{2}, 1_{3}, k_{5}\right) = 1 - \delta_{1}\left(k_{2}, 1_{3}, k_{5}\right)\\ &a^{0}_{1}\left(1_{2}, k_{3}, k_{5}\right) = 1 - \delta_{1}\left(1_{2}, k_{3}, k_{5}\right)\\ &a^{1}_{1}\left(0_{2}, 0_{3}, 1_{5}\right) = 1 - \delta_{1}\left(0_{2}, 0_{3}, 1_{5}\right). \end{aligned} $$

Now if we assume that the context does not affect the value of X ₁, i.e. the value of X ₁ can be fully determined by knowing the values of X ₂, X ₃, and X ₅, then we have the following equations:

$$\begin{array}{*{20}l} & a^{0}_{1}\left(k_{2}, k_{3}, 0_{5}\right) = P\left(X_{1}=0|X_{5}=0\right) = 1 \end{array} $$

(17a)

$$\begin{array}{*{20}l} & a^{0}_{1}\left(k_{2}, 1_{3}, k_{5}\right) = P\left(X_{1}=0|X_{3}=1\right) = 1 \end{array} $$

(17b)

$$\begin{array}{*{20}l} & a^{0}_{1}\left(1_{2}, k_{3}, k_{5}\right) = P\left(X_{1}=0|X_{2}=1\right)= 1 \end{array} $$

(17c)

$$\begin{array}{*{20}l} & a^{1}_{1}\left(0_{2}, 0_{3}, 1_{5}\right)= P(X_{1}=1|X_{2}=0,X_{3}=0, \\ & \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad X_{5}=1) = 1. \end{array} $$

(17d)

It can be seen from the equations above that for some setups of the regulator values, only a subset of them determines the value of X ₁, regardless of the other regulator values. If we assume that the value of X ₅ cannot be observed, for example X ₅ is an extracellular signal that cannot be measured in gene expression data and thereafter X ₅ is not in the features of our data, the only constraints relevant to the feature-label distribution that can be extracted from the regulating function knowledge will be

$$ \begin{aligned} &a^{0}_{1}\left(k_{2}, 1_{3}, k_{5}\right) = P\left(X_{1}=0|X_{3}=1\right) = 1\\ & a^{0}_{1}\left(1_{2}, k_{3}, k_{5}\right) = P\left(X_{1}=0|X_{2}=1\right)= 1. \end{aligned} $$

(18)

Special case of Dirichlet distribution

Fixing the value of a single gene, being ON or OFF (i.e. X _i =0 or X _i =1, respectively), corresponds to a partition of the states, $\mathcal {X }=\{1,\dots,b\}$. Here, the portions of $\mathcal {X}$ for which (X _i =k ₁,X _j =k ₂) and (X _i ≠k ₁,X _j =k ₂), for any k ₁,k ₂∈{0,1}, are denoted by $\mathcal {X}^{i,j}(k_{1},k_{2})$ and $\mathcal {X}^{i,j}(k_{1}^{c},k_{2})$, respectively. For the Dirichlet distribution, where θ=p and γ=α, the constraints on the expectation over the conditional probability in (15) can be explicitly written as functions of the prior probability parameters (hyperparameters). For the parameter of the Dirichlet distribution, a vector α indexed by $\mathcal {X}$, we denote the variable indicating the summation of its entities in $\mathcal { X}^{i,j}(k_{1},k_{2})$ by $\overline {\alpha }^{i,j}(k_{1},k_{2})=\sum _{k \in \mathcal {X}^{i,j}(k_{1},k_{2})}\alpha _{k}$. The notation can be easily extended for the cases having more than two fixed genes. In this setup, if the set of random variables corresponding to genes other than g _i and the vector of their corresponding values are shown by $\tilde {X}_{i}$ and $\tilde {X}_{i}$, respectively, the expectation over the conditional probability in (15) is [36]:

$$ \hspace{12pt} \begin{aligned} &E_{\boldsymbol{p}}\left[P\left(X_{i} = k_{i}| X_{1} = k_{1},\dots, X_{i-1} = k_{i-1},\right.\right.\\ &\quad\left. \left.X_{i+1} = k_{i+1},\dots, X_{m} = k_{m}\right)\right]\\ & =E_{\boldsymbol{p}}\left[\frac{\sum_{k\in \mathcal{X}^{i,\tilde{X}_{i}}\left(k_{i},\tilde{x}_{i}\right)}p_{k}}{\sum_{k\in \mathcal{X}^{i,\tilde{X}_{i}}\left(k_{i},\tilde{x}_{i}\right)}p_{k}+\sum_{k\in \mathcal{X}^{i,\tilde{X}_{i}}\left(k_{i}^{c},\tilde{x}_{i}\right)}p_{k}}\right] \\ &=\frac{\overline{\alpha }^{i,\tilde{X}_{i}}\left(k_{i},\tilde{x}_{i}\right)}{\overline{\alpha }^{i,\tilde{X}_{i}}\left(k_{i},\tilde{x}_{i}\right)+\overline{\alpha }^{i,\tilde{X}_{i}}\left(k_{i}^{c},\tilde{x}_{i}\right)}, \end{aligned} $$

(19)

where the summation in the numerator and the first summation in the denominator of the second equality is over the states (bins) for which ($X_{i} = k_{i}, \tilde {X}_{i} = \tilde {x}_{i}$), and the second summation in the denominator is over the states (bins) for which ($X_{i} = k_{i}^{c}, \tilde {X}_{i} = \tilde {x}_{i}$).

If there exists a set of genes that completely determines the value of gene g _i (or only a specific setup of their values that determines the value, as we had in our illustrative example in Eq. (17)), then the constraints on the conditional probability conditioned on all the genes other than g _i can be changed to be conditioned on that set only. Specifically, let R _i denote the set of random variables corresponding to such a set of genes/proteins and suppose there exists a specific setup of their values r _i that completely determines the value of gene g _i . If the set of all random variables corresponding to the genes/proteins other than X _i and R _i is denoted by $\boldsymbol {B}_{i}=\tilde {X}_{(i,\boldsymbol {R}_{i})}$, and their corresponding values by b _i , then the constraints on the conditional probability can be written as

$$ \begin{aligned} &E_{\boldsymbol{p}}\left[P\left(X_{i} = k_{i}| \boldsymbol{R}_{i} = \boldsymbol{r}_{i}\right)\right]\\ &=E_{\boldsymbol{p}}\left[\frac{\sum_{\boldsymbol{b}_{i}\in O_{\boldsymbol{B}_{i}}}\sum_{k\in \mathcal{X}^{i,\boldsymbol{R}_{i},\boldsymbol{B}_{i}}\left(k_{i},\boldsymbol{r}_{i},\boldsymbol{b}_{i}\right)}p_{k}}{\sum_{\boldsymbol{b}_{i}\in O_{\boldsymbol{B}_{i}}}\sum_{k\in \mathcal{X}^{i,\boldsymbol{R}_{i},\boldsymbol{B}_{i}}(k_{i},\boldsymbol{r}_{i},\boldsymbol{b}_{i})}p_{k}}\right.\\ & \quad \quad \quad\left.\frac{}{+ \sum_{\boldsymbol{b}_{i}\in O_{\boldsymbol{B}_{i}}}\sum_{k\in \mathcal{X}^{i,\boldsymbol{R}_{i},\boldsymbol{B}_{i}}\left(k_{i}^{c},\boldsymbol{r}_{i},\boldsymbol{b}_{i}\right)}p_{k}}\right] \\ &= \frac{\sum_{\boldsymbol{b}_{i}\in O_{\boldsymbol{B}_{i}}}\overline{\alpha }^{i,\boldsymbol{R}_{i},\boldsymbol{B}_{i}}(k_{i},\boldsymbol{r}_{i},\boldsymbol{b}_{i})}{\sum_{\boldsymbol{b}_{i}\in O_{\boldsymbol{B}_{i}}}\overline{\alpha }^{i,\boldsymbol{R}_{i},\boldsymbol{B}_{i}}(k_{i},\boldsymbol{r}_{i},\boldsymbol{b}_{i})} \\ & \quad \quad \frac{}{+\sum_{\boldsymbol{b}_{i}\in O_{\boldsymbol{B}_{i}}}\overline{\alpha }^{i,\boldsymbol{R}_{i},\boldsymbol{B}_{i}}(k_{i}^{c},\boldsymbol{r}_{i},\boldsymbol{b}_{i})}, \end{aligned} $$

(20)

where $\phantom {\dot {i}\!}O_{\boldsymbol {B}_{i}}$ is the set of all possible vectors of values for B _i .

For a multinomial model with a Dirichlet prior distribution, a constraint on the conditional probabilities translates into a constraint on the above expectation over the conditional probabilities (as in Eq. (15)). In our illustrative example and from the equations in Eq. (17), there are four of these constraints on the conditional probability for gene g ₁. For example, in the second constraint from the second line of Eq. (17) (Eq. 17b), X _i =X ₁, k _i =0, R _i ={X ₃}, r _i =[0], and B _i ={X ₂,X ₅}. One might have several constraints for each gene extracted from its regulatory function (more on extracting general constraints from regulating functions in the “Results and discussion” section).

Results and discussion

The performance of the proposed general prior construction framework with different types of objective functions and constraints is examined and compared with other methods on two pathways, a mammalian cell-cycle pathway and a pathway involving the gene TP53. Here we employ Boolean network modeling of genes/proteins (hereafter referred to as entities or nodes) [49] with perturbation (BNp). A Boolean Network with p nodes (genes/proteins) is defined as B=(V,F), where V represents the set of entities (genes/proteins) {v ₁,…,v _p }, and F is the set of Boolean predictor functions {f ₁,…,f _p }. At each step in a BNp, a decision is made by a Bernoulli random variable with the success probability equal to the perturbation probability, p _pert , as to whether a node value is determined by perturbation of randomly flipping its value or by the logic model imposed from the interactions in the signaling pathways. A BNp with a positive perturbation probability can be modeled by an ergodic Markov chain, and possesses a steady-state distribution (SSD) [50]. The performance of different prior construction methods can be compared based on the expected true error of the optimal Bayesian classifiers designed with those priors, and also by comparing these errors with some other well known classification methods. Another comparison metric of prior construction methods is the expected norm of the difference between the true parameters and the posterior mean of these parameters inferred using the constructed prior distributions. Here, the true parameters are the vectors of the true class-conditional SSDs, i.e. the vectors of the true class-conditional bin probabilities of the BNp.

Moreover, the performance of the proposed framework is compared with other methods on a publicly available gene expression dataset of non-small cell lung cancer when combined with the existing prior knowledge on relevant signaling pathways.

Mammalian cell cycle classification

A Boolean logic regulatory network for the dynamical behavior of the cell cycle of normal mammalian cells is proposed in [51]. Figure 3(a) shows the corresponding pathways. In normal cells, cell division is coordinated via extracellular signals controlling the activation of CycD. Rb is a tumor suppressor gene and is expressed when the inhibitor cyclins are not present. Expression of p27 blocks the action of CycE or CycA, and lets the tumor-suppressor gene Rb be expressed even in the presence of CycE and CycA, and results in a stop in the cell cycle. Therefore, in the wild-type cell-cycle network, expressing p27 lets the cell cycle stop. But following the proposed mutation in [51], for the mutated case, p27 is always inactive (i.e. can never be activated), thereby creating a situation where both CycD and Rb might be inactive and the cell can cycle in the absence of any growth factor.

Fig. 3

Signaling pathways corresponding to Tables 1 and 2. Signaling pathways for: 3(a) the normal mammalian cell cycle (corresponding to Table 1) and 3(b) a simplified pathway involving TP53 (corresponding to Table 2)

Table 1 Boolean regulating functions of normal mammalian cell cycle [51]. In the Boolean functions {AND, OR, NOT}={∧,∨,−}

Table 2 Boolean regulating functions corresponding to the pathway in Fig. 3(b) [54]. In the Boolean functions {AND, OR, NOT}={∧,∨,−}

The full functional regulations in the cell-cycle Boolean network are shown in Table 1.

Following [36], for the binary classification problem, y=0 corresponds to the normal system functioning based on Table 1, and y=1 corresponds to the mutated (cancerous) system where CycD, p27, and Rb are permanently down-regulated (are stuck at zero), which creates a situation where the cell cycles even in the absence of any growth factor. The perturbation probability is set to 0.01 and 0.05 for the normal and mutated system, respectively. A BNp has a transition probability matrix (TPM), and as mentioned earlier, with positive perturbation probability can be modeled by an ergodic Markov chain, and possesses a SSD [50]. Here, each class has a vector of steady-state bin probabilities, resulting from the regulating functions of its corresponding BNp and the perturbation probability. The constructed SSDs are further marginalized to a subset of seven genes to prevent trivial classification scenarios. The final feature vector is x=[E2F,CycE,CycA,Cdc20,Cdh1,UbcH10,CycB], and the state space size is 2⁷=128. The true parameters for each class are the final class-conditional steady-state bin probabilities, p ⁰ and p ¹ for the normal and mutated systems, respectively, which are utilized for taking samples.

Classification problem corresponding to TP53

TP53 is a tumor suppressor gene involved in various biological pathways [36]. Mutated p53 has been observed in almost half of the common human cancers [52], and in more than 90% of patients with severe ovarian cancer [53]. A simplified pathway involving TP53, based on logic in [54], is shown in Fig. 3(b). DNA double-strand break affects the operation of these pathways, and the Boolean network modeling of these pathways under this uncertainty has been studied in [53, 54]. The full functional regulations are shown in Table 2.

Following [36], two scenarios, dna-dsb=0 and dna-dsb=1, weighted by 0.95 and 0.05, are considered and the SSD of the normal system is constructed based on the ergodic Markov chain model of the BNp with the regulating functions in Table 2 by assuming the perturbation probability 0.01. The SSD for the mutated (cancerous) case is constructed by assuming a permanent down regulation of TP53 in the BNp, and perturbation probability 0.05. Knowing that dna-dsb is not measurable, and to avoid trivial classification situations, the SSDs are marginalized to a subset of three entities x=[ATM,Wip1,Mdm2]. The state space size in this case is 2³=8. The true parameters for each class are the final class-conditional steady-state bin probabilities, p ⁰ and p ¹ for the normal and mutated systems, respectively, which are used for data generation.

Extracting general constraints from regulating functions

If knowledge of the regulating functions exists, it can be used in the general constraint framework of the MKDIP, i.e. it can be used to constrain the conditional probabilities. In other words, the knowledge about the regulating function of gene i can be used to set ε _i (k ₁,…,k _i−1,k _i+1,…,k _m ), and $a^{k_{i}}_{i}(k_{1},\dots, k_{i-1}, k_{i+1},\dots, k_{m})$ in the general form of constraints in (15). If the true regulating function of gene i is known, and it is not context sensitive, then the conditional probability of its status, $a^{k_{i}}_{i}(k_{1},\dots, k_{i-1}, k_{i+1},\dots, k_{m})$, is known for sure, and δ _i (k ₁,…,k _i−1,k _i+1,…,k _m )=0. But in reality, the true regulating functions are not known, and are also context sensitive. The dependence on the context translates into δ _i (k ₁,…,k _i−1,k _i+1,…,k _m ) being greater than zero. The greater the context effect on the gene status, the larger δ _i is. Moreover, the uncertainty over the regulating function is captured by the slackness variables ε _i (k ₁,…,k _i−1,k _i+1,…,k _m ) in Eq. (15). In other words, the uncertainty is translated to the possible range of the slackness variable values in the prior construction optimization framework. The higher the uncertainty is, the greater the range should be in the optimization framework. In fact, slackness variables make the whole constraint framework consistent, even for cases where the conditional probability constraints imposed by prior knowledge are not completely in line with each other, and guarantee the existence of a solution.

As an example, for the classification problems of the mammalian cell-cycle network and the TP53 network, assuming the regulating functions in Tables 1 and 2 are the true regulating functions, the context effect can be observed in the dependence of the output of the Boolean regulating functions in the tables on the extracellular signals, non-measurable entities, and the genes that have been marginalized out in our setup. In the absence of quantitative knowledge about the context effect, i.e. $a^{k_{i}}_{i}(k_{1},\dots, k_{i-1}, k_{i+1},\dots, k_{m})$ for all possible setups of the regulator values, one can impose only those with such knowledge. For example, in the mammalian cell-cycle network, CycB’s regulating function only depends on the values included in the observed feature set; therefore the conditional probabilities are known for all regulator value setups. But for CycE the regulating function depends on Rb, which is marginalized out in our feature set, and also itself depends on an extracellular signal. Hence, the conditional probability constraints for CycE are known only for the setup of the features that determine the output of the Boolean regulating function independent of the other regulator values.

In our comparison analysis, $a^{k_{i}}_{i}(k_{1},\allowbreak \dots, k_{i-1},\allowbreak k_{i+1}, \allowbreak \dots,\allowbreak k_{m})$ for each gene/protein in Eq. (15) is set to one for the feature value setups that determine the Boolean regulating output regardless of the context. But since the observed data are not fully described by these functions, and the system has uncertainty, we let the possible range for the slackness variables in Eq. (15) be [0,1).

We now continue the examples on two of the mammalian cell-cycle network nodes, CycB and CycE. For CycB the following constraints on the prior distribution are extracted from its regulating function:

$$\begin{aligned} & E_{\boldsymbol{p}}[P(v_{10}=0|v_{8}=1)]\geq 1-\epsilon_{1} \\ & E_{\boldsymbol{p}}[P(v_{10}=0|v_{7}=1)]\geq 1-\epsilon_{2} \\ & E_{\boldsymbol{p}}[P(v_{10}=1|v_{7}=0,v_{8}=0)]\geq 1-\epsilon_{3}. \end{aligned} $$

For CycE, one of its regulators is Rb (v ₂), which is not included in the feature set, i.e. not observed, but is known to be down-regulated in the mutated (cancerous) case. Thus, the set of constraints extracted from the regulating function of CycE for the normal case includes only

$$\begin{aligned} & E_{\boldsymbol{p}}[P(v_{4}=0|v_{3}=0)]\geq 1-\epsilon_{1} \end{aligned} $$

and for the mutated case consists of

$$\begin{aligned} & E_{\boldsymbol{p}}[P(v_{4}=0|v_{3}=0)]\geq 1-\epsilon_{1} \\ & E_{\boldsymbol{p}}[P(v_{4}=1|v_{3}=1)]\geq 1-\epsilon_{2}. \end{aligned} $$

As another example, for the TP53 network, the set of constraints extracted from the regulating functions in Table 2 for the normal case are shown in the left panel of Table 3.

Table 3 The set of constraints extracted from the regulating functions and pathways for the TP53 network. Constraints extracted from the Boolean regulating functions in Table 2 corresponding to the pathway in Fig. 3(b) used in MKDIP-E, MKDIP-D, MKDIP-R (left). Constraints extracted based on [36] from the pathway in Fig. 3(b) used in RMEP, RMDIP, REMLP (right)

The first and second constraints for MKDIP in the left panel of Table 3 come from the regulating function of v ₂ in Table 2. Although v ₁ is an extracellular signal, the value of v ₄ imposes two constraints on the value of v ₂. But the regulating function of v ₄ in Table 2 only depends on v ₃, which is not included in our feature set, so we have no imposed constraints on the conditional probability from its regulating function. The other two constraints for MKDIP in the left panel of Table 3 are extracted from the regulating function of v ₅ in Table 2. Although v ₃ is not included in the observed features, for two setups of its regulators, (v ₂=1) and (v ₂=0,v ₄=1), the value of v ₅ can be determined, so the constraint is imposed on the prior distribution from the regulating function. For comparison, the constraints extracted from the pathway in Fig. 3(b) based on the method of [36] are provided in the right panel of Table 3.

Performance comparison in classification setup

For both the mammalian cell cycle and TP53 problems, the performance of 11 methods are compared for classification performance. OBC with the Jeffreys’ prior, OBC with our previous prior construction methods in [36] (RMEP, RMDIP, REMLP), OBC with our proposed general framework of constraints (MKDIP-E, MKDIP-D, MKDIP-R), and also well known methods including Histogram rule (Hist), CART [55], Random Forest (RF)[56], and Support Vector Machine classification (SVM) [57, 58]. Also, for all the Bayesian methods using OBC, the posterior mean of the parameters’ distance from the true parameters is calculated and compared. The samples from the true distributions are stratified fixing two different class prior probabilities. Following [36], we assume that $\max _{i}p_{i}^{y,true},~{for} ~y\in \{0,1\}$, is known within a +/−5% interval (can come from existing population statistics in practice). Two simulation scenarios are performed: one assuming the complete knowledge of the optimal precision factors [36] $\alpha _{0}^{y}=\sum _{i=1}^{b}\alpha _{i}^{y}, y\in \{0,1\}$ for prior construction methods (oracle precision factor); and the other estimating the optimal precision factor from the observed data itself. Two class prior probabilities, c=0.6 and c=0.5, are considered. Along with the true class-conditional SSDs of the two classes, the corresponding Bayes error corresponds to the best performance that any classification rule for that classification problem (feature-label distribution) can yield. Fixing c and the true class-conditional bin probabilities, n sample points by stratified sampling (n ₀=⌈c n⌉ sample points from class 0 and n ₁=n−n ₀ sample points from class 1) are taken for prior construction (if used by the method), classifier training, and posterior distribution calculations. Then the designed classifier’s true classification error is calculated for all classification methods. The posterior mean of parameter distance from the true parameter (true steady-state bin probabilities vector) is calculated based on $\sum _{y=0}^{1}||\boldsymbol {\alpha }^{y\ast }/\alpha _{0}^{y\ast }-\boldsymbol {p}^{y}||^{2}$, where α ^y∗ and p ^y represent the parameters of the posterior distribution and true bin probabilities vector for class y, respectively. For each fixed c and n, 800 Monte Carlo repetitions are done to calculate the expected classification errors and posterior distances from the true parameters for each parameter setup. For REMLP and MKDIP-R, which use a fraction of data in their prior construction procedure, 10 data points from each class are used for prior construction, and all for the inference and posterior calculation (here the number of data points used for prior construction is not fine-tuned, but a small number is chosen to avoid overfitting). The overall procedure taken for a fixed classification problem and a fixed sample size (fixed n) in each Monte Carlo repetition is as follows:

• The true bin probabilities p ⁰ and p ¹ are fixed.

• n ₀ and n ₁ are determined using c as n ₀=⌈c n⌉ and n−n ₀.

• Observations (training data) are randomly sampled from the multinomial distribution for each class, i.e. $(U^{y}_{1},\ldots,U^{y}_{b})\sim \mathcal {M}ult(\boldsymbol {p}^{y};n_{y})$, for y∈{0,1}.

• 10 data points are randomly taken from the training data points of each class to be used in the prior construction methods that utilize partial data (REMLP and MKDIP-R)

• All the classification rules are trained based on their constructed prior (if applicable to that classification rule) and the training data.

• The classification errors associated with the classifiers are computed using p ⁰ and p ¹. Also for the Bayesian methods, the posterior probability mass (mean) distance from the true parameters (true bin probabilities, p ⁰ and p ¹) is calculated.

The regularization parameter λ ₁ is set to 0.5, and λ ₂ is set to 0.25 and 0.5 for the mammalian cell cycle classification problem and the TP53 classification problem, respectively. The results of expected classification error and posterior mean distance from the true parameters for the mammalian cell-cycle network are shown in Tables 4 and 6, respectively. Tables 5 and 7 contain the results of expected classification error and posterior mean distance from the true parameters for the TP53 network.

Table 4 Expected true error of different classification rules for the mammalian cell-cycle network. The constructed priors are considered using two precision factors: optimal precision factor (left) and estimated precision factor (right), with c=0.5, and c=0.6, where the minimum achievable error (Bayes error) is denoted by E r r _Bayes

Table 5 Expected true error of different classification rules for the TP53 network. The constructed priors are considered using two precision factors: optimal precision factor (left) and estimated precision factor (right), with c=0.5, and c=0.6, where the minimum achievable error (Bayes error) is denoted by E r r _Bayes

The best performance (with the lowest error in Tables 4 and 5, and lowest distance in Tables 6 and 7) for each sample size, are written in bold. For the mammalian cell-cycle network, MKDIP methods show the best (or as good as the best) performance in all the scenarios in terms of both the expected classification error and posterior parameter estimates. For the TP53 network, MKDIP methods show the best performances in posterior parameter estimates, and are competitive with the previous knowledge-driven prior construction methods in terms of the expected classification error.

Table 6 Expected difference between the true model (for mammalian cell-cycle network) and estimated posterior probability masses. Optimal precision factor (left) and estimated precision factor (right), with c=0.5, and c=0.6

Table 7 Expected difference between the true model (for TP53 network) and estimated posterior probability masses. Optimal precision factor (left) and estimated precision factor (right), with c=0.5, and c=0.6

Performance comparison in mixture setup

The performance of the OBC with different prior construction methods, including OBC with the Jeffreys’ prior, OBC with prior constructions methods of [36] (RMEP, RMDIP, REMLP), and OBC with the general framework of constraints (MKDIP-E, MKDIP-D, MKDIP-R), are further compared in the mixture setup with missing labels, for both the mammalian cell-cycle and the TP53 systems. Also, the OBC with prior distribution centered on the true parameters with a relatively low variance (hereinafter abbreviated as PDCOTP method in Tables 8 and 9) is considered as the comparison baseline, though it is not a practical method. Similar to the classification problems, we assume that only two components (classes) exist, normal and mutated (cancerous). Here, c ₀ is fixed at 0.6 (c ₁=1−c ₀=0.4), but the sampling is not stratified. The component-conditional SSDs (bin probabilities) for the two components are as before in the classification problem, i.e. the same as the class-conditional SSDs in the classification problem.

Table 8 Expected errors of different Bayesian classification rules in the mixture model for the mammalian cell-cycle network. Expected true error (left) and expected error on unlabeled training data (right), with c ₀=0.6

Table 9 Expected errors of different Bayesian classification rules in the mixture model for the TP53 network. Expected true error (left) and expected error on unlabeled training data (right), with c ₀=0.6

For each sample point, first the label (y) is generated from a Bernoulli distribution with success probability c ₁, and then the bin observation is generated given the label, from the corresponding class-conditional SSD (class conditional bin probabilities vector, p ^y), i.e. the bin observation is a sample from a categorical distribution with parameter vector p ^y but the label is hidden for the inference chain and classifier training. n sample points are generated and fed into the Gibbs inference chain with different priors from the different prior construction methods. Then the OBC is calculated based on Eq. 9. For each sample size, 400 Monte Carlo repetitions are done to calculate the expected true error and the error of classifying the unlabeled observed data used for the inference itself.

To have a fair comparison of different methods’ class-conditional prior probability construction, we assume that we have a rough idea of the mixture weights (class probabilities). In practice this can come from existing population statistics. That is, the Dirichlet prior distribution over the mixture weights (class probabilities) parameters, ϕ in $\mathcal {D}(\boldsymbol {\phi })$, are sampled in each iteration from a uniform distribution that is centered on the true mixture weights vector +/−10% interval, and fixed for all the methods in that repetition. For the REMLP and MKDIP-R that need labeled data in their prior construction procedure, the predicted labels from using the Jeffreys’ prior are used and one fourth of the data points are used in prior construction for these two methods, and all for inference. The reason for using a larger number of data points in prior construction within the mixture setup compared to the classification setup is that in the mixture setup, data points are missing their true class labels, and the initial label estimates may be inaccurate. One can use a relatively larger number of data points in prior construction, which still avoids overfitting. The regularization parameters λ ₁ and λ ₂ are set as in the classification problem. Optimal precision factors are used for all prior construction methods. The results are shown in Tables 8 and 9 for the mammalian cell-cycle and TP53 models, respectively. The best performance (lowest error) for each sample size and the best performance among practical methods (all other than PDCOTP), if different, is written in bold. As can be seen from the tables, in most cases the MKDIP methods have the best performance among the practical methods. With larger sample sizes, MKDIP-R even outperforms PDCOTP in the mammalian cell-cycle system.

Performance comparison on a real data set

In this section the performance of the proposed methods are examined on a publicly available gene expression dataset. Here, we have considered the classification of two subtypes of non-small cell lung cancer (NSCLC), lung adenocarcinoma (LUA) versus lung squamous cell carcinoma (LUS). Lung cancer is the second most commonly diagnosed cancer and the leading cause of cancer death in both men and women in the United States [59]. About 84% of lung cancers are NSCLC [59] and LUA and LUS combined account for about 70% of lung cancers based on the American Cancer Society statistics for NSCLC. We have downloaded LUA and LUS datasets (both labeled as TCGA provisional) in the form of mRNA expression z-scores (based on RNA-Seq profiling) from the public database cBioPortal [60, 61] for the patient sets tagged as “All Complete Tumors", denoting the set of all tumor samples that have mRNA and sequencing data. The two datasets for LUA and LUS consist of 230 and 177 sample points, respectively. We have quantized the data into binary levels based on the following preprocessing steps. First, to remove the bias for each patient, each patient’s data are normalized by the mean of the z-scores of a randomly selected subset from the list of the recurrently mutated genes (half the size of the list) from the MutSig [62] (directly provided by cBioPortal). Then, a two component Gaussian mixture model is fit to each gene in each data set, and the normalized data are quantized by being assigned to one component, namely 0 or 1 (1 being the component with higher mean). We confine the feature set to {EGFR,PIK3CA,AKT,KRAS,RAF1,BAD,P53,BCL2} which are among the genes in the most relevant signaling pathways to the NSCLC [63]. These genes are altered, in different forms, in 86% and 89% of the sequenced LUA and LUS tumor samples on the cBioPortal, respectively. There are 256 bins in this classification setting, since the feature set consists of 8 genes. The pathways relevant to the NSCLC classification problem considered here are collected from KEGG [64, 65] Pathways for NSCLC and PI3K-AKT signaling pathways, and also from [63], as shown in Fig. 4. The corresponding regulating functions are shown in Table 10.

Fig. 4

Signaling pathways corresponding to NSCLC classification. The pathways are collected from KEGG Pathways for NSCLC and PI3K-AKT pathways, and from [63]

Table 10 Regulating functions corresponding to the signaling pathways in Fig. 4. In the Boolean functions {AND, OR, NOT}={∧,∨,−}

The informative prior construction methods utilize the knowledge in the pathways in Fig. 4, and the MKDIP methods also use the regulating relationships in Table 10 in order to construct prior distributions. The incidence rate of the two subtypes, LUA and LUS, varies based on demographic factors. Here, we approximate the class probability c=P(Y=LUA) as c≈0.57, based on the latest statistics of the American Cancer Society for NSCLC, and also based on a weighted average of the rates for 11 countries given in [66]. In each Monte Carlo repetition, n sample points by stratified sampling, i.e. n ₀=⌈c n⌉ and n ₁=n−n ₀ sample points, are randomly taken from preprocessed LUA (class 0) and LUS (class 1) datasets, respectively, for prior construction (if used by the method) and classifier training, and the rest of the sample points are held out for error estimation. For each n, 400 Monte Carlo repetitions are done to calculate the expected classification error. In the prior construction methods, first the optimization is solved for both classes with the precision factors $\alpha ^{y}_{0}=200, y\in \{0,1\}$, and then their optimal values are estimated using the training points. For REMLP and MKDIP-R, which use a fraction of the training data in their prior construction procedure, min(20,max(6,⌊0.25n _y ⌋)) sample points from the training data of each class (y∈{0,1}) are used for prior construction, and all the training data are used for inference. The regularization parameters λ ₁ and λ ₂ are set to 0.5 and 0.25, respectively. The results are shown in Table 11. In the table, the best performance among Hist, CART, RF and SVM is shown as Best Non Bayesian method. Best RM represents the best performance among RMEP, RMDIP, and REMLP. Best MKDIP denotes the best performance among the MKDIP methods.

Table 11 Expected error of different classification rules calculated on a real dataset. The classification is between LUA (class 0) and LUS (class 1), with c=0.57

The best performing rule for each sample size is written in bold. As can be seen from the table, OBC with MKDIP prior construction methods has the best performance among the classification rules. It is also clear that the classification performance can be significantly improved when pathway prior knowledge is integrated for constructing prior probabilities, especially when the sample size is small.

Implementation remarks

The results presented in this paper are based on Monte Carlo simulations, where thousands of optimization problems are solved for each sample size for each problem. Thus, the regularization parameters and the number of sample points used in prior construction are preselected for each problem. One can use cross validation to set these parameters in a specific application. It has been shown in [36] that by assuming precision factors greater than 1 ($\alpha _{0}^{y}>1, y\in \{0,1\}$), all three objective functions used are convex for the class of Dirichlet prior probabilities for multinomial likelihood functions. But unfortunately, we cannot guarantee the convexity of the feasible space due to the convolved constraints. Therefore, we have employed algorithms for nonconvex optimization problems and there is no guarantee of convergence to the global optimum. The method used for solving the optimization framework of the prior construction is based on the interior-point algorithm for nonlinear constrained optimization [67, 68] implemented in the fmincon function in MATLAB. In this paper, since the interest is in classification problems with small training sample sizes (which is often the case in bioinformatics) and also due to Monte Carlo simulations, we have only shown performance results on small networks with only a few genes. In practice, there would be no problem using the proposed method for larger networks, since there would then be a single one-time analysis. One should also note that with small sample sizes, one needs feature selection to keep the number of features small. In the experiments in this paper, feature selection is automatically done by focusing on the most relevant network by biological prior knowledge.

Conclusion

Bayesian methods have shown promising performance in classification problems in the presence of uncertainty and small sample sizes, which often occur in translational genomics problems. The impediment in using these methods is prior construction to integrate existing prior biological knowledge. In this paper we have proposed a knowledge-driven prior construction method with a general framework of mapping prior biological knowledge into a set of constraints. Knowledge can come from biological signaling pathways and other population studies, and be translated into constraints over conditional probabilities. This general scheme includes the previous approaches of using biological prior knowledge in prior construction. Here, the superior performance of this general scheme is shown on two important pathway families, the mammalian cell-cycle pathway and the pathway centering around TP53. In addition, prior construction and the OBC are extended to a mixture model, where data sets are with missing labels. Moreover, comparisons on a publicly available gene expression dataset show that classification performance can be significantly improved for small sample sizes when corresponding pathway prior knowledge is integrated for constructing prior probabilities.