Discovery of co-occurring driver pathways in cancer

Brief introduction of the maximum weight submatrix problem

where Γ(g)={i:A _{i g} =1} denotes the set of patients in which the gene g is mutated and Γ(M)=∪_g∈MΓ(g), |Γ(M)| measures the coverage of M and $\omega \left(M\right)=\sum _{g\in M}\left|\Gamma \right(g\left)\right|-\left|\Gamma \right(M\left)\right|$ measures the coverage overlap of M.

In addition to the stochastic MCMC search procedure [18], a BLP model has been introduced to solve this problem exactly [19]. Inspired by the BLP model, a binary linear programming model CoMDP to discover co-occurring driver pathways was developed here (Figure 1). The focus was placed on finding possible cooperative driver pathways in carcinogenesis. For example, in Figure 1, the gene set D not C can be detected using MCMC or BLP for k=2. This is because gene set D has higher mutation score W than that of gene set C. Mutations in the gene set B and the gene set C occurred simultaneously among a cohort of patients. CoMDP can successfully identify such co-occurring gene sets which may have been missed by previous approaches.

CoMDP: a binary linear programming model for the identification of co-occurring driver pathways

For the mutation matrix A, let us consider two submatrices M and N (which correspond to two gene sets or pathways S and T). Given the coverage Γ(M) and Γ(N) of the two gene sets (sometimes called individual coverage in this study), we define (1) the common coverage $c(M,N)=\left|\Gamma \right(M)\bigcap \Gamma (N\left)\right|$; (2) the union coverage $b(M,N)=\left|\Gamma \right(M)\bigcup \Gamma (N\left)\right|$. We further define the non-shared coverage d(M,N)=b(M,N)−c(M,N), which describes the extent of the mutation co-occurrence between the two gene sets: the smaller the value d, the larger the co-occurrence is. As suggested before, ω(M) and ω(N) reflect the exclusivity of M and N respectively.

where u _j and v _j are indicators whether column j of A falls into the submatrx M or N, so all the columns j’s with u _j =1 and v _j =1 constitute M and N respectively; x _i and y _i are indicators whether the entries of row i of M and N are not all zeros, so $\sum _{i=1}^m{x}_i$ and $\sum _{i=1}^m{y}_i$ represent the coverage of M and N (i.e., |Γ(M)| and |Γ(N)|) respectively; z _i is the indicator whether both x _i and y _i equal to 1, so $\sum _{i=1}^m{z}_i$ represents the overlap between the coverage of M and N (i.e., the common coverage c(M,N)). k is the total number of genes within S and T; and finally, λ and η are two parameters controlling the common coverage c(M,N) and the non-shared coverage d(M,N) of the two gene sets.

Note that $\sum _{i=1}^m{z}_i$ and $\sum _{i=1}^m(x_i+y_i-2z_i)$ in model (3) are always nonnegative according to the constraints in (4). One can properly set λ and η to be positive or negative to obtain gene sets with specific characteristics. For example, if λ<0 and η>0, the model tends to detect gene sets with large non-shared coverage but small common coverage under the maximization of G(x,y,z,u,v). Certainly, λ>0 and η<0 must be set if one wants to identify co-occurring driver pathways by maximizing the function H in (2), which is the main focus of this study. More discussion on the behavior of the model with λ and η can be referred to Simulation study below.

mod_CoMDP: Finding a pathway that co-occurs with a known one

where x _i and u _j are indicators whether the entries of row i in gene set C are not all zeros and whether the gene corresponding to column j of A falls into C, respectively; y _i ,z _i ,v _j and the parameters λ and η have the same meaning as in (3) and (4); r is the size of the desired gene set D.

Generally, a branch-and-bound algorithm or others can be used to produce an optimal exact solution for CoMDP (also for mod_CoMDP). In this study, an IBM ILOG CPLEX Optimizer was used to test the effectiveness of the model. The experiments were performed on a 2.50 GHz Core i5-2520M CPU PC. For each given k, CoMDP can automatically identify two gene sets when the sum of their sizes equals k. Although the problem was NP-hard, it can still be solved efficiently due to the sparsity of the mutation matrix.

Statistical significance

A permutation test was used to assess the significance of the results. As in a previous study [18], the weight W in (1) served as a statistic to test the significance of the exclusivity and coverage of each identified gene set (called individual significance). We employed the co-occurrence ratio, which is defined as the ratio of the common coverage to the union coverage, as the statistic to test the significance of the co-occurrence of these two gene sets (called co-occurrence significance).

