Integrated Analysis of Immune Infiltration Features for Cervical Carcinoma and Their Associated Immunotherapeutic Responses

Introduction

Cervical cancer is the fourth most prevalent cancer among women. In the developing countries, it is the leading cause of cancer-associated mortalities (Canfell, 2019). The median age for patients diagnosed with cervical cancer is 49 years. In general, cervical cancer lowers the quality of life of the affected persons (Canfell, 2019). It has been established that prolonged infection with human papillomavirus (HPV) type 16 and 18 is a risk factor for cervical cancer (Crosbie et al., 2013). Prophylactic vaccines against high-risk HPV types minimize the risk of developing cervical cancer. However, due to the limitations associated with HPV vaccines, reliable therapeutic options for cervical cancer, particularly recurrent or advanced tumors, are required (Shanmugasundaram and You, 2017). Current therapeutic options include surgical removal of the tumors based on the FIGO staging system, incorporated with chemo- or radiotherapies (Bhatla et al., 2019; Koh et al., 2019). As for recurrent cervical cancer, bevacizumab in combination with other therapies can significantly prolong patients’ survival time (Tewari et al., 2017). In addition, immunotherapy is a viable option for cervical cancer treatment.

Immunotherapy is effective against various solid tumors. Mechanistically, immunotherapy enhances immune responses by utilizing immune checkpoint inhibitors and adoptive cellular transfer (O’Donnell et al., 2019). However, immunotherapeutic depends on the tumor microenvironment (Gasser et al., 2017). Immunocyte infiltration degree, tumor mutational load, and T cell functions affect tumor sensitivity to immunotherapy. Programmed death ligand 1 (PD-L1) has been reported in over 90% of cervical cancer and tumors. Higher infiltration ratios of CD8+ T cells and CD4+ T regulatory cells confer better survival outcomes for tumor patients (Ramanathan et al., 2018; Otter et al., 2019). The immune system is also involved in HPV-induced tumorigenesis (Orbegoso et al., 2018). HPV has been known to trigger chronic inflammation, escape immune surveillance by hiding in keratinocytes, suppress cellar immunity, and wall itself with recruited immunocytes (Piersma, 2011). Based on these features, immunotherapy presents the best strategy for managing cervical cancer (Piersma, 2011; Smola, 2017; Kagabu et al., 2019).

Establishment of reliable biomarkers for the best choice of immunotherapy as well as an improved understanding of immune infiltration features with regard to cervical cancer are key to immunotherapy. Therefore, this study aimed at exploring cervical cancer-induced immune infiltration characteristics based on different clustering algorithms, with the aim of providing a strong foundation for research and rationale for immunotherapy. The ratio of immunocytes, overall survival outcomes, mutation burden, and immunotherapeutic responses between different groups were compared.

Materials and Methods

CESC Data and Preprocessing

A total of 291 samples from publicly available cell carcinoma (CESC) gene-expression datasets in The Cancer Genome Atlas (TCGA) were utilized in our analyses (Supplementary Table 10). The TCGA dataset was downloaded from UCSC Xena¹, and subsequent analysis was performed using the R software (version 3.6.1) and R Bioconductor packages.

Estimation of TME Infiltrating Cells

The respective proportions of the immune infiltrating cells in the cervical squamous CESC samples were quantified using the ssgsea algorithm (Hanzelmann et al., 2013). The gene sets comprised of 782 genes that could predict the abundance of 28 TIICs in individual tissue samples². CESC was selected because it allows for the determination of sensitivity and specificity of immune cell phenotypes. It can be used to discriminate up to 28 human infiltration immune cell phenotypes.

Unsupervised Consensus Clustering of TME-Infiltrating Cells

In order to generate more groups for further analyses, Partitioning Around Medoid (PAM) (Tahiri et al., 2018) was used to classify tumors with qualitatively diverse TME-infiltrating patterns. The optimal number of clusters in the TCGA cohort was determined using the ConsensuClusterPlus R package (Monti et al., 2003). The consensus ESTIMATE algorithm was performed to assess the infiltration of stromal and immune cells in CESC samples (Yoshihara et al., 2013).

Identification of TME-Associated Differentially Expressed Genes

The patients were grouped into two distinct TME clusters based on the expression of immune infiltrating cells. The R package limma was used to determine differentially expressed genes (DEGs) between the two TME cell-infiltrating clusters (Ritchie et al., 2015). Adjusted p-value < 0.01 and | logFC| > 1 were considered to be statistically significant for DEGs between the TME subtypes.

TME Gene Signatures Generations and Dimension Reduction

