Combination of Immune-Related Genomic Alterations Reveals Immune Characterization and Prediction of Different Prognostic Risks in Ovarian Cancer

Introduction

Ovarian cancer (OV) is the second leading cause of gynecological cancer and has the highest case-fatality rate among women, with 21,750 new cases and 13,940 deaths predicted for 2020 in the United States (Siegel et al., 2020). About 90% of patients suffering from OV have epithelial OV, which means it is of epithelial origin. High-grade serous ovarian cancer is the most common histological and the most aggressive subtype (Gloss and Samimi, 2014), with almost 80% of patients diagnosed as late stage and an approximately low 5-year survival rate of 30–40% due to concealing without effective characteristics (Dao et al., 2016; Torre et al., 2018). Despite the continuous progress in its diagnosis and treatment, the low sensitivity or specificity of the common OV biomarkers used for clinical diagnosis and recurrence surveillance, as well as the standard treatment that has no advanced improvement beyond cytoreductive surgery and platinum-based combination chemotherapy in decades, still makes it a significant threat to women’s lives. Therefore, it is of importance to understand the mechanisms of OV through its development and progression.

The development of OV is complex with several histopathological types and involves multiple alterations of oncogenes and tumor suppressor genes. Great efforts have been made to identify potential genomic alterations, either individually or jointly, many of which have been validated as major risk indicators for mortality; for instance, ERBB2 had been demonstrated as a poor prognostic predictor with elevated expression, and the combined expression of MANF and DOCK11 was identified as a novel risk factor (Qiu et al., 2014; Luo et al., 2018; Liu et al., 2020a; Tang et al., 2020). However, the functions of non-cancer cells such as stromal or immune cells and non-cellular components in a tumor microenvironment (TME) and their interactions are still poorly understood even though plenty of studies and clinical trials have been conducted for the purpose of improved survival rate and reduced chemotherapy resistance. Moreover, the TME has increasingly been shown to manipulate aberrant histological and cellular functions and plays a critical role in the subsequent evolution of malignancies, more progressive and resistant to chemotherapy (Mroue and Bissell, 2013). Accumulating evidence is uncovering the crucial roles of immunity in tumor immunosurveillance (Dunn et al., 2004; Koebel et al., 2007; Finn, 2008). Other studies of the TME during tumor development reveal multi-omics prognostic biomarkers that may be used for imaging or liquid biopsy analysis, both important to select the most suitable therapies and stratification of patients, including OV (Abadjian et al., 2017; Wu et al., 2017; Willumsen et al., 2018; Guo et al., 2019; Jiang et al., 2020). However, because of heterogeneity and developing drug resistance, consistent with low mutational burden, patients with OV often show a lower response to immunotherapy (Zhu et al., 2018). The lack of successful treatment leads us to measure comprehensive genomic and epigenomic alterations that affect outcomes and constitute therapeutic targets, and thus, further research studies are still needed urgently.

In this study, we employed high-throughput gene expression profiles with complete clinical pathological information offered in public databases to identify genes and features involved in immune-related processes and the prognosis of OV. Multiple machine-learning-based methods were employed to investigate and validate relative immune components and their interactions. An immune gene-based risk score model was constructed and verified using available clinical data. In summary, our findings may provide new ideas and targets for the precious medication of OV.

Materials and Methods

Data Collection and Processing

The fragments per kilobase million (FPKM) expression profile of TCGA RNA-sequencing data (level 3) for OV and the corresponding clinical pathological parameters were downloaded from UCSC Xena genome browser¹ (Goldman et al., 2020). Also, annotation information mapping probes to gene symbols was obtained from the GENCODE database², using the version for human release 22 (Frankish et al., 2019). To normalize both gene size and library size, FPKM values were then transformed to transcripts per million (TPM) (Wagner et al., 2012). Duplicated genes with the same stable ensemble ID were merged by their average values. Clinical data with paired expression data were then abstracted and summarized by the following criteria: (i) duplicated samples with both formalin-fixed and paraffin-embedded and frozen tissues subjected to sequencing analysis were removed, retaining one; (ii) patients without well-annotated clinical information were removed; and (iii) patients with overall survival time <90 days were also removed from the current research. The whole cohort was then stratified for training and testing the risk models using the methods below.

Two additional datasets – GSE9891 and GSE14764 – were downloaded using ‘‘GEOquery’’ package from the Gene Expression Omnibus (GEO) for external validation³. Additionally, a comprehensive immune gene dataset was obtained from the ImmPort database⁴ to filter genes enrolling in immune or inflammatory response for a more specific inspection (Bhattacharya et al., 2018).

Investigation of Tumor Infiltration Lymphocyte Subpopulations, Dimensionality Reduction, and Cluster Analysis

In the current analysis, the single-sample gene set enrichment analysis (ssGSEA) algorithm was employed to comprehensively identify immune cell types that are overrepresented in the TME calculating individual enrichment score (ES) based on weighted difference of the empirical cumulative distribution for each pairing of a sample and gene set (Barbie et al., 2009). Immune marker gene panels were collected from a literature resource representative of 28 subpopulations of tumor infiltration lymphocytes (TILs) related to both innate and adaptive immunity, and these genes are expressed neither in cancer nor normal tissues (Charoentong et al., 2017). These TILs were further classified into three categories based on their functional orientations. The enrichment score was then normalized by the min–max algorithm:

$z=\frac{x-\text{min}(x)}{\max (\mathrm{x})-\text{min}(\mathrm{x})}$

Here, x is the enrichment score calculated by ssGSEA.

Unsupervised hierarchical clustering was then performed on the basis of the Euclidean distance and complete linkage provided in “dist” and “hclust” functions, respectively, dividing OV patients into “Immune High” and “Immune Low” clusters for exploration of differentially expressed transcription patterns between these two clusters. t-Distributed stochastic neighbor embedding (t-SNE) is a popular nonlinear dimensionality reduction technique achieved via t-SNE modeling the probabilities as a Gaussian distribution computing low-dimensional coordinates of high-dimensional data embedding to a dimensionally reduced space using Cauchy distribution (Student’s t-distribution with 1 degree of freedom), often called a map (Belkina et al., 2019). These two clusters further confirmed the robustness using the “Rtsne” package.

Identification of Immune-Associated Components and Tumor Purity

We obtained the three scores of Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) described before to calculate the stromal and immune scores that represent the infiltration of stroma or immune cells in tumor tissues, as well as estimate scores, from which tumor purity can be inferred (Yoshihara et al., 2013). The algorithm was implemented based on “estimate” R package, and the Wilcoxon signed-rank test was used for comparisons between the two clusters.

The cell type identification by estimating relative subsets of RNA transcripts (CIBERSORT) uses a set of reference gene expression termed leukocyte gene signature matrix (LM22) containing 547 genes to normalize gene expression profiles, and quantifies either relative or absolute cell components with a support vector machine (Newman et al., 2015). An inference of 22 types of immune cell matrix following pairwise Pearson’s correlation coefficients and root mean square errors (RMSE) and empirical p values were also obtained at 1,000 permutations. These further determined the immune heterogeneity in different immune clusters and the correlation between genomic alterations and LM22.

Exploration of Differentially Expressed Patterns and Functional Enrichment Analysis

