Transcriptome and clinical data for OC patients were retrieved from the TCGA data portal using the “TCGAbiolinks” R package (Colaprico et al., 2016). A total of 353 OC samples were involved in the study after the initial quality control. The FPKM expression units were converted into TPM units for further analyses. The survival data of the TCGA-OV dataset was downloaded from the UCSC Xena platform (https://xena.ucsc.edu/). The raw gene expression profiles of GSE140082 and GSE32062 were acquired from the Gene Expression Omnibus (GEO) database by using the R package “GEOquery” (Davis and Meltzer, 2007). 380 samples were involved in the GSE140082 profile and 260 samples were involved in the GSE32062 after matching clinical data, respectively. GSE115978 and IMvigor210 were immunotherapy treated cohorts. The single-cell data set GSE115978 was available through the Single Cell Portal (https://singlecell.broadinstitute.org/single_cell). The complete annotated expression and clinical data of the IMvigor210 cohort were obtained from the study of Sanjeev Mariathasan under Creative Commons 3.0 License (Mariathasan et al., 2018). A total of 348 patients after immunotherapy treatment was involved after matching clinical data. More information on these datasets we utilized could be found in the supplementary materials.

The expression profiles of the TCGA-OV and GSE140082 cohorts were integrated with the relative proportions of CD8⁺ T cell populations to construct a weighted gene co-expression network by using the R package “WGCNA” (Langfelder and Horvath, 2008). The similarity matrix was first characterized given Pearson’s correlation value and then converted into an adjacency matrix, as selected by the weighting coefficient, ß. Subsequently, the adjacency matrix was transformed into a topological overlap matrix, followed by utilizing the dynamic tree cut method to recognize various modules with a module least size cutoff of 50.

The correlations between modules and clinical information (the infiltration level of CD8⁺ T cells, OS time, and OS statue) were investigated to determine the significance of modules by Pearson correlation’s analysis. A module was considered to have a significant correlation with CD8⁺ T cells while p-value was <0.05. The module with the highest correlation coefficient was then defined as the hub module.

The TCGA-OV patients were separated into different subtypes by performing non-negative matrix factorization (NMF) with the “brunet” standard and 50 iterations (Brunet et al., 2004; Gao and Church, 2005). The number of clusters, k, was set as 2 to 10, and the average contour width of the common member matrix was determined by the “NMF” package in R software. The minimum member of each subtype was set to 10. The optimal number of clusters was determined according to cophenetic, dispersion, and silhouette coefficients. The overall survival and progress-free survival of OC patients in the different clusters were determined by Kaplan-Meier (K-M) survival analysis. P values were calculated by the Log-rank test.

Quantification of immune infiltration was performed using three different methods. The CIBERSORT algorithm was used to calculate the proportion of infiltrating immune cell subsets (Newman et al., 2015). The ratio of the immune-stromal component in the TME under different clusters of OC patients was estimated by the Estimation of STromal and Immune cells in MAlignant Tumor tissues using the Expression data (ESTIMATE) algorithm (Computing, 2018), while the populations of eight immune cells and two stromal cells were calculated using the microenvironment cell populations counter (MCPcounter) method (Becht et al., 2016). Single-sample GSEA (ssGSEA) method was used to evaluate the functions of 12 immunotherapy-related gene sets in each sample of the TCGA-OV cohort (Hänzelmann et al., 2013). Wilcox test was used to analyze the statistical significance of differences between two groups, and the Kruskal test was used to compare differences among multiple groups. Spearman’s correlation test was used to perform the correlation analysis.

Gene ontology annotation, KEGG pathway, and GSEA analyses were conducted using the Clusterprofiler R package (Yu et al., 2012). For GO and KEGG, the enrichment analyses were performed on DEGs and only terms with P-value < 0.05 were considered significantly enriched. For GSEA, c2.cp.kegg.v7.4.symbols.gmt file was downloaded from the Molecular Signature Database as the target set to identify enriched KEGG pathways in distinct OC clusters. Only |NES| > 1 and FDR <0.05 were considered statistically significant.

Somatic mutations of the TCGA-OV patient cohort were retrieved from the online website cBioPortal (http://www.cbioportal.org/). For different mutational types, FRAME_SHIFT_DEL, FRAME_SHIFT_INS, In_FRAME_DEL, In_FRAME_INS, MUSSENSE, NOTHINE, NONSTOP, Splice_SITE, and TRANSPOT_START_SITE was classified as non-synonymous variants, while silent mutations and other mutation types were classified as synonymous variants, including Intron, 3′ UTR, 5′ UTR, 3′ Flank, 5’ Flank, IGR, RNA, and Splice_Region. Subsequently, the somatic mutation data were visualized by the “maftools” R package (Mayakonda et al., 2018), which was also used to evaluate mutated genes and calculate tumor mutational burden (TMB). The proportions of mutations were compared using a one-sided z test and a two-sided chi-squared test.

Univariate COX proportional hazards regression analysis was executed using the R package “survival”. Survival analysis and survival curves were conducted and plotted using R packages “survival” and “survminer” (Therneau and Grambsch, 2000). For certain gene expression, patients with values above and below the median were classified as “high expression” and “low expression” groups, respectively. Log-rank test was used to calculate the significance of differences in overall survival. Receiver operating characteristic (ROC) curves were plotted by combining normal tissue RNA sequencing data from the Genotype-Tissue Expression (GTEx) database with the data from the TCGA-OV cohort using R package pROC (Robin et al., 2011).

Detailed statistical analyses of bioinformatics were described above. The Association of the mutation rate between two clusters was evaluated by Fisher’s exact test. One-sided z test and chi-square test were used to compare the continuous and categorical variables between two clusters. For survival analysis, the Kaplan-Meier method and Log-rank test were employed. Spearman’s correlation test was implemented in correlation analysis. P values <0.05 were considered significant (*: p < 0.05; **: p < 0.01). Most analyses above were done by packages of R software (version 3.6.3), aside from a little piece of one that was performed joining with R version 4.1.0.

To quantify the relative proportions of CD8⁺ T cell populations in human OC samples, the CIBERSORT algorithm was applied to analyze the gene expression profiles of TCGA-OV and GSE140082 cohorts to infer immune infiltration. As shown in Supplementary Figure S1, the relative abundance of 22 distinct immune cell types was estimated in two cohorts. To address the complex regulatory processes involved in CD8⁺ T cell responses, we leveraged the WGCNA approach to reveal correlated modules in two cohorts and uncover associated clinical traits. To guarantee a scale-free network, the power of ß = 6 and 3 were selected as the soft-threshold value in TCGA-OV and GSE140082 cohorts, respectively (Figure 1A, Supplementary Figures S2, S3). In the TCGA-OV cohort, 23 regulatory gene modules were identified and the blue module showed the highest positive correlation coefficient with CD8⁺ T cell (Figure 1B). Highly correlated genes within the same module may represent similar expression patterns, biological processes, or mechanisms of regulation, therefore, 2663 genes in the blue module were selected for further screening and analysis in the TCGA-OV cohort. Similarly, in the GSE140082 cohort, 13 co-expression modules were obtained and 376 genes in the green module showed the highest correlation with CD8⁺ T cells (Figure 1C). After the intersection of CD8⁺ T cell-related genes from the two cohorts, a total of 313 candidate genes were obtained, which were further considered to be highly correlated with CD8⁺ T cell function (Figure 1D).

The mRNA levels of above 313 CD8⁺ T cell-related genes from the expression matrix of 353 OC patients in the TCGA-OV cohort were extracted to perform NMF, which was used to investigate a novel CD8⁺ T cell-based molecular classification of OC. The optimal number of clusters was set to k = 2 according to the cophenetic, dispersion, rss, and silhouette analyses (Figure 2A). Thereafter, two distinctive OC clusters were determined, including 283 patients in Cluster 1 and 70 patients in Cluster 2. The consensus map showed a clear and sharp boundary when k = 2, indicating a high correlation of OC patients in each sub-consensus (Figure 2B). Consensus maps were presented at K values of 3-10 (Supplementary Figure S4). We further used PCA and t-SNE to perform dimensionality reduction and the result revealed two main clusters, demonstrating the robustness of these clusters (Supplementary Figure S5). In addition, the K-M survival curve showed that overall survival (OS) and progress free survival (PFS) rates of Cluster 1 and Cluster 2 were significantly different (p < 0.05), and the Cluster 1 group showed a worse prognosis than the Cluster 2 group (Figures 2C,D).

The intercellular interactions in the TME could induce the changes in phenotype and biological features of many types of immune cells, especially antitumor T cell responses. To clarify the difference in microenvironment composition in the two clusters, we next applied the ESTIMATE algorithm to calculate the immune score, stromal score, ESTIMATE score, and tumor purity in the TCGA-OV cohort (Supplementary Figure S6). Cluster 2 showed elevated immune score (p < 0.001), stromal score (p = 0.012), and ESTIMATE score (p < 0.001), but lower tumor purity (p < 0.001), indicating a larger amount of immune and stromal components with more infiltrating normal cells. We next used the MCPcounter tool to quantify different immune cell populations and compare the immune scores of the two clusters in the TCGA-OV cohort. As was shown in Figure 3, several types of immune cells, including B lineage (p < 0.001), CD8⁺ T cell (p < 0.001), cytotoxic lymphocytes (p < 0.001), monocytic lineage (p < 0.001), myeloid dendritic cells (p < 0.001), natural killer (NK) cells (p < 0.001), and T cells (p = 0.019), showed higher abundance in Cluster 2. Of note, elevated tumor-infiltrating NK and CD8⁺ T killer lymphocytes indicated increasing anti-tumor functions of Cluster 2, which was consistent with a better prognosis (Figures 2C,D).

The chemokine system and other immunomodulators are required in recruiting T cell infiltration into the TME. We compiled information on 148 immunomodulators from Charoentong’s study, including chemokine, MHC, receptor, immune-stimulator, and immune-inhibitor (Charoentong et al., 2017). The comprehensive analysis of immune landscapes between clusters revealed more immune-activated Cluster 2 with significantly elevated chemokines with receptors, and immune activators (Supplementary Figure S7). We thus wondered if Cluster 2 would be more sensitive to certain immunotherapeutic strategies. We further collected 12 immunotherapy-related gene sets from previously published literature and the Molecular Signatures Database (MSigDB), including APC co-inhibition, APC co-stimulation, Checkpoint, Cytolytic activity, HLA, Inflammation-promoting, MHC class I, Para-inflammation, T cell co-inhibition, T cell co-stimulation, Type I IFN Response, and Type II IFN Response. ssGSEA was then employed to score each sample to evaluate the activity of these functions (Figure 4). Not surprisingly, Cluster 2 showed significantly elevated ssGSEA scores of these 12 gene sets, further indicating the potential effect of immunotherapies.

We conducted GO term and KEGG pathway enrichment analyses to explore the biological functions and pathways related to different OC clusters. As shown in Supplementary Figures S8A,B, differentially expressed genes (DEGs) between CD8⁺ T cell-based molecular subtypes were mainly enriched in T cell activation, external side of the plasma membrane, and chemokine receptor binding in the biological process, cellular component, and molecular function categories, respectively. KEGG pathway analysis revealed enrichment in cytokine-cytokine receptor interaction, chemokine signaling pathway, natural killer cell-mediated cytotoxicity, and PD-L1 expression and PD-1 checkpoint pathway in cancer. Further subtype characterization was performed with GSEA analysis (Supplementary Table S1; Supplementary Figures S8C,D). DEGs in Cluster 1 were mainly enriched in metabolic pathways, pathways in cancer, and pathways of neurodegeneration-multiple diseases, while DEGs in Cluster 2 were mainly enriched in natural killer cell-mediated cytotoxicity, T cell receptor signaling pathway, and Th1 and Th2 cell differentiation. These results demonstrated distinct tumor microenvironmental effects between CD8⁺ T cell-based molecular subtypes, which was consistent with the previous conclusion.

We next explored differences in terms of mutational landscapes between two clusters and the mutation profiles of OC patients were summarized and visualized using the “maftools” R package. The detailed mutation information in each sample was exhibited in a waterfall plot, where different mutation types were presented with various color annotations at the bottom (Figures 5A,B). A horizontal histogram revealed the genes with higher mutation frequency in two OC clusters, respectively. The top 10 mutated signature in Cluster 1 were TP53 (94%), TTN (38%), AHNAK2 (15%), FLG2 (13%), MUC16 (12%), CSMD3 (11%), FLG (11%), HMCN1 (10%), XIST (10%), and DST (10%), while the top 10 mutated signature in Cluster 2 were TP53 (96%), TTN (43%), CSMD3 (27%), FLG2 (24%), MUC16 (18%), OBSCN (18%), AHNAK (16%), MACF1 (16%), MUC17 (16%), and PDE4DIP (16%). Considering the evidence of a correlation between increased immune cell infiltrates and BRCA1/2 mutation was found (McAlpine et al., 2012), we also compared the mutated signature of BRCA1/2 in cluster 1 (Figure 5E) and cluster 2 (Figure 5F), respectively, which presented a slight alteration of BRCA1/2 in cluster 1 (6.74%) and in cluster 2 (10.2%). To get more intuition on the mutational differences, significantly mutated genes between two clusters were presented in Figure 5C. Four genes had dramatically higher mutation frequencies in Cluster 1, including CSMD3, MACF1, PDE4DIP, and OBSCN, while SYNE2 had a higher mutation frequency in Cluster 2. Meanwhile, the presence of significantly co-occurring associations was shown in Figure 5D and Supplementary Table S2, where CSMD3 had the highest correlation with MACF1 and PDE4DIP in their mutation pattern. It has been reported that TMB is closely related to immunotherapeutic response (Snyder et al., 2014; Rizvi et al., 2015). TMB scores of different mutation types were then calculated based on TCGA-OV somatic mutation data and compared between clusters. The results showed that all mutation burdens as well as non-synonymous mutation burdens were significantly higher in Cluster 2 (Figure 5G). In summary, we believe that patients in Cluster 2 had a better effect of immune activation and were more susceptible to immunotherapy.

To screen out core essential genes for CD8⁺ T cell-related genes, univariate COX proportional hazards regression analysis and K-M survival analysis were implemented to evaluate the prognostic ability of 313 candidate genes from the WGCNA analysis (Figure 1D). A total of 35 significant prognostic genes were identified by the univariate COX regression analyses, which were then subjected to the K-M survival analyses to get the last 11 genes with survival prognosis in the TCGA-OV cohort. Next, GSE140082 and GSE32062 were used to perform external validation. After conducting univariate COX regression and K-M survival analyses on these 11 genes within two GEO datasets, CD38 and CXCL13 were revealed to have statistical differences among three cohorts (Figures 6A–F; Table 1). We went on to assay the clinical significance of CD38 and CXCL13 and plotted ROC curves. The AUC values for both genes were more than 0.8 (CD38: 0.852, CXCL13: 0.828), indicating a promising predictive performance (Figure 6G). The survival curves of the remaining 9 genes in TCGA-OV, GSE140082, and GSE32062 cohorts were shown in Supplementary Figures S9–S11. We also applied ssGSEA to assess immune infiltrates in the TCGA-OV cohort and analyze their correlation with CD38 and CXCL13 expression (Supplementary Figure S12). Not surprisingly, CD38 and CXCL13 were both highly correlated with T cells and other effector cells with tumor-killing abilities, which was consistent with their better prognosis.

We next involved two immunotherapy treated cohorts (GSE115978 single-cell cohort and IMvigor 210 cohort) to determine whether CD38 and CXCL13 could be served as emerging predictive biomarkers for immune therapy. tSNE plot of various cell types in the GSE115978 cohort were provided in Figure 7A. As shown in Figures 7B,C, CD38 was significantly elevated not only in CD8⁺ T cells, but also in macrophages, NK cells, and CD4⁺ T cells. As such, the expression of CXCL13 in various immune cell types was provided in Figures 7D,E, revealing a higher expression not only in CD8⁺ T cells but also in B cells, and cancer-associated fibroblasts (CAF). We also performed survival analysis in the IMvigor210 cohort when mRNA expression data and clinical data were combined. The results showed that there was no significant relationship between CD38 expression and the prognosis, while low CXCL13 expression had a significant correlation with poor prognosis (Figures 7F,G). In addition, no statistically significant difference in CD38 expression was observed between patients with CR/PR and SD/PD (CR, complete response; PR, partial response; SD, stable disease; PD, progressive disease), while the CR/PR group had a significantly higher expression of CXCL13 than the SD/PD group (Figures 7H,I). The results above illustrated that CXCL13 was a prognostic biomarker in OC patients with immunotherapy, showing its potential ability as an immunotherapeutic target.