The DEGs between TME clusters were standardized for all the samples in the TCGA CESC cohort. Prognostic associated DEGs were filtered out by performing the Univariate Cox regression analysis (p-value < 0.05). The unsupervised clustering method (Hartigan and Wong, 1979) was used to classify patients into either of the TME gene clusters. Annotation of the TME gene pattern was performed using the clusterProfiler R package (Ghasemi and Zahediasl, 2012). The clustering algorithm (Monti et al., 2003) was used to define the gene clusters. Principal component one that served as the signature score was obtained using the principal component analysis (PCA). The TME score for each patient was determined based on the prognostic value of the gene signature (Sotiriou et al., 2006):

$\text{TME}\text{score}=\sum \mathrm{PC1i}-\sum \mathrm{SPC1j}$

where “i” is the signature score of clusters with HR > 1, while “j” represents the expression of genes with HR < 1.

Pathway Enrichment Analysis

The gene sets for pathway enrichment analysis were downloaded from the MSigDB database (Subramanian et al., 2005). Gene set variation analysis (GSVA) was performed on the TME score and the TME clusters using the clusterProfiler R package (Ghasemi and Zahediasl, 2012). Genes for the enriched Pathways in TME were identified using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG), with an adjusted p < 0.05.

Immunotherapeutic Response Prediction

The Tumor Immune Dysfunction and Exclusion (TIDE) algorithm was used to link individual responses to immunotherapeutic responses (Jiang et al., 2018; Lu et al., 2019). Differences in anti-PD-1 and CTAL-4 therapeutic response were evaluated using the Submap analysis. For the melanoma data set (GSE78220, N = 28), the expression profiles (FPKM normalized) of GSE78220 were transformed into TPM values which were then used to calculate the TME score (Wagner et al., 2012). With regard to the urothelial cancer data set (IMvigor, N = 298), the data package was downloaded from http://research-pub.gene.com/IMvigor210CoreBiologies. Quality control and trimming of the mean of M-values were performed using the R package arrayQualityMetrics to normalize the numerical data (Ritchie et al., 2015).

Statistical Analysis

The Shapiro–Wilk normality test was used to establish variable normality (Ghasemi and Zahediasl, 2012). For normally distributed variables, the unpaired Student t-test was used to compare differences between the two groups, whereas the Wilcoxon test was used to compare abnormally distributed variables. One-way analysis of variance (ANOVA) and Kruskal–Wallis tests were used for comparison of multiple groups.

Pearson and distance correlation analyses were performed to calculate correlation coefficients. The χ² contingency test was performed to determine the interrelationships between TME score and anti-PD-1 response. The overall survival and TME score were determined using the R package. The threshold for survival values was determined. Based on the dichotomized TME score, patients were grouped into either high or low TME clusters, while at the same time reducing the computational batch effect by the R package sva. The data were visualized using the ggplot2 for R package. In the analysis of differential gene expression, we used the Benjamini–Hochberg method that converts p-values to FDRs to identify significantly expressed genes (Schreiber et al., 2011). OncoPrint was used to delineate the mutation landscape in the TCGA dataset using the maftools R package (Wang et al., 2020). Survival curves for the subgroups were generated using the Kaplan–Meier method. Statistical significance between different data sets was determined using the log-rank test. The univariate and multivariate Cox proportional hazard regression models were performed using the R package to determine independent factors associated with prognosis. Survivorship curves were generated using the R package survminer. Heatmaps were generated based on pheatmap. All statistical analyses were performed using R³. The tests were two-sided, with p-values < 0.05 being considered to be statistically significant.

Results

The Landscape and Functional Annotation of CESC TME

The flowchart for this study is shown in Supplementary Figure 1A. Analysis on cluster stability performed on CESC in the TCGA dataset using ConsensusClusterPlus package to select the optimal cluster number is shown in Supplementary Figure 1B. PAM of the 291 tumors with corresponding TME cell expression profiles in the TCGA cohort on its part is shown in Figure 1A. Two TME phenotypes were established based on immune cell infiltration. They conferred significantly different OS for outcomes (log-rank test, p < 0.001) as shown in Figure 1B. PCA showed a clear separation between the two established groups in the TCGA dataset (Figure 1C). Figure 1D shows the distinct TME infiltration patterns for the two clusters. Based on the ESTIMATE algorithm, TME cluster 1 was strongly associated with the estimated, immune, and stromal scores compared with TME cluster 2 (Figure 1E). Furthermore, 20 immune-related signaling pathways and DNA regulation-related pathways in GO analysis (Supplementary Table 1) were identified in the TCGA data set (Supplementary Figure 2A). Intrinsic immune escape was attributed to the expression of seven immune checkpoint molecules including antigen presenters, co-stimulators, co-inhibitors, receptors, ligands, and cell adhesion proteins, among others (Schreiber et al., 2011; Wang et al., 2020). There was an elevated expression of immune checkpoint molecules around TME cluster 1 of the CESC that aid the respective tumor cell escape from immune killing in TCGA (Supplementary Figure 2B). Moreover, the correlation of TME clusters, hypoxia, and metabolism was explored. TME cluster 2 was found to correlate with more metabolism-related signaling pathways (Figure 2A). TME cluster 2 also positively correlated with hypoxia-related gene signatures, indicating a malignancy of TME cluster 2 (Figure 2B).

