Development and Validation of an Immune-Related Prognostic Signature in Cervical Cancer

Background: Cervical cancer is the fourth most frequent gynecological malignancy across the world. Immunotherapies have proved to improve prognosis of cervical cancer. However, few studies on immune-related prognostic signature had been reported in cervical cancer.

Methods: Raw data and clinical information of cervical cancer samples were downloaded from TCGA and UCSC Xena website. Immunophenoscore of immune infiltration cells in cervical cancer samples was calculated through the ssGSEA method using GSVA package. WGCNA, Cox regression analysis, LASSO analysis, and GSEA analysis were performed to classify cervical cancer prognosis and explore the biological signaling pathway.

Results: There were eight immune infiltration cells associated with prognosis of cervical cancer. Through WGCNA, 153 genes from 402 immune-related genes were significantly correlated with prognosis of cervical cancer. A 15-gene signature demonstrated powerful predictive ability in prognosis of cervical cancer. GSEA analysis showed multiple signaling pathways containing Programmed cell death ligand-1 (PD-L1) expression and PD-1 checkpoint pathway differences between high-risk and low-risk groups. Furthermore, the 15-gene signature was associated with multiple immune cells and immune infiltration in tumor microenvironment.

Conclusion: The 15-gene signature is an effective potential prognostic classifier in the immunotherapies and surveillance of cervical cancer.

Introduction

According to estimates from GLOBOCAN 2018, cervical cancer was the fourth most common cancer among women, with approximately 570,000 new cases and 311,000 deaths (1). Cure rate for earliest stages is more than 90%, whereas locally advanced lesions treated with multimodality therapy can only achieve 65% cure rate in stage II lesions and 55% in stage III lesions (2). Therapies to improve survivorship are in desperate need for cervical cancer especially locally advanced disease.

With rapid development of precision medicines, novel therapeutic strategies, especially immunotherapies, have been proposed to significantly improve clinical outcomes of cervical cancer (3, 4). Harnessing an antitumor immune response has been a fundamental strategy in cancer immunotherapy. A paradigm shift has appeared in cancer immunotherapy: from traditional immune enhancement with low-objective responses and frequent adverse events to more effective and less toxic reactions immune normalization (5, 6).

Cancer immunotherapies, such as administration of the cytokine Interleukin-2 (IL-2), adoptive cell transfer, and the checkpoint modulators CTLA-4 and PD-1, have proved effective in clinical practice (6). Blockade of the checkpoint modulators cytotoxic T-lymphocyte-associated antigen 4 (CTLA-4) and programmed cell death protein-1 (PD-1) starts the field of immune normalization in immunotherapy. Upregulated PD-1 in tumor microenvironment inhibits an effector T cell antitumor immune response, and therapies blocking this pathway have proven effective against multiple tumor types (5). Anti-PD therapy perform antitumor immunity mainly through the following three principles:1) targeting a tumor-induced immune escape mechanism, 2) selectively modulating immunity in the tumor microenvironment, and 3) resetting immunity in the tumor microenvironment (5). It has been proved that patients with high mutation burden and burden of potential neoepitopes benefit more from immunological checkpoint blockade (7–11). Immunotherapy plays a dispensable role in management of cervical cancer. In KEYNOTE-158, pembrolizumab has been approved by the US Food and Drug Administration for use in advanced cervical cancer with progressive disease either during or after chemotherapy (12). The objective response rate was 26.3% with a disease control rate of 68.4% for immunotherapy in cervical cancer (4). However, only small proportion patients benefited from immunotherapy, and this proportion will hopefully increase with better patient selection and combination therapy. Hence, it is urgent to find potential biomarkers for prediction of response to checkpoint immunotherapy and the rationale for the use of checkpoint immunotherapy.

In the present study, we qualify immune cells infiltration in cervical cancer and analyze the correlation between immune cells and cancer prognosis. Hub genes regulating prognosis through immune infiltration in cervical cancer were identified by weighted gene co-expression network analysis (WGCNA) and least absolute shrinkage and selection operator (LASSO). It was suggested by infiltrated immune cells and pathway enrichment analysis that our immune-related signature was closely related to tumor prognosis and could predict response of immunotherapy. A robust immune-related prognostic signature based on transcriptomics in cervical cancer was constructed and validated.

Materials and Methods

Data Source and Processing

