Introduction
Global cancer statistics released in 2020 indicated that endometrial cancer has become the sixth most common gynecological tumor, with an estimated 417,367 new cases and 97,370 deaths worldwide.1 The American Cancer Society also reported an increasing incidence and mortality trend for endometrial cancer in recent years, which underscores the urgent demand to develop more effective strategies for diagnosis and treatment.2 As endometrial cancer is frequently symptomatic and detectable at an early stage, such as abnormal vaginal bleeding, most cases are diagnosed at Stage I and have favorable outcomes after prompt surgically curable.3,4 However, worse clinical outcomes often occurred in patients with advanced-stage, recurrence, or metastatic endometrial cancers, partly direct to the abnormal molecular features of the highly aggressive tumor cells gained in the disease progress.5,6
Gremlin-1 (GREM1), a bone morphogenic protein (BMP) antagonist, has been known to induce carcinogenesis through multiple mechanisms, such as BMP dependent pathway or vascular endothelial growth factor receptor (VEGFR) signaling pathway.7 Previous studies have indicated elevated expression levels and carcinogenic roles of GREM1 in various malignancies, such as carcinomas of lung, ovary, kidney.8 Recently, GREM1 was identified as a critical regulator of cellular heterogeneity in pancreatic cancer whose over-expression could restrict epithelial-to-mesenchymal transition by interacting with BMP, inhibiting tumor metastasis.9 Cheng et al. showed that GREM1 promotes lineage plasticity and castration resistance by targeting fibroblast growth factor receptor 1, which is a candidate therapeutic target for prostate cancer.10 These inspiring findings provide a new perspective for understanding the function of GREM1 in cancer progression and treatment. Nevertheless, little evidence has illustrated the certain underlying role of GREM1 in endometrial cancer. Accordingly, more efforts should be urgently dedicated to elucidating the mechanism of GREM1 in endometrial cancer.
Herein, we conducted a comprehensive bioinformatics workflow, accompanied by machine learning-based integration, to explore the underlying mechanism of GREM1 in endometrial cancer. First, the gene expression profiles of endometrial cancer were obtained through The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases, which were used as the training cohort and test cohort respectively. The mRNA and protein expression patterns of GREM1 were then explored at the pan-cancer level, including endometrial cancer. Second, we classified the training cohort into low-GREM1 and high-GREM1 groups and identified the differentially expressed genes (DEGs) between the two groups. The Weighted Gene Co-expression Network Analysis (WGCNA) and Mfuzz clustering were conducted to discern gene expression patterns highly related to the high-GRME1 group. Third, we took the intersection of genes associated with high-GREM1 generated by the DEGs analysis, WCGNA, and Mfuzz. Subsequently, these genes were subjected to a machine learning-based integration procedure for further selection.11 Fourth, we assessed the classification performance of the selected key gene on low- and high-GREM1 groups. We then investigated the potential mechanism of GREM1 through the combination of pathway enrichment analysis and Bayesian network (BN) inference on key genes. To deepen the illustration of the underlying mechanism of GREM1 in immune regulation, we further explored the association of the key genes (GREM1 was also included) with immune-cell infiltration in endometrial cancer. Overall, our study revealed the potential pathways of GREM1 in endometrial cancer from integrative bioinformatics methods, which might expand our understanding of its function in tumor development.
Material and methods
Data gathering and processing
The gene expression profiles of endometrial cancer were collected from the TCGA database (https://portal.gdc.cancer.gov/ ) and the GEO database (http://www.ncbi.nlm.nih.gov/geo ). The training cohort [Cancer Genome Atlas-Uterine Corpus Endometrial Carcinoma (TCGA-UCEC)] was acquired from the TCGA database using TCGA biolinks R package, including 35 normal endometrial tissues and 554 endometrial cancer tissues.12 The test cohort was obtained from five GEO datasets using the GEOquery R package, including GSE2109, GSE17025, GSE36389, GSE106191, and GSE115810.13–16 The probes were then transformed to official gene symbols through the platform annotations. The raw count data from TCGA was converted into transcripts per kilobase million format and subjected to the log2 transformation process. The raw count data from GEO were processed with the limma and sva R packages for log2 transformation, normalization, and batch effect removal. The five GEO datasets were merged into the test cohort, which contained 216 endometrial cancer samples. The sample size of the included analysis of GREM1 expression in pan-caner and endometrial cancer datasets is shown in Table 1.
Table 1Sample size of datasets used in this study
| Dataset | Sample size |
|---|
| TCGA-UCEC | 589 samples (35 controls; 554 tumors) |
| GSE2109 | 17 tumors |
| GSE17025 | 91 tumors |
| GSE36389 | 20 tumors |
| GSE106191 | 64 tumors |
| GSE115810 | 24 tumors |
Leveraging the TCGA, genotype-tissue expression (GTEx), and clinical proteomic tumor analysis consortium (CPTAC) databases, we explored the mRNA and protein expression distribution of GREM1 across extensive cancer types. The datasets used in this part can be found in Table 2. The tumor immune estimation resource (TIMER) online tool (https://cistrome.shinyapps.io/timer/ ) was to investigate the GREM1 expression pattern in the TCGA database.17 The R package ggplot2 was used to characterize the GREM1 expression pattern of TCGA target GTEx datasets from The University of California Santa Cruz (UCSC) Genome Browser database (https://xenabrowser.net/ ).18 Subsequently, the protein expression distributions of GREM1 on the CPTAC database were evaluated by The University of ALabama at Birmingham CANcer (UALCAN) online tool (http://ualcan.path.uab.edu ).19 Additionally, the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/ ) was used to demonstrate the differential protein expression of GREM1 in normal endometrial tissue versus endometrial cancer tissue.
Table 2Datasets used for pan-cancer analysis
| Tool/database for pan-cancer analysis | Dataset list |
|---|
| TIMER | TCGA-ACC, TCGA-BLCA, TCGA-BRCA, TCGA-CESC, TCGA-CHOL, TCGA-COAD, TCGA-DLBC, TCGA-ESCA, TCGA-GBM, TCGA-HNSC, TCGA-KICH, TCGA-KIRC, TCGA-KIRP, TCGA-LAML, TCGA-LGG, TCGA-LIHC, TCGA-LUAD, TCGA-LUSC, TCGA-MESO, TCGA-OV, TCGA-PAAD, TCGA-PCPG, TCGA-PRAD, TCGA-READ, TCGA-SARC, TCGA-SKCM, TCGA-STAD, TCGA-TGCT, TCGA-THCA, TCGA-THYM, TCGA-UCEC, TCGA-UCS, TCGA-UVM. |
| UCSC | GBM, GBMLGG, LGG, UCEC, BRCA, CESC, LUAD, ESCA, STES, KIRP, KIPAN, COAD, PRAD, STAD, HNSC, KIRC, LUSC, LIHC, WT, SKCM, BLCA, THCA, READ, OV, PAAD, TGCT, UCS, LAML, PCPG, ACC, KICH, CHOL. |
| UALCAN | Breast cancer, colon cancer, uterine corpus endometrial carcinoma, lung adenocarcinoma, pancreatic adenocarcinoma, head and neck squamous carcinoma, glioblastoma. |
Analysis of GREM1 expression on prognosis and clinicopathological characteristics in different cancers
Herein, we focused on the prognosis and clinicopathological feature analyses of 12 cancer types: breast invasive carcinoma (BRCA), cholangiocarcinoma, head and neck squamous cell carcinoma (HNSC), lung adenocarcinoma, lung squamous cell carcinoma (LUSC), stomach adenocarcinoma (STAD), glioblastoma multiforme(GBM), kidney chromophobe (KICH), kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), thyroid carcinoma (THCA), and uterine corpus endometrial carcinoma (UCEC), which displayed the differential expression of GREM1 compared with the normal controls, in the TCGA database. To gain insight into the effect of GREM1 expression on the prognosis and clinicopathological characteristics of tumor patients, the gene expression profiling interactive analysis tool was first used to perform the prognosis analysis. Kaplan-Meier (K-M) curves for overall survival (OS) were performed to investigate the prognostic differences between the low-GREM1 and high-GREM1 expression groups. Clinicopathological characteristics of 12 cancer types, containing age, tumor grade, T/N/M stage, and tumor stage, were downloaded from the UCSC database. Samples with missing or incomplete clinicopathological information were excluded from the analysis. Then, we delineated and explored the GREM1 expression distributions in different clinicopathological characteristic groups using the ComplexHeatmap and ggpbur R packages.
Identification of DEGs in the low-GREM1 and high-GREM1 expression groups
The normal endometrial samples of the training cohort were excluded. The endometrial cancer samples were classified into low-GREM1 group and high-GREM1 group according to the medium value of GREM1 expression. Leveraging the limma R package, the differential expression analysis was performed to screen out the DEGs between the low-GREM1 and high-GREM1 groups. The screening criteria for DEGs were set as |log2Fold Change| > 1 and adjusted p < 0.05.20 Then, the pheatmap R package was used to showcase the differential expression pattern of DEGs among the low-GREM1 and high-GREM1 groups. Moreover, we conducted Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis on these screened DEGs, and the enriched pathways with q < 0.05 were considered significant.20 The CBNplot R package was used to construct a regulatory network of the 10 most enriched KEGG pathways with the Bayesian network and gene expression profile.21
Weighted gene co-expression network analysis
The WGCNA R package was used to generate a gene co-expression network using the expression profile of the training cohort.22 Sample labels in the cohort were previously grouped into low-GREM1 and high-GREM traits. The optimal soft threshold power (β = 1–20) was initially calculated to determine the scale-free topology of the network. Then, we constructed a weighted adjacency matrix and transformed it into a topological overlap matrix (TOM). Moreover, the dissTOM was obtained for hierarchical clustering. The dynamic tree-cutting method was adopted to identify various modules clustered by gene similarity. Herein, we set minModuleSize to 60 and MEDissThres to 0.3. Subsequently, we related the recognized modules to two traits (low-GREM1 and high-GREM1). Genes in the module that had high relevance with high-GREM1 were retained for subsequent analysis.
Mfuzz clustering analysis
The Mfuzz R package uses the fuzzy c-means clustering algorithm to cluster GREM1 expression patterns.23 We first obtained distinct clustering patterns according to the ascending GREM1 expression. To investigate the relationship of clustering patterns with GREM1, each sample of the expression profile was subjected to single sample Gene Set Enrichment Analysis (ssGSEA) to assign the scores of clustering patterns.24 We further explored the clustering characteristics between low-GREM1 and high-GREM1 expression. The correlation of clustering patterns and GREM1 was investigated, with the patterns closely related to GREM1 selected for further analysis.
Potential biological functions, pathways, and diseases of the hub genes
We intersected the genes in the WGCNA module and Mfuzz clustering patterns highly related to GREM1 with DEGs, thus gaining the hub genes linked to GREM1 expression. To explore the potential biological functions, pathways, and associated diseases of hub genes, gene ontology (GO), KEGG, and disease ontology (DO) enrichment analyses were conducted using the clusterProfiler and DOSE R packages.
Key genes generated from the machine learning-based combinations and protein-protein interaction network
To perform feature selection on hub genes, we integrated 12 machine-learning algorithms and 113 algorithm combinations. The machine learning algorithms for integration included absolute shrinkage and selection operator, ridge, elastic network, stepwise multiple generalized linear model, support vector machine, a generalized linear model with likelihood-based boosting (GLMBoost), linear discriminant analysis, partial least squares regression for generalized linear models, RandomForest, gradient boosting, eXtreme gradient boosting, and NaiveBayes. One algorithm was used to select features and another was used to construct a classification model using the leave-one-out cross-validation framework.11 For each model integration, area under curve (AUC) scores within the training and test cohort were calculated under cross-validation. The classification model performed with the highest AUC score was considered optimal, thus the selected features for the model establishment were then obtained. Furthermore, the String database was used to construct the feature gene-dominated Protein-Protein Interaction (PPI) network, and the top 10 key genes with high betweenness with GREM1 in the network were obtained by the Cytoscape software (version 3.9.1).25,26 We then visualized the chromosomal locations of the 10 key genes. Moreover, we assessed the principal component analysis (PCA) scores, expression distributions in low-GREM1 and high-GREM1 groups, and correlations of key genes within the training and test cohorts.
Validation of key genes for distinguishing from the low-GREM1 and high-GREM1 groups
To assess whether the screened key genes could accurately distinguish from the low-GREM1 and high-GREM1 groups, the glmnet R package was first used to construct a logistic regression model using key genes within the training cohort. The model was illustrated in a nomogram developed with the rms R package. The ability of the model to classify low-GREM1 and high-GREM1 groups was then systemically assessed by AUC scores, calibration curves, and decision curve analysis (DCA) in the training and test cohorts.27
Potential pathways and regulatory networks of key genes with GREM1
The regulatory relationship between GREM1 and key genes was inferred through the CBNplot R package. To further investigate the key gene-regulated pathways in endometrial cancer, the clusterProfiler R package was used to conduct KEGG enrichment analysis of the key genes. The enriched KEGG pathways with q values <0.05 were considered significant differences.20 Additionally, we depicted the regulatory network of pathways, and the genes enriched in the core pathways with more edges and interaction strengths using the CBNplot R package. To verify the in-silico results, the experimental evidence at the protein level was retrieved from the HPA database to show the expression of the key genes in endometrial cancer tissues.
Correlation of key genes with GREM1 with immune infiltration cells
Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) and ssGSEA were applied to estimate infiltration degrees of different immune cell types, in which CIBERSORT performs deconvolution and ssGSEA marks enriched scores on bulk RNA-seq data using immune-cell signature.28,29 The permutations of CIBERSORT were set as 1,000 and the calculated results with p < 0.05 were retained for subsequent analysis.30 We explored the infiltration patterns of immune cells in low-GREM1 and high-GREM1 groups. The relationship between key genes and immune-cell infiltration was then examined by Spearman’s correlation approach.
Statistical analysis
Statistical analysis was performed with R programming language (version 4.2.0). Wilcoxon test was performed to evaluate the mRNA expression of GREM1 in tumor tissue compared with normal tissues. Wilcoxon test was also carried out to investigate the expression distribution of GREM1 in different groups based on clinicopathological characteristics, including age, tumor grade, T/N/M stage, and tumor stage. An unpaired t-test was used to estimate the differential expression of GREM1 protein associated with key genes, and immune-cell infiltration levels in the two groups. Log-rank test was used to analyze differences in OS between the low-GREM1 and high-GREM1 expression groups. Spearman’s correlation method was adopted to analyze the association of GREM1 and Mfuzz clustering patterns and Pearson’s method for hub gene co-expression relationships. Spearman’s correlation analysis was performed to show whether there was a close relationship between extracellular matrix (ECM) signatures and the expression of GREM1 in different cancers, and p < 0.05 was considered statistically significant.
Results
GREM1 was downregulated in different types of cancer and endometrial cancer
First, the expression profiles of GREM1 across different cancer types were investigated via TIMER. Figure 1a shows that lower expression levels of GREM1 were observed in glioblastoma multiforme (GBM), KICH, KIRC, KIRP, THCA, and UCEC (all < 0.05). Next, we found that GREM1 was significantly under-expressed in 14 types of cancer compared with normal tissue in the TCGA target GTEx datasets, including GBM, glioma (GBMLGG), low-grade glioma (LGG), UCEC, cervical squamous cell carcinoma and endocervical adenocarcinoma, KIRP, pan-kidney, KIRC, high-risk Wilms tumor, skin cutaneous melanoma, THCA, uterine carcinosarcoma, adrenocortical carcinoma, and KICH (Fig. 1b). UALCAN was used to explore the protein expression levels of GREM1 in seven cancer types, which showed significant downregulation of GREM1 in UCEC and glioblastoma (all p < 0.05; Fig. 1c). We also observed significant downregulation of GREM1 in endometrial cancer compared with control tissue (all p < 0.05; Fig. 1d and f). Immunohistochemical staining patterns of GREM1 in normal endometrial tissue and endometrial cancer tissue were retrieved from the HPA database. The staining scores shown in Figure 1g and h show that GREM1 expression was significantly lower in endometrial cancer than in normal endometrial tissue. These results show that GREM1 was downregulated in endometrial cancer.
High expression of GREM is related to poor prognosis and clinicopathological characteristics
Leveraging the 12 tumors with differential expression patterns compared with normal controls in the TCGA database, we analyzed the effect of GREM1 expression on tumor prognosis and clinicopathological characteristics. To evaluate whether the expression of GREM1 is related to the prognosis of tumor patients, K-M survival curves and log-rank tests compared differences in OS in groups with low-GREM1 and high-GREM1 expression (Fig. 2). As shown in Figure 2e (LUSC), Figure 2i (KIRC), Figure 2j (KIRP), and Figure 2l (UCEC), OS was shorter in groups with high-GREM1 expression than it was in groups with low-GREM1 expression (all p < 0.05). These results indicate that the prognostic outcome of tumor patients with high GREM1 expression was poor. Additionally, we investigated the relationship between GREM1 expression and clinicopathological characteristics in 11 tumor types (excluding GBM datasets with insufficient clinicopathological features). Figure 3, shows the expression distribution of GREM1 in groups with diverse clinicopathological characteristics. Interestingly, we observed that the expression of GREM1 was closely linked to at least one characteristic in different tumor types, including BRCA (Fig. 3a), HNSC (Fig. 3c), LUSC (Fig. 3e), STAD (Fig. 3f), KIRC (Fig. 3h), KIRP (Fig. 3i), THCA (Fig. 3j), and UCEC (Fig. 3k). Compared with patients with lower N-stage tumors, those with advanced N-stage had significantly different GREM1 expression levels in BRCA (Fig. 3a), HNSC (Fig. 3c), LUSC (Fig. 3e), KIRC (Fig. 3h), and THCA (Fig. 3j). Additionally, the expression of GREM1 was upregulated in advanced T-stage patients with STAD (Fig. 3f), KIRC (Fig. 3h), KIRP (Fig. 3i), and THCA (Fig. 3j). Regarding tumor grading and staging, the over-expression patterns of GREM1 in advanced grade and stages were observed in STAD (Fig. 3f), KIRC (Fig. 3h), KIRP (Fig. 3i), THCA (Fig. 3j), and UCEC (Fig. 3k). It is noteworthy that endometrial cancer young patients with advanced tumor stages had higher GREM1 expression. Taken together, these results show that the high expression of GREM1 is relevant to poor prognosis and worse clinicopathological features of LUSC, KIRC, KIRP, and UCEC, which indicate its pivotal function in tumor development.
Identification of DEGs in low and high GREM1 groups
To identify the high GREM1-related genes, we first divided 554 endometrial cancer samples into low-GREM1 group (n = 277) and high-GREM1 group (n = 277) according to the medium expression value of GREM1 and performed DEGs analysis. We screened out a total of 65 DEGs, of which 64 genes were upregulated and one gene was downregulated in the high-GREM1 group (Fig. 4a). The heatmap shows a distinct expression pattern of 65 DEGs between the low-GREM1 and high-GREM1 groups (Fig. 4b). KEGG enrichment analysis showed that upregulated DEGs participated in the organization of the extracellular matrix, integrin cell-surface interactions, ECM proteoglycans, collagen chain trimerization, and degradation of the extracellular matrix (Fig. 4c). We further depicted the interaction network of the top 10 KEGG enriched pathways through the Bayesian inference (Fig. 4d), indicating the core roles of extracellular matrix organization, integrin cell-surface interactions, ECM proteoglycans, and degradation of the extracellular matrix with more edges and gene counts. These results showed that upregulation DEGs may highly cooperate with the extracellular matrix in endometrial cancer.
Crucial module with GREM1 by WGCNA
To further identify the high GREM1-related genes, WGCNA was carried out to develop a gene co-expression network and identify genes from the high GREM1-correlated module. First, the optimal soft threshold power (β = 2) was determined to generate a gene co-expression network (Fig. 5a). Then, seven modules were identified after the dynamic tree cut and merge (Fig. 5b). We then related the traits (low-GREM1 and high-GREM1 groups) to the modules, obtaining a brown module with the highest positive association to the high-GREM1 group (R = 0.4, p < 0.001; Fig. 5c). A significant positive correlation of gene significance (GS) and module membership (MM) was observed in the ME brown module (R = 0.51, p < 0.001; Fig. 5d), which suggested the genes from the ME brown module were associated with GREM1 and the most significant elements of the modules correlated with GREM1. Finally, we obtained 1,043 genes from the ME brown module for further analysis.
Critical patterns with GREM1 by Mfuzz clustering analysis
Mfuzz clustering analysis was first conducted according to the ascending GREM1 expression levels, combined with ssGSEA marked clustering scores for expression characteristics, to screen out GREM1 a highly positive expression pattern. A total of 50 expression patterns clustering on GREM1 expression were obtained (Fig. 6a and b). We then carried out ssGSEA to classify tumor samples into distinct patterns with the identified Mfuzz clustering results, and each sample was assigned to a certain clustering pattern. Subsequently, Spearman’s correlation between distinct clustering patterns and GREM1 was investigated (Fig. 6c). Three critical patterns significantly positive for GREM1 were then identified (R > 0.4, p < 0.001; Fig. 6d, e, and f), in which 901 genes obtained from these patterns were obtained for further selection.
Inference of biological functions, pathways, and diseases of the hub genes
We took the intersection of 64 upregulated DEGs, 1,043 genes in the ME brown module, and 901 genes in three critical clustering patterns as 57 hub genes with GREM1 (Fig. 7a). GO, DO, and KEGG enrichment analyses were then performed to explore the potential biological functions, pathways, and diseases of the selected hub genes. As shown in Figure 7b, the majority of the hub gene-enriched GO terms were related to extracellular matrix and collagen, mainly including extracellular matrix organization, collagen-containing extracellular matrix, and extracellular matrix structural constituents. The top 20 significant terms of DO enrichment analysis are shown in Figure 7c. Notably, the hub genes were enriched in several DO terms involving gynecologic tumors, such as uterine benign neoplasm, uterine fibroid, female reproductive organ benign neoplasm, and reproductive organ benign neoplasm, which indicated the linkage of hub genes and uterine tumors. Moreover, the KEGG enrichment analysis suggested that the upregulation of hub genes may participate in extracellular matrix-related pathways, including ECM-receptor interaction, focal adhesion, and proteoglycans in cancer (Fig. 7d). Figure 7e further depicts the interaction of enriched genes and the top five significant pathways, in which FN1 and THBS1 were the common genes enriched in ECM-receptor interaction, focal adhesion, and proteoglycans in cancer.
Machine learning-based combination and PPI network identified key genes
A total of 113 machine learning-based integration was fitted to the training and test cohorts via the leave-one-out cross-validation framework. As shown in Figure 8a, the combination of GLMBoost and RF was the optimal model performed with the highest average AUC score (0.698) in the training and test cohorts. Subsequently, we obtained 12 genes (FAP, THBS1, POSTN, INHBA, ASPN, COL3A1, DES, IGFBP5, COL8A1, FN1, COL5A2, EMILIN1) screened by the model. We then selected the 10 genes with the highest betweenness in the PPI network of GREM1 (Fig. 8b), thus identifying a final set of 10 key genes highly associated with GREM1 (FAP, THBS1, POSTN, INHBA, ASPN, COL3A1, IGFBP5, COL8A1, FN1, COL5A2). The chromosomal locations of GREM1 and 10 genes are shown in Figure 8c. The PCA score plots showing the 10 genes clearly distinguish the low-GREM1 and high-GREM1 groups within the training and test cohorts (Fig. 8d and e). The significant expression differences of THBS1, ASPN, COL3A1, COL5A2, COL8A1, INHBA, POSTN, and FAP between low-GREM1 and high-GREM1 groups in the training and test cohorts are shown in Figure 8f and g. Additionally, the positive correlation patterns of GREM1 and the 10 genes within the training and test cohorts are shown in Figures 8h and i.
Assessment of key genes for distinguishing the low-GREM1 and high-GREM1 groups
We constructed logistic regression models within the training and test cohorts to assess the classification performance of the 10 genes in low- and high-GREM1 groups. The nomograms intuitively visualized the logistic regression model established in the 10 genes within the training cohort (Fig. 9a). The excellent differentiation on low- and high-GREM1 groups of the models was evaluated by the AUC scores (all AUC > 0.8; Fig. 9b and c). In addition, the calibration curves showed a minor error between the predicted probability and actual probability, showing the high classification accuracies of the models (Fig. 9d and e). The DCA curves showed that the model deviated from the none and all baselines below the high-risk threshold, which ranged from 0 to 1 and further demonstrated the classification ability and certain application of the models in the training and test cohorts (Fig. 9f and g). These results suggest the critical role of 10 genes closely associated with GREM1 in endometrial cancer.
Inference of the pathways regulated by key genes and GREM1
The Bayesian inference-generated regulation network between GREM1 and key genes is shown in Figure 10a. The network demonstrated the intricate crosstalk between key genes that dominated through GREM1 as the upstream regulator. To reveal the potential pathways regulated by 10 genes and GREM1, we carried out KEGG enrichment analysis. Figure 10b shows the top 12 significant pathways of the 10 genes and GREM1. Among the enriched pathways, many extracellular matrix and collagen-related pathways were observed, such as integrin cell-surface interactions, extracellular matrix organization, nonintegrin membrane-ECM interactions, ECM proteoglycans, degradation of the extracellular matrix, collagen chain trimerization, assembly of collagen fibrils and other multimeric structures, and collagen degradation. We performed Bayesian inference of the enrichment pathways and gene expression profile to elucidate the pathway interactions. We obtained an interaction network of the enriched 12 pathways (Fig. 10c) in which integrin cell-surface interactions (Fig. 10d), extracellular matrix organization (Fig. 10e), ECM proteoglycans (Fig. 10f), degradation of the extracellular matrix (Fig. 10g), and collagen degradation (Fig. 10h) were cores with more edges and similarities in the network. The enriched gene interactions in the core pathway are also shown. Generally, five enriched genes (FN1, COL3A1, THBS1, COL5A2, and COL8A1) were observed in the networks. Of these genes, FN1 and THBS1 proved to be significant in the previous KEGG enrichment analysis, suggesting they have potential crucial roles along with GREM1 in endometrial cancer. The HPA database was searched to obtain experimental evidence of the expressions of these five common enriched key genes at the protein level. Immunohistochemical analysis of four genes, COL3A1, COL8A1, FN1, and THBS1, were obtained (Fig. 10i–l). According to the staining and intensity levels (Fig. 10i, k and l), upregulation of COL3A1, FN1, and THBS1 was observed in endometrial cancer tissue and COL8A1 protein was not detected (Fig. 10j). The experimental evidence verified the in-silico results and further indicate transactivation of these key genes in endometrial cancer, particularly over-expression of COL3A1, FN1, and THBS1 at the protein level. To gain insight into the GREM1-regulated mechanism, we also investigated the relationship of fibrotic status and GREM1 expression in 12 types of cancer. Leveraging the ECM signatures reflecting the fibrotic status, we observed that it was positively correlated with GREM1 expression of (all R > 0, p < 0.001) in these cancer types (Fig. 11). These results further showcased the role of GREM1 in tumor fibrosis of various cancers.
Expression profiles of immune infiltration cells and correlation analysis
To investigate whether the immune infiltration cells varied in the low-GREM1 and high-GREM1 groups, we performed CIBERSORT and ssGSEA of the gene expression profiles. The results showed the expression profiles of immune infiltration cells (Fig. 12a and b), in which CD8+ T cells, CD4+ memory resting T cells, and M0 macrophages were the main cell types that infiltrated endometrial cancer (Fig. 12a). The infiltration of 10 immune cells in the low-GREM1 and high-GREM1 groups was evaluated by CIBERSORT (Fig. 12c). Regarding the ssGSEA results, a total of 16 immune cells had significantly different expression in the low- and high-GREM1 groups (Fig. 12d). Additionally, correlation of immune-cell infiltration and 11 genes (including GREM1) was explored (Fig. 12e and f). Interestingly, infiltration of mast cells was upregulated in the high-GREM1 groups and was positively correlation with ASPN, COL3A1, FAP, GREM1, IGFBP5, POSTN, THBS1, which suggested the mast cells were highly influenced by GREM1 in endometrial cancer.
Discussion
GREM1 involvement in carcinogenesis has been investigated and is seen as a prospective therapeutic target in various malignancies.7–10,31 Nevertheless, the role of GREM1 in endometrial cancer and its association with immune cells is not well understood. Herein, a comprehensive bioinformatics analysis combined with a machine learning-based combination was applied to reveal the potential mechanism of GREM1 in endometrial cancer. Inspired by evidence that GREM1 over-expression inhibits epithelial-to-mesenchymal transition in pancreatic cancer,9 we wondered whether a similar mechanism occurred in endometrial cancer. We focused the analysis on the expression profiles of cancer tissues with high-GREM1 expression in the TCGA and GEO databases. We first found downregulation of GREM1 mRNA and protein in endometrial cancer, which was consistent with the expression pattern previously reported in pancreatic cancer.9 Downregulation of GREM1 expression was observed in glioblastoma. Although we found overall downregulation of GREM1 in endometrial cancer, the heterogeneity of cancer exists, suggesting a certain cell population with a higher GREM1 expression pattern. Additionally, we performed prognosis and clinicopathological analysis of GREM1 in 12 cancer types in the TCGA database. The expression of GREM1 was significantly different from that in normal controls (up/downregulation pattern). The results indicated that high expression of GREM1 predicted poor prognosis (shorter OS) and clinicopathological features (advanced tumor grade/stage) in lung squamous cell carcinoma, kidney renal clear cell carcinoma, kidney renal papillary cell carcinoma, and endometrial carcinoma. We thus suggest that the patients with tumors that have high GREM1 expression have poor clinical outcomes, indicating that GREM1 contributes to tumor progression.
Then, we classified the endometrial cancer samples into low-GREM1 and high-GREM1 groups. We screened out a total of 65 DEGs (64 upregulated genes and one downregulated gene) between these two groups. Subsequently, WGCNA was performed to discern the ME brown module highly positively linked to high-GREM1 trait, of which 1,043 genes were included. Together with ssGSEA scoring and correlation analysis, Mfuzz clustering was also conducted to screen out three high-GREM1-related gene expression patterns containing 901 genes. We took the intersection of the 64 upregulated DEGs, 1,043 genes identified by WGCNA, and Mfuzz-selected 901 genes. We observed a set of 57 overlapping genes, indicating that the prefeature selection procedure using multiple bioinformatics methods was reliable and accessible. In addition, we conducted DO, GO, and KEGG enrichment analysis on the selected 57 hub genes to investigate the biological processes and pathways in endometrial cancer. Interestingly, the hub genes were enriched in uterine benign neoplasms, uterine fibroids, female reproductive organ benign neoplasms, and reproductive organ benign neoplasms, suggesting an association between the hub genes and uterine carcinomas. Furthermore, GO and KEGG analysis indicated that the hub genes were enriched in extracellular matrix-related terms, including extracellular matrix organization, collagen-containing extracellular matrix, ECM-receptor interaction, focal adhesion, and proteoglycans in cancer. We speculate that the hub GREM1 genes participate in ECM-related pathways. Deregulated ECM dynamics, such as collagen deposition or an increase in ECM stiffness,32 are signatures of tumor proliferation and metastasis. It is well known that ECM is a complex and dynamic environment mainly composed of collagen, proteoglycans, laminin, and fibronectin.33 Thus, we also wondered which specific ECM composition was a great contribution or highly correlated to the higher level of GREM1, further exploring the underlying regulatory relationship with GREM1 in endometrial cancer.
To define the underlying molecular mechanism of GREM1 more precisely in endometrial cancer, feature selection on the hub genes was still needed. Herein, the machine learning-based integrative procedure was taken for feature selection. We obtained a total of 12 genes (FAP, THBS1, POSTN, INHBA, ASPN, COL3A1, DES, IGFBP5, COL8A1, FN1, COL5A2, and EMILIN1) via best-performance algorithm integration (GLMBoost target RF) with the highest AUC scores. To explore whether the functional associations between the 12 genes and GREM1 exist, we constructed a PPI network interacting with 12 genes and GREM1. Using Cytoscape software, we calculated the betweenness of 12 genes and GREM1 in the network. DES and EMILIN1, with zero betweenness were excluded, and a PPI network composed of 10 genes and GREM1 was then visualized. We thus obtained a final set of 10 genes (FAP, THBS1, POSTN, INHBA, ASPN, COL3A1, IGFBP5, COL8A1, FN1, and COL5A2), which are considered as key GREM1 genes. Using three-dimensional PCA score plots, we observed the classification performances of key genes in low-GREM1 and high-GREM1 groups. The elevated expression of key genes in the high-GREM1 groups was then demonstrated, following their positive expression correlations with GREM1, which implies that the co-expressions of key genes and GREM1 exist in endometrial cancer.
To systemically assess the classification abilities of key genes from the low-GREM1 and high-GREM1 groups, we established logistic regression models, following ROC curves, calibration curves, and DCA curves as evaluation indexes. We noted robust AUC scores, strong calibrations, and distinct utility of the discrimination models within the training and test cohorts. These results show that the key genes were clearly distinguished in the low-GREM1 and high-GREM1 groups, indicating their expression was similar to that of GREM1. The key genes were subjected to KEGG enrichment analysis, together with conducting BN inference on the enriched pathways and gene expression patterns. Interestingly, the key genes were also enriched in extracellular matrix-related pathways, especially those related to ECM structural constituent and dynamic changes. Through the pathway interaction based on BN inference, we found integrin cell-surface interactions, extracellular matrix organization, ECM proteoglycans, degradation of the extracellular matrix, and collagen degradation were core pathways. We then investigated the gene interactions of the core pathways and observed five common genes (FN1, COL3A1, THBS1, COL5A2, and COL8A1) enriched in the interactions. Notably, the collagen genes COL3A1, COL5A2, and COL8A1 almost occurred in the enriched pathways, suggesting the critical role of collagen in ECM-related pathways. Moreover, we obtained immunohistochemical staining results of FN1, COL3A1, THBS1, and COL8A1 in endometrial cancer from the HPA database, which was as external experimental evidence supporting our analysis. These results showcased the elevated expression of FN1, COL3A1, and THBS1 proteins in endometrial cancer, which strengthened the in-silico conclusions.
Recent studies have demonstrated that COL3A1 over-expression leads to poor prognosis or promotes cancer cell progression in several malignancies, including bladder cancer, colorectal cancer, and triple-negative breast cancer.34–37 Similarly, recent evidence has revealed that COL5A2 or COL8A1 expression was related to poor prognosis or tumor progression,38–45 also highlighting their potential diagnostic values for different malignancies.38,46,47 FN1, acting as an upstream or downstream regulator, was previously shown to promote or inhibit tumor progression in various malignancies.48–51 THBS1 was found to positively associate with aggressive tumor development and progression.52–54 Additionally, elevated microvessel counts were detected in endometrial cancer patients with high THBS1 expression, suggesting a role of THBS1 in tumor angiogenesis.55 But few studies reported roles of FN1, COL3A1, THBS1, COL5A2, and COL8A1 in endometrial cancer, more studies are warranted. Taken together, we hypothesize that GREM1 is involved in the modulation of ECM-related pathways, especially those mainly regulated by collagen genes. Modulation of the ECM may allow for tumor invasion or metastasis, which has been shown in studies of estrogen receptor-negative breast cancer.56
The tumor microenvironment, especially immune-cell infiltration, was shown to be linked to ECM modulation in tumor progression.57,58 Leveraging the CIBERSOT and ssGSEA algorithms, the infiltration of immune cells in low-GREM1 and high-GREM1 groups was measured to explore the linkage of GREM1 and the tumor microenvironment of endometrial cancer. We aimed to identify the specific immune cells correlated with GREM1 and its interactive key genes, with the upregulation of infiltration in high-GREM1 groups. Differential expression analysis showed that the infiltration of mast cells was elevated in high-GREM1 groups. Analysis of correlations with immune-cell subtypes revealed that mast cell infiltration levels were positively associated with GREM1 and key genes. Previous evidence supports the linkage between increased mast cell count and the progression of endometrial cancer, such as angiogenesis.59 Furthermore, the high density of mast cells is significantly associated with myometrial invasion in endometrial cancer, which shows the function of mast cells in tumor progression.60 However, because of the limited number of reports on how mast cells contribute to the progression of endometrial cancer, more functional studies are required to clarify their role. Overall, our data indicate that mast cell infiltration in endometrial cancer may be influenced by ECM modulation, which is correlated with GREM1 and key genes.