FIGURE 1

Figure 1. The landscape of CESC TME, and the characteristics of TME subtypes in TCGA. (A) Unsupervised clustering of TME cells for 291 patients in TCGA. (B) Kaplan–Meier curves for 291 patients in TCGA, showing the association between TME infiltration patterns and OS (log-rank test, P < 0.001). (C) PCA for the two TME clusters. (D) The distribution of immunocytes in the two TME clusters. Within each group, the scattered dots represent values for cellular expression at TME whereas the thick line represents the median value. The bottom and top (lines) in the boxes are the 25th and 75th percentiles (interquartile range). The whiskers encompass 1.5 times the interquartile range. The statistical difference of the two TME clusters was compared using the Kruskal–Wallis test. **, P < 0.01; ****, P < 0.0001; ns, not statistically significant. (E) The differential expression of Estimate, Immune, and Stromal Score for TME clusters in the TCGA dataset.

FIGURE 2

Figure 2. Functional annotation of TME clusters. (A) Heatmap depicting the correlation of TME clusters and metabolic pathways. (B) Heatmap depicting the correlation of TME clusters and hypoxia-related gene signatures. *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001; ns, not statistically significant.

Generation of TME Gene Signatures and Functional Annotation

A total of 383 DEGs (Supplementary Table 2) were identified and used to classify the patients into genomic types, to further investigate the potential biological characteristics of each TME infiltration cell pattern. The clustering stability established by ConsensusClusterPlus package for the optimal number of clusters (Supplementary Figure 1C) was in tandem with the two CESC gene clusters (gene clusters 1 and 2) generated in TCGA (Figure 3A). Survival analysis of the two clusters revealed that the expression of gene cluster 1 was associated with better survival outcome (Figure 3B). Compared to cluster 2 genes, cluster 1 gene expressions were also associated with immune, stromal, and estimate scores (Figure 3C). Furthermore, cluster 1 genes were correlated with higher expression levels of infiltrating immune cells (Figure 3D). Compared to cluster 2 genes, CESC in gene cluster 1 were significantly associated with the immunosuppressive process that was mediated by a higher expression of immune checkpoint molecules (Figure 4A).

FIGURE 3

Figure 3. TME signatures and functional annotation constructs. (A) Unsupervised analysis and hierarchical clustering of common DEGs based on expression data of CESC derived from the TCGA: Gene clusters 1 and 2. (B) Kaplan–Meier curves for the two TME gene clusters (log-rank test showed an overall P < 0.001). (C) The differential expression of Estimate, Immune, and Stromal Score in TME clusters in the TCGA dataset. (D) The distribution of cells in TME gene clusters. *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001; ns, not statistically significant.

FIGURE 4

Figure 4. (A) The expression pattern of seven types of immune checkpoints in TME gene clusters for the TCGA dataset. (B) The GO enrichment analysis of the 98 DEG for TME signatures. (C) KEGG enrichment analysis of the 98 DEG for TME signatures. *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001; ns, not statistically significant.

Univariate Cox regression analysis for TME scores and the corresponding transcriptome traits as well as clinical characteristics for the top98 DEGs are shown in Supplementary Table 3, respectively. GO enrichment analysis revealed that 98 genes were associated with T cell activation and proliferation, regulation of macrophage activation, positive regulation of tumor necrosis factors, and neutrophil-mediated immunity. Combined, these pathways regulate the immune system (Figure 4B and Supplementary Table 4). KEGG enrichment analysis was consistent with GO analysis, with 98 genes found to be associated with immune system regulation (Figure 4C and Supplementary Table 5). TME scores for patients with CESC in the TCGA dataset are shown in Supplementary Table 6. GO analysis revealed that high TME scores were significantly associated with immune-related pathways, including T cell selection, T cell activation, regulation of macrophage activation, negative regulation of lymphocyte-mediated immunity, mast cell activation, regulation of myeloid dendritic cell activation, regulatory T cell differentiation, and immune response activation (Figure 5A and Supplementary Table 7). Furthermore, CESC with high TME scores exhibited higher immune checkpoint expression levels (Figure 5B). High TME scores were also associated with macrophage, mast cells, MDSC, and regulatory T cell infiltration, all of which induce an immunosuppressive environment. Moreover, a high TME score was also a predictor for a more activated immune environment. This is because high TME scores were correlated with a higher infiltration of multiple T cells and natural killer cells (Figure 6A). The two-sided role of high TME score could be attributed to the complexity of multiple immune cell-infiltrated tumor microenvironment. High TME scores were also positively correlated with Estimate, Immune, and Stromal Score (Figure 6B).

