Classification and clinical value of three immune subtypes of ovarian cancer based on transcriptome data

: Background : Abundant evidence suggests that tumor immune infiltration was involved in the occurrence of ovarian cancer (OvCa). Current studies have demonstrated the effect of tumor-infiltrating immune cells (TIICs) on OvCa development; few studies have found the immune genomic profile and immune subclasses of OvCa based on transcriptome data, which may help to optimally stratify patients who respond to immunotherapy. Methods : Based on the two publicly available OvCa transcriptomics data, three immunogenomic subgroups were classified using unsupervised hierarchical clustering. The GO and KEGG were analyzed in each subtype. Response to immunotherapy and anti-cancer drug targets was predicted by the TIDE, Submap algorithm, and GDSC dataset. Results: The three types of immunogenomic OvCa subsets were classified based on immune signatures. We identified three OvCa subtypes, termed hyperimmunogenic (Immunity_H), moderately immunogenic (Immunity_M), and hypoimmunogenic (Immunity_L). Each subtype has specific pathways. In the Immunity_H subtype, a number of cancer-related and immune-related pathways are overactivated. In contrast, the Immunity_L subtype is predominantly enriched in lipid metabolism. Immunity_H subtype has higher immune cell infiltration, antitumor immunoreactivity, and better survival prognosis compared to other subtypes. Predicted clinical response to immune checkpoint blockade was used by Submap and TIDE algorithm and screened potential chemotherapeutic drug targets for OvCa was employed using GDSC. After the prediction for potential drug targets, we identified several potential drug targets for the treatment of OvCa, including Parthenolide and AKT inhibitor VIII, Paclitaxel. Also, Immunity_H subgroup has an early FIGO stage, and was susceptible to respond to immunotherapy. Conclusions: The characterization of immune-based OvCa subgroups possessed potential clinical implications for OvCa treatment and has the potential to guide personalized treatment of OvCa patients through immunotherapy. (serous (SC), endometrioid (EC), clear cell (CCC), and mucinous (MC) carcinomas). OvCa can also be classified into type I (low grade) and type II (high grade) on histology, pathogenesis, tumor growth, prognosis, and response to treatment 41 . Meanwhile, Previous studies have identified four phenotype-based expression-based high-grade serous tuboovarian cancer (HGSOC) molecular subtypes based on the transcriptome data, including C1/mesenchymal (C1.MES) subtype, C2/immunoreactive (C2.IMM) subtype, C4/Differentiated (C4.DIF) subtype, and C5/Proliferative (C5.PRO) subtype 42, 43 . Currently, there is a lack of reliable molecular prognostic biomarkers for laryngeal cancer stratification based on DNA methylation profiles. Therefore, we identified subgroups of immune-related molecules based on the transcriptome data and evaluated the specific mechanisms in each subgroup in the work. In the work, the immunogenomic OvCa subsets were classified based on immune signatures using transcriptome data. Our findings show that OvCa can be divided into three subsets: Immunity_H, Immunity_M, and Immunity_L. Each subtype has specific characteristics. For example, Immunity_H had higher immune cell infiltration and stromal score, and a better prognosis, while Immunity_L contains higher tumor purity. Immunity_H has an early FIGO stage, while Immunity_L has an advanced FIGO stage. Moreover, each cluster possessed its own distinct functional enrichment terms. According to the algorithm prediction of the TIDE algorithm, the results show that the anti-PD1 immunotherapy response of Immunity_H patients is better than that Studies have also that Immunity_H can most HLA genes, indicates that it more than


Introduction
As the second most common type of cancer in women, Ovarian cancer (OvCa) has a high risk of leading mortality worldwide. OvCa has an estimated incidence of 22,530 new cases and 13,980 estimated deaths in 2019 1, 2 . Considering that the ovaries are located in the pelvic cavity, the early symptoms are not obvious, and about half of the patients are diagnosed in the late stage. Therefore, OvCa has the worst prognosis among all gynecological tumors. In most countries, the 5-year survival rate is still less than 45% 3 . Despite significant advances in surgery and chemotherapy, there is still a considerable gap to improve the prognosis. Many studies have found OvCa to be an immunogenic disease, so it was important to understand how immune cells was involved in the tumor microenvironment during tumor development [4][5][6] .
Tumor inflammation exerts a vital role in the occurrence and development of cancer. Multiple studies have demonstrated that tumor-infiltrating immune cells (TIIC) can help the host resist cancer cells and promote the development of solid tumors [7][8][9] , especially in the OvCa 10, 11 . Recently, Immunotherapy promises to be a promising option for cancer treatment, and it can also prevent drug resistance. For some types of cancer (e.g. malignant melanoma), it has yielded satisfactory results but has not worked well for OvCa 12, 13 . Therefore, It's worth considering immunotherapy in OvCa, as the treatment options for this disease are very limited. However, at present, Immunotherapy has shown beneficial effects in less than 20 percent of cancer patients. This suggests that not all OvCa patients respond to immunotherapy 14,15 . In the present work, Three immunogenomic OvCa subtypes were identified based on immune signatures in two publicly datasets. We have also identified OvCa subtypes that are sensitive to immunotherapy, which will help to explore the reasons for the poor response of OvCa to immunotherapy.

