Landscape of Immune Microenvironment in Epithelial Ovarian Cancer and Establishing Risk Model by Machine Learning

Background Epithelial ovarian cancer (EOC) is an extremely lethal gynecological malignancy and has the potential to benefit from the immune checkpoint blockade (ICB) therapy, whose efficacy highly depends on the complex tumor microenvironment (TME). Method and Result We comprehensively analyze the landscape of TME and its prognostic value through immune infiltration analysis, somatic mutation analysis, and survival analysis. The results showed that high infiltration of immune cells predicts favorable clinical outcomes in EOC. Then, the detailed TME landscape of the EOC had been investigated through “xCell” algorithm, Gene set variation analysis (GSVA), cytokines expression analysis, and correlation analysis. It is observed that EOC patients with high infiltrating immune cells have an antitumor phenotype and are highly correlated with immune checkpoints. We further found that dendritic cells (DCs) may play a dominant role in promoting the infiltration of immune cells into TME and forming an antitumor immune phenotype. Finally, we conducted machine-learning Lasso regression, support vector machines (SVMs), and random forest, identifying six DC-related prognostic genes (CXCL9, VSIG4, ALOX5AP, TGFBI, UBD, and CXCL11). And DC-related risk stratify model had been well established and validated. Conclusion High infiltration of immune cells predicted a better outcome and an antitumor phenotype in EOC, and the DCs might play a dominant role in the initiation of antitumor immune cells. The well-established risk model can be used for prognostic prediction in EOC.


Introduction
Ovarian cancer is the second leading cause of cancer death in adult women, and the epithelial ovarian cancer (EOC) is the most common histological subtype characterized by complex adjacent anatomical structures, high-degree malignancy, and heterogeneity [1,2]. e surgical excision and chemotherapy are effective for EOC patients, but a large fraction of EOC patients will subsequently relapse and develop to chemoresistance, which seriously shortens the patients' clinical survival [3].
Nowadays, the immune checkpoint blockades (ICBs) have received wide attention and have emerged as efficient agents for tumor therapy. It can specifically block the immune checkpoints and retrieve the antitumor immunity [4]. Nevertheless, due to the lack of comprehensive analysis on the tumor microenvironment (TME) and EOC, there is a limited application of ICBs in EOC therapy. e TME not only plays a critical role in carcinogenesis but also has various regulatory functions in tumor growth and metastasis [5]. It has been proposed to be valuable in diagnosis and prognosis in a wide range of tumors, but the utility on EOC has not been investigated in detail [6]. e immune cells are a master component in TME, composed of CD8 + cytotoxic T lymphocytes (CTLs), B cells, plasma cells, macrophages, and dendritic cells (DCs) [7][8][9]. e class II major histocompatibility complex-(MHC-I-) restricted CTLs implicated as critical components in antitumor immunity, contributing to tumor killing [10]. And CD4 + T-helper ( ) cells are an important helper in the activation of CTLs. CTLs also can indirectly mediate tumor killing by the secretion of lymphokines such as gamma-interferon (IFN-c), lymphotoxin [11]. In addition, DCs, a professional antigen presenting cells (APCs), are central regulators in the T-cell-mediated antitumor immunity [12]. Conversely, tumor-intrinsic immune checkpoints' ligand like programmed death ligand-1 (PD-L1) can inhibit the activation of CTLs via binding to its receptor PD-1, resulting in the tumor immune evasion [8,13].
In this study, we intended to investigate the potential role and underlying mechanisms of TME in EOC. Immune infiltration analysis has indicated that the EOC patients with high-intensity infiltrating immune cells companied with better clinical outcomes, higher TMB, CSMD3, and MUC16 mutation. In addition, immune score was significantly related to antitumor immunity and the efficacy of ICBs in EOC. We emphasized that infiltrating dendritic cells plays a leading role in this antitumor immunity. Furthermore, we identified six DC-related prognostic genes for establishing the DC-related risk model using machine-learning Lasso regression, support vector machines (SVMs), and random forest. e well-established DC-related risk model in this study can accurately stratify patients into subgroups with different clinical survival.

Immune Infiltration Analysis.
e "ESTIMATE" (Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data) algorithm (R package "estimate") was carried out to calculate the immune score, which refers to the infiltration level of total immune cells [18]. e populations of major types of infiltrating immune cells were evaluated through "xCell" (R package "xCell") [19].