FIGURE 5

Figure 5. Immune-related functional annotation for the TME score. (A) GSVA for the TME score based on GO for TCGA. (B) The expression pattern of seven types of immune checkpoints in the TME score-based clusters for TCGA. *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001; ns, not statistically significant.

FIGURE 6

Figure 6. Genomic profiles associated with the TME score. (A) The distribution of immune cells in TME score clusters for TCGA. (B) The differential expression of Estimate, Immune, and Stromal Score for TME score in the TCGA dataset. (C) GISTIC 2.0 amplifications and deletions in CESC with high TME score. Chromosomal locations of peaks of significantly recurring focal amplification (red) and deletions (blue). (D) Differential somatic mutations in CESC with a high TME score. *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001; ns, not statistically significant.

TME Score Is Associated With Unique Patterns of Genomic Alterations

CNA and somatic mutation analysis performed on the TCGA dataset to determine the association between TME score and CESC genomic profiles revealed that samples with high TME scores frequently amplified several genomic regions particularly drivers of oncogenesis and immune regulatory genes including NRAS (1p13.27), DUP3Q29 (3q29), LYZ (5p11), HLA-DQA1 (6p21.32), CHEK2P2 (15q11.1), STAT3 (17q21.2), and KLK3 (19q13.33). These gene sets were associated with COL11A1 (1p21.1), MCL1 (1q21.2), UGT2B7 (4q13.2), ANGPT2 (8p23.1), PTEN (10q23.31) TNFRSF13B (17p11.2), TNNI3 (19q13.42), and GSTT1 (22q11.23) gene deletions as shown in Figure 6C. Different genomic profiles were observed in low TME score as shown in Supplementary Figure 3A. Furthermore, somatic mutation profiling revealed a high frequency of mutations in TTN (51%), MUC4 (35%), PIK3CA (32%), MUC17 (25%), MUC16 (24%), and SYNE1 (23%) among genes with high TME score (Figure 6D), whereas TTN (32%), MUC4 (29%), PIK3CA (29%), and MUC16 (27%) were the frequently mutated genes in the low TME score cluster (Supplementary Figure 3B).

TME Score Predicts Therapeutic Responses

In the analysis of the association between prognosis and TME score for the CESC cohort, it was found that high TME scores and targeted therapy were respective markers and intervention for positive prognosis of CESC in the TCGA dataset. Univariate and multivariate analyses revealed that patient ages, tumor stage, and radiation therapy were correlated with poor CESE prognosis (Figure 7A). A high TME score was associated with better survival outcomes for patients with CESC, BRCA, and OV in the TCGA dataset (Figure 7B). TME scores for BRCA and OV patients in the TCGA dataset included in this study are shown in Supplementary Tables 8, 9, respectively. Analysis of the probable response to immunotherapeutic response for CESE in the TCGA based on the TIDE algorithm revealed that low TME score tumors responded better to immunotherapy compared to high TME score tumors (Figure 7C). Subsequently, anti-PD-1 and anti-CTLA-4 therapeutic responses were analyzed. Tumors in the two TME score clusters were found to respond differently to immunotherapies. Tumors in the high TME score clusters exhibited a better response to anti- than those with low TME scores (Figure 7D). In the melanoma dataset, GSE78220 patients with high TME scores exhibited significantly longer OS outcomes than patients with lower TME scores (Figure 7E). High TME scores were also correlated with complete and partial anti-PD-1 responses (Figure 7F). In the urothelial cancer dataset, patients with high TME scores exhibited significantly longer OS outcomes than those with low TME scores in the IMvigor210 cohort (Figure 7G). Additionally, high TME scores were also associated with strong anti-PD-1 responses (Figure 7H).

FIGURE 7