Data source and processing
RNA-seq data were acquired from two publicly available datasets: The Cancer Genome Atlas (TCGA) and the International Cancer Genome Consortium (ICGC) cohorts. The TCGA dataset included 380 primary tumor patients, ICGC dataset enrolled 85 patients. We simultaneously obtained the OvCa mutation data from the TCGA website. The clinical data was also downloaded.

Unsupervised clustering analysis
Firstly, the single-sample gene set enrichment analysis (ssGSEA) was used to quantify the enrichment level based on the 29 immune signals in each OvCa sample 16,17 . Secondly, the unsupervised hierarchical clustering among OvCa samples was performed based on the immune enrichment scores(Supplementary Excel).

Estimation of immune cell infiltration level, tumor purity, and stromal content in OvCa
The ESTIMATE method 18 was employed to estimate the level of immune cell infiltration including immune score, tumor purity, and stromal score for each OvCa sample. We used the Kruskal-Wallis method to examine the difference between OvCa subtypes.

Survival analysis
We compared the survival analysis of three OvCa subsets in the ICGC and TCGA datasets separately. Kaplan-Meier curves were drawn to exhibit differences in survival rate. Log-rank tests were employed to examine the significance of survival differences.
A P value less than 0.05 was considered statistically significant.

The proportion of immune cell subpopulations in OvCa subtypes
We used the CIBERSORT algorithm 19 for inferring the proportion of LM22 human immune cells in TCGA and ICGC. We used 1000 permutations and P < 0.05 as inclusion criteria for tumor samples. The LM22 human cells among the three subclasses were compared. The Kruskal-Wallis test was used to examine the difference among OvCa subtypes.

Gene Set Enrichment Analysis
The gene enrichment analysis (GSEA) was employed to calculate the overall enrichment score in the TCGA data set 20 . The signature gene sets curated from c2/c5 was downloaded from Broad Institute's Molecular Signature Database (MSigDB). Then we compared the KEGG pathway between Immunity_H and Immunity_L. According to a false discovery rate (FDR) <0.05, a significant enrichment pathway was identified.

Mutation analysis
We downloaded the mutation profiling in MAF of OvCa samples in the TCGA data set. We performed the mutation analysis using the R package "maftools" 21 , which calculate the mutation frequency of given genes and compare the most mutant genes between Immunity_H and Immunity_L. A P value <0.05 was considered statistically significant.