Simulation data

Three datasets were constructed to illustrate the properties of the proposed method. The first set of simulated data, Sim _data1, was generated as in a previous study [19]. First, an empty m (samples) ×n (genes) matrix was given (m=500,n=1,000 were used). Then, gene sets M _i (i=1,⋯,I; and each set has 10 genes) with a mutation probability p _i were embedded in the matrix (p _i =1−i·Δ, Δ=0.05, and I=10 were used here). For each sample, a gene uniformly chosen from M _i with p _i was mutated, and once one gene was mutated, the other genes in M _i had a probability p₀ to be mutated (here p₀=0.04 was used). Finally, the genes not in M _i were mutated in at most three samples. The second dataset, Sim _data2, was generated using the strategy described above for noisy probability p₀ from 0.04 to 0.24 in steps of 0.02.

The third dataset, Sim _data3, was generated as follows. Starting with an empty r×s matrix (here r=600,s=1,000 were used), we embedded J gene sets N₁,N₂,⋯,N _J (J≥1, here J=9 was used) into it. N _i has size m _i ×n₀ (in this study n₀=5 and m _i = [ m/2]+2ⁱ⁻¹ were used where [m/2] denotes the integer part of m/2). N _i was constructed according to the strategy like that for M _i stated above. Similarly, we mutated the genes not in N _i at most in three samples.

We note that the average mutation rate for genes in a dataset in current simulation study is comparable to those of real datasets. For example, for Sim _data2 with the noisy probabilities of 0.04 and 0.24, each gene has an average mutation rate 0.0142 and 0.0274 respectively. For the four biological datasets (GBM1, GBM2, lung cancer and ovarian cancer) introduced in the following subsection, each dataset has an average mutation rate of 0.0658, 0.0416, 0.0206, and 0.0134 for each gene, respectively.

Biological data

To assess the proposed methods for practical applications, four biological datasets were collected due to their popularity and abundant prior knowledge. Note that the CoMDP can be easily applied to other cancer mutation datasets.

The glioblastoma multiforme data 1 (GBM1) and the lung adenocarcinoma dataset were obtained directly from a previous study [18]. These sets contain mutations in 178 genes across 84 GBM patients (samples) and 356 genes in 163 lung cancer patients, respectively. The GBM2 and ovarian carcinoma datasets were obtained from another previous study [16]. These sets contain CNAs for 1269 genes spanning 169 GBM patients, somatic mutations for 343 genes across 135 GBM patients, CNAs in 966 genes across 559 ovarian patients, and somatic mutations in 8431 genes across 320 ovarian cancer patients. For the last two datasets, somatic mutations and CNAs were first integrated by merging the genes on the common patients. Finally, a binary mutation matrix A was obtained for each of the four datasets. The genes that are mutated in the same samples were combined into a gene set which was named as a metagene in this study. Note that the definition of a metagene differs from that defined based on the matrix factorization method.

Simulation study

CoMDP can get the optimal solution of the original maximum weight submatrix problem

Like the BLP model [19], CoMDP can detect all the embedded gene sets in Sim _data1 with k=10, η=1 and λ<0 (here λ=−10 was used). In the current situation, CoMDP usually degenerates to find one gene set which corresponds to the optimal solution of BLP. Sometimes it can produce two gene sets with no common coverage. For example, with Sim _data1 where ten gene sets were embedded in the 500 (samples) × 1000 (genes) matrix, we identified two sets with two and eight genes respectively, which are mutated in 74 and 340 samples respectively. But their common coverage is 0. In fact, these two sets constitute one of the embedded gene sets which can be found using the BLP method, so they can be viewed as the same driver gene set as obtained using BLP directly.

Note that as p₀ increases, the exclusivity among the genes in M _i decreases, so the detection of the embedded gene sets M _i becomes more and more difficult. Let k=10. CoMDP (λ=−10, η=1) and BLP were applied on Sim _data2. Both were able to precisely identify all ten embedded gene sets when p₀≤0.10. For BLP, the average number of detected embedded gene sets decreased sharply as p₀ increased. However, by properly choosing the parameter η, CoMDP can obtain more accurate and robust results than BLP at high values of p₀. For example, CoMDP with λ=−10 and η=2 has much higher identification accuracy than BLP for p₀≥0.20 (Figure 2A).