The “limma” package was used to screen out the differentially expressed genes (DEGs) between clusters. Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were carried out using the “clusterProfiler” package to obtain classified annotations and functional enrichment of DEGs (Yu et al., 2012). Additionally, gene set enrichment analysis for biological pathways was conducted by GSEA, a Java-based software (Subramanian et al., 2005). The annotated gene sets were specified by “hallmarker gene sets” and “c7 immunologic signature gene sets” downloaded from MSigDB as input files (Liberzon et al., 2015; Godec et al., 2016). The significance for the corresponding enriched terms was statistically set as p < 0.05 adjusted by the Benjamini–Hochberg (BH) method and visualized using the “enrichplot” package.

Least Absolute Shrinkage and Selection Operator Cox Regression

Based on the intersection of DEGs and genes offered in the ImmPort database, we explored the potential interactions among immune genes indicating prognosis. Least absolute shrinkage and selection operator (LASSO) Cox regression provided in “glmnet” package was used for a linear combination, performing continuous shrinkage and also feature selection (Tibshirani, 1997). Currently, LASSO is widely used for the survival analysis of high-dimensional data (Jiang et al., 2018). Ten-fold cross-validation was used in this study to derive the best-fit lambda value while minimizing the mean cross-validated error. A LASSO Cox model was constructed through the formula:

$\mathrm{Risk}\mathrm{score}=\sum _{i=1}^n\text{coefi}*\text{expi}$

Here, n represents the number of selected prognostic genes, and coef_i is the coefficient of each non-zero gene_i, while exp_i represents the expression value of each screened gene_i contributing to the model. Each sample enrolled in the study was then calculated using the formula and grouped for subsequent analyses.

Survival Analysis and Model Judgment

Analyses were performed in each cohort independently. The patients in the two clusters and the whole TCGA cohort were then categorized into high group and low group, respectively, after calculating the optimal cutoff value provided in “surv_cutpoint” function based on the risk score model. The overall survival (OS) rates between each of the two subgroups were compared using Kaplan–Meier survival curves, and statistical significance was implemented by the log-rank test. Patients with censored values were marked as “+,” and their survival curves were also plotted. The univariate and multivariate Cox regression analyses were also implemented among clinical pathological variables, and the risk score model was also implemented using “ezcox” package (Wang et al., 2019). Finally, time-dependent receiver operating characteristic (ROC) curve analysis was performed to assess OS prediction of sensitivity and specificity. Additionally, patients’ risk scores and survival status, as well as expressed patterns of identified prognostic genes, were also explored to illustrate their distributions with different stratifications of clinical parameters.

Meanwhile, two external cohorts from GEO were employed to validate our risk score model. Samples only with the same pathological diagnoses consistent with the TCGA cohort remained, consisting of 251 and 68 samples, respectively. Besides, our risk score model was compared among several prognostic biomarkers identified previously, revealing stability and reliability in predicting OS.

Statistical Analysis

Statistical analyses were all performed on R software version 3.6.3⁵. Comparisons between two variables were performed by the Wilcoxon signed-rank sum test. For comparisons of more than two variables, the Kruskal–Wallis test was performed. Hazard ratio (HR) and 95% confidence intervals (CI) for each variable were also calculated where needed. A two-sided, p value < 0.05 adjusted by the BH or false discovery rate (FDR) method and | log2-fold change (FC)| > 1 were regarded as statistically significant thresholds.

Results

Summary of Expression and Clinical Pathological Data

After obtaining the expression of all probes as well as clinical pathological variables from the UCSC Xena database, we set up the criteria for more rigorous strength, as samples with OS < 90 days excluded meant more evidence about treatment and medication. Here, a total of 309 samples expressing 19,711 mRNAs were enrolled. An overview of the patients included in the whole TCGA cohort and of each cluster is shown in Table 1.

TABLE 1

Table 1. Clinical pathological characteristics of OV patients in TCGA (n = 309).

Construction of Immune-Related Clusters and Exploration Heterogeneity of Components

Using the ssGSEA method, we estimated 28 subpopulations of TILs including major types that participated in antitumor and promoted tumor procession closely linked with adaptive immunity and innate immunity functions or pathways, some of which are vital components of the tumor tissue. The whole TCGA cohort was split into two different clusters based on their normalized ES (NES) and unsupervised clustering analysis: “Immune High” (n = 194, 62.78%) and “Immune Low” (n = 115, 37.22%) (Figure 1A). Furthermore, we explored the expression of GZMA and PRF1, whose geometric mean value represents immune infiltration and immune cytolytic activity (Rooney et al., 2015). These two genes both showed higher expression in the “Immune High” cluster (Figures 1B,C). We also applied another unsupervised dimensionality reduction algorithm t-SNE confirming that two clusters possessed robust assignments, in accordance with former results (Supplementary Figures 1A,B). To further explore tumor purity and heterogeneity of components between two clusters, three scores according to the ESTIMATE algorithm were assessed. We found that immune scores, stromal scores, and estimate scores in the “Immune High” cluster were all significantly increased when compared with those in the “Immune Low” cluster, meaning higher infiltrations of immune and stromal cells in the “Immune High” cluster (Wilcoxon signed-rank test, p < 0.0001) (Figures 1D–F). However, tumor purity inferred from these three scores showed a significantly opposite trend between two clusters (Wilcoxon signed-rank test, p < 0.0001), indicating declined components of tumor cells comprised as integrated TME with non-tumor cells (Figure 1G). These results suggested the presence of intratumoral heterogeneity in OV, and stratification was observed even compared with different methods.

FIGURE 1

Figure 1. Construction of immune-related clusters and exploration heterogeneity of components. (A) Clusters were constructed in The Cancer Genome Atlas (TCGA) cohort determined by ssGSEA and unsupervised hierarchical clustering analysis, and heterogeneity was explored by the ESTIMATE algorithm. Violin plots showed the immune cytolytic activity between the expression of GZMA (B) and PRF1 (C) between two clusters. Distribution of immune score (D), stromal score (E), ESTIMATE score (F), and tumor purity (G) in the “Immune High” cluster and the “Immune Low” cluster indicated heterogeneity between two clusters. NES, normalized enrichment score; ****p < 0.0001.

Differentially Expressed Genes and Potential Mechanisms Underlying Immune-Related Roles

