Integration of single-cell RNA-seq and bulk RNA-seq data to construct and validate a cancer-associated fibroblast-related prognostic signature for patients with ovarian cancer

Gene Expression Omnibus (GEO) data set GSE184880↗ was obtained, which contained scRNA-seq information from seven OC samples and five ovarian tissues. Initially, single cells were assessed to ensure that each gene was expressed in at least three cells and that each cell expressed a minimum of 250 genes to generate scRNA-seq data. The percentages of mitochondria and rRNA were subsequently calculated utilizing the PercentageFeatureSet tool in the Seurat R package. Additional screening was performed on the single cells by setting them to express a minimum of 6,000 genes with UMI > 100. Finally, 40,810 cells were retained. From The Cancer Genome Atlas (TCGA), we retrieved the transcriptomic data, copy number variations (CNV) data of the Masked Copy Number Segment, single-nucleotide variant (SNV) data, and relevant OC clinical data. We accessed the OV project of the TCGA database (http://cancergenome.nih.gov/↗) to extract the transcriptomes and clinical data of 379 OC patients, whereas the information on 88 normal tissues was gathered from the GTEx database. Patients without adequate clinical data were excluded from the analyses. After excluding samples of normal tissue and tumors derived from the GEO database lacking data on follow-up and outcomes, the GSE140082↗ cohort containing 380 OC samples was retrieved for use as a validation dataset. The literature was searched for ten cancer-associated pathways, including HIPPO, TP53, NOTCH, PI3K, TGF-Beta, RAS, NRF1, WNT, MYC, and Cell Cycle [15].

We performed a re-analysis of the OC scRNA-seq data with the aid of the Seurat program [16] to characterize the CAFs’ signature fully. The first step was to exclude the cells that had either > 5000 or < 250 expressed genes, which was accompanied by log normalization of these genes. Next, the uniform various approximation and projection approach was used for the non-linear dimensionality reduction, with 15 principal components and a resolution of 0.1. With the use of the FindNeighbors and FindClusters functions, single cells were organized into a variety of distinct subgroups (dim = 30 and resolution = 0.1). Afterward, the RunTSNE function was adopted to execute t-distributed stochastic neighbor embedding (TSNE) dimensionality reduction. The four marker genes, namely ACTA2, FAP, PDGFRB, and NOTCH3 were identified as being specifically expressed in fibroblasts and were annotated accordingly. The FindClusters and FindNeighbors tools from the original method were employed to re-cluster the fibroblasts. Clusters of fibroblasts were then subjected to dimensionality reduction using TSNE. With the FindAllMarkers tool, marker genes for each CAFs cluster were determined by performing pairwise comparisons across clusters based on an adjusted p-value < 0.05, minpct = 0.35, and logFC = 0.5. Employing the clusterProfiler program, an enrichment analysis was completed with the Kyoto Encyclopedia of Genes and Genomes (KEGG) on the CAF clusters’ marker genes [17]. Subsequently, the CopyKAT in R program was utilized to discriminate between tumors, and normal cells present per sample by analyzing the CNV features across the CAFs clusters [18].

The limma program was employed to search for differentially expressed genes (DEGs) between tumor tissues and normal tissues with a |log2(Fold change) |>1 and false discovery rate (FDR) < 0.05 [19]. We next used Pearson’s correlation to determine which DEGs were most strongly associated with each CAF cluster. We used this information to determine the most important genes involved in CAF with p < 0.05 and cor > 0.4. In addition, genes associated with prognosis were discovered via univariate Cox regression analysis in the survival package at a p-value < 0.05. We conducted a least absolute shrinkage and selection operator (Lasso) cox regression analysis to minimize the number of genes, after which a multivariate Cox regression analysis was conducted using a stepwise regression technique. The equation below describes the risk signature generated from the output values of the multivariate Cox model: risk score=, whereby the risk signature gene is denoted by i, the expression profile of gene i is denoted by Expi, and bi denotes the gene i coefficients in the multivariate model. Zero-mean normalization was then performed to classify the patients into low- and high-risk categories. The timeROC program was utilized to conduct a receiver operating characteristic (ROC) analysis of the risk signature’s prediction accuracy. The validation cohort was also subjected to similar analyses.

To develop a nomogram model suitable for clinical application, we initially conducted univariate and multivariate analyses on clinicopathological and risk signature parameters. With the rms program, the factors in the multivariate model with a p-value of < 0.05 were utilized to construct a nomogram for estimating the prognosis of OC individuals [20]. The generation of the calibration curve aided in evaluating of the model’s ability to make accurate predictions. Decision curve analysis (DCA) was conducted to assess the model’s reliability.

The CIBERSORT algorithm [21], a method for evaluating immune cell infiltration, was employed to probe the distributions of 22 subtypes of immunocytes in the TCGA cohort. To additionally evaluate the TME, we used the ESTIMATE method to calculate stromal and immune scores [22].

The “maftools” R software was utilized for the purpose of constructing the genomic landscape of CAF-associated genes with SNV and CNV from the TCGA datasets.

To investigate the direct predictive value of the risk score on PD-1 therapy response, we utilized the IMvigor210 cohort, which comprised 298 patients with urothelial carcinoma and included transcriptomic data along with treatment response to immunotherapy [23].

To validate the markers at the cellular level, we retrieved the mRNA expression patterns of those markers in 41 fibroblasts and 69 OC cell lines from the CCLE platform (https://portals.broadinstitute.org/ccle↗) [24]. We used a heat map to analyze the differences in the expression of these markers between fibroblasts and OC cell lines.

The R software (v4.2.2) was applied to conduct all statistical data analyses. Pearson or Spearman correlations were utilized to create the matrices of correlation. All pairwise comparisons between the two groups were made via the Wilcoxon test. K-M curves and the Log-rank test were conducted to evaluate the significance of survival differences. The significance criterion was established at P < 0.05.

After preliminary screening, the scRNA-seq data yielded 40,810 cells. Figure S1 displays the detailed results of data preprocessing. Additionally, dimensionality reduction and log-normalization yielded 16 subpopulations, and using four marker genes-NOTCH3, PDGFRB, FAP, and ACTA2-five CAFs groups were found (Figures S2A, B). The cells of 5 CAFs clusters were extracted for further clustering and dimensionality reduction. The same clustering algorithm was applied to the CAFs clusters, identifying of five CAFs clusters (Figures S2C, D). None of the five CAFs subpopulations showed expression of the epithelial cell-specific gene, thereby validating the reliability of CAFs detection (Figure S3). The TSNE plot for the whole set of 12 distributions is presented in Fig. 1A. Consequently, five CAFs clusters were established and then analyzed (Fig. 1B). Overall, 1476 DEGs were found across the five CAFs clusters, and Fig. 1C displays the distribution of the five leading DEGs (detected as the gene markers of CAF clusters) throughout the five clusters. Moreover, Fig. 1D depicts the distribution of the five clusters per cohort. The KEGG analysis findings in Fig. 1E illustrated the enrichment of these DEGs in several different pathways, including the cGMP-PKG signaling pathway, PI3K-Akt signaling pathway, focal adhesion, etc. Furthermore, 1533 cancer and normal cells are distributed throughout the five CAFs clusters based on their CNV features (Fig. 1F).

We analyzed the features of 10 tumor-related pathways in the five CAFs clusters to clarify the links between the clusters and tumor development. Figure 2A depicts the GSVA scores of the 10 tumor-related pathways across the various CAFs clusters. Notably, CAF_1, CAF_2, and CAF_3 had a significantly greater proportion of cancer cells than the other two clusters (Fig. 2B). Additionally, we compared cancerous and non-cancer cells within each CAF cluster using GSVA scores for the 10 tumor-related pathways and found some modest variations (Fig. 2C–G). We initially computed the ssGSEA score of the marker genes of each CAF cluster using the TCGA dataset (Fig. 1C displays the five most significant DEGs from CAF clusters) to identify any associations of the CAFs clusters with prognosis. As shown by the data, tumor samples scored considerably higher on the CAF_1 and CAF_3 clusters in contrast with normal samples, whilst other CAFs clusters showed the reverse pattern, with greater scores in normal tissues compared to tumors (Fig. 3A). The OC samples from the TCGA dataset were classified into high- and low-CAF score groups as per the optimal cut-off value analyzed by the survminer R package. Samples with higher CAFs scores fared poorly in both the CAF_0 and CAF_1 clusters in contrast with those in the low-CAF score group, whereas there was no correlation between the CAF_2, CAF_3, and CAF_4 clusters and a poor outcome in OC individuals (Fig. 3B–F).

We began by screening DEGs between tumor and normal tissues to develop a risk signature. Figure 5A demonstrates the identified 5808 DEGs, with 2769 DEGs showing an upregulation and 3039 DEGs showing a downregulation. In addition, 530 genes, in particular, were strongly correlated with the prognosis-associated CAFs clusters. Moreover, univariate analysis was undertaken to determine the prognostic significance of each gene and 66 genes was shown to be prognostic genes (Fig. 4A, B). After conducting Lasso Cox regression analysis to minimize the number of genes, only 14 remained with a lambda value of 0.0412 (Fig. 4C, D). Following the execution of multivariate analysis using the stepwise regression approach, we ultimately selected seven genes to construct the risk signature, including WD repeat domain 77 (WDR77), V-set, and immunoglobulin domain containing 4 (VSIG4), selectin L (SELL), mannosidase alpha class 2 A member 1 (MAN2A1), C-X-C motif chemokine ligand 9 (CXCL9), calcium voltage-gated channel subunit alpha1 C (CACNA1C), and ETS transcription factor ELK3 (ELK3) (Fig. 4E). The ultimate equation for the 7-gene signature is as follows: risk score = (0.174*ELK3) + (0.397*CACNA1C) +(-0.181*CXCL9) + (0.262*MAN2A1) + (-0.271*SELL) + (0.262*VSIG4) + (-0.158*WDR77). After applying the z-mean normalization, we determined each sample’s risk score and categorized them as either high or low risk predicated on that score. In the TCGA dataset, the AUC values for the model for 1- to 5-year survival varied from 0.65 to 0.69. However, in the GEO dataset, they varied from 0.65 (Fig. 4F, G). The Kaplan-Meier (KM) survival analyses illustrated that high-risk patients exhibited considerably worse survival status than those at low-risk in both the TCGA and GEO datasets (Fig. 4H, I).

The next step involved analyzing the SNV mutations in each of the seven genes used to develop the risk signature. A larger number of samples were found to contain SNV mutations in CACNA1C, VSIG4, CXCL9, ELK3, and MAN2A1, but no SNV variants were found in SELL or WDR77 (Figure S4A). We examined the probability that these important genes would co-occurrence along with the top 10 highly mutated genes. Figure S4B shows that whereas mutations in the aforementioned 5 genes did not show a statistically significant likelihood of co-occurrence, mutations in CSMD3 and CACNA1C did. We discovered that relatively few samples showed gain/loss of CNV in any of the 7 genes (Figure S4C). We investigated the links between the risk genes and various molecular signatures of OC to comprehend the nature of these links. The findings indicated that Homologous Recombination Defects, Fraction Altered, and the Number of Segments were strongly negatively correlated with ELK3, while CXCL9 and SELL were considerably positively associated with these same measures (Figure S4D). We also examined the possible pathways linked to each risk gene. Finally, as illustrated in Fig. 5A and B, the JAK-STAT, the Toll-like receptors receptor, the chemokine signaling, etc., were among the 28 pathways substantially linked to these seven genes.

Based on our findings, we found that ELK3, CXCL9, MAN2A1, SELL, and VSIG4 all had positive associations with the stromal, immune, and estimate scores, whereas WDR77 had negative associations with all three (FigureA). We then evaluated the three scores among groups that were differentiated by their median expression levels. Results demonstrated a considerable variation between the high-expression and low-expression groups regarding the immune score for the ELK3, CACNA1C, CXCL9, MAN2A1, SELL, and VSIG4 genes, with the former exhibiting significant upregulation (FigureB). Strong inverse correlations were observed between M0 macrophages and CXCL9, MAN2A1, SELL, and VSIG4 (FigureC). Furthermore, analysis of correlations indicated that CXCL9, SELL, and VSIG4 exhibited a substantial positive correlation with most T cells (FigureD). S5 S5 S5 S5

T-cell immunotherapy has made remarkable progress and is promising in treating cancer [25]. As a consequence, we evaluated the IMvigor210 risk signature’s predictive usefulness for immune-checkpoint therapy. Among the 348 patients in the IMvigor210 cohort, response to anti-PD-L1 receptor blockers varied, with outcomes ranging from complete response (CR) and partial response (PR) to stable disease (SD) and progressive disease (PD). As depicted in Fig. 6A, risk scores were higher in SD/PD patients as compared to CR/PR patients. A more significant proportion of people with SD/PD were found in the high-risk group compared to those in the low-risk category (Fig. 6B). Compared to the high-risk group, the low-risk group in the IMvigor210 cohort saw more favorable clinical outcomes and had substantially improved overall survival (OS) (Fig. 6C). Figure 6D shows a substantial difference in survival duration across various risk groups for patients in Stages I and II but not for those in Stages III and IV (Fig. 6E). Evidence from this study revealed that the risk score was more accurate for individuals at an earlier stage.

Through univariate and multivariate analyses, we incorporated the clinicopathological features and risk score to enhance the risk signature’s predictive performance power. The risk signature emerged as the most influential prognostic predictor of OC in a multivariate model [hazard ratio (HR) = 1.632, 95% confidence interval (CI): 1.424–1.870, P < 0.001] (Fig. 7A, B). Consequently, the variables displayed in Fig. 7C (age, lymphatic invasion, residual disease, venous invasion, and risk score) were used to develop a nomogram. The nomogram’s ability to accurately predict real-world survival rates was proved via a calibration plot (Fig. 7D). Figure 7E further demonstrates that DCA found the nomogram to be more discriminative than the risk score and stage in identifying high-risk patients. Among the TCGA dataset, the AUC values of the risk score and nomogram were shown to be elevated in contrast with those of any other indicator, as evidenced by timeROC analysis (Fig. 7F).

We used the CCLE repository to confirm that fibroblast cell lines exhibited elevated mRNA expression levels of numerous genes (MAN2A1, CACNA1C, and ELK3) in contrast with OC cell lines (Wilcoxon test, all p < 0.001; Fig. 8A and B).

The authors declare no competing interests.

Not applicable.

Not applicable.

The datasets used and/or analyzed in this study are available from the corresponding author on reasonable request.