Integrative analysis of genomic and epigenetic regulation of endometrial cancer

Molecular subtypes based on CNV-G and MET-G genes

Next, we determined whether the expression profiles of CNV-G and MET-G genes could predict prognosis. Cluster analysis showed an optimal clustering number of 4 for both CNV-G and MET-G (Figure 2A and 2B). Kaplan-Meier (KM) plots, for overall survival (OS), revealed significant differences in prognostic outcomes between the groups (Figure 2C), with marked differences observed in the MET-G subclass (Figure 2D). In addition, there was a significant overlap among the four subclasses of both CNV-G and MET-G (Figure 2E, 2F).

Figure 2. Molecular subtypes based on CNV-G and MET-G genes. (A) NMF-based clustering of CNV-G. (B) NMF-based clustering of MET-G. (C) KM survival curve of CNV-G subtype. (D) KM survival curve of MET-G subtype. (E) The overlap between the CNV-G subtype and the MET-G subtype. (F) The overlap between the CNV-G subtype, the MET-G subtype and the histological subtype.

CNV, MET, and EXP datasets were integrated for cluster analysis

We repeated the clustering 20 times in K=2 (class 3) and K=3 (class 4) in order to optimize clusters created by iCluster. The results indicated a more stable clustering in rank=2 (class 3) than rank=3 (class 4). We finally concluded that the iCluster was aggregated into three subclasses: iC1 (95 samples), iC2 (128 samples), and iC3 (198 samples). Based on the CNV level distribution of iCluster, we found a higher CNV in iC2 compared to iC1 and iC3 (Figure 3A), although their methylation levels were comparable between the groups (Figure 3B). A further comparison of iCluster and existing EC subtypes revealed that iC2 mainly corresponded to CNV high subtypes, whereas CNV low and MSI subtype samples were mainly concentrated in iC3 (Figure 3A). KM survival analysis indicated significant differences in OS between the three groups (Figure 3C). Further comparisons among the three groups revealed significantly different prognosis among iC1, iC3 and iC2 subtypes (Figure 3D, 3E). However, there was no significant difference in prognosis between iC1 and iC3 subtypes (Figure 3F). Moreover, progression-free survival (PFS) was significantly different among the three subtypes (Supplementary Figure 1).

Figure 3. (A) CNV levels of subtype CNV-G identified by iCluster. (B) Methylation levels of MET-G subtype identified by iCluster. (C) KM curve for the subtypes identified by iCluster. (D) KM curve for iC1 and iC2 subtypes. (E) KM curve for iC2 and iC3 subtypes. (F) KM curve for the iC1 and iC3 subtypes.

Abnormality in DNA copy number is consistent with methylation abnormality

To study the relationship between CNV and MET abnormalities, we defined the β value of CNV > 0.3 as CNV Gain; β value < -0.3 as Loss; β value of MET > 0.8 as MetHyper (hypermethylation); and β value < 0.2 denoted as MetHypo (demethylation). We counted the numbers of CNV Gain, Loss, MetHyper and MetHypo for each sample and found a significant correlation between Gain, Loss and MetHypo (Figure 4A, 4C), not with MetHyper (Figure 4B). Additionally, we did not record a significant correlation between Loss and MetHyper (Figure 4D), although both were significantly correlated with MetHypo (Figure 4E). MetHyper and MetHypo showed a strong negative correlation (Figure 4F).

Figure 4. (A) Frequency distribution of CNV Gain and 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.

Characteristics of the tumor microenvironment among the molecular subtypes

We categorized EC into three subgroups according to multi-group data, then compared the differences in clinical characteristics of iC subtypes in stage, grade, age, and BMI (body mass index). Results revealed significant differences in the distribution of iC subtypes among samples with different clinical characteristics (Table 1). High-grade and advanced samples were more likely to be distributed in the worst prognostic iC2 subtype (Supplementary Figure 2). We used tumor immune estimation resource (TIMER) to compare immune scores across the three subtypes, and found that six immune cell scores were lower in iC2 subtype and had the worst prognosis compared to the other subtypes (Figure 5A, 5B), indicating that the iC2 subtype may represent an immunosuppressive state. Comparative analysis further showed a significantly lower macrophage regulation and lymphocyte infiltration score in the iC2 subtype relative to the other subtypes, whereas the wound healing and inflammation (IFN-gamma response)-related score was significantly higher in iC2 than other subtypes (Figure 5C). This further suggested that immune status may affect the prognosis of EC.

