Identification of m6A-related lncRNAs
Pearson correlation analysis was used to assess the relationship between 14086 lncRNAs and 34 m6A regulators. Total 112lncRNAswith absolute correlation coefficient > 0.40 and P value < 0.001were considered as m6A-related lncRNAs. Univariate Cox regression analysis was used to explore the prognostic roles of m6A-related lncRNAs. Of the 112 m6A-related lncRNAs, 7 were associated with the OS (Table1). These results indicated that the 7 m6A-related lncRNAs, including AC024270.4, AC008124.1, AL109811.2, AC015922.2, AC099850.4,AC025176.1, and RPP38-DT,might be potential prognostic biomarkers of CC, termed as m6A-relatedprognosticlncRNAs.
Expression profiles of m6A-relatedprognosticlncRNAs
To explore the potential biological function of m6A-related lncRNAs in the occurrence of CC, we compared the expression profiles of 7 m6A-related prognostic lncRNAs between CC samples and normal tissues. Notably, the tumor samples showed significantly lower expression levels of AC024270.4, AC008124.1, AL109811.2, and AC015922.2, but higher levels of AC099850.4, AC025176.1, and RPP38-DT, compared with the normal samples (Supplementary Fig.1).These findings suggested that the 7 m6A-related prognostic lncRNAs might possess important biological roles in the development of CC.
Consensus clusteringpatterns of m6A-relatedprognostic lncRNAs
The “ConsensusClusterPlus” package, using the7 m6A-related prognostic lncRNAs, was utilized to explore the molecular subtypes of patients.According to the cumulative distribution function (CDF), the area under the CDF curve, and the tracking plot from k = 2 to 9 (SupplementaryFig. 2a-2c), the k = 2 was identified as the most appropriate choice to divide patients into two different m6A-relatedlncRNAclusteringpatterns(Fig.1a), including209 casesinm6AlncRNAclusterAand 64 cases in m6AlncRNAclusterB.m6AlncRNAclusterA had a notably better outcome compared withclusterB(Fig.1b). In addition, the heatmap revealed that m6AlncRNA clusterB was preferentially related to a highFIGOstage (Fig. 1c).
TME immuneinfiltrationcharacteristics ofdifferent m6AlncRNAclusteringpatterns
GSVAwas applied toexplore the biological behaviors betweendifferentm6AlncRNAclusteringpatterns. m6AlncRNA clusterA presented enrichment pathways related to oxidative phosphorylation,cardiac muscle contraction, histidine catabolism, and arachidonic acid metabolism (Fig.1d). m6AlncRNAclusterB was significantly enriched in immune evasion, stromal, and tumorigenic activationpathways such as TGF beta signaling pathway, ubiquitin mediated proteolysis, focal adhesion, and pathways in cancer. Subsequently, we further compared TME cell infiltrates between two m6AlncRNAclusters. ClusterA showed higher infiltration levels of multiple immune cells such asactivated B cellandactivated CD8 T cellthan clusterB (Fig.1e). The TMEcell-infiltratingcharacteristic of clusterA was consist with its matching survival advantage. As expected, clusterA exhibited higher immune score(Fig. 1f)and ESTIMATE score (Fig. 1g),suggesting that clusterA had a significantly higher immune cell content and lower tumor purity. However, no significant difference of stromal score was displayed between two clusters (Fig.1h). These results indicated that the two distinct m6AlncRNAclustering patterns had markedly different TME.
Generation ofm6AlncRNAgenesand KEGGpathwayenrichmentanalysis
To further explore the potentialbiological behaviorsof each m6AlncRNA clustering pattern,we determined 786DEGsbetween two m6AlncRNAclustersusing the limma package, named as m6AlncRNAgenes. Then, we used the clusterProfiler package to perform KEGGenrichment analysis for the DEGs. Fig.2a showed the pathways with significant enrichment. To our surprise, these genes presented enrichment of pathways associated with PD-L1 expression and PD-1 checkpoint pathway and infection-related pathways such as viral carcinogenesis and Epstein-Barr virus infection. Afterward, we utilized univariate Cox regression analysis to explore the effect of DEGs on the survival of patients. Among the 786 genes, 140 were positively or negatively related to the OS with Pvalue < 0.05, regarded as m6AlncRNAprognosticgenes(Supplementary Table1).
Consensus Clustering of m6AlncRNAprognosticgenes
To further assess the regulation mechanism of m6AlncRNAclustering patternin CC, we subsequently performed consensus clustering analysis based on the 140 prognosticDEGsso as to divide patients. The consensus clustering of the 140 m6AlncRNAprognostic genesclassified patients into two different genomic subtypes, considered as geneclusterA(n = 224) and geneclusterB (n = 49), respectively (Supplementary Fig.3a-3cand Fig.2b). We found that 4 out of the 7 m6A-related prognostic lncRNAs showed significantly different expression levels in the two geneclusters (Supplementary Fig.4a).Gene clusterB had significantly better prognosis than gene clusterA (Fig. 2c). Moreover, the heatmap showed gene clusterA was preferentially associated with m6AlncRNA clusterB (Fig. 2d).
Generation of m6AlncRNAscore and prognosticvalue
To reveal the role of 140 prognostic DEGs in CC, we used PCA toconstruct a scoring system to quantify them6AlncRNA clustering pattern in each patient,termed as m6AlncRNAscore. We then divided patients into thehigh-m6AlncRNAscoregroup (n = 147) and the low-m6AlncRNAscore group(n = 126) according to the cutoff value -0.85 determined by the survminer package. Of the 7 m6A-related prognostic lncRNAs, 6 displayed significantly different levels between two different m6AlncRNAscore groups (Supplementary Fig. 4b).A better prognosis was observed in the high-m6AlncRNAscore subgroup (Fig.3a). The alluvial diagram showed the corresponding relationship between m6AlncRNA cluster grouping, gene cluster grouping,m6AlncRNAscoregrouping, and survival outcomes (Fig. 3b). The matching rates of m6AlncRNA clusterA with high-m6AlncRNAscore and m6AlncRNA clusterB with low-m6AlncRNAscore were 67.5% and 90.6%, respectively. To assess the accuracy of m6AlncRNAscore in predicting the OS, we performed ROC analysis and found that the 3-year AUC value was 0.708, implying that m6AlncRNAscore had a good prognostic discrimination performance (Fig. 3c). Subsequently, we compared the m6AlncRNAscore between different clustering subtypes. The m6AlncRNAscore in m6AlncRNA clusterA, as expected, was dramatically higher than that in m6AlncRNA clusterB (Fig. 3d). Similarly, gene cluster B had significantly higher m6AlncRNAscore compared with gene cluster A (Fig. 3e).Further stratified survival analysis results showed thatthe OS time in the low-m6AlncRNAscore group was dramatically shorter compared with the high-m6AlncRNAscore group, no matter for patients with grade1/2, grade3/4, age ≤ 60 years, age > 60 years, stage I/II, or stageIII/IV(Fig. 3f-3k).
Independentprognosticvalueof m6AlncRNAscore in theprognosisof CC
As shown in Fig. 4a, the univariate Cox analysis resultsshowed thatm6AlncRNAscore, age, and FIGO stageweresignificantly related to the OSof CC patients. Subsequent multivariate Cox analysisresultsdisplayedthatage was an independent risky factor (HR = 2.157, Pvalue = 0.015), but m6AlncRNAscore was an independent protective factor (HR = 0.918, Pvalue < 0.001) for the OS of CC patients (Fig. 4b).
Considering the significance ofPFSandDSSin tumor prognosis, we further validated theprognostic value ofm6AlncRNAscore in the PFS and DSS. In the univariate analysis, high-m6AlncRNAscore was significantly associated with better PFS (Fig. 4c) and DSS (Fig.4d). Moreover, further multivariate Cox regression analysis results showed that m6AlncRNAscore was not only an independent prognostic factor for the PFS (Fig.4e), but also an independent prognostic factor for the DSS (Fig. 4f).Our results strongly indicated that m6AlncRNAscore had goodprognostic value in CC.
TMEimmuneinfiltrationcharacteristics ofdifferentm6AlncRNAscoregroups
Toverify the biological behaviors ofm6AlncRNA clustering patterns in TME, we performed GSVA and ssGSEA analyses in two m6AlncRNAscore groups. The high-m6AlncRNAscore groupwas characterized by enrichment of hallmark pathways such asoxidative phosphorylation andcardiac muscle contraction (Fig.5a)and infiltration of activated CD8 T cell, CD56dim natural killer cell, and monocyte (Fig. 5b). The low-m6AlncRNAscore group was characterized by enrichment of immunosuppressive, stromal, and carcinogenic activation pathways such as wnt signaling pathway, TGF beta signaling pathway, MAPK signaling pathway, ERBB signaling pathway, focal adhesion, extracellular matrix (ECM)-receptor interaction, and pathways in cancer (Fig. 5a). Besides, the low-m6AlncRNAscore group was rich inT helper 2 (Th2)and regulatory T (Treg) cells,two types of tumorimmunosuppressiveT cells(Fig.5b).We then explored the expression profiles of immunosuppressive factors IL-10 (Fig.5c) and TGF-beta1 (Fig.5d) and found that their levels in the low-m6AlncRNAscore group were significantly higher. The above results indicated again that m6A-related lncRNA clustering patterns played a vital role in shaping TME landscape.
Clinical andsomaticmutationcharacteristicsofdifferent m6AlncRNAscoregroups
As expected, the m6AlncRNAscore was significantly higher in patients with stage I/II than those with stageIII/IV(Fig. 6a). However, no significant m6AlncRNAscore differencewasobserved in different age or grade subgroups. TMB quantification analysis results showed that the low-m6AlncRNAscore group presented no significant TBM difference in relative to the high-m6AlncRNAscore group (Fig.6b). Next, the survminer package was applied to classify patients with information of somatic mutation and survival into the high-TMB group (n = 28) and the low-TMB group (n = 226) according to the cutoff value 6.32.A better prognostic tendency was observed in the high TMB group, while no significant difference was displayed betweenthe high- and low-TMB groups(Fig. 6c). Moreover, we found that patients with low-m6AlncRNAscore and low-TMB had the worst prognosis, and the prediction power of m6AlncRNAscore was not disturbed by TMB during the first 5 years of follow-up (Fig.6d).
Patient’sresponse to ICIsindifferent m6AlncRNAscoregroups
It has been reported that IPS values could predict the response of patients toICIs. Whether patients received anti-CTLA-4 (Fig. 7a), anti-PD-L1 (Fig.7b)or anti-CTLA-4 and anti-PD-L1 combination treatments (Fig.7c), the IPS values of the high-m6AlncRNAscoregroup were dramatically higher compared with the low-m6AlncRNAscore group, suggesting that the corresponding ICI therapy responses in thehigh-m6AlncRNAscore groupwere significantly better than those of the low-m6AlncRNAscore group. These results indicated thatpatients with high-m6AlncRNAscore were more likely to benefit from ICIs.
Validation of the expression levels of four m6A-related lncRNAs in CC samples
qRT-PCR assay was used to detect the expression ofAC024270.4, AC008124.1,AC025176.1andRPP38-DTin 6tumor tissuesand 6normalsamples. As shown inSupplementary Fig. 5, compared with normal tissues, cervical cancer tissues had higher AC024270.4 and AC008124.1, but lower AC025176.1. There was no difference in AC024270.4 expression between tumor samples and normal samples.