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 (4–6). 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.
Materials and Methods
Expression profile data extraction
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.
Sources of immune-related genes
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.
TCGA data preprocessing
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.
Identification of ISs and immune gene modules and functional analysis
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.
Assessment of immune-related clinical, molecular, and cellular features
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.
Construction of the immune landscape
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 (17–19). 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).
Real-time quantitative PCR and immunohistochemistry
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).
Ethics approval and consent to participate
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.
Availability of data and materials
The datasets used or analyzed during this study are available from the corresponding author on reasonable request.
Molecular typing based on immune-related genes
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.
Relationship between ISs and the tumor mutational burden and common gene mutations
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.
Expression of classic markers of the chemotherapy-induced immune response and expression of immune checkpoint genes in ISs
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.
Expression differences in tumor marker–related genes among the three ISs
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.
Sensitivity of ISs to immunotherapy and chemotherapy
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.
Immune characteristics in different ISs
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.
Immune landscape of EC
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.
Identification of immune gene coexpression modules
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).
Function and prognosis analyses of immune gene coexpression modules
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.
Expression of SIRPG and SLAMF1 as new immune-related biomarkers in EC
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.
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:
cumulative distribution function
fragments per kilobase million
myeloid-derived suppressor cell
programmed cell death protein 1
The Cancer Genome Atlas
tumor immune microenvironment
tumor mutational burden
University of California Santa Cruz
The authors have no financial conflicts of interest.