Comprehensive single-cell and bulk RNA-seq analyses reveal a novel CD8⁺ T cell-associated prognostic signature in ovarian cancer

Single-cell transcriptome analysis

We obtained scRNA-seq data from seven ovarian tumors from the Gene Expression Omnibus (GEO) database (GSE184880) for this investigation [9]. The Seurat package [10] (v4.1.1) in R (v4.1.3) was used to perform unsupervised clustering of individual cells using the read count matrix as input. Rigorous quality control measures were applied, focusing on the number of detected genes and the percentage of mitochondrial gene counts per cell. Cells with fewer than 200 detected genes and those with more than 20% mitochondrial gene counts were systematically excluded from the analysis. We also omitted genes detected in fewer than 3 cells to reduce potential extraneous signals. To address batch effects, we employed the Harmony algorithm to integrate multi-sample data and correct any inconsistencies. Following the Seurat-guided tutorial, we performed dimension reduction clustering and differential expression analysis. Principal component analysis (PCA) and uniform manifold approximation and projection (UMAP) dimension reduction were conducted based on the top 20 principal components. Cell cluster annotations were determined using canonical gene markers.

Collection of public datasets

The RNA-sequencing expression matrix and clinical information of ovarian cancer samples and para-cancerous tissues were downloaded from TCGA database, available at UCSC Xena (https://xena.ucsc.edu/). Two additional independent datasets (GSE30161 and GSE63885) were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) [11, 12]. The somatic mutation data were retrieved from the Genomic Data Commons (GDC, https://portal.gdc.cancer.gov/). The somatic mutation data were then analyzed and sorted in the form of Mutation Annotation Format (MAF), followed by using the R package maftools [13] to calculate tumor mutation burden (TMB).

Construction and validation of a CD8⁺ T cell-related prognostic signature

Firstly, we collected 1,020 genes associated with CD8⁺ T cells through single-cell transcriptome analysis. We then analyzed the RNA expression matrices of breast cancer samples in the TCGA, GSE30161, and GSE63885 datasets to identify co-expressed genes. To develop a prognostic signature using CD8⁺ T cell-related genes, we performed univariate Cox regression and LASSO regression analyses. In the end, we identified 17 genes: SLFN5, CD40LG, EMB, ISG20, CD226, CCR7, S1PR4, RNF167, MLLT3, PLEKHF1, RNPEPL1, STK39, FXYD5, BTN3A3, WNK1, and SUSD3. Based on these genes, we constructed a CD8⁺ T cell-related prognostic model.

To classify ovarian cancer patients, we calculated the risk score for each patient in the training set using the following formula:

$$\text{Risk score}={\displaystyle \sum \text{ni}}={\displaystyle \sum (\text{Coefi}\ast \text{xi})}$$

Next, we divided the ovarian cancer patients into high-risk and low-risk groups based on the median of the risk score. We used the R package survivalROC [14] to estimate the predictive sensitivity of the risk score. Finally, we evaluated the effectiveness of the model in the validation set using the same coefficients and cutoff values as the training set.

Biological functional analysis between high/low-risk group patients

The Deseq2 [15] R package was used to analyze differentially expressed genes (DEGs). DEGs were identified using a cutoff of an adjusted p-value of < 0.05 and a |Log2 fold change| > 1. Gene set enrichment analysis (GSEA) was conducted with the clusterProfiler R package [16]. Gene sets with false discovery rate (FDR)-corrected p-values < 0.05 were considered significant using Fisher’s exact test. Gene set variation analysis (GSVA) was performed using the R package GSVA [17]. Gene signatures representing recurrent cancer cell states were obtained from a previous study [18].

Tumor immune microenvironment in ovarian cancer patients

We used the “estimate” algorithm to calculate immune scores for ovarian cancer samples [19]. To analyze the infiltration of immune cells, we used TIMER2.0, an algorithm that efficiently predicts immune cell infiltration based on bulk tumor gene expression data (http://timer.cistrome.org/). Using CIBERSORT, we quantified the relative abundance of 22 immune cells for each sample. Additionally, we gathered a collection of tumor immunomodulators from literature and calculated the correlation between the risk score and these modulators.

Predicting drug responses and immunotherapy sensitivity

The tumor immune exclusion score was obtained from TIDE (http://tide.dfci.harvard.edu/). The Immunophenoscore (IPS) was calculated using The Cancer Immunome Atlas (https://tcia.at/). To evaluate the predictive ability of the risk score for chemotherapeutic agents, we used the R package oncoPredict [20] to calculate patients’ half maximal inhibitory concentration (IC50) for various common chemotherapeutic agents. The Wilcoxon rank test was then employed to compare the difference in IC50 between the high-risk and low-risk groups. Additionally, we utilized the IMvigor210 dataset (specifically for uroepithelial carcinoma) to further validate the predictive value of the aforementioned immunotherapy responses.

Univariate and multivariable Cox regression

We performed univariate Cox regression analysis to investigate the correlation between gene expression and OS in patients with ovarian cancer. In the same group, we utilized multivariate Cox regression to identify independent risk factors. Genes and factors that displayed a FDR < 0.05 were deemed statistically significant in terms of patient survival. The outcomes of both univariate and multivariate Cox regression analyses were obtained and visually represented using the R package forestplot.

Establishment of the nomogram

This study utilized the Cox regression model and the R package “rms” to develop a nomogram that predicts OS at 1-, 3-, and 5-year time frames. The accuracy of the nomogram was assessed using the C-index, while calibration plots were employed to visually compare the predicted and observed OS at 1-, 3-, and 5-year intervals.

Statistical analysis

All statistical analyses were performed using R version 4.1.3 (https://www.r-project.org/) and its adequate packages. Statistical significance was set at p ≤ 0.05.

Data availability

The source of datasets supporting the conclusions are included in this article.

Interrogating the cellular constitution of ovarian cancer at single-cell resolution

In order to minutely investigate the cellular constitution of ovarian cancer at single-cell resolution and identify markers of CD8⁺ T cells, we re-analyzed the scRNA-seq data of tumors from seven treatment-naive ovarian cancer samples. Firstly, we integrated these data and corrected the potential batch effects using Harmony algorithm. After strict quality control and data filtering, data of 21785 genes within 27045 cells was obtained. We identified three major compartments within the TME of ovarian cancer, including the epithelial subset, immune subset and stromal subset (Figure 1A). UMAP visualization showed that seven scRNA-seq data were integrated and mixed uniformly (Figure 1B). By interrogating a list of classical lineage markers, we annotated the subpopulations among three main cellular subsets. For example, the stromal subset consisted of endothelial cells identified by expressions of VWF, PECAM1 and KDR, pericytes with high RGS5, ACTA2, and COL4A1 expressions, and cancer-associated fibroblasts (CAFs) identified by significant expressions of COL1A1, DCN and LUM (Figure 1C). Additionally, we identified several immune subpopulations including myeloid cells, NK cells, T cells, B cells, and plasma cells in the ovarian TME. Interestingly, a cycling T cell subpopulation which showed markers of both T cell lineage (CD2, CD3D and CD3G) and cell proliferation (HMGB2, MKI67 and TOP2A) was identified. All major cell types were represented across all seven tumors, but showed diverse composition patterns (Figure 1D). For example, the CD8⁺ T cells were enriched in OV-1, OV-3, OV-5 and OV-7 samples, while the CAFs were mainly enriched in the OV-6 sample. These observed variations in cellular composition might partly explain the inter-tumoral heterogeneity. Considering the boundary of NK cells and T cells on the UMAP plot was equivocal, we devoted to minutely discriminate distinct cellular subsets. The NK cell subpopulation was then identified based on the negative expression of canonical T cell lineage markers including CD2, CD3D, CD3G, CD4, CD8A and CD8B, along with specifically high expression of the chemokine XCL1, as well as GNLY and TRDC (Figure 1C, 1E). Furthermore, CD4⁺ T cells and CD8⁺ T cells were distinguished based on the expression of CD4, CD8A and CD8B (Figure 1C–1E). Taken together, we minutely interrogated the cellular constitution of ovarian cancer at single-cell resolution and identified cellular markers of each subpopulation, such as the CD8⁺ T cell subpopulation.

Figure 1. Interrogating the cellular constitution of ovarian cancer at single-cell resolution. (A, B) UMAP plot showing the major cell subpopulations in ovarian cancer. (C) Bubble heatmap showing expression levels of selected signature genes in ovarian cancer. Dot size indicates fraction of expressing cells, colored based on normalized expression levels. (D) Relative proportions of diverse cell types across each sample. (E) Feature plots to further identify various immune cells, based on the expression levels of marker genes.

Construction and validation of a prognostic model based on specific CD8⁺ T cell markers in ovarian cancer

Conducting differential analyses on the scRNA-seq data, we calculated a list of 1020 genes as the markers of the CD8⁺ T cell subpopulation (Supplementary Table 1). As expected, the Gene Ontology (GO) enrichment analysis exhibited that these markers were enriched in multiple T cell-associated pathways, including T cell differentiation, T cell receptor signaling pathway and alpha-beta T cell activation (Supplementary Figure 1A). Moreover, the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis showed a similar result, including enrichment of signaling pathways such as T cell receptor signaling pathway, Th1 and Th2 cell differentiation, Th17 cell differentiation as well as PD-L1 expression and PD-1 checkpoint pathway in cancer (Supplementary Figure 1A).

In order to investigate the prognostic role of CD8⁺ T cells in ovarian cancer, we firstly carried out a univariate Cox regression analysis to identify potential prognostic markers associated with CD8⁺ T cells for ovarian cancer in the TCGA cohort. A total of 64 genes was found to be correlated with the ovarian cancer prognosis (Supplementary Table 2). Subsequently, we performed the LASSO regression analysis to reduce the number of prognostic genes and finally identified 17 prognostic signature genes, including SLFN5, CD40LG, EMB, ISG20, CD226, CCR7, S1PR4, RNF167, MLLT3, PLEKHF1, RNPEPL1, STK39, FXYD5, BTN3A3, WNK1 and SUSD3 (Figure 2A, 2B). Based on the coefficients and expression levels of these prognostic signature genes, we computed the risk score for each sample in the TCGA ovarian cancer cohort, which served as our training set. Subsequently, we categorized the ovarian cancer patients in the TCGA training cohort into high-risk and low-risk groups using the median risk score as the threshold. Our analysis revealed that patients in the high-risk group exhibited markedly inferior outcomes (Figure 2C–2E). To assess the effectiveness of our CD8⁺ T cell-associated risk model, we presented receiver operator characteristic (ROC) curves and determined the area under the ROC curve (AUC) values at 1, 2, 3, and 5 years, which were 0.63, 0.67, 0.68, and 0.77, respectively (Figure 2F).

Figure 2. Construction and validation of a prognostic model based on specific CD8⁺ T cell markers in ovarian cancer. (A) Coefficient profiles in the LASSO regression model. (B) Cross-validation for tuning parameter selection in the LASSO regression. (C) Kaplan-Meier survival analysis was performed on the relationship between the risk score and OS using the TCGA training cohort. (D) The rank of risk score in the TCGA training cohort. (E) Survival status in the TCGA training cohort. (F) Time-dependent ROC curve analysis of the prognostic model (1, 2, 3, and 5 years) in the TCGA training cohort. (G) Kaplan-Meier survival analysis was performed on the relationship between the risk score and OS using the GSE30161 validation cohort. (H) The rank of risk score in the GSE30161 validation cohort. (I) Survival status in the GSE30161 validation cohort. (J) Time-dependent ROC curve analysis of the prognostic model (1, 2, 3, and 5 years) in the GSE30161 validation cohort.

To validate the reliability of the CD8⁺ T cell-related prognostic signature, we conducted additional validation tests in several independent datasets. Similarly, in the GSE30161 cohort, we stratified ovarian cancer patients based on the same risk score, consistent with the approach used in the TCGA training set and observed that patients in the high-risk group had poorer prognoses (Figure 2G–2I). The AUC values for the risk score in the GSE30161 dataset were 0.81 for 1 year, 0.73 for 2 years, 0.76 for 3 years, and 0.66 for 5 years (Figure 2J). Additionally, ovarian cancer patients with higher risk score in the GSE63885 dataset also exhibited significantly inferior prognoses (Supplementary Figure 2A–2C). The AUC values for the CD8⁺ T cell-related risk score in the GSE63885 dataset were 0.81 for 1 year, 0.73 for 2 years, 0.76 for 3 years, and 0.66 for 5 years (Supplementary Figure 2D). In summary, we have developed and validated a novel CD8⁺ T cell-related prognostic signature for predicting ovarian cancer outcomes.

Functional and genomic features of high/low-risk groups

In our quest to elucidate the underlying mechanisms explaining the prognostic significance of the CD8⁺ T cell-related risk score, we initially set out to explore its functional and genomic attributes. GO enrichment analysis revealed that the high-risk group patients showed significant enrichment in multiple immune-related pathways, such as the B cell receptor signaling pathway, T cell receptor complex and antigen binding pathways (Figure 3A). In contrast, the intermediate filament cytoskeleton organization, catenin complex and G protein-coupled peptide receptor activity pathways were found to be enriched in the low-risk group patients (Figure 3B). In addition, the results of GSEA indicated that patients in the high-risk group had higher enrichment scores in pathways associated with epithelial-mesenchymal transition (EMT), Hedgehog signaling, and estrogen response (Figure 3C and Supplementary Figure 3A–3C), while those in the low-risk group showed elevated enrichment scores in pathways related to interferon-gamma response, DNA repair, and MYC targets V1 (Figure 3C and Supplementary Figure 3D–3F). In order to obtain further insights into the transcriptional heterogeneity among ovarian cancers, we conducted the GSVA algorithm to calculate 16 recurrent cancer cell states, which interacting with the TME to take shape organized systems qualified to promote immune escape, metastasis and drug resistance [18, 21]. We found the high-risk group patients harbored higher signature scores of astrocyte (AC)-like, alveolar and partial epithelial-mesenchymal transition (pEMT) modules (Figure 3D). Nevertheless, the cell cycle, interferon and oxidative phosphorylation modules were enriched among low-risk group ovarian cancer patients (Figure 3D).

Figure 3. Functional and genomic features of high/low-risk group patients. (A, B) GO enrichment of DEGs in high/low-risk groups. (C) Bar plot showing different pathways enriched in high/low-risk groups of ovarian cancer calculated by GSEA. (D) Boxplots showing the signature score of 16 cancer cell states in high/low-risk groups of ovarian cancer scored by GSVA. Paired two-sided Wilcoxon test. The asterisks represent the statistical P-value (*p<0.05; **p<0.01; ***p < .001; ****p < 0.0001; ns p>0.05).

Moreover, we illustrated and compared the mutation profiles of ovarian cancer patients in the high-risk and low-risk groups (Supplementary Figure 4A). Surprisingly, we found that TMB scores, which have vital potential to predict immunotherapy responses, were markedly elevated in the low-risk group patients, and negatively correlated with the CD8⁺ T cell-associated risk score (Supplementary Figure 4B–4C).

Dissection of tumor immune microenvironment features based on CD8⁺ T cell-related prognostic signature

To unveil the disparities in the immune milieu between the high-risk and low-risk group ovarian cancer patients, we first calculated the immune score using the ESTIMATE method. The results indicated that the immune score was dramatically higher in the low-risk group patients compared to high-risk patients, and it was negatively associated with this risk score (Figure 4A, 4B). Moreover, we interrogated a series of immune signature score in two risk group patients, and found that the low-risk group showed higher scores in signatures related to antigen processing and presentation, B cell receptor (BCR) signaling pathway, T cell receptor (TCR) signaling pathway, and interferons (Figure 4C). Subsequently, the CIBERSORT algorithm was performed to speculate the infiltration levels of diverse immune cell types in the ovarian cancer immune microenvironment. It was obviously indicated that high-risk patients had lower fractions of the M1-like macrophages, CD8⁺ T cells, regulatory T cells and follicular helper T cells (Figure 4D). In addition, the risk score showed a significantly correlation with the infiltration levels of the M1-like macrophages, CD8⁺ T cells (Figure 4E, 4F). Finally, we explored the potential association between our CD8⁺ T cell-related risk score and different immunomodulators. As shown in the bar plot, the risk score exhibited markedly positive correlation with C10orf54 and ICOSLG, but significantly negative correlation with CXCL10, CD40LG and IFNG (Figure 4G). Taken together, these results indicated us the high/low-risk patients had distinct immune cell infiltration and immune features, which might potentially contribute to different outcomes of ovarian cancers.

Figure 4. Dissection of tumor immune microenvironment features based on CD8⁺ T cell-related prognostic signature. (A) Immune score in high/low-risk groups. (B) Spearman correlation between Immune score and risk score. (C) Boxplots showing the signature score of 17 immune pathways in high/low-risk groups of ovarian cancer scored by GSVA. Paired two-sided Wilcoxon test. (D) Boxplots showing the proportion of 22 immune cells in high/low-risk groups of ovarian cancer estimated by CIBERSORT. Paired two-sided Wilcoxon test. (E, F) Scatter plots showing the correlation between the risk score and the proportion of M1-like macrophages and CD8⁺ T cells (G) Bar plot of the correlation between immunomodulators and the risk score. The asterisks represent the statistical P-value (*p<0.05; **p<0.01; ***p < .001; ****p < 0.0001; ns p>0.05).

High/low-risk group patients differ in response to immunotherapy

Considering the evident correlation between the CD8⁺ T cell-associated risk score and diverse immune checkpoints, we subsequently devoted to figuring out whether the risk score correlated with the tumor immune exclusion score, which was recognized as one of immunotherapy predictors [22]. As expected, we found that ovarian cancer patients in the high-risk group harbored a higher tumor immune exclusion score, which positively associated with the risk score (Figure 5A, 5B). To evaluate the response of ovarian cancers in different risk groups to immune checkpoint inhibitors (ICIs), we calculated the IPS score (including CTLA4 blocker and PD-1/PD-L1/PD-L2 blocker) was much significantly higher in the low-risk group patients (Figure 5C–5E). This suggested that ICIs targeting CTLA4 or PD-1/PD-L1/PD-L2 are more suitable for low-risk patients rather than the high-risk group patients. Finally, to further validate the value of our CD8⁺ T cell-related prognostic model in predicting the immunotherapy response, we analyzed RNA-seq data of pretreatment samples of 298 bladder cancer patients from the IMvigor210 cohort before anti-PD-L1 treatment. We found high-risk patients were mainly in the non-response group, and patients with higher risk score also exhibited inferior outcomes as showed in the Kaplan-Meier analysis (Figure 5F, 5G). Collectively, the above findings indicate that patients in the low risk-group are more likely to benefit from immunotherapy, and the CD8⁺ T cell-related risk score could act as a potential biomarker to distinguish ovarian cancer patients who may benefit from immunotherapy.

Figure 5. High/low-risk group patients differ in response to immunotherapy. (A) Tumor Immune Exclusion score in high/low-risk groups. (B) Spearman correlation between Tumor Immune Exclusion score and risk score. (C–E) Immunophenoscore (IPS) in high/low-risk groups. (F) Proportion of immunotherapy responses between high/low-risk groups. (G) Kaplan-Meier survival analysis was performed on the relationship between the risk score and OS in the IMvigor210 immunotherapy cohort.

Distinct drug sensitivities were observed among high/low-risk group patients

The substantial implication of tumor heterogeneity lies in the fact that distinct groups of patients exhibit varying responses to anti-cancer drug treatments, thereby exacerbating the risk of treatment failure and recurrence. Therefore, we determined to explore the potential relationship between our prognostic model and diverse anti-cancer drug sensitivity. The IC50 values of multiple drugs in ovarian cancer patients were calculated and compared between the two risk groups as shown in the bubble plot (Figure 6A). Interestingly, the results suggested that the IC50 values of cisplatin, cyclophosphamide and docetaxel were significantly higher in the high-risk patients, and positively associated with the CD8⁺ T cell-related risk score (Figure 6B–6D). It indicated that ovarian cancer patients in the high-risk group might be more resistant to standard chemotherapy regimens. Additionally, we also assessed the correlation between the model genes and these anti-cancer drugs. Specifically, the expression levels of WNK1 and RNPEPL1 were both conspicuously positively related to the IC50 of GSK591. On the other hand, the IC50 values of cyclophosphamide negatively correlated with the expression levels of BTN3A3, S1PR4, CCR7, CD226 and CD40LG (Figure 6A). Summarily, our findings suggest that this novel CD8⁺ T cell-related prognostic model can be a reliable predictor for screening efficient drugs in ovarian cancers.

Figure 6. High/low-risk group patients differ in drug sensitivity and response to immunotherapy. (A) Bubble plot showing the relationship between drugs, risk score, and model genes. (B) Violin plot showing the comparison of IC50 of cisplatin between high/low-risk groups, and scatter plot showing the correlation between the IC50 of drugs and the risk score. (C) Violin plot showing the comparison of IC50 of cyclophosphamide between high/low-risk groups, and scatter plot showing the correlation between the IC50 of drugs and the risk score. (D) Violin plot showing the comparison of IC50 of docetaxel between high/low-risk groups, and scatter plot showing the correlation between the IC50 of drugs and the risk score.

Construction of a nomogram to forecast survival for ovarian cancer

Subsequently, we conducted univariate and multivariate Cox analyses to determine whether the CD8⁺ T cell-associated risk score was a robust independent prognostic factor. The univariate Cox regression analysis revealed that the risk score was identified as a risk factor for ovarian cancer patients (Figure 7A). When adjusting for other confounding factors, the multivariate analysis also revealed that the risk score was still a prognostic factor independent of age, tumor stage and tumor grade (Figure 7B). Moreover, a predictive nomogram was established to improve the prognosis efficacy of our prognostic model and provide a quantitative and visual tool for estimating outcomes of ovarian cancer patients in the TCGA cohort (Figure 7C). As demonstrated by the calibration plots, the nomogram exhibited robust stability in quantifying the survival probabilities of ovarian cancer patients at 1, 2, 3, and 5 years (Figure 7D). These above results indicated that this constructed nomogram exhibited excellent prediction efficacy for ovarian cancer patients.

Figure 7. Construction of a nomogram to forecast survival for ovarian cancer. (A) Univariate analysis for the clinicopathologic characteristics and the risk score in TCGA cohort. (B) Multivariate analysis for the clinicopathologic characteristics and the risk score in TCGA cohort. (C) A nomogram was established to predict the prognostic of ovarian cancer patients. (D) Calibration analysis of the nomogram. (E) ROC curve analysis of the nomogram. The asterisks represent the statistical P-value (*p<0.05; **p<0.01; ***p < .001; ****p < 0.0001; ns p>0.05).