Integrated bioinformatics data analysis reveals a risk signature and PKD1 induced progression in endometrial cancer patients with postmenopausal status

Data collection

The Cancer Genome Atlas (TCGA) UCSC XENA (https://xena.ucsc.edu/) and the Gene Expression Omnibus (GEO) database were used to obtain clinical data and gene expression profiles. Affymetrix Human Gene 2.0 ST Array was used to process the GSE17025 dataset. The level of TCGA gene expression was measured as Transcripts per million reads (TPM). The data of this study are from GEO and TCGA public databases.

Identifications of differentially expressed genes (DEGs)

The GSE17025 expression profile was normalized and analyzed using R software and the limma package, while the TCGA EC dataset was normalized and analyzed using R software and the edge R package. Cut-off criteria were considered the |log₂ Fold Change (log₂FC)| >1 and adjusted p-value <0.05.

Construction of WGCNA and module preservation

We applied the WGCNA to construct the gene co-expression network and identify the co-expression modules in R software [7]. In this study, we selected the minimum size (genome) 30 for the gene tree, the cutting line of 0.25 for the module tree, and combined some modules. We assessed how similar the module eigengenes (MEs) were, identified a cut line for the module dendrogram, and merged certain modules to better understand the module.

Identification of clinically significant modules

Real hub genes were selected by drawing a Venn diagram which combining the module, DEGs in GEO database (p < 0.05, log₂FC>1), and DEGs in TCGA database (p < 0.05, log₂FC>1) together. The expression similarity of different samples is used to identify the WGCNA, while the relationship between the external clinicopathological information and the gene modules is used to identify clinically significant modules. Finally, gene modules that were highly correlated with certain clinical features were chosen as modules of interest for further analysis. Functional and pathway enrichment analysis of significant module was conducted according to the methods previously reported [13]. P < 0.05 was set as the cut off criterion.

Construction and evaluation of the postmenopausal related prognostic model

LASSO regression analysis was performed to select essential prognostic genes by the “glmnet” R package. We listed the formula of the risk score for the predicting patients’ survival as follow:

$$\text{Risk score}=\sum {\text{n}}_{\text{i}}=\sum (\text{Coefi}\times {\text{x}}_{\text{i}})$$

Coefi represented the coefficient and x_i represented the expression level of hub genes. Subsequently, risk scores were calculated by the formula and gene expression for each patient. Then patients were divided into two groups (high- and low-risk) according to the median value of risk score. In these two subgroups, the clinicopathological features and gene expression profiles of each patient were displayed by "pheatmap" and "survival" R packages. The Kaplan-Meier survival analysis was conducted to compare the overall survival (OS) rate of the two subgroups by “survival” package of R (P < 0.05). The accuracy of the risk model was further evaluated by receiver operating characteristic (ROC) curve.

We validated possible uses of the predictive risk model and performed a survival analysis within our cohorts, and we obtained RNA-seq expression and clinical data from 30 patients that underwent surgery in the Department of Obstetrics and Gynecology, Nanjing First Hospital. Thirty samples without neoadjuvant therapy and who underwent surgical resection between January 2019 and December 2021 from patients in our center were selected. RNA isolation and reverse transcription-quantitative PCR were both conducted according to previously described methods [13]. The Institutional Ethics Committee (Human Research) of our hospital approved our research, and we obtained informed consent from all participants.

In vitro validation of hub genes

Total protein of different types of tissues were extracted and used for western blot as described earlier [6]. Total RNA was extracted from tissues using TRIzol reagent (Tiangen, China). We used the PrimeScript RT-polymerase (Vazyme) to reverse-transcribe cDNAs to the mRNA of interest. Real-time quantitative RT-PCR (qRT-PCR) was conducted using SYBR-Green Premix (Vazyme) with specific PCR primers (Thermo Fisher Scientific, USA). GAPDH was taken as an internal control. The relative expression of the target gene was calculated using the 2^−ΔΔCt method. qRT-PCR was performed according to the manufacturer’s instructions.

Statistical analysis

Statistical differences in the expression of hub genes in the normal and tumor samples were analyzed using a two-tailed student’s t-test in both the TCGA and GEO databases. All statistical tests and graphing were performed using R software and GraphPad prism 7.0. All analyses were conducted three times and represented data form three separate experiments. P value <0.05 was considered statistically significant.

Identification of DEGs in GEO and TCGA databases

A brief study design was shown in Figure 1. The DEGs in both of the GEO and TCGA datasets were analyzed respectively using p < 0.05 and |log₂FC| >1 as the cutoff criteria. The volcano maps of the DEGs in the two groups were constructed using the R package. Based on the screening guidelines, 4163 DEGs were acquired, consisting of 1900 up-regulated genes and 2263 down-regulated genes from EC tissues compared with normal endometrial tissues in TCGA dataset (Figure 2A, 2B). What’s more, we found a total of 2461 DEGs in the normal endometrial tissue in GSE17025, including 1487 downregulated and 974 upregulated genes (Figure 2C, 2D). In total, 1240 common genes are acquired by Venn diagram (Figure 2E).

Figure 1. The flowchart of the study design.

Figure 2. Screening of common differentially expressed genes (DEGs) of TCGA and GEO databases. (A) Volcano plot of DEGs between EC and normal endometrial samples in TCGA database. (B) Bar plot for DEGs of dysregulated genes in TCGA database. Red bar, up-regulated mRNA, blue bar, down-regulated mRNA. (C) Volcano plot of DEGs between EC and normal endometrial samples in GSE17025. (D) Bar plot for DEGs of dysregulated genes in GSE17025. (E) Venn diagram of DEGs between TCGA and GEO databases, the blue circle represents for total number of DEGs in TCGA, and the pink circle represents for total number of DEGs in GSE17025, the overlapped part is used for further analysis.

WGCNA and identification of key modules

We constructed co-expression network by TCGA dataset including 545 EC samples associated with whole clinical information. When constructing co-expression network with “WGCNA” software package, the expression value of the first 25% DEG is included. In our study, β = 4 (scale free R² = 0.89) was selected as the soft-thresholding power to ensure a scale-free network (Figure 3A–3C). All selected genes were classified using the dissimilarity measure based on topological overlap matrix (TOM), which is based on the dynamic tree cutting algorithm and divides the tree into eight modules marked with different colors (Figure 3D). We identified the relationship between each module and the EC clinical information, with a particular emphasis on postmenopausal status (Figure 3E). Pearson’s correlation coefficient was used to assess the association between the co-expression modules in this study. Modules with various gene clusters were labeled in the topological overlap heatmap using different colors, with blue indicating a negative relationship and red indicating a positive relationship (Figure 3F). Then 1000 genes were randomly selected for the heatmap. Results of the clustering analysis indicated that the hierarchical clustering of the module eigengenes was representative of the modules. The dendrogram branches were clustered according to the eigengene relationships (Figure 3G).

Figure 3. Weighted gene correlation network (WGCNA) on the RNA-seq database and selection of hub genes. (A, B) Screening and validation of the soft threshold. (C) Checking the scale free topology when β = 4. The x-axis demonstrates the logarithm of whole network connectivity, while the y-axis shows the logarithm of the corresponding frequency distribution. (D) Clustering dendrogram of common DEGs in EC tissues. (E) Correlation between modules and postmenopausal status. (F) Correlations between different modules. (G) Heatmap depicts the Topological Overlap Matrix (TOM) of genes selected for WGCNA. Light color represents lower overlap and deep color represents higher overlap.

First, modules with a greater MS were considered to have more connection with patients in postmenopausal status. However, most of the correlations were between low and moderate (R² < 0.5), and the MS of the brown module (R² = 0.83, P = 4e-18) was found to be higher than that of the rest of the modules. Therefore, the brown module with postmenopausal status was identified as the clinically significant module, which was extracted for further analysis.

GO and pathway enrichment analysis of hub genes in brown module

GO and KEGG enrichment analysis were used to explore the function and pathways of the hub genes. The brown module, which consisted of 471 genes, was mainly associated with the following subclasses of GO classification (Figure 4A): calcium ion binding, translational initiation, integral component of plasma membrane, and cell adhesion. KEGG pathway analysis showed that top enriched terms were P53 signaling pathway, cell adhesion molecules, T cell receptor signaling pathway, and chemical carcinogenesis (Figure 4B). These suggested that brown module genes were closely related to the calcium channel and cell adhesion.

Figure 4. Functional analysis of common DEGs from TCGA and GSE17025 datasets. (A) GO analysis showing the differentially expressed postmenopausal related genes. (B) The significantly enriched pathways of the DEGs determined by KEGG analysis. Abbreviations: GO: gene ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes.

Establishment and validation of a postmenopausal related prognostic signature

For screening out the potentially prognostic biomarkers for patients with postmenopausal status, LASSO regression was conducted. LASSO analyses identified 11 genes significantly associated with prognosis: PKD1, PCDHB15, PRND, TBC1D8B, ZNF578, SLC26A4, CACNA1D, CPEB1, CYP46A1, EPPK1, and KBTBD12 (Figure 5A, 5B). According to the results of LASSO regression analysis, we used the coefficients (Table 1) to construct the prognostic model as following: risk score = (PKD1 × 0.13192805) − (PCDHB15 × 0.082172079) − (PRND × 0.026889625) − (TBC1D8B × 0.093931462) − (ZNF578 × 0.048951829) + (SLC26A4 × 0.209421425) + (CACNA1D × 0.372518929) + (CPEB1 × 0.045936749) + (CYP46A1 × 0.039213661) + (EPPK1 × 0.019827739) − (KBTBD12 × 0.025230047). The patients were then ranked, in ascending order, according to the parameter, after which they were classified into low-risk and high-risk groups based on the median risk values. We further analyzed the relationship between the 11 genes (Figure 5C). We found that they were significantly relevant, especially between PKD1 and ZNF578, PKD1 and SLC26A4, CACNA1D and ZNF578, SLC26A4 and PRND. In addition, the expression of 11 genes in TCGA dataset in low-risk and high-risk patients was also confirmed in the heatmap. The expression profiles of the prognostic genes showed that PKD1, SLC26A4, CACNA1D, CPEB1, CYP46A1, and EPPK1 were highly expressed in the high-risk subgroup. We observed significant differences in the high- and low-risk groups related to peritoneal cytology, histology, menopause, tumor LNM, recurrence, grade, age, stage, and death status (Figure 5D). The dot plot displays the ranked risk score and survival status of each individual, with notable differences between the groups in terms of overall survival (OS) (Figure 5E, 5F). Kaplan-Meier curve analysis indicated that the OS of the high-risk group was significantly shorter than that of the low-risk group (P = 1.827e−05, Figure 5G). Analysis of the ROC curve indicated that the area under the ROC curve (AUC) of the prognostic HRGs model for OS was 0.729 (Figure 5H). These results demonstrated that the menopause-related risk signature had a high precision in the prediction of EC patients.

Figure 5. Identification of postmenopausal related prognostic signature in EC patients. (A) Plots of the ten-fold cross-validation error rates. (B) LASSO coefficient profiles of the eleven postmenopausal related genes. (C) Spearman correlation analysis with the selected 11 genes. (D) Heatmap and clinicopathological features of high- and low-risk groups. The samples are ordered by risk score, and the score decreases from left to right. (E, F) Risk score distribution in low- and high-risk groups. (G) Kaplan-Meier survival analysis of the low and high-risk group. (H) Time-dependent ROC curves for overall survival prediction for EC patients.

The risk signature mentioned above was further validated in the clinical cohort in our center (Supplementary Table 1). We conducted the RNA level sequencing of 30 EC patients and calculated the risk score of each patient with the formula mentioned above. EC patients in the cohort were divided into low- and high-risk groups based on the median risk score as before. Expression of PKD1 was compared between patients with pre- and postmenopausal status (Supplementary Figure 1A). Kaplan-Meier survival curves showed that patients with low risk scores had a longer OS, which was consistent with the predicted survival results of the risk model (p = 0.031) (Supplementary Figure 1B).

Verification of expression of PKD1 in database and in vitro

We further validated the expression of PKD1 in different clinicopathological characteristics of EC patients in TCGA database. The results indicated that PKD1 was highly expressed in cancer tissue, grade 3, positive lymph node metastasis, positive peritoneal cytology, and recurrence groups (Figure 6A–6E). Patients were divided into two cohorts according to median value of PKD1. The results indicated that the expression of PKD1 was significantly associated with patient prognosis (Figure 6F). The analysis showed that PKD1 was statistically different in worse outcomes of EC patients. GSEA results showed that high expression of PKD1 was enriched in apoptosis, and low expression of PKD1 was enriched in B cell receptor signaling pathway, endometrial cancer, mTOR signaling pathway, and VEGF signaling pathway (Supplementary Figure 2). The immune infiltration was significant diverse in different expression of PKD1 (Supplementary Figure 3). To further explore the expression of PKD1 in different menopause status, we performed the qPCR and western blot (WB) validation in clinical specimens following the steps described above. As shown in Figure 7A, 7B, the expression levels of PKD1 in postmenopausal patients was more than that in the premenopausal patients in EC group. These results indicated that PKD1 may play a pivotal role in EC especially in postmenopausal women.

Figure 6. Expression and survival validation of PKD1 in TCGA database. Expression of PKD1 in different clinicopathological features. (A) Normal and cancer tissues. (B) Different grade. (C) Positive and negative lymph node metastasis. (D) Positive and negative peritoneal cytology. (E) Different recurrence status. (F) Kaplan-Meier survival plot of low and high expression of PKD1 according to the median value.

Figure 7. Relative expression of PKD1 in 5 pairs of pre- and post-menopausal patients in our center for (A) protein level and (B) RNA level. Pre, premenopausal status, Post, postmenopausal status.