Figure 7. Prognosis potential of TME score for immunotherapeutic response. (A) Univariate and multivariate Cox regression model estimating prognostic potential of TME score, patient age, tumor stage, tumor grade, pregnancy count, radiation therapy, and targeted therapy in TCGA. (The length of horizontal line represents the 95% confidence interval for each group. The vertical dotted line represents the HR of all patients. HR < 1.0 indicates that a high TME score is a biomarker for positive prognosis). (B) Kaplan–Meier curves for the two groups of patients classified along the TME score for CESC, BRCA, and OV from the TCGA. (Log-rank test showed an overall P < 0.001). (C) The TIDE value and response to immunotherapy for the TME score clusters. (D) Submap analysis based on the TIDE algorithm for differential response to CTAL-4 and anti-PD-1 therapy with regard to the TME score for the TCGA dataset. (E) Kaplan–Meier curves for the two groups of melanoma patients classified based on TME score for GSE78220. (Log-rank test showed an overall P < 0.001). (F) TME scores for groups with different anti–PD-1 clinical response status (CR/PR and SD/PD) in the GSE78220 dataset. (Wilcoxon, P = 0.019). (G) Kaplan–Meier curves for the two groups of patients classified based on the TME score in the IMvigor cohort. (Log-rank test: P < 0.001). (H) Distribution of TME scores in groups with different anti–PD-L1 clinical response statuses in the IMvigor cohort.

Discussion

The component of the cervical cancer microenvironment can affect cervical cancer progression, but potential mechanisms are still elusive. In this study, we analyzed the influence of immunocyte infiltration shed on tumor response to immunotherapy by applying different algorithms and proposed the subtype of cervical cancer based on immunocyte infiltration. CESC in the TCGA dataset were grouped into two clusters based on the differential expression of immunocytes, and lower immunocyte infiltration ratio samples showed worse survival outcome. In the meantime, several metabolic related pathways and hypoxia-associated genes were also differentially activated or expressed between the low-risk group (TME cluster 1) and the high-risk group (TME cluster 2). Therefore, the immune subtype of cervical cancer, TME cluster 1, and TME cluster 2 can affect the formation of the tumor microenvironment and cervical cancer progression.

To step further, 383 DEGs in total were identified and used to deeply explore this immune subtype of cervical cancer. High TME score group samples are usually accompanied with higher infiltration ratios of immunocytes like NK cells, T cells, and macrophages and higher immunotherapy-related gene expression (like HLA-A, HLA-B, HLA-C, CCL5, CXCL10, CD40, CTLA4, and PDCD1). In addition, immune activation-related pathways are also differentially activated in the low and high TME score groups. The previous study reported that chemokines like CCL5 and CXCL10 can modulate tumor sensitive to immunotherapy (Vilgelm and Richmond, 2019). Biomarkers like CTLA4, CD40, and the HLA family have been confirmed to be involved in immune surveillance (Waldman et al., 2020). Therefore, tumors with a low TME score group may be more sensitive to immunotherapy than high-scoring samples. However, its prediction ability still requires being examined with more clinical samples.

Previous studies also discussed the association between cervical cancer and its immune landscape. A previous study analyzed the proportion of immunocytes in cervical cancer and identified prognostic related immunocytes (Wang et al., 2019). Another study stratified samples based on the expression profile of differentially expressed immune-related genes and suggested that samples with higher infiltrated CD8 T cells and mast cells are more sensitive to immune checkpoint inhibitors (Yang et al., 2019). Moreover, tumor mutation loads are also critical intrinsic factors that affect tumor response to immunotherapy (Gasser et al., 2017; Havel et al., 2019). For instance, the amplification of HLA-DQA1 and STAT3 (Bae et al., 2020; Zou et al., 2020) and the deletion of ANGPT2 and TNFRSF13B (Rotolo et al., 2016; Lauret Marie Joseph et al., 2020) from the high TME score group have been proved to be immunotherapy-associated factors. Therefore, the TME score may serve as a potential tool to evaluate the tumor sensitivity to immunotherapy.

Previous studies proposed the “hot” and “cold” tumor analogy to describe tumor sensitivity to immunotherapy (Galon et al., 2006; Galon and Bruni, 2019). Given that tumors with high TME scores exhibited higher infiltrations of activated immunocytes and inflammatory related cells, tumors in this group may be referred as “hot” tumors. Moreover, high TME score tumors exhibited better PD-L1 receptor therapeutic response than low TME score tumors. Thus, adopting different strategies may improve patients’ clinical outcome. For instance, T-cell-targeted therapy (Buchbinder and Desai, 2016; Hellmann et al., 2016) or microbiome modulation (Snyder et al., 2015; Routy et al., 2018) was recommended to “hot” tumors. Chemotherapy, in combination with T cell enhancement or stimulatory signals, can improve “cold” tumor sensitivity (Whiteside et al., 2016). Taking the TME score into consideration in the choice of cervical cancer treatment may improve patients’ survival outcome.

In this study, infiltration of activated immunocytes was preferentially enhanced in the high TME score group. A higher immunocyte infiltration and enhanced immune checkpoint gene expression in the high TME score cluster implies that these tumors are more sensitive to immunotherapy. In conclusion, this study highlights the impact of the tumor microenvironment to immunotherapeutic sensitivity in tumors.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author Contributions