Somatic Mutation Analysis.
e mutation MAF (minor allele frequency) file of EOC was downloaded through R package "TCGAbiolinks". e total number of somatic mutations in each sample from TCGA was extracted through R package "maftools". In TCGA, GRCh38 is the reference genome with a length of about 35 Mb, and the formula for calculating the tumor mutation burden (TMB) is as follows [20,21]: TMB � Sn/35 (where Sn represents the total number of somatic mutations). e R package "maftools" was also used to plotting.

Gene Set Variation Analysis.
All gene sets in this study for gene set variation analysis (GSVA) were obtained from the molecular signatures database (MSigDB, https://www. gsea-msigdb.org/gsea/msigdb, Supplementary Table 1). e GSVA score of each gene set was calculated using R package "GSVA" with ssGSEA (single-sample gene set enrichment analysis) method [22,23].

Differential Expression
Analysis. EOC patients were classified into high-DC and low-DC infiltration groups based on the median level of infiltrating DCs. R package "Limma" was conducted for differential expression analysis. Absolute values of log 2 fold change (log FC) > 1 and adjusted p (adj.p) value < 0.05 were used as thresholds to identify differentially expressed genes (DEGs) [24].

Statistics.
All statistics were performed in R program (version 4.0.0). Student's t-test or Wilcox test (according to the Shapiro test and Bartlett test) was applied to calculate the p value between two groups with continuous variables. Kaplan-Meier analysis was used for survival analysis, and the statistical significance was tested using the log-rank test. e correlation between two variables was assessed by the spearman correlation analysis. p < 0.05 was considered statistically significant.

High Infiltration of Immune Cells Predicts Favorable
Clinical Outcomes in Epithelial Ovarian Cancer. To better understand the role of the infiltrating immune cells in EOC progression, the immune score based on the gene expression profiles was calculated using the "ESTIMATE" algorithm. As shown in results from GSE105437 and GSE4122, EOC patients had significantly higher immune score than normal tissues (Figure 1(a) and 1(b)). is result was consistent with our previous findings that immune cells have been increasingly infiltrated into tumor microenvironment [26]. In GSE32062, the alive patients had higher immune score than the dead (Figure 1(c), p < 0.001). In contrast, patients with high grade (Grades 3 and 4) demonstrated significantly higher immune score (Figure 1(d), p < 0.001), while the immune score in patients with FIGO stage III and stage IV did not differ significantly (Figure 1(e)).
Regarding the most common mutation like TP53 and BRCA1/2 in EOC [27], we plotted the distribution of the immune score across other frequently mutated gene TTN, CSMD3, and MUC16 in TCGA mutation data (Figure 1(f )). Patients with high immune score had higher TMB, CSMD3, and MUC16 mutation than in wild type (Figure 1 To evaluate the potential prognostic utility of the infiltrating immune cells in EOC, the patients were initially classified into high and low groups according to the median immune score. In GSE32062 cohort, the median OS of the high and low groups were 69 and 50 months (Figure 2(a), p � 0.01). And the immune score also was an indicator of longer progression-free survival (PFS) (Figure 2(b)). Moreover, GSE63885 also confirmed that the high immune score was correlated with better OS (Figure 2(c)). e prognostic accuracy of immune score to predict the OS in GSE32062 and GSE63885 were assessed via receiver operating characteristic (ROC) curve curves, as shown in Supplementary Figure 2, e area under the curve (AUC) value indicated that the immune score can be used to predict OS in EOC patients (Supplementary Figure 2).

EOC Patients with High Infiltrating Immune Cells Have an Antitumor Phenotype.
e dynamic counterbalance between immunity and evasion in TME contributes to diverse immune phenotype, which is affected by lots of factors, including the infiltrating immune cells, cytokines, and immune checkpoints [28]. "xCell" algorithm was applied to quantify the different types of infiltrating immune cells. We observed high density of antitumor immune cells in highimmune-score group, including CD8 + T cells, B cells, DCs, and macrophages M1 (Figure 2(d), all p < 0.001). Conversely, immunosuppressive cells Tregs and macrophages M2 were more accumulated in the low group ( Figure 2(e), all p < 0.01). Besides, 16 immunostimulatory-related cytokines were overexpressed in patients with high infiltrating immune cells, which include chemokines and receptors (CXCL10, CCL11, CXCL13, CXCL9, CXCL11, CXCR3, CCL5, CCL4, CCR1, and CCL8), interferons and receptors (IL2RB, IL32, IL2RG, and IL10RA), and other cytokines (IDO1 and CSF1) ( Figure 2(f ), |log FC| > 1, p < 0.05). Moreover, the GSVA scores of innate immunity and adaptive immunity were significantly correlated with infiltrating immune cells (Figure 2(g), all p < 0.001). We defined that EOC patients with high infiltrating immune cells trend to form an antitumor phenotype with more antitumor immune cells and immunostimulatory-related cytokines.
We next sought to investigate the correlation between infiltrating immune cells and immune checkpoints TIGIT, CD48, PDCD1 (PD-1), LAG3, CD274, HAVCR2, LAIR1, and PDCD1LG2(PD-L2). e result elucidated that infiltrating immune cells were positively correlated with all above immune checkpoints ( Figure 2(h), cor > 0.6). Notably, the correlation coefficients with CD48, HAVCR2, and PDCD1LG2 were more than 0.8, suggesting that EOC patients with high infiltrating immune cells may be more sensitive to ICBs therapy.

Dendritic Cells Shape the Immune Phenotype of EOC Microenvironment.
e correlation between infiltrating immune cells and different types of immune cell was assessed by Spearman's correlation analysis, and the results showed that the DCs were strongly positively correlated with immune score (Figure 3(a), cor � 0.896, p < 0.001). Meanwhile, the enrichment levels of DC-related biological processes were high in patients with high immune score ( Figure 3(b)). It is worth mentioning that cytokines CCL4, CXCL9, and CXCL10, which were highly expressed in the high infiltrating immune cells group, can attract DCs (Figure 2(f )). In addition, K-M survival analysis indicated that the patients with high infiltrating DCs presented significantly better clinical survival (Figure 3(c)).

Construction and Validation of Dendritic Cell-Related
Risk Model. As shown in volcano plot, a total of 120 genes (including 119 upregulated and 1 downregulated genes) were found to be differently expressed in high infiltrating DC group (Figure 4(a), |log FC| > 1, p < 0.05). en, 15 out of 120 DEGs were screened as the OS-related DEGs via univariate Cox analysis (Figure 4(b)). Afterward, the machinelearning Lasso regression, SVM, and random forest were employed to identify the important factors (Figure 4(c)-4(f )). Collectively, six DC-related prognostic genes CXCL9, VSIG4, ALOX5AP, TGFBI, UBD, and CXCL11 were identified to establish a risk model for EOC (Figure 4(g)). ree of them (CXCL9, UBD, and CXCL11) were protective genes with hazard ratio (HR) < 1, and others (VSIG4, ALOX5AP, and TGFBI) were risky genes (Figure 4(b)). e DC-related risk model had been applied to TCGA cohort, which classified EOC patients into high-and lowrisk groups based on median risk score ( Figure 5(a)). Different OS was noticed between two subgroups ( Figure 5 We also confirmed that the risk score was significantly negatively correlated with infiltrating immune cells and DCs but not related to tumor purity ( Figure 5(d)-5(f )).

Nomogram.
To analyze the relationship between the DC-related risk model and clinical parameters, the nomogram was developed based on clinical parameters (FIGO stage and Age) integrated with the risk score. e result illustrated that the risk score we constructed was an independent predictor in EOC (Figure 6(a)). Calibration curve showed that the nomogram has a good reliability in predicting 3-and 5-year survival (Figure 6(b)).

Journal of Oncology
In this study, we found that the immune cells were highly infiltrated into TME in carcinogenesis and predict a better clinical outcome in EOC, which are consistent with the study by Hao et al. [32], that the immune score can be used as a powerful predictive tool for both prognosis and chemotherapeutic sensitivity of EOC.
As we all know, TME phenotype was determined by a complex regulatory network involving antitumor immune responses, tumor escape, and cellular and molecular characteristics. In the current work, an antitumor phenotype with more antitumor immune cells and immunostimulatoryrelated cytokines was observed in high infiltrating lymphocytes group. In addition, consistent with Nishino et al. [4,33], we also observed a strong positive correlation between infiltrating immune cells and immune checkpoints, such as CD48, HAVCR2 and PDCD1LG2, which is an important intrinsic immune surveillance escape mechanism; the PD-1, CTLA-4, BTLA, HAVCR2, and Lag3 have been proved to be the potent   Journal of Oncology immune checkpoints to mediate tumor immunosuppression [34,35]. Among which, predominately centered on PD-L1 that overexpressed in tumor cells, leading to inhibition of CTLs' proliferation [36]. Single-agent PD-1 blockade and PD-1/PD-L1 dual-regimen therapy both exhibit longer PFS in persistent and recurrent EOC [37]. HAVCR2 is significantly overexpressed in CD4 ＋ and CD8 ＋ T cells and suppress the antitumor immune responses in primary ovarian cancer, the blockade of HAVCR2 leads to sustained antitumor reactions [38]. is study suggests that EOC patients with high infiltrating immune cells may be able to achieve better efficacy in ICB's therapy.
In addition, we observed that DCs might highly express MHC-I and MHC-II molecules, then initiating MHC-Irestricted CTLs responses and MHC-II-restricted CD4 + 1 responses, which enables T-cell immune response [39]. All in all, DCs may play a dominant role in promoting the infiltration of immune cells into TME and forming an antitumor immune phenotype. e MHC-I/II CTLs and 1  Name  IL2RG  CD2  CD3E  CXCL9  CD3D  SLAMF7  VSIG4  GZMB  ALOX5AP  CXCL13  TGFBI  UBD  CXCL10  immune responds were activated by DCs, both can further enhance the antitumor immunity to kill the tumor cells. Immunological epigenetic alterations have already become innovative and precise cancer biomarkers in urologic, brain, lung, breast, and colorectal cancers [40][41][42]. However, in EOC with highly heterogeneous, the traditional tumor, node, and metastasis (TNM) classification are not efficient for prognostic assessment [43]. In this study, six DC-related prognostic biomarkers (CXCL9, VSIG4, ALOX5AP, TGFBI, UBD, and CXCL11) were identified to construct risk model, which could accurately stratify the EOC patients into two subtypes with different survival outcomes, providing an informative prognostic assessment for EOC patients. We hypothesize these molecules may play important regulatory role in EOC. CXCL9 and CXCL11 are known for their tumor suppressive properties and are expressed on the DCs; CXCL9 and CXCL11 have been described to enhance antitumor immunity by activating 1 [46]. UBD (ubiquitin D) also known as FAT10, is induced by mature DCs and contributes to antigen presentation [47]. As for risky genes, VSIG4 (V-set and Ig containing 4), a transmembrane receptor, was expressed specifically in DCs and macrophages, which inhibited the T cells' proliferation and immune response [48,49]. e activation of ALOX5AP (Arachidonate 5-lipoxygenase activating protein) correlates  with the HER2 (human epidermal growth factor receptor 2) and promotes the growth and migration of breast cancer [50]. TGFBI (transforming growth factor beta induced) has been reported to be an oncogene in a variety of cancers including prostate, ovarian, and breast cancer [51][52][53]. Fico and Santamaria-Martínez found that TGFBI induced breast cancer metastasis by regulating the TME and hypoxia [51]. e combination of interstitial TGFBI and intratumoral CD8 + T cells can be used as a predictor for PD-1/L1 blockade nivolumab [54].

Conclusion
In summary, our results indicate that the heavily infiltrated immune cells were positively related to a better outcome and antitumor phenotype in EOC. ICBs therapy should be considered for EOC patients with high infiltrating immune cells. In addition, we provided the mechanistic insights in this study, showing that the DCs may play an important role in the initiation of antitumor immune. Finally, we verified that the established DC-related risk model is able to accurately stratify   patients into subgroups with different survival outcomes. Nevertheless, further validation studies of these molecular mechanisms are still warranted in the future.

Data Availability
Previously reported sequencing data were used to support this study and are available at TCGA (https://www.cancer. gov/) and GEO (http://www.ncbi.nlm.nih.gov/geo). ese prior studies and datasets are cited at relevant places within the text as references [14][15][16].

Conflicts of Interest
e authors declare that they have no conflicts of interest.

Authors' Contributions
Shi-yi Liu and Rong-hui Zhu contributed equally to this work.

Acknowledgments
is work has benefited from e Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO). e authors would like to express our gratitude to these programs for its generous sharing of large amounts of data. is work was supported by the National Natural Science Foundation of China (grant no. 81860276) and the Key Special Fund Projects of National Key R&D Program (grant no. 2018YFC1003200).

Supplementary Materials
Supplementary Figure 1: the box plot shows that the BRCA1/ 2 mutation has no significant correlation with immune score (t-test, p � 0.271). Supplementary Figure 2: time-dependent ROC curves in (left) GSE32062 and (right) GSE63885 indicating high accuracy of immune score in OS prediction. ROC: receiver operating characteristic; OS: overall survival. Supplementary Figure 3: the work flow of this study. Supplementary Table 1: gene sets for gene set variation analysis. (Supplementary Materials)