Some patients with endometrial cancer (EC) suffer from limited survival benefits after immunotherapy, suggesting that there may be a specific pattern associated with immunotherapy. Immune-related genes were extracted from The Cancer Genome Atlas databases. We analyzed the differences among immune subtypes (ISs) in the distribution of the tumor mutational burden, chemotherapy-induced immune response markers, immune checkpoint-related genes, immunotherapy, and chemotherapy. We applied dimensionality reduction and defined the immune landscape of EC. Then, we used the Weighted Gene Co-Expression Network Analysis package to identify the coexpression modules of these immune genes. Finally, hub genes were selected and detected by quantitative PCR and immunohistochemistry. We obtained three ISs. There were differences in the distribution of the tumor mutational burden, chemotherapy-induced immune response markers, and immune checkpoint–related genes among the ISs. Regarding immunotherapy and chemotherapy, the IS2 subtypes were more sensitive to programmed cell death protein 1 inhibitors. In addition, different positions in the immune landscape map exhibited different prognostic characteristics, providing further evidence of the ISs. The IS2 subtypes were significantly positively correlated with yellow module gene list, indicating a good prognosis with high score. SIRPG and SLAMF1 were identified as the final characteristic genes. The quantitative PCR and immunohistochemistry results showed that the expression levels of SIRPG and SLAMF1 were low in human EC tissue. In this study, we identified three reproducible ISs of EC. The immune landscape analysis further revealed the intraclass heterogeneity of the ISs. SIRPG and SLAMF1 were identified to be associated with progression, suggesting that they may be novel immune-related biomarkers of EC.

Endometrial cancer (EC) is the fourth most commonly diagnosed cancer among women (1), and its incidence has continued to increase over the past few years; it is projected to continue to increase in the future (2, 3). The increase in morbidity is attributed to various factors, including the obesity epidemic, the widespread reduction in the use of menopausal hormone therapy (including progesterone), changes in reproductive behavior, and the increase in the prevalence of diabetes (46). The 5-year survival rate in localized stage EC is good, reaching 95%, and the rate of patients with distant metastasis is 17% (7). Paclitaxel plus carboplatin, a standard first-line treatment, has limited efficacy in advanced, recurrent, and metastatic EC. Recently, pembrolizumab has demonstrated efficacy in patients with metastatic high microsatellite instability EC after front-line chemotherapy failure (8, 9). In the KEYNOTE-158 study (9), 49 patients with high microsatellite instability and recurrent EC treated with single-agent pembrolizumab had an overall response rate of 57%, with 16% of the women having complete responses and 41% having partial responses. These results prompted us to explore the immune microenvironment of EC to identify patients who could benefit from immunotherapy.

The tumor immune microenvironment (TIM) is composed of immune cells, mesenchymal cells, endothelial cells, inflammatory mediators, and extracellular matrix molecules (10, 11). The occurrence and development of EC are closely related to the regulation of the TIM. Through a series of mechanisms, cancer cells eventually escape immune surveillance and inhibit the cytotoxic function of tumor-antagonizing immune cells (12). In recent years, several immunotherapy regimens for therapeutic use have been developed based on the immune evasion features of cancer cells. The use of immune cells to attack cancer cells introduces a new strategy for cancer treatment. However, few studies show the immune landscape of EC, and the understanding of the TIM of EC is limited, largely hindering the selection of ideal populations suitable for immunotherapy. Furthermore, the inability to distinguish the difference in the tumor microenvironment between those patients who cannot benefit from immunotherapy and those who benefit from immunotherapy hinders the development of new targets for overcoming immune resistance.

In this study, we identified three reproducible EC immune subtypes (ISs) through a multicohort retrospective study and used independent data for subtype verification and comprehensive molecular identification. The results revealed that each IS was related to different gene expression profiles, showing different patterns in tumor genetic aberrations, infiltrating immune cell composition, functional orientation, and cytokine profile. In addition, different ISs had significant differences in clinical prognosis. This study has clinical significance for selecting populations with advantages in immunotherapy, designing new immunotherapy strategies, and overcoming immune resistance, providing a conceptual framework for understanding the TIM of EC.

The gene expression profile data of EC samples were downloaded from The Cancer Genome Atlas (TCGA) Genomic Data Commons Application Programming Interface database and further processed and transformed into transcripts per million data. The data of a total of 542 tumor samples and 23 normal samples were obtained. The clinical follow-up information data of EC were also downloaded from TCGA Genomic Data Commons Application Programming Interface, and samples lacking survival time and survival status data were eliminated.

According to previous reports (13, 14), the following categories of genes were collected as immune-related genes for the subsequent analysis: immune cell–specific genes derived from single-cell RNA sequencing (RNA-seq) data, genes of costimulatory and cosuppressive molecules, cytokines and cytokine receptor genes, genes involved in Ag processing and presentation, and other immune-related genes. In total, 2006 immune-related genes were collected.

The following steps were applied to the RNA-seq data from TCGA-University of California Santa Cruz (TCGA-UCSC) to obtain 542 ECs: (1) samples without clinical data were removed; (2) normal tissue sample data were removed; (3) genes with an expression level fragments per kilobase million (FPKM) equal to 0 in >50% of the samples were removed; and (4) the expression profile of immune cell–related genes was retained, and log transformation log2 (FPKM +1) was performed.

The following steps were performed to preprocess the RNA-seq data from TCGA-UCSC to obtain 56 samples of uterine carcinoma: (1) samples without clinical data were removed; (2) normal tissue sample data were removed; (3) genes with an expression level (FPKM) equal to 0 in >50% of the samples were removed; and (4) the expression profile of immune cell–related genes was retained, and log transformation log2 (FPKM +1) was performed.

A consensus matrix was constructed through consensus clustering (ConsensusClusterPlus), and the samples were classified by clustering (15). Using the expression data of 1967 immune-related genes obtained through screening, the ISs of the samples were obtained. We used the algorithm and the “1 − Pearson correlation coefficient” as the measurement distances, and we performed 500 bootstraps. Each bootstrap process included 80% of the patients in the training set.

