Skip to main content
Download PDF

Integrated analysis of N6-methyladenosine- and 5-methylcytosine-related long non-coding RNAs for predicting prognosis in cervical cancer

Hereditas volume 161, Article number: 34 (2024) Cite this article

Abstract

Background

N6-methyladenosine (m6A) and 5-methylcytosine (m5C) play a role in modifying long non-coding RNAs (lncRNAs) implicated in tumorigenesis and progression. This study was performed to evaluate prognostic value of m6A- and m5C-related lncRNAs and develop an efficient model for prognosis prediction in cervical cancer (CC).

Methods

Using gene expression data of TCGA set, we identified m6A- and m5C-related lncRNAs. Consensus Clustering Analysis was performed for samples subtyping based on survival-related lncRNAs, followed by analyzing tumor infiltrating immune cells (TIICs). Optimal signature lncRNAs were obtained using lasso Cox regression analysis for constructing a prognostic model and a nomogram to predict prognosis.

Results

We built a co-expression network of 23 m6A-related genes, 15 m5C-related genes, and 62 lncRNAs. Based on 9 m6A- and m5C-related lncRNAs significantly associated with overall survival (OS) time, two molecular subtypes were obtained, which had significantly different OS time and fractions of TIICs. A prognostic model based on six m6A- and m5C-related signature lncRNAs was constructed, which could dichotomize patients into two risk subgroups with significantly different OS time. Prognostic power of the model was successfully validated in an independent dataset. We subsequently constructed a nomogram which could accurately predict survival probabilities. Drug sensitivity analysis found preferred chemotherapeutic agents for high and low-risk patients, respectively.

Conclusion

Our study reveals that m6A- and m5C-related lncRNAs are associated with prognosis and immune microenvironment of CC. The m6A- and m5C-related six-lncRNA signature may be a useful tool for survival stratification in CC and open new avenues for individualized therapies.

Background

Cervical cancer (CC) is a highly heterogeneous cancer and ranks fourth in incidence and mortality in females worldwide, with cervical squamous cell carcinoma (CESC) as the most common type [1]. CC incidence has experienced an evident increase in recent years because of the availability of effective screening programs [2]. Human papilloma virus (HPV) is a well-documented etiological factor of CC, and thus receiving HPV vaccine remains a primary method for CC prevention [3]. CC patients are commonly treated with combinations of surgical resection, chemotherapy, radiotherapy or other therapies [4]. Nonetheless, the long-term survival and prognosis of advanced CC remains unsatisfactory [5]. Therefore, there is an evident interest in gaining a deeper understanding of molecular mechanisms behind CC carcinogenesis and discovering promising prognostic signatures capable of predicting clinical outcomes of patients reliably and accurately.

Long non-coding RNAs (LncRNAs), a group of non-coding RNAs longer than 200 nucleotides, participate in a variety of biological processes involved in tumorigenesis and progression [6, 7]. A significant body of literature has supported important biological roles of lncRNAs in CC progression, invasion and metastasis, and their potentials as prognostic biomarkers [8, 9]. N6-methyladenosine (m6A) and 5-methylcytosine (m5C) are two important forms of modifications to messenger RNA (mRNAs) and non-coding RNAs (ncRNAs), playing a role in diverse functions, such as RNA splicing and translation. Either m6A or m5C modification has three classes of components: intracellular methyltransferases (“writers”), demethylases (“erasers”), and signal transducers (“readers”) [10]. m6A methyltransferase METTL3 up-regulation has been observed in CC, and may serve as prognostic biomarkers [11]. m6A reader YTHDF1 plays a oncogenic role and is related to poor prognosis of CC patients [12]. Moreover, growing evidences are highlighting m6A-related genes and m6A-related lncRNAs as potential prognostic biomarkers in CC [13, 14]. Substantial evidences have revealed that m5C modification is implicated in tumorigenesis, cancer migration, and metastasis of a large number of cancers including CC [15,16,17]. However, there is a lack of studies on characterization of biological roles of m5C modification and m5C-related lncRNAs in the pathophysiology of CC.

Biological significance of m6A and m5C modifications implies the potential of m6A-and m5C-related lncRNAs to be used for prognostic purposes in CC. To decipher the regulatory networks of m6A-and m5C-related lncRNAs in tumorigenesis and progression of CC and explore their prognostic value, we identified m6A-and m5C-related lncRNAs in CESC samples downloaded from The Cancer Genome Atlas (TCGA). Out of them, we determined survival-related lncRNAs and subsequently, revealed two molecular subtypes. Furthermore, with prognostic lncRNAs we constructed a risk score model and a nomogram to predict prognosis in CESC. In addition, we unraveled associations of m6A-and m5C-related lncRNAs with tumor infiltrating immune cells (TIICs) in tumor microenvironment (TME) as well.

Methods

Data sources

