Development and verification of a 7-lncRNA prognostic model based on tumor immunity for patients with ovarian cancer

Identification and preliminary evaluation of two subtypes of OC

According to the immunological characteristics of 379 tumor samples in the TCGA-OC cohort, we divided them into high immunity group (Immunity_H) and low immunity group (Immunity_L) based on 29 immune gene sets along with ssGSEA algorithm. R packages (“estimate”, “limma”) were used to calculate the Immune score, Stromal score, ESTIMATE score and Tumor Purity of the two subtypes (Fig. 1A). The heatmap of immune responses based on the ESTIMATE algorithms and single-sample GSEA (ssGSEA) is depicted in Fig. 1B. We further used tSNE algorithm for clustering analysis of TCGA-OC and obtained similar classification results (Fig. 1C). In addition, the results revealed the Immune Score, Stromal Score, and ESTIMATE Score of the Immunity_H was higher than that of the Immunity_L. Moreover, the violin plot also showed significant differences in Immune Score, Stromal Score and ESTIMATE Score between the two subtypes (Fig. 1D). We further explored the expressions of HLA genes between the two subtypes and discovered that the expressions of all HLA genes in Immunity_H were significantly higher than that in Immunity_L (Fig. 1E). These results illustrated the significance of our classification of OC into two subtypes, which could largely distinguish the characteristics of OC.

Fig. 1figure 1

Immune subtypes and clustering in OC patients. A Based on the results of ssGSEA, OC patients were divided into Immunity_H and Immunity_L by hierarchical clustering algorithm. B Immune infiltration and tumor microenvironment landscape of TCGA-OC patients. C Verification of immune subtypes by tSNE. D The comparison of Immune Score, Stromal Score and ESTIMATE Score in Immunity_H and Immunity_L groups. E Comparison of the expression levels of HLA genes between two subtypes

Detection of immunity‐related module and hub genes by WGANA

In WGCNA analysis, we identified 13 co-expression modules and analyzed their association with the immune-related cluster from ssGSEA. Based on the lncRNAs, a co-expression network was established by R package “WGCNA”, which could reveal the modules and genes that were significantly associated with the immunity cluster (Fig. 2A). In this study, β = 5 was the best choice for soft thresholds to construct a scale-free network (Fig. 2B). We next visualized the gene network with the meta-modules (Fig. 2C). After adjusting the parameters of WGCNA, we classified the DEGs into 13 modules (Fig. 2D). The results indicated that purple module was the most correlated module of immunity-cluster (r = 0.69, P = 7e-34, Fig. 2E). There were 141 genes in the purple module (Table S1). In the module-trait analysis, 8 genes with GS value > 0.3 and MM value > 0.8 were defined as hub genes: CCDC69, CLMP, FAM110B, FAM129A, GUCY1B3, PALLD, PLEKHO1, and STY11. Afterward, we defined genes in the purple module as stemness-related hub genes (Fig. 2F). These results suggested that the genes in the purple module was significantly related to the stemness of OC cells.

Fig. 2figure 2

Detection and validation of immunity-related module by WGCNA. A The cluster was based on the transcriptome data from TCGA. The color intensity represents the immunity cluster. B Analysis of the scale-free fit index for various soft-thresholding powers and the mean connectivity for various soft-thresholding powers. C The heatmap identified groups of correlated eigengenes termed meta-modules. D TOM cluster dendrogram of WGCNA: Branches with different colors corresponding to different modules. Dynamic Tree Cut represents the original module, while Merged Dynamic represents the final module. E Heatmap of the correlation between gene modules and the immunity cluster of ovarian cancer. The purple module was the most significant module with immunity. F Scatter plot of module eigengenes in the purple module

Function and DEGs of lncRNA related genes

We next conducted the correlation analysis for the key lncRNAs resulted from the WGCNA. Finally, 1269 genes were identified. Functional annotation analyses of the selected genes were then performed. GO enrichment showed that the lncRNA related genes were mainly involved in “autophagy”, “process utilizing autophagic mechanism”, “transcription regulator complex”, “focal adhesion”, and “protein serine/threonine kinase activity” (Fig. 3A). KEGG pathway enrichment analysis showed that the target genes were mainly involved in “MAPK signaling pathway”, “MicroRNAs in cancer”, and “Human cytomegalovirus infection” (Fig. 3B). Meanwhile, we compared the DEGs between normal and cancer tissues of the related genes. A total of 72 DEGs, were identified including 31 down-regulated genes and 41 up-regulated genes (Fig. 3C). The heatmap of DEGs were shown in Fig. 3D. These outcomes indicated that the key genes were functional in the progression of OC cells.

Fig. 3figure 3

Enrichment analysis and differentially expressed genes of immune-related lncRNAs. A GO enrichment analysis and B KEGG pathway enrichment analysis was performed. C Volcano plot of DEGs. D Heatmap of the DEGs

