3.1 Identification of mRLs in OV
Firstly, we abstracted the expression of lncRNAs from TCGA-OV dataset and identified 419 mRLs (p < 0.05, absolute correlation coefficient > 0.3) through Pearson correlation analysis with m6A regulators. The relational network between m6A regulators and mRLs was shown in Figure 1B. After collecting the expression of lncRNAs in the GSE9891 dataset, we found 61 mRLs shared in both datasets. When combined with clinical information, we screened five mRLs related to prognosis through univariate regression analysis in TCGA-train (Table S2).
3.2 Establish and verify of a prognostic signature based on mRLs in OV patients
Firstly, the 374 OV patients in were classified into train and test sets at a 1:1 ratio randomly, including 187 OV patients, respectively. Secondly, we performed the LASSO (Figure S1A-B) and multivariate regression analysis (Figure S1C) in train set to build a reliable prognostic signature. The formula was as following:
Risk core= (-0.068088177 * expr (WAC-AS1)) + (-0.276777737 * expr (LINC00997)) + (0.09155959 * expr (DNM3OS)) + (-0.137781856 * (FOXN3-AS1), expr means the expression value. WAC-AS1, LINC00997, DNM3OS, and FOXN3-AS1 were finally identified as key mRLs. Both the outcomes of the expression analyses in the GEPIA dataset and the qRT-PCR results indicated that the expression of DNM3OS and FOXN3-AS1 were higher in normal compared with tumor tissue (p <0.05) (Figure S2A-H). Afterward, we assessed the risk score of OV patients in TCGA-train and separated them into two groups. The distributions of risk score, survival status of each sample (Figure 2A), and the heatmap of key gene expression patterns (Figure 2E) in TCGA-train were performed in Figure 2. People in high risk had poor survival prognosis than low-risk groups in TCGA-train (Figure 2I). Time-dependent receiver operating characteristic curves (ROC)analysis based on the signature revealed that the AUC of 1-, 2-, and 3-years survival was 0.624, 0.694, and 0.630, respectively (Figure 2M). To verify the mRL prognostic signature, we calculated the risk score of the patients in TCGA-test (N=187), TCGA-entire (N=374), and GEO datasets (N=278) and divided patients into either group in terms of the median risk score of the TCGA-train (Figure 2B-D). The heatmap all showed that except for DNM3OS, the expression of WAC-AS1, LINC00997, and FOXN3-AS1 was higher in low-risk groups (Figure 2E-H). The Kaplan-Meier analyses of the TCGA-test, TCGA-entire, and GEO datasets revealed the same trend that the OV patients with low risk had higher OS than high-risk group (p < 0.01) (Figure 2J-L). Furthermore, time-dependent ROC curve indicated the predict capability of the prognostic signature based on mRLs. The AUC of 1-, 2- and 3-years survival in three sets were shown in Figure 2N-2P. In addition, the distributions of the groups were different in PCA plots (Figure S3). These results all demonstrated the accurate and reliability of the mRL prognostic signature in OV prediction.
3.3 The relationships between the model and clinical factors
The correlations between risk core and clinical characteristics are performed in Figure S4A-F. The risk score in age≤60 group was lower than group which age > 60 (p < 0.05, Figure S4D), and the proportion of older patients with high-risk score is 52%, while that people with low-risk score was 40% (Figure S4A). But the difference of risk score in the groups assigned by tumor grade and tumor stage was not significant (Figure S4B-C, E-F). What’s more, the risk score of high-immune score group was markedly higher than lower score patients (p < 0.05, Figure S4G), and the significant difference could be found in the risk core of Cluster1 and Cluster2 subgroups (p < 0.05, Figure S4H). Stratification analysis grouped by patient age, tumor grade, and tumor stage performed that the OS of patients with high risk was worse than that low-risk group (Figure S4I-N). By comparing the patients’ state between the groups, we also found a higher mortality of patients with high risk (Figure S5A), and the risk scores of dead patients were statistically higher than alive patients (p =0.002, Figure S5B). K–M curves of OS for patients subjected to chemotherapy performed that people with high risk had poorer prognoses than those with lower risk (p <0.001, Figure S5C). In patients with BRCA1, low-risk groups also had a favourable survival outcome (Figure S5D).
3.4 Clinical application of the signature
Univariate and multivariate regression analyses in TCGA-train, TCGA-test, TCGA-entire, and GEO datasets all indicated the independence of the mRL model (Table S3). Moreover, patient age and tumor stage were also crucial poor prognostic factors for OV. To advance the predict performance of the signature, we built a nomogram comprising the risk score, age, grade, and stage in both TCGA (Figure 3A) and GEO datasets (Figure S6). Figure 3B-D revealed that the predicted rates of OS were highly consistent with the observed rates. In addition, the prognostic signature showed superior predictive ability than clinical characteristics, and when combined with clinical characteristics, the model showed better predictive power than the signature used only (Figure 3E).
3.5 The correlations between the mRL prognostic signature and TME
GSEA indicated that pathway related to immunity were enriched in people with high-risk score (Figure 4A-B). Firstly, we implemented the ESTIMATE algorithm which indicated that the stromal score, immune score, and ESTIMATE score were higher in the high-risk group (p < 0.05), and the risk score had positive association with them (Figure 4C-H). The same trends were found in the GEO dataset (Figure 4I-N). Then through comparison of 22 TIICs between the groups, we identified that Dendritic cells resting and Macrophages M2 were positively related to risk score and (p < 0.05) and the infiltration level of T cells follicular helper and T cells regulatory (Tregs) had negative correlation with risk score (p < 0.05) in TCGA dataset (Figure 5A-I). In the GEO dataset, the risk score had positive correlation with T cells gamma delta and T cells CD4 memory resting, while risk score has negative relationship with Tregs, NK cells resting, and Dendritic cells activated (Figure S7). The differences in immune score and abundance of TIICs between the groups showed that patients with high risk had a repressive immune phenotype, which could partly explain the poorer OS of people with high risk.
To identify the immunological role of the four mRLs, we assessed their association with the TIICs through CIBERSORT algorithm. We found DNM3OS had significantly association with T cells follicular helper, Macrophages M2 and B cells memory (p < 0.001); WAC-AS1 was correlated with T cells follicular helper and Macrophages M2 (p < 0.001) (Figure S8A). Besides, WAC-AS1 (R = -0.15, p = 0.0041) and LINC00997 (R = -0.23, p < 0.001) had negative relationships with immune score, whereas DNM3OS (R = 0.13, p =0.008) and FOXNS-AS1 (R = 0.19, p < 0.001) were positively related to immune score (Figure S8B). Next, to clarify the potential mechanisms of the lncRNAs in tumor progression, we used GSEA analysis to find the enriched pathways of the mRLs in OV (Figure S8C-F). DNM3OS was highly enriched in the calcium signaling pathway, the focal adhesion pathway, the hematopoietic cell lineage pathway, the neuroactive ligand-receptor interaction pathway, and cancer pathways.
The expression level of PD-L2 was higher in tumour tissues than in normal tissues (p < 0.05, Figure 6A), and the it was markedly higher in people with high risk score in both TCGA and GEO dataset (p < 0.05, Figure 6C-F), suggesting that people in high risk promise to benefit less from immune checkpoint inhibitor (ICI). In addition, LINC00997 showed a high correlation to PD-L2 (Figure 6B), so the function of LINC00997 in the tumour immune microenvironment deserved further research. Comparing the IPSs between the groups, low-risk groups had significantly higher IPSs, suggesting a more immunogenic phenotype (Figure 6G-J). These results all demonstrated that the signature promised to predict the efficacy to immunotherapy for OV patients.
3.6 Consensus clustering for mRLs related to OV prognosis and TME
After the consensus clustering, OV patients in TCGA dataset were classified into Cluster1 and Cluster2 subgroups in the light of the expression of the significant prognostic mRLs (Figure 7A-D). The heatmap of the expression pattern between cluster1 and cluster2 subgroups was performed in Figure 7E. In addition, the stroma, immune, and ESTIMATE scores were markedly higher and the tumor purity was substantially lower (p < 0.05) in Cluster2 (Figure 7G-J). Besides, the infiltration of B cells naive (p < 0.05), T cell CD4 memory activated (p < 0.05), T cells follicular helper (p < 0.05), Macrophages M1 (p < 0.05) and Mast cells resting (p < 0.05) were higher in Cluster1 subgroup, whereas the abundance of Macrophages M0 (p < 0.05), Macrophages M2 (p < 0.05), Mast cells activated (p < 0.05) and Neutrophils (p < 0.05) were higher in Cluster2 subgroup (Figure 7F), so Cluster2 subgroup tended to behave a more immune-suppressive phenotype. Moreover, the results of Kaplan-Meier analyses performed that the OS of the Cluster2 was dramatically worse than that of Cluster1 (Figure 7K), and the expression of PD-L2 was higher in the Cluster2 (Figure 7L); this indicates patients in the Cluster2 may respond more sensitively to immunotherapy. The above results all demonstrated that the expression of the mRLs influences the TME, leading to different prognoses.
3.7 The expression and immune status of lncRNAs in pan-cancer
We achieved the pan-cancer data of 33 cancers from TCGA to further study the relationships between the four mRLs and immune features, stem-like properties, and patients' prognoses in pan-cancer. The expression pattern of lncRNAs in 33cancers (Figure 8A-D) and the heatmap (Figure S9A) showed inter-tumor heterogeneity, but in most cancers, the expression level of WAC-AS1 was higher in tumor tissues, whereas the expression level of DNM3OS was opposite, which were consistent with the expression trend in OV (Figure S2). The expression distribution of the four mRLs across all 33 cancer types showed the same trend as the above results (Figure 8E). Moreover, there is no significant correlation between the key mRLs (Figure S9B). The Kaplan-Meier analyses (Figure S9C) and univariate regression (Figure S10A) analyses showed that the role of the mRLs in the survival of patients with different cancers varied.
To identify the relationships between the key mRLs and immune features in pan-cancer, we investigated the immune subtype and the landscape of correlation with stromal score, immune score, stemness scores based on DNA-methylation (DNAss)and stemness scores based on mRNA (RNAss). Through an extensive immunogenomic analysis of pan-cancer, six immune subtypes, including wound healing (C1), IFN-γ dominant (C2), inflammatory (C3), lymphocyte depleted (C4), immunologically quiet (C5), and TGF-β dominant (C6), were identified according to differences in immune, genetic and clinical features, in which C3 had great OS, C2 and C1 had poor prognosis, while C4 and C6 had the least favorable outcomes(32). We evaluated the correlations between the mRLs expression level and the six immune subtypes, and found that all of them were related to pan-cancer immune subtypes, and elevated C1, C2, and C6 subtypes with upregulated expression of DNM3OS (Figure 9A). Additionally, DNM3OS had positive correlation with stromal score and immune score, whereas WAC-AS1 and LINC00997 were negatively related to stromal and immune score (Figure 9B-D), which showed the similar correlation in OV. All the four key mRLs were negatively related to RNAss (Figure 9E). Therefore, the four mRLs were took part in the formation of specific TIME. Finally, we applied the CellMiner database to evaluate the correlation between the mRLs and drug sensitivity. A higher Z-score means higher sensitivity to one specific drug. Significantly, the expression of FOXNS-AS1 influenced drug sensitivity in various cell lines (Figure S10B).