Prognostic Risk Model of Immune-Related Genes In Ovarian Cancer

Aim: The prognosis of patients with ovarian cancer (OC) is highly heterogeneous which hinders to make an appropriate clinical decision. The study aimed to stratify patients’ prognoses by establishing a risk assessment model in the context of mRNA levels of immune-related genes (IRGs). Methods: Comprehensive bioinformatics analyses were done using datasets from The UCSC Xena platform, ICGC Data Portal, The Cancer Genome Atlas (TCGA), and the Genotype ‐ Tissue Expression (GTEx) project. LASSO regression was done to determine the independence of associations of specic factors with overall survival (OS). Nomogram that combined the independent prognostic factors was constructed to predict the OS of OC patients. The tumor microenvironment and immune response were estimated by cell type identication via estimating the relative subset of known RNA transcripts (CIBERSORT) and immunophenoscore (IPS). Results: Overall ten IRGs were signicantly associated with the OS probability and were used for the prognostic model construction of OC patients. According to the prognostic model, ovarian cancer samples were identied as high- or low-risk group. A nomogram containing risk score, stage, histological grade and age could signicantly predict 1-year, 3-years, and 5-years OS probability of OC patients, as well as a higher IPS score and a higher immunoreactivity phenotype, which were probably correlated with better immunotherapy response and good prognosis. Conclusions: We established a reliable IRGs-based risk model with potential prognostic value for patients with ovarian cancer. Further studies are needed to conrm these prognosis-associated biomarkers.


Introduction
As one of the most common malignancies of the female reproductive system, ovarian cancer (OC) carries the highest mortality [1]. Due to the unclear etiology and no effective early diagnosis method, 70% of OC patients are diagnosed at an advanced stage, leading to an extremely low 5-year survival rate [2]. At present, some clinical features such as residual tumor after surgery, FIGO stage and histological type are accepted to be associated with OC prognosis [3]. However, due to the high heterogeneity of ovarian cancer, it is urgent to explore more precise and practical prognostic markers and models at molecular levels for clinical OC diagnosis and therapy.
As we know, the immune system plays an important role in the development of various cancers including ovarian cancer [4]. A recent study found that immunological data (such as the type, density and location of immune cells in tumor samples) can predict the survival of OC patients [5]. Likewise, immune cells could play signi cant parts in the tumor microenvironment and affect the progression and metastasis of ovarian cancer [6]. Tumor-in ltrating macrophages and T cells in ovarian cancer were also related to the response of OC to immunotherapy and chemotherapy, and high generalized tumor in ammatory in ltrating could be an excellent marker for OC prognosis [7]. Moreover, immunotherapy such as immune checkpoints inhibitors was also proposed as a potential effective method for OC treatment [5,6]. In addition, studies have shown that immune-related genes (IRGs) in the ovary were closely related to the occurrence and development of ovarian cancer [8,9]. However, there is currently no effective prognostic model based on IRGs to assess the overall OC prognosis. Therefore, it remains an important prospect to construct an immune-based model to predict OC prognosis.
In this study, we identi ed differentially expressed immune-related genes in ovarian cancer. Then, the IRGs were integrated to construct immune-related prognostic models. In addition, we veri ed that the risk model can be used as an effective independent prognostic indicator for ovarian cancer.

Differential Gene Analysis and Functional Enrichment of IRGs
The difference analysis was carried out through the limma package of the R language, FDR<0.05, |log2FC|> 1 were gured out as the differentially expressed genes (DEGs). We downloaded 2,483 immune genes from the Immport database and screened the differential immune genes to explore the immune differences between normal tissues and tumor tissues. To explore the biological functions of differential immune genes, we used the Cluster Pro le package for Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) for analysis. In addition, ggplot2 was used to produce static visualization.

Construction and Validation of the Prognostic Risk Model
Thirty immune genes related to prognosis were selected using univariate Cox proportional hazard regression analysis and Kaplan-Meier (K-M) analysis, (p<0.05). LASSO Cox analysis was performed to determine the most important prognostic genes, then a multivariate Cox regression analysis was performed to establish a prognostic risk model for ovarian cancer patients. The formula for calculating the risk score was as follows: The K-M curve and receiver operating characteristic curve (AUC) were used to evaluate the predictive ability of the prognostic risk model. The calibration curve of the prediction model was used to evaluate the difference between the predicted probability and the actual probability. Subsequently, ICGC's gene expression pro les of ovarian cancer and squamous cervical cancer were used to verify the accuracy of the prognostic risk model.