Figure 5. (A) Immune cell scores obtained from all samples. (B) A comparison of all immune cell scores among the three subtypes of iCluster. (C) A comparison of 5 immune signatures scores.

Molecular characteristics of the subtypes

Based on the results from iCluster, we compared the differentially expressed genes (DEGs) between iC1/iC3 and iC2. We identified a total of 207 DEGs in the three groups, after removing the low expression levels. Gene ontology (GO) analysis indicated a significant enrichment of terms related to immune regulation, such as leukocyte migration and adaptive immune response. The CNV frequency of 207 DEGs in iC2 was significantly higher than that in iC1 and iC3, suggesting that CNV influenced the prognosis of EC (mean CNV: 8005/28579/5899) (Figure 6A). However, no significant differences were observed in methylation levels between the molecular subtypes (mean methylation: 134682/143148/140185) (Figure 6B). A correlation between expression level, methylation and CNV, revealed a high expression of DEGs in demethylated samples (Figure 6C), but this was not observed in CNV. This indicated that the effect of methylation on the expression of DEGs was stronger than the effect of CNV on the expression of DEGs. Univariate survival analysis identified 24 genes that were significantly associated with prognosis. We also analyzed expression of three genes: KIAA1324, nonexpresser of pathogenesis-related genes1 (NPR1) and idiopathic hypogonadotropic hypogonadism (IHH), and the relationship between methylation and CNV. We also found a significant negative correlation between expression of these genes and their methylation status, but not in the CNV (Figure 7). Studies have shown that KIAA1324 is activated by estrogenase, which suggests that estrogenase may play a role in the occurrence of EC.

Figure 6. (A) Distribution pattern for CNV in iCluster. (B) Distribution for methylation level in iCluster. (C) Heatmap of differentially expressed genes in iCluster subtypes.

Figure 7. (A–C) A correlation of methylation of KIAA1324 gene with its expression, expression of iC subtypes, and the KM curve of the high/low expression groups. (D–F) A correlation of methylation of NPR1 gene with its expression, expression of iC subtypes, and the KM curve of the high/low expression groups. (G–I) The relationship between methylation of IHH gene with its expression, expression of iC subtypes, and the KM curve of the high/low expression group.

Mutation spectrum of molecular subtypes

We further mapped the mutation spectrum of various molecular subtypes to identify differentially expressed genes in iC subtypes. Using the Fisher’s test (with FDR < 0.001), we obtained a total of 48 genes. Mutation spectrum analysis showed a significantly lower mutation frequency of PTEN, ARID1A, CTNNB1 in the iC1/iC3 subtype with better prognosis than that in the iC2 subtype (FDR < 0.0001). Moreover, the mutation frequency was significantly higher in subtype iC2 than in iC1/iC3 (Figure 8A). Overall, there were fewer silent/nonsilent mutations and neoantigens in iC2 than in iC1/iC3 (Figure 8B). However, the number of CNVs in iC2 was significantly higher than that of iC1/iC3 (Figure 8C), suggesting that the effect of gene copy number variation on prognosis was stronger than that of genomic mutations. Additionally, methylated MetHyper/MetHypo levels varied significantly among the molecular subtypes (Figure 8D).

Figure 8. (A) Profiles of significant mutations in 48 genes across iC subtypes. (B) Distribution pattern of mutation number of the 48 genes with significant mutations. (C) Distribution of silent/nonsilent and neoantigens among iC subtypes. (D) Distribution of CNV Gain/Loss and methylated MetHyper/MetHypo among iC subtypes.

Data collection