Consensus cluster analysis for selected key gene expression profiles

Then, we performed the consensus clustering analysis to investigate the relationship between these prognostic genes and OC subtypes. According to the CDF value, we classified the 379 OC patients into three clusters (k = 3, Fig. 4A–D). Cluster 1 (n = 226), cluster 2 (n = 80), and cluster 3 (n = 73) were generated from a total of 379 patients. We used principal component analysis (PCA) to display differences in gene expression levels among the three subgroups (Fig. 4E). We also found that the patients from cluster 2 tended to survive longer than the patients from cluster 1 and cluster 3 (Fig. 4F), implying a significant prognostic value of these DEGs.

Fig. 4figure 4

Identification of the molecular subtypes of the OC patients using the DEGs associated with prognosis. A The OC patients were stratified into 3 clusters based on the consensus clustering matrix (k = 3). B-D Consensus clustering model with cumulative distribution function (CDF) by k from 2 to 9. E The results of PCA analysis among the three subtypes. F Survival curves of patients in the three clusters

Establishment and validation of the risk signature based on lncRNA expression

We then constructed a risk model by the lncRNA resulting from the purple module. First, we selected these genes to conduct an additional LASSO regression analysis on 136 lncRNAs (Fig. 5A-B). Table 1 listed the genes and coefficients used to calculate each subject’s risk score. The formula was as follows: Risk score = (AL391832.3*0.547455232)-(AC078788.1*0.219104082)- (AC104971.3*0.240874313) –(LINC00892*0.060063563) + (LINC02207*0.115620759)- (LINC02416*0.160019691) –(PSMB8.AS1*0.015774278). The risk scores of OC patients in TCGA were evaluated, and all patients were divided into high-risk group and low-risk group according to the median risk score (Fig. 5C). There was no doubt that the mortality rate in the high-risk group was considerably higher than that in the low-risk group (Fig. 5D). Differential-expression levels of the 7 lncRNA and clinicopathological features in the high- and low-risk groups are shown in heatmaps (Fig. 5E). The results showed that living status, tumor residual disease, tumor status, recurrence, grade, stage, and neoplasm subdivision were differentially distributed in the two risk groups. The correlation analyses were also performed to check the expression correlation between the hub genes. (Fig. 5F). To evaluate the role of the 7-lncRNA signature in OC, we drew K-M curves for the high- and low-risk groups of the TCGA-OC cohort (Fig. 5G). These two subgroups significantly differed in OS (P < 0.01). Thereafter, we used a time-dependent ROC curve to predict the efficacy of the risk signature. The area under the ROC curve (AUC) of the prediction model was 0.72 of the OS (Fig. 5H). The contents of seven lncRNAs in different immunity groups were also compared. As shown in supplementary Fig. S1, the expressions of seven lncRNAs were all elevated in immunity_H group. We also validated the function of PSMB8-AS1 in SKOV3 cell line in vitro and in vivo. As shown in supplementary Fig. S2, the metastatic capability significantly decreased after knocking down PSMB8-AS1. These results suggested that the 7 lncRNAs play essential roles in the progression of ovarian cancer.

Fig. 5figure 5

Construction of the immune-related lncRNA risk model in TCGA cohort. A LASSO coefficient profiles of the common genes. B Cross-validation for tuning parameter screening in the LASSO regression model. C, D Distribution of risk score, survival time and survival status. E Heatmap of the 7 lncRNA expression pattern in clinicopathologic characteristics and risk score in the TCGA database. F Correlation between the 7 lncRNAs. G Survival analysis for different risk groups in the combined TCGA-OC patients. H Time-dependent ROC curve analysis of the risk model

Table 1 Seven immune cluster associated genes and corresponding coefficient valueFunctional analysis of the risk score model

We further investigated the correlation between the risk score and ESTIMATE related score including immune score, stromal score, and estimate score. We found a low positive relationship between Immune Score, Stromal Score and risk score with r = 0.13 and 0.28, respectively (both p < 0.01, Fig. 6A-B), which pointed out that stromal and immune cell was higher in the high risk group. However, the relationship between Tumor Purity and risk score was negatively correlated (Fig. 6C). These results indicated that patients with an unfavorable prognosis in the high risk group associated with the variation in tumor immune microenvironment of OC. To clarify the important pathway of signature enrichment related to the risk signature, we conducted GSEA. Finally, 55 enrichment pathways with significant variations between low and high-risk groups were identified at the criteria of FDR < 0.25, P-value < 0.05. The top five signaling pathways in the high-risk group were calcium signaling pathway, cell cycle, fatty acid metabolism, GnRH signaling pathway, and mismatch repair. On the other hand, the top five signaling pathways in the low-risk group were P53 signaling pathway, pyrimidine metabolism, regulation pf actin cytoskeleton, TGF-β signaling pathway, and tight junction (Fig. 6D). Furthermore, we stratified the patients into four subgroups according to the immune score and risk score. The result indicated that patients with high immune score and low risk score had the most favorable prognosis. However, patients with high immune score and high risk score had the worst prognosis (Fig. 6E). These results illustrated the relationship between ESTIMATE score and risk score was significant, and the potential function of the risk signature was also meaningful.