Gene expression data (log2 (FPKM + 1), Illumina HiSeq 2000 RNA Sequencing platform) of 291 CC samples and corresponding clinicopathological characteristics were downloaded from TCGA database (https://portal.gdc.cancer.gov/) and were used as the training set. The validated set of this study used GSE44001 dataset [18] downloaded from gene expression omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), including gene expression data and survival information of 300 CC samples (Illumina HumanHT-12 WG-DASL V4.0 R2 expression beadchip).

Identification of m6A- and m5C-related lncRNAs

According to a recently published study [19], 23 m6A-related genes and 15 m5C-related genes were identified (Supplemental Table 1). Expression levels of genes associated with m5C and m6A were extracted from the TCGA dataset. Using the cor function in R3.6.1 (http://77.66.12.57/R-help/cor.test.html), Pearson correlation coefficients (PCC) were calculated between the expression levels of 23 m6A-related genes, 15 m5C-related genes, and all annotated lncRNAs. The m6A- and m5C-related lncRNAs were identified using a threshold of |PCC| > 0.3 and a significance p-value < 0.05. Finally, we identified lncRNAs significantly associated with both m5C and m6A genes and selected those that were significantly correlated with both. A co-expression network between m5C and m6A genes and their commonly associated lncRNAs was constructed and visualized using Cytoscape 3.6.1 [20] (https://cytoscape.org/).

Consensus clustering analysis

Based on the expression levels of m6A- and m5C-related lncRNAs in CESC tumor samples from the TCGA training dataset, univariate Cox regression analysis was performed using the survival package Version 2.41-1 (http://bioconductor.org/packages/survivalr/) in R 3.6.1 to identify lncRNAs significantly associated with survival prognosis, with a significance threshold set at P < 0.05. Based on expression levels of these survival-related lncRNAs, CC samples in the TCGA set were separated into different subtypes by performing Consensus Clustering Analysis using ConsensusClusterPlus package [21] (version 1.54.0, http://www.bioconductor.org/packages/release/bioc/html/ConsensusCluster Plus.html) in R language. The parameters were set as d, maxK = 10, reps = 50, pItem = 0.8, pFeature = 1, clusterAlg="hc”, distance="pearson”. Kaplan- Meier (KM) survival curves were plotted for different subtypes using survival package.

Analysis of immune infiltration

R3.6.1 GSVA (http://www.bioconductor.org/packages/release/bioc/html/GSVA.html) Version Version 1.36.3, which was based on single sample gene set enrichment analysis (ssGSEA) algorithm was adopted to evaluate the proportion of TIICs in TCGA samples. Subsequently, we performed differential distribution comparisons of TIIC proportions across different subtypes using group t-tests in R3.6.1, with a significance threshold set at a p-value of less than 0.05. Finally, we analyzed the correlation between the expression levels of prognostically significant lncRNAs and the TIIC types that showed significant distribution differences.

Establishment of m6A- or m5C-related lncRNA signature for survival prediction

Briefly, out of the survival-related lncRNAs, the optimal prognostic lncRNAs were selected by performing Lasso Cox regression analysis [22] using lars package (https://cran.r-project. org/web/packages/lars/index.html) in R. The optimal lambda value was determined by performing a 10-fold cross-validation. Based on LASSO coefficients and expression levels of the optimal signature lncRNAs, prognostic model was constructed using the following formula:

Risk score (RS) = ∑CoeflncRNAs×ExplncRNAs.

Where CoeflncRNAs represents the estimated LASSO coefficient of lncRNAs; ExplncRNAs represents expression level of lncRNAs.

With median risk score as the cutoff, samples were divided into two risk subgroups in the TCGA set. KM survival curves were plotted for different risk subgroups and compared using log-rank tests. Accuracy, sensitivity and specificity of the prognostic model were assessed using receiver operating characteristic (ROC) curve. Predictive performance of the prognostic score model was validated using GSE44001 dataset.

Identification of independent prognostic factors and nomogram construction

Clinical factors of patients in the TCGA set were subjected to uni-variable and multi-variable Cox regression analysis to identify independent prognostic factors using survival package in R (log-rank p-value < 0.05). Nomogram [23] that could provide a predicted probability was generated using rms package (version 5.1-2) in R (https://cran.r-project.org/web/packages/rms/index.html).

Drug sensitivity analysis

Based on gene expression profiles of CC samples in the TCGA set and chemotherapeutic agents downloaded from Genomics of Drug Sensitivity in Cancer (GDSC) [24] database (https://www.cancerrxgene.org/), sensitivity of CESC patients in high-risk and low-risk groups to chemotherapeutic drugs was assessed by predicting half-maximal inhibitory concentration (IC50) using pRRophetic [25] package in R (http://127.0.0.1:22402/doc/html/Search?objects=1&port=22402). IC50 values were compared between different risk groups using Students’ t test to assess differential therapeutics effects of chemotherapeutic drugs.

Function annotation analyses

Differentially expressed genes (DEGs) were screened between different risk groups of the TCGA set using limma [26] package in R (https://bioconductor.org/packages/release/bioc/html/limma.html) with FDR < 0.05 and |log2FC|>0.5 as the significance cutoff. Gene ontology (GO) function and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were conducted using The Database for Annotation, Visualization and Integrated Discovery (DAVID) tool (https://david.ncifcrf.gov/, version 6.8) with the strict cutoff at FDR < 0.05.

Results

A co-expression network was constructed including 23 m6A genes, 15 m5C genes and 343 lncRNAs

As mentioned in Method section, we recognized 23 m6A-related genes and 15 m5C-related genes, and extracted 343 common lncRNAs from both the TCGA set and GSE44001 dataset. By performing Pearson correlation analysis, we obtained 202 m6A genes-lncRNAs pairs and 264 m5C genes-lncRNAs pairs, which shared 62 common lncRNAs. With these m6A genes, m5C genes and 62 m6A and m5C-related lncRNAs, a co-expression network was constructed (Figure S1). Out of the 62 m6A and m5C-related lncRNAs, 9 lncRNAs were significantly related to OS time in uni-variate Cox regression analysis (Fig. 1A).

Fig. 1
figure 1

Identification of two subtypes with signficantly different survival time. (A) Forest plot of 9 survival-related lncRNAs in uni-variable Cox regression analysis; (B) samples are classified into two subtypes (subtype 1 and 2) by Consensus Clustering Analysis; (C) Kaplan- Meier curvers of two subtypes

Full size image

Two molecular subtypes were obtained based on nine survival-related lncRNAs

Based on expression levels of the 9 survival-related lncRNAs in the TCGA set, samples were categorized into subtype 1 (N = 139) and subtype 2 (N = 152) through performing Consensus Clustering Analysis (Fig. 1B). KM survival curves showed that subtype 1 showed significantly longer OS time than subtype 2 (p = 0.0011, HR = 2.287 (1.373–3.810), Fig. 1C).

Characterization of tumor infiltrating immune cells in different subtypes

In order to investigate the characteristics of TME and TIICs of the two molecular subtypes, we analyzed and compared fractions of different TIICs between them. As shown in Fig. 2A, subtype 1 had significantly lower levels of activated CD4 + T cells, regulatory T cells, myeloid derived suppressor cells, memory B cells, natural killer T cells and neutrophils, and higher levels of activated B cells, plasmacytoid dendritic cells, CD56 + bright natural killer cells, immature dendritic cells, and monocytes (p-value < 0.05). Moreover, PCCs between the 9 survival-related lncRNAs and 11 immune cells were shown in a heatmap (Fig. 2B). Fractions of activated CD4 + T cells, and neutrophils were negatively correlated with expression levels of all nine lncRNAs. On the contrary, fractions of plasmacytoid dendritic cells, and monocytes were positively correlated with expression levels of all nine lncRNAs.

Fig. 2
figure 2

Characterization of tumor infiltrating immune cells in different subtypes. (A) Comparative analysis of fractions of tumor infiltrating immune cells between two subtypes; (B) Pearson correlation coefficients of 9 survival-related lncRNAs with 11 infiltrating immune cells

Full size image

A prognostic model based on six m6A and m5C-related lncRNAs was constructed and successfully validated in an independent set

LASSO Cox regression was performed based on the aforementioned nine survival-related lncRNAs. As a result, six optimal lncRNAs were obtained, consisting of ASB16-AS1, DANCR, FBXL19-AS1, LINC01089, ST7-AS1, and ZNF213-AS1 (Fig. 3A-B). Based on LASSO coefficients and expression levels of the six optimal lncRNAs in the TCGA set, risk score was calculated for each sample using the following formula:

Fig. 3
figure 3

LASSO Cox regression analysis. (A) Performance indicators of LASSO regression analysis (B) LASSO coefficients of 6 signature lncRNAs. (C) Kaplan-Meier curves for patients stratified by each signature lncRNA. Patients are divided into high-expression group and low-expression groups according to median expression level of a signature lncRNA (ASB16-AS1, DANCR, LINC01089, ST7-AS1, FBXL19-AS1 or ZNF213-AS1)

Full size image

RS = (-0.0587601)*ExpLINC01089 + (-0.03925616)*Exp DANCR +(-0.03465494)*Exp ST7−AS1+ (-0.02335107)*Exp ASB16−AS1 +(-0.02032213)*Exp ZNF213−AS1+(0.11831476)*Exp FBXL19−AS1.

With regard to ASB16-AS1 (p = 0.0053), DANCR (p = 0.0077), LINC01089 (p = 0.013), ST7-AS1 (p = 0.025), or ZNF213-AS1 (p = 0.035), high-expressed patients had significantly longer OS time compared to low-expressed patients (Fig. 3C). Conversely, patients with low-expressed FBXL19-AS1 had longer OS time in comparison with patients with high-expressed FBXL19-AS1 (p = 0.019, Fig. 3C). In the TCGA set, with median risk score (-0.2913) as cutoff, tumor samples were separated into a high-risk group and a low-risk group (Fig. 4A-B). As depicted in Fig. 4C-D, OS time was significantly longer in low-risk patients compared to high-risk patients (p-value = 3.494e-06, HR = 3.310 (1.937–5.565)), with an AUC of 0.872 (0.841, 0.873). Similarly, the risk score model was applied to GSE44001 dataset. As shown in Fig. 4E-H, patients were dichotomized into high-risk and low-risk patients with significantly different OS time (p-value = 0.023, HR = 2.145 (1.094–4.209)) and an AUC of 0.795 (0.873, 0.705).

Fig. 4
figure 4

Performance of the risk model in TCGA set and validation set. Risk score distribution (A), relationship of risk score and survival time (B), KM curves of high-risk and low risk patients (C) and ROC curve (D) of patients in high-risk and low-risk subgroups in the TCGA set. Risk score distribution (E), relationship of risk score and survival time (F), KM curves of high-risk and low risk patients (G) and ROC curve (H) of patients in high-risk and low-risk subgroups in the validation set (GSE44001)

Full size image

Sankey diagram was used to analyze the relationships of the above-mentioned two subtypes, two risk subgroups and vital status of CC patients in the TCGA set. The majority of subtype 1 samples were classified into the low-risk group and had better prognosis, while the majority of subtype 2 samples were classified into the high-risk group and had poor prognosis (Fig. 5A). Risk scores were significantly decreased in subtype 1 samples compared to subtype 2 samples (Fig. 5B, p-value < 2.22e-16). These results reveal that the subtyping results based on the 9 survival-related lncRNAs were consistent with the risk stratification results based on the six-m6A and m5C-related lncRNAs signature.

Fig. 5
figure 5

Associations of molecular subtypes with risk score and nomogram construction. (A) Sankey diagram of molecular subtypes, risk score status and vital status of patients; (B) The risk scores of two subtypes were significantly different; (C) a nomogram for predicting 1-, 3- and 5- year survival probability in the TCGA set; (D) Calibration plots of the nomogram; (E) ROC curves for predicting 1-, 3- and 5- year survival probability

Full size image

A nomogram was constructed integrating risk score with prognostic factors

Uni-variable and multi-variable Cox regression analysis was carried out for clinical factors of patients in the TCGA set. As shown in Table 1, Pathologic N (p-value = 2.19E-03, 95%CI = 1.409–5.596), Pathologic T (p-value = 2.31E-05, 95%CI = 1.371–2.452), Pathologic Stage (p-value = 2.86E-04, 95%CI = 1.195–1.852), Tobacco smoking history (p-value = 6.30E-03, 95%CI = 1.081–1.904) and RS model (p-value = 3.49E-06, 95%CI = 1.937–5.656) were significantly associated with OS time of CC patients. Furthermore, RS model was an independent prognostic factor (p-value = 4.34E-02, 95%CI = 1.027–5.861).

Table 1 Uni-variable and multi-variable Cox regression analysis of clinical factors
Full size table

Combining risk score model with pathologic N, pathologic T, pathologic Stage, and tobacco smoking history, a nomogram was built to predict survival probabilities in CC patients (Fig. 5C). Calibration curves showed that 1-year C-index was 0.771, with a p-value of 1.258e-02; 3-year C-index = 0.743, with a p-value of 2.295e-05; 5-year C-index = 0.772, with a p-value of 6.036e-09 (Fig. 5D). Moreover, in ROC curves, 1-year, 3-year and 5-year AUC was 0.856 (0.750, 0.875), 0.875 (0.906, 0.833), 0.869 (0.900, 0.762), separately (Fig. 5E). These results collectively suggest that the nomogram model has high accuracy in predicting 1-, 3-, 5-year survival in CC patients.

Differential sensitivity of high-risk and low-risk patients to chemotherapeutic agents

We compared sensitivities of high-risk and low-risk patients to 21 chemotherapeutic agents (Axitinib, AZD6482, AZD7762, AZD8055, Bortezomib, Camptothecin, Cisplatin, Cytarabine, Dasatinib, Erlotinib, Gefitinib, Gemcitabine, GSK269962A, Lapatinib, Nilotinib, Paclitaxel, Rapamycin, Sorafenib, Vinblastine, Vinorelbine, and Vorinostat) in the TCGA set. Low-risk patients were more sensitive to AZD8055 (p-value = 0.0012), Rapamycin (p-value = 6.2e-07) and Vorinostat (p-value = 0.046), whereas high-risk patients were more sensitive to Dasatinib (p-value = 0.00021, Fig. 6A).

Fig. 6
figure 6

Drug sensitivity and functional enrichment analysis. (A) IC50 values of AZD8055, Rapamycin, Vorinostat and Dasatinib between high-risk and low-risk patients in the TCGA set. (B) Top 10 GO biological processes and top 10 KEGG signaling pathways

Full size image

Functional annotation of the DEGs between two risk groups in the TCGA set

A total of 624 DEGs (FDR < 0.05 and |log2FC| > 0.5) were identified between the two risk subgroups of the TCGA set. These genes were significantly enriched in 25 GO biological processes and 10 KEGG signaling pathways, such as inflammatory response, positive regulation of interleukin-6 production, Calcium signaling pathway, MAPK signaling pathway, and IL-17 signaling pathway. Top 10 BP terms and KEGG signaling pathways are showcased in Fig. 6B.

Discussion

LncRNAs play critical regulatory roles in progression, metastasis, and prognosis of CC [27]. As two well-studied forms of methylation modifications, m6A and m5C modifications that contribute to cancer occurrence and development are gaining increasing attention [28]. Our study underlined the regulatory mechanisms of m6A and m5C modifications-related lncRNAs in CC. We identified 62 m6A- and m5C-related lncRNAs in CC samples, and constructed a co-expression network with these lncRNAs, 23 m6A-related genes and 15 m5C-related genes. Based on nine survival-related lncRNAs we obtained two molecular subtypes with distinctive TME characteristics. Noticeably, we constructed a prognostic model based on six m6A- and m5C-related signature lncRNAs, which separated patients into high-risk and low-risk groups with significantly different OS time. Robustness of the prognostic model was successfully validated in an independent cohort. Low-risk patients showed higher sensitivity to AZD8055, Rapamycin and Vorinostat, whereas high-risk patients showed higher sensitivity to Dasatinib. Furthermore, a nomogram based on risk score and prognostic clinical factors was built and showed robust and accuracy performance in predicting prognosis in CC patients. These results suggest that the prognostic model based on six m6A- and m5C-related lncRNA is a promising tool for prognosis stratification of CC patients and provide new hints to direct individualized therapies.

Previous studies suggest that m6A RNA methylation regulators play critical roles.

in the malignant progression of CC, and have prognostic implications [29, 30]. Moreover, m6A-related lncRNAs may be promising prognostic biomarkers for CC [14]. The m6A- and m5C-related lncRNAs signature obtained by our study contained lncRNA ASB16-AS1, DANCR, FBXL19-AS1, LINC01089, ST7-AS1, and ZNF213- AS1. LncRNA ASB16-AS1 is a promising pan-cancer prognostic biomarker, with an association with immune infiltration [31]. There is evidence that lncRNA ASB16-AS1 expression is up-regulated in CC, strengthening cell proliferation, and migration via Wnt/β-catenin signal pathway [32]. LncRNA DANCR has prognostic significance in human cancers, and correlates with worse prognosis [33]. Substantial evidences have supported an oncogenic role of lncRNA DANCR in CC, which promotes CC proliferation, metastasis and progression [34, 35]. Emerging studies show that lncRNA FBXL19-AS1 promotes CC proliferation and metastasis [36, 37]. Besides, prognostic potential of lncRNA FBXL19-AS1 has been suggested for CC [38]. LncRNA LINC01089 exerts suppressive effects on the development of CC [39]. LncRNA ST7-AS1 overexpression exhibits positive correlations with shorter OS time, a higher frequency of lymph node metastasis and deeper cervical invasion [40]. Our study consistently found involvement of the 5 m6A- and m5C-related lncRNAs in the progression of CC. Elevated expression of lncRNA ZNF213-AS1 plays a role in differentiation and proliferation of acute myeloid leukemia and low-grade gliomas, with an insignificant association with poor prognostic outcome [41]. However, biological roles of lncRNA ZNF213-AS1 in CC remain undefined. Our study indicates that these m6A- and m5C-related lncRNAs may serve as prognostic biomarkers for CC.

Our study revealed two molecular subtypes based on express levels of prognostic lncRNAs. The majority of subtype 1 samples had low-risk scores and longer OS time, whereas the majority of subtype 2 samples had high-risk scores and shorter OS time. TME has been recognized as a principal component of CC tumorigenesis and development, and influences prognosis and treatment of patients [42, 43]. Plasmacytoid dendritic cells play an immunosuppressive role in TME and promote tumor growth [44]. Our study unveiled distinctive TME characteristics of the two subtypes. Subtype 1 had significantly higher levels of activated B cells, CD56 + bright natural killer cells, plasmacytoid dendritic cells, and lowers levels of myeloid derived suppressor cells. It implies that stronger anti-tumor immune function is an important contributor to better prognosis of subtype 1 compared to subtype 2. These findings indicate that m6A- and m5C-related lncRNAs could serve as effective prognostic biomarkers and predictors for clinical outcomes and immunotherapeutic responses in CC patients. Additionally, our study found that the low-risk patients were more sensitive to AZD8055, Rapamycin and Vorinostat, whereas high-risk patients were more sensitive to Dasatinib. It offers valuable data to facilitate individualized clinical treatments for CC patients. Our study also unraveled that various inflammation-related biological processes and pathways, such as inflammatory response, positive regulation of interleukin-6 production, and IL-17 signaling pathway, as well as Calcium signaling pathway and MAPK signaling pathway might participate in the regulatory mechanisms of m6A and m5C modifications-related lncRNAs in CC in concordance with previous findings [45,46,47].

This study has several limitations. First, since the data were analyzed from the TCGA and GEO databases, there is a lack of validation through experiments. Second, the prognostic value of the RiskScore model requires validation with additional external datasets before it can be applied clinically. Finally, the potential molecular mechanisms of these m6A/m5C-related lncRNAs in cervical cancer remain unclear, and we plan to conduct in vitro or in vivo experiments in future research to validate our findings.

Conclusion

Taken together, our study comprehensively analyzed interactions of m6A and m5C modifications with lncRNAs in CC samples, and established a prognostic model based on a signature of six m6A- and m5C-related lncRNAs, which could be used to predict outcome of patients. These signature lncRNAs had close associations with TIICs and cancer prognosis, highlighting great promises as prognostic biomarkers and therapeutic targets for CC. Our study expands the knowledge concerning involvement of m6A- and m5C-related lncRNAs in CC pathogenesis and has clinical implications for risk stratification and personalized therapeutics of patients.

Data availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Abbreviations

(m6A):

N6 methyladenosine

(m5C):

5 methylcytosine

(lncRNAs):

long non-coding RNAs

(CC):

cervical cancer

(TIICs):

tumor infiltrating immune cells

(OS):

overall survival

(CESC):

cervical squamous cell carcinoma

(HPV):

Human papilloma virus

(mRNAs):

messenger RNA

(ncRNAs):

non-coding RNAs

(TCGA):

The Cancer Genome Atlas

(TIICs):

tumor infiltrating immune cells

(GEO):

gene expression omnibus

(KM):

Kaplan- Meier

(ROC):

receiver operating characteristic

(GDSC):

Genomics of Drug Sensitivity in Cancer

(DEGs):

Differentially expressed genes

(GO):

Gene ontology

(KEGG):

Kyoto Encyclopedia of Genes and Genomes

(DAVID):

Database for Annotation, Visualization and Integrated Discovery

References

  1. Sung H, Ferlay J, Siegel RLL, Soerjomataram M, Jemal I, Bray A. Global Cancer statistics 2020: GLOBOCAN estimates of incidence and Mortality Worldwide for 36 cancers in 185 countries. Cancer J Clin. 2021;71(3):209–49.

    Article  Google Scholar 

  2. Pimple SA, Mishra GA. Global strategies for cervical cancer prevention and screening. Minerva Ginecol. 2019;71(4):313–20.

    Article  PubMed  Google Scholar 

  3. Shami S, Coombs J. Cervical cancer screening guidelines: an update. JAAPA: Official J Am Acad Physician Assistants. 2021;34(9):21–4.

    Article  Google Scholar 

  4. Hill EK. Updates in Cervical Cancer Treatment. Clin Obstet Gynecol. 2020;63(1):3–11.

    Article  PubMed  Google Scholar 

  5. Ward ZJ, Scott AM, Hricak H, Atun R. Global costs, health benefits, and economic benefits of scaling up treatment and imaging modalities for survival of 11 cancers: a simulation-based analysis. Lancet Oncol. 2021;22(3):341–50.

    Article  PubMed  PubMed Central  Google Scholar 

  6. Lin W, Zhou Q, Wang CQ, Zhu L, Bi C, Zhang S, Wang X, Jin H. LncRNAs regulate metabolism in cancer. Int J Biol Sci. 2020;16(7):1194–206.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  7. Taniue K, Akimitsu N. The functions and unique features of LncRNAs in Cancer Development and Tumorigenesis. Int J Mol Sci. 2021;22(2):1–20.

    Article  Google Scholar 

  8. Yang Q, Al-Hendy A. The Regulatory functions and the mechanisms of Long non-coding RNAs in Cervical Cancer. Cells 2022, 11(7).

  9. He J, Huang B, Zhang K, Liu M, Xu T. Long non-coding RNA in cervical cancer: from biology to therapeutic opportunity. Biomed Pharmacotherapy = Biomedecine Pharmacotherapie. 2020;127:1–10.

    Google Scholar 

  10. Wang T, Kong S, Tao M, Ju S. The potential role of RNA N6-methyladenosine in Cancer progression. Mol Cancer. 2020;19(1):1–18.

    Article  CAS  Google Scholar 

  11. Wang Q, Guo X, Li L, Gao Z, Su X, Ji M, Liu J. N(6)-methyladenosine METTL3 promotes cervical cancer tumorigenesis and Warburg effect through YTHDF1/HK2 modification. Cell Death Dis. 2020;11(10):1–10.

    Article  PubMed  PubMed Central  Google Scholar 

  12. Wang H, Luo Q, Kang J, Wei Q, Yang Y, Yang D, Liu X, Liu T, Yi P. YTHDF1 aggravates the progression of Cervical Cancer through m(6)A-Mediated Up-Regulation of RANBP2. Front Oncol. 2021;11:1–13.

    Google Scholar 

  13. Condic M, Ralser DJ, Klümper N, Ellinger J, Qureischi M, Egger EK, Kristiansen G, Mustea A, Thiesler T. Comprehensive analysis of N6-Methyladenosine (m6A) writers, erasers, and readers in Cervical Cancer. Int J Mol Sci. 2022;23(13):1–11.

    Article  Google Scholar 

  14. Zhang H, Kong W, Zhao X, Han C, Liu T, Li J, Song D. N6-Methyladenosine-related lncRNAs as potential biomarkers for predicting prognoses and immune responses in patients with cervical cancer. BMC Genomic data. 2022;23(1):1–13.

    Article  PubMed  PubMed Central  Google Scholar 

  15. Nombela P, Miguel-López B, Blanco S. The role of m(6)A, m(5)C and Ψ RNA modifications in cancer: novel therapeutic opportunities. Mol Cancer. 2021;20(1):1–30.

    Article  Google Scholar 

  16. Guo G, Pan K, Fang S, Ye L, Tong X, Wang Z, Xue X, Zhang H. Advances in mRNA 5-methylcytosine modifications: detection, effectors, biological functions, and clinical relevance. Mol Therapy Nucleic Acids. 2021;26:575–93.

    Article  CAS  Google Scholar 

  17. Wang L, Zhang J, Su Y, Maimaitiyiming Y, Yang S, Shen Z, Lin S, Shen S, Zhan G, Wang F, et al. Distinct roles of m(5)C RNA methyltransferase NSUN2 in major gynecologic cancers. Front Oncol. 2022;12:1–15.

    Google Scholar 

  18. Zuo Z, Xiong J, Zeng C, Jiang Y, Xiong K, Tao H, Guo Y. Exploration of a robust and prognostic Immune Related Gene signature for cervical squamous cell carcinoma. Front Mol Biosci. 2021;8:1–15.

    Article  Google Scholar 

  19. Song W, Ren J, Xiang R, Yuan W, Fu T. Cross-talk between m(6)A- and m(5)C-Related lncRNAs to Construct a Novel signature and predict the Immune Landscape of Colorectal Cancer patients. Front Immunol. 2022;13:1–16.

    Google Scholar 

  20. Doncheva NT, Morris JHG, Jensen J. Cytoscape StringApp: Network Analysis and Visualization of Proteomics Data. J Proteome Res. 2019;18(2):623–32.

    Article  PubMed  CAS  Google Scholar 

  21. Zhang X, Ren L, Yan X, Shan Y, Liu L, Zhou J, Kuang Q, Li M, Long H, Lai W. Identification of immune-related lncRNAs in periodontitis reveals regulation network of gene-lncRNA-pathway-immunocyte. Int Immunopharmacol. 2020;84:1–12.

    Article  Google Scholar 

  22. Goeman JJ. L1 penalized estimation in the Cox proportional hazards model. Biometrical J Biometrische Z. 2010;52(1):70–84.

    Article  Google Scholar 

  23. Balachandran VP, Gonen M, Smith JJ, DeMatteo RP. Nomograms in oncology: more than meets the eye. Lancet Oncol. 2015;16(4):e173–180.

    Article  PubMed  PubMed Central  Google Scholar 

  24. Yang W, Soares J, Greninger P, Edelman EJ, Lightfoot H, Forbes S, Bindal N, Beare D, Smith JA, Thompson IR, et al. Genomics of Drug Sensitivity in Cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res. 2013;41(Database issue):D955–961.

    PubMed  CAS  Google Scholar 

  25. Zhu Y, Meng X, Ruan X, Lu X, Yan F, Wang F. Characterization of Neoantigen Load Subgroups in gynecologic and breast cancers. Front Bioeng Biotechnol. 2020;8:1–11.

    Article  Google Scholar 

  26. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W. Smyth GK: limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):1–13.

    Article  Google Scholar 

  27. Aalijahan H, Ghorbian S. Long non-coding RNAs and cervical cancer. Exp Mol Pathol. 2019;106:7–16.

    Article  PubMed  CAS  Google Scholar 

  28. Zhu H, Zhu H, Tian M, Wang D, He J, Xu T. DNA methylation and hydroxymethylation in Cervical Cancer: diagnosis, prognosis and treatment. Front Genet. 2020;11:1–12.

    Google Scholar 

  29. Pan J, Xu L, Pan H. Development and validation of an m6A RNA methylation Regulator-based signature for prognostic prediction in cervical squamous cell carcinoma. Front Oncol. 2020;10:1–8.

    Article  Google Scholar 

  30. Ma X, Li Y, Wen J, Zhao Y. m6A RNA methylation regulators contribute to malignant development and have a clinical prognostic effect on cervical cancer. Am J Translational Res. 2020;12(12):8137–46.

    CAS  Google Scholar 

  31. Wu L, Liao W, Wang X, Zhao Y, Pang J, Chen Y, Yang H, He Y. Expression, prognosis value, and immune infiltration of lncRNA ASB16-AS1 identified by pan-cancer analysis. Bioengineered. 2021;12(2):10302–18.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  32. Liu W, Zhuang R, Feng S, Bai X, Jia Z, Kapora E, Tan W. Long non-coding RNA ASB16-AS1 enhances cell proliferation, migration and invasion via functioning as a ceRNA through miR-1305/Wnt/β-catenin axis in cervical cancer. Biomed Pharmacotherapy = Biomedecine Pharmacotherapie. 2020;125:1–10.

    Google Scholar 

  33. Liu W, Wang QP, Guo J. Prognostic significance of lncRNA DANCR expression in human cancers: a systematic review and meta-analysis. Biosci Rep. 2021;41(8):1–9.

    Article  Google Scholar 

  34. Hu C, Han Y, Zhu G, Li G, Wu X. Krüppel-like factor 5-induced overexpression of long non-coding RNA DANCR promotes the progression of cervical cancer via repressing microRNA-145-3p to target ZEB1. Cell Cycle (Georgetown Tex). 2021;20(14):1441–54.

    Article  PubMed  CAS  Google Scholar 

  35. Cao L, Jin H, Zheng Y, Mao Y, Fu Z, Li X, Dong L. DANCR-mediated microRNA-665 regulates proliferation and metastasis of cervical cancer through the ERK/SMAD pathway. Cancer Sci. 2019;110(3):913–25.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  36. Huang X, Shi H, Shi X, Jiang X. LncRNA FBXL19-AS1 promotes proliferation and metastasis of cervical cancer through upregulating COL1A1 as a sponge of miR-193a-5p. J Biol Res (Thessalonike Greece). 2021;28(1):1–14.

    Google Scholar 

  37. Wan S, Ni G, Ding J, Huang Y. Long noncoding RNA FBXL19-AS1 expedites cell growth, Migration and Invasion in Cervical Cancer by miR-193a-5p/PIN1 signaling. Cancer Manage Res. 2020;12:9741–52.

    Article  CAS  Google Scholar 

  38. Dai S, Yao D. An immune-associated ten-long noncoding RNA signature for predicting overall survival in cervical cancer. Translational cancer Res. 2021;10(12):5295–306.

    Article  CAS  Google Scholar 

  39. Li S, Han Y, Liang X, Zhao M. LINC01089 inhibits the progression of cervical cancer via inhibiting miR-27a-3p and increasing BTG2. J Gene Med. 2021;23(1):1–10.

    Article  Google Scholar 

  40. Qi H, Lu L, Wang L. Long noncoding RNA ST7-AS1 Upregulates TRPM7 expression by sponging microRNA-543 to promote Cervical Cancer Progression. OncoTargets Therapy. 2020;13:7257–69.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  41. Ramilowski JA, Yip CW, Agrawal S, Chang JC, Ciani Y, Kulakovskiy IV, Mendez M, Ooi JLC, Ouyang JF, Parkinson N, et al. Functional annotation of human long noncoding RNAs via molecular phenotyping. Genome Res. 2020;30(7):1060–72.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  42. Xu F, Shen J, Xu S. Integrated Bioinformatical Analysis identifies GIMAP4 as an Immune-related prognostic Biomarker Associated with Remodeling in Cervical Cancer Tumor Microenvironment. Front cell Dev Biology. 2021;9:1–13.

    CAS  Google Scholar 

  43. Chen Q, Qiu B, Zeng X, Hu L, Huang D, Chen K, Qiu X. Identification of a tumor microenvironment-related gene signature to improve the prediction of cervical cancer prognosis. Cancer Cell Int. 2021;21(1):1–16.

    PubMed  PubMed Central  CAS  Google Scholar 

  44. Zhou B, Lawrence T, Liang Y. The role of Plasmacytoid dendritic cells in cancers. Front Immunol. 2021;12:1–10.

    Article  Google Scholar 

  45. Liu Y, Li L, Li Y, Zhao X. Research Progress on Tumor-Associated macrophages and inflammation in Cervical Cancer. Biomed Res Int. 2020;2020:1–6.

    CAS  Google Scholar 

  46. Zhou B, Li T, Xie R, Zhou J, Liu J, Luo Y, Zhang X. CircFAT1 facilitates cervical cancer malignant progression by regulating ERK1/2 and p38 MAPK pathway through miR-409-3p/CDK8 axis. Drug Dev Res. 2021;82(8):1131–43.

    Article  PubMed  CAS  Google Scholar 

  47. Xu J, Pan Q, Ju W. Ras inhibition by zoledronic acid effectively sensitizes cervical cancer to chemotherapy. Anticancer Drugs. 2019;30(8):821–7.

    Article  PubMed  CAS  Google Scholar 

Download references

Acknowledgements

None.

Funding

None.

Author information

Authors and Affiliations

  1. Department of Laboratory Medicine, The Sixth Affiliated Hospital of Harbin Medical University, 998 Aiying Street, Songbei District, Harbin, 150028, Heilongjiang Province, China

    Jie Gao, Xiuling Zhang, Anqi Xu & HaiYan Gao

  2. Department of Rheumatism and Immunology, The 2nd Affiliated Hospital of Harbin Medical University, 246 Xuefu Road, Nangang District, Harbin, 150001, Heilongjiang Province, China

    Wei Li

Authors
  1. Jie Gao
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Xiuling Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Anqi Xu
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Wei Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. HaiYan Gao
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

Jie Gao carried out the Conception and design of the research, Xiuling Zhang participated in the acquisition of data. Anqi Xu and Wei Li carried out the Analysis and interpretation of data. HaiYan Gao participated in the design of the study and performed the statistical analysis. Jie Gao drafted the manuscript and revision of manuscript for important intellectual content. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Wei Li or HaiYan Gao.

Ethics declarations

Ethical approval and informed consent statement

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic supplementary material

Below is the link to the electronic supplementary material.

Supplementary Material 1

Supplementary Material 2

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Gao, J., Zhang, X., Xu, A. et al. Integrated analysis of N6-methyladenosine- and 5-methylcytosine-related long non-coding RNAs for predicting prognosis in cervical cancer. Hereditas 161, 34 (2024). https://doi.org/10.1186/s41065-024-00336-w

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s41065-024-00336-w

Keywords

Hereditas

ISSN: 1601-5223

Contact us