The gene expression profiles and clinical information of cervical cancer were downloaded from The Cancer Genome Atlas (TCGA) Genomic Data Commons Data Portal (https://portal.gdc.cancer.gov/). Patients with pathologically confirmed cervical cancer and complete information about transcriptomics overall survivals (OSs) were included in this study. Finally, a total of 304 primary cervical cancer samples, two metastatic cervical cancer samples, and three normal cervix samples from TCGA were analyzed in our study. Details on clinical information of the included samples were summarized in Table 1. A workflow of this study was indicated in Figure 1.

TABLE 1

Table 1 Clinical characteristics of 309 samples from TCGA.

FIGURE 1

Figure 1 Flowchart of the 15-gene signature prognostic model construction process.

Infiltration of Immune Cells

Tumor-infiltrating immune cells can be quantified from RNA sequencing (RNA-seq) data of human tumors using bioinformatics approaches. Single-sample gene set enrichment analysis (ssGSEA) calculated and qualified the infiltration of immune cells through RNA-seq data. The infiltration of 28 immune cells was obtained. Immunophenoscore (normalized enrichment score) of each TCGA cervical cancer sample was calculated through ssGSEA. The ssGSEA ranks the genes by their absolute expression in a sample and computes enrichment score by integrating the differences between the empirical cumulative distribution functions of the gene ranks (13).

Survival Analysis

Univariate Cox regression analysis was performed to identify the association between cervical cancer survival and immunophenoscore of infiltrated immune cells. Forest plot was drawn to demonstrate the influence of infiltrated immune cells on survival. The best separation statistic was performed by use of “survminer” package that divides gene expression into high groups and low groups based on best separation. Next, Kaplan–Meier curve was made to further analyze relationship between survival and infiltration immune cells.

Construction of Weighted Gene Co-Expression Network Analysis

WGCNA is an algorithm for finding genetic interactions in a weighted manner. It is used to build a gene co-expression networks to mine network modules closely associated with clinical traits through systematic biological method (14). In this study, immunophenoscore of the infiltrated immune cells was regarded as target clinical traits. Genes expressing not available (NA) were removed. The top 25% genes with most median absolute deviation used as a robust measure of variability were selected for WGCNA analysis (15, 16).

Using the WGCNA function adjacency, an adjacency matrix is constructed by computing the Pearson correlation between all pairs of genes in the selected sample. Genes were divided into different gene modules based on the dissimilarity measure. A hierarchical clustering tree was constructed with different branches of the tree representing different gene modules (15–17).

The WGCNA network was built. There were eight immune cells associated with cervical cancer survival. The magenta module was the module that most related to infiltration of immune cells. A total of 402 genes contained in magenta module were analyzed for survival using the standardized expression data FPKM. Finally, we found that 153 genes were associated with survival (p < 0.05) in magenta module.

Screening Hub Genes by LASSO and Survival Analysis

The 153 genes with the highest correlation associated with survival were picked out from WGCNA analysis. Then, hub genes were further screened from the 153 genes by the use of LASSO. Survival analysis of hub genes was performed using “survival” and “survminer” packages to verify whether differently expressed genes affected tumor prognosis. Data of cervical cancer from University of California, Santa Cruz (UCSC) Xena were downloaded and applied to compare the expression of hub genes in normal cervix and cervical cancer.

Construction of Prognostic Scoring Model

Patients with cervical cancer in TCGA were randomly divided into two groups by stage: the training group of 70% and the test group of 30% by the use of caret R package. Finally, a total of 199 patients were used as training group and 74 patients were regarded as test group. The selected key genes for support vector machine analysis were used to fit a LASSO Cox–proportional hazards (Cox-PH) model for selecting an optimal panel of predictive genes with penalized package of R (18). Optimal lambda value was computed through a 10 cross-validations. Next, Cox-PH coefficients and infiltration of immune cells levels of these selected genes were used to calculate prognostic score as follows: Risk Score = Σ (coef_gene × immunophenoscore_gene), where Coef_gene and immunophenoscore_gene suggest Cox-PH coefficient and immunophenoscore level of a gene, respectively.

Validation of Prognostic Scoring Model

We calculated risk score of every patient in the training and test group by the model. We separated the training group and test group into a high-risk group and a low-risk group with the median risk score as cutoff, respectively. Kaplan–Meier curve was applied to obtain and compare the OS time of two risk groups. ROC (operating characteristic curve) was performed to evaluate the predictive accuracy of the model.

Comparison of Eight Immune Cell Subtypes Between High-Risk and Low-Risk Groups

To explore the differences of immune cell subtypes between high-risk and low-risk groups, the eight immune cell subtypes associated with OS in cervical cancer were assessed in test and train cohorts. Mann–Whitney U-test was used to compare differences in immune cell subtypes in the high-risk and low-risk groups.

Gene Set Enrichment Analysis (GSEA)

All 304 patients with primary cervical cancer were evaluated by prognostic scoring model and then divided into high-risk group and low-risk group based on the standard cutoff. The global gene expression was analyzed and displayed with volcanic maps using the limma package in R. GSEA was conducted respectively to search “all gene sets” enriched in the samples with high-risk group and low-risk group. The differentially expressed genes were enriched in PD-L1 pathway functional sets (19).

Results

Quantify Immune Cell Infiltration and Survival Analysis

We used ssGSEA to quantify mRNA data for immune cell infiltration. Finally, 28 infiltrating immune cells were included. Immunophenoscores of 304 primary cervical cancer samples, two metastatic cervical cancer samples, and three normal cervix samples for 28 immune cells were calculated and demonstrated in Figure 2A. Univariate Cox regression analysis was performed to identify association between immune cells and OS of cervical cancer, and the results are shown in Figure 2B. It indicated that eight immune cells were correlated with OS significantly. Moreover, these eight immune cells (activated B cell, activated CD8 T cell, eosinophil, monocyte, activated CD4 T cell, effector memory CD8 T cell, immature B cell, and plasmacytoid dendritic cell) were all protective factors for OS. Kaplan–Meier curve validated that high expression of these eight immune cells were related to longer OS time suggested worse prognosis (Figure 3).

FIGURE 2

Figure 2 Immunophenoscores and hazard ratios of 28 immune cells in cervical cancer. (A) Boxplot of immunophenoscores from 28 immune cells TCGA cervical cancer and normal samples. (B) Forest plot of hazard ratios in each infiltrated immune cell for OS. Infiltrated immune cell subtypes significantly associated with overall survival were underlined with a red dashed line. *: P<0.05, **: P<0.01, ns, no significance.

FIGURE 3

Figure 3 Kaplan–Meier curves for overall survival of eight infiltrated immune cells. Activated B cell, activated CD8 T cell, eosinophil, monocyte, activated CD4 T cell, effector memory CD8 T cell, immature B cell, and plasmacytoid dendritic cell were associated with overall survival of cervical cancer.

The Weighted Gene Co-Expression Network Analysis Construction and Key Module Identification

The top 25% of the gene expression of variance was screened by the quartile of gene expression level, and 609 genes were screened out to construct co-expressed gene networks and the sample dendrogram and trait heatmap are constructed (Figure 4A). In this study, the power of β = 8 was selected to ensure a scale-free network (Figures 4B, C) (scale-free R² = 0.9, slope = −1.65).

FIGURE 4

Figure 4 Sample tree and trait heatmap of 28 immune cells. (A) Sample dendrogram and trait heatmap of 28 immune cells. The clinical trait information is immunophenoscores of 28 infiltrated immune cells in each sample. (B) Analysis of the scale-free fit index for various soft-thresholding powers (β) and the mean connectivity for soft threshold powers. The soft-thresholding power in the WGCNA was determined on the basis of a scale-free R² (R² = 0.90). The left panel presents the relationship between the soft-threshold and scale-free R². The right panel presents the relationship between the soft-threshold and mean connectivity. (C) Histogram of connectivity distribution when β = 8 and checking the scale-free topology when β = 8.

After merging the modules with the high similarity of feature genes in the gene cluster dendrogram through a cutline (0.25) (Figure 5A), nine modules were identified by the dynamic tree cut method. The clustering dendrograms of genes were shown in Figure 5B. A heatmap illustrating the correlation between immunophenoscores of infiltrated immune cells and key genes in the module was created (Figure 5C).

FIGURE 5

Figure 5 Weighted gene co-expression network in cervical cancer. (A) Clustering of module eigengenes. A cutline (0.25) was selected for module dendrogram and merged some modules according to dissimilarity of estimated module eigengenes. (B) A dendrogram of the differentially expressed genes clustered based on different metrics. Each branch in the figure represents one gene, and every color below represents one co-expression module. (C) A heatmap showing the correlation between the gene module and clinical traits. The correlation coefficient in each cell represented the correlation between gene module and the clinical traits. The magenta module showed the highest correlated with the infiltrated immune cells.

It was obviously that magenta module was the most correlated with immune infiltration in the heatmap. Hence, magenta module was selected as the clinical significant module for further analysis. We found that the correlation coefficient between magenta module and 11 immune infiltration cells was ≥ 0.7, which suggested strong correlation. This indicates that the genes in this module are most relevant to tumor OS. All 402 genes in magenta module were analyzed for survival of cervical cancer. Through univariate Cox regression survival analysis, 153 genes found to be significantly correlated with prognosis of cervical cancer were selected for further analysis.

Screening the Key Genes and Survival Analysis

To develop a prognostic scoring model based on WGCNA, the 153 key genes were used to fit the LASSO Cox-PH model. Using parameter lambda (0.038) obtained upon performing 1000 cross-validations, a combination of 15 genes was obtained. These 15 genes were as follows: LAG3, CD74, CCL22, CH25H, OLR1, MIAT, BATF, IKZF3, TRARG, ACSL6, C11orf21, GTSF1, APOL1, CD1C, and LINC00158 (Table 2 and Figure 6). Risk score of each sample in the training set was calculated. Then, the samples in the training cohort were divided into high- and low-risk groups according to median risk score. As demonstrated in Figure 6, the high-risk group patients in the training cohort have shorter OS time than the low-risk group (p < 0.001, HR: 5.1, 95% CI: 3.4–7.5), with an AUC of 0.803 in 1 year, 0.809 in 3 years, and 0.800 in 5 years. Next, the 15-gene signature was then verified by the validation cohort from TCGA. Consistently, the high-risk group patients in the validation cohort had worse prognosis than the low-risk group (p < 0.001, HR: 1.3, 95% CI: 1.1–1.4), and the AUC was 0.845 in 1 year, 0.703 in 3 years, and 0.761 in 5 years (Figure 7). All these suggested that the 15-gene signature associated with immune infiltration is able to predict prognosis in patients with cervical cancer.

TABLE 2

Table 2 Fifteen genes from LASSO analysis.

FIGURE 6

Figure 6 Construction of the prognostic classifier. (A, B) Determination of the number of factors by the LASSO analysis. (C) The distribution of risk score in test cohort. (D) The survival duration and status of patients in test cohort. (E) A heatmap of immune related genes in the classifier in test cohort. (F) Kaplan–Meier curve for patients with cervical cancer in test cohort. (G) ROC curve for cervical cancer in test cohort.

FIGURE 7

Figure 7 The distribution of time-dependent ROC curves and Kaplan–Meier survival based on the integrated classifier in the train cohort. (A) The distribution of risk score in train cohort. (B) The survival duration and status of patients in train cohort. (C) A heatmap of immune related genes in the classifier in train cohort. (D) Kaplan–Meier curve for patients with cervical cancer in train cohort. (E) ROC curve for cervical cancer in train cohort.

Immune Cell Subtypes Between High-Risk and Low-Risk Groups

The expression levels of the 15 genes in test and validation cohorts are shown in Figure 8A. Different immune scores had differential OS in patients with cervical cancer. In both test cohorts and train cohorts, these eight immune cell subtypes (activated B cell, activated CD8 T cell, eosinophil, monocyte, activated CD4 T cell, effector memory CD8 T cell, immature B cell, and plasmacytoid dendritic cell) expressed differentially in high-risk and low-risk groups. (Figures 8B, C).

FIGURE 8

Figure 8 Expression profile of 15 genes and immunophenoscores of eight immune cell subtypes in test and train sets. (A) Expression profile of 15 genes. The test cohort is shown in the left, and the train cohort is shown in the right. (B, C) Immunophenoscores of eight infiltrated immune cells in high-risk and low-risk groups based on the integrated classifier. The test cohort is shown in the left, and the train cohort is shown in the right. *: P<0.05, **: P<0.01, ***: P<0.001****: P<0.0001, ns, no significance.

GSEA Analysis

All 304 primary cervical cancer samples were divided into high-risk group and low-risk group based on the 15-gene signature. To identify potential function of the 15 key genes, GSEA was conducted respectively to search “all gene sets” enriched in all 304 samples. GSEA showed 33 significant Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway associated with risk score, including PD-L1 expression and PD-1 checkpoint pathway in cancer (Figure 9 and Table 3). These results validated that the key genes in clinical significant module were mainly involved in the regulation of immune system. From GSEA, the 15-gene signature obviously participated in regulation of PD-1/PD-L1 pathway in cervical cancer.

FIGURE 9

Figure 9 GSEA analysis and PD-L1 pathway. (A) GSEA analysis demonstrated 33 KEGG pathways associated with risk score in cervical cancer. (B) PD-L1 expression and PD-1 checkpoint pathway in cancer was one of the 33 pathways that significantly related with risk score.

TABLE 3

Table 3 Significant KEGG pathways from GSEA analysis.

Discussion

Cervical cancer represents a major public health problem even in developed countries such as United States, with 14,000 new cases and 4,285 deaths estimated in 2022 (20). Increasing incidence of cervical adenocarcinoma and attenuation of earlier declines for cervical squamous cell carcinoma emphasize the importance of the need for improved therapeutic options to reduce the burden of cervical cancer (21). The most important risk factors affecting prognosis of cervical cancer are stage, status of the lymph nodes, tumor volume, depth of tumor invasion into the cervical stroma, and Lymphovascular space invasion (LVSI) (22). The 5-year OS decreased significantly with rising stages (23). A substantial percentage of patients with advanced cervical cancer will undergo recurrence and poor prognosis (2, 24), and the recurrence rate fluctuates from 10% to 74% (25). Tumor stage was one of the most pivotal factors related to recurrence: 10% for stage IB, 17% for stage IIA, 23% for stage IIB, 42% for stage III, and 74% for stage IVA (26). Prognosis of metastatic and recurrent cervical cancer was extremely poor with a median survival time of 12 months (27). Therefore, improved therapeutic options in management of cervical cancer, especially locally advanced cervical cancer is in urgent need.

High-risk subtypes of the human papilloma virus (HPV) are the cause of cervical cancer (28, 29). Viral oncoproteins E6 and E7 leads to dysregulation of p53 and Hypoxia-inducible factor-1 (HIF-1) alpha, thus affecting cell cycle proteins and VEGF expression (28–32). HPV is highly immunogenic and elicits immune responses in humans; thus, immune response might play important roles in carcinogenesis of cervical cancer. Incidence of cervical cancer was substantial declined in countries with high HPV vaccine coverage (33–35). Immunotherapy fights against tumor cells through activating endogenous immune response, which seems to switch on the new frontier of the anticancer treatment (36). Immunotherapy includes different approaches, such as active immunotherapy (vaccine), passive immunotherapy (adoptive cellular transfer, antibodies, and cytokines), and immunomodulation (cyclooxygenase 2 inhibitor).

Discovery of immune checkpoint such as CTLA-4 and PD-1 plays an indispensable role in the development of cancer immunotherapy. It was surprising that immune checkpoint inhibitors anti–CTLA-4 and anti–PD-1 displayed enormous success in solid tumors (37). Similar to other solid tumors, novel immunotherapeutic approaches, such as immune checkpoint inhibitors, have shown encouraging results in cervical cancer (12). Implementing immunotherapeutic approaches earlier in advanced cervical cancer would seem to be most appropriate. However, the objective response rate with anti–PD-1/PD-L1 monotherapy is hovering at 20%. Moreover, immune-related toxicities and severe adverse effects can occur during PD-1/PD-L1 blockade therapy (38). PD-1/PD-L1 blockade therapy does not demonstrate efficacy in almost 80% of patients with cervical cancer, suggesting that the potential mechanism of PD-1/PD-L1 in immunotherapy remains to be further clarified. Thus, new immune checkpoint inhibitors or comprehensive understanding of specific mechanism underlying PD-1/PD-L1 regulation in carcinogenesis is in urgent need.

Immune infiltration of cervical cancer determines the immune activation of tumor microenvironment and is related with clinical outcome of patients. In this study, immunophenoscore of 28 infiltrated immune cells in TCGA cervical cancer was calculated. Univariate Cox regression analysis showed that activated B cell, activated CD8 T cell, eosinophil, monocyte, activated CD4 T cell, effector memory CD8 T cell, immature B cell, and plasmacytoid dendritic cell were strongly associated with OS of cervical cancer. WGCNA revealed that magenta module was the most relevant module to immune infiltration. In total, we identified 10 kinds of infiltrated immune cells with strong correlation with magenta module. These were as follows: activated B cell (R = 0.88), activated CD8 T cell (R = 0.82), activated dendritic cell (R = 0.79), central memory CD4 T cell (R = 0.78), effector memory CD8 T cell (R = 0.85), immature B cell (R = 0.93), macrophage (R = 0.72), Myeloid-derived suppressor cells (MDSC) (R = 0.85), regulatory T cell (R = 0.71), T follicular helper cell (R = 0.72), and type 1 T helper cell (R = 0.83).

By the use of WGCNA, magenta module showed the highest correlation with the immune infiltration of cervical cancer. The magenta module contained 609 immune-related genes, and 153 genes were picked up for further analysis after taking intersection of prognosis and immune infiltration. Through LASSO analysis, 15 immune-related genes (LAG3, CD74, CCL22, CH25H, OLR1, MIAT, BATF, IKZF3, TRARG, ACSL6, C11orf21, GTSF1, APOL1, CD1C, and LINC00158) were included in prognostic classifier. In addition to CTLA-4– and PD-1/PD-L1–targeted cancer immunotherapy, LAG3 (lymphocyte activation gene-3, CD223) is the third clinically targeted inhibitory receptor (39). CD74 (invariant chain) plays a dispensable role in process of immune systems that it participates in antigen presentation, B-cell differentiation, and inflammatory signaling. CD74 has the potential to be a therapeutic target in cancer and autoimmune disease (40). The chemokine CCL22 promoted regulatory T cell communication with dendritic cells to control immunity and was associated with poor prognosis (41, 42). CH25H produces 25-hydroxycholesterol, which inhibited tumor-derived extracellular vesicles uptake and correlated with prognosis in patients with melanoma (43). OLR1 (oxidized low-density lipoprotein (LDL) receptor 1) was a possible link between obesity, dyslipidemia, and cancer. OLR1 played carcinogenic role by activating Nuclear factor kappa B (NF-κB) pathway to promote proliferation and migration and to inhibit apoptosis and de novo lipogenesis (44). MIAT Long non-coding RNAs (lncRNA) was overexpressed in a number of malignancies and caused poor prognosis (45–47). Basic leucine zipper transcription factor, ATF-like (BATF) was an important transcription factor regulating differentiation of early effector CD8⁺ T cells (48) and was a prognostic indicator for patients with colon cancer (49). IKZF3 promoted growth of multiple tumors to cause poor prognosis (50).

Immune infiltration played important roles in the survival of cervical cancer. In this study, we identified 15 immune infiltration associated genes in cervical cancer and built a prognostic signature. Immune scores depended on the expression of these 15 genes and were associated with the survival of cervical cancer. High immune scores meant good prognosis. GSEA analysis showed that the 15-gene prognostic signature was obviously associated with PD-L1 expression and PD-1 checkpoint pathway in cancer. The prognostic signature could provide basis for potential immunotherapy in the future. Similarly, other researchers also constructed prognostic signatures for cervical cancer based on immune-related genes (51, 52). Compared with our study, although these studies use different types and different quantity of genes, all the signatures can well predict the prognosis of cervical cancer. However, the study has several limitations. First, no in vitro or in vivo molecular experiment was performed to verify our analysis. Second, our study was a retrospective study. Thus, prospective study is in need to validate the findings of our study in the future.

Conclusion

In conclusion, we successfully constructed a 15-gene prognostic signature with powerful predictive function. Differences in the OS of high- and low-risk groups are implicated in immune infiltration, tumor microenvironment, PD-L1 expression, and PD-1 checkpoint pathway. These findings revealed the underlying mechanism of immunotherapy and provided basis for cervical cancer pathogenesis and clinical treatment.

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material. Further inquiries can be directed to the corresponding authors.

Author Contributions

RS is responsible for experimental design. CWJ and HB are responsible for instrument operation. JX is responsible for data analysis. LZ and CJJ are for providing overall ideas. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the National Natural Science Foundation of China (grant no. 82002736 to CJJ) and the clinical characteristic medical technology cultivation plan of Shanghai General Hospital (no. 02.DY12.06.22.07 to LZ and no. 02.DY12.06.22.12 to JX).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank Dr. Jianming Zeng (University of Macau) and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes, especially Xiaojie Sun.

References

1. Arbyn M, Weiderpass E, Bruni L, de Sanjosé S, Saraiya M, Ferlay J, et al. Estimates of Incidence and Mortality of Cervical Cancer in 2018: A Worldwide Analysis. Lancet Global Health (2020) 8(2):e191–203. doi: 10.1016/S2214-109X(19)30482-6

CrossRef Full Text | Google Scholar

2. Monk BJ, Tewari KS, Koh WJ. Multimodality Therapy for Locally Advanced Cervical Carcinoma: State of the Art and Future Directions. J Clin Oncol Off J Am Soc Clin Oncol (2007) 25(20):2952–65. doi: 10.1200/JCO.2007.10.8324

CrossRef Full Text | Google Scholar

3. Chung HC, Ros W, Delord JP, Perets R, Italiano A, Shapira-Frommer R, et al. Efficacy and Safety of Pembrolizumab in Previously Treated Advanced Cervical Cancer: Results From the Phase II KEYNOTE-158 Study. J Clin Oncol Off J Am Soc Clin Oncol (2019) 37(17):1470–8. doi: 10.1200/JCO.18.01265

CrossRef Full Text | Google Scholar

4. Cohen PA, Jhingran A, Oaknin A, Denny L. Cervical Cancer. Lancet (2019) 393(10167):169–82. doi: 10.1016/S0140-6736(18)32470-X

CrossRef Full Text | Google Scholar

5. Sanmamed MF, Chen L. A Paradigm Shift in Cancer Immunotherapy: From Enhancement to Normalization. Cell (2018) 175(2):313–26. doi: 10.1016/j.cell.2018.09.035

CrossRef Full Text | Google Scholar

6. Tran E, Robbins PF, Rosenberg SA. 'Final Common Pathway' of Human Cancer Immunotherapy: Targeting Random Somatic Mutations. Nat Immunol (2017) 18(3):255–62. doi: 10.1038/ni.3682

CrossRef Full Text | Google Scholar

7. Asaoka Y, Ijichi H, Koike K. PD-1 Blockade in Tumors With Mismatch-Repair Deficiency. N Engl J Med (2015) 373(20):1979. doi: 10.1056/NEJMc1510353

CrossRef Full Text | Google Scholar

8. Hugo W, Zaretsky JM, Sun L, Song C, Moreno BH, Hu-Lieskovan S, et al. Genomic and Transcriptomic Features of Response to Anti-PD-1 Therapy in Metastatic Melanoma. Cell (2016) 165(1):35–44. doi: 10.1016/j.cell.2016.02.065

CrossRef Full Text | Google Scholar

9. Powles T, Eder JP, Fine GD, Braiteh FS, Loriot Y, Cruz C, et al. MPDL3280A (Anti-PD-L1) Treatment Leads to Clinical Activity in Metastatic Bladder Cancer. Nature (2014) 515(7528):558–62. doi: 10.1038/nature13904

CrossRef Full Text | Google Scholar

10. Rizvi NA, Hellmann MD, Snyder A, Kvistborg P, Makarov V, Havel JJ, et al. Cancer Immunology. Mutational Landscape Determines Sensitivity to PD-1 Blockade in Non-Small Cell Lung Cancer. Science (2015) 348(6230):124–8. doi: 10.1126/science.aaa1348

CrossRef Full Text | Google Scholar

11. van Rooij N, van Buuren MM, Philips D, Velds A, Toebes M, Heemskerk B, et al. Tumor Exome Analysis Reveals Neoantigen-Specific T-Cell Reactivity in an Ipilimumab-Responsive Melanoma. J Clin Oncol Off J Am Soc Clin Oncol (2013) 31(32):e439–442. doi: 10.1200/JCO.2012.47.7521

CrossRef Full Text | Google Scholar

12. Marabelle A, Fakih M, Lopez J, Shah M, Shapira-Frommer R, Nakagawa K, et al. Association of Tumour Mutational Burden With Outcomes in Patients With Advanced Solid Tumours Treated With Pembrolizumab: Prospective Biomarker Analysis of the Multicohort, Open-Label, Phase 2 KEYNOTE-158 Study. Lancet Oncol (2020) 21(10):1353–65. doi: 10.1016/S1470-2045(20)30445-9

CrossRef Full Text | Google Scholar

13. Finotello F, Trajanoski Z. Quantifying Tumor-Infiltrating Immune Cells From Transcriptomics Data. Cancer Immunology Immunotherapy CII (2018) 67(7):1031–40. doi: 10.1007/s00262-018-2150-z

CrossRef Full Text | Google Scholar

14. Langfelder P, Horvath S. WGCNA: An R Package for Weighted Correlation Network Analysis. BMC Bioinf (2008) 9:559. doi: 10.1186/1471-2105-9-559

CrossRef Full Text | Google Scholar

15. Dong J, Horvath S. Understanding Network Concepts in Modules. BMC Syst Biol (2007) 1:24. doi: 10.1186/1752-0509-1-24

CrossRef Full Text | Google Scholar

16. Horvath S, Dong J. Geometric Interpretation of Gene Coexpression Network Analysis. PloS Comput Biol (2008) 4(8):e1000117. doi: 10.1371/journal.pcbi.1000117

CrossRef Full Text | Google Scholar

17. Langfelder P, Zhang B, Horvath S. Defining Clusters From a Hierarchical Cluster Tree: The Dynamic Tree Cut Package for R. Bioinf (Oxford England) (2008) 24(5):719–20. doi: 10.1093/bioinformatics/btm563

CrossRef Full Text | Google Scholar

18. Ferlay J, Soerjomataram I, Dikshit R, Eser S, Mathers C, Rebelo M, et al. Cancer Incidence and Mortality Worldwide: Sources, Methods and Major Patterns in GLOBOCAN 2012. Int J Cancer (2015) 136(5):E359–386. doi: 10.1002/ijc.29210

CrossRef Full Text | Google Scholar

19. Subramanian A, Kuehn H, Gould J, Tamayo P, Mesirov JP. GSEA-P: A Desktop Application for Gene Set Enrichment Analysis. Bioinformatics (2007) 23(23):3251–3. doi: 10.1093/bioinformatics/btm369

CrossRef Full Text | Google Scholar

20. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer Statistics, 2022. CA Cancer J Clin (2022) 72(1):7–33. doi: 10.3322/caac.21708

CrossRef Full Text | Google Scholar

21. Islami F, Fedewa SA, Jemal A. Trends in Cervical Cancer Incidence Rates by Age, Race/Ethnicity, Histological Subtype, and Stage at Diagnosis in the United States. Prev Med (2019) 123:316–23. doi: 10.1016/j.ypmed.2019.04.010

CrossRef Full Text | Google Scholar

22. Peters WA 3rd, Liu PY, Barrett RJ 2nd, Stock RJ, Monk BJ, Berek JS, et al. Concurrent Chemotherapy and Pelvic Radiation Therapy Compared With Pelvic Radiation Therapy Alone as Adjuvant Therapy After Radical Surgery in High-Risk Early-Stage Cancer of the Cervix. J Clin Oncol Off J Am Soc Clin Oncol (2000) 18(8):1606–13. doi: 10.1200/JCO.2000.18.8.1606

CrossRef Full Text | Google Scholar

23. Sturdza A, Pötter R, Fokdal LU, Haie-Meder C, Tan LT, Mazeron R, et al. Image Guided Brachytherapy in Locally Advanced Cervical Cancer: Improved Pelvic Control and Survival in RetroEMBRACE, a Multicenter Cohort Study. Radiotherapy Oncol J Eur Soc Ther Radiol Oncol (2016) 120(3):428–33. doi: 10.1016/j.radonc.2016.03.011

CrossRef Full Text | Google Scholar

24. Liontos M, Kyriazoglou A, Dimitriadis I, Dimopoulos MA, Bamias A. Systemic Therapy in Cervical Cancer: 30 Years in Review. Crit Rev Oncology/Hematology (2019) 137:9–17. doi: 10.1016/j.critrevonc.2019.02.009

CrossRef Full Text | Google Scholar

25. Peiretti M, Zapardiel I, Zanagnolo V, Landoni F, Morrow CP, Maggioni A. Management of Recurrent Cervical Cancer: A Review of the Literature. Surg Oncol (2012) 21(2):e59–66. doi: 10.1016/j.suronc.2011.12.008

CrossRef Full Text | Google Scholar

26. Perez CA, Grigsby PW, Nene SM, Camel HM, Galakatos A, Kao MS, et al. Effect of Tumor Size on the Prognosis of Carcinoma of the Uterine Cervix Treated With Irradiation Alone. Cancer (1992) 69(11):2796–806. doi: 10.1002/1097-0142(19920601)69:11<2796::AID-CNCR2820691127>3.0.CO;2-O

CrossRef Full Text | Google Scholar

27. Elit LM, Hirte H. Management of Advanced or Recurrent Cervical Cancer: Chemotherapy and Beyond. Expert Rev Anticancer Ther (2014) 14(3):319–32. doi: 10.1586/14737140.2014.866041

CrossRef Full Text | Google Scholar

28. Wardak S: Human Papillomavirus (HPV) and Cervical Cancer. Medycyna doswiadczalna i mikrobiologia (2016) 68(1):73–84.

Google Scholar

29. Brianti P, De Flammineis E, Mercuri SR. Review of HPV-Related Diseases and Cancers. New Microbiol (2017) 40(2):80–5. doi: 10.1016/S0140-6736(13)60022-7

CrossRef Full Text | Google Scholar

30. Burd EM: Human Papillomavirus and Cervical Cancer. Clin Microbiol Rev (2003) 16(1):1–17. doi: 10.1128/CMR.16.1.1-17.2003

CrossRef Full Text | Google Scholar

31. Monk BJ, Willmott LJ, Sumner DA. Anti-Angiogenesis Agents in Metastatic or Recurrent Cervical Cancer. Gynecologic Oncol (2010) 116(2):181–6. doi: 10.1016/j.ygyno.2009.09.033

CrossRef Full Text | Google Scholar

32. Nakamura M, Bodily JM, Beglin M, Kyo S, Inoue M, Laimins LA. Hypoxia-Specific Stabilization of HIF-1alpha by Human Papillomaviruses. Virology (2009) 387(2):442–8. doi: 10.1016/j.virol.2009.02.036

CrossRef Full Text | Google Scholar

33. Brotherton JM, Fridman M, May CL, Chappell G, Saville AM, Gertig DM. Early Effect of the HPV Vaccination Programme on Cervical Abnormalities in Victoria, Australia: An Ecological Study. Lancet (2011) 377(9783):2085–92. doi: 10.1016/S0140-6736(11)60551-5

CrossRef Full Text | Google Scholar

34. Drolet M, Bénard É, Pérez N, Brisson M. Population-Level Impact and Herd Effects Following the Introduction of Human Papillomavirus Vaccination Programmes: Updated Systematic Review and Meta-Analysis. Lancet (2019) 394(10197):497–509. doi: 10.1016/S0140-6736(19)30298-3

CrossRef Full Text | Google Scholar

35. Gargano JW, Zhou F, Stokley S, Markowitz LE. Human Papillomavirus Vaccination in Commercially-Insured Vaccine-Eligible Males and Females, United States, 2007-2014. Vaccine (2018) 36(23):3381–6. doi: 10.1016/j.vaccine.2018.03.045

CrossRef Full Text | Google Scholar

36. Ventriglia J, Paciolla I, Pisano C, Cecere SC, Di Napoli M, Tambaro R, et al. Immunotherapy in Ovarian, Endometrial and Cervical Cancer: State of the Art and Future Perspectives. Cancer Treat Rev (2017) 59:109–16. doi: 10.1016/j.ctrv.2017.07.008

CrossRef Full Text | Google Scholar

37. Balar AV, Weber JS. PD-1 and PD-L1 Antibodies in Cancer: Current Status and Future Directions. Cancer Immunology Immunotherapy CII (2017) 66(5):551–64. doi: 10.1007/s00262-017-1954-6

CrossRef Full Text | Google Scholar

38. Topalian SL, Hodi FS, Brahmer JR, Gettinger SN, Smith DC, McDermott DF, et al. Safety, Activity, and Immune Correlates of Anti-PD-1 Antibody in Cancer. N Engl J Med (2012) 366(26):2443–54. doi: 10.1056/NEJMoa1200690

CrossRef Full Text | Google Scholar

39. Andrews LP, Marciscano AE, Drake CG, Vignali DA. LAG3 (CD223) as a Cancer Immunotherapy Target. Immunol Rev (2017) 276(1):80–96. doi: 10.1111/imr.12519

CrossRef Full Text | Google Scholar

40. Borghese F, Clanchy FI. CD74: An Emerging Opportunity as a Therapeutic Target in Cancer and Autoimmune Disease. Expert Opin Ther Targets (2011) 15(3):237–51. doi: 10.1517/14728222.2011.550879

CrossRef Full Text | Google Scholar

41. Rapp M, Wintergerst MWM, Kunz WG, Vetter VK, Knott MML, Lisowski D, et al. CCL22 Controls Immunity by Promoting Regulatory T Cell Communication With Dendritic Cells in Lymph Nodes. J Exp Med (2019) 216(5):1170–81. doi: 10.1084/jem.20170277

CrossRef Full Text | Google Scholar

42. Wang Q, Schmoeckel E, Kost BP, Kuhn C, Vattai A, Vilsmaier T, et al. Higher CCL22+ Cell Infiltration is Associated With Poor Prognosis in Cervical Cancer Patients. Cancers (2019) 11(12). doi: 10.3390/cancers11122004

CrossRef Full Text | Google Scholar

43. Ortiz A, Gui J, Zahedi F, Yu P, Cho C, Bhattacharya S, et al. An Interferon-Driven Oxysterol-Based Defense Against Tumor-Derived Extracellular Vesicles. Cancer Cell (2019) 35(1):33–45.e36. doi: 10.1016/j.ccell.2018.12.001

CrossRef Full Text | Google Scholar

44. Khaidakov M, Mitra S, Kang BY, Wang X, Kadlubar S, Novelli G, et al. Oxidized LDL Receptor 1 (OLR1) as a Possible Link Between Obesity, Dyslipidemia and Cancer. PloS One (2011) 6(5):e20277. doi: 10.1371/journal.pone.0020277

CrossRef Full Text | Google Scholar

45. Huang X, Gao Y, Qin J, Lu S. lncRNA MIAT Promotes Proliferation and Invasion of HCC Cells via sponging miR-214. Am J Physiol Gastrointestinal liver Physiol (2018) 314(5):G559–65. doi: 10.1152/ajpgi.00242.2017

CrossRef Full Text | Google Scholar

46. Sha M, Lin M, Wang J, Ye J, Xu J, Xu N, et al. Long Non-Coding RNA MIAT Promotes Gastric Cancer Growth and Metastasis Through Regulation of miR-141/DDX5 Pathway. J Exp Clin Cancer Res CR (2018) 37(1):58. doi: 10.1186/s13046-018-0725-3

CrossRef Full Text | Google Scholar

47. Liu W, Wang Z, Wang C, Ai Z. Long Non-Coding RNA MIAT Promotes Papillary Thyroid Cancer Progression Through Upregulating LASP1. Cancer Cell Int (2019) 19:194. doi: 10.1186/s12935-019-0913-z

CrossRef Full Text | Google Scholar

48. Kurachi M, Barnitz RA, Yosef N, Odorizzi PM, DiIorio MA, Lemieux ME, et al. The Transcription Factor BATF Operates as an Essential Differentiation Checkpoint in Early Effector CD8+ T Cells. Nat Immunol (2014) 15(4):373–83. doi: 10.1038/ni.2834

CrossRef Full Text | Google Scholar

49. Dai L, Cui X, Zhang X, Cheng L, Liu Y, Yang Y, et al. SARI Inhibits Angiogenesis and Tumour Growth of Human Colon Cancer Through Directly Targeting Ceruloplasmin. Nat Commun (2016) 7:11996. doi: 10.1038/ncomms11996

CrossRef Full Text | Google Scholar

50. Liu XP, He L, Zhang QP, Zeng XT, Liu SQ. Baicalein Inhibits Proliferation of Myeloma U266 Cells by Downregulating IKZF1 and IKZF3. Med Sci monitor Int Med J Exp Clin Res (2018) 24:2809–17. doi: 10.12659/MSM.907058

CrossRef Full Text | Google Scholar

51. Nie H, Bu F, Xu J, Li T, Huang J. 29 Immune-Related Genes Pairs Signature Predict the Prognosis of Cervical Cancer Patients. Sci Rep (2020) 10(1):14152. doi: 10.1038/s41598-020-70500-5

CrossRef Full Text | Google Scholar

52. Mei J, Xing Y, Lv J, Gu D, Pan J, Zhang Y, et al. Construction of an Immune-Related Gene Signature for Prediction of Prognosis in Patients With Cervical Cancer. Int Immunopharmacol (2020) 88:106882. doi: 10.1016/j.intimp.2020.106882

CrossRef Full Text | Google Scholar