Materials and method

Dataset
A global gene expression (GE) dataset (GSE38376) from 1) cells sensitive to lapatinib (said to be under "parental conditions") and 2) cells with acquired resistance to lapatinib was obtained from Komurov et al. [17]. Expression values were measured using Illumina HumanHT-12 V3.0 expression beadchip (GPL6947). Samples include SKBR3 parental and resistant (SKBR3-R) each under basal conditions and in response to 0.1 μM and 1 μM lapatinib after 24 hours, where the resistant cell line variant (SKBR3-R) showed 100-fold more resistance to lapatinib treatment than the parental SKBR3 cell line, as reported by Komurov et al. [17]. These gene expression datasets used probe-level annotation, which we converted into gene-level annotation. To obtain gene-level GE values, probes were mapped to gene symbols using the corresponding annotation file (GPL6947). While mapping, the average GE values were calculated across all probes if the same gene symbol was annotated to multiple probes. Two GE data matrices were constructed for parental SKBR3 cell lines and resistant SKBR3-R cell lines, respectively, where rows were labelled with gene symbols and columns were labelled with different treatment conditions (0, 0.1 μM and 1 μM of lapatinib).

Construction of a gene-gene relationship network
We define the gene-gene relationship network as GGR:= (S,R) for each GE data matrix. Here, S is a set of 370 cancer related genes collected from the Cancer Gene Census [23]. R is defined as the set of pair-wise relationships among seed genes. A gene pair (genei, genej) is included in R if the corresponding absolute Pearson Correlation Coefficient (PCC) is above some threshold, and defined as a pair-wise relationship. These threshold values were empirically chosen for parental and resistant conditions individually, based on the corresponding distributions of all pairwise absolute PCC values. Note PCC values resulting from probes mapped to the same gene were trivially ignored.

Bayesian statistical modeling of GGR network

Network model
For statistical modeling of networks, exponential families of distributions offer robust and flexible parametric models [24]. These probabilistic models can be used to evaluate the probability that an edge is present in the network. They can also be used to quantify topological properties of networks by summarizing them in a parametric form and associating sufficient statistics with those parameters [19,24]. In this study, we use a special class of exponential family distributions known as ERGM (Exponential Random Graph Models), also known as the p1-model, which was introduced by Holland and Leinhardt [24].
A gene-gene relationship network with g genes can be regarded as a random variable X taking values from a set G containing all 2g(g−1) possible relationship networks [24,25]. Let u be a generic point of G which can alternatively be denoted as the realization of X by X = u. Let the binary outcome uij = 1 if genei interacts with genej, or uij = 0 otherwise. Then u is a binary data matrix [19]. Let Pr(u) be the probability function on G given by (1) Pr(u)=Pr(X=u)=1κθexp∑pθpzpu
where zp(u) is the network statistic of type p, θp is the parameter associated with zp(u) and κ(θ) is the normalizing constant that ensures Pr(u) is a proper probability distribution (sums to 1 over all u in G) [26]. The parameter θ is a vector of model parameters associated with network statistics and needs to be estimated. See [24] for further details.
A major limitation of the p1-model is the difficulty of calculating the normalizing constant, κ(θ), since it is a sum over the entire graph space. Estimating the maximum likelihood of this model becomes intractable as there are 2g(g−1) possible directed graphs (or 2g(g−1)2 undirected graphs), each having g nodes (genes). A technique called maximum pseudolikelihood estimation has been developed to address this problem [27]. This technique employs MCMC methods such as Gibbs or Metropolis-Hastings sampling algorithms [28].
The construction of the p1-model for a directed network is described in an Appendix Additional file 1: Appendix I. For the gene-gene relationship network with undirected edges, the description of the p1-model can be simplified by using only two Bernoulli variables Yij0 and Yij1 instead of four as follows: Yijk=1ifuij=k,0otherwise
The simplified p1-model can then be defined using the following two equations to predict the probability of an edge being present between genei and genej: (2) logPrYij1=1=λij+θ+αi+αj
(3)  log Pr Y  ij 0   = 1     = λ  ij
for i<j. Note that λij is chosen to ensure Pr(Yij0=1)+Pr(Yij1=1)=1. In this formulation, the expansiveness and attractiveness parameters were reduced to a single parameter, α, which represents the propensity of a gene to be connected in an undirected network. Hence, the p1-model seeks to find the probabilities of edge formation in a network considering its structural features explicitly.