To identify DEGs between the “Immune High” cluster and the “Immune Low” cluster, differential expression analysis was conducted, revealing and identifying a total of 381 DEGs, all of which were explored for functional analyses. The top 10 of the enriched GO terms significantly shown in the bar plot indicated that most of them were associated with immunological and tumor-associated processes, such as regulation of leukocyte activation, leukocyte migration, granulocyte activation, T-cell activation, and neutrophil activation in biological process (BP); receptor regulator activity, receptor-ligand activity, cytokine activity, and cytokine receptor binding in molecular function (MF); and extracellular matrix, plasma membrane protein complex, cytoplasmic vesicle lumen, and secretory granule membrane in cellular component (CC) (Figure 2A). The top 30 KEGG pathways enriched significantly also indicated immune-related terms such as cytokine–cytokine receptor interaction, phagosome, chemokine signaling pathway, cell adhesion molecules, NF-kappa B signaling pathway, Th17 cell differentiation, Th1 and Th2 cell differentiation and, unexpectedly, coronavirus disease 2019 (COVID-19) (Figure 2B). Recent studies have indicated that female cancer patients have an increased infection risk and develop more severe forms of COVID-19, and overexpression of CTSL pivotal for COVID-19 infection is a marker of invasion and metastasis in ovarian cancer (Sui et al., 2016; Liang et al., 2020; Rugge et al., 2020; Montopoli et al., 2021). These results may suggest a variety of coping strategies during a pandemic between inflammation and tumorigenesis. We also downloaded “hallmarker gene sets” and “c7 immunologic signature gene sets” as background gene sets for GSEA. These results indicated that immune signatures, such as “EPITHELIAL MESENCHYMAL TRANSITION” (p adjusted = 0.0044) and “INFLAMMATORY RESPONSE” (p adjusted = 0.0044), were most significantly enriched in patients in the “Immune High” cluster (Figures 3A,B and Supplementary Figures 2A–F). Additionally, a significant enrichment in the “Immune Low” cluster was “WNT_BETA_CATENIN_SIGNALING” (p adjusted = 0.0167) (Figure 3C). As for immunologic gene sets, we also found a positive enrichment in the “Immune High” cluster corresponding to immune cells and other relative terms (all p adjusted < 0.05) (Figures 3D–F and Supplementary Tables 1, 2). Importantly, these intimate relationships between clusters and immune-related gene sets were confirmed without a doubt. It could be suggested that exploration of pathways and signatures aforementioned in OV development and revealing inherent molecular mechanisms involved may be urgent.

FIGURE 2

Figure 2. Functional analyses based on differentially expressed genes (DEGs). (A) Top 10 results of the GO enrichment analysis in BP, CC, and MF. The bright yellow-to-salmon pink chromatograms indicate the corresponding p values corrected by the BH method. (B) Top 30 results of the KEGG enrichment analyses. The bright yellow-to-salmon pink chromatograms indicate the corresponding p values corrected by the BH method. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; BP, biological processes; CC, cell components; MF, molecular functions.

FIGURE 3

Figure 3. Significant enrichment terms for GSEA in “hallmarker gene sets” (A–C) and in “c7 immunologic signature gene sets” (D–F). The significance of NES was all adjusted on behalf of the false discovery rate. NES, normalized enrichment score; FDR, false discovery rate.

Extraction of Differentially Expressed Immune Genes and Construction of a Prognostic Model

Focusing on functional analyses, we next intended to extract immune-related genes from DEGs, including 322 upregulated genes and 59 downregulated genes (Figure 4A and Supplementary Figure 3A). Based on the ImmPort database, 122 genes were extracted and their expression patterns were explored as differentially expressed immune genes (DEIGs), consisting of 117 upregulated and 5 downregulated genes, respectively (Figures 4B,C). To acquire genes with the greatest potential prognostic values, we used the “Immune High” cluster as the training set. LASSO Cox regression analysis was performed with 10-fold cross-validation to evaluate and eliminate variables which contributed less to the model. Finally, a total of 11 mRNAs, namely MSR1, FPR1, RNASE2, GBP2, CXCL9, CXCL11, C5AR1, CCL13, FGF17, CXCL14, and PI3, related to OS with non-zero coefficients were selected as candidate predictors contributing to a linear model, and then they were validated and assessed, and significant differences in OS were observed (Figures 4D,E; Supplementary Figure 3B; Supplementary Table 3). The expression levels and regression coefficients were integrated, and therefore, a risk score model was developed. For further survival analyses, a risk score based on the model was calculated for each sample.

FIGURE 4

Figure 4. Screening immune-related genes and construction of a prognostic model. (A) A volcano plot of all DEGs was shown combined with | log2-FC| and adjusted p value. (B) Extraction of immune-related genes from DEGs. (C) The landscape of expression pattern for all DEIGs between two immune clusters. Ten-fold cross-validation for tuning parameter selection (D) and coefficients extraction (E) in the LASSO Cox model. DEGs, differentially expressed genes; DEIGs, differentially expressed immune genes.

Stratification of Samples and Verification of Independent Prognostic Model

The whole TCGA cohort and the “Immune Low” cluster were all enrolled for validation, and these three cohorts were stratified into high-risk group and low-risk group followed by the cutoff point. The Kaplan–Meier plot showed significant differences in terms of patients’ OS between the high- and low-risk groups: high-risk group (n = 124) versus low-risk group (n = 70) in the “Immune High” cluster (log-rank test, p < 1.0E-4) (Figure 5A); high-risk group (n = 21) versus low-risk group (n = 94) in the “Immune Low” cluster (log-rank test, p < 1.0E-4) (Figure 5B); and high-risk group (n = 122) versus low-risk group (n = 187) in the whole TCGA cohort (log-rank test, p < 1.0E-4) (Figure 5C). Though no significant difference in patients’ OS between two immune clusters was observed (log-rank test, p = 0.2, Figure 5D), we could also distinguish the significance stratified by the median cutoff point of each cluster, meaning that a higher risk score may more likely belong to the “Immune High” cluster and may predict a worse prognosis (log-rank test, p = 4.6E-5, Figure 5E).

FIGURE 5

Figure 5. Assessment of candidate immune genes and prognostic capacity. A log-rank test was set up for statistics. Kaplan–Meier survival analysis based on the risk score model in the “Immune High” cluster (A), the “Immune Low” cluster (B), and the whole TCGA cohort (C). Patients grouped by best cutoff values were listed in the risk score chart below at different time points. (D) Kaplan–Meier survival analysis between two immune clusters. (E) Kaplan–Meier survival analysis based on the risk score model between two immune clusters. Patients were grouped by median value in each cluster. ROC and the corresponding AUC of the risk score model and other clinical pathological characteristics for the “Immune High” cluster (F) and the whole TCGA cohort (G). Time–ROC of the risk score model for the “Immune High” cluster (H) and the whole TCGA cohort (I). The ROC and AUC of the predictions for 1, 3, and 5 years are shown, respectively. (J) Correlation heatmap of immune genes extracted by LASSO, risk score model, and clinical pathological characteristics. The lower triangular was the correlation coefficients between two variables and the upper triangular matrix was the significant adjusted p values. RD, tumor residual disease; ROC, receiver operating characteristic; AUC, area under the curve; *p < 0.05, **p < 0.01, ***p < 0.001.

Univariate Cox regression analysis was performed for these studies. The risk score model was independently a negative prognostic factor for the training cohort (HR: 7.02, 95% CI: 4.06–12.12, p < 0.001). These results also showed that age (HR: 1.02, 95% CI: 1.01–1.04, p = 0.010), tumor residual disease (HR: 1.34, 95% CI: 1.10–1.63, p = 0.004), and stage (HR: 1.61, 95% CI: 1.05–2.46, p = 0.027) served as independent prognostic risk factors (Table 2A). Multivariate Cox regression analysis was performed using the significant prognostic factors identified in the univariate analysis. The risk score model was further indicated to possess predictive performance ability for OS, owning the most significant prediction (HR: 5.24, 95% CI: 2.92–9.43, p < 0.001). Similar results were also obtained in the whole cohort, as the risk score model performed best in both univariate analysis (HR: 3.90, 95% CI: 2.61–5.81, p < 0.001) and multivariate analysis (HR: 3.01, 95% CI: 1.98–4.57, p < 0.001, Table 2B). Time–ROC and area under the curve (AUC) of each factor and all combined were also presented, and our risk score model showed a high AUC (Figures 5F–I). Pairwise Pearson’s correlation analysis was employed, and an adjusted p value was set up for a threshold to derive the regulation of genes enrolled in the risk score model and clinical pathological characteristics (Figure 5J). Besides CXCL14 and PI3, the other nine genes showed significant positive or negative correlations, indicating co-expression patterns between them. A significant positive relationship between venous invasion and lymphatic invasion was also observed, indicating an underlying cooperation between them in the long-term survival period in patients suffering from OV.

