Exploration of a Robust and Prognostic Immune Related Gene Signature for Cervical Squamous Cell Carcinoma

Abstract

Background: Cervical squamous cell carcinoma (CESC) is one of the most frequent malignancies in women worldwide. The level of immune cell infiltration and immune-related genes (IRGs) can significantly affect the prognosis and immunotherapy of CESC patients. Thus, this study aimed to identify an immune-related prognostic signature for CESC.

Methods: TCGA-CESC cohorts, obtained from TCGA database, were divided into the training group and testing group; while GSE44001 dataset from GEO database was viewed as external validation group. ESTIMATE algorithm was applied to evaluate the infiltration levels of immune cells of CESC patients. IRGs were screened out through weighted gene co-expression network analysis (WGCNA). A multi-gene prognostic signature based on IRGs was constructed using LASSO penalized Cox proportional hazards regression, which was validated through Kaplan–Meier, Cox, and receiver operating characteristic curve (ROC) analyses. The abundance of immune cells was calculated using ssGSEA algorithm in the ImmuCellAI database, and the response to immunotherapy was evaluated using immunophenoscore (IPS) analysis and the TIDE algorithm.

Results: In TCGA-CESC cohorts, higher levels of immune cell infiltration were closely associated with better prognoses. Moreover, a prognostic signature was constructed using three IRGs. Based on this given signature, Kaplan–Meier analysis suggested the significant differences in overall survival (OS) and the ROC analysis demonstrated its robust predictive potential for CESC prognosis, further confirmed by internal and external validation. Additionally, multivariate Cox analysis revealed that the three IRGs signature served as an independent prognostic factor for CESC. In the three-IRGs signature low-risk group, the infiltrating immune cells (B cells, CD4/8 + T cells, cytotoxic T cells, macrophages and so on) were much more abundant than that in high-risk group. Ultimately, IPS and TIDE analyses showed that low-risk CESC patients appeared to present with a better response to immunotherapy and a better prognosis than high-risk patients.

Conclusion: The present prognostic signature based on three IRGs (CD3E, CD3D, LCK) was not only reliable for survival prediction but efficient to predict the clinical response to immunotherapy for CESC patients, which might assist in guiding more precise individual treatment in the future.

Introduction

Cervical cancer ranks as the fourth most frequent malignancy and the fourth leading cause of cancer-related death in women worldwide, with approximately 570,000 new cases and 311,000 deaths reported by the International Agency for Research on Cancer (IARC) in 2018 alone (Bray et al., 2018). Cervical squamous cell carcinoma (CESC) comprises the primary histological type of cervical cancer, accounting for more than 80% of cases. Although recent advances in comprehensive treatments, including surgical resection, chemotherapy, and radiotherapy, have been made, the overall mortality has increased annually around the world (Jemal et al., 2011; Torre et al., 2015; Bray et al., 2018). For instance, in China, the number of deaths increased at an annual percentage change of 5.9% from 2001 to 2011 (Chen et al., 2016). Notably, the high risk of metastasis and recurrence makes it challenging for CESC patients with traditional therapies to maintain a good prognosis; thus, the development of a reliable prognostic assessment and novel therapeutic strategies is urgently needed (Pfaendler and Tewari, 2016).

In recent years, immunotherapy, especially the use of immune checkpoint inhibitors (ICIs), has made appreciable progress in antitumor practice and is becoming the new first-line therapeutic option in oncology (Fillon, 2020); e.g., for treatment of bladder cancer (Patel et al., 2020), melanoma (Zaretsky et al., 2016), renal cancer (Curti, 2018). Distinct from traditional treatments, the clinical benefits for immunotherapy patients are achieved by stimulating the persistent antitumor immune response (Kennedy and Salama, 2020), which depends on immunomodulation between the tumor microenvironment (TME) and cancer cells. Numerous studies have reported that infiltrating immune cells, a vital part of the TME, are closely related to cancer progression and the efficacy of immunotherapy (Hanahan and Coussens, 2012; Byrne et al., 2020). For instance, the mechanism of immune checkpoint blockade (ICB) therapy is to enhance a patient’s anti-tumor immune response by blocking the inhibitory effect of tumor cells on immune cells. In this process, programmed cell death 1 (PD1) and cytotoxic T lymphocyte antigen 4 (CTLA-4) are the most common targets of ICIs, which generally inhibit the activation and amplification of T cells that render the anti-tumor response ineffective. Although ICIs were expected to exhibit great potential in the immunotherapy of CESC patients, the clinical outcomes and prognosis were far from satisfactory, as only less than 20% patients achieved a partial or complete response where most had a PD-L1-positive tumors (Frenel et al., 2017; Chung et al., 2019). Meanwhile, the heterogeneity of increased levels of tumor infiltrating lymphocytes and tumor mutational burden in CESC participants also gives the additional explanation for it (Otter et al., 2019). Therefore, to achieve precise personalized decision-making for ICI treatment, a robust prognostic signature will be important to determine the prognosis and predict the sensitivity of immunotherapy in CESC patients.

Although several risk models have been developed to predict the prognosis of CESC patients (Yang et al., 2019; Liu et al., 2020), these have been based on differentially expressed genes, which ignored the complicated interactions of genes, as well as hub genes that showed no significant differences between tumor and normal tissues, but were closely linked to clinical features. Additionally, a promising prognostic signature based on immune-related genes (IRGs) has been confirmed great potential to predict prognosis and immunotherapy responsiveness in cancers but were less reported in CESC. In the present study, weighted gene co-expression network analysis (WGCNA) was the first time to be applied to identify significant IRGs for CESC patients, which were highly associated with infiltration levels of immune cells calculated by the use of ESTIMATE algorithm. Furthermore, hub IRGs were selected to construct a promising prognostic signature through the least absolute shrinkage and selection operator (LASSO) penalized Cox proportional hazards regression. Ultimately, internal and external validations further demonstrated that this robust prognostic signature based on IRGs in this study was not only reliable for survival prediction but efficient to predict the clinical response of ICIs for CESC patients, which might facilitate personalized counseling for immunotherapy in the future.

