Data

The Level 3 dataset of gene expression, copy number and DNA methylation for the same set of ovarian cancer samples (Table 1) were obtained from the publicly available TCGA website (https://tcga-data.nci.nih.gov/tcga/). Gistic2.0 was used to analyze the copy number dataset (Level 3) for the identification of recurrent regions of copy number alteration and the copy number of genes. The beta values of DNA methylation are continuous, ranging from 0 (unmethylated) to 1 (completely methylated). The probe IDs were mapped to Gene symbols with the annotation table for Illumina Human-Methylation27 platform, which detected the methylation level of 27,578 CpG loci located within the proximal promoter regions of transcription start sites of 14,495 genes. If there were multiple probes corresponding to the same gene, we adopted the averaged intensity of these probes as the beta value of the gene and removed the probes with no value or corresponding gene. We selected a K-nearest neighbor-based method that imputes missing values in gene expression profiles, which was implemented by an R package (impute). In addition, we have added a list of the samples into supplementary material (see S1 Table).

To integrate HPRD[17], Reactome[18], MSKCC Cancer Cell Map, and the NCI/Nature Pathway Interaction Database[19], Pathway interaction data and protein-protein interaction data were used to establish the initial network. Pathway data sets for Reactome, the NCI/Nature Pathway Interaction Database, and the MSKCC Cancer Cell Map were downloaded in the Simple Interaction Format (SIF) format from Pathway Commons, protein-protein interaction data was downloaded from HPRD. The Human Background Network (HBN) was the unified set of the four dataset. Simultaneously, redundant edges and self-connected edge were removed (Table 2).

The HBN we built consists of genes and interactions in the forms of nodes and edges. The interaction reflect the functional associations between two genes, such as a physical interaction, or an indirect interaction via the common pathway.

We acquired 973 seed genes (S2 Table) from four well-established cancer- and disease-related gene databases: Cosmic[20], GAD[21], OMIM[22] and phenopedia[23]. Ovarian cancer seed genes were defined as known oncogenes or tumor suppressor genes associated with cancer in the well-known databases. The workflow of our approach is depicted in Fig. 1 and further details are provided in the next section.

Difference analysis of genes in a single level

Gistic2.0[24] was used to analyze the copy number dataset to identify recurrent regions of copy number alteration and the copy number of genes. We identified a number of recurrent focal somatic copy number alteration (SCNA) events, including 55 significant amplifications and 48 deletion peaks. The SAM[25] algorithm was applied to two sets of ovarian samples (tumor/normal) to identify differentially expressed genes: we identified 549 highly expressed genes and 805 low-expressed genes that were differentially expressed in cancer (fold change > = 2 and false discovery rate (FDR) <0.05). For DNA methylation data, we identified highly significant (FDR<0.005) differentially methylated genes in tumor samples compared to normal samples using the Mann-Whitney-Wilcoxon test, including 1445 hypermethylated genes and 1219 hypomethylated genes.

The construction of the integrated co-alteration network and performance comparison

To simultaneously use multiple features of genes and establish the correlation between genes at the genome, epigenome and transcriptome level, we designed a framework based on CCA, a statistical method used to analyze the degree of correlation between two sets of random variables. CCA can turn the ordinary correlation between two variables into the canonical correlation between two sets of variables. The purpose of CCA is to seek maximization of the correlation between two linear combinations of the variables[26, 27].

In this work, the features of genes were seen as random variables; the possibility of two genes being co-altered on all levels was then measured by the following procedure.

We defined two genes: g₁, g₂. Suppose that G₁ = [g₁⁽¹⁾, g₁^{(2) …},g₁^(p)]^T, G₂ = [g₂⁽¹⁾, g₂^(2)…,g₂^(p)]^T, and the two vectors consist of p types of information of g₁ and g₂. In this study, we set p = 3.Take G₁ for example: g⁽¹⁾ denoted the expression values of g1 in samples, g₁⁽²⁾ denoted the copy number values of g1 in samples, and g₁⁽³⁾ denoted the methylation values of g1 in samples. Similarly, we can define G₂.

Let,

Then the covariance matrix is defined as:, in which each element is calculated by formula (1).

(1)

We use the correlation of linear combination of vectors (namely a^TG₁, b^TG₂) to measure the linear relationship between G₁ and G₂.

The construction of ICan was implemented by seeking the maximum correlation coefficient between U = a^TG₁ and V = b^TG₂

(2)

Solutions to the optimization problem (2) satisfied the conditions: Var(a^TG₁) = 1, Var(b^TG₂) = 1.

Our purpose was to seek the most suitable a and b such that corr(U,V) was the largest. The first pair of linear combinations was called the first pair of canonical variables; their largest correlation ρ(U₁,V₁) was called the first canonical correlation. Next, if there exists a_k and b_k such that the following conditions were satisfied:

1. was uncorrelated with initially K-1 pair canonical variables;
2. 
3. The correlation coefficient between and is the largest.

were called the first K pair of canonical variables and ρ(U_k,V_k)was called the first K canonical correlation. In this study, we set K = 3. The Rayleigh quotient matrix:.

The first correlation coefficient is equal to the square root of the largest eigenvalue λ₁ of the matrix R. Similarly, the first K correlation coefficient is equal to the square root of the largest eigenvalue λ_k of the matrix R. After that, the linear correlation coefficient (ρ₁,ρ₂,ρ₃) was calculated between every gene pair in the data set.

Canonical correlation is an extension of ordinary correlation; it can measure the correlation between two sets of variables[28]. Compared with using a single data type, it showed more accuracy in the quantification of the linear relationships between genes using their different features[29]. Next, similar to previous works[29], we used the chi-squared test to measure whether the canonical correlation coefficient (ρ₁,ρ₂,ρ₃)[30] was significant.

The null hypothesis is H₀: λ_k = … = λ_p = 0

Let P_k be the p-value of the K-th test statistic T^k, with:, and T^k ~ [29], where n is the number of samples. Finally, we used a combination of weights (3) to assign a weight to the edges connecting two genes, (3) Where

The final weight, ω, represents the correlation between genes more precisely. ω measures the possibility of two genes being co-altered on the level of copy number, DNA methylation and gene expression. We then assigned the weight to the HBN and constructed the integrated co-alteration network referred to as ICan. The method can measure the strength of association between genes on multiple levels. In this work we implemented the CCA method and chi-square-based statistical significance test by the library "CCA" and “Chi-square test” in the R statistical software.

Meanwhile, we computed the Pearson correlation coefficient of the expression profiles (copy number profiles and methylation profiles) between every pair of genes and established a co-expression network(GCE), a co-copy number network(GCC) and a co-methylation network(GCM). This process was also implemented in the R statistical software. To better reflect the performance of our network, we compared ICan and CNAmet, and between three single data networks.

Identifying candidate ovarian cancer-related genes

Random Walk with Restarts[31] is a sorting algorithm. It simulates the process of walking step by step from seed nodes to direct neighbor nodes; nodes in the network are ranked by the probabilities of reaching the node. Assuming W is the adjacency matrix of the ICan and P_t is a vector whose i-th element holds the probability of arriving at node i at step t, the random walk was computed by

(4)

The distribution of values of seed nodes in the initial probability vector P₀ was set as uniform, with the sum of the probabilities equal to 1; r represents the probability to restart at seed nodes, which was set to 0.7. After N steps, this probability will reach a steady state, which was determined by the difference between P_t and P_t+1. We performed the iteration until the L1 norm between them fell below 1E-10. The Random Walk with Restarts probability for all the genes in the network was calculated. We then analyzed the differential alteration of the top 20% genes in the various levels.

Kaplan-Meier survival analysis for candidate cancer-related genes

A non-parametric Kaplan-Meier estimator was applied to estimate the influence of different factors on survival time. In this work, to explore the possible prognostic value of identified candidate genes, we used the “survival” package in the R statistics software. A p-value <0.05 and an FDR < 0.25 were used as a cutoffs for statistical significance by the log-rank test.

We investigated the alteration of each gene in the samples, and discretized the three datasets according to the features of oncogenes and tumor suppressor genes, i.e., amplification, overexpression, hypomethylation; and the reverse: deletion, low expression and hypermethylation, respectively. For copy number data, we adopted the results of GISTIC2.0 discrete copy number calls. The samples were classified as gene homozygous deletion (-2) or amplification (1/2). For the gene expression data, we calculated the mean value and standard deviation (SD) for each gene: the values that were higher than mean + SD were considered overexpression. Conversely, the values that were lower than the mean—SD were considered low expression. For the DNA methylation data, we set the threshold based on empirical analysis of the beta value distributions: a beta value less than 0.2 was regarded as hypomethylation; a value more than 0.8 was regarded as hypermethylation.

Identifying functional modules for ICan

We identified functional modules from ICan and constructed three single-level networks using MCODE[32]. The use of MCODE was preferred for an easier comparison of ICan and the three single-factor networks, as the same modules were identified from the unweighted network. The edge-weighting procedure was performed separately for each network, and the M scores of each module were calculated according to a scoring formula (see Additional file S4 Table for details). A functional enrichment analysis was performed on the candidate cancer-related gene set and the genes inside the module using the DAVID tool[33] (http://david.abcc.ncifcrf.gov/).

ICan has the properties of complex networks

The integrated co-alteration network is represented as an undirected weighted graph, where nodes represent genes and edges connecting the nodes represent the correlations of co-alteration between genes. First, making use of human interaction data and pathway knowledge, we established an HBN that comprised 9,195 nodes and 65,720 edges.

In 574 ovarian cancer tumor samples, there are 11,384 genes that are present in all three profiles of copy number, promoter methylation and gene expression. According to CCA, we then calculated the weight between every two genes to measure their linear correlation by the three features. Next, the edges in the network were assigned weights and the genes not contained in molecular profiles were removed. Eventually, we constructed ICan, which comprised 6,345 nodes and 40,125 edges. The closer ω is to 1, the higher the correlation between the two genes. In addition, we used the Pearson correlation coefficient for the levels of gene expression, copy number, and DNA methylation to construct three same sized networks.

Network topology plays an important role in the biological functions and information transmission in the network. After analyzing the properties of the network topology, we found that ICan showed a scale-free structure, with a power-law distribution of node degrees. This means that ICan includes only a small number of nodes whose degree is high, suggesting the importance of the hub nodes. We then applied the weighted random walking method to identify hub nodes. This method can effectively optimize candidate disease genes and accurately predict candidate key genes of cancer.

ICan improves the accuracy of prioritizing candidate cancer-related genes

ICan contains 604 known ovarian cancer-related genes, which were used as the gold standard to plot receiver operator characteristic curves, and to calculate the area under the curve (AUC). Based on five-fold cross validation, we selected 80% of the genes as seed genes; the remaining 20% were reserved for final validation. To prove the accuracy of our method, using the same data set, we applied the CNAmet method to predict oncogenes and tumor suppressor genes, and compared the outcomes with the ICan outcome. As a result, the AUC value of CNAmet was significantly less than the AUC value of ICan (ICan: the Max AUC = 0.8179; CNAmet: AUC = 0.5183, p-value = 3.158e-14, the first two sheets in S5 Table) (Fig. 2). The significance of the difference of the AUC for two ROC curves was determined by DeLong's test in “pROC package”[34].

To more accurately predict the cancer-related genes in ovarian cancer, we used a weighted random walking method to calculate the proximity between other nodes and seed genes to determine correlations with oncogenes. This method is often referred to as the “guilt-by-direct-association” principle, by which the genes that are associated with disease genes tend to have similar functions. We randomly chose genes in ICan as seed genes, and compared them with the original results. This process was repeated 1000 times; an adjusted p-value below 0.05 was considered significant for cancer-related genes. On the other hand, we compared the difference in the degree[35] and gene length between candidate genes and the other genes. Recent research has shown that a greater gene length often results in more domains in the translated proteins, thus leading to greater interactivity, which means a greater possibility of the gene being cancer gene[36]. The results showed that not only were there significant differences in the gene length of candidate cancer-related genes compared with the other genes (p-value = 2.64E-02, Fig. 3, S6 Table), but also the results were similar in terms of gene degree (p-value = 6.176E-07).

Finally, we identified 155 candidate cancer-related genes (S7 Table), and analyzed the co-alteration events of these genes in detail. CHEK1, IGF1R and MSH3 were co-altered in common at all three levels; CHEK1, IGF1R, MSH3 and FANCA were co-altered at the copy number and expression levels; and CHEK1, FGF18, IGF1R, IGFBP1, IGFBP2, MSH3, PLAU, RAD51 and EIF2AK2 were co-altered at the level of DNA methylation and expression.

CHEK1, FANCA and RAD51 are involved in the inspection of breakpoints in the cell cycle regulation and repair process, and play important roles either in the p53 signaling pathway or the MAPK signaling pathway. The MAPK signaling pathway is an important cancer pathway; activation of this pathway can promote endothelial cell proliferation and angiogenesis. The newly generated blood vessels could provide more nutrients to tumor cells, accelerating tumor growth and promoting proliferation of cancer cells[37]. MSH3 and IGF1R have important roles in DNA replication, recombination, and repair. Deficiency of mismatch repair, especially loss of expression of the seven main genes (MSH2, MSH3, MSH6, MLH1, MLH3, PMS1 and PMS2), can increase the risk of ovarian cancer[38].

In addition, we analyzed the differential proportion of the top 20% genes in ICan by random walking. Fig. 4 shows that the proportion of differential methylation was the highest in each bar among the top 100; however, only two genes have simultaneous differential changes on all three levels. The numbers of genes with only one type of alteration (CNA, differential methylation or differential expression) were 13, 19 and 18, respectively. We found that the number of genes that were differentially altered on multiple levels tended to stabilize after the top 600, which indicated that the probability of these genes is much higher, suggesting a closer relationship with known seed genes.

The alteration of a gene on a single level represented a copy number abnormality, differential expression or differential methylation, respectively (S3 Table, sheet 1–3).

Novel cancer-related genes of ovarian cancer may affect survival

To estimate the impact of candidate genes on patient survival, and look for genomic and epigenetic genomic features related to patients’ prognosis, we applied survival analysis to estimate the contribution of 6 features for each of the 155 genes (930 total features) on survival time. We identified six significant oncogenic risk factors and 11 significant tumor suppressor factors (S8 Table).

Interestingly, the impact of homozygous deletions of candidate genes on survival was not significant. We speculated that it might result from heterogeneity of the tumor samples. Although the high expression of PDPN did not have a particularly significant impact on poor prognosis (p-value = 7.80E-04, FDR = 0.12, Fig. 5). Cancer cells with high PDPN expression have higher malignant potential because of enhanced platelet aggregation, which promotes alteration of cell motility, metastasis and epithelial-mesenchymal transition[39]. Previous studies have shown that overexpression of PDPN in fibroblasts is significantly correlated with a poor prognosis in ovarian carcinoma[40].

We also noted that the overexpression of EphA2 was associated with a shorter patient survival time (Fig. 5). Increased expressions of Eph receptor tyrosine kinases have been implicated in tumor progression in a number of malignancies[41, 42]. It was also observed that abnormal expression of EphA2 could lead to survival of patients with ovarian cancer[43].

ICan achieved a better modularity

The cutoff is the key parameter of MCODE that influences the size of the module. To select a more appropriate cutoff, we chose 0.01, 0.02, 0.03, 0.04, and 0.05, respectively, for parameter optimization. We found that when the cutoff was 0.02 (Fig. 6), the number of nodes in each module tended to stabilize. Eventually, 133 modules were identified. MCODE does not consider the weight of edges; therefore, we calculated the score M for each module using formula (4), and ranked the modules by the scores.

(5)

where a represents the ID of the module, E(a) represents the number of edges in module a, and N represents the number of nodes in module a.

To further explore the biological functions of ICan, we also compared ICan with the Pearson correlation coefficient method with only a single data type (copy number, methylation or gene expression). The results showed that the weights in the same module were significantly lower than that of the CCA method (p-value = 2.2e-16, Fig. 7, S4 Table). The modules of ICan were tighter in structure.

More precisely, the average weight of M7 (Fig. 8) was 0.7210 (CCA), 0.1176 (gene expression), 0.1113 (copy number) and 0.2305 (DNA methylation). We noted that the average weight of ICan was the highest, achieving more than three times the weight of the single-level networks. For instance, CTNNB1 and RPA2 were not only co-altered at the level of the copy number, but also at the level of DNA methylation. With respect to a single data type, we uncovered fine correlations between genes. Recent research has also shown that alterations in the genome region of CTNNB1 and RPA2 were closely related to the occurrence of ovarian cancer[44, 45]. The M7 module involved 20 genes (NIPSNAP1, ACTN4, ACTB, PIPOX, ACTG2, ACADVL, KRT8, LMNA, KRT1, KRT18, HSPA9, HSPD1, SET, HSPA5, XRCC5, CTNNB1, GBAS, C1QBP, CDH1 and ANP32B) that were significantly enriched to GO terms, including regulation of programmed cell death (2.2E-4), negative regulation of apoptosis (6.5E-4) and adherens junction (6.7E-4) (S9 Table). Among these, CTNNB1, CDH1, XRCC5 are important ovarian cancer genes. They are mainly involved in tissue invasion, metastasis and proliferation of cancer.