TABLE 2

Table 2. Univariate (A) and multivariate (B) Cox regression analysis of the risk score model and other clinical variables.

Two independent cohorts – GSE9891 and GSE14764 – were then employed as external validations for the model to confirm prognostic accuracy, and each was grouped into high and low groups using the same algorithm aforementioned. Kaplan–Meier curves showed significant differences in GSE9891 (log-rank test, p = 0.0035) and GSE14764 (log-rank test, p = 0.042), similar to our previous results, and time–ROC indicated that our immune risk score model had high sensitivity and specificity to predict survival probability (Figures 6A–D). Also, we verified our model to possess stable and reliable ability compared with several biomarkers identified previously (Figure 6E). Hence, these results indicated that our model is reliable in making a precise prediction.

FIGURE 6

Figure 6. Assessment of prognostic capacity using the two external cohorts GSE9891 (A,B) and GSE14764 (C,D), respectively. (E) ROC curves show the sensitivity and specificity of the immune risk model signature and other known biomarkers in predicting the OS of patients. ROC, receiver operating characteristic; AUC, area under the curve.

Stratification and Validation Among Clinical Pathological Factors as a Predicted Indicator

We explored the distribution of risk scores and survival time in the whole TCGA cohort, annotated with the distribution of expression pattern of immune genes between two immune clusters (Figure 7A). Patients with higher risk scores tended to have shorter survival time and more likely remained in the “Immune High” cluster, whereas patients with lower risk scores tended to have longer survival time and also more likely remained in the “Immune Low” cluster. Moreover, we specifically explored the distribution of risk scores and expression levels of each gene in the whole cohort stratified by different clinical characteristics. As previously described, risk scores in the older group (median age > 58 years old) gained higher levels, accompanied by higher expression levels of CXCL14 and FGF17 and lower expression levels of GBP2 (Wilcoxon signed-rank test, p < 0.05, Figure 7B). Patients with increasing tumor residual disease diameters similarly gained higher values of risk scores (Kruskal–Wallis test, p < 0.05, Figure 7C). While the number of patients in the current study with either venous invasion or lymphatic invasion was comparatively small, we did observe some notable and significant intergroup differences among immune genes using the prognostic model, indicating its potential prognostic role in diverse pathological situations (Wilcoxon signed-rank test, p < 0.05, Figures 7D,E).

FIGURE 7

Figure 7. Risk score distribution and clinical stratification for the whole TCGA cohort. (A) The distribution of risk scores, survival status, and expression pattern of immune genes in the whole TCGA dataset. (B) OV patients were divided into two groups according to median age, and then the relationships among risk score distribution and expression levels of immune genes between age groups were explored. (C) Relationships among risk score distribution, expression levels of immune genes, and different status of tumor residual disease. (D) Relationships among risk score distribution, expression levels of immune genes, and different status of venous invasion. (E) Relationships among risk score distribution, expression levels of immune genes, and different status of lymphatic invasion. ^∗p < 0.05, ^∗∗p < 0.01, ^∗∗∗p < 0.001, ^****p < 0.0001.

Tumor Immune Landscape and Upstream Regulatory Mechanisms

Three hundred and nine samples were further quantified for a view of tumor immune landscape using the CIBERSORT method, as well as the different identification between the two clusters, which pointed out proportion changes in the immune microenvironment (Supplementary Figure 4). We inspected the association between risk score distribution and each of the 22 leukocyte cells. Pairwise Pearson’s correlation analysis indicated that the risk score model had an intimate relationship with immune component changes, suggesting diversities of systematic antitumor treatment on changes derived from cell-mediated immune response prospectively. Among neutrophils (R = 0.39), M2 macrophages (R = 0.31) were positively associated with risk score (p < 0.05), indicating a worse prognosis, while M1 macrophages (R =−0.15) and CD8 T cells (R =−0.12) presented negative associations, hence, a relative better prognosis (p < 0.05) (Figures 8A–D and Supplementary Figure 5). Additionally, we examined the upstream regulatory mechanisms of the genes contributing to the prognostic model. Transcription factor (TF) datasets were downloaded from TRRUST (version 2) (Han et al., 2018) and LncMAP databases (Li et al., 2018). Combined with the expression of all the TFs and immune-related genes, Pearson’s correlation analysis was applied and then a Sankey plot was shown after integrating a significant pairwise comparison of the TF–immune gene regulatory network (Figure 8E). A total of 29 TF–gene pairs targeting 6 immune genes significantly were extracted from 63 pairs (Supplementary Table 4). While CXCL11 targeted by JUN owned a negative regulatory mechanism (R = −0.142, adjusted p = 0.0293), all the other TFs were observed to have positive relationships while targeting immune genes, respectively.

FIGURE 8

Figure 8. (A–D) Pearson’s correlation analysis of the risk score distribution and immune cell infiltration significantly identified by CIBERSORT in OV patients. All p values were corrected by the BH method. (E) Upstream regulatory mechanisms revealed the network among TFs and immune genes using the Sankey plot. The width of each alluvial stripe in the TF column represents -log2 (p values) of pairwise Pearson’s correlation adjusted. TFs, transcription factors.

Discussion

Using TCGA dataset in an unbiased manner, we systematically evaluated the TME in OV by multiple approaches to investigate immune activity and proposed prognostic analyses based on bulk immune genes. A risk score model was then constructed and further analyzed by machine learning through different immune clusters and clinical pathological variables. Association with the model-stratified groups identified a significant correlation with OS, suggesting that it provided an independent and reliable measurement for gyneco-oncologists to ensure appropriate pre- or post-surgery treatment and chemotherapy management which is, to some extent, crucial and urgent, due to increasing incidence, lower 5-year survival rate, higher rate of recurrence, and resistance to chemotherapy in OV (Dao et al., 2016; Torre et al., 2018; Siegel et al., 2020).