GBM1 dataset

For k=4, two gene sets were detected: {CDKN2A, MG ₁} and {MTAP, CYP27B1} (MG ₁ is a metagene consisting of CDK4, FAM119B, MARCH9, TSFM, CENTG1, METTL1 and TSPAN31) with individual significance p₁=0.0207, p₂=0.0058, co-occurrence significance p_1,2<0.0001, and co-occurrence ratio r_1,2=0.9412 (Table 1). The two genes MTAP and CDKN2A were found to be frequently co-deleted [27, 28]. They are both located on chromosome 9p21, a typical tumor suppressor region whose deletion is related to many different types of cancers. CYP27B1 and the metagene MG ₁ were mutated in the same patients with one exception: a single-nucleotide mutation was recorded in one additional patient for CYP27B1 (Figure 3A). Previous studies have suggested that CDK4 is the target of a common CNA in the corresponding patients [29]. Two protein products of CDKN2A, INK4A (also known as p16) and ARF (also known as p14), are involved in the p53 and RB tumor suppressor pathways (Figure 3B). It has been shown that any error disrupting these pathways causes tumor formation [30]. CDKN2A and CDK4 are considered part of the RB pathway. Both MTAP and CYP27B1 encode important enzymes. The enzymes encoded by MTAP play a major role in polyamine metabolism and those encoded by CYP27B1 play a role in calcium metabolism and tissue differentiation.

k	Gene set 1	Gene set 2	p 1	p 2	n 1	n 2	r 1,2	p 1,2
4	CDKN2A, M G1	MTAP, CYP27B1	0.0207	0.0058	50	49	0.9412	<0.0001
	CDKN2A, TP53,
5	M G 1	CDKN2B, CYP27B1	0.0003	0.0018	68	57	0.7606	<0.0001
	CDKN2A, PTEN,	CDKN2B, TP53,
6	CYP27B1	M G 1	0.0002	0.0003	69	71	0.8182	<0.0001
	CDKN2A, PTEN,	CDKN2B, RB1,
7	CYP27B1	TP53, M G1	0.0002	0.0001	69	74	0.8571	<0.0001
	CDKN2A, PTEN,	CDKN2B, RB1,
8	NF1, CYP27B1	M G1, ERBB2	0.0015	<0.0001	72	70	0.8933	<0.0001
	CDKN2A, PTEN, NF1,	CDKN2B, RB1,
9	CYP27B1, KDR	M G1, ERBB2	0.0006	<0.0001	74	70	0.9200	<0.0001
	CDKN2A, PTEN, CYP27B1,	CDKN2B, NF1, RB1,
10	KDR, M G2	M G1, ERBB2	<0.0001	<0.0001	72	73	0.9333	<0.0001

Here p₁ and p₂ are the p-values of the individual significance of two identified gene sets, p_1,2 represents the p-value of their co-occurrence significance, n₁ and n₂ denote their respective coverage, and r_1,2 is the ratio of the common coverage to their union coverage (i.e., co-occurrence ratio). There are same meanings in the following tables. MG₁ is a metagene including seven genes: CDK4, FAM119B, MARCH9, TSFM, CENTG1, METTL1, TSPAN31. MG₂ is a metagene including four genes: WT1, SLC1A2, PAX6, ABCC4.

For k=5, two gene sets including {CDKN2A, TP53, MG ₁} and {CDKN2B, CYP27B1} were detected with p₁=0.0003, p₂=0.0018, p_1,2<0.0001 and r_1,2=0.7606. CDKN2B encoding INK4B (also known as p15) also locates in chromosome 9p21 homozygous deletion region, and CDKN2B is usually co-deleted with CDKN2A. This disrupts the p53 and RB pathways. For this reason, combinatorial inactivation of CDKN2A and CDKN2B is frequently observed in these tumors. The cross-talk between the p53 and RB pathways (Figure 3B) suggests that CDKN2A, TP53 and CDK4 are in the same pathway (Figure 3C(a) or 3C(b)).

