Integrated multi-omics analysis of genomics, epigenomics, and transcriptomics in ovarian carcinoma

Comparison of CNVcor and METcor gene sets

A total of 3,990 CNVcor genes and 9,651 METcor genes were identified at a significance level of p < 0.01. In the z-value distribution, CNVcor gene correlations were significantly shifted to the right, while METcor gene correlations were significantly shifted to the left (Figure 1A) (D'Agostino test: p < 1e-5). These results indicate a positive correlation between CNVcor genes and gene expression and a negative correlation between METcor genes and gene expression. Due to the large number of genes in these two sets, only those genes significantly related to prognosis in each set (p < 0.05, 413 CNVcor genes and 103 METcor genes) were included in subsequent analyses. There was no significant overlap between the two gene sets (Figure 1B), indicating a possible lack of interaction between CNVcor and METcor genes. Further analysis of the genomic distribution of CNVcor and METcor genes revealed that CNVcor is more inclined to appear on chr14, chr11, chr12 and chr1 chromosomes (FDR< 0.05), while METcor is more inclined to appear on chr12 and chr15, chr14 and chr16 (p < 0.05) (Figures 1C, 1D, S1, Table 1) and consist primarily of protein-coding genes (Figure 1E). The MET loci were primarily located in the S-shore, N-shore, and island regions (Figure 1F).

Figure 1. Identification of DNA copy number-correlated (CNVcor) and DNA methylation-correlated (METcor) genes in ovarian carcinoma. (A) Correlation z-value distributions for CNVcor genes and METcor genes. Density distribution is shown on the y-axis; the dotted line represents the median CNV and MET correlation coefficient values. (B) Venn diagram of overlapping CNVcor genes and METcor genes. The green area represents the 413 CNVcor genes. The yellow area represents the 103 METcor genes. The overlapping region contains 46 genes. (C) Box plot of CNVcor gene chromosomal distribution (upper figure) and correlations (lower figure). In the top image, the proportion of CNVcor genes on each chromosome (adding to a total of 1) is shown on the y-axis; in the bottom image, the correlation coefficient of CNVcor genes on each chromosome is shown on the y-axis. (D) Box plot of METcor gene chromosomal distribution. The proportion of METcor genes on each chromosome is shown on the y-axis. (E) METcor gene types. The proportion of METcor genes by type is shown on the x-axis. F: Distribution of MET loci. The proportion of types of methylated loci is shown on the x-axis.

Molecular subtypes were identified based on CNVcor and METcor genes

The CNVcor and METcor gene sets were clustered using non-negative matrix factorization (NMF). The standard “brunet” was selected by NMF and then subjected to 50 iterations. The clustering number k was set at 2-10. The average contour width of the sharing membrane matrix was determined by NMF in the R package, and the minimal member number for each subtype was set at 10. The optimal clustering number was determined using cophenetic, dispersion, and silhouette. The optimal clustering number was two for CNVcor genes (Figure 2A, Figure S2) and three for METcor genes (Figure 2B, Figure S3). Because no significant differences were observed when the CNVcor genes were clustered into three subtypes rather than two, and because the optimal clustering number for METcor genes was three, the optimal clustering number for CNVcor genes was also set at three for subsequent analyses. There were significant differences in OS for both the CNVcor and METcor genes in the three subtypes (Figure 2C, 2D, Figure S4, Table S4). Considerable overlap was observed between the CNVcor and METcor genes clustered into subtype three (Figure 2E and 2F, p < 0.05).

Figure 2. Identification of ovarian carcinoma molecular subtypes based on CNVcor and METcor genes. (A) NMF clustering results for CNVcor genes. (B) NMF clustering results for METcor genes. (C) KM survival curve for CNVcor gene clustering subsets. Survival time is shown on the x-axis, and survival rate determined by log rank P test is shown on the y-axis. (D) KM survival curve for METcor gene clustering subsets. Survival time is shown on the x-axis and survival rate determined by log rank P test is shown on the y-axis. (E) Overlap between CNVcor and METcor gene clustering subsets. F: Overlap test for CNVcor and METcor gene clustering subsets. Blue circles represent the proportion of overlapping samples between two the clusters; significance was determined using the Kolmogorov-Smirnov test.

CNV, MET, and EXP data were integrated to divide samples into four categories

Because METcor genes were clustered into three subtypes, K= 2-3 was chosen when calculating lambda values in the multi-omics clustering analysis. The final lambda values obtained at K=2 and 3 (3 and 4 clustering subtypes) were 0.004950495, 0.391089109, and 0.836633663, respectively. A considerable difference in sample size distribution among the three clustering subtypes was noted at K=2. Consequently, K=3 and 4 clustering subtypes with sample sizes of 1, 65, 128, and 156 were selected. Because the minimal clustering number cannot be less than 10, the subtypes containing one and 65 samples were merged into a single subtype. The final three resulting subtypes were iC1 (66 samples), iC2 (128 samples), and iC3 (156 samples). The clustering results for the three datasets are presented in Figure 3A and 3B, and clustering information is shown for each sample in Table S5.

Figure 3. OV sample molecular typing results in multiple datasets. (A) Expression heatmap for CNVcor gene subsets identified by iCluster. (B) Expression heatmap for METcor gene subsets identified by iCluster. (C) OS and PFS curves for subtypes identified by iCluster. (D) Overlap between METcor and CNVcor gene subsets with iCluster gene subsets. Blue circles represent the proportion of overlapping samples between two clusters; significance was determined using the Kolmogorov-Smirnov test.

OS (Figure 3C, p < 0.001), but not PFS (Figure 3C, p>0.05), differed significantly among the 3 subtypes. In addition, there was considerable overlap between the iCluster gene clustering results and the CNVcor and METcor gene clustering results (Figure 3C, 3D, p < 0.001).

Correlations between CNV and MET variation

Next, we examined the relationship between CNV and MET variation. For CNVs, CNV Gain was defined by a β>0.3 and CNV Loss by a β<-0.3; for MET, hypermethylation (MetHyper) was defined by a β>0.8 and hypomethylation (MetHypo) by a β<0.2. The number of genes showing CNV Gain, CNV Loss, MetHyper, and MetHypo was counted separately for each sample. CNV Gain and CNV Loss were positively correlated (R=0.15, p=0.0037) (Figure 4A), but there was no correlation between CNV Gain and CNV MetHyper (Figure 4B). There was also a strong positive correlation between CNV Gain and MetHypo (R=0.29, p=3e-08) (Figure 4C); CNV Loss and MetHyper were positively correlated as well (R=0.2, p=0.00022) (Figure 4D). There was no significant correlation between CNV Loss and MetHypo (Figure 4E). Finally, MetHyper and MetHypo were negatively correlated (R=-0.29, p=3.1e-08) (Figure 4F, Table S6).

Figure 4. Associations between aberrations in DNA copy number and DNA methylation in ovarian carcinoma. (A) Frequency distribution of CNV Gain and CNV Loss. (B) Frequency distribution of CNV Gain and MetHyper. (C) Frequency distribution of CNV Gain and MetHypo. (D) Frequency distribution of CNV Loss and MetHyper. (E) Frequency distribution of CNV Loss and MetHypo. (F) Frequency distribution of MetHyper and MetHypo. Larger correlation coefficients (R values) indicate stronger correlations; the log rank P test was used.

Immune scores differ among Ovarian Cancer (OV) gene subtypes

All OV genes included in this study were clustered into three subtypes based on available multi-omics data. No significant differences in clinical characteristics (gender and age) were found among the three subtypes (Table 2). Immune scores were determined for samples grouped by subtype using tumor immune estimation resource (TIMER) tool. Compared to the other two subtypes, B cell scores were lower, while macrophage cell scores were higher, in iC1 subtype samples, which were also associated with the poorest prognoses (Figure 5A, p < 0.01). The immune scores of iC1 subtypes in CD4 T cells, neutrophils and dendritic cells are significantly lower than that of iC3 (P<0.05). We speculate that these molecular subtypes have different cellular immune functions and affect the survival and prognosis of patients to a certain extent.

Figure 5. Infiltration of iC molecular subtypes in different immune cells. (A) Comparison of scores for six types of immune cells among the three iCluster gene subtypes. Immune cell scores are shown on the y-axis, and different molecular subtypes are shown on the x-axis. The Kruskal-Wallis rank test was used to identify significant differences in immune cell scores between the iC molecular subtypes. (B) Scores for six types of immune cells in all samples. iCluster indicates multi-omics molecular subtypes, METcor_C indicates METcor gene molecular subtypes, and CNVcor_C indicates CNVcor gene molecular subtypes.

Scores for the six immune cell types examined in the three sample subtypes are shown in Figure 5B; the associated data are shown in Table S7.

Molecular characteristics of OV gene subtypes

Differences in CNV, MET, and mRNA expression between the iCluster iC1 and iC3 subtypes, which also differed significantly in prognosis, were then analyzed (Figure S5). Samples were separated into the following three types based on CNV and MET data as previously described: CNV Gain (MetHyper), CNV Loss (MetHypo), and CNV Normal (MET Normal). Genes that differed significantly in CNV and MET data between the iC1 and iC3 subtypes were identified using the Fisher-exact test. Genes with differential expression based on gene expression (EXP) data between iC1 and iC3 subtypes were identified by DESeq2 (differential genes with p < 0.05). The genes that differed in CNV, MET, and expression spectrum are shown in Figure 6.

Figure 6. The 100 genes with the most significant differences in CNV, Met, and gene expression among iC1 molecular subtypes. (A) CNV distribution in iCluster gene subtypes; (B) MET distribution in iCluster gene subtypes; (C) Heatmap of differential genes for iCluster gene subtypes.

Associations between gene expression, CNV, and MET

To investigate the relationship between gene expression, CNV, and MET, genes with both differential expression in the iC1 and iC3 subtypes and with a significant difference in CNV Gain/Loss and MetHypo/MetHyper (p < 0.05) were examined in a prognostic survival analysis. A total of 1,338 genes differing in CNV Gain/Loss, 11 genes differing in MetHyper/MetHypo, and 8,195 genes with differential expression between the iC1 and iC3 subtypes were included. Of the seven genes which differed in all three measures between the two subtypes (URI1, AKT2, ZHX3, RAB34, FBXO6, IL18BP, and UBB), two genes (UBB and IL18BP) were significantly associated with prognosis according to one-factor survival analysis (Figure 7. Table S8 and Table S9, p < 0.05). Specifically, low expression of both UBB and IL18BP was associated with poorer prognosis. Furthermore, in samples associated with poor prognosis, expression of these two genes was lower in iC1 subtype samples than in iC3 subtype samples, perhaps due to higher MET levels and lower CNV levels in the iC1 subtype.

Figure 7. CNV and MET in UBB and IL18BP were significantly correlated with prognosis. (A) Comparison of CNV, MET, and EXP in the UBB gene between iC1 and iC3 subtypes and prognostic survival curve for the TIMM50 gene. From left to right: the proportion of UBB Gain and Loss (range: 0-1) in iC1 and iC3 samples, the proportion of UBB hypermethylation and hypomethylation (range: 0-1) in iC1 and iC3 samples, UBB expression in iC1 and iC3 samples, and survival curve with samples divided into high and low UBB groups based on median gene expression level. (B) Comparison of CNV, MET, and EXP in the IL18BP gene between iC1 and iC3 subtypes and prognostic survival curve for the TIMM50 gene. From left to right: the proportion of IL18BP Gain and Loss (range: 0-1) in iC1 and iC3 samples, the proportion of IL18BP hypermethylation and demethylation (range:0-1) in iC1 and iC3 samples, IL18BP expression in iC1 and iC3 samples, and survival curve with samples divided into high and low IL18BP groups based on median gene expression level. Survival time is shown on the x-axis and overall survival is shown on the y-axis; differences were identified using the log rank P test.

Survival analysis for UBB and IL18BP

Associations between the two genes of interest and prognosis were validated using 10 independent datasets from the KMplot website www.kmplot.com/lung. KM survival curves for OS, PFS (p < 0.05), and post-progression survival (PPS) for these genes are shown in Figure 8.

Figure 8. UBB and IL18BP gene survival curves in the KMplot dataset. (A) OS, PFS, and PPS curves for the UBB gene in a dataset from the KMplot database. Survival time is plotted on the x-axis; OS, PFS, and PPS are plotted on the y-axis. (B) OS, PFS, and PPS curves for the IL18BP gene in the KMplot dataset. Survival times are plotted on the x-axes; OS, PFS, and PPS are plotted on the y-axes. Differences were identified using the log rank P test.

Comparison of mutation spectrum in the OV gene subtypes

Based on the iCluster gene clustering results, we analyzed mutation spectra in the different subtypes and screened a set of genes that differed significantly between the iC1 and iC3 subtypes (Figure 9). In total, 83 genes with p < 0.05 were selected based on ranked Fisher test p values; associated data are shown in Table S10, and Fisher test results for SNV loci by group are presented in Table S11.

Figure 9. Genes with differential mutations in the iC1 and iC3 molecular subtypes. Different colors represent different numbers of mutations in the corresponding genes. iCluster indicates multi-omics molecular subtype, METcor_C indicates METcor gene molecular subtype, and CNVcor_C indicates CNVcor gene molecular subtype.

To investigate mutations in and expression of the UBB and IL18BP genes, we calculated the correlation between expression of these two genes at each SNV locus; the top ten loci are listed in Tables 3 and 4, and all results are shown in Tables S12 and S13. It was inferred from these results that the mutation loci significantly related to expression levels might have specific effects on UBB and IL18BP expression.

Download of TCGA data

The most recent clinical follow-up data were obtained from the TCGA Genetic Disease Control (GDC) API on January 24, 2019; CNV, MET, and RNA-seq (including read count) data were also obtained for subsequent analysis of differential gene expression in different patient subsets. In addition, SNV data (mutect version) were downloaded from TCGA. Data from 351 patients in three datasets were included in the analysis; sample information for all three datasets is shown in Table S1.

Profiling of DNA copy numbers, DNA methylation, mRNA expression, and SNV data

The CNV data were pre-processed as follows. Two regions with 50% overlap were considered identical. Regions covering <5 probes were deleted. The CNV region was mapped to corresponding genes using the GRCh38 release 22 (https://www.gencodegenes.org/human/release_22.html). Multiple CNV regions in a gene were merged into a single region, and CNV values were averaged to provide a merged CNV value. MET data were pre-processed by deleing absent loci in >70% of samples. Missing data were imputed using the KNN (k-Nearest Neighbor) algorithm. Probes in the TSS region from 2kb upstream to 200bp downstream were preserved using GRCh38 release 22 and mapped to the corresponding genes. RNA-seq data were pre-processed by deleting genes with low expression levels (FPKM = 0 in <0.5% of all samples). SNV data were pre-processed by deleting mutations in intron regions and silent mutations.

Identification of CNVcor and METcor gene sets

The Pearson correlation coefficients for associations between CNV and RNA-seq and between MET and RNA-seq were calculated separately and converted into z-values using the formula ln((1+r)/(1-r)). Genes with p < 0.05 in the correlation coefficient test were included in the CNVcor and METcor gene sets. CNVcor and METcor gene data are shown in Table S2 and S3, respectively.

Sample clustering via integration of CNV, MET, and gene expression data (EXP) data

The “iCluster” R package was used to conduct multi-omics clustering analysis by integrating CNV data from CNVcor genes, MET data from METcor genes, and EXP data from both CNVcor and METcor genes. Optimal weights for CNV, MET, and EXP datasets were determined based on lambda values. After completing 20 iterations to optimize lambda values, a total of 101 lambda sample points valued 0-1 were selected.

Survival analysis

The KMplot website (http://kmplot.com/analysis/) was used to validate the data [30]. This database system contains integrated data from 8 independent datasets consisting of a total of 1,657 TCGA Ovarian Cancer (TCGA-OV) samples. Ovarian cancer patients were divided into 2 groups (high and low expression) based on the expression of the gene of interest. Kaplan-Meier plots were used to analyze overall survival in ovarian cancer patients. Hazard ratios (HR), 95% confidence intervals (CIs), and log rank P-values were evaluated.