The unsupervised hierarchical clustering analysis based on ssGSEA identified two clusters that presented obvious differences of immune activity and TME heterogeneity estimated by the ESTIMATE algorithm. Inverse correlations with tumor purity were seen for most expression of immune genes (Li et al., 2016; Rhee et al., 2018), as stromal cells and immune cells are the main non-tumor components co-existing with malignant cells in bulk tumor specimen, by which the tumor purity is deduced. These non-tumor cells dilute the purity of OV cells in TME. The cytolytic activity defined as the geometric mean of GZMA and PRF1 confirmed this robust stratification for the two clusters, for their ultimate effect mechanisms in the cancer immunity cycle. However, no significant difference in overall survival time was observed between these two different clusters, suggesting that genomic alterations may play dominant roles in affecting the functionality of immune cells via modification of TME and, finally, the immunotherapy response (Abdalla et al., 2014). Genomic expression data were further utilized for expanding the depth of understanding. We identified 381 DEGs, followed by GO and KEGG analyses, and also GSEA utilizing two gene sets as genetic backgrounds revealed that significantly enriched terms were observed in immune or inflammatory pathways, underlying their associations with immune activation or inhibition mechanisms in the progression of OV. For instance, it has been proven that high levels of IL-6, a major immunoregulatory cytokine, are associated with TME alterations by binding to its specific receptor IL-6R, whose increased expression, as well as its soluble spliced variant, is regulated by tumor-associated inflammation, leading to death (Rath et al., 2010; Lane et al., 2011). Epithelial–mesenchymal transition (EMT) is a well-known mechanism involved in the biological process of tumorigenesis and resistance to adjuvant therapies, such as paclitaxel-resistant OV cells (Kajiyama et al., 2007; Savagner, 2010). It shed light on blocking the EMT pathways for preventing tumor migration and invasion, at the same time remodeling to chemotherapy and immunity therapy (Rosano et al., 2011; Du et al., 2013). By activating NF-κB, OV cells overexpressed Her2/neu to induce the expression of VEGF which substantially increases vascular permeability, suggesting its involvement in the formation of malignant ascites (Hsieh et al., 2004). Previous studies have already approved the immune-related mechanisms in pan-cancer through TCGA and GEO databases (Charoentong et al., 2017) and the immunotherapeutic strategy available for multiple cancer types such as melanoma (Pandolfi et al., 2008; Achkar and Tarhini, 2017) and non-small cell lung cancer (Herbst et al., 2018), narrowly but successfully, except for OV as a result of low mutational burden or other reasons. In the current study, aiming to genotype and immunophenotype relations involved in OV, we filtered the expression of immune genes from the ImmPort database. LASSO Cox regression analysis identified MSR1, FPR1, RNASE2, GBP2, CXCL9, CXCL11, C5AR1, CCL13, FGF17, CXCL14, and PI3 as hub genes. Among these, FPR1 participates in tumorigenicity of human cervical cancer cells via activation of immune cells induced by N-formyl peptide (Cao and Zhang, 2018; Minopoli et al., 2019). CXCL9 and CXCL11 have been associated with activation of Th1 immunity within TME and a favorable response to chemotherapy and immunotherapy in melanoma (Harlin et al., 2009; Hong et al., 2011). CCL13 can be expressed by M2 macrophages with both anti-inflammatory and tissue repair functions (Grage-Griebenow et al., 2001; Murray and Wynn, 2011). It could be speculated that these biomarkers might play vital roles in the carcinogenesis of OV. Here, we chose the parameter family = “cox” as a response type to filter genes associated with OS. Moreover, considering smaller variables, we also applied other methods to narrow target genes based on AIC or multivariable Cox regression. However, these results failed to reach higher AUCs, suggesting an optimum balance employing genes. Then, a risk score model was constructed as a linear fit form and validated on testing sets as well as on two GEO datasets. Our results predominantly indicated that the risk score model could predict prognosis, as a higher score accompanies a worse prognosis. Besides, focusing on the whole TCGA cohort, a higher risk score may more likely incline to the “Immune High” cluster with a truncated survival time. With more genome-wide annotations’ acquisition, bipartite regulation, especially TF-regulating networks, is highly specific to different cell types (Neph et al., 2012). Upstream regulatory mechanisms were then explored and more molecular interactions were gained, which means a complex co-regulation network exists while considering the immune-related events in OV that need further exploration.

Univariate and multivariate Cox regression analyses were explored and identified independent clinical predictors additionally. As a result, the risk score model performed well entirely, indicating a reverse association with prognosis. Unexpectedly, age was screened out as an independent predictor in the whole TCGA cohort. Various changes occur during aging, while aging stimulates senescence in vivo (Boulestreau et al., 2020). Accompanying senescence, hematopoietic and immune health both decline significantly, contributing to an impaired immune function in the elderly (Beerman et al., 2010). However, previous studies often prove age as a non-significant predictor in OV or other cancers (Liu et al., 2020b). In the current study, we filtered out OS < 90 days for the purpose of long-term influence on OS and treatment management. Tumor residual disease (RD), that is, lesion diameters after cytoreductive surgery, is a consistently important factor across molecular subtypes and disease patterns (Aletti et al., 2007; Chi et al., 2009; Wallace et al., 2017; Wang et al., 2017). Minimizing RD remains an important goal to improve OS when feasible. A multivariate analysis controlling for age, preoperative albumin, stage, disease dissemination pattern, molecular subtypes, and RD posed that RD is the only variable independently associated with OS (Torres et al., 2018). Here, we confirmed the above statements. We observed a significant positive relationship between venous invasion and lymphatic invasion, suggesting their roles involved in immune regulation for a long-term survival incidence. Tumoral lymphovascular space invasion (LVSI), defined as the presence of tumor cells inside capillary lumens of either a lymphatic or a microvascular system within OV, has been reported as a new biomarker of progression (Matsuo et al., 2012, 2014). Moreover, we have set up and proved an enlightening and independent method for gyneco-oncologists and pathologists detecting immune-related genes to stratify prognosis of patients with OV.

Again, we employed CIBERSORT inferring quantitatively infiltrating lymphocytes from tumor transcriptomes. These results not only demonstrated diversity between immune clusters but also showed significant Pearson’s correlations between the risk score model and leukocytes. Previous literature shows that the complexity of the TME determines the functions of immune cells, especially those, such as neutrophils, with dual functions (Medina-Echeverz et al., 2014). Polarizing from M0 macrophages, the distinct immunoregulatory functions of activated M1 and M2 macrophages are antitumoral and protumoral, respectively (Noy and Pollard, 2014). Moreover, lymphocytes can strengthen cancer immune surveillance to suppress tumor cell proliferation, invasion, and metastasis (Dunn et al., 2004). The levels of tumor-infiltrated CD8+ T cells in the HGSOC tumors reveal a positive correlation with the patients’ survival regardless of the extent of residual disease, therapy, or BRCA1 mutation (Ovarian Tumor Tissue Analysis Consortium et al., 2017). Combining all the results mentioned above, we characterized the immune landscape with a risk score model and introduced a novel biomarker to predict prognosis which can further guide treatment decisions in patients with OV.

To conclude, our study has proposed multiple methods to investigate the TME landscape. Genes and the risk score model based on immune clusters were analyzed. We found that the risk score model was significantly associated with OV prognosis. Further analyses indicated that the risk score model was independent and more sensitive and specific than other clinical characteristics. Thus, we strongly believe that the immune-related model represents an important contribution and will enhance the identification of complex mechanistic insights of heterogeneity in OV.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: https://xena.ucsc.edu/; https://www.gencodegenes.org/human/release_22.html; https://www.ncbi.nlm.nih.gov/geo/; and https://www.immport.org.

Author Contributions

GZ designed and directed all the research. XZ, SC, QG, YC, TL, and JW performed the data processing and experimental analysis. GZ and XZ drafted the manuscript. All authors reviewed and approved the final version of the manuscript.

Funding

This work was supported by the National Natural Science Foundation of China (81772780 and 81902646), Heilongjiang Provincial Natural Science Foundation (LH2020H046), Postdoctoral Science Foundation of China (2020M670922), and Postdoctoral Foundation of Heilongjiang Province (LBH-Z19077).

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.653357/full#supplementary-material

Footnotes

1. ^ https://xena.ucsc.edu/
2. ^ https://www.gencodegenes.org/human/release_22.html
3. ^ https://www.ncbi.nlm.nih.gov/geo/
4. ^ https://www.immport.org
5. ^ https://www.r-project.org/