Bayesian modeling
We used a fully Bayesian approach for modeling our gene-gene relationship network. Parameter estimation is a crucial step in statistical modeling, for which a classical approach is maximum likelihood estimation (MLE). However, unlike MLE, Bayesian techniques involve calculation of posterior probabilities of model parameters by training the model with given data. We assume that the data  follows the generative model , and assign a prior probability Pθ|ℳ to the parameter vector θ under the model . Then Bayes’ rule for calculating posterior probability is as follows: (4) Prθ|ℳ,D=PrD|θ,ℳ×Prθ|ℳZ
where PrD|θ,ℳ is the likelihood function. Now, the marginal likelihood  can be expressed as (5) Z=PrD|ℳ=∫PrD|ℳ,θ×Pθ|ℳdθ,
Computing the exact solution for the marginal likelihood  is often intractable since it is prone to the curse of dimensionality. Fortunately, Markov Chain Monte Carlo (MCMC) methods such as Gibbs sampling and Metropolis-Hastings methods do not require  to be explicitly computed. In general, MCMC methods are stochastic simulation techniques which generate samples from the joint distribution Pℳ,θ|D for calculating the posterior probabilities of parameters. Here we used Gibbs sampling methods, which sample iteratively, one parameter at a time, from the full conditional distribution given the current and previous values of all other parameters. To implement Gibbs sampling, we employed WinBUGS [29], which is a high-level software package providing an easy interface for implementing complex Bayesian models. In WinBUGS, users are free from background lower-level programming details, and only have to express the model precisely.
We hypothesized that gene-pairs involved in drug resistance are likely to be found with high probabilities in the resistant network but low probabilities in the parental network. Therefore, we built two networks, one from resistant datasets and the other from parental datasets. In this Bayesian approach, the model likelihood is defined in Equations (2) and (3), where Yk is the data matrix calculated from the observed data u. Here we have two Yk data matrices, namely a gene-gene relationship network YkR derived from resistant samples and YkP derived from parental samples.
Our approach is a hierarchical Bayesian model in that model parameters are in turn dependent on hyperparameters. We assign the density parameter θ in Equation (2) a normal prior distribution with mean 0 and standard deviation σθ. (6) θ∼N0,σθ2
Note, in WinBUGS the parameter τ, called the precision, replaces the standard deviation parameter σ of the normal distribution, where, τ=σ−2. For the hyperparameter τθ we specify a gamma prior distribution as follows, since it is a conjugate prior for the normal distribution: (7) τθ∼Gammaa0,b0
We set a0 = 0.001 and b0 = 0.001 to make the prior for θnoninformative, making its standard deviation wide to express large uncertainty [19]. For attractiveness/ expansiveness parameters αi and αj, we followed the approach used by Adams et al. [30]. (8) αiRαiP∼N00,Σ
(9)  Σ  − 1   ∼ Wishart 1  0   0  1      , 2
Here, αiR and αiP represent the expansiveness/attractiveness parameters for the network model of resistant and parental conditions, respectively.

Drug resistant cross-talk prediction
Since, Lapatinib is an EGFR and ErbB inhibitor, we considered the cross-talks between the EGFR/ErbB signaling pathway and other signaling pathways. Here cross-talks can be defined as any gene-pair (genei,genej) in which genei ∈ {genes in EGFR/ErbB signaling pathway} and genej ∈ {genes in other signaling pathways}, or vice versa [31]. Thus if both genes in any gene-pair were found in the same signaling pathway, that particular gene-pair was trivially ignored. For that purpose, we collected 24 signaling pathways from Reactome [32] (downloaded at 19/05/2014), 35 signaling pathways from KEGG [33,34] (downloaded at 21/10/2014), and 63 signaling pathways from WikiPathway [35] (downloaded at 16/10/2014) databases. Each signaling pathway downloaded from these databases was encoded as tab-delimitated lists of gene symbols.
To determine whether a given gene-pair is involved in drug resistance, we calculated a simple odds ratio of the corresponding two posterior probabilities: (10) odds=PrYij1R=1PrYij1P=1
where, Yij1R and Yij1P are gene-gene relationships defined over resistant and parental networks, respectively, and the probabilities are estimated using MCMC sampling. We then selected only those gene-pairs for which the odds score and PruijR=1 are greater than conservative thresholds, and identified these as the gene-pairs which are potentially involved in drug-resistance.