We downloaded recent clinical follow-up information from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/) using the TCGA GDC API (https://gdc.cancer.gov/developers/gdc-application-programming-interface-api) at 2019.01.24. CNV, Methylation, RNA-seq data and the SNP data processed by the mutect software were also downloaded. The RNA-seq data also includes the UCEC count for subsequent group differential expression analysis. A total of 161 samples of all three sets of data used in subsequent analyses.

Data preprocessing

The following processing procedures were performed on the CNV, Methylation, RNA-seq, and SNV datasets in the EC samples from the TCGA database:

CNV data preprocessing

The CNV intervals were merged using the criteria shown below:

1. 50% regional overlap between two intervals was considered as the same interval.

2. The number of coverage probes <5 intervals were removed.

3. The CNV interval was mapped to the corresponding gene using the GRh38 version of 22.

4. Multiple CNV regions in one gene region were combined into one, and the combined CNV value was averaged.

Preprocessing of methylation data

1. Missing sites in more than 70% of the samples were removed.

2. KNN (k-nearest Neighbour) algorithm was applied to fill in missing values.

3. Probes 2 kb upstream and 200 bp downstream of the TSS intervals were retained by the annotated version of gencode.v22 and mapped to the corresponding genes.

Preprocessing RNA-seq data

1. Lowly expressed genes (samples with FPKM of 0 accounting for <0.5 of all samples) were removed.

SNV data preprocessing

1. Mutations in the intron interval were removed.

2. Mutations annotated as silence were removed.

Identification of CNV-G and MET-G

We calculated Spearman correlation coefficients for each gene corresponding to CNV and expression profile (RNA-seq) and methylation and expression profile, respectively. These coefficients were then converted into z-value using the formula ln(1+r)/(1-r). Genes with p < 1e-5, which were tested for correlation coefficients, constituted a CNV-G and a MET-G.

Identification of molecular subtypes of CNV-G and the MET-G

Nonnegative matrix factorization (NMF) is an unsupervised clustering method widely used to identify molecular subtypes of tumors based on genomics [19, 20]. To further explore the association between CNV-G and MET-G expressions with clinical phenotypes, samples were clustered using the NMF method based on the expression profiles of the CNV-G and MET-G gene sets. Briefly, we selected the standard "brunet", using the NMF method 50 iterations. The number of clusters k was set between 2 to 10, and the average profile width of the common member matrix calculated using the NMF package implemented in R [21], with the minimum member of each subclass was set to 10.

Identification of molecular subtypes

We employed the R package ‘iCluster’ [13], to perform multi-group data integration cluster analyses and integrate the copy number variation (CNV) data of the CNV-G gene, methylation data (MET) of the MET-G gene, as well as the expression profile data (EXP) of the genes in CNV-G and MET-G. Subsequently, 20 iterations and10 lambda sample points between 0-1 were used for optimal lambda value screening to identity optimal CNV, MET, and EXP data weight values (lambda values). Considering the number of molecular subtypes identified by CNV-G and MET-G, we chose 2-4 as the number of clustering K.

Relationship between molecular subtypes and tumor microenvironment

Tumor immune estimation resource (TIMER) is a platform used for systematic assessment of the clinical impact of different immune cells in various types of cancer [22]. This method was used to estimate the abundance of six immune cell types, namely: B, CD4 T, CD8 T, and neutral cells, macrophages and dendritic cells. The abundance of those cells in the tumor microenvironments were analyzed in different molecular subtypes.

Analysis of genetic differences in molecular subtypes

To examine differences in gene expression among the molecular subtypes, we employed DESeq2 [23] tool using 2-fold differences and FDR < 0.05 as thresholds for identifying differentially expressed genes between molecular subtypes.

Relationship between molecular subtypes and tumor genomic variation

To assess the differences in genomic variation between molecular subtypes, we downloaded SNP data from TCGA, then removed introns and silent mutations. We then used the Fisher's exact test to compare differences in mutations between the two samples. A threshold of p<0.05 was used to identify mutated genes.

Functional enrichment analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using the R package cluster profiler 28 for genes. We then identified over-represented GO terms in three categories namely; biological process, molecular function, and cellular component. For these analyses, a FDR < 0.05 was considered for the determination of statistical significance.

Statistical analysis

Kaplan-Meier was used to visualize the differences in subtype prognosis, while univariate survival analysis was performed to estimate overall survival. The log-rank test was used to test prognostic differences at a significance of p < 0.05. All the analyses were performed in R software version 3.4.3.