For k=10, two gene sets {CDKN2A, PTEN, CYP27B1, KDR, MG ₂} and {CDKN2B, NF1, RB1, MG ₁, ERBB2} with p₁,p₂ and p_1,2 less than 0.0001 and r_1,2=0.9333 were identified (Table 1 and Figure 4A). The first gene set was found to be involved in the p53 and PI3K/Akt signaling pathways and the second in the RB and RTK/RAS/ERK signaling pathways. The RTK/RAS/PI3K signaling pathway can also be induced by the mutations in these two gene sets (Figure 4B). These pathways are implicated in biological processes associated with cell survival, cell cycle, protein synthesis, and cell proliferation. p53, RB, and RTK/RAS/PI3K have been previously reported to contribute to GBM pathogenesis in original TCGA GBM studies [1]. Five well-known tumor suppressors (CDKN2A, CDKN2B, PTEN, NF1, and RB1) are involved in these two gene sets. Besides the co-occurrence of CDKN2A and CDKN2B, NF1 and RB1 in the second gene set have exclusive mutations, which are co-occurrent with mutations of PTEN in the first set (Figure 4A). Recently, several studies have shown the cooperativity of tumor suppressors in carcinogenesis [33–35]. For example, Rahrmann et al. demonstrated that co-occurring mutations in PTEN and NF1 cooperate in the development of grade 3 PNSTs (peripheral nerve sheath tumors) in mice, suggesting that they may cooperate in human MPNST (malignant PNST) progression [33]. Another study by Chow et al.[34] showed that cooperativity among PTEN, TP53, and RB1 can cause high-grade astrocytoma in mouse adult brain, in which the majority of glioblastomas arise. For another two almost simultaneously mutated oncogenes CYP27B1 and CDK4 (Figure 4A), Beckner et al. have demonstrated their cooperative amplification and co-expression for potential modulation of vitamin D in glioblastomas [36].

On the other hand, WT1 encodes a transcription factor that plays an essential role in cellular development and cell survival [37]. It regulates the expression of numerous target genes, including the famous tumor suppressor TP53 and the Wnt signaling pathway [38]. KDR encodes a VEGF (vascular endothelial growth factor) receptor. VEGF plays a crucial role in angiogenesis and progression of malignant brain tumors. ERBB2 encoding the protein HER2 (human epidermal growth factor receptor 2) is a member of the epidermal growth factor receptor (EGFR) family, and it has been shown to play an important role in the pathogenesis and progression of many different types of cancer. WT1, KDR and ERBB2 may drive the carcinogenesis of GBM, indicating that CoMDP can identify low-frequency candidate driver genes that play important roles in cancer initiation and development.

Lung cancer

For k=10, {STK11, ATM, TP53, PAK4} and {KRAS, NTRK3, EGFR, GNAS, EPHA3, NRAS} were identified. All four genes STK11, ATM, TP53, PAK4 have been demonstrated to be closely related to the p53 signaling pathway in lung cancer [40, 41]. Two members of the RAS subfamily, KRAS and NRAS function as binary molecular switches controlling the intracellular signaling networks that regulate several key cancer-related processes, such as proliferation, differentiation, cell adhesion, apoptosis, and cell migration. GNAS is a guanine nucleotide-binding protein (G protein). It acts as a modulator or transducer in various transmembrane signaling systems. GNAS may interact with MDM2, which may lead to MDM2-mediated degradation of TP53. Solomon et al. found that many kinds of TP53 mutations can regulate RAS in different ways, inducing a cancer-related gene signature [42]. Kosaka et al. demonstrated that TP53, EGFR and KRAS may cooperatively determine the prognosis of the patients in lung adenocarcinoma [43].

GBM2 dataset

We observed that some new co-occurring gene sets were identified for GBM2 compared to GBM1 (Table 3). For k=9, we identified {CDKN2A, TP53, MG ₃, PIK3R1, TAF1} and {CDKN2B, CYP27B1, RB1, SYNE1} (MG ₃ is a metagene including CDK4, MARCH9 and TSPAN31) with p₁=0.0002, p₂<0.0001, p_1,2<0.0001 and r_1,2=0.8642 (Figure 5A and Table 3). In addition to the cooperative effects between CDKN2A and CDKN2B, TP53 and RB1 have been reported to have frequently co-occurring mutations related to several cancers, including the central nervous system tumor [13]. Recently, the collaboration of TP53 and CDKN2B was also studied with respect to cell apoptosis and aneurysm formation [44]. On the other hand, for the two detected low-frequency mutated genes TAF1 (2/170) and SYNE1 (3/170), TAF1 encoding a transcription initiation factor phosphorylates TP53 during G1 cell-cycle progression, so TAF1 may be a member of the p53 signaling pathway; SYNE1 was found to be associated with the GBM patients’ lifetime, and was therefore considered to be an important biomarker of glioblastoma survival [45]. Our studies indicated that the p53, RB, and the PI3K-related signaling pathways may collaboratively contribute to carcinogenesis in GBM via combined genetic alterations (Figure 4B and Figure 5B).