The number of clusters was set from 2 to 10, and the best classification was determined by calculating the consistency matrix and the consistency cumulative distribution function (CDF). Furthermore, uniform clustering was used to group immune genes, and the same settings and parameters were used to obtain immune gene modules.

We annotated the biological functions of the immune gene modules. With the collected immune-related genes as the background, we used the DAVID (v6.8) tool to annotate the biological processes of the genes in each module in Gene Ontology (GO). We used the ANOVA algorithm to evaluate the association between the ISs and 57 immune-related molecular and cellular characteristics.

We first used stage, age, and grade as covariates in the training set samples, and we used overall survival as the end point. A log-rank test and univariate and multivariate Cox regression methods were used to evaluate the prognostic value of the ISs. In the validation set samples, an ANOVA was used to evaluate the correlation between the ISs and various immune-related molecules and cell characteristics.

Considering the dynamic characteristics of the immune system, we used graph-based learning methods to perform a dimensionality reduction analysis to reveal the internal structure of the immune system and visualize the distribution of individual patients. This analysis projects high-dimensional gene expression data into a tree structure in a low-dimensional space by retaining local geometric information (16). This method has previously been used to simulate the progression of cancer and define the development trajectory using large amounts of single-cell gene expression data (1719). In this study, we extended the analysis to immune gene expression profiles. This immune landscape reflects the relationship between patients in a nonlinear manifold, which may complement the discrete ISs defined in the linear Euclidean space.

Fresh EC tissues were provided by the Cancer Hospital of the Chinese Academy of Medical Sciences. Experienced pathologists performed the EC diagnosis and disease classification. The study was approved by the Ethics Committee of the Cancer Hospital of the Chinese Academy of Medical Sciences (17-099/1355).

The total RNA extraction was performed using RNA-easy Isolation Reagent (No. RC112-01; Vazyme). Quantitative PCR (qPCR) was performed using a HiScript III 1st Strand cDNA Synthesis Kit (No. R312-01; Vazyme) and ChamQ Universal SYBR qPCR Master Mix (No. Q712-02; Vazyme) according to the manufacturer’s instructions. The total RNA extraction, qPCR analysis, and immunohistochemistry were performed as previously described (20). The following Abs were used for immunohistochemistry: anti-SIRPG Ab (NBP2-33512; Bio-Techne, Toronto, ON, Canada) and anti-SLAMF1 Ab (ab228978; Abcam, Cambridge, UK).

The difference in overall survival between the groups was evaluated using a Kaplan–Meier analysis with a log-rank test. Student t test was used to determine the significance of the difference between two groups. All statistical analyses were performed using R software (version 3.5.3) and SPSS 22.0 software (SPSS, Chicago, IL) (*p < 0.05, **p < 0.01, ***p < 0.001, and ****p < 0.0001).

The study was approved by the Ethics Committee of the Cancer Hospital of the Chinese Academy of Medical Sciences (17-099/1355). Informed consent was obtained from all participants in the study.

The datasets used or analyzed during this study are available from the corresponding author on reasonable request.

First, the EC expression profile of immune-related genes was extracted from the RNA-seq data of TCGA-UCSC, and 1967 genes were obtained for the subsequent analysis. In total, 542 EC samples were clustered through ConsensusClusterPlus, and the optimal number of clusters was determined according to the CDF. The CDF delta area curve indicated that when the cluster was selected as 3, there was a relatively stable clustering result (Fig. 1A, 1B). We selected k = 3 to obtain three ISs (IS1, IS2, and IS3) (Fig. 1C). Further analysis of the prognostic characteristics of these three ISs revealed significant prognostic differences. In general, IS3 had a good prognosis, whereas the IS1 subtypes had a poor prognosis (Fig. 1D). In addition, we compared the relationship among these three ISs by stage and grade (no age information, and all patients were female). We observed that there were significant differences among these three subtypes in different stages. There was a significantly high proportion of the IS1 subtypes (Fig. 1E) in stage IV patients and a significantly high proportion of the IS1 subtypes (Fig. 1F) in the G3 samples, indicating that these three ISs may achieve a more accurate prognostic classification of patients. In addition, we used the same method to molecularly type TCGA patients with uterine sarcoma, and we observed that the prognosis of these three ISs also had significant differences (Fig. 1G), which was consistent with the training set. Similarly, comparing the relationship among these three ISs by stage and age (no grade information, and all patients were female), we observed that there were significant differences among the three ISs in different stages (Fig. 1H). There was a significantly high proportion of the IS1 subtypes in the stage III and stage IV samples (Fig. 1I), indicating that these three immune-based molecular subtypes are transplantable in different research cohorts.

FIGURE 1.

ISs in uterine corpus endometrial carcinoma (UCEC). (A and B) The CDF curve of TCGA-UCEC cohort samples (A) and CDF delta area curve of TCGA-UCEC cohort samples (B) indicated that when the cluster was selected as 3, there was a relatively stable clustering result. (C) We selected k = 3 to obtain three ISs (IS1, IS2 and IS3). (D) Kaplan–Meier curves of three ISs showed that IS3 had a good prognosis, whereas the IS1 subtypes had a poor prognosis. There was a significantly high proportion of the IS1 subtypes (E) in stage IV patients and a significantly high proportion of the IS1 subtypes (F) in G3 samples. (G) The prognosis of these three ISs also significantly differs in uterine sarcoma. (H) There were significant differences among the three ISs in different stages. (I) There was a significantly high proportion of the IS1 subtypes in the stage III and stage IV samples.

FIGURE 1.

ISs in uterine corpus endometrial carcinoma (UCEC). (A and B) The CDF curve of TCGA-UCEC cohort samples (A) and CDF delta area curve of TCGA-UCEC cohort samples (B) indicated that when the cluster was selected as 3, there was a relatively stable clustering result. (C) We selected k = 3 to obtain three ISs (IS1, IS2 and IS3). (D) Kaplan–Meier curves of three ISs showed that IS3 had a good prognosis, whereas the IS1 subtypes had a poor prognosis. There was a significantly high proportion of the IS1 subtypes (E) in stage IV patients and a significantly high proportion of the IS1 subtypes (F) in G3 samples. (G) The prognosis of these three ISs also significantly differs in uterine sarcoma. (H) There were significant differences among the three ISs in different stages. (I) There was a significantly high proportion of the IS1 subtypes in the stage III and stage IV samples.

