Identification of ovarian cancer associated genes using an integrated approach in a Boolean framework

Identification of differentially expressed genes

We extracted and analysed TCGA (The Cancer Genome Atlas) level 3 (Batch 9) ovarian serous cystadenocarcinoma data from the Affymetrix platform [17]. TCGA gene expression data are normalised, annotated and validated for expression variation relevant to the tissue types and with the type of array platforms, thus increasing the robustness in analysing expression data. Rather than a fold-change, we have calculated the differential expression of each gene by considering the percentage of false prediction (pfp) ≤5% using the RankProd R package [18]. RankProd uses the rank product non-parametric method to indentify up/down-regulated genes under one condition against the other (in our case tumor vs. normal ovarian samples). This is based on the null hypothesis that the order of all items is random and the probability of finding a specific item among the top r of n items in a list is p = r/n. Multiplying these probability leads to the identification of the rank product $\mathit{RP}={\displaystyle \prod _i\frac{r_i}{n_i}}$, where r _i is the rank of the item and n _i is the total number of items in the i ^th list. The smaller the RP value, the smaller the probability that the observed placement of the item at the top of the list is due to chance.

Relevant functional attributes in the disease condition

Although microarrays measure the relative abundance of mRNA transcripts, their translated proteins are likely to be differentially present in diseased tissue. Moreover, the extent of differential protein concentration under the disease condition is quite difficult to estimate due to the heterogeneity of cells in the tumor sample. Therefore, we considered a Boolean combination of six proteins functional attributes for searching genes associated with ovarian cancer, where the causative effects are not additive but combinatorial as well as non-linear. These functional attributes are tissue specificity (TS), transcription factors (TFs), post-translation modifications (PTMs), protein kinases (PKs), secreted proteins (SPs) and whether the protein is a hub in the interactome, with node connectivity greater than four (i.e. node connectivity > = 5) along with the gene attribute of methylation (METH), in cancer vs. non-cancer associated genes.