Tumor-In ltrating Immune Cells Fraction Calculation
We calculated the tumor-in ltrating immune-cell fraction of each OC patient by the ESTIMATE CIBERSORT algorithm [10]. ESTIMATE is an algorithm that infers the ratio of stromal cells and immune cells in tumor samples based on gene expression characteristics. CIBERSORT is a universal calculation method for quantifying cell fractions from bulk tissue gene expression pro les, which was used to estimate the relative proportion of immune cells. Using 22 leukocyte expression data as a reference, the relative proportion of 22 in ltrating immune cells (including B cells) was derived through the CIBERSORT algorithm.

Statistical Analysis
ESTIMATE was calculated by the "ESTIMATE" package. CIBERSORT was implemented through code (https://rdrr.io/github/singha53/amritr/src/R/supportFunc_cibersort.R). GraphPad Prism 8.0.2 was used for paired t-test and correlation analysis. All graphs were drawn using R 3.5.3 and GraphPad Prism 8.0.2. We took two-sided P < 0.05 as the threshold of statistical signi cance.

Data source and processing
We obtained 5777 DEGs between 379 OC patients and 88 normal controls. Overall 419 immune genes were included in these DEGs (Figure 1a), including 156 down-regulated genes and 263 up-regulated genes. Then GO and KEGG pathway enrichment analyses were carried out to determine the abnormally regulated biological processes and signal pathways of the 419 differential IRGs. Speci cally, GO analysis showed that differential immune genes are involved in cytokine and growth factors related functions and pathways ( Figure 1b). KEGG pathway enrichment analysis showed that differential immune genes were enriched in "chemokine signaling pathway" and "Th17 cell differentiation" and "antigen presentation pathway" (Figure 1c). Furthermore, univariate Cox regression and Kaplan-Meier survival analysis were performed to select characteristic genes for predicting patient prognosis, and 30 differential immune genes with p<0.05 were retained for further analysis (Figure 1d). In order to avoid over tting, Lasso Cox regression was used to further screen variables (Figure 1e, f). A total of 22 genes were screened for multivariate Cox analysis, and nally 10 genes were identi ed to construct an immune prognosis model. The 10 genes were IL27RA, OBP2A, PI3, HSP90AB1, AKT2, CSPG5, CD3D, TEC, ISG20 and FGF7 ( Figure   1g). Risk score = (0.5723×expression of IL27RAl) + (0.2461×expression of OBP2A) + (0.4379×expression of PI3) + (-1.7327×expression of HSP90AB1) + (0.9677×expression of AKT2) + (-0.7023×expression of CSPG5) + (-0.6142×expression of CD3D) + (-0.5767×expression of TEC) + (-0.4948×expression of ISG20) + (0.4801×expression of FGF7). In addition, IL27RA, OBP2A, PI3, AKT2, and FGF7 were highly expressed in the high-risk group, CSPG5, CD3D, TEC and ISG20 were low expressed in the high-risk group (Pall<0.001), and HSP90AB1 had a low expression trend in the high-risk group (Figure 1h, P=0.0726).

Construction of immune-related Risk Model
In order to evaluate the ability of the 10-mRNA model to predict OC survival, we calculated the risk score of each case based on the expression of 10 mRNAs and divided the cases into high-risk groups (n = 229) and low-risk groups (n = 131) through the "maxstat" package (

Validation of the immune-related risk model
We further analyzed whether the risk score can be used as an independent predictor of OC patients. Univariate Cox regression analysis showed that age (P <0.001, HR = 1.024), stage (P = 0.040, HR = 1.371) and risk score (P <0.001, HR = 2.736) were associated with OS in OV patients (Figure 3a). Furthermore, we used the above clinical parameters and risk scores to perform a multivariate Cox regression analysis. The results showed that age (P <0.001, HR = 1.025), grade (P = 0.036, HR = 1.579) and prognostic risk model (P <0.001, HR = 2.759) can be used as independent prognostic indicators to predict the patient's OS (Figure 3b). Next, we develop a nomogram using risk score, pathological staging, histological grade, and age to predict 1-year, 3-year, and 5-year overall survival rates (Figure 3c). The consistency index (C-index) of the nomogram is 0.706, which increased the predictive ability of OS compared with the 10-mRNA classi er (C-index = 0.675). The calibration curves of 1-year, 3-year, and 5-year overall survival indicated that OS corresponded to the actual outcome closely (Figure 3d-f).

Pathway enrichment analysis of different groups
In order to study the underlying molecular mechanism of the 10-mRNA classi er in OC, we performed Gene Set Enrichment Analysis (GSEA) on the high-risk and low-risk groups. The results showed that in the high-risk group immune-related pathways such as adaptive immune response and cell differentiation of Th1 and Th2 cells were inactivated (Figure 4a), yet the MAPK pathway and Ras pathway were activated, indicating that these pathways are related to OC progression ( Figure 4b). Next, we selected the top 500 genes related to risks core characteristics for enrichment analysis and found that these genes were enriched in some important immune-in ammatory related pathways, including transendothelial migration of leukocytes pathway and TNF signaling pathway (Figure 4c).

Correlation between 10-mRNA characteristics and immune cell in ltration
Subsequently, the "ESTIMATE" algorithm was performed to calculate the stromal score, immune score, and estimate score for each OC sample. Compared with the low-risk group, the high-risk group had a higher stromal score (-343.5 vs -532.2, p=0.0146) and a lower immune score (260.9 vs 432.5, p=0.0564) (Figure 5a, b). Furthermore, correlations between 22 kinds of tumor-in ltrating immune cells' abundance were discussed in the two groups. The results exhibited higher in ltration of plasma cells, CD8+ T cells, T follicular helper cells, M1 macrophage and resting mast cells in the low-risk group, while the high-risk group was presented with higher naive B cells, CD4+ memory resting and M0 macrophage in ltration (Figure 5c).
In addition, we found that higher plasma cells, T follicular helper cells, M0 macrophage, M1 macrophage and resting mast cells in ltration tended to have a better prognosis, while low CD4+ memory resting immune cell in ltration suggested good prognosis (Figure 5d-g). As shown in Figure 5h-j, clinicopathological features including tumor grade (p=0.6259) and pathological stage (p=0.1298) showed no signi cant difference between the low-and high-risk groups, but patients with high pathological stage tend to have a higher risk score. These ndings demonstrated that the speci c changes of tumor-in ltrating immune cells might be associated with the OS of OC patients.

The potential of 10-mRNA labeling as an indicator of immunotherapy response
To explore the immunogenicity of patients, the expression of immune genes (PD1, PDL1, and CTLA4) in different groups was analyzed. As shown in Figure 6a-c, the low-risk group exhibited higher expression of PD1 (P<0.0001), PDL1 (P=0.0031) and CTLA4 (P=0.0001), suggesting higher immunogenicity of the lowrisk group. We further assessed the relation between IPS scores and risk groupings. As a result, the scores of IPS, IPS+PD1, IPS+CTLA4, and IPS+PD+1+CTLA4 increased signi cantly in the low-risk group (Figure   6d-g). What's more, we performed evaluated the correlation between mRNA characteristics and PD1, PDL1, and CTLA4. The results indicated a signi cant negative correlation between risk score and the expression of PD1 (r=-0.2461, p=0.0009), PDL1 (r=-0.2775, p=0.0002) and CTLA4 (r=-0.2534, p=0.0006) in the high-risk group, indicating that OC patients in the high-risk group were less likely would bene t from immunotherapy (Figure 6h-j).

Discussion
Ovarian cancer is the most lethal gynecologic malignancy, with a 5-year survival rate of only 20%-30% [1,2]. The low survival rate of OC is related to the low early detection rate and the lack of relevant molecular markers for effective targeted therapy [2].
Gene expression levels and mutation characteristics could provide insight into the likelihood of tumorigenesis and prognosis prediction [11,12]. Changes in mRNA expressions and mutations of OCrelated genes were likely to be crucial determinants for OC development and progression [12,13]. On the other hand, the immune system played a critical role in the development of various cancers, including ovarian cancer [4,14]. The immunological status could provide similar survival predictors compared to the histopathological features used for OC patients [5][6]. Immune cells were important components of the tumor microenvironment and involved in the development and metastasis of OC [7,14]. It suggested that the aberrantly expressed OC-related IRGs were correlated to the development of ovarian cancer and these genes might serve as biomarkers for diagnosis and prognosis.
Therefore, we conducted this study to nd a valid prognostic model based on IRGs mRNAs to assess the overall prognosis of OC patients. Our study was based on mRNA expression pro le data from TCGA and GTEx databases, a total of 10 mRNAs signi cantly associated with OC prognosis were identi ed for further analysis, including IL27RA, OBP2A, PI3, HSP90AB1, AKT2, CSPG5, CD3D, TEC, ISG20 and FGF7. Among them, IL27RA, HSP90AB1, AKT2 and FGF7 have been reported in OC development through multiple molecular pathways. Interleukin receptors IL27RA were involved in the regulation of the microenvironment of OC [15]; HSP90AB1, an ATP-dependent molecular chaperone protein, could affect the maturation, stability, and activation of several diverse client proteins, and then promoted persistent activation of various OC-related cellular kinases and transcription factors [16]; AKT2 might promote OC invasion through increased expression of β1 integrins and activation PI3K-dependent pathway [17]; FGF7 modulated growth and platinum sensitivity in ovarian cancer by decreasing cisplatin 50% inhibiting concentration (IC50) [18]. In addition to these four mRNAs, this experiment also screened out six less studied mRNAs that are also closely related to ovarian cancer prognosis, including OBP2A, PI3, CSPG5, CD3D, TEC and ISG20. The underlying mechanisms of abnormal mRNAs expression of these genes still need to be con rmed by further studies.
The study further constructed a risk model via lasso regression analysis based on the ten mRNAs. The nomogram included in the model could more accurately predict 1-year, 3-year and 5-year overall survival of OC patients than a single clinical variable, which helped take personalized anti-tumor treatment.
Further analysis for potential molecular mechanisms of mRNA showed that these immune-related gene mRNAs are involved in such as MAPK pathway and Ras pathway, which are directly involved in ovarian cancer tumor progression [19,20]; also involved in immune in ammation-related pathways in cancer [5][6], such as adaptive immune response, Th1 and Th2 cell differentiation, leukocyte transendothelial migration and TNF signaling pathway. These all suggest that tumor development is not only related to abnormal gene expression and variation, but also tumor immunity.
The tumor microenvironment plays an important role in tumor immunity. Tumor immune in ltrating lymphocytes (TILs), a component of the tumor micro-environment, have been shown to correlate with the overall survival of ovarian cancer [21]. In our study, the composition of tumor-in ltrating immune cells was estimated by ESTIMATE and CIBERSORT for each sample. Abundance of plasma cells, CD8+ T cells, T follicular helper cells, M1 macrophages and resting mast cells were found in the low-risk group and was correlated with a favorable prognosis, which was consistent with previous investigations. Meanwhile, the biomarkers of ICI were explored. Low-risk patients displayed up-regulation of PD-1, PDL1, and CTLA4. Besides, patients in the low-risk group had remarkably increased IPS, IPS+PD1, IPS+CTLA4, and IPS+PD+1+CTLA4 scores. The increased risk score and PD1, PDL1 and CTLA4 expression showed a signi cant negative correlation. Collectively, these results demonstrated that our risk model held the capacity to determine tumor response to immunotherapy in OC patients. The low-risk group had a higher immunoreactivity. The key ten mRNAs in the present study affected prognosis by interfering with immune cell in ltration in OC, as well as the presence of tumor tissue immune-related genes PD1, PDL1, and CTLA4, which were correlated with the prognosis of OS. All of these further con rming that immune factors play a crucial role in the development of ovarian cancer.
In summary, we identi ed a novel 10-immune related gene-based signature that may be applied as a prognostic predictor of ovarian cancer, thereby bene ting clinical decision-making and the personalized treatment of patients with OC. This nding provides new insights into the relationship between the immune system and OC, offering useful clinical value for OC prognosis and therapy. In the future, more data are needed to con rm this conclusion.

Declarations
Ethics approval and consent to participate The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved.

Consent for publication
Not applicable.

Competing interests
The authors (Bo Ding, Wenjing Yan, Shizhi Wang, Yang Shen) declare that they have no competing interests.    Pathway enrichment analysis of different groups (a, b) Immune-related pathways was analysed by GSEA.
(c) Enrichment analysis of IRGs to nd related pathways. Correlation between 10-mRNA characteristics and immune cell in ltration (a, b) Stromal scores were calculate with "ESTIMATE" algorithm in high-and low-risk group. (c) Tumor-in ltrating immune cells' abundance were discussed in the two groups. (d, e, f, g) Relation between different tumor-in ltrating immune cells status and the survival. (h, i , j) Relation between clinicopathological features stromal score.

Figure 6
Estimation of the Immunoreactivity (a, b, c) The expression PD1, PDL1 and CTLA4 in the high-and lowrisk group. (d, e, f, g) The scores of IPS, IPS+PD1, IPS+CTLA4, and IPS+PD+1+CTLA4 in the high-and lowrisk group. (h, i, j) The relation of risk score and the expression of PD1, PDL1 and CTLA4 in high-risk group.