Close modal

We downloaded the mutation dataset of TCGA-UCSC processed by mutect2 software, calculated the tumor mutational burden (TMB), and analyzed the distribution of the TMB in the three ISs (Supplemental Fig. 1A). The TMB in the IS1 subtype was lower than that in the IS2 and IS3 subtypes. In addition, we counted the differences in the number of gene mutations in the samples in the different ISs and observed that the number of gene mutations in the IS1 subtype was lower than that in the IS2 and IS3 subtypes (Supplemental Fig. 1B). Subsequently, we screened a total of 17,785 genes with a mutation frequency >3 in each immune molecular subtype. In each subtype, a χ2 test was used to screen genes with significantly high frequency mutations (p < 0.05), and 9671 genes were obtained. The characteristics of the top 10 genes with the highest mutation frequency among the different subtypes are shown in Supplemental Fig. 1C. The proportion of TP53 mutations in the IS1 subtypes was significantly higher than that in the IS2 and IS3 subtypes, which may be related to the poor prognosis of tumors caused by TP53 mutations. In addition, the proportion of TP53 mutations in the IS3 subtype was significantly lower than that in the IS1 subtype, which may be related to the good prognosis.

To observe the differences in the expression distribution of the classic markers of chemotherapy-induced immune responses among the three ISs, we calculated the differences in these genes in TCGA-UCSC cohort. In total, there were 26 genes in the TCGA-UCSC cohort, and 22 (84.6%) genes had significant differences in each subtype (Supplemental Fig. 1D). For instance, CALR, HMGB1, PANX1, and P2RX7 were significantly upregulated in the IS1 tumors in the TCGA cohort, while TLR4 and HGF were significantly upregulated in the IS3 tumors, suggesting that there are differences in chemotherapy-induced immune response markers among the different ISs. These differences may lead to different clinical prognoses of the disease. In addition, we obtained 47 immune checkpoint–related genes from previous studies (21) and analyzed the differences in these genes among the different ISs. In TCGA-UCSC cohort, 43 (91.5%) genes were significantly different (Supplemental Fig. 1E). For instance, VTCN1 and CD40 were significantly upregulated in the IS1 tumors in the TCGA cohort, while CD200, NRP1, and CD44 were significantly upregulated in the IS3 tumors. These results indicate that the differential expression of each IS may lead to differences in disease progression and provide potential therapeutic targets for immunotherapy.

We extracted the expression profiles of genes corresponding to the CA125, CA153, CA199, YKL-40, and DJ-1 proteins from TCGA-UCSC cohort (the results of HE4 were not significant) and analyzed their distribution differences in the different ISs. The results showed that the expression profiles of these proteins significantly differed among the three ISs. The expression levels of CA125 (Supplemental Fig. 2A), CA153 (Supplemental Fig. 2B), CA199 (Supplemental Fig. 2C), and YKL-40 (Supplemental Fig. 2D) in the IS1 IS were significantly higher than those in the IS3 subtype. However, the expression level of DJ-1 (Supplemental Fig. 2E) did not significantly differ between the two subtypes. These results suggested that immunotyping could be used as a precise prognostication strategy for OC, and it may be superior to assessing CA125, CA153, and CA199 expression in predicting patient prognosis.

Immunotherapy and chemotherapy are important strategies in the treatment of EC; thus, we analyzed the differences in the response to immunotherapy and chemotherapy among the three ISs. In this study, we used the subclass mapping method in which lower p values indicate higher similarity to compare the similarities among these three ISs and immunotherapy patients in the GEO GSE91061 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi) dataset (47 patients who received programmed cell death protein 1 [PD-1] immune check point inhibitors or cytotoxic T lymphocyte–associated protein-4 immune check point inhibitors). The results showed that the IS2 subtype was more sensitive to PD-1 inhibitors than the other two subtypes (Supplemental Fig. 2F). Furthermore, we analyzed the response of different ISs to traditional chemotherapy drugs, such as paclitaxel, cisplatin, doxorubicin, and docetaxel. We found that the IS2 subtype was more sensitive to paclitaxel and cisplatin than the other subtypes, while the IS1 subtype was more sensitive to doxorubicin and docetaxel (Supplemental Fig. 2G–J). These ISs have different potential to benefit from PD-1 inhibitors and various chemotherapeutic agents through immunophenotyping analysis, providing a theoretical basis and guidance for immunotherapy and chemotherapeutic drug selection for EC.

To compare the distribution of immune cell components among the different ISs, we obtained 28 immune cell marker genes from a previous study (22). The single sample gene set enrichment analysis method was used to score each type of immune cell to determine the score of 28 types of immune cell in each patient. We calculated 28 immune scores in the patients in TCGA-UCSC cohort, and these immune cells were mainly divided into five categories (Supplemental Fig. 3A). In addition, most immune cell components varied among the different subtypes. For example, the effector memory CD4 T cell frequency was significantly higher in the IS1 subtype than in the IS3 subtype. Th2 cells, immature dendritic cells, eosinophils, mast cells, NKT cells, and plasmacytoid dendritic cells were significantly more abundant in the IS3 subtype than in the IS1 subtype (Supplemental Fig. 3B). These findings suggest that the poor prognosis of EC may be related to the activation of effector memory CD4 T cells and the inhibition of Th2 cells, immature dendritic cells, eosinophils, mast cells, NKT cells, and plasmacytoid dendritic cells.