Fig. 6figure 6

The association between tumor microenvironment and risk score. Discrepancy of A Immune Score, B Stromal Score and C Tumor Purity in two groups. D Functional enrichment analysis between low- and high-risk groups of top five signaling pathways in the high- and low-risk subgroup. E Survival analysis for four groups stratified by combining the immune signature and the risk score characteristic in the TCGA-OC cohort

Construction and validation of the prognostic-nomogram model

Next, we performed univariate and multivariate Cox regression analyses in the TCGA-OC patients to assess the independent prognostic value of the lncRNA related risk signature. We observed that in univariate analysis, age, stage, tumor status, tumor residual, and risk score were significantly correlated with prognosis (Fig. 7A). Furthermore, multivariate analysis indicated that age, tumor status, tumor residual, and risk score were independent prognostic factors in the TCGA-OC patients (Fig. 7B; both P < 0.05). A nomogram model based on four independent risk factors was established to evaluate the prognostic significance of the risk signature in OC patients (Fig. 7C). The corresponding score of each variable is shown in Table 2. The calibration curves revealed a favorable consistency between expected and observed survival rates (Fig. 7D). Then patients with OC were divided into three subgroups evenly according to the total points from the nomogram namely low-, moderate-, and high-score group. The overall survival curve of the three groups was shown in Fig. 7E. The results showed that patients with high scores had the worst prognosis. What’s more, the ROC showed that nomogram could accurately predict the survival outcome of patients, and the AUC values of 1, 3 and 5 years were 0.780, 0.823 and 0.837 respectively (Fig. 7F). Taken together, the results described above suggested that the nomogram model had good reliability in predicting OS in OC patients.

Fig. 7figure 7

Construction and validation of the nomogram based on immune-related lncRNA signature and clinicopathological features. A-B Univariate and multivariate Cox regression analyses in the TCGA cohort. C The nomogram was established using age, tumor status, tumor residual, and the risk signature in the TCGA-OC cohort. D Calibration diagram of the nomogram for predicting the probability of OS at 1, 3, and 5-years. E Survival curve of patients in low, moderate, and high score groups. F Prediction of the nomogram indicated by AUC based on clinical traits and risk score

Table 2 Corresponding risk score for each variable and total scoreValidation of the lncRNA-related risk signature in GEO database

To assess the predictive value of the risk model, we used the risk score algorithm in the GSE datasets. The results in the validation cohort revealed that OC patients in the high-risk group had worse OS and PFS rates in GSE17260 (Fig. 8A-B), and OS in GSE14764 (Fig. 8C) than those in the low-risk group. The AUCs for survival were 0.774, 0.759, and 0.786, respectively (Fig. 8D-F). These findings suggested that the 7-lncRNA risk model could accurately predict the prognosis of patients with OC.

Fig. 8figure 8

Verification of the immune-related signature in two independent cohorts. A OS plot of the lncRNA-related signature in the GSE17260. B PFS plot of the signature in the GSE17260. C OS plot of the signature in the GSE14764. D Time-dependent ROC curve of the lncRNA-related signature in GSE17260 of OS. E Time-dependent ROC curve of the signature in GSE17260 of PFS. F Time-dependent ROC curve of the signature in GSE14764 of OS

Chemotherapeutic drug sensitivity analysis and validation

To observe the differences in drug sensitivity of commonly used chemotherapeutic agents between the different risk groups, drug selection was used. “pRRophic” package is a method used to predict sensitivity of some kinds of chemotherapy drugs. By using the “pRRophetic” package for drug sensitivity analysis, we observed that patients in the high-risk group were more sensitive to rapamycin (Fig. 9A). SKOV3 is derived from ascites isolated cells from patients with ovarian cancer. It has resistance to some chemotherapy drugs including cisplatin and adriamycin. A2780 is similar with SKOV3 cell line. Therefore, we further validated the function of rapamycin by in vivo experiments with SKOV3 and A2780 cell lines. Transwell and wound healing experiments indicated that rapamycin inhibited invasion and metastasis in OC cell lines (Fig. 9B-C).

Fig. 9figure 9

In vitro validation of rapamycin in OC cell lines. A Differences in drug sensitivity of rapamycin between high- and low-risk groups. B MTT assay, C transwell assay, and D gap closure analysis were performed in SKOV3 cell line after treated with rapamycin at the concentration of 20um. E MTT assay, F transwell assay, and G gap closure analysis were also performed in A2780 cell line after treated with rapamycin at the concentration of 20um. NC, Normal Control

留言 (0)

沒有登入
gif