Prediction of immunotherapy response and anti-cancer drug target
We predicted the immunotherapy response using the tumor immune dysfunction and rejection (TIDE) algorithm 22 and subclass mapping 23 . We also used the public pharmacogenomics database [Genomics of Drug Sensitivity in Cancer (GDSC) https://www.cancerrxgene.org/] for predicting drug response in subtypes of samples using the R package of pRRophetic 24 . Half-peak inhibition concentration (IC50) and prediction accuracy were demonstrated by ridge regression and 10-fold cross-validation depending on the GDSC training set, respectively.

Statistical analysis
All statistical tests were conducted using R software version 4.0.0. The Wilcox test was used to compare two groups, while the Kruskal-Wallis test was used for more than two groups 25 . Kaplan-Meier curve was performed to identify prognostic related subtypes for survival data 26 . Statistically significant differences were assessed by the log-rank test using the R package "survival". All P values were set as bilateral, and a Pvalue <0.05 was regarded as statistically significant.

Identification of Immunogenomic OvCa subtypes
Firstly, a set of 29 immune-related genes that represent multiple immune cell types, functions, and pathways was downloaded. Secondly, we used the ssGSEA method to estimate the enrichment score of immune cells in tumor samples. Thirdly, immunogenomic OvCa subtypes were identified by unsupervised cluster analysis in the OvCa dataset. The two datasets exhibited similar clustering results, with three clusters separated. The three clusters: hyperimmunogenic (Immunity_H), moderately immunogenic (Immunity_M), and hypoimmunogenic (Immunity_L) were observed ( Figure 1A, B).
In all data sets, stroma scores and the immune scores in the Immunity_H subset were significantly higher, while the tumor purity in the Immunity_L subset was significantly higher (Kruskal-Wallis test, P<0.05) (Figure 2A, B). These findings indicated that the Immunity_H subtype had a higher percentage of immune and stromal cells, while the Immunity_L subtype contained a higher percentage of tumor cells.
Survival analysis results demonstrated that these OvCa subtypes have different clinical outcomes. The Immunity_H exhibited a better survival prognosis than the Immunity_M and Immunity_L in both datasets, but no statistically significant difference between the Immunity_M and Immunity_L subtypes in the ICGC dataset was detected ( Figure 3A, B).

Comparisons of HLA and immune checkpoint among three OvCa subtypes
In thymic T-cell development, T-cell responses to autoantigens are eliminated, while T-cells, for exogenous antigenic responses by their human leukocyte antigenpresenting (HLA) molecules, are preserved. The HLA locus is located on chromosome 6 and covers a 7.6 Mb region that contains more than 250 highly polymorphic genes.
The region is divided into three subregions: class I, II, and III, all of which have distinct functions. It is reported that the HLA complex is one of the important components of the immune system 27,28 . In the current study, HLA-associated genes showed significantly higher expression in Immunity_H and lower expression in Immunity_L (Kruskal-Wallis test, P < 0.05, Figure 4A (cytotoxic T cells) 29 showed the same trend, which was highest in Immunity_H and lowest in Immunity_L. These findings also proved that Immunity_H obtains a higher concentration of immune cell infiltration.
We analyzed gene expression of PD-L1 (programmed death-ligand 1), PD1 (benzoate dehydratase 1), and PD-L2 (programmed death-ligand 2) for three OvCa subtypes in two datasets. The findings demonstrated that Immunity_H exhibited the highest expression of PD-L1, PD1, and PD-L2, while the lowest expression of PD-L1, PD1, and PD-L2 was observed in the Immunity_L (Kruskal-Wallis test, P <0.05) ( Figure 5A, B). Based on the above results, we speculate that the OvCa subtype Immunity_H may respond better to PD-L1 immunotherapy because PD-1 / PD-L1 expression is usually positively correlated with more or less in vivo immune response 30 .

Identification of subtype-specific immune cell infiltration patterns among the three OvCa subsets
The CIBERSOFT method using gene expression profiling was applied to impute the proportion of LM22 infiltration in OvCa samples. The results demonstrated that significant differences between OvCa subtypes in 11 of 22 immune cells; In the Immunity_H subgroup, activated dendritic cells, macrophage M1, T-cell CD8, T-cell follicular helper cells, and T-cell CD4 memory were activated, whereas dendritic cells were activated and monocytes were significantly reduced in the subgroup Immunity_L of the TCGA dataset ( Figure 6A). In the ICGC dataset, differences between OvCa subtypes were observed for only 5 of the 22 immune cells. Dendritic cells activated, plasma cells in the Immunity_H subgroup were elevated, while in the Immunity_L subset monocyte and T-cell regulatory (Treg) was reduced ( Figure 6B). This result is lined with previous results that the Immunity_H subtype was rich in immune cells.

clinical characteristics of OvCa subtypes
In terms of clinical features, Immunity_H had a young population (Age less than 65) and early FIGO TNM stage compared to Immunity_L in the TCGA data set.
Immunity_L had advanced the FIGO TNM stage ( Figure 7A, B, C). In the ICGC data set, we found similar results that Immunity_H had a young population (Age less than 65) and early FIGO TNM stage compared to Immunity_L. During the follow-up, Immunity_H had a higher complete remission, while Immunity_L had a higher relapse.
The patients in the Immunity_L subset were mainly FIGO stage II, while patients with FIGO stage IV were enriched in the Immunity_H subset ( Figure 7D, E, F).

Identification of subtype-specific GO and KEGG Pathway
In the TCGA dataset, GSEA analysis reveals subtype-specific GO and KEGG Pathway categories between Immunity_H and Immunity_L subtype ( Figure 8A, B).
Several immune-associated GO categories were enriched in the Immunity_H subtype, including MHC class II protein complexes, regulation of natural killer cell chemotaxis, and neutrophil activation. Compared with the Immunity_H subtype, iron−sulfur cluster assembly, and intestinal lipid absorption, triglyceride lipase activity was activated in the Immunity_L subset. Similarly, the immune-associated KEGG pathways were hyperactivated in the Immunity_H subclass and included Th17 cell differentiation, B cell receptor signaling pathway, PD-1 checkpoint pathway in cancer, and B cell and T cell receptor signaling pathways. These findings reconfirmed that immunoreactivity was elevated in the subtype Immunity_H. In addition, several cancer-related pathways are overactivated in the Immunity_H subtype, such as the JAK-STAT signaling pathway, the NOD-like receptor signaling pathway, the TNF signaling pathway, and the Toll-like receptor signaling pathway. In contrast, the Immunity_L subtype is predominantly associated with enrichment in pathways related to porphyrin and chlorophyll metabolism, protein digestion, and insulin secretion.

Comparison of OvCa subtype-specific gene mutations
The OvCa subtype-specific gene mutations were explored between the Immunity_H and Immunity_L subtype. Highly mutated genetic profiles were exhibited in Figure 9A, B. The TP53, TTN, and CSMD3 were the most frequently mutated genes in the Immunity_H subtype, while TP53, TTN, and RYR2 were observed in the Immunity_L subtype. The USH2A and SYNE2 showed high mutation frequency in the Immunity_L subtype, with cutoffs less than 0.05 ( Figure 9C).

Prediction for Immunotherapy and Anti-cancer Drug Response
Immunotherapy utilizing monoclonal antibodies against the immune checkpoint inhibitor (ICPC) shows a potential application for OvCa. Immunity_H responded better to treatment with ICPC compared with Immunity_L according to the TIDE algorithm (P <0.05). The expression profiles of OvCa subtypes were compared with another published data set using the Submap algorithm for the specific immune checkpoint. The findings suggested that the patients in Immunity_H were highly sensitive to programmed cell death 1 (PD1) immunotherapy (Bonferroni correction P <0.05) ( Figure 10A).
Anti-tumor drug therapy is one of the basic treatments for patients with OvCa. The prediction model was employed to assess the response to anti-cancer drugs between two molecular subgroups based on the GDSC cell line dataset. The prediction accuracy of the approach was confirmed with 10-fold cross-validation, and the response sensitivity was estimated with IC50. Among the 138 types of anti-cancer drugs, 85 exhibited the difference between Immunity_H and Immunity_L subtype (Supplementary Table 2). The result suggested that Immunity_H was more sensitive to these kinds of drugs than Immunity_L. Parthenolide and AKT inhibitor VIII, Paclitaxel exhibit great potential treatment value for OvCa (P <0.05) ( Figure 10B). In particular, paclitaxel as the first-line chemotherapy has been widely used in clinical practice.

Discussion
Since OvCa remains the worst prognosis among gynecological tumors, it is necessary to find new therapeutic targets to improve the prognosis of OvCa and make more efforts in this field to reveal the role of immune infiltration in OvCa 11, 31, 32 .
Normally, immune cells recognize and attack cancer cells to inhibit their growth.  42,43 . Currently, there is a lack of reliable molecular prognostic biomarkers for laryngeal cancer stratification based on DNA methylation profiles. Therefore, we identified subgroups of immune-related molecules based on the transcriptome data and evaluated the specific mechanisms in each subgroup in the work. In the work, the immunogenomic OvCa subsets were classified based on immune signatures using transcriptome data. Our findings show that OvCa can be divided into three subsets: Immunity_H, Immunity_M, and Immunity_L. Each subtype has specific characteristics. For example, Immunity_H had higher immune cell infiltration and stromal score, and a better prognosis, while Immunity_L contains higher tumor purity. Immunity_H has an early FIGO stage, while Immunity_L has an advanced FIGO stage. Moreover, each cluster possessed its own distinct functional enrichment terms. According to the algorithm prediction of the TIDE algorithm, the results show that the anti-PD1 immunotherapy response of Immunity_H patients is better than that of Immunity_L patients. Studies have also confirmed that Immunity_H can more express most HLA genes, which indicates that it is more immunogenic than other subtypes. Although the specific mechanism of immune escape is unknown, this finding provides a possible target for tailored patient therapy. The above discussion shows that the heterogeneity of the tumor immune microenvironment will produce different responses to immunotherapy or anti-cancer drugs. However, bioinformatics and further clinical studies are needed to confirm these results.
In conclusion, the identification of immune-associated OvCa subtypes has potential clinical implications for OvCa therapy.          that Immune-H is more sensitive to the immunotherapy (Bonferroni-corrected P < 0.05).