To observe the relationship among the three ISs identified in this study and six previously identified pan-cancer immunotypes, we extracted the molecular subtype data of these samples from a previous study (17). The IS1 subtype was mainly composed of the C1 and C2 subtypes, and the proportion of the C3 subtype in the IS3 subtype was higher than that in the IS1 subtype (Supplemental Fig. 3C). Previous pan-cancer IS studies have demonstrated that the C3 and C5 subtypes have better prognoses, whereas the C1, C2, C4, and C6 subtypes have worse prognoses, which is consistent with our results. In addition, we assessed the correlation between the immunophenotype and 56 previously defined immune molecular characteristics with a false discovery rate < 0.01. In total, 37 of the most significant immune-related features were identified as shown in Supplemental Fig. 3D. The IS1 subtype had the highest intratumor heterogeneity, proliferation, wound healing, IFN-γ response, number of segments, B cell memory, activated dendritic cells, resting mast cells, mast cells, and dendritic cells. The IS3 subtype had a significantly higher TGF-β response and activation of naive B cells and CD4 memory T cells than the IS1 and IS2 subtypes. These findings are consistent with the cellular signature results. Therefore, immunophenotyping can be used to reflect the immune status of EC and identify those who are suitable for immunotherapy.

To facilitate visualization and uncover the underlying structures of the distribution of individual patients, we applied the dimensionality reduction method based on the graph learning technique to the immune gene expression profiles. This analysis places a single patient into a graph with a sparse tree structure and defines the immune landscape of EC, with IS1, IS2, and IS3 generally distributed in different branches of the tree (Fig. 2A). The position of the patient represents the overall characteristics of the TIM of the corresponding subtype. Principal component 1 (PCA1) (horizontal axis) was negatively correlated with activated B cells, activated CD8 T cells, immature B cells, Th1 cells, myeloid-derived suppressor cells (MDSCs), and NK cells. Interestingly, PCA2 (vertical axis) was positively correlated with effector memory CD8 T cells, immature B cells, activated dendritic cells, macrophages, and MDSCs (Fig. 2B). These results suggest that there was significant intraclass heterogeneity among the subtypes. The immune landscape further revealed significant intracluster heterogeneity within each subtype. We observed that certain subtypes appeared to be more diverse and heterogeneous than others. For instance, the IS1 subtype was distributed at the opposite ends of the immune landscape, indicating that there was significant intraclass heterogeneity in this subtype, and IS1 could be further divided into two subgroups (IS1A and IS1B) based on their location in the immune landscape, showing different immune gene expression profiles in specific modules. Similar results were obtained for IS2 (IS2A, IS2B, and IS2C) and IS3 (IS3A and IS3B) (Fig. 2C). Furthermore, we analyzed the specific immune expression patterns in various subgroups. For example, IS1A showed lower scores for activated CD4 T cells, type 2 Th cells, eosinophils, and immature dendritic cells than IS2B. The scores of activated B cells, activated CD4 T cells, activated CD8 T cells, MDSCs, and neutrophils in the IS2C subtype were higher than those in IS2A and IS2B. IS3A had higher activated B cell scores, activated CD4 T cell scores, activated CD8 T cell scores, immature B cells, and MDSCs than IS3B (Fig. 2D). Therefore, tumor Ags may be more effective in IS1A, IS2A, and IS3A than in the other subgroups. In addition, the samples with extreme distribution positions in the immune landscape were subjected to a prognostic comparison (groups 1 and 5) (Fig. 2E), and the patients in group 1 showed a better survival probability (Fig. 2F). These results indicate that our immune landscape analysis provides complementary value to previously identified ISs.

FIGURE 2.

Immune landscape of EC. (A) The immune landscape of EC with IS1, IS2, and IS3 generally distributed in different branches of the tree. (B) The position of the patient represents the overall characteristics of the TIM of the corresponding subtype. PCA1 was negatively correlated with activated B cells, activated CD8 T cells, immature B cells, Th1 cells, MDSCs and NK cells. PCA2 was positively correlated with effector memory CD8 T cells, immature B cells, activated dendritic cells, macrophages, and MDSCs. (C) IS1 could be further divided into two subgroups (IS1A and IS1B) based on their location in the immune landscape, which showed different immune gene expression profiles in specific modules. Similar results were obtained for IS2 (IS2A, IS2B, and IS2C) and IS3 (IS3A and IS3B). (D) IS1A had lower scores in activated CD4 T cells, Th2 cells, eosinophils, and immature dendritic cells than IS2B. The scores of activated B cells, activated CD4 T cells, activated CD8 T cells, MDSCs, and neutrophils in the IS2C subtype were higher than those in IS2A and IS2B. IS3A had a higher activated B cell score, activated CD4 T cell score, activated CD8 T cell score, immature B cells, and MDSCs than IS3B. (E and F) Samples with extreme distribution positions in the immune landscape were subjected to a prognostic comparison (groups 1 and 5) (E), and the patients in group 1 showed a better survival probability (F).

FIGURE 2.

Immune landscape of EC. (A) The immune landscape of EC with IS1, IS2, and IS3 generally distributed in different branches of the tree. (B) The position of the patient represents the overall characteristics of the TIM of the corresponding subtype. PCA1 was negatively correlated with activated B cells, activated CD8 T cells, immature B cells, Th1 cells, MDSCs and NK cells. PCA2 was positively correlated with effector memory CD8 T cells, immature B cells, activated dendritic cells, macrophages, and MDSCs. (C) IS1 could be further divided into two subgroups (IS1A and IS1B) based on their location in the immune landscape, which showed different immune gene expression profiles in specific modules. Similar results were obtained for IS2 (IS2A, IS2B, and IS2C) and IS3 (IS3A and IS3B). (D) IS1A had lower scores in activated CD4 T cells, Th2 cells, eosinophils, and immature dendritic cells than IS2B. The scores of activated B cells, activated CD4 T cells, activated CD8 T cells, MDSCs, and neutrophils in the IS2C subtype were higher than those in IS2A and IS2B. IS3A had a higher activated B cell score, activated CD4 T cell score, activated CD8 T cell score, immature B cells, and MDSCs than IS3B. (E and F) Samples with extreme distribution positions in the immune landscape were subjected to a prognostic comparison (groups 1 and 5) (E), and the patients in group 1 showed a better survival probability (F).