YK: data collection, analysis, and manuscript preparation. JH, YL, and NZ: manuscript revision. QC and YZ: study design. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by the China Postdoctoral Science Foundation (No. 2018M63302), the Hunan Provincial Natural Science Foundation of China (No. 2018JJ3838), the Hunan Provincial Health and Health Committee Foundation of China (C2019186), and Xiangya Hospital Central South University Postdoctoral Foundation.

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.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.573497/full#supplementary-material

Supplementary Figure 1 | Flow chart and cluster analysis related information. (A) The diagrammatical flow chart for this study. (B) TME clusters; TCGA matrixes for each k (k = 2–4), displaying the clustering stability using 1,000 iterations of hierarchical clustering. (C) TME gene clusters; TCGA matrixes for each k (k = 2–4), displaying the clustering stability using 1000 iterations of hierarchical clustering.

Supplementary Figure 2 | GSVA analysis and immune landscape based on TME clusters. (A) GSVA of TME clusters in TCGA-GO (B). The expression pattern of seven types of immune checkpoints in TME clusters for TCGA.

Supplementary Figure 3 | Genomic profiles associated with TME score. (A) GISTIC 2.0 amplifications and deletions in CESC with high TME score. Chromosomal locations of peaks of significantly recurring focal amplification (red) and deletions (blue). (B) Differential somatic mutations in CESC with high TME score.

Supplementary Table 1 | GO pathways for TME clusters based on GSVA in TCGA.

Supplementary Table 2 | 383 differentially expressed genes in TCGA.

Supplementary Table 3 | 98 representative TME signature genes.

Supplementary Table 4 | GO enrichment analysis for TME gene clusters based on 98 genes.

Supplementary Table 5 | KEGG enrichment analysis for TME gene clusters based on 98 genes.

Supplementary Table 6 | TME score for CESC patients in TCGA.

Supplementary Table 7 | GSVA of TME score.

Supplementary Table 8 | TME score for BRCA patients.

Supplementary Table 9 | TME score for OV patients.

Supplementary Table 10 | Characteristics of all the included patient samples.

Footnotes

1. ^ https://xenabrowser.net/
2. ^ http://software.broadinstitute.org/gsea/msigdb/index.jsp
3. ^ https://www.r-project.org/, version 3.6.1

References

Bae, J. Y., Choi, K. U., Kim, A., Lee, S. J., Kim, K., Kim, J. Y., et al. (2020). Evaluation of immune-biomarker expression in high-grade soft-tissue sarcoma: HLA-DQA1 expression as a prognostic marker. Exp. Ther. Med. 20:107. doi: 10.3892/etm.2020.9225

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhatla, N., Berek, J. S., Cuello Fredes, M., Denny, L. A., Grenman, S., Karunaratne, K., et al. (2019). Revised FIGO staging for carcinoma of the cervix uteri. Int. J. Gynaecol. Obstet. 145, 129–135. doi: 10.1002/ijgo.12749

PubMed Abstract | CrossRef Full Text | Google Scholar

Buchbinder, E. I., and Desai, A. (2016). CTLA-4 and PD-1 Pathways: similarities, differences, and implications of their inhibition. Am. J. Clin. Oncol. 39, 98–106. doi: 10.1097/COC.0000000000000239

PubMed Abstract | CrossRef Full Text | Google Scholar

Canfell, K. (2019). Towards the global elimination of cervical cancer. Papillomavirus Res. 8:100170. doi: 10.1016/j.pvr.2019.100170

PubMed Abstract | CrossRef Full Text | Google Scholar

Crosbie, E. J., Einstein, M. H., Franceschi, S., and Kitchener, H. C. (2013). Human papillomavirus and cervical cancer. Lancet 382, 889–899. doi: 10.1016/S0140-6736(13)60022-7

CrossRef Full Text | Google Scholar

Galon, J., and Bruni, D. (2019). Approaches to treat immune hot, altered and cold tumours with combination immunotherapies. Nat. Rev. Drug Discov. 18, 197–218. doi: 10.1038/s41573-018-0007-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Galon, J., Costes, A., Sanchez-Cabo, F., Kirilovsky, A., Mlecnik, B., Lagorce-Pages, C., et al. (2006). Type, density, and location of immune cells within human colorectal tumors predict clinical outcome. Science 313, 1960–1964. doi: 10.1126/science.1129139

PubMed Abstract | CrossRef Full Text | Google Scholar