Methods

Data Acquisition and Preprocessing

The RNA-seq datasets (FPKM profiles) and clinical information of the CESC cohorts were retrieved from The Cancer Genome Atlas (TCGA) website (https://portal.gdc.cancer.gov/). Then, the gene expression matrix was formed via Perl script based on the CMD command. IDs were converted to gene symbols according to the Ensembl database (http://asia.ensembl.org/index.html). Meanwhile, mRNA expression profiles (GSE44001) were also downloaded from Gene Expression Omnibus (GEO) database to be viewed as the external validation group, where the matrix profile was conducted by Robust Multi-Array Average (RMA) algorithm in R affyPLM package (Irizarry et al., 2003). But of note, any patients with the follow-up time less than 1 day were excluded. In addition, if matched with multiple IDs in two matrix profiles, the gene expression values were averaged. Finally, log2 processing of the matrix data was performed for further analyses.

Evaluation of Tumor Microenvironment in The Cancer Genome Atlas-Cervical Squamous Cell Carcinoma Cohort

ESTIMATE algorithm (Yoshihara et al., 2013) has been popularly applied to predict infiltration levels of immune cells and stromal cells in the TME based on specific gene expression, of which the final results were evaluated using the immune score and stromal score, respectively. The ESTIMATE score was also calculated to measure the overall tumor immunity of CESC. Then, all patients were divided into two groups (high score and low score) based on the mean score, and survival analyses were performed to explore whether the overall survival were related with the immune infiltration levels.

Construction of Weighted Gene Co-expression Networks

In the present study, the top 25% most-varying genes calculated by the mean absolute deviation (MAD) algorithm were selected for subsequent analysis. The outlier samples were firstly identified and removed using the goodSamplesGenes function in the WGCNA package. Next, Pearson correlation coefficients were calculated between any two genes to construct the gene expression similarity matrix with the following formula: , where X_i and X_j represent the expression values of gene i and j. Furthermore, the cut-off value of scale-free R2 was set at 0.9 to obtain the lowest soft-thresholding power β that would be used to build an adjacency matrix so that gene distribution conformed to scale-free networks. Based on the topological overlap matrix (TOM) converted from the adjacency matrix, dissimilarity between genes was employed to cluster genes with similar expression profiles by the hierarchical clustering method, which were then cut into different modules by dynamic shear module recognition and visualized by the dendrogram with colored assignments (minimum module size was 80). Meanwhile, a cut-height of 0.35 was considered as the cut-off value to merge similar modules.

Identification of Significant Modules and Genes Related to Clinical Features

The clinical traits of each sample were combined with module eigengenes (MEs) to construct the relevant clinical modules. MEs were defined as the first principal component of the module, representing the gene expression profile of the entire module. The clustering modules that were the most closely associated with immune cell infiltration were selected as the significant modules. Gene significance (GS) and module membership (MM) were also calculated to screen significant genes. GS represents the correlation between gene expression and clinical traits, and MM reflects the correlation between gene expression profile and genes within given modules. Ultimately, GS > 0.4 and MM > 0.8 were set as the criteria to identify hub genes that were strongly related to the clinically significant traits. (Zhang and Horvath, 2005).

To further screen potential genes that play an essential role in immune cell infiltration of CESC, a total of 1811 IRGs were obtained from the ImmPort database (https://immport.niaid.nih.gov) (Bhattacharya et al., 2014), and these overlapping genes were selected as hub IRGs for subsequent analysis.

Functional Enrichment and Interactions Analyses

The functional enrichment analyses including Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed using the R clusterProfiler package (Yu et al., 2012). GO terms were further divided into biological process (BP), molecular function (MF), and cellular component (CC). Adj. p. value <0.05 were considered significant. To identify the interactive relationships of hub genes, a protein–protein interaction (PPI) network was constructed using the STRING database (https://string-db.org/) with the minimum required interaction score being 0.7. In addition, the R corrplot package was applied to conduct Pearson correlation analysis between hub IRGs.

Development and Validation of the Prognostic Signature

CESC patients were randomly divided into training and test groups with a 7:3 ratio (Table 1). In the training set, univariate Cox regression analysis was performed to explore the relationship between the expression of each key IRG and overall survival (OS). Then, these genes with a p-value < 0.05 were selected as hub IRGs, which were further analyzed by using LASSO penalized Cox proportional hazards regression to find the best risk model in the R package “glmnet” (Wang et al., 2019a). The risk score was calculated using the following formula: , where β is the coefficient of each prognostic hub gene, G represents the expression value of hub genes, and n denotes the number of hub genes. Patients were classified into low-risk (<mean risk score) and high-risk (>mean risk score) groups depending on the mean risk score. Moreover, the survival curve was determined using the Kaplan–Meier method in the R survminer package, where the differences between low- and high-risk groups were calculated by the log-rank test. Meanwhile, a time-dependent receiver operating characteristic (ROC) curve was adopted using R ROC package (Kamarudin et al., 2017), of which the area under the curve (AUC) was calculated to assess the accuracy of the prognostic risk model.

To further verify the predictive performance of the prognostic model, the risk scores were also calculated in the testing group and external groups using the same prognostic formula, and the Kaplan–Meier survival curve and ROC curve were generated with a cutoff value of the mean risk score.

Immune Infiltration Patterns with Prognostic Signature

Immune infiltrating cells, such as T cells, B cells, and monocytes, have been shown to be important in the TME and to significantly affect cancer progression. Thus, ImmuCellAI (http://bioinfo.life.hust.edu.cn/web/ImmuCellAI/) (Miao et al., 2020), a novel and efficient tool based on the ssGSEA algorithm that can estimate the abundance of 24 immune infiltrate cells (18 T-cell subtypes and B cells, natural killer cells, monocytes, macrophages, neutrophils, and dendritic cells) from the gene expression data set (Bindea et al., 2013), was applied to calculate the abundance of immune cells and to compare the differences in infiltrated patterns between low- and high-risk groups. Meanwhile, we also employed TIMER database to explore the correlation of the expression levels of hub IRGs with the expression of immune checkpoint molecules and infiltrating levels of immune cells (Li et al., 2017).

Immunogenicity and Immunotherapeutic Sensitivity with Prognostic Signature

To further validate the predictive performance of the given prognostic signature for the response to ICIs, two independent methods including immunophenoscore (IPS) analysis and the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm were employed to assess the immunogenicity and immunotherapeutic sensitivity of CESC patients, respectively.

The immunogenicity of a patient was determined by four main parts (effector cells, immunosuppressive cells, MHC molecules, and immunomodulators), which could be calculated without bias using machine learning methods by IPS analysis (Charoentong et al., 2017). Higher IPS scores (range 0–10) represent increased immunogenicity. The IPSs of patients with CESC were obtained from The Cancer Immunome Atlas (TCIA) (https://tcia.at/home).

The TIDE algorithm, designed by Dr. Shirley Liu and colleagues, was considered as highly reliable method to predict the clinical response of patients to ICB therapy (PD1 and CTLA-4 inhibitor) in recent studies (Jiang et al., 2018). The results were measured by TIDE score. According to the default settings, a patient with a TIDE value <0 was defined as a responder (positive sensitivity to immunotherapy), whereas a patient with a TIDE value >0 was defined as a non-responder (negative sensitivity to immunotherapy). Furthermore, a correlation analysis of CESC samples was performed to explore the correlation between the given prognostic signature and immune function.

Statistical Analysis

The WGCNA method and all statistical analyses were performed in R software (3.6.1) and GraphPad Prism (8.0). The Mann–Whitney test was applied to compare differences of continuous data between two groups, whereas ANOVA was used for more than two groups. A Chi square test was used to test for differences between categorical variables. p values <0.05 were considered to be statistically significant.

Results

Correlation of ESTIMATE Score and Clinical Characteristics

The flow diagram of the present study is shown in Figure 1. A total of 291 and 301 CESC samples were enrolled from TCGA and GEO database according to the inclusion criteria, respectively. The ESTIMATE score of each patient, which reflects the landscape of the TME and overall immune-infiltration degree was calculated with the ESTIMATE algorithm. As shown in Figure 2, patients with high immune and ESTIMATE scores demonstrated a better OS than those with low scores (p < 0.05). Nevertheless, all scores showed no statistical significance between different race, age, and International Federation of Gynecology and Obstetrics (FIGO) stage (Supplementary Figure S1 and Supplementary Table S1).

FIGURE 1

FIGURE 2

Identification of Immune-Related Modules and Hub Genes by Weighted Gene Co-Expression Network Analysis

First, three outlier samples were excluded by the goodSamplesGenes function. The dendrogram and heatmap (Supplementary Figure S2) display the distribution and clinical traits of the remaining 288 samples in detail. The soft-threshold power β = 4, the first parameter to meet the requirements of scale-free R2 > 0.9, was selected to construct the scale-free networks (Supplementary Figure S3). Furthermore, after similar modules were merged with a cut-height of 0.35, a total of 21 modules was ultimately identified for further analyses (Figure 3). Interestingly, the green and dark turquoise modules were found to be most associated with immune cell infiltration levels of CESC in the TME (Figure 4A), in which the correlation coefficients were 0.97 and 0.85, respectively. Therefore, the green (Figure 4B) and dark turquoise (Figure 4C) modules were considered immune-related key modules, and significant genes were selected with cut-off criteria of GS > 0.4 and MM > 0.8.

FIGURE 3

FIGURE 4

Function Analyses of Immune-Related Modules and Hub Genes

As shown in Figure 5A, most genes within immune-related modules were highly enriched in the inflammatory response, immune response, and proteolysis in the BP category, while chemokine activity, immunoglobulin receptor binding, and antigen binding; and plasma membrane, extracellular exosome, and cytosol were enriched in the MF and CC categories, respectively. Moreover, KEGG pathway analysis (Figure 5B) indicated that virus infection, cytokine−cytokine receptor interaction, antigen processing and presentation, and chemokine signaling pathway were significant. Therefore, the above results suggested that most defined genes in immune modules positively participated in immune-related biological processes in CESC.

FIGURE 5

Establishment and Validation of the Prognostic Signature

A total of 31 common genes were identified between significant module genes and IRGs obtained from the ImmPort database (Figure 6A). Subsequently, TCGA-CESC datasets were randomly divided into a training set (n = 201) and testing set (n = 87); basic clinical information is provided in Table 1. Using univariate Cox regression analysis in the training group, 20 genes were correlated with CESC patient survival (p < 0.05) (Table 2). The PPI network is shown in Figure 6B, and the results of Pearson correlation analysis indicated that 20 hub IRGs were significantly co-expressed in CESC samples (Figure 6C). GO analysis (Figure 6D) illustrated that the biological process terms were highly focused on T cell activation, positive regulation of leukocyte activation, leukocyte cell−cell adhesion, lymphocyte differentiation, and T cell differentiation, which further confirmed the important role of identified hub IRGs in the activities of immune infiltration and response.

FIGURE 6

After the optimal model of prognostic prediction, the individualized risk scores were calculated with coefficient values extracted by LASSO Cox regression analysis. The formula was as follows: risk score = expression of CD3E*(−0.11337) + CD3D*(−0.01026) + LCK*(−0.0523). The LASSO plot is shown in Supplementary Figure S4, and the distribution of patient risk scores and survival status is shown in Figure 7A. The survival analysis indicated that the prognosis of high-risk patients was significantly worse than that of low-risk patients (Figure 7B). Additionally, the time-dependent ROC showed that the AUC values of 1-, 3-, and 5-year OS were 0.705, 0.641, and 0.662, respectively (Figure 7C). Ultimately, the 3-IRGs prognostic signature was validated using OS data from the testing set, of which the results remained consistent (Figure 8D), with the AUC values of 1-, 3-, and 5-year OS of 0.767, 0.770, and 0.702, respectively. Interestingly, the AUC values of 1-, 3-, and 5-year OS were 0.651, 0.648, 0.612 in external group (Figure 9), indicated the reliable predictive potential of the given prognostic signature. Meanwhile, the FIGO stage and prognostic model (low/high) were confirmed as independent risk factors for survival of CESC patients by the univariate and multivariate Cox regression analyses (Table 3).

FIGURE 7

FIGURE 8

FIGURE 9

Immune Infiltration Patterns with Prognostic Signature

As shown in Figure 10, the expression levels of CD3D, CD3E and LCK were significantly positively correlation with the expression of three immune checkpoint molecules, including PD1(PDCD1), PDL1(CD274), and CTLA4, respectively. In addition, we also found that the CD3D, CD3E and LCK expression was significantly negatively related to tumor purity while positively associated with the infiltrating levels of immune cells, such as B cell, CD8 + T cell, and CD4 + T cell (Supplementary Figure S5). Thus, these results suggested that CESC patients with higher expression of three IRGs may had a better response to immunotherapy.

FIGURE 10

The correlation of prognostic signature and ESTIMATE score was shown in Supplementary Figure S6. To further investigate whether there were significant differences in immune cell infiltration between low- and high-risk groups, the abundance of 24 immune cells of each patient was determined using the ssGSEA algorithm in the ImmunCellAI database. As shown in Figure 11, the abundance of immune infiltrating cells, such as CD4/8 naïve T cells, neutrophils, monocytes, and NK T cells, was increased (p < 0.05) in the high risk group compared to the low risk group; whereas the infiltrating levels of B cells, CD4/8+ T cells, Th1/2 cells, iTregs, macrophages, and Cytotoxic T cells in the low-risk group were markedly higher than those of the high-risk group. These results are shown in detail in Supplementary Table S2.

FIGURE 11

Immunogenicity and Immunotherapeutic Sensitivity with Prognostic Signature

In this work, the IPS of each CESC patient was evaluated to explore the association between immunogenicity and the two prognostic risk subgroups. As shown in Figures 12A–D, the IPS, IPS-CTLA4, IPS-PD1-PD-L1-PD-L2, and IPS-PD1-PD-L1-PD-L2-CTLA4 scores in the low-risk group were significantly higher than those of the high-risk group (p < 0.05), indicating a more immunogenic phenotype in the low-risk group. In addition, according to the TIDE algorithm, the response to ICIs of CESC patients was measured by the TIDE value (Supplementary Figure S7). The results also indicated that risk scores were significantly negatively correlated with IFNG, CD274 (PD-L1), CD8+ T cells, and dysfunction, and positively correlated with exclusion and MDSC (Figure 12E), which further confirmed our findings; the details are provided in Supplementary Table S3. Moreover, the number of CESC patients that exhibited a positive response to ICIs in the low-risk group was higher than that in the high-risk group (p < 0.05) (Figure 12F). The Kaplan–Meier survival curve showed improved survival in the responder group compared with the non-responder group (p = 0.036) (Figure 12G). Overall, the results demonstrate that the low-risk group identified by the three IRGs appeared to present with a more positive response to immunotherapy and a better prognosis.

FIGURE 12

Discussion

Although advanced diagnostic methods and molecular anticancer therapies have been rapidly developed, the overall prognosis of CESC patients remains poor (McLachlan et al., 2017). ICB therapy was expected to be a breakthrough in CESC treatment (Hamanishi et al., 2017). In several early-phase randomized clinical trials, the application of pembrolizumab and nivolumab (PD-L1/PD-1 blockage) suggested promising clinical outcomes for CESC patients with metastasis and recurrence (Frenel et al., 2017; Chung et al., 2019). However, only a small portion of CESC patients showed a positive response to ICI treatments. Therefore, we identified a robust prognostic signature based on IRGs to forecast prognosis and immunotherapy sensitivity for CESC patients, which may facilitate personalized counseling.

In the present study, the abundance of immune cells and stromal cells in the TME was first investigated using the ESTIMATE algorithm, which was considered a clinical trait for further analysis. The results also indicated that the high infiltration of immune cells showed a better prognosis for CESC patients, indicating that the TME plays an essential role in the prognosis of CESC patients. GO enrichment and KEGG pathway analyses demonstrated that most genes identified as closely related to immune cell infiltration by WGCNA were enriched in immunomodulatory activities, such as the inflammatory response, immune response, and chemokine signaling pathway. By performing univariate and LASSO Cox regression analyses, three hub IRGs (CD3E, CD3D, and LCK) were selected to establish an immune-related prognostic signature for CESC patients, where the AUC values of 1-, 3-, and 5-year OS were 0.705, 0.641, and 0.662 in the training group, and 0.767, 0.770, and 0.702 in the testing group and 0.651, 0.648, 0.612 in the external validation, respectively, indicating the reliably predictive capacity for CESC patient prognosis.

Additionally, functional annotation further suggested that CD3E, CD3D, and LCK were involved in positive regulation of T cell activation and leukocyte cell-cell adhesion that were known as the chief determinant for the efficacy of tumor immunotherapy (Ding and Chen, 2019; Alvarez et al., 2020). Meanwhile, existing evidence has accumulated demonstrating the important role of three hub IRGs in the regulation of immune responses of tumor tissues. The promising tumor immunotherapy mainly depends on the recognition of T-cell receptor (TCR) to special tumor antigens to stimulate the activation of self T cells in order to attack cancer cells (Kennedy and Salama, 2020). Notably, the CD3 co-receptor complex is vital for signal transduction after specific binding of TCR, of which the integrity is considered as the crucial factor for cytotoxic T cell responses to tumor antigens (Fuehrer et al., 2014). CD3E usually participates in encoding the CD3ε chain, one of the major components (γ-, δ-, ε- and ζ-chain) of the CD3 co-receptor complex, whose deficiency will cause the severe combined immunodeficiency (Firtina et al., 2017; Erman et al., 2020). Hart et al. has reported that a reduced cell surface abundance of CD3E could lead to a significant inhibition of T cell killing capacity (Hart et al., 2019). In contrast, the increased expression of CD3E was found markedly related to the effective response in 31 cancer types patients who received anti-PD1 immunotherapy (Gaffney et al., 2019). CD3D has been reported as a potential biomarker for the response to ICIs and prognosis in cancers, including colon cancer and muscle-invasive bladder cancer (Klintman et al., 2016; Shi et al., 2019; Yang et al., 2020). Homozygous mutations in the CD3D gene can lead to markedly abnormal T-cell development, and thus, to early-onset immunodeficiency (Fischer et al., 2005; Gil et al., 2011). Moreover, high expression of CD3E and CD3D gene were reported significantly related to positive OS in CESC (Wang et al., 2019b). LCK, also known as lymphocyte-specific protein tyrosine kinase p56, was found as a key molecule in T cell activation by phosphorylating the TCR/CD3 complex to initiate TCR signaling (Wei et al., 2020). Recent studies showed that the improved LCK activity was contributed to improve the efficacy of chimeric antigen receptors (CARs) immunotherapy in cancers (Gulati et al., 2018; Bommhardt et al., 2019; Suryadevara et al., 2019). In contrast, the inhibited targeted drugs of LCK was reported that could cause the loss of T-cell immune response and result in immunosuppression for patients cancer (Zhao et al., 2008). Therefore, these defined three hub IRGs were expected as immunotherapeutic biomarkers and potential therapeutic targets for CESC patients; this will be the subject of future studies.

The number and proportion of infiltrating immune cells in the TME are recognized as important factors affecting cancer progression and immunotherapy response. To further elucidate the role of the TME associated with this given prognostic signature, the ImmunCellAI database based on ssGSEA algorithm was the first time to be employed to analyze the immune cell infiltrating landscape in CESC patients. We found that most immune cells in low risk groups were more abundant than in high risk group, such as cytotoxic T cell, exhausted T cells, Th1/2 cells, and NK cells, suggesting a more potent immune defense in low risk CESC patients. Additionally, the increased B cells, CD4+ T cell, and CD8+ T cell infiltration indicated that better therapeutic outcomes may be achieved with ICI immunotherapy in low risk groups compared to high risk groups (Liu et al., 2019b; Matsuzaki et al., 2019). Dr. Joy et al. (Hsu et al., 2018) reported that NK cells, in addition to T cells, could enhance the effect of ICI immunotherapy, whereas that efficacy may be weakened due to competition between NK cells and T cells shown in a recent study (Alvarez et al., 2020). In contrast, Tregs and tumor-related macrophages have been reported to be immunosuppressive cells that can form an immunosuppressive atmosphere to facilitate tumor progression by disrupting the adaptive immune response (Milowsky et al., 2016).

Ultimately, to further explore the predictive value of the prognostic signature in ICI immunotherapy, two independent methods, IPS analysis and TIDE algorithm, were applied to calculate the response of CESC patients to ICIs. The results showed that IPS, IPS-CTLA4, IPS-PD1/PD-L1/PD-L2, and IPS/PD1/PD-L1/PD-L2 + CTLA4 scores were significantly increased in this prognostic signature low-risk group. And TIDE algorithm indicated that low risk patients appeared to present with more positive response to anti-ICIs immunotherapy. Both of these founding further support the potential of this given immune-related prognostic signature to determine the immunotherapy sensitivity for CESC patients.

This study represents the first application of WGCNA to identify hub IRGs linked to immune cell infiltration in an effort to develop a prognostic signature for predicting CESC patients prognosis. And it is the first study that employs ImmunCellAI and TIDE algorithm to analyze the immune cell infiltrating landscape and predict immunotherapy sensitivity for CESC, respectively. Compared with individual biomarkers (e.g., PD1 and PDL1) that were more susceptible to interference, the multiple genes signature showed a more reliable stability. In addition, the performance of our three IRGs signature to predict the progression and immune response was better than that in previous study (Liu et al., 2020). Importantly, the versatility of this prognostic model was further verified using external validation, which was less reported in previous study. Nevertheless, some limitations in this work still remain. Firstly, the relative values of gene expression and the difference of FIGO stages between TCGA and GEO database may contribute to the discrepancies in verification groups (e.g., the limited performance in external validation) or in further clinical trails. Secondly, the prognostic signature of CESC patients was developed based on the OS rate, but some external factors such as the TMN degrees were not extensively evaluated due to a lack of relevant data. Thirdly, although bioinformatics tools are helpful to display the interaction of hub IRGs, the external experiments are also important to further elucidate the molecular mechanisms. Finally, TCGA-CESC patients have not received relevant immunotherapy, and the response to ICI treatment was calculated by cutting-edge bioinformatics technologies. Although the potential of this prognostic signature to stratify CESC patients with different immune response was verified reliable by the consistent results of two independent and powerful approaches, a multicenter and large-scale study is still needed to evaluate its practicality in clinical tests and to strengthen its the clinical evidence.

Conclusion

In the present study, we have identified hub genes related to the immune infiltration in CESC microenvironment and constructed a robust three IRGs (CD3E, CD3D, and LCK) signature to predict the prognosis of CESC patients. Meanwhile, the versatility of the signature was verified by using internal and external validation groups. In the further exploration, we also found that this model had a reliable potential to forecast the sensitivity to ICI immunotherapy for CESC patients, which was able to facilitate personalized counseling for immunotherapy. Further testing of this model in clinical practice will be necessary for prognostic stratification and treatment management in the future.

References

• 1

  AlvarezM.SimonettaF.BakerJ.MorrisonA. R.WenokurA. S.PieriniA.et al (2020). Indirect impact of PD-1/PD-L1 blockade on a murine model of NK cell exhaustion. Front. Immunol.11, 7. 10.3389/fimmu.2020.00007

  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  BhattacharyaS.AndorfS.GomesL.DunnP.SchaeferH.PontiusJ.et al (2014). ImmPort: disseminating data to the public for the future of immunology. Immunol. Res.58 (2–3), 234–239. 10.1007/s12026-014-8516-1

  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  BindeaG.MlecnikB.TosoliniM.KirilovskyA.WaldnerM.ObenaufA. C.et al (2013). Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity39 (4), 782–795. 10.1016/j.immuni.2013.10.003

  • CrossRef
  • Google Scholar

• 4

  BommhardtU.SchravenB.SimeoniL. (2019). Beyond TCR signaling: emerging functions of lck in cancer and immunotherapy. Int. J. Mol. Sci.20 (14), 3500. 10.3390/ijms20143500

  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  BrayF.FerlayJ.SoerjomataramI.SiegelR. L.TorreL. A.JemalA. (2018). Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J. Clin.68 (6), 394–424. 10.3322/caac.21492

  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  ByrneA.SavasP.SantS.LiR.VirassamyB.LuenS. J.et al (2020). Tissue-resident memory T cells in breast cancer control and immunotherapy responses. Nat. Rev. Clin. Oncol.17 (6), 341–348. 10.1038/s41571-020-0333-y

  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  CharoentongP.FinotelloF.AngelovaM.MayerC.EfremovaM.RiederD.et al (2017). Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep.18 (1), 248–262. 10.1016/j.celrep.2016.12.019

  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  ChenW.ZhengR.BaadeP. D.ZhangS.ZengH.BrayF.et al (2016). Cancer statistics in China, 2015. CA Cancer J. Clin.66 (2), 115–132. 10.3322/caac.21338

  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  ChungH. C.RosW.DelordJ. P.PeretsR.ItalianoA.Shapira-FrommerR.et al (2019). Efficacy and safety of pembrolizumab in previously treated advanced cervical cancer: results from the phase II KEYNOTE-158 study. J. Clin. Oncol.37 (17), 1470–1478. 10.1200/JCO.18.01265

  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  CurtiB. D. (2018). Immunotherapy in advanced renal cancer–is cure possible?N. Engl. J. Med.378 (14), 1344–1345. 10.1056/NEJMe1801682

  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  DingL.ChenF. (2019). Predicting tumor response to PD-1 blockade. N. Engl. J. Med.381 (5), 477–479. 10.1056/NEJMcibr1906340

  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  ErmanB.FırtınaS.FışgınT.BozkurtC.ÇipeF. E. (2020). Biallelic form of a known CD3E mutation in a patient with severe combined immunodeficiency. J. Clin. Immunol.40 (3), 539–542. 10.1007/s10875-020-00752-3

  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  FillonM. (2020). Immune checkpoint inhibitors may be safe for patients with preexisting autoimmune disease. CA Cancer J. Clin.70 (1), 3–4. 10.3322/caac.21587

  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  FirtinaS.NgY. Y.NgO. H.NepesovS.YesilbasO.KilercikM.et al (2017). A novel pathogenic frameshift variant of CD3E gene in two T-B+ NK+ SCID patients from Turkey. Immunogenetics69 (10), 653–659. 10.1007/s00251-017-1005-7

  • CrossRef
  • Google Scholar

• 15

  FischerA.de Saint BasileG.Le DeistF. (2005). CD3 deficiencies. Curr. Opin. Allergy Clin. Immunol.5 (6), 491–495. 10.1097/01.all.0000191886.12645.79

  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  FrenelJ. S.Le TourneauC.O’NeilB.OttP. A.Piha-PaulS. A.Gomez-RocaC.et al (2017). Safety and efficacy of pembrolizumab in advanced, programmed death ligand 1-positive cervical cancer: results from the phase Ib KEYNOTE-028 trial. J. Clin. Oncol.35 (36), 4035–4041. 10.1200/JCO.2017.74.5471

  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  FuehrerM.PannickeU.SchuetzC.JacobsenE. M.SchulzA.FriedrichW.et al (2014). Successful haploidentical hematopoietic stem cell transplantation in a patient with SCID due to CD3ε deficiency: need for IgG-substitution 6 years later. Klin Padiatr226 (3), 149–153. 10.1055/s-0033-1361142

  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  GaffneyS. G.PerryE. B.ChenP. M.GreensteinA.KaechS. M.TownsendJ. P. (2019). The landscape of novel and complementary targets for immunotherapy: an analysis of gene expression in the tumor microenvironment. Oncotarget10 (44), 4532–4545. 10.18632/oncotarget.27027

  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  GilJ.BustoE. M.GarcillánB.CheanC.García-RodríguezM. C.Díaz-AldereteA.et al (2011). A leaky mutation in CD3D differentially affects αβ and γδ T cells and leads to a Tαβ-Tγδ + B + NK + human SCID. J. Clin. Invest.121 (10), 3872–3876. 10.1172/JCI44254

  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  GulatiP.RühlJ.KannanA.PircherM.SchuberthP.NytkoK. J.et al (2018). Aberrant lck signal via CD28 costimulation augments antigen-specific functionality and tumor control by redirected T cells with PD-1 blockade in humanized mice. Clin. Cancer Res.24 (16), 3981–3993. 10.1158/1078-0432.CCR-17-1788

  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  HamanishiJ.MandaiM.KonishiI. (2017). Immune checkpoint (PD-1 and CTIA-4) signal inhibitors for gynecologic oncology; up to date. Nihon Rinsho75 (2), 234–244.

  • Google Scholar
  • View reference in article

• 22

  HanahanD.CoussensL. M. (2012). Accessories to the crime: functions of cells recruited to the tumor microenvironment. Cancer Cell21 (3), 309–322. 10.1016/j.ccr.2012.02.022

  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  HartM.Walch-RückheimB.FriedmannK. S.RheinheimerS.TänzerT.GlombitzaB.et al (2019). miR-34a: a new player in the regulation of T cell function by modulation of NF-κB signaling. Cell Death Dis.10 (2), 46. 10.1038/s41419-018-1295-1

  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  HsuJ.HodginsJ. J.MaratheM.NicolaiC. J.Bourgeois-DaigneaultM. C.TrevinoT. N.et al (2018). Contribution of NK cells to immunotherapy mediated by PD-1/PD-L1 blockade. J. Clin. Invest.128 (10), 4654–4668. 10.1172/JCI99317

  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  IrizarryR. A.HobbsB.CollinF.Beazer-BarclayY. D.AntonellisK. J.ScherfU.et al (2003). Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics4 (2), 249–264. 10.1093/biostatistics/4.2.249

  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  JemalA.BrayF.CenterM. M.FerlayJ.WardE.FormanD. (2011). Global cancer statistics. CA Cancer J. Clin.61 (2), 69–90. 10.3322/caac.20107

  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  JiangP.GuS.PanD.FuJ.SahuA.HuX.et al (2018). Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat. Med.24 (10), 1550–1558. 10.1038/s41591-018-0136-1

  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  KamarudinA. N.CoxT.Kolamunnage-DonaR. (2017). Time-dependent ROC curve analysis in medical research: current methods and applications. BMC Med. Res. Methodol.17 (1), 53. 10.1186/s12874-017-0332-6

  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  KennedyL. B.SalamaA. K. S. (2020). A review of cancer immunotherapy toxicity. CA Cancer J. Clin.70 (2), 86–104. 10.3322/caac.21596

  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  KlintmanM.BuusR.CheangM. C.SheriA.SmithI. E.DowsettM. (2016). Changes in expression of genes representing key biologic processes after neoadjuvant chemotherapy in breast cancer, and prognostic implications in residual disease. Clin. Cancer Res.22 (10), 2405–2416. 10.1158/1078-0432.CCR-15-1488

  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  LiT.FanJ.WangB.TraughN.ChenQ.LiuJ. S.et al (2017). TIMER: a web server for comprehensive analysis of tumor-infiltrating immune cells. Cancer Res.77 (21), e108–e110. 10.1158/0008-5472.CAN-17-0307

  • CrossRef
  • Google Scholar
  • View reference in article

• 32

  LiuJ.WuZ.WangY.NieS.SunR.YangJ.et al (2020). A prognostic signature based on immune-related genes for cervical squamous cell carcinoma and endocervical adenocarcinoma. Int. Immunopharmacol.88, 106884. 10.1016/j.intimp.2020.106884

  • CrossRef
  • Google Scholar
  • View reference in article

• 33

  LiuM.SunQ.WangJ.WeiF.YangL.RenX. (2019b). A new perspective: exploring future therapeutic strategies for cancer by understanding the dual role of B lymphocytes in tumor immunity. Int. J. Cancer144 (12), 2909–2917. 10.1002/ijc.31850

  • CrossRef
  • Google Scholar
  • View reference in article

• 34

  MatsuzakiJ.TsujiT.ChodonT.RyanC.KoyaR. C.OdunsiK. (2019). A rare population of tumor antigen-specific CD4 + CD8 + double-positive αβ T lymphocytes uniquely provide CD8-independent TCR genes for engineering therapeutic T cells. J. Immunother. Cancer7 (1), 7. 10.1186/s40425-018-0467-y

  • CrossRef
  • Google Scholar
  • View reference in article

• 35

  McLachlanJ.BoussiosS.OkinesA.GlaessgenD.BodlarS.KalaitzakiR.et al (2017). The impact of systemic therapy beyond first-line treatment for advanced cervical cancer. Clin. Oncol. (R Coll. Radiol.)29 (3), 153–160. 10.1016/j.clon.2016.10.002

  • CrossRef
  • Google Scholar
  • View reference in article

• 36

  MiaoY. R.ZhangQ.LeiQ.LuoM.XieG. Y.WangH.et al (2020). ImmuCellAI: a unique method for comprehensive T-cell subsets abundance prediction and its application in cancer immunotherapy. Adv. Sci. (Weinh)7 (7), 1902880. 10.1002/advs.201902880

  • CrossRef
  • Google Scholar
  • View reference in article

• 37

  MilowskyM. I.RumbleR. B.BoothC. M.GilliganT.EapenL. J.HaukeR. J.et al (2016). Guideline on muscle-invasive and metastatic bladder cancer (European association of urology guideline): American society of clinical oncology clinical practice guideline endorsement. J. Clin. Oncol.34 (16), 1945–1952. 10.1200/JCO.2015.65.9797

  • CrossRef
  • Google Scholar
  • View reference in article

• 38

  OtterS. J.ChatterjeeJ.StewartA. J.MichaelA. (2019). The role of biomarkers for the prediction of response to checkpoint immunotherapy and the rationale for the use of checkpoint immunotherapy in cervical cancer. Clin. Oncol. (R Coll. Radiol.)31 (12), 834–843. 10.1016/j.clon.2019.07.003

  • CrossRef
  • Google Scholar

• 39

  PatelV. G.OhW. K.GalskyM. D. (2020). Treatment of muscle-invasive and advanced bladder cancer in 2020. CA Cancer J. Clin.70 (5), 404–423. 10.3322/caac.21631

  • CrossRef
  • Google Scholar

• 40

  PfaendlerK. S.TewariK. S. (2016). Changing paradigms in the systemic treatment of advanced cervical cancer. Am. J. Obstet. Gynecol.214 (1), 22–30. 10.1016/j.ajog.2015.07.022

  • CrossRef
  • Google Scholar

• 41

  ShiM. J.MengX. Y.WuQ. J.ZhouX. H. (2019). High CD3D/CD4 ratio predicts better survival in muscle-invasive bladder cancer. Cancer Manag. Res.11, 2987–2995. 10.2147/CMAR.S191105

  • CrossRef
  • Google Scholar

• 42

  SuryadevaraC. M.DesaiR.FarberS. H.ChoiB. D.SwartzA. M.ShenS. H.et al (2019). Preventing lck activation in CAR T cells confers treg resistance but requires 4-1BB signaling for them to persist and treat solid tumors in nonlymphodepleted hosts. Clin. Cancer Res.25 (1), 358–368. 10.1158/1078-0432.CCR-18-1211

  • CrossRef
  • Google Scholar

• 43

  TorreL. A.BrayF.SiegelR. L.FerlayJ.Lortet-TieulentJ.JemalA. (2015). Global cancer statistics, 2012. CA Cancer J. Clin.65 (2), 87–108. 10.3322/caac.21262

  • CrossRef
  • Google Scholar

• 44

  WangH.LengerichB. J.AragamB.XingE. P. (2019a). Precision Lasso: accounting for correlations and linear dependencies in high-dimensional genomic data. Bioinformatics35 (7), 1181–1187. 10.1093/bioinformatics/bty750

  • CrossRef
  • Google Scholar

• 45

  WangQ.LiP.WuW. (2019b). A systematic analysis of immune genes and overall survival in cancer patients. BMC Cancer19 (1), 1225. 10.1186/s12885-019-6414-6

  • CrossRef
  • Google Scholar

• 46

  WeiQ.BrzostekJ.SankaranS.CasasJ.HewL. S.YapJ.et al (2020). Lck bound to coreceptor is less active than free Lck. Proc. Natl. Acad. Sci. U S A.117 (27), 15809–15817. 10.1073/pnas.1913334117

  • CrossRef
  • Google Scholar

• 47

  YangS.WuY.DengY.ZhouL.YangP.ZhengY.et al (2019). Identification of a prognostic immune signature for cervical cancer to predict survival and response to immune checkpoint inhibitors. Oncoimmunology8 (12), e1659094. 10.1080/2162402X.2019.1659094

  • CrossRef
  • Google Scholar

• 48

  YangY.ZangY.ZhengC.LiZ.GuX.ZhouM.et al (2020). CD3D is associated with immune checkpoints and predicts favorable clinical outcome in colon cancer. Immunotherapy12 (1), 25–35. 10.2217/imt-2019-0145

  • CrossRef
  • Google Scholar

• 49

  YoshiharaK.ShahmoradgoliM.Mart¨nezE.VegesnaR.KimH.Torres-GarciaW.et al (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun.4, 2612. 10.1038/ncomms3612

  • CrossRef
  • Google Scholar

• 50

  YuG.WangL. G.HanY.HeQ. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. Omics16 (5), 284–287. 10.1089/omi.2011.0118

  • CrossRef
  • Google Scholar

• 51

  ZaretskyJ. M.Garcia-DiazA.ShinD. S.Escuin-OrdinasH.HugoW.Hu-LieskovanS.et al (2016). Mutations associated with acquired resistance to PD-1 blockade in melanoma. N. Engl. J. Med.375 (9), 819–829. 10.1056/NEJMoa1604958

  • CrossRef
  • Google Scholar

• 52

  ZhangB.HorvathS. (2005). A general framework for weighted gene co-expression network analysis. Stat. Appl. Genet. Mol. Biol.4, Article17. 10.2202/1544-6115.1128

  • CrossRef
  • Google Scholar

• 53

  ZhaoW.GuY. H.SongR.QuB. Q.XuQ. (2008). Sorafenib inhibits activation of human peripheral blood T cells by targeting LCK phosphorylation. Leukemia22 (6), 1226–1233. 10.1038/leu.2008.58

  • CrossRef
  • Google Scholar