References

Abadjian, M. Z., Edwards, W. B., and Anderson, C. J. (2017). Imaging the tumor microenvironment. Adv. Exp. Med. Biol. 1036, 229–257. doi: 10.1007/978-3-319-67577-0_15

CrossRef Full Text | Google Scholar

Abdalla, D. R., Aleixo, A. A., Murta, E. F., and Michelin, M. A. (2014). Innate immune response adaptation in mice subjected to administration of DMBA and physical activity. Oncol Lett. 7, 886–890. doi: 10.3892/ol.2013.1774

PubMed Abstract | CrossRef Full Text | Google Scholar

Achkar, T., and Tarhini, A. A. (2017). The use of immunotherapy in the treatment of melanoma. J. Hematol. Oncol. 10:88. doi: 10.1186/s13045-017-0458-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Aletti, G. D., Dowdy, S. C., Podratz, K. C., and Cliby, W. A. (2007). Relationship among surgical complexity, short-term morbidity, and overall survival in primary surgery for advanced ovarian cancer. Am. J. Obstet. Gynecol. 197:676.e1-7. doi: 10.1016/j.ajog.2007.10.495

PubMed Abstract | CrossRef Full Text | Google Scholar

Barbie, D. A., Tamayo, P., Boehm, J. S., Kim, S. Y., Moody, S. E., Dunn, I. F., et al. (2009). Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature 462, 108–112. doi: 10.1038/nature08460

PubMed Abstract | CrossRef Full Text | Google Scholar

Beerman, I., Maloney, W. J., Weissmann, I. L., and Rossi, D. J. (2010). Stem cells and the aging hematopoietic system. Curr. Opin. Immunol. 22, 500–506. doi: 10.1016/j.coi.2010.06.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Belkina, A. C., Ciccolella, C. O., Anno, R., Halpert, R., Spidlen, J., and Snyder-Cappione, J. E. (2019). Automated optimized parameters for T-distributed stochastic neighbor embedding improve visualization and analysis of large datasets. Nat. Commun. 10:5415. doi: 10.1038/s41467-019-13055-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhattacharya, S., Dunn, P., Thomas, C. G., Smith, B., Schaefer, H., Chen, J., et al. (2018). ImmPort, toward repurposing of open access immunological assay data for translational and clinical research. Sci. Data 5:180015. doi: 10.1038/sdata.2018.15

PubMed Abstract | CrossRef Full Text | Google Scholar

Boulestreau, J., Maumus, M., Rozier, P., Jorgensen, C., and Noel, D. (2020). Mesenchymal stem cell derived extracellular vesicles in aging. Front. Cell Dev. Biol. 8:107. doi: 10.3389/fcell.2020.00107

PubMed Abstract | CrossRef Full Text | Google Scholar

Cao, G., and Zhang, Z. (2018). FPR1 mediates the tumorigenicity of human cervical cancer cells. Cancer Manag. Res. 10, 5855–5865. doi: 10.2147/CMAR.S182795

PubMed Abstract | CrossRef Full Text | Google Scholar

Charoentong, P., Finotello, F., Angelova, M., Mayer, C., Efremova, M., Rieder, D., et al. (2017). Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep. 18, 248–262. doi: 10.1016/j.celrep.2016.12.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Chi, D. S., Eisenhauer, E. L., Zivanovic, O., Sonoda, Y., Abu-Rustum, N. R., Levine, D. A., et al. (2009). Improved progression-free and overall survival in advanced ovarian cancer as a result of a change in surgical paradigm. Gynecol. Oncol. 114, 26–31. doi: 10.1016/j.ygyno.2009.03.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Dao, F., Schlappe, B. A., Tseng, J., Lester, J., Nick, A. M., Lutgendorf, S. K., et al. (2016). Characteristics of 10-year survivors of high-grade serous ovarian carcinoma. Gynecol. Oncol. 141, 260–263. doi: 10.1016/j.ygyno.2016.03.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, F., Wu, X., Liu, Y., Wang, T., Qi, X., Mao, Y., et al. (2013). Acquisition of paclitaxel resistance via PI3Kdependent epithelialmesenchymal transition in A2780 human ovarian cancer cells. Oncol. Rep. 30, 1113–1118. doi: 10.3892/or.2013.2567

PubMed Abstract | CrossRef Full Text | Google Scholar

Dunn, G. P., Old, L. J., and Schreiber, R. D. (2004). The immunobiology of cancer immunosurveillance and immunoediting. Immunity 21, 137–148. doi: 10.1016/j.immuni.2004.07.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Finn, O. J. (2008). Cancer immunology. N. Engl. J. Med. 358, 2704–2715. doi: 10.1056/NEJMra072739

PubMed Abstract | CrossRef Full Text | Google Scholar

Frankish, A., Diekhans, M., Ferreira, A. M., Johnson, R., Jungreis, I., Loveland, J., et al. (2019). GENCODE reference annotation for the human and mouse genomes. Nucleic Acids Res. 47, D766–D773. doi: 10.1093/nar/gky955

PubMed Abstract | CrossRef Full Text | Google Scholar

Gloss, B. S., and Samimi, G. (2014). Epigenetic biomarkers in epithelial ovarian cancer. Cancer Lett. 342, 257–263. doi: 10.1016/j.canlet.2011.12.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Godec, J., Tan, Y., Liberzon, A., Tamayo, P., Bhattacharya, S., Butte, A. J., et al. (2016). Compendium of immune signatures identifies conserved and species-specific biology in response to inflammation. Immunity 44, 194–206. doi: 10.1016/j.immuni.2015.12.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldman, M. J., Craft, B., Hastie, M., Repecka, K., McDade, F., Kamath, A., et al. (2020). Visualizing and interpreting cancer genomics data via the Xena platform. Nat. Biotechnol. 38, 675–678. doi: 10.1038/s41587-020-0546-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Grage-Griebenow, E., Flad, H. D., and Ernst, M. (2001). Heterogeneity of human peripheral blood monocyte subsets. J. Leukoc. Biol. 69, 11–20.

Google Scholar

Guo, Q., He, Y., Sun, L., Kong, C., Cheng, Y., Wang, P., et al. (2019). Identification of potential prognostic TF-associated lncRNAs for predicting survival in ovarian cancer. J. Cell. Mol. Med. 23, 1840–1851. doi: 10.1111/jcmm.14084

PubMed Abstract | CrossRef Full Text | Google Scholar

Han, H., Cho, J. W., Lee, S., Yun, A., Kim, H., Bae, D., et al. (2018). TRRUST v2: an expanded reference database of human and mouse transcriptional regulatory interactions. Nucleic Acids Res. 46, D380–D386. doi: 10.1093/nar/gkx1013

PubMed Abstract | CrossRef Full Text | Google Scholar

Harlin, H., Meng, Y., Peterson, A. C., Zha, Y., Tretiakova, M., Slingluff, C., et al. (2009). Chemokine expression in melanoma metastases associated with CD8+ T-cell recruitment. Cancer Res. 69, 3077–3085. doi: 10.1158/0008-5472.CAN-08-2281

PubMed Abstract | CrossRef Full Text | Google Scholar

Herbst, R. S., Morgensztern, D., and Boshoff, C. (2018). The biology and management of non-small cell lung cancer. Nature 553, 446–454. doi: 10.1038/nature25183