Gasser, S., Lim, L. H. K., and Cheung, F. S. G. (2017). The role of the tumour microenvironment in immunotherapy. Endocr. Relat. Cancer 24, T283–T295. doi: 10.1530/ERC-17-0146

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghasemi, A., and Zahediasl, S. (2012). Normality tests for statistical analysis: a guide for non-statisticians. Int. J. Endocrinol. Metab. 10, 486–489. doi: 10.5812/ijem.3505

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 14:7. doi: 10.1186/1471-2105-14-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Hartigan, J. A., and Wong, M. A. (1979). Algorithm AS 136: A K-Means clustering algorithm. J. Royal Stat. Soc. Ser. C (Appl. Stat.) 28, 100–108. doi: 10.2307/2346830

CrossRef Full Text | Google Scholar

Havel, J. J., Chowell, D., and Chan, T. A. (2019). The evolving landscape of biomarkers for checkpoint inhibitor immunotherapy. Nat. Rev. Cancer 19, 133–150. doi: 10.1038/s41568-019-0116-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hellmann, M. D., Friedman, C. F., and Wolchok, J. D. (2016). Combinatorial cancer immunotherapies. Adv. Immunol. 130, 251–277. doi: 10.1016/bs.ai.2015.12.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, P., Gu, S., Pan, D., Fu, J., Sahu, A., Hu, X., et al. (2018). Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat. Med. 24, 1550–1558. doi: 10.1038/s41591-018-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Kagabu, M., Nagasawa, T., Fukagawa, D., Tomabechi, H., Sato, S., Shoji, T., et al. (2019). Immunotherapy for uterine cervical cancer. Healthcare (Basel) 7:108. doi: 10.3390/healthcare7030108

PubMed Abstract | CrossRef Full Text | Google Scholar

Koh, W. J., Abu-Rustum, N. R., Bean, S., Bradley, K., Campos, S. M., Cho, K. R., et al. (2019). Cervical cancer, version 3.2019, NCCN clinical practice guidelines in oncology. J. Natl. Compr. Canc. Netw. 17, 64–84. doi: 10.6004/jnccn.2019.0001

PubMed Abstract | CrossRef Full Text | Google Scholar

Lauret Marie Joseph, E., Laheurte, C., Jary, M., Boullerot, L., Asgarov, K., Gravelin, E., et al. (2020). Immunoregulation and clinical implications of ANGPT2/TIE2(+) M-MDSC signature in non-small cell lung cancer. Cancer Immunol. Res. 8, 268–279. doi: 10.1158/2326-6066.CIR-19-0326

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, X., Jiang, L., Zhang, L., Zhu, Y., Hu, W., Wang, J., et al. (2019). Immune signature-based subtypes of cervical squamous cell carcinoma tightly associated with human papillomavirus type 16 expression, molecular features, and clinical outcome. Neoplasia 21, 591–601. doi: 10.1016/j.neo.2019.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Monti, S., Tamayo, P., Mesirov, J., and Golub, T. (2003). Consensus clustering: a resampling-based method for class discovery and visualization of gene expression microarray data. Mach. Learn. 52, 91–118.

Google Scholar