k	Gene set 1	Gene set 2	p 1	p 2	n 1	n 2	r 1,2	p 1,2
4	CDKN2A, M G 3	CDKN2B, CYP27B1	0.0045	0.0056	60	63	0.9219	<0.0001
	CDKN2A,
5	CYP27B1, COL1A2	CDKN2B, M G 3	0.0072	0.0018	62	63	0.9531	<0.0001
	CDKN2A,	CDKN2B,
6	CYP27B1, COL1A2	M G3, ERBB2	0.0073	0.0003	62	65	0.9538	<0.0001
	CDKN2A, TP53,	CDKN2B, RB1,
7	M G 3 , TAF1	CYP27B1	0.0001	0.0001	77	70	0.8375	<0.0001
	CDKN2A, TP53,	CDKN2B, RB1,
8	M G 3 , TAF1	CYP27B1, PRNP	0.0002	<0.0001	77	71	0.8500	<0.0001
	CDKN2A, TP53, TAF1,	CDKN2B, RB1,
9	M G 3 , PIK3R1	CYP27B1, SYNE1	0.0002	<0.0001	79	72	0.8642	<0.0001
	CDKN2A, TP53, TAF1,	CDKN2B, RB1, SYNE1,
10	M G 3 , PIK3R1	CYP27B1, M G 4	0.0001	<0.0001	79	73	0.8765	<0.0001

MG₃ is a metagene including three genes: CDK4, MARCH9, TSPAN31. MG₄ is a metagene including 168 genes.

Ovarian cancer

The mutation distribution among genes in the ovarian carcinoma data is quite nonuniform. Among all the 314 samples, TP53 is mutated in 251 of them and all the other genes are mutated in less than 26% of samples. This indicates that TP53 plays a crucial role in the carcinogenesis of ovarian cancer (TTN was removed in the present analysis because of the possible artifacts of its mutations [46]). Determining whether there are other driver genes or pathways collaborating with TP53 will be helpful for understanding the pathogenesis of this cancer.

To identify other driver gene sets coupled to TP53, we applied mod_CoMDP with r=1∼10 to the ovarian cancer data and significant results were obtained for r=3∼10 (Additional file 1: Table S1). For example, for r=10 we identified {MYC, CCNE1, NINJ2, M G₅, USH2A, NF1, HMCN1, ZNF596, USP35, M G₆} (M G₆ is a metagene including four genes: STMN3, SLC2A4RG, ZGPAT, RTEL1) with r a=0.6563 (the co-occurrence ratio with TP53). Frequent somatic mutations in NF1 have been previously shown to co-occur with TP53 mutations in ovarian carcinomas [52, 53]. STMN3 and NF1 have been demonstrated to be involved in the MAPK signaling pathway [19]. Furthermore, to discover possible collaborations of multiple driver pathways with T P53, we combined TP53 and the aforementioned 10-gene set into one nominal gene, which was considered mutated in a sample if both sets were mutated in that sample. Then we applied mod_CoMDP to identify gene sets significantly co-occurring with the nominal gene. For r=1∼10 we identified PPP2R2A, which is generally implicated in the negative control of cell growth and division. Kalev et al. revealed that PPP2R2A plays a critical role in DNA double-strand break repair through modulation of ATM phosphorylation [54]. Youn and Simon recently studied mutator alterations relevant to ovarian cancer [55]. Besides the well-known mutator gene TP53, they identified PPP2R2A and the chromosomal region 22q13.33 as the new mutator candidates. We find that these so called mutator genes, which increase genomic instability when altered, may be collaboratively involved in the processes of DNA synthesis and repair, chromosome segregation, damage surveillance, cell cycle checkpoints, and apoptosis. The discovered driver patterns here may provide new information to enhance our understanding of the ovarian carcinoma pathogenesis, and further explorative analysis is needed to verify their biological relevance.