Close modal

An immune gene coexpression module was used to classify hub immune-related genes that influence the effectiveness of immunotherapy. We used the Weighted Gene Co-Expression Network Analysis R software package to identify the coexpression modules of these immune genes. Based on the scale independence and mean connectivity, a soft threshold of 4 was selected to screen coexpression modules (Fig. 3A, 3B). The coexpression network conforms to the scale-free network, i.e., the logarithm log(k) of a node with a connection degree of k is negatively correlated with the logarithm log(p(k)) of the probability that the node appears, resulting in a correlation coefficient >0.85. To ensure that the network was a scale-free network, we selected β = 4. Next, we converted the expression matrix into an adjacency matrix, which was then converted into a topological matrix. Based on the topological overlap matrix, we used the average-linkage hierarchical clustering method to cluster genes according to the standard of the hybrid dynamic shearing tree and set each minimum number of genes in a gene network module to 30. After determining the gene modules using the dynamic shear method, we calculated the eigengenes of each module, performed a cluster analysis of the modules, merged the modules that were closer to each other into a new module, and set height = 0.25, DeepSplit = 4, and minModuleSize = 30. In total, 10 modules (Fig. 3C) were obtained. Notably, the gray module represented a gene list that could not be aggregated into other modules. Then, we analyzed the correlation between each module and age, histological type, stage and grade, IS1, IS2, and IS3 (Fig. 3D). The turquoise and yellow modules were significantly correlated with the IS1, IS2, and IS3 subtypes. We found that the IS2 subtype was significantly positively correlated with the genes in the turquoise (Fig. 3E) and yellow modules (Fig. 3F, Supplemental Table I).

FIGURE 3.

Identification of immune gene coexpression modules. Based on the scale independence (A) and mean connectivity (B), a soft threshold of 4 was selected to screen the coexpression modules. (C) In total, 10 modules were obtained by the coexpression network analysis. (D) We analyzed the correlation between each module and age, histological type, stage, and grade. IS1, IS2, IS3, the turquoise module, and the yellow module were significantly correlated with the IS1, IS2 and IS3 subtypes. The IS2 subtype was significantly positively correlated with the genes in the turquoise (E) and yellow (F) modules.

FIGURE 3.

Identification of immune gene coexpression modules. Based on the scale independence (A) and mean connectivity (B), a soft threshold of 4 was selected to screen the coexpression modules. (C) In total, 10 modules were obtained by the coexpression network analysis. (D) We analyzed the correlation between each module and age, histological type, stage, and grade. IS1, IS2, IS3, the turquoise module, and the yellow module were significantly correlated with the IS1, IS2 and IS3 subtypes. The IS2 subtype was significantly positively correlated with the genes in the turquoise (E) and yellow (F) modules.

Close modal

Of the nine modules, the turquoise and yellow modules were significantly associated with EC prognosis, revealing that the expression of genes in the turquoise and yellow modules was significantly associated with the prognosis of EC patients (Fig. 4A). The low score of the turquoise module indicated a poor prognosis, while the high score of the yellow module indicated a good prognosis. The GO analysis showed that the turquoise module was associated with leukocyte migration and myeloid leukocyte migration (Fig. 4B), which were highly negatively correlated with PCA1 in the immune landscape (Fig. 4C). The yellow module was related to T cell activation and regulation of T cell activation (Fig. 4D). Similarly, the yellow module was highly negatively correlated with the PCA1 component in the immune landscape (Fig. 4E). Furthermore, we extracted the gene expression profiles of the turquoise and yellow modules that were most correlated with prognosis from TCGA-UCSC dataset. Then, we used the mean value of expression as the sample feature to classify patients with high and low features, and we analyzed the prognostic differences between the high-scoring and low-scoring patients. The prognosis did not significantly differ between the groups (Fig. 4F) in the turquoise module, while the prognosis of the high-scoring group was significantly better than that of the low-scoring group in the yellow module (Fig. 4G). Finally, we selected the genes in the yellow module with a correlation >0.9 as hub genes. We identified two genes, namely, SIRPG and SLAMF1, as the final characteristic genes, and a low expression of these genes was correlated with a poor prognosis. Therefore, the hub genes can act as biomarkers for predicting the prognoses of EC patients and identifying patients suitable for immunotherapy.

FIGURE 4.

Analysis of immune gene coexpression modules. (A) The low score of the turquoise module indicated a poor prognosis, while the high score of the yellow module indicated a good prognosis. GO analysis showed that (B) the turquoise module was associated with leukocyte migration and myeloid leukocyte migration and was (C) highly negatively correlated with PCA1 in the immune landscape. (D) The yellow module was related to T cell activation and regulation of T cell activation. (E) The yellow module was highly negatively correlated with the PCA1 component in the immune landscape. The difference in prognosis between the groups was not significant (F) in the turquoise module, while (G) the prognosis of the high-scoring group was significantly better than that of the low-scoring group in the yellow module.

FIGURE 4.

Analysis of immune gene coexpression modules. (A) The low score of the turquoise module indicated a poor prognosis, while the high score of the yellow module indicated a good prognosis. GO analysis showed that (B) the turquoise module was associated with leukocyte migration and myeloid leukocyte migration and was (C) highly negatively correlated with PCA1 in the immune landscape. (D) The yellow module was related to T cell activation and regulation of T cell activation. (E) The yellow module was highly negatively correlated with the PCA1 component in the immune landscape. The difference in prognosis between the groups was not significant (F) in the turquoise module, while (G) the prognosis of the high-scoring group was significantly better than that of the low-scoring group in the yellow module.

Close modal

We applied qPCR and immunohistochemistry to detect the differences in the expression of SIRPG and SLAMF1 between 10 low-risk tumor tissues (stage I+II) and 10 high-risk tumor tissues (stage III+IV). The qPCR results showed that the expression levels of SIRPG (Fig. 5A) and SLAMF1 (Fig. 5B) mRNA were low in the high-risk tumor tissues. The protein levels of SIRPG (Fig. 5C) and SLAMF1 (Fig. 5D) were low in the high-risk tumor tissues as detected by immunohistochemistry. The experimental results further validated our bioinformatics analysis.

