Xiaoyan Shen
Xiao Gu
Ruiqiong Ma1†
Xiaoping Li
Jianliu Wang- 1Department of Gynecology, Peking University People’s Hospital, Beijing, China
- 2Department of Oncology, Shengjing Hospital of China Medical University, Shenyang, China
Ovarian cancer (OV) is a deadly gynecological cancer. The tumor immune microenvironment (TIME) plays a pivotal role in OV development. However, the TIME of OV is not fully known. Therefore, we aimed to provide a comprehensive network of the TIME in OV. Gene expression data and clinical information from OV patients were obtained from the Cancer Genome Atlas Program (TCGA) database. Non-negative Matrix Factorization, NMFConsensus, and nearest template prediction algorithms were used to perform molecular clustering. The biological functions of differentially expressed genes (DEGs) were identified using Metascape, gene set enrichment analysis (GSEA), gene ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. The copy number variations (CNVs), single nucleotide polymorphisms (SNPs) and tumor mutation burden were analyzed using Gistic 2.0, R package maftools, and TCGA mutations, respectively. Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data and CIBERSORT were utilized to elucidate the TIME. Moreover, external data from the International Cancer Genome Consortium (ICGC) and ArrayExpress databases were used to validate the signature. All 361 samples from the TCGA OV dataset were classified into Immune Class and non-Immune Class with immune signatures. By comparing the two classes, we identified 740 DEGs that accumulated in immune-related, cancer-related, inflammation-related biological functions and pathways. There were significant differences in the CNVs between the Immune and non-Immune Classes. The Immune Class was further divided into immune-activated and immune-suppressed subtypes. There was no significant difference in the top 20 genes in somatic SNPs among the three groups. In addition, the immune-activated subtype had significantly increased proportions of CD4 memory resting T cells, T cells, M1 macrophages, and M2 macrophages than the other two groups. The qRT-PCR results indicated that the mRNA expression levels of RYR2, FAT3, MDN1 and RYR1 were significantly down-regulated in OV compared with normal tissues. Moreover, the signatures of the TIME were validated using ICGC cohort and the ArrayExpress cohort. Our study clustered the OV patients into an immune-activated subtype, immune-suppressed subtype, and non-Immune Class and provided potential clues for further research on the molecular mechanisms and immunotherapy strategies of OV.
1 Introduction
Ovarian cancer (OV) is the third most commonly diagnosed cancer and the second leading cause of reproductive system death-related mortality in 2020 (Sung et al., 2021). Last year, 313,959 new cases of OV were reported, and 207,252 women died of this cancer worldwide. Despite the advances in surgery, chemotherapy, and target therapy for OV in recent decades, the 5-year survival rate for OV remains around 25%–35% (Lee et al., 2018). Therefore, research into potential molecular mechanisms and biomarkers to develop new anti-tumor therapy targets and to further improve the therapy response of OV patients is urgently needed.
Avoiding immune destruction is widely recognized as a key hallmark of cancer (Hanahan and Weinberg, 2011). The cancer cells can escape immune surveillance during tumor development (Chew et al., 2012; Liu and Cao, 2016) and evade immunological destruction by host immune systems (Mantovani and Sica, 2010; Hart et al., 2011; Veglia et al., 2018). Thus, many types of cancer, such as lung (Xia et al., 2021) and kidney cancer (Koh et al., 2021), are characterized by immune dysfunction. Recently, immunotherapy has emerged as a promising anti-tumor treatment strategy. Several studies have shown that immunotherapies are effective in clinical practice for cancer patients. For example, the patients with metastatic renal cell carcinoma responded to nivolumab (a programmed cell death protein 1 (PD-1) immune checkpoint inhibitor antibody) plus ipilimumab (a cytotoxic T-lymphocyte antigen-4 (CTLA-4) immune checkpoint inhibitor antibody) with acceptable toxicity (Gul et al., 2020). In patients with metastatic nonsquamous non-small cell lung cancer, pembrolizumab (an anti-PD-1 monoclonal antibody) plus pemetrexed-platinum improved overall survival (OS) and progression-free survival (PFS) with manageable safety (Gadgeel et al., 2020). In patients with advanced or metastatic esophageal squamous cell carcinoma, camrelizumab significantly improved OS, with a manageable safety (Huang et al., 2020).
Recently, evidence of a spontaneous anti-tumor immune response, satisfactory clinical response to immunotherapy, and immune evasion mechanisms have indicated that OV is an immunogenic cancer (Rodriguez-Garcia et al., 2017). Several studies have confirmed that tumor-infiltrating lymphocytes (TILs) are associated with a good prognosis of OV. Zhang et al. (Zhang et al., 2003) demonstrated that intratumoral T cells were observed in 54.8% of OV patients and were associated with a higher 5-year OS rate (38% vs. 4.5%). Other investigators showed that the patients with intraepithelial CD4+ and CD8+ TILs had better OS (CD4+ TILs, hazard ratio [HR] = 0.260; CD8+ TILs, HR = 0.503; both p < 0.05) and better PFS (CD4+ TILs, HR = 0.389; CD8+ TILs, HR = 0.478, both p < 0.005) than those without them (Pinto et al., 2018). Moreover, another study showed that 44.7% of OV patients had a high expression level of PD-L1 (PD-L1high), and patients with PD-L1high had a worse OS (HR = 2.877; p = 0.001) and PFS (HR = 1.843; p = 0.021) than those with low expression of PD-L1 (Zhu et al., 2017). M1 and M2 macrophages represent an important component of the tumor microenvironment (TME). M1 macrophages promoted OV metastasis through the activation of NF-κB pathway (Cho et al., 2018). Infiltration M2 macrophages was associated with poor prognosis of OV patients (Badmann et al., 2020). A high ratio of M1/M2 predicted a higher OS and PFS in OV patients (Maccio et al., 2020). Although these results offered the possibility that immunotherapy could be used for OV, the efficacy of immune checkpoint blockers in the clinic was still far from optimal. For example, the KEYNOTE-028 trial (Varga et al., 2019) showed that the objective response rate (ORR) in PD-L1-expressing advanced OV patients treated with pembrolizumab monotherapy was only 11.5%. On the other hand, another anti-PD-1 antibody, nivolumab, resulted in a 15% ORR in OV patients (Hamanishi et al., 2015). Based on the dilemma, the potential immunologic molecular mechanisms and the predictive biomarkers for immunotherapy efficacy in OV are urgently explored.
Emerging evidence suggested that the tumor immune microenvironment (TIME) changed and played a pivotal role in OV development and response to immunotherapies. For example, Huo et al. (2021) found ten TME genes related to the prognosis of OV patients. Olalekan et al. (2021) revealed the immune cell types and their roles in TME of metastatic OV by single-cell transcriptomics. Collagen type XI alpha 1 promoted OV growth and invasion by activating CAF (Wu et al., 2021). Multiple chemical agents, such as platinum derivatives, taxanes, and PARP inhibitors, regulate the interaction between tumor and stromal cells bidirectionally and extensively affect the TME (Eckert et al., 2021). Although many studies focus on the TIME of OV, the molecular mechanisms of TIME regulation of OV development remain unclear and require extensive research. In addition, a thorough understanding of the tumor immune microenvironment will guide the development of more effective immunotherapy targets for OV patients.
Therefore, we aimed to provide a comprehensive network of the immune microenvironment in OV. To achieve this goal, we divided the OV cohort of the Cancer Genome Atlas (TCGA) into Immune Class and non-Immune Class with immune signatures and analyzed the biological functions in each group. The Immune Class was further divided into immune-activated and immune-suppressed subtypes, and the gene signatures and TIME status were assessed in the three groups. Moreover, the Immune signatures of TIME were validated using the International Cancer Genome Consortium (ICGC) cohort and the ArrayExpress cohort.
2 Methods
2.1 Data Processing
Single-nucleotide polymorphisms (SNPs), clinical data, gene expression, and survival data of OV from the Cancer Genome Atlas Program (TCGA) were downloaded from the UCSC Xena platform. Gene expression data and clinical information of OV patients were obtained from the Australian Ovarian Cancer Study cohort (OV-AU). This study used TCGA data as the training data set and two OV datasets from International Cancer Genome Consortium (ICGC) and ArrayExpress as the external validation data sets.
The R package idmap1 annotated the expression data based on the Human Genome Organization’s Gene Nomenclature Committee (HGNC). We removed the mRNA expression data of the normal samples and, therefore, only expressions of tumor samples were retained to construct expression profiles for the following analyses.
We obtained the raw chip data and the corresponding annotations of the E-MTAB-62 dataset from ArrayExpress. The R packages affy (Gautier et al., 2004) and makecdfenv were used to read the raw chip data and generate the cdf package and environment. The normalized expression data using Robust Multi-array Average (RMA) method were annotated to the corresponding platform (GPL20967-3976), resulting in a gene symbol expression profile for subsequent analyses.
2.2 NMF Clustering
We used Non-negative Matrix Factorization (NMF) to perform molecular clustering and find the function features for each group. NMF factorizes a non-negative matrix V into two non-negative matrices, the base matrix W and the coefficient matrix H, so that V = W × H. In subsequent analyses, the coefficient matrices represent dimension-reduced matrices.
The TCGA mRNA expression profiles were subjected to NMF using the R package NMF in R environment version 4.0.5 (Gaujoux and Seoighe, 2010). The rank K of the maximum variation in cophenetic values was used as the optimal number of clustering.
2.3 Identification of Immune Microenvironment
We utilized the Estimation of STromal and Immune cells in MAlignant Tumor tissues using the Expression data (ESTIMATE; Becht et al., 2016) algorithm to predict the content of stromal and immune cells in each tumor sample. The stromal scores, immune scores, and ESTIMATE scores represent the ratio of stromal, immune, and the sum of both, respectively. We integrated the results of NMF and ESTIMATE to elucidate the immune microenvironments in each group.
2.4 Classification by NMFConsensus and NTP
We obtained and calculated the Immune Factors of the groups based on NMF and ESTIMATE. The top 150 genes of the Immune Factor were selected as Immune Genes to perform functional analysis using the Metascape database (Zhou et al., 2019). The Immune Gene expression profiles of all samples were analyzed by the NMFConsensus module of Genepattern (Reich et al., 2006) with default parameters. Samples were further divided into the immune group and non-immune group according to the result of NMFConsensus.
Nearest template prediction (NTP) is a convenient method for predicting clinical disease subtypes. To further study the immune features in OV, the GenePattern NTP module with default parameters was performed based on the activated stroma signature (Moffitt et al., 2015) of the immune group expression profile. The immune group was further divided into two subgroups: immune-suppressed and immune-activated subtypes.
2.5 Identification and Functional Analyses of Differentially Expressed Genes
The paired t-test of R package limma (Ritchie et al., 2015) was used to identify DEGs between the immune and non-immune groups. The absolute values of log2 [fold-change (FC)] > 1 and adjusted p-values < 0.01 were used to identify DEGs.
To identify the functions and relevant pathways of DEGs, we performed gene ontology (GO), and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis on DEGs using the Database for Annotation, Visualization and Integrated Discovery (DAVID) provides (Dennis et al., 2003). The GO terms of the biological process (BP) and KEGG pathways were visualized with the bubble plots. We carried out Gene Set Enrichment Analysis (GSEA) using the GSEA software to determine whether DEGs exhibit significant and consistent differences between the immune and non-immune groups (Subramanian et al., 2005).
2.6 Polymorphism and Tumor Mutation Burden Analyses
To analyze the copy number variation (CNV) mutations in the immune-suppressed, immune-activated and non-immune groups, we obtained CNV data for TCGA OV from GDAC Firehose. The Gistic 2.0 module from GenePattern with default parameters was conducted to identify regions with CNVs in immune-suppressed, immune-activated, and non-immune groups. The CNVs analysis results were collated and visualized using the R package maftools (Mayakonda et al., 2018). Analysis of Variance (ANOVA) tests were used to compare the CNVs among the three groups.
SNP mutation annotation format (MAF) files were downloaded from the UCSC Xena platform and processed using the R package maftools. We explored the top 20 genes of the somatic SNPs in the immune-suppressed, immune-activated, and non-immune groups by plotting the waterfall plots to identify the mutations of genes related to OV in the three groups.
We performed TMB analysis to identify the number of nonsynonymous somatic mutations in a specific genomic region in each group with TCGA OV mutation data using the R package TCGAmutations. The student’s t-test was used to compare TMB between the immune and non-immune groups.
2.7 TME Analysis
The analyses of the immune microenvironment are essential to identify the proportions of immune cells in tumor tissues. CIBERSORT (Newman et al., 2015) is used for the deconvolution of the expression matrix of human immune cell subtypes. This method provides gene expression signatures consisting of 22 immune cell subtypes (LM22) to estimate the abundances of immune cells. We applied CIBERSORT to analyze the immune microenvironment of expression matrices and obtained the proportion of 22 distinct immune cell types in OV tissues of different groups. We used ANOVA to compare the immune cell abundances among the three groups.
2.8 Sample Collection
18 OV tissues and 5 normal tissues were collected at People’s Hospital of Peking University between January 2021 and January 2022. OV tissues obtained from patients with serous adenocarcinoma. Normal ovarian tissues obtained from patients with benign gynecologic diseases who undergone bilateral salping-oophenrectomy. The protocol for this study was approved by the Ethical committee of People’s Hospital of Peking University (2021PHB239-001). Informed consent was obtained from all patients before surgery for residual tissue use.
2.9 Cell Lines and Cell Culture
OV cell lines SKOV3, A2780, ES2, and CAOV3 were used in the study. The human OV cells SKOV3 were incubated in McCoy’s 5A medium (Gibco, USA) contained with 10% fetal bovine serum (FBS) (Gibco, USA). The human OV cells A2780 and ES2 were incubated in RPMI 1640 medium (Gibco, USA) contained with 10% FBS. The human OV cells CAOV3 were incubated in Dulbecco’s modified Eagle medium (DMEM) medium (Gibco, USA) contained with 10% FBS. All cells were cultured at 37°C in 5% CO2.
2.10 Quantitative Real-Time PCR
Total RNA was collected with the TRIzol Reagent (Invitrogen, USA). The RNA was reversely transcribed into cDNA using a reverse transcription kit (Takara, Japan). qRT-PCR was conducted with SYBR Green PCR Master Mix (ABI, USA). GAPDH was used as a control. The expression level of genes was normalized to GAPDH. The relative expression of genes was calculated by the 2−ΔΔCT method. All results are presented as the mean ± standard deviation (SD). The primer sequences used were presented in Supplementary Table S1.
3 Results
3.1 Classification of OV Based on NMF and ESTIMATE
The total TCGA-OV expression profile data (361 samples) were grouped into 4 clusters by NMF. We estimate K from 2 to 15 and found that the cophenetic value changed the most when K varied from 4 to 5 (Figure 1A). Consequently, K = 4 was chosen as the optimal rank for NMF decomposition, and the decomposition results showed that the samples had stable and clear clusters (Figure 1B). We then conducted ESTIMATE analysis on all samples, and the results were integrated with the NMF results. Since cluster 2 (red module) from NMF analysis corresponded to a higher Immune Score from ESTIMATE analysis, we took cluster 2 as the Immune Factor for subsequent analysis (Figure 1C).
FIGURE 1. NMF and ESTIMATE analysis of OV dataset of TCGA mRNA expression profile. (A) The cophenetic values varied with K = 2 to 15 in NMF analysis. (B) The consensus map of NMF analysis results from the dimensional reduction of the original matrix and consensus analysis between modules. (C) The heat map of the NMF clustering results and ESTIMATE analysis results of OV. (D) Functional analysis of the top 150 Immune Genes of Immune Factor.
To ensure the accuracy of the results of Immune Factor and Immune Score from NMF analysis combined with ESTIMATE, we extracted the top 150 genes corresponding to cluster 2 from NMF analysis as Immune Genes to perform functional analysis using Metascape. The results showed that the top 150 Immune Genes from Immune Factor were significantly enriched in many immune-related biological functions, such as negative regulation of humoral immune response function and myeloid leukocyte activation function, further illustrating the accuracy of the results of Immune Factor and Immune Genes (Figure 1D).
All the samples were classified as Immune Class and non-Immune Class based on Immune Genes using NMFConsensus analysis. Figure 2A showed that the Immune Class contained almost all Immune Factors. Furthermore, the Immune Score of the Immune Class was higher than that of the non-Immune Class. In addition, principal component analysis (PCA) also showed that the Immune Class and the non-Immune Class had obvious characteristics (Figures 1B,C).
FIGURE 2. Classification of OV dataset into Immune Class and non-Immune Class. (A) NMFConsensus analysis with the top 150 Immune Genes. (B,C) The results of the Immune Class and the non-Immune Class based on principal component analysis (PCA). Each single point represents a sample.
3.2 Functional Analyses of DEGs Between Immune Class and Non-Immune Class
We conducted functional enrichment analyses to identify the biological functions involving the Immune and non-Immune Classes. We identified 740 DEGs by comparing Immune Class with non-Immune Class (adj. p-value < 0.01 and |log2(FC)| > 1), including 230 upregulated and 510 downregulated DEGs (Figures 3A,B). We conducted biological functions and pathways enrichment analysis of the 740 DEGs using DAVID. The top 15 GO terms were enriched with biological processes, such as inflammatory response, immune response regulation, cell-cell signaling, leukocyte migration, signal transduction, cellular defense response, chemical synaptic transmission, and adaptive immune response, shown in Figure 3C. The top 15 entries enriched in KEGG pathways include Staphylococcus aureus infection, cytokine-cytokine receptor interaction, protein digestion and absorption, ECM-receptor interaction, phagosome, primary immunodeficiency, and complement and coagulation cascades (Figure 3D).
FIGURE 3. DEGs and functional analysis of Immune Class and non-Immune Class. (A) The distribution of DEGs is shown in the volcano plot. (B) NMF clustering among all the DEGs into the Immune Class and the non-Immune Class. (C) The top 15 enriched GO items in BP were enriched among all DEGs. (D) The top 15 enriched KEGG pathways among all DEGs.
The identified DEGs, involved in many biological processes and pathways related to the immune system, suggested the reliability of the immune classification of Immune Genes by NMFConsensus. The DEGs were also involved in some interesting terms, such as extracellular matrix organization, cell adhesion, endodermal cell differentiation, signal transduction in biological processes, and neuroactive ligand-receptor interaction, nicotine addiction, protein digestion and absorption, ECM-receptor interaction, and cell adhesion molecules in KEGG pathways.
GSEA analysis showed that numerous pathways were related to immune cell and immune response enrichment in the Immune Class, such as B cell receptor signaling pathway, leukocyte transendothelial migration pathway, natural killer cell-mediated cytotoxicity pathway, T cell receptor signaling pathway, antigen processing and presentation pathway, chemokine signaling pathway, cytokine-cytokine receptor interaction pathway (Figure 4; Table 1). Some pathways related to proinflammatory were enriched in the Immune Class, such as NOD-like receptor signaling pathway and Toll-like receptor signaling pathway. There were several pathways related to tumor promotion enriched in the Immune Class, such as MAPK signaling pathway. The results of GO and KEGG pathway analyses of DEGs identified more function and pathways related to the immune system in the Immune Class than in the non-Immune Class.
FIGURE 4. The immune-related gene set enrichment in the GSEA analysis of Immune Class and non-Immune Class.
TABLE 1. Gene sets enriched in the Immune Class and non-Immune Class.
3.3 TCGA OV Dataset Characteristics Analysis
We applied the NTP algorithm to identify immune-suppressed and immune-activated subtypes in the Immune Class identified by integrated DEGs and classification analyses. The TCGA OV dataset was divided into the non-Immune Class, immune-suppressed, immune-activated subtypes according to their immune microenvironments. The decrease in immune-related scores and the increase in Tumor Purity from the Immune activated subtype and the Immune suppressed subtype to non-Immune subtype were observed (Figure 5A), which is consistent with the features of the three subtypes.
FIGURE 5. NTP, CNVs, and SNPs analyses among immune-suppressed subtype, immune-activated subtype, and non-Immune Class. (A) TCGA OV dataset was divided into the immune-suppressed subtype, the immune-activated subtype, and the non-Immune Class. (B) The TMB was estimated in the Immune Class and the non-Immune Class. (C) Statistical analyses of CNV amplification among the immune-suppressed subtype, the immune-activated subtype, and the non-Immune Class. (D) Statistical analysis of CNV deletion among the immune-suppressed subtype, the immune-activated subtype, and the non-Immune Class. (E–G) CNV analyses in the non-Immune Class (E), the immune-activated subtype (F) and the immune-suppressed subtype (G) (red for CNV amplification and blue for CNV deletion, some significant genes were marked). (H) Waterfall diagram of SNP analyses.
We analyzed CNVs, SNPs, and TMB among the subgroups to identify genetic variations among subgroups. CNV amplifications differ significantly between the Immune-suppressed subtype and non-Immune Class (p = 0.017) and Immune-activated subtype (p = 0.02) (Figure 5C). A difference in CNV deletion was observed between Immune Activated subtype and non-Immune Class (p = 0.013) (Figure 5D). There were significant differences in the amplification and deletion of genes between the non-Immune Class and the Immune Class. The significantly amplified genomic regions in both subtypes included 8q24.21 (with the highest G score), 8q24.3, 3q26.2, 19q12, 3q29. Meanwhile, the genomic regions with a significant deletion in both subtypes included 5q12.1, 19p13.3. The amplification of 19q12, 3q27.1, and 8q24.22 was frequently associated with the non-Immune Class, and the deletion of 18q23 was more frequently occurring in this class (Figures 5E–G). However, the non-Immune Class did not differ from the Immune Class in TMB (p = 0.55) (Figure 5B) and not from the Immune-activated subtype in CNVs amplification (p = 0.77). In addition, no significant difference of top 20 genes in TCGA OV somatic SNPs among the three groups was found (Figure 5H), which agrees with the TMB results (Figure 5B).
3.4 Different Tumor Immune Microenvironments Among Subtypes
We analyzed the tumor immune microenvironments in each subgroup with CIBERSORT. We found that the most common immune cell type was CD4 memory resting T cells, followed by M2 macrophages and M0 macrophages (Figure 6A). We further estimated the proportions of the specific 22 immune cells in different groups from the NMF analysis (Figure 6B).
FIGURE 6. Estimation of immune microenvironments in TCGA OV samples. (A) The proportions of specific 22 immune cells in the non-Immune Class, the Immune-suppressed subtype, and the Immune-activated subtype. (B) The proportions of specific 22 immune cells in each group from NMF analysis. (C–F) The differential proportions of memory CD4 resting T cells (C), CD8 T cells (D), M1 macrophages (E), and M2 macrophages (F) among the non-Immune Class, the Immune-suppressed subtype, and the Immune-activated subtype.
The proportions of the main immune cells were different among the non-Immune Class, the Immune-suppressed subtype, and the Immune-activated subtype. Compared with the Immune-suppressed subtype, we found higher proportions of CD4 memory resting T cells in the Immune-activated subtype (p = 0.025). Compared with the non-Immune Class, we found higher proportions of D8 T cells in the Immune-activated subtype (p = 0.0096). The Immune-activated subtype had a higher proportion of M1 macrophages than the Immune-suppressed subtype (p = 0.013) and the non-Immune Class (p = 0.021), and increased proportions of M2 macrophages than the Immune-suppressed subtype (p = 0.011) and the non-Immune Class (p = 0.00038) (Figures 6C–F).
3.5 External Validation of the Mutation Genes
To further validate the mutation genes at the mRNA level, qRT-PCR was used in on 18 OV tissues, 5 normal tissues, and 4 OV cells lines. The results indicated that the mRNA expression of RYR2, FAT3, MDN1 and RYR1 was significantly down-regulated in OV compared with normal tissues. The mRNA expression of MUC16, FLG, MYH4, HMCN1, FCGBP, LRP1B, LRP2 and AHNAK2 level tended to increase in the OV compared with the normal tissues but without statistical significance. The mRNA expression of TP53, TTN, FLG2, KMT2C, S1 and MACF1 level tended to decrease in the OV compared with the normal tissues but without statistical significance. The mRNA expression of USH2A was not detected in OV tissues. The expressions of genes in OV and normal tissues are shown in Figure 7.
FIGURE 7. The mRNA expression level of the genes in OV (n = 18) and normal ovarian tissues (n = 5) (*p < 0.05, **p < 0.01, ***p < 0.001).
The expressions of the mutation genes were validated by qRT-PCR in four OV cells (SKOV3, A2780, ES2 and CAOV3). The results showed that TTN, MUC16, MYH4, FAT3, KMT2C, LRP1B, RYR1, S1 and AHNAK2 were highly expressed in SKOV3 cell lines. TP53, FLG2, RYR2, FLG, HMCN1 and FCGBP were highly expressed in A2780 cell lines. CSMD3 and MDN1 were highly expressed in ES2 cell lines. LRP2 and MACF1 were highly expressed in CAOV3 cell lines. USH2A mRNA expression was not detected in all OV cells (Figure 8).
FIGURE 8. The mRNA expression level of the genes in OV ovarian cells (SKOV3, A2780, ES2 and CAOV3) (*p < 0.05, **p < 0.01, ***p < 0.001).
3.6 Validation of the Immune Microenvironment Analysis of TCGA OV Samples in ICGC
We used a total of 81 samples from ICGC OV-AU data as a validation set to confirm our results of the TCGA OV tumor microenvironment analyses.
The total ICGC OV-AU expression profiles were clustered into five groups by NMF. Then, we estimated the optimal K and found that the cophenetic value changed the most when K varied from 5 to 6 (Figure 9A). Therefore, the optimal rank of 5 was chosen for NMF decomposition (Figure 9B) with 13 samples in factor-1, 34 samples in factor-2, 13 samples in factor-3, 15 samples in factor-4, and six samples in factor-5, where factor-2 was considered to be the Immune factor.
FIGURE 9. Validation of the immune microenvironment analysis of TCGA OV samples using ICGC. (A) The cophenetic values varied with K = 2 to 15 in NMF analysis. (B) The consensus map of NMF analysis results from the dimensional reduction of original matrices and consensus analysis among modules. (C) The heat map of the NMF clustering results and ESTIMATE analysis results of OV. (D) All the samples were classified as Immune Class and non-Immune Class using NMFConsensus analysis. (E) ICGC OV-AU dataset was divided into the immune-suppressed subtype, the immune-activated subtype, and the non-Immune Class using NTP analysis. (F) The proportions of specific 22 immune cells in each group from NMF analysis. (G) The proportions of specific 22 immune cells in the non-Immune Class, the Immune-suppressed subtype, and the Immune-activated subtype. (H) The differential proportions of memory CD4 resting T cells among the non-Immune Class, the Immune-suppressed subtype, and the Immune-activated subtype.
ESTIMATE analysis, including all samples, was conducted and integrated with the NMF results to estimate Immune Factor (Figure 9C). We classified the validation set into an Immune Class with 54 samples and a non-Immune Class with 27 samples. Consistent with the TCGA OV analysis results, the Immune Class contained almost all Immune Factors and higher Immune Scores than the non-Immune Class (Figure 9D). The Immune Class was further grouped into the Immune-activate subtype with 26 samples and Immune-suppressed subtype with 28 samples by NTP analysis (Figure 9E).
To identify the immune microenvironment of OV samples in ICGC, we calculated the proportions of 22 immune cells in OV samples from ICGC. The most common immune cell type was M2 macrophages, followed by memory CD4 resting T cells (Figure 9G). We calculated the proportions of the specific 22 immune cells in each subgroup from (Figure 9F). The Immune-activated subtype had higher proportions of memory CD4 resting T cells than the Immune-suppressed subtype (p = 0.038) and the non-Immune Class (p = 0.00029) (Figure 9H), which is consistent with the results from TCGA OV analysis.
3.7 Validation of the Immune Microenvironment Analysis of TCGA OV Samples in ArrayExpress
The ArrayExpress gene symbol expression profiles were further used to validate our TCGA OV tumor microenvironment analysis results. We estimated the optimal K and found that the cophenetic value changed the most when K varied from 4 to 5 by NMF (Figure 10A). Therefore, the optimal rank of 4 was chosen for NMF decomposition, and 65 samples in factor-1, 62 samples in factor-2, 69 samples in factor-3, and 69 samples in factor-4 were identified by NMF analysis of 264 samples (Figure 10B), where factor-1 was considered to be the Immune factor. ESTIMATE analysis was conducted integrated with the NMF results to estimate Immune Factor (Figure 10C). We grouped 170 samples in the Immune Class and 95 samples in the non-Immune Class by NMFconsensus analysis. The Immune Class contained almost all Immune Factors and higher Immune Scores of the Immune Class than those of non-Immune Class (Figure 10D), which agrees with the results from TCGA OV analysis.
FIGURE 10. Validation the immune microenvironment analysis of TCGA OV samples with the gene symbol expression profiles from ArrayExpress. (A) The cophenetic values varied with K = 2 to 15 in NMF analysis. (B) The consensus map of NMF analysis results from the dimensional reduction of original matrices and consensus analysis among modules. (C) The heat map of the NMF clustering results and ESTIMATE analysis results of OV. (D) All the samples were classified as Immune Class and non-Immune Class using NMFConsensus analysis. (E) ArrayExpress gene symbol expression profiles were divided into the immune-suppressed subtype, the immune-activated subtype, and the non-Immune Class using NTP analysis. (F) The proportions of specific 22 immune cells in each group from NMF analysis. (G) The proportions of specific 22 immune cells in the non-Immune Class, the Immune-suppressed subtype, and the Immune-activated subtype. (H) The differential proportions of memory CD4 resting T cells among the non-Immune Class, the Immune-suppressed subtype, and the Immune-activated subtype.
The 170 samples in the Immune Class were further grouped into 92 samples as the Immune-activate subtype and 78 samples as the Immune-suppressed subtype (Figure 10E).
We calculated the proportions of 22 immune cells in OV samples for the ArrayExpress gene symbol expression profiles to identify the immune microenvironments. The most common immune cell type was memory CD4 resting T cells, followed by CD8 T cells (Figure 10G). We calculated the proportions of the specific 22 immune cells in each group from NMF analysis (Figure 10F). The Immune-suppressed subtype had higher proportions of M0 macrophages than the Immune-activated subtype (p = 0.000013) and the non-Immune Class (p = 0.0000057) (Figure 10H), consistent with the results from TCGA OV analysis.
4 Discussion
The current study aimed to comprehensively analyze the heterogeneous OV immune microenvironment subtypes underlying total immune genes. For this purpose, we analyzed 361 samples based on the TCGA OV dataset. We validated them using two external datasets, 81 samples from the ICGC OV dataset and 264 samples from the ArrayExpress OV dataset. We found the optimal rank as 4 in NMF analysis with the training set containing 361 samples in TCGA OV data. We obtained four factors: factor-1 had 74 samples, factor-2 had 80 samples, factor-3 had 69 samples, and factor-4 had 138 samples. Factor-2 was considered an Immune factor. We further grouped samples into two immune-related classes, among which the Immune Class had 244 samples and the non-Immune Class had 117 samples. The Immune Class was further divided into two subtypes by NTP analysis, of which the Immune-activate subtype had 121 samples, and the Immune-suppressed subtype had 123 samples. We identified 740 DEGs by comparing the two classes. The functional enrichment analysis of the 740 DEGs showed that immune-related, cancer-related, inflammation-related biological functions and pathways were enriched. There were significant differences in the amplification and deletion of genes between the Immune Class and non-Immune Class. The immune-activated subtype had increased proportions of memory CD4 resting T cells, CD8 T cells than the non-Immune Class, M1 macrophages, and M2 macrophages than the Immune-suppressed subtype and the non-Immune Class.
The Metascape database found the top 150 Immune Genes enriched in many immune-related GO items; the most significant enrichment item was the negative regulation of the humoral immune response. This result was consistent with cancer mediating immunosuppression and also suggested that our classification was accurate. Several studies demonstrated the biological function of humoral immunity in OV. A recent study found that CD19+ B cells within the total tumor were associated with better overall survival in high-grade serous OV than without CD19+ B cells. Tumor B-cell-derived IgA in the OV microenvironment induced OV cell death by redirecting myeloid cells against extracellular oncogenic drivers (Biswas et al., 2021). Most B cells presented in lymphoid structures in the stroma of high-grade serous OV metastases. They could improve the cytotoxic immune response to tumor cells through the secretion of cytokines and chemokines, which helped recruit and support antigen-presenting cells. Besides, immunoglobulins IgGs produced by B cells could target tumor antigens and form immune complexes that were helpful to activate antigen-presenting cells in the TME (Montfort et al., 2017). Thus, these results suggested that immunobiological pathways might be responsible for OV progression, and immunotherapy is a promising treatment strategy.
There were also many significant cancer-related GO enrichment elements in the top 150 Immune Genes, including integrin-mediated cell adhesion, histone H3-K9 modification, double-strand break repair by homologous recombination. Integrin plays an extensive role in cancer progression (Gu et al., 2017; Hou et al., 2020). In OV, it was reported that integrin is associated with proliferation, migration, invasion, chemoresistance, stemness, TIME, etc (Wu et al., 2020; Yan et al., 2021; Yin et al., 2021). The integrin-mediated adhesion was considered part of the metastasis. The methylation of histone H3K9 has been implicated in the development of various cancers due to its involvement in the transcriptional inactivation of chromatin and the induction of expression of cancer suppressor genes (Casciello et al., 2015). Homologous recombination is a key pathway involved in repairing DNA double-strand breaks. Konstantinopoulos et al. (2015) reported that around 50% of epithelial ovarian cancers are deficient in DNA repair via homologous recombination. New studies revealed that deficient DNA repair via homologous recombination is associated with immunotherapy response (Mouw et al., 2017). In the Keynote-162 (NCT02657889) phase I/II trial, 62 platinum-resistant recurrent OV patients were treated with pembrolizumab plus the inhibitor niraparib. The study reached an ORR (primary endpoint of the study) of 25% and a disease-control rate (DCR) of 68% (Konstantinopoulos et al., 2019). Therefore, the immune genes in OV might play various biological functions rather than immunity, indicating that cross-talks existed in immunobiological and other biological pathways.
Through the GO enrichment analysis of DEGs by comparing the Immune Class with the non-Immune Class, we found that many biological functions were highly enriched in immune-related processes, including immune response, immune response regulation, cell-cell signaling, leukocyte migration, signal transduction, cellular defense response, and adaptive immune response. Besides, KEGG enrichment analysis showed that the DEGs were significantly enriched in the hematopoietic cell lineage, the cytokine-cytokine receptor interaction, the primary immunodeficiency, and allograft rejection. These results indicated that DEGs played a potential role in regulating the TIME and greatly supported that the TIME was involved in OV malignant progression, consistent with previous studies (Westergaard et al., 2020; Shen et al., 2021). Moreover, it provided supporting evidence that both the Immune Class and the non-Immune Class existed in the OV immune microenvironment and differed significantly in their biological functions.
Inflammation is a hallmark of cancer, and the inflammatory state of premalignancy promotes tumor progression by various immune cells (Hanahan and Weinberg, 2011). In the present study, we performed GO analysis and identified inflammatory response as the most significant GO enrichment term of DEGs. KEGG pathway analysis showed that DEGs were significantly enriched in Staphylococcus aureus infection, phagosomes, asthma, viral myocarditis, and malaria. Numerous studies have shown that systemic inflammatory response markers, such as the neutrophil to lymphocyte ratio, and the platelet to lymphocyte ratio, could provide useful prognostic information among patients with OV (Zhang et al., 2017; Kwon et al., 2018; Yoshida et al., 2019). Interestingly, we also found that extracellular matrix organization and cell adhesion were enriched in GO analysis, and ECM-receptor interaction and cell adhesion molecules were enriched in KEGG analysis. These terms reflected the cross-talks that existed between inflammation and angiogenesis (Szewczyk et al., 2019) and between inflammation and metastasis (DiGiacomo and Gilkes, 2018). Angiogenesis is a common and critical biomarker of the inflammation-to-cancer transition in cancers. The infiltrated immune cells and their secreted cytokines were partly responsible for the inflammation-to-cancer transition. Furthermore, inflammation-to-cancer cytokines and angiogenesis genes might serve as predictors of survival and immune therapy response (Chen et al., 2019). Accordingly, it is suggested that the combination therapies of immunotherapy with anti-inflammatory therapy and anti-angiogenesis therapy might be possible and efficient.
In GSEA analysis, we discovered the NOD-like receptor signaling pathway, the Toll-like receptor signaling pathway, and many immune cells and immune response associated pathways significantly enriched in the Immune Class. A recent study demonstrates that activation of Toll-like receptor 8 reversed the immunosuppression function of regulatory T cells (Tregs) among TME by regulating glucose metabolism of Tregs in OV (Xu et al., 2021). Toll-like receptor 4 signaling pathway promoted OV cell proliferation and metastasis by activating osteopontin (Xu et al., 2017). Thus, it provided new insights for OV treatment and predictive prognosis.
The CNV-based risk score is an independent and discriminatory biomarker for the survival of OV patients (Graf et al., 2021). We also explored the CNV differences between the Immune Class and non-Immune Class. The amplification of 8q24.21 was detected in both Immune-activated and Immune-suppressed subtypes. The amplification of MYC, one of the genes located on 8q24.21, increased the sensitivity of OV cells to PARP inhibitors (Papp et al., 2018). TP53 was the most common mutated gene among all OV samples, and there was no difference among the three subtypes, likely due to the mutations in the TP53 gene in a high proportion (70%) of OV. This observation was consistent with a previous study (Cancer Genome Atlas Research, 2011). These results suggested the CNVs of the Immune Class might be closely associated with immunobiological pathways, further providing new targets for immunotherapy.
Notably, it has been previously reported that high levels of various immune cells are infiltrated into the TME of OV (Wang et al., 2018). Our independent analysis consistently discovered the immune cells expressed in tumors, especially memory CD4 resting T cells, macrophages M0, and macrophages M2. Memory CD4 resting T cells expression was associated with the prognosis in patients with OV (An and Yang, 2020). Lampert et al. (Lampert et al., 2020) found that memory CD4 resting T cells increased after treatment with a checkpoint 1 inhibitor of the cell cycle, perhaps due to the adaptive immune response activation. M0 macrophages conferred favorable OS, whereas the M2 macrophages were correlated with a poor OS in OV (Liu et al., 2020). M2 Macrophages inhibited apoptosis and increased the proliferation, invasion, and migration of OV cells (Yi et al., 2020). Recently, one study revealed that M2 macrophages controlled the vascular barrier by regulating the VCAM1/RAC1/ROS/p-PYK2/p-VE-cad axis. According to our results, the immune cells in the TIME played critical roles in OV progression, indicating the possibility of immunotherapy for OV.
Additionally, we also validated the aforementioned immune microenvironment analysis of TCGA OV samples in ICGC and ArrayExpress. All the results were as expected, suggesting that our results provided an accurate and repeatable classification of OV samples, as well as the potential targets for immunotherapy and biomarkers to predict outcomes.
OV is an aggressive malignancy that is still one of the most lethal gynecological cancers worldwide (Sung et al., 2021). The advanced stage was diagnosed in seventy-five percent of OV patients (Ferlay et al., 2015). The etiology of OV is unclear, and the major challenge in the clinical management of OV is the lack of effective treatment options. OV is currently stratified into different subtypes based on clinical and pathological characteristics. In contrast, the OS of different subtypes varies, indicating that biological heterogeneity still exists within each subtype. Therefore, it is necessary to investigate the underlying molecular subtypes according to specific gene patterns to evaluate and improve individualized medical decisions for OV patients. Taking advantage of the development of bioinformatics in recent years, especially the high-throughput sequencing technology, it is possible to systematically research the underlying mechanism of OV at the genome level.
OV is considered to be immunogenic, however, several clinical studies did not show a promising benefit for immunotherapy in OV recently (Moore et al., 2021; Pujade-Lauraine et al., 2021). Despite some OV patients show response to immunotherapy, there remains a subset of patients with PD-L1 expression who do not respond. Several markers have been found to be related to the efficacy of immunotherapy, such as TMB, PD-1, PD-L1, homologous repair deficient and proficient, TME and TILs (Conway et al., 2018; Keenan et al., 2019; Pellegrino et al., 2020; Plesca et al., 2020; Paijens et al., 2021). These factors were also been found in the present study. There are still lots of problems about the complex network of the immunogenicity of the cancer and the TIME needs to be further investigated, for example, how the immune system accesses the tumor, how the immune cells perform the killing functions, which is the ideal markers for response to immunotherapy. The present study improved understanding of immune interactions of OV and TIME.
We also performed external validation with the OV tissues and cells. The qRT-PCR results showed that the mRNA expression levels of FAT3, a putative tumor suppressor gene which codes for an atypical cadherin, were down-regulated in OV tissues. A study showed that FAT3 was enriched with strong mutations in metastases lesions in OV patients (Ojasalu et al., 2020). A recent study showed that the lung cancer patients with co-mutation of FAT3 and LRP1B had significantly prolonged immunotherapy PFS, which indicated that co-mutation of FAT3 and LRP1B is a promising biomarker to predict the efficacy of immunotherapy (Zhu et al., 2021).
However, our study may have some limitations. First, although three independent datasets were involved in the present study, more samples are still needed for comprehensive analysis to combat bias. Secondly, more clinical data should be included to analyze the immune signature and clinical characteristics comprehensively. Thirdly, more experimental evidence is needed to explain the molecular mechanisms and biological significance of our immunogenomic analysis.
In conclusion, we clustered the OV patients into three subtypes. The three subtypes of OV patients showed distinct gene signatures and TIME status, and the classification was efficient and repeatable. Moreover, novel functional genes and pathways that might contribute to the TME of OV were identified. Thus, the current study provided potential clues for further research on the molecular mechanisms and immunotherapy strategies of OV.
Data Availability Statement
Publicly available datasets were analyzed in this study. This data can be found here: TCGA (https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga), ICGC (https://dcc.icgc.org/), and ArrayExpress (https://www.ebi.ac.uk/arrayexpress/). The R code for this article is available at https://github.com/wangjianliu1203/Immune-signatures-for-ovarian-cancer.git.
Ethics Statement
The current work was supported by the Ethics Committee of People’s Hospital of Peking University. The patients/participants provided their written informed consent to participate in this study.
Author Contributions
JW designed the study. XS performed the data analysis and drew the pictures. XG wrote the manuscript. RM collected the OV samples, performed the qRT-PCR assay, and wrote the methods and results of corresponding portions. XL contributed to download the data. All authors have read and approved the final manuscript.
Funding
Project (RDY2018-02) supported by Peking University People’s Hospital Scientific Research Development Fund.
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.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2022.772701/full#supplementary-material
References
An, Y., and Yang, Q. (2020). Development and Validation of an Immune-Related Prognostic Signature for Ovarian Cancer Based on Weighted Gene Coexpression Network Analysis. Biomed. Res. Int. 2020, 1–17. doi:10.1155/2020/7594098
Badmann, S., Heublein, S., Mayr, D., Reischer, A., Liao, Y., Kolben, T., et al. (2020). M2 Macrophages Infiltrating Epithelial Ovarian Cancer Express MDR1: A Feature that May Account for the Poor Prognosis. Cells 9 (5), 1224. doi:10.3390/cells9051224
Becht, E., Giraldo, N. A., Lacroix, L., Buttard, B., Elarouci, N., Petitprez, F., et al. (2016). Estimating the Population Abundance of Tissue-Infiltrating Immune and Stromal Cell Populations Using Gene Expression. Genome Biol. 17 (1), 218. doi:10.1186/s13059-016-1070-5
Biswas, S., Mandal, G., Payne, K. K., Anadon, C. M., Gatenbee, C. D., Chaurio, R. A., et al. (2021). IgA Transcytosis and Antigen Recognition Govern Ovarian Cancer Immunity. Nature 591 (7850), 464–470. doi:10.1038/s41586-020-03144-0
Cancer Genome Atlas Research, N. (2011). Integrated Genomic Analyses of Ovarian Carcinoma. Nature 474 (7353), 609–615. doi:10.1038/nature10166
Casciello, F., Windloch, K., Gannon, F., and Lee, J. S. (2015). Functional Role of G9a Histone Methyltransferase in Cancer. Front. Immunol. 6, 487. doi:10.3389/fimmu.2015.00487
Chen, X., Xu, C., Hong, S., Xia, X., Cao, Y., McDermott, J., et al. (2019). Immune Cell Types and Secreted Factors Contributing to Inflammation-To-Cancer Transition and Immune Therapy Response. Cell Rep. 26 (7), 1965–1977. doi:10.1016/j.celrep.2019.01.080
Chew, V., Toh, H. C., and Abastado, J.-P. (20122012). Immune Microenvironment in Tumor Progression: Characteristics and Challenges for Therapy. J. Oncol. 2012, 1–10. doi:10.1155/2012/608406
Cho, U., Kim, B., Kim, S., Han, Y., and Song, Y. S. (2018). Pro-inflammatory M1 Macrophage Enhances Metastatic Potential of Ovarian Cancer Cells through NF-Κb Activation. Mol. Carcinog 57 (2), 235–242. doi:10.1002/mc.22750
Conway, J. R., Kofman, E., Mo, S. S., Elmarakeby, H., and Van Allen, E. (2018). Genomics of Response to Immune Checkpoint Therapies for Cancer: Implications for Precision Medicine. Genome Med. 10 (1), 93. doi:10.1186/s13073-018-0605-7
Dennis, G., Sherman, B. T., Hosack, D. A., Yang, J., Gao, W., Lane, H. C., et al. (2003). DAVID: Database for Annotation, Visualization, and Integrated Discovery. Genome Biol. 4 (5), P3. doi:10.1186/gb-2003-4-5-p3
DiGiacomo, J. W., and Gilkes, D. M. (2018). Tumor Hypoxia as an Enhancer of Inflammation-Mediated Metastasis: Emerging Therapeutic Strategies. Targ Oncol. 13 (2), 157–173. doi:10.1007/s11523-018-0555-4
Eckert, M. A., Orozco, C., Xiao, J., Javellana, M., and Lengyel, E. (2021). The Effects of Chemotherapeutics on the Ovarian Cancer Microenvironment. Cancers 13 (13), 3136. doi:10.3390/cancers13133136
Ferlay, J., Soerjomataram, I., Dikshit, R., Eser, S., Mathers, C., Rebelo, M., et al. (2015). Cancer Incidence and Mortality Worldwide: Sources, Methods and Major Patterns in GLOBOCAN 2012. Int. J. Cancer 136 (5), E359–E386. doi:10.1002/ijc.29210
Gadgeel, S., Rodríguez-Abreu, D., Speranza, G., Esteban, E., Felip, E., Dómine, M., et al. (2020). Updated Analysis from KEYNOTE-189: Pembrolizumab or Placebo Plus Pemetrexed and Platinum for Previously Untreated Metastatic Nonsquamous Non-small-cell Lung Cancer. Jco 38 (14), 1505–1517. doi:10.1200/JCO.19.03136
Gaujoux, R., and Seoighe, C. (2010). A Flexible R Package for Nonnegative Matrix Factorization. BMC Bioinformatics 11, 367. doi:10.1186/1471-2105-11-367
Gautier, L., Cope, L., Bolstad, B. M., and Irizarry, R. A. (2004). affy--analysis of Affymetrix GeneChip Data at the Probe Level. Bioinformatics 20 (3), 307–315. doi:10.1093/bioinformatics/btg405
Graf, R. P., Eskander, R., Brueggeman, L., and Stupack, D. G. (2021). Association of Copy Number Variation Signature and Survival in Patients with Serous Ovarian Cancer. JAMA Netw. Open 4 (6), e2114162. doi:10.1001/jamanetworkopen.2021.14162
Gu, X., Hua, Z., Dong, Y., Zhan, Y., Zhang, X., Tian, W., et al. (2017). Proteome and Acetylome Analysis Identifies Novel Pathways and Targets Regulated by Perifosine in Neuroblastoma. Sci. Rep. 7, 42062. doi:10.1038/srep42062
Gul, A., Stewart, T. F., Mantia, C. M., Shah, N. J., Gatof, E. S., Long, Y., et al. (2020). Salvage Ipilimumab and Nivolumab in Patients with Metastatic Renal Cell Carcinoma after Prior Immune Checkpoint Inhibitors. Jco 38 (27), 3088–3094. doi:10.1200/JCO.19.03315
Hamanishi, J., Mandai, M., Ikeda, T., Minami, M., Kawaguchi, A., Murayama, T., et al. (2015). Safety and Antitumor Activity of Anti-PD-1 Antibody, Nivolumab, in Patients with Platinum-Resistant Ovarian Cancer. Jco 33 (34), 4015–4022. doi:10.1200/JCO.2015.62.3397
Hanahan, D., and Weinberg, R. A. (2011). Hallmarks of Cancer: the Next Generation. Cell 144 (5), 646–674. doi:10.1016/j.cell.2011.02.013
Hart, K. M., Byrne, K. T., Molloy, M. J., Usherwood, E. M., and Berwin, B. (2011). IL-10 Immunomodulation of Myeloid Cells Regulates a Murine Model of Ovarian Cancer. Front. Immun. 2, 29. doi:10.3389/fimmu.2011.00029
Hou, J., Yan, D., Liu, Y., Huang, P., and Cui, H. (2020). The Roles of Integrin α5β1 in Human Cancer. Ott 13, 13329–13344. doi:10.2147/OTT.S273803
Huang, J., Xu, J., Chen, Y., Zhuang, W., Zhang, Y., Chen, Z., et al. (2020). Camrelizumab versus Investigator's Choice of Chemotherapy as Second-Line Therapy for Advanced or Metastatic Oesophageal Squamous Cell Carcinoma (ESCORT): a Multicentre, Randomised, Open-Label, Phase 3 Study. Lancet Oncol. 21 (6), 832–842. doi:10.1016/S1470-2045(20)30110-8
Huo, X., Sun, H., Liu, S., Liang, B., Bai, H., Wang, S., et al. (2021). Identification of a Prognostic Signature for Ovarian Cancer Based on the Microenvironment Genes. Front. Genet. 12, 680413. doi:10.3389/fgene.2021.680413
Keenan, T. E., Burke, K. P., and Van Allen, E. M. (2019). Genomic Correlates of Response to Immune Checkpoint Blockade. Nat. Med. 25 (3), 389–402. doi:10.1038/s41591-019-0382-x
Koh, M. Y., Sayegh, N., and Agarwal, N. (2021). Seeing the forest for the Trees-Single-Cell Atlases Link CD8+ T Cells and Macrophages to Disease Progression and Treatment Response in Kidney Cancer. Cancer Cell 39 (5), 594–596. doi:10.1016/j.ccell.2021.03.008
Konstantinopoulos, P. A., Ceccaldi, R., Shapiro, G. I., and D'Andrea, A. D. (2015). Homologous Recombination Deficiency: Exploiting the Fundamental Vulnerability of Ovarian Cancer. Cancer Discov. 5 (11), 1137–1154. doi:10.1158/2159-8290.CD-15-0714
Konstantinopoulos, P. A., Waggoner, S., Vidal, G. A., Mita, M., Moroney, J. W., Holloway, R., et al. (2019). Single-Arm Phases 1 and 2 Trial of Niraparib in Combination with Pembrolizumab in Patients with Recurrent Platinum-Resistant Ovarian Carcinoma. JAMA Oncol. 5 (8), 1141–1149. doi:10.1001/jamaoncol.2019.1048
Kwon, B. S., Jeong, D. H., Byun, J. M., Lee, T. H., Choi, K. U., Song, Y. J., et al. (2018). Prognostic Value of Preoperative Lymphocyte-Monocyte Ratio in Patients with Ovarian clear Cell Carcinoma. J. Cancer 9 (7), 1127–1134. doi:10.7150/jca.24057
Lampert, E. J., Cimino-Mathews, A., Lee, J. S., Nair, J., Lee, M.-J., Yuno, A., et al. (2020). Clinical Outcomes of Prexasertib Monotherapy in Recurrent BRCA Wild-type High-Grade Serous Ovarian Cancer Involve Innate and Adaptive Immune Responses. J. Immunother. Cancer 8 (2), e000516. doi:10.1136/jitc-2019-000516
Lee, J.-Y., Kim, S., Kim, Y. T., Lim, M. C., Lee, B., Jung, K.-W., et al. (2018). Changes in Ovarian Cancer Survival during the 20 Years before the Era of Targeted Therapy. BMC Cancer 18 (1), 601. doi:10.1186/s12885-018-4498-z
Liu, R., Hu, R., Zeng, Y., Zhang, W., and Zhou, H.-H. (2020). Tumour Immune Cell Infiltration and Survival after Platinum-Based Chemotherapy in High-Grade Serous Ovarian Cancer Subtypes: A Gene Expression-Based Computational Study. EBioMedicine 51, 102602. doi:10.1016/j.ebiom.2019.102602
Liu, Y., and Cao, X. (2016). Immunosuppressive Cells in Tumor Immune Escape and Metastasis. J. Mol. Med. 94 (5), 509–522. doi:10.1007/s00109-015-1376-x
Macciò, A., Gramignano, G., Cherchi, M. C., Tanca, L., Melis, L., and Madeddu, C. (2020). Role of M1-Polarized Tumor-Associated Macrophages in the Prognosis of Advanced Ovarian Cancer Patients. Sci. Rep. 10 (1), 6096. doi:10.1038/s41598-020-63276-1
Mantovani, A., and Sica, A. (2010). Macrophages, Innate Immunity and Cancer: Balance, Tolerance, and Diversity. Curr. Opin. Immunol. 22 (2), 231–237. doi:10.1016/j.coi.2010.01.009
Mayakonda, A., Lin, D.-C., Assenov, Y., Plass, C., and Koeffler, H. P. (2018). Maftools: Efficient and Comprehensive Analysis of Somatic Variants in Cancer. Genome Res. 28 (11), 1747–1756. doi:10.1101/gr.239244.118
Moffitt, R. A., Marayati, R., Flate, E. L., Volmar, K. E., Loeza, S. G. H., Hoadley, K. A., et al. (2015). Virtual Microdissection Identifies Distinct Tumor- and Stroma-specific Subtypes of Pancreatic Ductal Adenocarcinoma. Nat. Genet. 47 (10), 1168–1178. doi:10.1038/ng.3398
Montfort, A., Pearce, O., Maniati, E., Vincent, B. G., Bixby, L., Böhm, S., et al. (2017). A Strong B-Cell Response Is Part of the Immune Landscape in Human High-Grade Serous Ovarian Metastases. Clin. Cancer Res. 23 (1), 250–262. doi:10.1158/1078-0432.CCR-16-0081
Moore, K. N., Bookman, M., Sehouli, J., Miller, A., Anderson, C., Scambia, G., et al. (2021). Atezolizumab, Bevacizumab, and Chemotherapy for Newly Diagnosed Stage III or IV Ovarian Cancer: Placebo-Controlled Randomized Phase III Trial (IMagyn050/GOG 3015/ENGOT-OV39). Jco 39 (17), 1842–1855. doi:10.1200/JCO.21.00306
Mouw, K. W., Goldberg, M. S., Konstantinopoulos, P. A., and D'Andrea, A. D. (2017). DNA Damage and Repair Biomarkers of Immunotherapy Response. Cancer Discov. 7 (7), 675–693. doi:10.1158/2159-8290.CD-17-0226
Newman, A. M., Liu, C. L., Green, M. R., Gentles, A. J., Feng, W., Xu, Y., et al. (2015). Robust Enumeration of Cell Subsets from Tissue Expression Profiles. Nat. Methods 12 (5), 453–457. doi:10.1038/nmeth.3337
Ojasalu, K., Brehm, C., Hartung, K., Nischak, M., Finkernagel, F., Rexin, P., et al. (2020). Upregulation of Mesothelial Genes in Ovarian Carcinoma Cells Is Associated with an Unfavorable Clinical Outcome and the Promotion of Cancer Cell Adhesion. Mol. Oncol. 14 (9), 2142–2162. doi:10.1002/1878-0261.12749
Olalekan, S., Xie, B., Back, R., Eckart, H., and Basu, A. (2021). Characterizing the Tumor Microenvironment of Metastatic Ovarian Cancer by Single-Cell Transcriptomics. Cell Rep. 35 (8), 109165. doi:10.1016/j.celrep.2021.109165
Paijens, S. T., Vledder, A., de Bruyn, M., and Nijman, H. W. (2021). Tumor-infiltrating Lymphocytes in the Immunotherapy Era. Cell Mol Immunol 18 (4), 842–859. doi:10.1038/s41423-020-00565-9
Papp, E., Hallberg, D., Konecny, G. E., Bruhm, D. C., Adleff, V., Noë, M., et al. (2018). Integrated Genomic, Epigenomic, and Expression Analyses of Ovarian Cancer Cell Lines. Cell Rep. 25 (9), 2617–2633. doi:10.1016/j.celrep.2018.10.096
Pellegrino, B., Musolino, A., Llop-Guevara, A., Serra, V., De Silva, P., Hlavata, Z., et al. (2020). Homologous Recombination Repair Deficiency and the Immune Response in Breast Cancer: A Literature Review. Translational Oncol. 13 (2), 410–422. doi:10.1016/j.tranon.2019.10.010
Pinto, M. P., Balmaceda, C., Bravo, M. L., Kato, S., Villarroel, A., Owen, G. I., et al. (2018). Patient Inflammatory Status and CD4+/CD8+ Intraepithelial Tumor Lymphocyte Infiltration Are Predictors of Outcomes in High-Grade Serous Ovarian Cancer. Gynecol. Oncol. 151 (1), 10–17. doi:10.1016/j.ygyno.2018.07.025
Plesca, I., Tunger, A., Müller, L., Wehner, R., Lai, X., Grimm, M.-O., et al. (2020). Characteristics of Tumor-Infiltrating Lymphocytes Prior to and during Immune Checkpoint Inhibitor Therapy. Front. Immunol. 11, 364. doi:10.3389/fimmu.2020.00364
Pujade-Lauraine, E., Fujiwara, K., Ledermann, J. A., Oza, A. M., Kristeleit, R., Ray-Coquard, I.-L., et al. (2021). Avelumab Alone or in Combination with Chemotherapy versus Chemotherapy Alone in Platinum-Resistant or Platinum-Refractory Ovarian Cancer (JAVELIN Ovarian 200): an Open-Label, Three-Arm, Randomised, Phase 3 Study. Lancet Oncol. 22 (7), 1034–1046. doi:10.1016/S1470-2045(21)00216-3
Reich, M., Liefeld, T., Gould, J., Lerner, J., Tamayo, P., and Mesirov, J. P. (2006). GenePattern 2.0. Nat. Genet. 38 (5), 500–501. doi:10.1038/ng0506-500
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 (7), e47. doi:10.1093/nar/gkv007
Rodriguez-Garcia, A., Minutolo, N. G., Robinson, J. M., and Powell, D. J. (2017). T-cell Target Antigens across Major Gynecologic Cancers. Gynecol. Oncol. 145 (3), 426–435. doi:10.1016/j.ygyno.2017.03.510
Shen, J., Liu, T., Bei, Q., and Xu, S. (2021). Comprehensive Landscape of Ovarian Cancer Immune Microenvironment Based on Integrated Multi-Omics Analysis. Front. Oncol. 11, 685065. doi:10.3389/fonc.2021.685065
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 (43), 15545–15550. doi:10.1073/pnas.0506580102
Sung, H., Ferlay, J., Siegel, R. L., Laversanne, M., Soerjomataram, I., Jemal, A., et al. (2021). Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA A. Cancer J. Clin. 71 (3), 209–249. doi:10.3322/caac.21660
Szewczyk, G., Maciejewski, T. M., and Szukiewicz, D. (2019). Current Progress in the Inflammatory Background of Angiogenesis in Gynecological Cancers. Inflamm. Res. 68 (4), 247–260. doi:10.1007/s00011-019-01215-1
Varga, A., Piha-Paul, S., Ott, P. A., Mehnert, J. M., Berton-Rigaud, D., Morosky, A., et al. (2019). Pembrolizumab in Patients with Programmed Death Ligand 1-positive Advanced Ovarian Cancer: Analysis of KEYNOTE-028. Gynecol. Oncol. 152 (2), 243–250. doi:10.1016/j.ygyno.2018.11.017
Veglia, F., Perego, M., and Gabrilovich, D. (2018). Myeloid-derived Suppressor Cells Coming of Age. Nat. Immunol. 19 (2), 108–119. doi:10.1038/s41590-017-0022-x
Wang, W., Zou, W., and Liu, J. R. (2018). Tumor-infiltrating T Cells in Epithelial Ovarian Cancer: Predictors of Prognosis and Biological Basis of Immunotherapy. Gynecol. Oncol. 151 (1), 1–3. doi:10.1016/j.ygyno.2018.09.005
Westergaard, M. C. W., Milne, K., Pedersen, M., Hasselager, T., Olsen, L. R., Anglesio, M. S., et al. (2020). Changes in the Tumor Immune Microenvironment during Disease Progression in Patients with Ovarian Cancer. Cancers 12 (12), 3828. doi:10.3390/cancers12123828
Wu, A., Zhang, S., Liu, J., Huang, Y., Deng, W., Shu, G., et al. (2020). Integrated Analysis of Prognostic and Immune Associated Integrin Family in Ovarian Cancer. Front. Genet. 11, 705. doi:10.3389/fgene.2020.00705
Wu, Y.-H., Huang, Y.-F., Chang, T.-H., Chen, C.-C., Wu, P.-Y., Huang, S.-C., et al. (2021). COL11A1 Activates Cancer-Associated Fibroblasts by Modulating TGF-Β3 through the NF-Κb/igfbp2 axis in Ovarian Cancer Cells. Oncogene 40 (26), 4503–4519. doi:10.1038/s41388-021-01865-8
Xia, Y., Ying, S., Jin, R., Wu, H., Shen, Y., Yin, T., et al. (2021). Application of a Classifier Combining Bronchial Transcriptomics and Chest Computed Tomography Features Facilitates the Diagnostic Evaluation of Lung Cancer in Smokers and Nonsmokers. Int. J. Cancer 149 (6), 1290–1301. doi:10.1002/ijc.33675
Xu, C., Li, H., Yin, M., Yang, T., An, L., and Yang, G. (2017). Osteopontin Is Involved in TLR4 Pathway Contributing to Ovarian Cancer Cell Proliferation and Metastasis. Oncotarget 8 (58), 98394–98404. doi:10.18632/oncotarget.21844
Xu, R., Wu, M., Liu, S., Shang, W., Li, R., Xu, J., et al. (2021). Glucose Metabolism Characteristics and TLR8-Mediated Metabolic Control of CD4+ Treg Cells in Ovarian Cancer Cells Microenvironment. Cell Death Dis 12 (1), 22. doi:10.1038/s41419-020-03272-5
Yan, L., He, Z., Li, W., Liu, N., and Gao, S. (2021). The Overexpression of Acyl-CoA Medium-Chain Synthetase-3 (ACSM3) Suppresses the Ovarian Cancer Progression via the Inhibition of Integrin β1/AKT Signaling Pathway. Front. Oncol. 11, 644840. doi:10.3389/fonc.2021.644840
Yi, J., Lin, Y., Yicong, W., Chengyan, L., Shulin, Z., and Wenjun, C. (2020). Effect of Macrophages on Biological Function of Ovarian Cancer Cells in Tumor Microenvironment In Vitro. Arch. Gynecol. Obstet. 302 (4), 1009–1017. doi:10.1007/s00404-020-05719-8
Yin, H., Wang, J., Li, H., Yu, Y., Wang, X., Lu, L., et al. (2021). Extracellular Matrix Protein-1 Secretory Isoform Promotes Ovarian Cancer through Increasing Alternative mRNA Splicing and Stemness. Nat. Commun. 12 (1), 4230. doi:10.1038/s41467-021-24315-1
Yoshida, K., Yoshikawa, N., Shirakawa, A., Niimi, K., Suzuki, S., Kajiyama, H., et al. (2019). Prognostic Value of Neutrophil-To-Lymphocyte Ratio in Early-Stage Ovarian clear-cell Carcinoma. J. Gynecol. Oncol. 30 (6), e85. doi:10.3802/jgo.2019.30.e85
Zhang, H., Lu, J., Lu, Y., Zhou, J., Wang, Z., Liu, H., et al. (2017). Prognostic Significance and Predictors of the System Inflammation Score in Ovarian clear Cell Carcinoma. PLoS One 12 (5), e0177520. doi:10.1371/journal.pone.0177520
Zhang, L., Conejo-Garcia, J. R., Katsaros, D., Gimotty, P. A., Massobrio, M., Regnani, G., et al. (2003). Intratumoral T Cells, Recurrence, and Survival in Epithelial Ovarian Cancer. N. Engl. J. Med. 348 (3), 203–213. doi:10.1056/NEJMoa020177
Zhou, Y., Zhou, B., Pache, L., Chang, M., Khodabakhshi, A. H., Tanaseichuk, O., et al. (2019). Metascape Provides a Biologist-Oriented Resource for the Analysis of Systems-Level Datasets. Nat. Commun. 10 (1), 1523. doi:10.1038/s41467-019-09234-6
Zhu, J., Wen, H., Bi, R., Wu, Y., and Wu, X. (2017). Prognostic Value of Programmed Death-Ligand 1 (PD-L1) Expression in Ovarian clear Cell Carcinoma. J. Gynecol. Oncol. 28 (6), e77. doi:10.3802/jgo.2017.28.e77
Zhu, M., Zhang, L., Cui, H., Zhao, Q., Wang, H., Zhai, B., et al. (2021). Co-Mutation of FAT3 and LRP1B in Lung Adenocarcinoma Defines a Unique Subset Correlated with the Efficacy of Immunotherapy. Front. Immunol. 12, 800951. doi:10.3389/fimmu.2021.800951
Glossary
ANOVA analysis of Variance
BP biological process
CNVs copy number variations
CTLA-4 cytotoxic T-lymphocyte antigen-4
DAVID Database for Annotation, Visualization and Integrated Discovery
DCR disease-control rate
DEGs differentially expressed genes
DMEM Dulbecco’s modified Eagle medium
ESTIMATE Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data
FBS fetal bovine serum
FC fold-change
GO gene ontology
GSEA gene set enrichment analysis
HGNC Human Genome Organization’s Gene Nomenclature Committee
HR hazard ratio
ICGC International Cancer Genome Consortium
KEGG Kyoto Encyclopedia of Genes and Genomes
LM22 22 immune cell subtypes
MAF mutation annotation format
NMF Non-negative Matrix Factorization
NTP nearest template prediction
ORR objective response rate
OS overall survival
OV ovarian cancer
OV-AU Australian Ovarian Cancer Study cohort
PCA principal component analysis
PD-1 programmed cell death protein 1
PFS progression-free survival
qRT-PCR Quantitative Real-Time PCR
RMA Robust Multi-array Average
SD standard deviation
SNPs single nucleotide polymorphisms
TCGA the Cancer Genome Atlas Program
TILs tumor-infiltrating lymphocytes
TIME tumor immune microenvironment
TMB tumor mutation burden
TME tumor microenvironment
Keywords: ovarian cancer, immune, tumor microenvironment, molecular subtype, bioinformatics analysis
Citation: Shen X, Gu X, Ma R, Li X and Wang J (2022) Identification of the Immune Signatures for Ovarian Cancer Based on the Tumor Immune Microenvironment Genes. Front. Cell Dev. Biol. 10:772701. doi: 10.3389/fcell.2022.772701
Received: 08 September 2021; Accepted: 02 March 2022;
Published: 17 March 2022.
Edited by:
Xianquan Zhan, Shandong First Medical University, ChinaReviewed by:
Na Li, Central South University, ChinaEswari Dodagatta-Marri, University of California, San Francisco, United States
Peng-Chan Lin, National Cheng Kung University, Taiwan
Copyright © 2022 Shen, Gu, Ma, Li and Wang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jianliu Wang, wangjianliu@pkuph.edu.cn
†These authors have contributed equally to this work and share first authorship