O’Donnell, J. S., Teng, M. W. L., and Smyth, M. J. (2019). Cancer immunoediting and resistance to T cell-based immunotherapy. Nat. Rev. Clin. Oncol. 16, 151–167. doi: 10.1038/s41571-018-0142-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Orbegoso, C., Murali, K., and Banerjee, S. (2018). The current status of immunotherapy for cervical cancer. Rep. Pract. Oncol. Radiother. 23, 580–588. doi: 10.1016/j.rpor.2018.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Otter, S. J., Chatterjee, J., Stewart, A. J., and Michael, A. (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. 31, 834–843. doi: 10.1016/j.clon.2019.07.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Piersma, S. J. (2011). Immunosuppressive tumor microenvironment in cervical cancer patients. Cancer Microenviron. 4, 361–375. doi: 10.1007/s12307-011-0066-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramanathan, P., Dhandapani, H., Jayakumar, H., Seetharaman, A., and Thangarajan, R. (2018). Immunotherapy for cervical cancer: can it do another lung cancer? Curr. Probl. Cancer 42, 148–160. doi: 10.1016/j.currproblcancer.2017.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

Rotolo, A., Caputo, V., and Karadimitris, A. (2016). The prospects and promise of chimeric antigen receptor immunotherapy in multiple myeloma. Br. J. Haematol. 173, 350–364. doi: 10.1111/bjh.13976

PubMed Abstract | CrossRef Full Text | Google Scholar

Routy, B., Le Chatelier, E., Derosa, L., Duong, C. P. M., Alou, M. T., Daillere, R., et al. (2018). Gut microbiome influences efficacy of PD-1-based immunotherapy against epithelial tumors. Science 359, 91–97. doi: 10.1126/science.aan3706

PubMed Abstract | CrossRef Full Text | Google Scholar

Schreiber, R. D., Old, L. J., and Smyth, M. J. (2011). Cancer immunoediting: integrating immunity’s roles in cancer suppression and promotion. Science 331, 1565–1570. doi: 10.1126/science.1203486

PubMed Abstract | CrossRef Full Text | Google Scholar

Shanmugasundaram, S., and You, J. (2017). Targeting persistent human papillomavirus infection. Viruses 9:229. doi: 10.3390/v9080229

PubMed Abstract | CrossRef Full Text | Google Scholar

Smola, S. (2017). Immunopathogenesis of HPV-Associated cancers and prospects for immunotherapy. Viruses 9:254. doi: 10.3390/v9090254

PubMed Abstract | CrossRef Full Text | Google Scholar

Snyder, A., Pamer, E., and Wolchok, J. (2015). IMMUNOTHERAPY. could microbial therapy boost cancer immunotherapy? Science 350, 1031–1032. doi: 10.1126/science.aad7706

PubMed Abstract | CrossRef Full Text | Google Scholar

Sotiriou, C., Wirapati, P., Loi, S., Harris, A., Fox, S., Smeds, J., et al. (2006). Gene expression profiling in breast cancer: understanding the molecular basis of histologic grade to improve prognosis. J. Natl. Cancer Inst. 98, 262–272. doi: 10.1093/jnci/djj052

PubMed Abstract | CrossRef Full Text | Google Scholar

Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U.S.A. 102, 15545–15550. doi: 10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

Tahiri, N., Willems, M., and Makarenkov, V. (2018). A new fast method for inferring multiple consensus trees using k-medoids. BMC Evol. Biol. 18:48. doi: 10.1186/s12862-018-1163-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Tewari, K. S., Sill, M. W., Penson, R. T., Huang, H., Ramondetta, L. M., Landrum, L. M., et al. (2017). Bevacizumab for advanced cervical cancer: final overall survival and adverse event analysis of a randomised, controlled, open-label, phase 3 trial (Gynecologic Oncology Group 240). Lancet 390, 1654–1663. doi: 10.1016/S0140-6736(17)31607-0

CrossRef Full Text | Google Scholar

Vilgelm, A. E., and Richmond, A. (2019). Chemokines modulate immune surveillance in tumorigenesis, metastasis, and response to immunotherapy. Front. Immunol. 10:333. doi: 10.3389/fimmu.2019.00333

PubMed Abstract | CrossRef Full Text | Google Scholar

Wagner, G. P., Kin, K., and Lynch, V. J. (2012). Measurement of mRNA abundance using RNA-seq data: RPKM measure is inconsistent among samples. Theory Biosci. 131, 281–285. doi: 10.1007/s12064-012-0162-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Waldman, A. D., Fritz, J. M., and Lenardo, M. J. (2020). A guide to cancer immunotherapy: from T cell basic science to clinical practice. Nat. Rev. Immunol. 20, 651–668. doi: 10.1038/s41577-020-0306-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Li, Z., Gao, A., Wen, Q., and Sun, Y. (2019). The prognostic landscape of tumor-infiltrating immune cells in cervical cancer. Biomed. Pharmacother. 120:109444. doi: 10.1016/j.biopha.2019.109444

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, S., Zhang, Q., Yu, C., Cao, Y., Zuo, Y., and Yang, L. (2020). Immune cell infiltration-based signature for prognosis and immunogenomic analysis in breast cancer. Brief Bioinform. 6:bbaa026. doi: 10.1093/bib/bbaa026

PubMed Abstract | CrossRef Full Text | Google Scholar

Whiteside, T. L., Demaria, S., Rodriguez-Ruiz, M. E., Zarour, H. M., and Melero, I. (2016). Emerging opportunities and challenges in cancer immunotherapy. Clin. Cancer Res. 22, 1845–1855. doi: 10.1158/1078-0432.CCR-16-0049

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, S., Wu, Y., Deng, Y., Zhou, L., Yang, P., Zheng, Y., et al. (2019). Identification of a prognostic immune signature for cervical cancer to predict survival and response to immune checkpoint inhibitors. Oncoimmunology 8:e1659094. doi: 10.1080/2162402X.2019.1659094

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoshihara, K., Shahmoradgoli, M., Martínez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun. 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

Zou, S., Tong, Q., Liu, B., Huang, W., Tian, Y., and Fu, X. (2020). Targeting STAT3 in Cancer immunotherapy. Mol. Cancer 19:145. doi: 10.1186/s12943-020-01258-7

PubMed Abstract | CrossRef Full Text | Google Scholar