FIGURE 5.

Expression of SIRPG and SLAMF1 in EC. SIRPG and SLAMF1 mRNA levels in tumor tissues from patients with and without risk. The qPCR results showed that the expression of SIRPG (A) and SLAMF1 (B) mRNA was low in the high-risk tumor tissues. SIRPG and SLAMF1 protein levels in tumor tissues from patients with and without risk. The protein levels of SIRPG (C) and SLAMF1 (D) were low in the high-risk tumor tissues as detected by immunohistochemistry. *p < 0.05.

FIGURE 5.

Expression of SIRPG and SLAMF1 in EC. SIRPG and SLAMF1 mRNA levels in tumor tissues from patients with and without risk. The qPCR results showed that the expression of SIRPG (A) and SLAMF1 (B) mRNA was low in the high-risk tumor tissues. SIRPG and SLAMF1 protein levels in tumor tissues from patients with and without risk. The protein levels of SIRPG (C) and SLAMF1 (D) were low in the high-risk tumor tissues as detected by immunohistochemistry. *p < 0.05.

Close modal

The incidence of EC has been gradually increasing worldwide. Traditional treatments, such as surgery, radiation therapy, and chemotherapy, are applied, but the mortality incidence of EC is still increasing (23, 24). Immunotherapy is an available option to treat a growing number of cancers, including melanoma (25) and lung cancer (26). However, EC patients suffer from low response rates and limited survival benefits after immunotherapy (27), suggesting that there may be a specific pattern associated with immunotherapy benefit.

In this multicohort retrospective study, we identified three ISs of EC, and we used independent data for subtype validation and comprehensive molecular characterization. The results were as follows. We obtained three ISs based on immune-related genes and analyzed the distribution of the TMB among the ISs. The proportion of TP53 mutations in the IS1 subtype was significantly higher than that in the IS2 and IS3 subtypes, which may be related to the poor prognosis. These results are consistent with previous studies (28, 29), thereby demonstrating the feasibility of the IS analysis. Each IS was associated with a different gene expression profile. Regarding chemotherapy-induced immune responses and immune checkpoint–related genes, 22 (84.6%) chemotherapy-induced immune response markers and 43 (91.5%) immune checkpoint–related genes significantly differed among the different ISs. Tumor marker–related genes were also analyzed, and the expression levels of CA125, CA153, CA199, and YKL-40 in the IS1 subtype were significantly higher than those in the IS3 subtype. In addition, we analyzed the response of the different ISs to immunotherapy and chemotherapy. We found that the IS2 subtypes were more sensitive to paclitaxel and cisplatin than the other subtypes, while the IS1 subtypes were more sensitive to doxorubicin and docetaxel. These results show that the ISs have clinical significance for the design and reasonable combination strategy of new immunotherapies. We also analyzed the immune characteristics of the ISs. The poor prognosis of EC may be related to the activation of effector memory CD4 T cells and the inhibition of Th2 cells, immature dendritic cells, eosinophils, mast cells, NKT cells, and plasmacytoid dendritic cells in the IS3 and IS1 subtypes. To visualize and reveal the potential structure of the individual patient distribution, we applied the dimensionality reduction method based on graph learning to the immune gene expression profiles. This analysis placed a single patient into a graph with a sparse tree structure and defined the immune landscape of EC. Furthermore, the Weighted Gene Co-Expression Network Analysis package was used to identify the coexpression modules among the immune genes. We found that the IS2 subtype was significantly positively correlated with the yellow module gene list, which was associated with a good prognosis. Finally, we selected the genes in the yellow module with a correlation >0.9 as hub genes and identified two genes, namely, SIRPG and SLAMF1, as the final characteristic genes. The qPCR and immunohistochemistry results showed that the expression levels of SIRPG and SLAMF1 were low in high-risk tumor tissues.

In this study, we established three ISs and analyzed the immune landscape of EC, which has rarely been mentioned in previous studies. SIRPG and SLAMF1 were identified as new immune-related biomarkers, which may provide a new strategy for immunotherapy for EC. Previous studies have regarded SIRPG as an immune target for improving therapeutic regimens in patients with head and neck squamous cell carcinoma (30) and autoimmune diseases in the thymus (31). SLAMF1 is a hub gene affecting immune infiltration in ovarian cancer (32), leukemias, and lymphomas (33), and it enhances antitumor activity in vitro and in vivo by T cell activation (34). These studies suggest that SIRPG and SLAMF1 are important markers associated with immunotherapy, but these genes have not been studied in EC.

This study had several limitations. The number of samples used for PCR and immunohistochemistry was small, and there was a lack of exploration of the gene function of SIRPG and SLAMF1 in vitro and in vivo. We aim to expand the sample size and further explore the functions of the hub gene in future studies.

In conclusion, this study provided a conceptual framework for understanding the TIM of EC, and our findings suggest that SIRPG and SLAMF1 may be immune-related biomarkers for immunotherapy for EC. These findings may provide new immune characteristics of EC for immunotherapy strategies.

We thank TCGA and GEO databases for the availability of the data.

This work was supported by the National Natural Science Foundation of China (8217113735).

The online version of this article contains supplemental material.

G.Y. and L.W. conceived and designed the study. L.L. and Y.Z. performed the experiments. J.L. and J.Z. analyzed the data and drafted the manuscript. All authors read and approved the final manuscript.

Abbreviations used in this article:

CDF

cumulative distribution function

EC

endometrial cancer

FPKM

fragments per kilobase million

GO

Gene Ontology

IS

immune subtype

MDSC

myeloid-derived suppressor cell

PCA

principal component

PD-1

programmed cell death protein 1

qPCR

quantitative PCR

RNA-seq

RNA sequencing

TCGA

The Cancer Genome Atlas

TIM

tumor immune microenvironment

TMB

tumor mutational burden

UCSC

University of California Santa Cruz