We hope to capture the underlying enabling factors for cancer, by considering the above protein functional attributes. Large-scale data analysis supports the fact that disease genes are generally tissue-specific and are over-expressed in those tissues where changes in gene expression result in pathology [19]. TFs are DNA-binding proteins regulating gene expression and thereby control cell development, differentiation and growth [20] and their aberrant activity has been implicated in the cancer disease condition [21]. Oncogenic conversion of normal cells into cancerous cells involves changes in transcription factor, e.g. c-Fos component of TF c-Jun/JUN/AP-1 is crucial for the estrogen receptor α (ERα) mediated transcription in breast cancer [22]. PTMs of key regulatory or structural proteins are known to play an important role in the progression of cancer by activation of signalling pathways, enhanced proliferation and impaired cell division and death [23]. PTMs contributing to tumorigenesis include phosphorylation, acetylation, methylation, glycosylation, prolyl isomerisation, hydroxylation, oxidation, glutathionylation, sumolyation and ubiquitination. For example, clinical evidence suggests that phosphorylation, acetylation and sumolyation of ERα lead to prostate and breast cancer in humans [24]. PKs are important signalling molecules for maintaining normal tissue architecture and function, hence mutation in these genes are a common cause of human cancer (http://www.sanger.ac.uk/genetics/CGP/Kinases/) [25]. Recent developments in proteomic analyses suggest an increasingly large number of genes overexpressed in ovarian cancer, of which several encode secreted proteins [26]. For example, the high expression of prostasin [27] and osteopontin [28] are recorded in the serum of ovarian cancer patients. Highly connected proteins, i.e. hubs are shown to be essential in connecting diverse functional modules inside the cell [29]. Also, epigenetic inactivation of tumor-suppressor genes due to methylation is well known in carcinogenesis [30].

Data integration from multiple experiments

These filters were used to define a HC protein interaction set to study the network properties of molecular functions and biological processes of interacting proteins. In this study, scoring schema for interactions were considered for those protein nodes with more than four interactions, as this is the empirical value of hubs suggested in gene co-expression stability in the analysis of protein interaction networks [46]. Therefore, we weighted such highly connected protein nodes encoded by the known cancerous genes.

Weighting schema for Boolean-based probability calculation

where N is the total number of genes.

Scoring schema on the weighted functional attributes for ranking genes

We used the Boolean algorithm proposed by Nagaraj and Reverter [15] for ranking the differentially expressed genes in ovarian samples, with our own set of Boolean variables representing relevant functional attributes in the disease condition. The particular combination across the seven Boolean variables i.e. functional attributes for a given differentially and non-differentially expressed genes, was decomposed into its root. For example, if a given gene has four known functional attributes, then 2⁴ Boolean states are known to exist containing (2⁴-1) roots, i.e. all possible combinations of Boolean states at the positions of known functional attributes, excluding the Boolean value with all zero status. The probability of each root is simply the average sum of all the weights associated with known functional attributes calculated via r _Φ . These root probabilities are then used to rank the differentially and non-differentially expressed genes by summing up all the probability values associated with the individual roots.

Validation set

We retrieved the raw expression data for 153 ovarian tumor samples from the Gene Expression Omnibus entry GSE1349, containing samples in four tumor stages [9]. Raw expression values for each probe were transformed to log-scale with base 2. Probe IDs were converted to Entrez Gene IDs using AILUN [47]. For genes with multiple probes, the probes with the highest variance across the samples were used to describe the expression value for the genes. Probes with multiple or without Gene IDs were removed from the analysis. Pearson’s correlation coefficients were calculated based on the co-expression values alone, to define the pairwise gene co-expressions. We have taken a Pearson’s coefficient > 0.5 to define a link between co-expressed genes in the gene expression network.

Network analysis of human signalling

We performed network analysis using the manually curated human signalling network [48]. The signalling network was pruned to contain associations between proteins alone and hence, small molecules were removed from the network, resulting in 1522 protein nodes and 4276 edges. The R package, igraph [49] was used for the network analysis. The Ingenuity Pathway Analysis system, (IPA, http://www.ingenuity.com) was further considered, to interpret the interaction of cancerous genes in humans.

Clinical characterization using survival analysis

For the high scoring Boolean-based differentially expressed genes, we performed Kaplan-Meier survival analysis using the Cox-Mantel log-rank test, implemented as an R package. The significance of these genes in the normal and tumor samples were evaluated from the Welch two sample t-test. A less conservative P- value < 0.1 was considered for the statistical significance of genes in the ovarian tumor sample classification. The clinical data were downloaded from the TCGA data portal.

Protein kinases

Protein kinases are important regulators of cell function and belong to a functionally diverse gene family. They affect the activity, localisation and overall function of other proteins by adding a phosphate group and thereby control the activity of cellular processes. Kinases are particularly important in signal transduction and co-ordination of complex functions such as cell cycle and pathological conditions. Identification of IRAK1 as a differentially expressed gene in ovarian cancer suggests its important role in this disease. It is a putative Ser/Thr kinase known to partially interact with transcription factor, NF-κB. Activation of NF-κB leads to cell proliferation, survival and migration [50]. Over-expression of this gene suggests indirect cell survival and proliferation in the ovarian tumor condition. Similarly, IGF1R is a receptor with tyrosine kinase activity, which binds an insulin-like growth factor. It is over-expressed in most malignant tissue, acting as an anti-apoptotic agent by enhancing cell survival [51, 52]. LYN is a non-receptor tyrosine kinase, phosphorylating caspase 8, rendering it inactive and thereby assisting apoptosis of the inflammatory cell [53]. In the absence of the normal expression of LYN, active caspase 8 may prevent the tumor cells from undergoing apoptosis.

Other important kinases in cell survival and proliferation during tumorigenesis are associated with key cell cycle proteins. CDC7 (cell-division cycle 7 homolog of S. cerevisiae) and BUB1 (budding uninhibited by benzimidazoles 1 homolog of S. cerevisiae) encode protein kinases which induce G1/S transition and are involved with the spindle checkpoint function, respectively during cell mitosis. CDC7 is known to be overexpressed in the epithelial ovarian carcinoma, resulting in tumor progression, genomic instability and accelerated cell division [54]. On the other hand, BUB1 overexpression induces aneuploidy and tumor formation [55]. CHEK1 (checkpoint kinase 1) is an another important cell-cycle molecule of Ser/Thr protein kinase family mediating signals from ATM and ATR cell cycle proteins involved in the DNA damage response and associated with chromatin in the meiotic prophase I. The importance of this protein in tumor invasiveness has been suggested by researchers in lung, bladder, liver, prostate, gastric, brain, cervical and colorectal cancers and B-cell lymphoma [56–58]. CHEK2 (checkpoint kinase 2) is yet another important cell cycle protein which regulated key proteins during cell division. It interacted with BRCA1 (breast cancer 1) to restore survival in response to DNA damage with known association with endometrial cancer risk [59]. We observed overexpression of IRAK1, BUB1, CDC7, CHEK1 and CHEK2 genes in TCGA samples at a high Boolean probability score of 0.607561, together with the co-expression of other key cell-cycle molecules in an independent validation expression set GSE1349 suggesting their association in ovarian cancer (Figure 4 and Additional file 5). The presence of high probability up-regulated genes in the co-expression network (Pearson’s correlation coefficient > 0.5) is shown in Additional file 6. The co-expression network of downregulated genes is available from Additional file 7.

Serine proteases

Serine proteases are proteolytic enzymes, hydrolysing the peptide bond of protein substrates via a nucleophilic serine residue in the active site [60]. Serine proteases play diverse roles in human health, from non-specific digestion to highly regulated functions like embryonic development, immune response and blood coagulation. Moreover, insufficient or excess protease activity can promote significant pathologies like cancer, inflammation, hemophilia, heart attack, stroke, pancreatitis and parasite infection [61]. We suggest the potential use of KLK6 (kallikrein-related peptidase 6) as a potential biomarker for ovarian cancer based on its high Boolean probability score (0.697808) (Figure 4). KLK6 is a serine protease with diverse functional roles inside the cell. It has been suggested that overexpression of this protein leads to the loss of cell-cell adhesion in skin cancer (melanoma) [62]. Moreover, a recent study reports the up-regulation of KLK6 in colon cancer and its use as a potential biomarker and therapeutic agent [63].

Secreted proteins

Secreted proteins are secreted from the cell into the extracellular space and have important biological regulatory roles with the potential for therapeutics. STC2 (Stanniocalcin 2) is a secreted homodimeric glycoprotein that is expressed in a variety of tissues. STC2 is known to promote the epithelial-mesenchymal transition and invasiveness in human ovarian cancer under inadequate oxygen supply to the tissue [64]. Our results show that STC2 is a significant up-regulated gene, promoting ovarian cancer. On the other hand, CLU (clusterin) and LCN2 (lipocalin 2) are down-regulated genes in our analysis. CLU encodes a protein which is secreted under stress conditions, that functions as a strong anti-migratory and anti-invasive agent by inducing the destruction of the actin cytoskeleton inside the cell [65]. The decreased expression of CLU thus promotes the cancerous disease condition. LCN2 encodes a 25 kDa secretory protein involved with iron-transportation and contributes to endometrial carcinoma [66]. Moreover, it is a key molecule in various signalling pathways (Additional file 4). Down-regulation of LCN2 due to epigenetic inactivation may lead to ovarian carcinoma.

Other types of proteins

We observed down-regulation of genes with high probability associated with phosphoproteins, transcription factors and receptors due to epigenetic inactivation. Phosphoprotein DAB2 is a mitogen-responsive agent, acting as tumor suppressor in normal ovaian epithelial cells and down-regulation of this gene modulates the TGF-β signalling pathway [67]. FOXL2 (forkhead box L2) encodes a transcription factor which helps in the normal development of ovarian tissue. IGFBP7 (insulin-like growth factor binding domain) is known as the tumor suppressor gene, leading to lung cancer due to the epigenetic inactivation [68]. PGR (progesterone receptor) encodes a protein playing a central role in the reproductive system by maintaining progesterone levels and ensuring normal pregnancy. AR (androgen receptor) encodes a protein which functions as a steroid hormone-activated transcription factor and has been shown to be involved in prostate cancer [69] as well as in ovarian cancer in association with p44[70]. VIM (vimentin) encodes a protein that is responsible for maintaining cell shape, integrity of the cytoplasm and stabilizing cytoskeleton interaction. Thus, the decreased expression of these genes could be indicative of ovarian cancer.

Relevance to cancer

The complexity of cellular regulation is encapsulated in the signalling network. Moreover, signalling network also helps to understand the dynamics of oncogenic cross-talk by determining the sites at which oncogenic signals occur and through which oncogenic signals are transduced. Therefore, differentially and non-differentially expressed genes with high Boolean scores identified in the TCGA dataset (listed in Table 2 and Additional file 3, respectively) were mapped to the manually curated human signalling network [48], shown in Figure 5. Nine of the 17 differentially expressed genes and 32 of 48 non-differentially expressed genes could be mapped to the signalling network (Figure 5A). We have computed the importance of a node in the signalling network using betweenness centrality (the number of times a node acts as a bridge along the shortest path between two other nodes) and the influence of a node in the network using eigenvector centrality (relative scores to all nodes in the network with connections to high-scoring nodes receiving higher scored compared to low-ranking nodes). A plot of the two signalling network centrality measures (betweenness centrality and eigenvector centrality; Figure 5B) suggests the importance of the differentially expressed AR and the non-differentially expressed MAP3K7 (mitogen-activated protein kinase kinase kinase 7, involved in stress response) genes, which have ten shared interactions, from BioGraph [73], ranked 636 out of 18180 ‘Gene’ concepts (top 3.50%) in BioGraph’s knowledge base. Moreover, the nine differentially expressed genes (AR, CHEK1, CHEK2, PGR, VIM, LYN, IRAK1, IGF1R and DAB2) mapped to the signalling network were further identified using the Ingenuity Pathway Analysis system to visualize the interaction of these genes with the known oncogenes (Figure 6). The central role played by CHEK1 in the DNA damage response signalling network, has been confirmed by Dai and Grant [74], where CHEK2, CDC7 and BUB1 have also been identified from the 17 differentially expressed genes reported here.

Clinical characterization

Table 2 lists 17 genes, of which seven are up-regulated and ten are down-regulated in ovarian cancer patients. The expression patterns of these genes suggest that the sum of the up-regulated gene expression values minus the sum of the down-regulated gene expression values should be maximized in ovarian cancer patients compared to controls without ovarian cancer (Additional file 10). Figure 7 shows that this is indeed the case for the 38 ovarian clinical samples and seven normal samples in this dataset and that this simple formula for the 17 genes identified here can be used to successfully distinguish between normal and ovarian cancer patients (p-value < 1.2E-06).

Survival analysis was carried out to suggest if any of above 17 genes or their combinations, can be used in the classification and prognosis of ovarian cancer, to classify good and poor prognostic tumors. To demonstrate the survival analysis across the 38 ovarian clinical samples in this dataset, expression levels of each of the 17 genes were ranked from lowest to highest expression. Tumor samples associated with the lower 50% of the expression values for a given gene were labelled as “low-expression” for that gene; otherwise, they were labelled as a “high-expression” sample for that gene. Log-rank tests were then performed to suggest the difference between expected vs. observed survival outcomes for the low- and high-expression tumor samples for each of the genes. As there were only 38 ovarian tumor samples with clinical data, we chose the less stringent log-rank P-value of 0.1 and discovered three genes, CHEK1, AR and LYN exhibit a prognostic value, based on this cut-off level (see Figure 8).

In Figure 8, the lower of the two curves in each of the four survival analysis plots indicates tumor samples associated with poor prognosis. Interestingly, though the survival curves associated with gene AR indicate poor prognosis is expected for tumor samples within the high expression range of AR, from Table 2 we note that AR is down-regulated in ovarian cancer. From Figure 8, it is seen that high expression for up-regulated CHEK1 and down-regulated AR and low expression for LYN leads to poor prognosis. The clinical data thus suggests a preference for limited down-regulation of AR. Therefore, combining the expression levels of these three genes as CHEK1+AR-LYN (Figure 8), then ranking this score from lowest to highest values and associating the patients into low and high expression groups, as before, gave greater significance in the prognostic outcome for classifying good and poor tumour outcomes than did the individual genes. Biologically, this combination represents increased cell cycle control, particularly for entry into mitosis (CHEK1), decreased expression of the androgen receptor (AR), whose expression levels have controversial reports as a favourable prognostic factor in epithelial ovarian cancer [75, 76] and moderately decreased expression of LYN, resulting in apoptosis of tumor cells.