PubMed Abstract | CrossRef Full Text | Google Scholar

Hong, M., Puaux, A. L., Huang, C., Loumagne, L., Tow, C., Mackay, C., et al. (2011). Chemotherapy induces intratumoral expression of chemokines in cutaneous melanoma, favoring T-cell infiltration and tumor control. Cancer Res. 71, 6997–7009. doi: 10.1158/0008-5472.CAN-11-1466

PubMed Abstract | CrossRef Full Text | Google Scholar

Hsieh, C. Y., Chen, C. A., Chou, C. H., Lai, K. P., Jeng, Y. M., Kuo, M. L., et al. (2004). Overexpression of Her-2/NEU in epithelial ovarian carcinoma induces vascular endothelial growth factor C by activating NF-kappa B: implications for malignant ascites formation and tumor lymphangiogenesis. J. Biomed. Sci. 11, 249–259. doi: 10.1007/BF02256568

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, Y., Wang, C., and Zhou, S. (2020). Targeting tumor microenvironment in ovarian cancer: premise and promise. Biochim. Biophys. Acta Rev. Cancer 1873:188361. doi: 10.1016/j.bbcan.2020.188361

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, Y., Zhang, Q., Hu, Y., Li, T., Yu, J., Zhao, L., et al. (2018). ImmunoScore signature: a prognostic and predictive tool in gastric cancer. Ann. Surg. 267, 504–513. doi: 10.1097/SLA.0000000000002116

PubMed Abstract | CrossRef Full Text | Google Scholar

Kajiyama, H., Shibata, K., Terauchi, M., Yamashita, M., Ino, K., Nawa, A., et al. (2007). Chemoresistance to paclitaxel induces epithelial-mesenchymal transition and enhances metastatic potential for epithelial ovarian carcinoma cells. Int. J. Oncol. 31, 277–283.

Google Scholar

Koebel, C. M., Vermi, W., Swann, J. B., Zerafa, N., Rodig, S. J., Old, L. J., et al. (2007). Adaptive immunity maintains occult cancer in an equilibrium state. Nature 450, 903–907. doi: 10.1038/nature06309

PubMed Abstract | CrossRef Full Text | Google Scholar

Lane, D., Matte, I., Rancourt, C., and Piche, A. (2011). Prognostic significance of IL-6 and IL-8 ascites levels in ovarian cancer patients. BMC Cancer 11:210. doi: 10.1186/1471-2407-11-210

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, B., Severson, E., Pignon, J. C., Zhao, H., Li, T., Novak, J., et al. (2016). Comprehensive analyses of tumor immunity: implications for cancer immunotherapy. Genome Biol. 17:174. doi: 10.1186/s13059-016-1028-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Li, L., Wang, Z., Pan, T., Sahni, N., Jin, X., et al. (2018). LncMAP: pan-cancer atlas of long noncoding RNA-mediated transcriptional network perturbations. Nucleic Acids Res. 46, 1113–1123. doi: 10.1093/nar/gkx1311

PubMed Abstract | CrossRef Full Text | Google Scholar

Liang, W., Guan, W., Chen, R., Wang, W., Li, J., Xu, K., et al. (2020). Cancer patients in SARS-CoV-2 infection: a nationwide analysis in China. Lancet Oncol. 21, 335–337. doi: 10.1016/S1470-2045(20)30096-6

CrossRef Full Text | Google Scholar

Liberzon, A., Birger, C., Thorvaldsdottir, H., Ghandi, M., Mesirov, J. P., and Tamayo, P. (2015). The molecular signatures database (MSigDB) hallmark gene set collection. Cell Syst. 1, 417–425. doi: 10.1016/j.cels.2015.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., Meng, H., Nie, S., Sun, Y., Jiang, P., Li, S., et al. (2020a). Identification of a prognostic signature of epithelial ovarian cancer based on tumor immune microenvironment exploration. Genomics 112, 4827–4841. doi: 10.1016/j.ygeno.2020.08.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., Xu, F., Cheng, W., and Gao, L. (2020b). Identification and verification of a ten-gene signature predicting overall survival for ovarian cancer. Exp. Cell Res. 395:112235. doi: 10.1016/j.yexcr.2020.112235

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, H., Xu, X., Ye, M., Sheng, B., and Zhu, X. (2018). The prognostic value of HER2 in ovarian cancer: a meta-analysis of observational studies. PLoS One 13:e0191972. doi: 10.1371/journal.pone.0191972

PubMed Abstract | CrossRef Full Text | Google Scholar

Matsuo, K., Sheridan, T. B., Yoshino, K., Miyake, T., Hew, K. E., Im, D. D., et al. (2012). Significance of lymphovascular space invasion in epithelial ovarian cancer. Cancer Med. 1, 156–164. doi: 10.1002/cam4.31

PubMed Abstract | CrossRef Full Text | Google Scholar

Matsuo, K., Yoshino, K., Hiramatsu, K., Banzai, C., Hasegawa, K., Yasuda, M., et al. (2014). Effect of lymphovascular space invasion on survival of stage I epithelial ovarian cancer. Obstet. Gynecol. 123, 957–965. doi: 10.1097/AOG.0000000000000240

PubMed Abstract | CrossRef Full Text | Google Scholar

Medina-Echeverz, J., Aranda, F., and Berraondo, P. (2014). Myeloid-derived cells are key targets of tumor immunotherapy. Oncoimmunology. 3:e28398. doi: 10.4161/onci.28398

PubMed Abstract | CrossRef Full Text | Google Scholar

Minopoli, M., Polo, A., Ragone, C., Ingangi, V., Ciliberto, G., Pessi, A., et al. (2019). Structure-function relationship of an Urokinase Receptor-derived peptide which inhibits the Formyl Peptide Receptor type 1 activity. Sci. Rep. 9:12169. doi: 10.1038/s41598-019-47900-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Montopoli, M., Zorzi, M., Cocetta, V., Prayer-Galetti, T., Guzzinati, S., Bovo, E., et al. (2021). Clinical outcome of SARS-CoV-2 infection in breast and ovarian cancer patients underwent anti-estrogenic therapy. Ann. Oncol. S0923-7534(21)00097-1. doi: 10.1016/j.annonc.2021.01.069

PubMed Abstract | CrossRef Full Text | Google Scholar

Mroue, R., and Bissell, M. J. (2013). Three-dimensional cultures of mouse mammary epithelial cells. Methods Mol. Biol. 945, 221–250. doi: 10.1007/978-1-62703-125-7_14

CrossRef Full Text | Google Scholar

Murray, P. J., and Wynn, T. A. (2011). Protective and pathogenic functions of macrophage subsets. Nat. Rev. Immunol. 11, 723–737. doi: 10.1038/nri3073

PubMed Abstract | CrossRef Full Text | Google Scholar

Neph, S., Stergachis, A. B., Reynolds, A., Sandstrom, R., Borenstein, E., and Stamatoyannopoulos, J. A. (2012). Circuitry and dynamics of human transcription factor regulatory networks. Cell 150, 1274–1286. doi: 10.1016/j.cell.2012.04.040

PubMed Abstract | CrossRef Full Text | Google Scholar

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, 453–457. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