1.
Siegel
R.
,
D.
Naishadham
,
A.
Jemal
.
Cancer statistics.
2013
.
CA Cancer J. Clin.
63
:
11
30
.
2.
Cote
M. L.
,
J. J.
Ruterbusch
,
S. H.
Olson
,
K.
Lu
,
R.
Ali-Fehmi
.
2015
.
The growing burden of endometrial cancer: a major racial disparity affecting Black women.
Cancer Epidemiol. Biomarkers Prev.
24
:
1407
1415
.
3.
Sheikh
M. A.
,
A. D.
Althouse
,
K. E.
Freese
,
S.
Soisson
,
R. P.
Edwards
,
S.
Welburn
,
P.
Sukumvanich
,
J.
Comerci
,
J.
Kelley
,
R. E.
LaPorte
,
F.
Linkov
.
2014
.
USA endometrial cancer projections to 2030: should we be concerned?
Future Oncol.
10
:
2561
2568
.
4.
Duncan
M. E.
,
V.
Seagroatt
,
M. J.
Goldacre
.
2012
.
Cancer of the body of the uterus: trends in mortality and incidence in England, 1985-2008.
BJOG
119
:
333
339
.
5.
Cowie
C. C.
,
K. F.
Rust
,
D. D.
Byrd-Holt
,
E. W.
Gregg
,
E. S.
Ford
,
L. S.
Geiss
,
K. E.
Bainbridge
,
J. E.
Fradkin
.
2010
.
Prevalence of diabetes and high risk for diabetes using A1C criteria in the U.S. population in 1988-2006.
Diabetes Care
33
:
562
568
.
6.
Friberg
E.
,
N.
Orsini
,
C. S.
Mantzoros
,
A.
Wolk
.
2007
.
Diabetes mellitus and risk of endometrial cancer: a meta-analysis.
Diabetologia
50
:
1365
1374
.
7.
Friedenreich
C. M.
,
L. S.
Cook
,
Q.
Wang
,
R. L.
Kokts-Porietis
,
J.
McNeil
,
C.
Ryder-Burbidge
,
K. S.
Courneya
.
2020
.
Prospective cohort study of pre- and postdiagnosis physical activity and endometrial cancer survival. [Published erratum appears in 2021 J. Clin. Oncol. 39: 3650.]
J. Clin. Oncol.
38
:
4107
4117
.
8.
O’Malley
D. M.
,
G. M.
Bariani
,
P. A.
Cassier
,
A.
Marabelle
,
A. R.
Hansen
,
A.
De Jesus Acosta
,
W. H.
Miller
Jr.
,
T.
Safra
,
A.
Italiano
,
L.
Mileshkin
, et al
2022
.
Pembrolizumab in patients with microsatellite instability-high advanced endometrial cancer: results from the KEYNOTE-158 study.
J. Clin. Oncol.
40
:
752
761
.
9.
Marabelle
A.
,
D. T.
Le
,
P. A.
Ascierto
,
A. M.
Di Giacomo
,
A.
De Jesus-Acosta
,
J. P.
Delord
,
R.
Geva
,
M.
Gottfried
,
N.
Penel
,
A. R.
Hansen
, et al
2020
.
Efficacy of pembrolizumab in patients with noncolorectal high microsatellite instability/mismatch repair-deficient cancer: results from the Phase II KEYNOTE-158 study.
J. Clin. Oncol.
38
:
1
10
.
10.
Hanahan
D.
,
L. M.
Coussens
.
2012
.
Accessories to the crime: functions of cells recruited to the tumor microenvironment.
Cancer Cell
21
:
309
322
.
11.
Chen
P.
,
Y.
Yang
,
Y.
Zhang
,
S.
Jiang
,
X.
Li
,
J.
Wan
.
2020
.
Identification of prognostic immune-related genes in the tumor microenvironment of endometrial cancer.
Aging (Albany NY)
12
:
3371
3387
.
12.
Lei
X.
,
Y.
Lei
,
J. K.
Li
,
W. X.
Du
,
R. G.
Li
,
J.
Yang
,
J.
Li
,
F.
Li
,
H. B.
Tan
.
2020
.
Immune cells within the tumor microenvironment: biological functions and roles in cancer immunotherapy.
Cancer Lett.
470
:
126
133
.
13.
Kobayashi
S. D.
,
F. R.
Deleo
.
2013
.
Systems biology and innate immunity.
J. Innate Immun.
5
:
97
99
.
14.
Huang
X.
,
G.
Zhang
,
T.
Tang
,
T.
Liang
.
2021
.
Identification of tumor antigens and ISs of pancreatic adenocarcinoma for mRNA vaccine development.
Mol. Cancer
20
:
44
.
15.
Wilkerson
M. D.
,
D. N.
Hayes
.
2010
.
ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking.
Bioinformatics
26
:
1572
1573
.
16.
Thorsson
V.
,
D. L.
Gibbs
,
S. D.
Brown
,
D.
Wolf
,
D. S.
Bortone
,
T. H.
Ou Yang
,
E.
Porta-Pardo
,
G. F.
Gao
,
C. L.
Plaisier
,
J. A.
Eddy
, et al
Cancer Genome Atlas Research Network
.
2018
.
The immune landscape of cancer. [Published erratum appears in 2019 Immunity 51: 411–412.]
Immunity
48
:
812
830.e14
.
17.
Zhang
H.
,
Y.
Chen
.
2021
.
Identification of glioblastoma ISs and immune landscape based on a large cohort.
Hereditas
158
:
30
.
18.
Wang
L.
,
Q.
Mao
.
2019
.
Probabilistic dimensionality reduction via structure learning.
IEEE Trans. Pattern Anal. Mach. Intell.
41
:
205
219
.
19.
Trapnell
C.
,
D.
Cacchiarelli
,
J.
Grimsby
,
P.
Pokharel
,
S.
Li
,
M.
Morse
,
N. J.
Lennon
,
K. J.
Livak
,
T. S.
Mikkelsen
,
J. L.
Rinn
.
2014
.
The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells.
Nat. Biotechnol.
32
:
381
386
.
20.
Liang
L.
,
J.
Yu
,
J.
Li
,
N.
Li
,
J.
Liu
,
L.
Xiu
,
J.
Zeng
,
T.
Wang
,
L.
Wu
.
2021
.
Integration of scRNA-seq and bulk RNA-seq to analyse the heterogeneity of ovarian cancer immune cells and establish a molecular risk model.
Front. Oncol.
11
:
711020
.
21.
Danilova
L.
,
W. J.
Ho
,
Q.
Zhu
,
T.
Vithayathil
,
A.
De Jesus-Acosta
,
N. S.
Azad
,
D. A.
Laheru
,
E. J.
Fertig
,
R.
Anders
,
E. M.
Jaffee
,
M.
Yarchoan
.
2019
.
Programmed Cell Death Ligand-1 (PD-L1) and CD8 expression profiling identify an immunologic subtype of pancreatic ductal adenocarcinomas with favorable survival.
Cancer Immunol. Res.
7
:
886
895
.
22.
Charoentong
P.
,
F.
Finotello
,
M.
Angelova
,
C.
Mayer
,
M.
Efremova
,
D.
Rieder
,
H.
Hackl
,
Z.
Trajanoski
.
2017
.
Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade.
Cell Rep.
18
:
248
262
.
23.
Makker
V.
,
M. H.
Taylor
,
C.
Aghajanian
,
A.
Oaknin
,
J.
Mier
,
A. L.
Cohn
,
M.
Romeo
,
R.
Bratos
,
M. S.
Brose
,
C.
DiSimone
, et al
2020
.
Lenvatinib Plus pembrolizumab in patients with advanced endometrial cancer.
J. Clin. Oncol.
38
:
2981
2992
.
24.
Soslow
R. A.
,
C.
Tornos
,
K. J.
Park
,
A.
Malpica
,
X.
Matias-Guiu
,
E.
Oliva
,
V.
Parkash
,
J.
Carlson
,
W. G.
McCluggage
,
C. B.
Gilks
.
2019
.
Endometrial carcinoma diagnosis: use of FIGO grading and genomic subcategories in clinical practice: recommendations of the International Society of Gynecological Pathologists.
Int. J. Gynecol. Pathol.
38
(
Suppl. 1
):
S64
S74
.
25.
Harel
M.
,
R.
Ortenberg
,
S. K.
Varanasi
,
K. C.
Mangalhara
,
M.
Mardamshina
,
E.
Markovits
,
E. N.
Baruch
,
V.
Tripple
,
M.
Arama-Chayoth
,
E.
Greenberg
, et al
2019
.
Proteomics of melanoma response to immunotherapy reveals mitochondrial dependence.
Cell
179
:
236
250.e18
.
26.
Alessi
J. V.
,
M. M.
Awad
.
2020
.
Immunotherapy in lung cancer: effective for patients with poor performance status?
Lancet Respir. Med.
8
:
838
839
.
27.
Di Tucci
C.
,
C.
Capone
,
G.
Galati
,
V.
Iacobelli
,
M. C.
Schiavi
,
V.
Di Donato
,
L.
Muzii
,
P. B.
Panici
.
2019
.
Immunotherapy in endometrial cancer: new scenarios on the horizon.
J. Gynecol. Oncol.
30
:
e46
.
28.
Sherman
M. E.
,
M. E.
Bur
,
R. J.
Kurman
.
1995
.
p53 in endometrial cancer and its putative precursors: evidence for diverse pathways of tumorigenesis.
Hum. Pathol.
26
:
1268
1274
.
29.
Levine
R. L.
,
C. B.
Cargile
,
M. S.
Blazes
,
B.
van Rees
,
R. J.
Kurman
,
L. H.
Ellenson
.
1998
.
PTEN mutations and microsatellite instability in complex atypical hyperplasia, a precursor lesion to uterine endometrioid carcinoma.
Cancer Res.
58
:
3254
3258
.
30.
Wang
J.
,
Y.
Tian
,
G.
Zhu
,
Z.
Li
,
Z.
Wu
,
G.
Wei
,
L.
Zhuang
,
Z.
Li
,
X.
Chen
,
X.
Zhang
, et al
2021
.
Establishment and validation of immune microenvironmental gene signatures for predicting prognosis in patients with head and neck squamous cell carcinoma.
Int. Immunopharmacol.
97
:
107817
.
31.
Gabrielsen
I. S.
,
S. S.
Amundsen
,
H.
Helgeland
,
S. T.
Flåm
,
N.
Hatinoor
,
K.
Holm
,
M. K.
Viken
,
B. A.
Lie
.
2016
.
Genetic risk variants for autoimmune diseases that influence gene expression in thymus.
Hum. Mol. Genet.
25
:
3117
3124
.
32.
Su
R.
,
C.
Jin
,
L.
Zhou
,
Y.
Cao
,
M.
Kuang
,
L.
Li
,
J.
Xiang
.
2021
.
Construction of a ceRNA network of hub genes affecting immune infiltration in ovarian cancer identified by WGCNA.
BMC Cancer
21
:
970
.
33.
Gomes
A. Q.
,
D. V.
Correia
,
A. R.
Grosso
,
T.
Lança
,
C.
Ferreira
,
J. F.
Lacerda
,
J. T.
Barata
,
M. G.
Silva
,
B.
Silva-Santos
.
2010
.
Identification of a panel of ten cell surface protein antigens associated with immunotargeting of leukemias and lymphomas by peripheral blood gammadelta T cells.
Haematologica
95
:
1397
1404
.
34.
Mehrle
S.
,
J.
Schmidt
,
M. W.
Büchler
,
C.
Watzl
,
A.
Märten
.
2008
.
Enhancement of anti-tumor activity in vitro and in vivo by CD150 and SAP.
Mol. Immunol.
45
:
796
804
.

The authors have no financial conflicts of interest.

This article is distributed under The American Association of Immunologists, Inc., Reuse Terms and Conditions for Author Choice articles.

Supplementary data