Noy, R., and Pollard, J. W. (2014). Tumor-associated macrophages: from mechanisms to therapy. Immunity 41, 49–61. doi: 10.1016/j.immuni.2014.06.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Ovarian Tumor Tissue Analysis Consortium, Goode, E. L., Block, M. S., Kalli, K. R., Vierkant, R. A., Chen, W., et al. (2017). Dose-response association of CD8+ tumor-infiltrating lymphocytes and survival time in high-grade serous ovarian cancer. JAMA Oncol. 3:e173290. doi: 10.1001/jamaoncol.2017.3290

PubMed Abstract | CrossRef Full Text | Google Scholar

Pandolfi, F., Cianci, R., Lolli, S., Dunn, I. S., Newton, E. E., Haggerty, T. J., et al. (2008). Strategies to overcome obstacles to successful immunotherapy of melanoma. Int. J. Immunopathol. Pharmacol. 21, 493–500. doi: 10.1177/039463200802100302

PubMed Abstract | CrossRef Full Text | Google Scholar

Qiu, J. J., Lin, Y. Y., Ye, L. C., Ding, J. X., Feng, W. W., Jin, H. Y., et al. (2014). Overexpression of long non-coding RNA HOTAIR predicts poor patient prognosis and promotes tumor metastasis in epithelial ovarian cancer. Gynecol. Oncol. 134, 121–128. doi: 10.1016/j.ygyno.2014.03.556

PubMed Abstract | CrossRef Full Text | Google Scholar

Rath, K. S., Funk, H. M., Bowling, M. C., Richards, W. E., and Drew, A. F. (2010). Expression of soluble interleukin-6 receptor in malignant ovarian tissue. Am. J. Obstet. Gynecol. 203:230.e1-8. doi: 10.1016/j.ajog.2010.03.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Rhee, J. K., Jung, Y. C., Kim, K. R., Yoo, J., Kim, J., Lee, Y. J., et al. (2018). Impact of tumor purity on immune gene expression and clustering analyses across multiple cancer types. Cancer Immunol. Res. 6, 87–97. doi: 10.1158/2326-6066.CIR-17-0201

PubMed Abstract | CrossRef Full Text | Google Scholar

Rooney, M. S., Shukla, S. A., Wu, C. J., Getz, G., and Hacohen, N. (2015). Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell 160, 48–61. doi: 10.1016/j.cell.2014.12.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosano, L., Cianfrocca, R., Spinella, F., Di Castro, V., Nicotra, M. R., Lucidi, A., et al. (2011). Acquisition of chemoresistance and EMT phenotype is linked with activation of the endothelin A receptor pathway in ovarian carcinoma cells. Clin. Cancer Res. 17, 2350–2360. doi: 10.1158/1078-0432.CCR-10-2325

PubMed Abstract | CrossRef Full Text | Google Scholar

Rugge, M., Zorzi, M., and Guzzinati, S. (2020). SARS-CoV-2 infection in the Italian Veneto region: adverse outcomes in patients with cancer. Nat. Cancer 1, 784–788. doi: 10.1038/s43018-020-0104-9

CrossRef Full Text | Google Scholar

Savagner, P. (2010). The epithelial-mesenchymal transition (EMT) phenomenon. Ann. Oncol. 21(Suppl. 7), vii89–vii92. doi: 10.1093/annonc/mdq292

PubMed Abstract | CrossRef Full Text | Google Scholar

Siegel, R. L., Miller, K. D., and Jemal, A. (2020). Cancer statistics, 2020. CA Cancer J. Clin. 70, 7–30. doi: 10.3322/caac.21590

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

Sui, H., Shi, C., Yan, Z., and Wu, M. (2016). Overexpression of Cathepsin L is associated with chemoresistance and invasion of epithelial ovarian cancer. Oncotarget 7, 45995–46001. doi: 10.18632/oncotarget.10276

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, W., Li, J., Chang, X., Jia, L., Tang, Q., Wang, Y., et al. (2020). Construction of a novel prognostic-predicting model correlated to ovarian cancer. Biosci. Rep. 40:BSR20201261. doi: 10.1042/BSR20201261

PubMed Abstract | CrossRef Full Text | Google Scholar

Tibshirani, R. (1997). The lasso method for variable selection in the Cox model. Stat. Med. 16, 385–395. doi: 10.1002/(sici)1097-0258(19970228)16:4<385::aid-sim380<3.0.co;2-3

CrossRef Full Text | Google Scholar

Torre, L. A., Trabert, B., DeSantis, C. E., Miller, K. D., Samimi, G., Runowicz, C. D., et al. (2018). Ovarian cancer statistics, 2018. CA Cancer J. Clin. 68, 284–296. doi: 10.3322/caac.21456

PubMed Abstract | CrossRef Full Text | Google Scholar

Torres, D., Wang, C., Kumar, A., Bakkum-Gamez, J. N., Weaver, A. L., McGree, M. E., et al. (2018). Factors that influence survival in high-grade serous ovarian cancer: a complex relationship between molecular subtype, disease dissemination, and operability. Gynecol. Oncol. 150, 227–232. doi: 10.1016/j.ygyno.2018.06.002

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

Wallace, S., Kumar, A., Mc Gree, M., Weaver, A., Mariani, A., Langstraat, C., et al. (2017). Efforts at maximal cytoreduction improve survival in ovarian cancer patients, even when complete gross resection is not feasible. Gynecol. Oncol. 145, 21–26. doi: 10.1016/j.ygyno.2017.01.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, C., Armasu, S. M., Kalli, K. R., Maurer, M. J., Heinzen, E. P., Keeney, G. L., et al. (2017). Pooled clustering of high-grade serous ovarian cancer gene expression leads to novel consensus subtypes associated with survival and surgical outcomes. Clin. Cancer Res. 23, 4077–4085. doi: 10.1158/1078-0432.CCR-17-0246

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, S., Zhang, J., He, Z., Wu, K., and Liu, X. S. (2019). The predictive power of tumor mutational burden in lung cancer immunotherapy response is influenced by patients’ sex. Int. J. Cancer 145, 2840–2849. doi: 10.1002/ijc.32327

PubMed Abstract | CrossRef Full Text | Google Scholar

Willumsen, N., Thomsen, L. B., Bager, C. L., Jensen, C., and Karsdal, M. A. (2018). Quantification of altered tissue turnover in a liquid biopsy: a proposed precision medicine tool to assess chronic inflammation and desmoplasia associated with a pro-cancerous niche and response to immuno-therapeutic anti-tumor modalities. Cancer Immunol. Immunother. 67, 1–12. doi: 10.1007/s00262-017-2074-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, X., Giobbie-Hurder, A., Liao, X., Connelly, C., Connolly, E. M., Li, J., et al. (2017). Angiopoietin-2 as a biomarker and target for immune checkpoint therapy. Cancer Immunol. Res. 5, 17–28. doi: 10.1158/2326-6066.CIR-16-0206

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoshihara, K., Shahmoradgoli, M., Martinez, 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

Yu, G., Wang, L. G., Han, Y., and He, Q. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16, 284–287. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, J., Armstrong, A. J., Friedlander, T. W., Kim, W., Pal, S. K., George, D. J., et al. (2018). Biomarkers of immunotherapy in urothelial and renal cell carcinoma: PD-L1, tumor mutational burden, and beyond. J. Immunother. Cancer 6:4. doi: 10.1186/s40425-018-0314-1

PubMed Abstract | CrossRef Full Text | Google Scholar