Gene signatures, immune infiltration, and drug sensitivity based on a comprehensive analysis of m6a RNA methylation regulators in cervical cancer
Journal of Translational Medicine volume 20, Article number: 385 (2022)
Cervical cancer is the fourth most common cancer in women. N6-dimethyladenosine (m6A) mRNA methylation is closely associated with cervical cancer.
Using TCGA database, we studied the expression and mutation of m6A-related genes in cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC) and obtained genetic characteristics based on an m6A risk model and prognostic value of m6A. We studied the effects of the m6A risk score on immune features and genomic changes of patients with CESC, evaluated the sensitivity of patients with CESC to different small-molecule drugs based on the m6A risk score, and established a clinical prediction model.
Ten m6A-related genes were differentially expressed between CESC and normal tissues. High-risk patients had a low overall survival (OS) and significantly low immune scores but showed no significantly altered stromal scores. The tumor mutation burden (TMB) and tumor neoantigen levels significantly differed between the high- and low-risk groups. In the high-risk group, copy number variation (CNV) changes mainly led to gene amplification, while in the low-risk group, CNV changes primarily manifested as gene copy number deletions. ZC3H13 expression was low in CESC tissues. ZC3H13 knockdown promoted CESC cell proliferation, migration, and invasion, reducing the RNA methylation levels. Rapamycin suppressed the CESC cell proliferation, migration, and invasion abilities, increasing the m6A levels.
m6A mRNA methylation is closely related to the occurrence, development, immune invasion, drug sensitivity, and prognosis of cervical cancer. The prognostic m6A feature model of m6A signature genes can accurately predict the OS of patients with CESC. Drugs targeting factors regulating m6A mRNA methylation might offer a good prospect for treating cervical cancer.
Cervical cancer is the fourth most common cancer in women [1, 2]. More than 570,000 new cases of cervical cancer are diagnosed worldwide each year, with approximately 311,000 deaths . Presently, at least 170 RNA modifications have been discovered in cervical cancer, for instance, N6-dimethyladenosine (m6A), inosine, pseudouridine, 5-methylcytidine (m5C), 5-hydroxymethylcytidine (hm5C), and N1-methyladenosine (m1A). Furthermore, most RNA species contain at least one chemical modification and are widely associated with physiopathological processes [4,5,6,7]. m6A is one of the most abundant mRNA modifications [8, 9]. m6A modifies approximately 0.1–0.4% of all adenosines in RNA, accounting for ~ 50% of the total methylated ribonucleotides . m6A is a dynamic and reversible epigenetic modification. The cellular m6A status is regulated by groups of proteins called m6A methyltransferases (“writers”), m6A demethylases (“erasers”), and m6A-binding proteins (“readers”), which add, remove, or recognize m6A-modified sites, respectively, modulating the stability, splicing, intracellular distribution, and translational changes of mRNA while affecting certain biological processes [11,12,13,14]. As a potential tumor biomarker, m6A plays a key role in various biological processes and malignancies . Changes in m6A-modifying enzyme levels affect the expression of downstream oncogenes or tumor suppressor genes by altering mRNA methylation . Furthermore, epigenetic markers potentially serve as diagnostic, prognostic, and predictive biomarkers and might be used as novel targets for cancer precision medicine [17,18,19,20].
Here, we studied the expression and mutation of m6A-related genes in patients with cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), obtained the gene characteristics of patients with CESC based on an m6A-risk model, and constructed a prognostic m6A feature model based on m6A signature genes. The effects of the m6A risk score on biological function, immune characteristics, and genomic changes of patients with CESC were analyzed. The sensitivity of patients with CESC to different small molecule drugs was evaluated based on the m6A risk score, and a clinical prediction model was constructed. Subsequently, we studied the expression and function of ZC3H13 in cervical cancer tissues according to the screening results. We selected rapamycin to study the effects of the screened drug on m6A and cervical cancer phenotype. Through this research, we attempted to study the gene signatures, immune infiltration, and drug sensitivity based on a comprehensive analysis of m6A RNA methylation regulators in cervical cancer.
Materials and methods
Data were downloaded from the TCGA GDC website (https://portal.gdc.cancer.gov/), and the cervical squamous carcinoma and adenocarcinoma (cervical squamous cell carcinoma and endocervical adenocarcinoma, CESC) expression profile data (FPKM results) of gene expression sequencing in patients were obtained. The FPKM values were converted into the TPM values as follows:
Clinical data of corresponding patients, including age, TNM stage, and survival prognosis, were downloaded. After excluding patients with CESC with missing clinical information, data regarding 279 tumor tissues and three paracancer tissues were obtained. The copy number variation (CNV) of somatic cells of patients with CESC was downloaded. RCircos package was used to map the gene CNV in 23 pairs of chromosomes . “Masked mutation” was selected as the somatic mutation data, and R’s MAfTools package was used to visualize somatic mutations. The tumor mutation burden (TMB) of each patient was obtained . The sequencing results of 19 cervical tissues were downloaded from the GTEx database and converted into TPM values.
In addition, gene expression data of the GSE52903 chip, clinicopathological features, and prognostic information of patients were downloaded from samples in the GEO database, including 55 tumor tissue and 17 normal cervical tissue samples . Among them, all chip data samples were from Homo Sapiens, and the platform was mainly based on the GPL6244 [Hugene-1_0-ST] Affymetrix Human Gene 1.0 ST Array [Transcript (Gene) version]. Next, TCGA data, GTEx data, and GEO chip data were combined, and the limma and sva packages of R were used for homogenization and batch removal [24, 25]. The total process is composed of two steps. The first is data consolidation; only common genes and corresponding expression values are retained in the three data sets. The second step involves removing the batch effects. The first step adopts the merge function in R to merge the three data sets. The second step uses the combat function in the SVA package to remove the batches. The demographic and clinical characteristics of the cervical cancer patients are shown in Additional file 1 materials.
Risk model construction based on m6A-related genes
To analyze the expression of m6A-related genes in cervical cancer, we first analyzed the differential expression of m6A-related genes in cervical cancer and normal tissues, the correlation of gene expression, and its impact on the prognosis of patients with CESC. Subsequently, expression data of patients with CESC were pooled using TCGA-CESC and GEO data. m6A-related genes were included in the model, and the least absolute shrinkage and selection operator (LASSO) algorithm was used to analyze the dimension reduction and obtain characteristic genes associated with prognosis. A risk score formula was established using the penalty coefficient weighted normalized gene expression values obtained by LASSO Cox analysis. The patients were divided into high- and low-risk groups based on the median risk score.
Gene set enrichment analysis (GSEA)
Gene Ontology (GO) analysis is a common method for large-scale functional enrichment studies that includes three categories: biological process (BP), molecular function (MF), and cellular component (CC). KEGG is a widely used database that stores information regarding genomes, biological pathways, diseases, and drugs. R's clusterProfiler package was used for GO annotation analysis and KEGG pathway enrichment analysis of signature genes. A value of FDR < 0.05 was considered statistically significant .
To study the differences in BPs between different subgroups, we performed GSEA based on the gene expression profile dataset of patients with CESC. GSEA is a computational method to analyze whether a specific gene set has statistically significant differences between two biological states and is usually used to estimate changes in pathways and BP activity among samples . The “h.all.v7.2.symbols.gmt” gene set was downloaded from MSigDB for GSEA analysis . An adjusted P-value of less than 0.05 was considered statistically significant.
Assessment of patients’ biological characteristics between risk groups
We further analyzed the correlation between different subgroups and some biologically related processes using the GSVA method . Mariathasan et al. constructed a gene set for storing genes related to some BPs, including (1) immune checkpoint; (2) antigen processing; (3) characteristics of CD8+ T cells; (4) epithelial-mesenchymal transformation (EMT) markers, including EMT1, EMT2, and EMT3; (5) angiogenesis; (7) characteristics of TGF-β response in pan-FTBRs; (8) WNT characteristics; (9) DNA damage repair; (10) mismatch repair; (11) nucleotide excision repair; (12) DNA replication; and (13) antigen handling and presentation [28,29,30]. Gene sets corresponding to different biological characteristics were downloaded to calculate the corresponding enrichment scores of patients and to compare the differences between two groups.
Analysis of m.6A-related differentially expressed genes (DEGs)
To identify genes associated with the m6A risk model, THE limma package of R was used to analyze DEGs between high and low subgroups in patients with CESC, and the DEGs with significant differences were defined as genes with absolute log (Fold change) > 0.5 and FDR < 0.05 . Hierarchical clustering was used to divide the tumors into different gene groups based on the Euclidean distance of differential gene expression and named them Geneclusters A, B, and C. Among them, R's ConsensusClusterPlus package was used to determine the number of clusters in the dataset and was repeated 1000 times to ensure the stability of classification . Meanwhile, based on the expression changes of specific genes, they were divided into signature-A and -B gene groups.
Calculation of dimensionality reduction and m6A score
First, unsupervised clustering was used to classify TCGA data of patients according to DEG values. m6A signature-A and -B gene groups were reduced in dimension according to the Boruta algorithm. The principal component PC1 was extracted by the principal component analysis (PCA) algorithm as the A score. Finally, we applied a method similar to the gene expression grade index to define each patient’s immune checkpoint inhibitor (ICI) score:
Identification and correlation analysis of tumor immune infiltrating cells
To further quantify the proportion of different immune cells in CESC samples, the CIBERSORT algorithm, and LM22 gene set were used to investigate the phenotypes of 22 human immune cells (including B cells, T cells, and natural killer cells) in the tumor microenvironment (TME) . Macrophages are highly sensitive and specific. CIBERSORT is a deconvolution algorithm that uses a set of reference gene expression values (with 547 characteristic genes) considered to be the minimum representative of each cell type. Based on these values, the proportion of cell types can be inferred in data from a large tumor sample with mixed cell types.
At the same time, R's ESTIMATE package was used to assess tumor immune activity . ESTIMATE an immune score for each tumor sample obtained by quantifying the immune activity (level of immune invasion) in the tumor sample based on its gene expression profile. Then, differences in the characteristics of immune infiltration in patients with CESC between the high- and low-risk groups were obtained.
Analysis of CNV
To analyze the changes in the copy number in different risk score groups of TCGA-patients with CESC, the TCGAbiolinks package of R was used to download the masked copy number segment data of patients. GISTIC 2.0 analysis of the downloaded CNV fragments was performed by GenePattern5. During analysis, except for a few parameters (e.g., the confidence is 0.99; Excluding X chromosomes prior to analysis), GISTIC 2.0 analysis was used with the default settings. Finally, R's MAfTools package was used to visually display the results of the GISTIC 2.0 analysis.
Sensitivity analysis of anticancer drugs
Genomics of Drug Sensitivity in Cancer (GDSC; https://www.cancerrxgene.org/), which is used for molecular cancer therapy and mutation, was used to explore public databases . R's pRRophetic package was used to download cell line gene mutation data and IC50 values of different anticancer drugs to analyze the correlation between patients with high- and low-risk scores and sensitivity to different anticancer drugs .
Building a clinical prediction model based on the m6A risk model
To prove that the risk score combined with the clinicopathological features can help in the personalized assessment of patient prognosis, univariate and multivariate Cox analyses were conducted to analyze the ability of the risk score combined with clinicopathological features to predict the overall survival (OS). Then, a nomogram was constructed by incorporating the risk scoring model and clinicopathological features into the model. To quantify the differentiated performance, Harrell’s Consistency index (C-index) was measured. A calibration curve was generated to evaluate the performance of the rosette by comparing the predicted values of the rosette with the actual observed survival rates.
Human tissue samples
Ten cervical cancer tissues and ten normal cervical tissues were obtained from patients with cervical cancer who underwent ovariectomy prior to chemotherapy and radiotherapy. The resected cervical cancer tissues and normal tissues were immediately stored at − 80 °C for further study. All patients signed informed consent forms. This study was performed in accordance with the Helsinki Declaration and approved by the ethics committee of the Second Affiliated Hospital of Zhengzhou University.
Cell lines and culture conditions
The human cervical cancer cell lines HeLa and SiHa were purchased from Procell Life Science & Technology Co., Ltd. (Wuhan, P.R. China). To prepare the complete growth medium, the cell culture media were supplemented with fetal bovine serum (Gibco, Grand Island, USA) at a final concentration of 10%.
Total RNAs of human cells were extracted using the Trizol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions and treated with RQ1 DNase (Promega, Madison, WI, USA) to remove DNA. The quality and quantity of purified RNA were determined by measuring the absorbance at 260 nm and 280 nm (A260 and A280, respectively) using a SmartSpec Plus Spectrophotometer (Bio-Rad Laboratories, Inc., Hercules, CA, USA). RNA integrity was further verified by electrophoresis using a 1.5% agarose gel. All RNA samples were stored at -80 °C for future analysis. Reverse transcription reactions were carried out using the ReverTra Ace qPCR RT Kit (TOYOBO Life Science, Shanghai, P.R. China), according to the manufacturer’s instructions.
Quantitative real-time PCR (qRT-PCR)
Expression levels of the ZC3H13 gene were detected by qRT-PCR. The human (species) ACTB gene was used as a control. Specific primers were designed based on cDNA sequences. Primer sequences of ZC3H13 were as follows: 5′- ACATTCATTAGGCTCTGGTGC -3′ forward), 5′—TTCTCCTCATCCTGTTGGTCC—3′ (reverse). qRT-PCR was performed on a Bio-Rad S1000 with Bestar SYBR GreenRT-PCR Master Mix (TOYOBO). PCR conditions consisted of denaturing at 95 °C for 1 min, and 40 cycles of denaturing at 95 °C for 15 s followed by annealing and extension at 60 °C for 30 s. The relative gene expression was calculated using the Livak and Schmittgen 2−ΔΔCt method (Livak and Schmittgen 2001), normalized with the reference gene actin. PCR amplifications were performed in triplicate for each sample.
The tissue slides (4-μm-thick sections) were initially treated for deparaffinization, rehydration, and antigen-retrieval using 3% H2O2. The sections were incubated with anti-ZC3H13 (Bioss, P.R. China) and then with horseradish peroxidase (HRP)-labeled IgG secondary antibodies (Beyotime Institute of Biotechnology, Shanghai, P.R. China). Fields from each slide were examined and photographed using a light microscope (Olympus, Japan).
Western blot analysis
Cells were lysed in a radioimmunoprecipitation assay buffer containing 1 nM phenylmethanesulfonyl fluoride, and protein concentrations were determined using a bicinchoninic acid (BCA) protein assay kit (Biosharp, Guangzhou, China) according to the manufacturer’s instructions. Protein samples (40 μg) of each group were boiled with 6 × sodium dodecyl sulfate loading buffer for 10 min before electrophoresis on 12% sodium dodecyl sulfate–polyacrylamide gel electrophoresis gels. The resolved proteins were electrotransferred onto polyvinylidene difluoride membranes (Millipore, Billerica, MA, USA) in a transferring buffer (25 mM Tris, 0.2 M glycine, and 25% methanol). After blocking with 5% skimmed milk, the membranes were incubated with anti-ZC3H13 (Bioss, P.R. China) and anti-β-actin antibodies (Bioss, P.R. China), followed by incubation with appropriate HRP-conjugated secondary antibodies (Abcam, England).
ZC3H13-siRNA (5′-GAAGACAUCUGCAGUAUCU-3′, antisense 5′-AGAUACUGCAGAUGUCUUC-3′) was synthesized by GeneCreate Bioengineering Co., Ltd. Gene Create (Wuhan, P.R. China). According to the manufacturer’s instructions, Lipofectamine 2000 (Invitrogen, USA) was used for transfecting siRNAs into HeLa and SiHa cells.
Cells were digested with trypsin and inoculated into 96-well plates. The culture plate was placed in a 5% CO2 incubator for 0, 24, and 48 h at 37℃. Next, rapamycin (GlpBio, USA) and/or FBS (Gibco, USA) were added. Then, the CCK-8 (Solarbio, Beijing, China) solution was added, and the samples were incubated for 0.5 h. The absorbance of each well was measured at 450 nm.
Transwell chamber assay
The cell invasion ability was evaluated using Transwell chambers precoated with Matrigel. A total of 1 × 104 cells were inoculated into the top chamber and incubated at 37 °C with 5% CO2 for 48 h. After cells on top of the filter were removed, cells on the bottom of the filter were fixed in 4% paraformaldehyde and stained with 1% crystal violet (Beyotime, P.R. China). The invading cells were counted under a microscope.
Approximately 5 × 105 HeLa cells were seeded in 6-well plates. The cells were washed with PBS three times, and fresh medium was added. Cells were cultured in a 37 ℃ 5% CO2 incubator, photographed, and recorded under the microscope at 0, 24, and 48 h.
The total RNA extracted was detected by EpiQuik ™ m6A RNA Methylation Quantification Kit (Epigentek, USA), following the manufacturer’s instructions.
All data processing and analysis were completed in R (version 3.6.2). To compare two groups of continuous variables, the statistical significance of the normally distributed variables was estimated using the independent student’s t-test, and the differences between non-normally distributed variables were analyzed using the Mann–Whitney U test (i.e., Wilcoxon rank-sum test). The Chi-square test or Fisher’s exact test was used to compare and analyze the statistical significance between two groups of categorical variables. The correlation coefficients of different genes were calculated using Pearson correlation analysis. R’s survival package was used for survival analysis, the Kaplan–Meier survival curve was used to highlight the survival difference, and the log-rank test was used to evaluate the significance of the survival time difference between two groups . R’s pROC package was used to plot the receiver operating characteristic (ROC) curve and calculate the area under the curve (AUC) to assess the accuracy of the risk score in estimating prognosis. Univariate and multivariate Cox analyses were used to determine independent prognostic factors . All statistical P-values were bilateral, and P < 0.05 was considered statistically significant.
A flow chart of the study is shown in Fig. 1A.
Expression and mutation profile of m6A-related genes in patients with CESC
To analyze the overall expression of m6A-related genes in patients with CESC, we analyzed the genomic mutations and mRNA expression, including single nucleotide polymorphism, CNV, and gene expression level. SNP analysis showed that 37 of 289 samples had single nucleotide mutations in m6A-related genes (Fig. 1B). Moreover, studies on the change frequency of CNV showed that the changes of m6A-related genes in CNV levels were widespread in patients with CESC, and most were focused on the loss of copy number. Fig. 1C and D show the position and frequency of CNV changes on chromosomes, respectively.
Next, we comprehensively analyzed the expression of m6A-related genes in cervical cancer and cervical tissue in the TCGA database, GTEx database, and GSE52903 dataset, and batch effects were removed. PCA results showed significant differences in overall expression levels between tumor samples and normal tissues after batch removal (Fig. 1E). We focused on m6A-related genes and compared the differences of m6A-related genes between the two groups by the Wilcoxon rank-sum test. Differential analysis results showed that many m6A-related genes were differentially expressed between cervical cancer tissues and normal tissues, including METTL14, METTL16, ZC3H13, YTHDC1, and YTHDC2 (Fig. 1F). After further correlation analysis, a heat map was constructed, which showed the correlation of m6A-related gene expression. Most genes showed a positive correlation in CESC tissues (Fig. 1G).
We further selected and analyzed the influence of m6A-related genes on the prognosis of patients with CESC based on the TCGA and GEO databases. The gene network depicted a comprehensive picture of m6A-related gene interactions, source grouping, and their impact on OS in patients with CESC (Fig. 2A).
Subsequently, we constructed a risk scoring system based on the expression of m6A-related genes to quantitatively evaluate the impact of m6A-related genes on the prognosis of each patient with CESC. First, 20 m6A-related genes were included in the LASSO Cox analysis, and 12 genes with the best prognostic value were selected and identified (Fig. 2B, C). Simultaneously, based on the penalty coefficient of important characteristic genes calculated by LASSO Cox analysis, the gene expression was multiplied and added to the corresponding coefficient to establish the risk score, and the final risk score of each sample was calculated. The significance of the LASSO Cox analysis here was to eliminate the multicollinearity between genes (resulting from the use of LASSO), and ultimately, to ensure that the genes included in the multivariate Cox regression are independent. The distribution of risk score, survival status, and expression patterns of characteristic genes were shown in Fig. 2D. Kaplan–Meier analysis showed that the OS of patients with high-risk scores was relatively poor (log-rank P < 0.001; Fig. 2E). The time-dependent ROC analysis of the risk score showed that the risk score had a good predictive ability for the OS of patients with CESC, and the AUC of 1-, 3-, and 5-year OS were 0.734 (0.642–0.825), 0.702 (0.648–0.797), and 0.723 (0.665–0.829), respectively (Fig. 2F).
Influence of m6A risk score on the biological function of patients with CESC
We then divided patients into high- and low-risk groups based on the median m6A risk score for patients with CESC and assessed the changes in biological function between the two groups. To analyze the variation in different pathways in patients with CESC, the enrichment scores of patients with CESC were analyzed using the GSVA method, and the related signaling pathways with significant differential expression in the high- and low-risk patients were shown through heat maps (Fig. 2G). In addition, the analysis revealed significant differences in the enrichment of related biological pathways, such as CD8+ T cell effectors, immune checkpoints, EMT pathways, and angiogenesis, between high- and low-risk patients (P < 0.05; Fig. 2H).
Construction of genetic characteristics of patients with CESC based on the m6A risk model
To determine the underlying biological characteristics of different m6A-associated phenotypes, the LIMMA package was used to analyze differential genes between different risk models in the dataset and yielded 207 DEGs. The DEGs are shown in Additional file 2 materials. Subsequently, based on the DEGs, patients with CESC were divided into three different subtypes named Geneclusters A, B, and C (Fig. 3A) using unsupervised clustering. According to the expression and correlation of DEGs between different groups, genes were divided into m6A signature gene-A and B groups. At the same time, survival analysis results showed significant differences in the prognosis of patients with three different gene subtypes. Patients in Genecluster A had the worst prognosis (log-rank P < 0.001; Fig. 3B).
DEGs were divided into two groups, m6A signature-A and m6A signature-B, according to the correlation of gene expression. m6A signature-A gene set had 101 DEGs, and m6A signature-B gene set included 106 DEGs. To explore differences in biological functions between the two groups of genes, functional enrichment analysis was performed (Tables 1 and 2). GO and KEGG enrichment analysis showed that signature-A and -B gene sets showed different unique BPs. Signature-A gene set involved leukotriene catabolic processes, leukotriene B4 catabolic processes, leukotriene B4 metabolic processes, glycosphingolipid biosynthesis-lacto and neolacto series, nicotine addiction, and other pathways (Fig. 3C, D). In contrast, the gene set of overexpressed signature gene-B mainly showed heat generation, positive regulation of angiogenesis, positive regulation of vasculature development, AGE-RAGE signaling pathway in diabetic complications, pertussis, cellular senescence, and other pathways (Fig. 3E, F). The pathways are shown in Additional file 3 materials.
Construction of prognostic m6A feature model based on m6A signature gene
To better predict the impact of m6A characteristics on patient prognosis, we constructed a new prognostic-related risk scoring system. According to the expression of m6A Signature-A and -B gene sets in patients with CESC, the corresponding PCA1 of m6A signature-A and -B gene sets in each patient was calculated by PCA, and the corresponding m6A score of each patient was obtained by subtraction, which was named the m6A group. Similarly, patients were divided into high- and low-risk groups based on the median prognostic model scores. The Sankey diagram showed the corresponding gene group (Geneclusters A, B, and C) of each patient with CESC, the prognosis model m6A group, and the corresponding relationship among the survival statuses of patients (Fig. 3G). Simultaneously, survival analysis showed that the prognostic score model could accurately predict the OS of patients with CESC (log-rank P < 0.001; Fig. 3H).
Subsequently, we analyzed the effect of the high- and low-risk m6A group on patients’ biology-related functions. GSEA showed that epithelial-mesenchymal transition, angiogenesis, Myc target V1, and other pathways were mainly enriched in the high-risk group (Fig. 3I), whereas bile acid metabolism, KRAS signaling DN, estrogen response late, and other pathways were significantly enriched in low-risk patients (Fig. 3J) (Table 3).
Influence of m6A risk score on immune characteristics of patients with CESC
Next, we assessed the impact of the m6A risk score on the overall immune characteristics and different levels of immune cell infiltration in patients with CESC. The immune score was significantly lower in the high-risk group than in the low-risk group (P = 0.025, Fig. 4A), whereas the stromal score was not significantly different (P = 0.64, Fig. 4B). Meanwhile, the CIBERSORT algorithm was further used to evaluate the invasion levels of 22 different immune cells (Fig. 4C). The difference analysis showed that the infiltration levels of various immune cell subsets differed significantly between the high- and low-risk groups (Fig. 4D), including T cells CD8, T cells, CD4 memory activated T cells, follicular helper cells, and M1 macrophages.
Influence of m6A risk score on genomic changes in patients with CESC
Subsequently, we further assessed the impact of the m6A risk score on changes in levels of genetic variation in patients with CESC, including single nucleotide polymorphism and CNV. Analysis of single nucleotide mutations in driver genes during common tumorigenesis revealed that SNP levels in different driver genes were different between the high- and low-risk subgroups (Fig. 5A). Simultaneously, the overall level analysis showed that the TMB (P = 0.035, Fig. 5B) and tumor neoantigen levels (P = 0.0079, Fig. 5C) significantly differed between the high- and low-risk groups. Moreover, studies on the change frequency of CNV showed that CNV changes in high-risk patients were mainly focused on gene amplification (Fig. 5D), whereas in low-risk patients, CNV changes were mainly manifested as gene copy number deletion (Fig. 5E).
Sensitivity analysis of patients with CESC to different small molecule drugs based on m6A risk score
To analyze the sensitivity of the m6A risk score to different drugs and small molecule substances in patients, we downloaded the cell line gene mutation data and IC50 values of different anticancer drugs from the GDSC database. IC50 values for patients with CESC were predicted based on the reactivity of cell lines to 138 different chemotherapeutic agents and small molecule anticancer agents. The results showed that the IC50 values of multiple chemotherapeutic agents and small molecule anticancer agents significantly differed between patients with high and low m6A risk scores (P < 0.001; Fig. 6), especially, KIN001.135, Akt.inhibitor.viii, and rapamycin.
Construction of the clinical prediction model based on m6A risk score
Subsequently, we assessed the impact of the m6A risk score on the outcome of patients with CESC. Univariate and multivariate Cox analysis showed that the m6A risk score was an independent risk factor for predicting the prognosis of patients with CESC (Fig. 7A). Subsequently, the m6A risk score was combined with different clinicopathological features to construct a predictive rosette to predict the OS of patients with CESC (Fig. 7B). The C-index was used to calculate the distinguishing ability of the lipogram, which showed a high degree of differentiation [0.736 (0.680–0.792)]. Moreover, the calibration curves showed a good agreement between the 1-, 3-, and 5-year OS estimates from the rosettes and the actual patient observations (Fig. 7C).
ZC3H13 expression in tissues
We used qPCR to study ZC3H13 expression in ten normal cervical tissues and ten cervical cancer tissues and found that ZC3H13 was significantly downregulated in cervical cancer tissues (Fig. 8A). The IHC analysis showed that ZC3H13 was localized in the nucleus and had a low expression in cervical cancer tissues (Fig. 8B). We performed gray level difference analysis of IHC results using the Image-Pro Plus software. The results showed that the mean integral optical density (IOD) of ZC3H13 in normal cervical tissues was significantly higher than that of ZC3H13 in cervical cancer tissues. The WB analysis also demonstrated the low expression of ZC3H13 in cervical cancer tissues (Fig. 8C).
Function of ZC3H13 in cervical cancer and effects of rapamycin on cell phenotypes
The CCK8 assay showed that the proliferation of HeLa and SiHa cells was significantly enhanced after the ZC3H13 knockdown (Fig. 8D). Further, the CCK8 assay showed a significant reduction in the proliferation of HeLa and SiHa cells after rapamycin treatment (Fig. 8E). The wound-healing assay showed that ZC3H13 knockdown significantly enhanced the invasion ability of HeLa and SiHa cells (Fig. 8F). Transwell experiments showed that ZC3H13 knockdown significantly enhanced the cell migration ability (Fig. 8G). The wound-healing assay showed a significant reduction in the invasion abilities of HeLa and SiHa cells after rapamycin treatment (Fig. 8F). The Transwell assay showed a significant reduction in the migration abilities of HeLa and SiHa cells after rapamycin treatment (Fig. 8G). The m6A ELISA showed that the m6A levels decreased significantly after the ZC3H13 knockdown (Fig. 8H) and that the m6A levels increased significantly in both HeLa and SiHa cells after rapamycin treatment (Fig. 8I).
As one of the most common RNA modifications, m6A mRNA methylation is closely associated with cervical cancer. m6A mRNA methylation might promote cervical cancer development. The m6A level was significantly reduced in cervical cancer compared to adjacent normal tissue. The reduction in m6A levels significantly correlated with the FIGO stage, tumor size, differentiation, lymph invasion, and cancer recurrence . m6A methyltransferase methyltransferase-like 3 (METTL3) enhanced the stability of FOXD2-AS1, and its expression was maintained. METTL3/FOXD2-AS1 accelerated cervical cancer progression via an m6A-dependent modality . METTL3 enhanced the HK2 stability through YTHDF1-mediated m6A modification, thereby promoting the Warburg effect in cervical cancer . GAS5-AS1 interacted with the tumor suppressor GAS5 and increased its stability by interacting with RNA demethylase ALKBH5 and decreasing GAS5 m6A modification. m6A-mediated GAS5 RNA degradation relied on the m6A reader protein YTHDF2-dependent pathway . FTO can control the m6A modification of E2F1 and Myc transcripts to regulate the proliferation and migration of cervical cancer cells . The FTO-mediated stabilization of HOXC13-AS epigenetically upregulated FZD6 and activated Wnt/β-catenin signaling, driving CC proliferation, invasion, and EMT . FTO and its substrate m6A may be critical factors for regulating chemo-radiotherapy resistance . YTHDF1 regulated the translation of RANBP2, which potentiated the growth, migration, and invasion of cervical cancer cells in an m6A-dependent manner without any effect on its mRNA expression . circARHGAP12 exerted the oncogenic role in cervical cancer progression through the m6A-dependent IGF2BP2/FOXM1 pathway . KCNMB2-AS1 and IGF2BP3 formed a positive regulatory circuit that increased the tumorigenic effect of KCNMB2-AS1 in cervical cancer . m6A-associated downregulation of miR-193b promoted cervical cancer aggressiveness by targeting CCND1 . ZFAS1 and its m6A modification may be a promising target for cervical cancer treatment .
According to our study, differential analysis results showed differences in the expression of multiple m6A-related genes between cervical cancer tissues and normal tissues; furthermore, these genes were related to prognosis. The established m6A risk scoring system is closely related to the biological function of cervical cancer, immune invasion of cervical cancer, and sensitivity to small molecule drugs. The model based on the m6A signature genes and the model based on the m6A risk score can accurately predict the prognosis of patients with cervical cancer. Studies showed that m6A was closely related to epithelial–mesenchymal transition, angiogenesis, Myc targets V1, KRAS signal pathway, and estrogen in the occurrence and development of tumors [43, 50,51,52,53,54]. Our results showed that epithelial-mesenchymal transition, angiogenesis, Myc targets V1, and other pathways are mainly enriched in the high-risk groups, while bile acid metabolism, KRAS signaling DN, late estrogen response, and other pathways are significantly enriched in low-risk patients. The results of the present study may help understand the relationship between these pathways and m6A, revealing their mechanism in cervical cancer occurrence and development. We also showed that the m6A regulator correlates with the survival and clinicopathological characteristics of patients with CESC. The m6A regulator-based prognostic signature may predict the prognosis of CESC .
According to our analysis, the IC50 of multiple chemotherapy drugs differed significantly between patients with high and low m6A risk scores. Meanwhile, previous studies have also shown a close relationship between m6A, mTOR pathway, and rapamycin. mTORC1 could stimulate oncogenic signaling and control anabolic cell growth via m6A [56, 57]. m6A has been reported to play a fundamental role in the function of the mTOR pathway in gastrointestinal cancer . Reductions in the m6A levels increased the expression of mTORC2 in endometrial cancer . The FTO protein was reported to participate in tumorigenesis via interactions with the target of the mammalian target protein Rapamycin (mTOR) . METTL3 was shown to cause the activation of mTORC1 signaling and colorectal cancer development in an m6A-dependent manner . Rapamycin may regulate the m6A levels and affect the biological characteristics of OSCC cells . METTL3-high tumors showed more sensitivity to everolimus, a rapamycin analog, in gastric cancer .
Our study showed that knocking down ZC3H13, an m6A methyltransferase, can promote the proliferation, migration, and invasion abilities of cervical cancer cells and reduce the m6A levels. Rapamycin suppresses the proliferation, migration, and invasion abilities of cancer cells and can enhance the m6A levels, demonstrating the role of the m6A methylation modulator in the development of cervical cancer and the effectiveness of the drugs targeting the modulators of m6A methylation in treating cervical cancer.
In conclusion, our study shows that m6A regulatory factors are closely related to the occurrence, development, immune invasion, drug sensitivity, and prognosis of cervical cancer. Drugs targeting the factors regulating m6A offer good prospects for treating cervical cancer.
Availability of data and materials
The datasets used or analyzed during the current study are available from the TCGA GDC website (https://portal.gdc.cancer.gov/).
Bray F, Ferlay J, Soerjomataram I, et al. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2018;68:394–424. https://doi.org/10.3322/caac.21492.
Sung H, Ferlay J, Siegel RL, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2021;71:209–49. https://doi.org/10.3322/caac.21660.
Gultekin M, Ramirez PT, Broutet N, Hutubessy R. World Health Organization call for action to eliminate cervical cancer globally. Int J Gynecol Cancer. 2020;30:426–7. https://doi.org/10.1136/ijgc-2020-001285.
Delaunay S, Frye M. RNA modifications regulating cell fate in cancer. Nat Cell Biol. 2019;21:552–9. https://doi.org/10.1038/s41556-019-0319-0.
Helm M, Motorin Y. Detecting RNA modifications in the epitranscriptome: predict and validate. Nat Rev Genet. 2017;18:275–91. https://doi.org/10.1038/nrg.2016.169.
Harcourt EM, Kietrys AM, Kool ET. Chemical and structural effects of base modifications in messenger RNA. Nature. 2017;541:339–46. https://doi.org/10.1038/nature21351.
Roundtree IA, Evans ME, Pan T, He C. Dynamic RNA modifications in gene expression regulation. Cell. 2017;169(7):1187–200. https://doi.org/10.1016/j.cell.2017.05.045.
He L, Li H, Wu A, et al. Functions of N6-methyladenosine and its role in cancer. Mol Cancer. 2019;18:176. https://doi.org/10.1186/s12943-019-1109-9.
Boccaletto P, Machnicka MA, Purta E, et al. MODOMICS: a database of RNA modification pathways. Nucleic Acids Res. 2018;46:D303-d7. https://doi.org/10.1093/nar/gkx1030.
Sun T, Wu R, Ming L. The role of m6A RNA methylation in cancer. Biomedicine & pharmacotherapy = Biomedecine & pharmacotherapie. 2019;112:108613. https://doi.org/10.1016/j.biopha.2019.108613.
Zhou Y, Kong Y, Fan W, et al. Principles of RNA methylation and their implications for biology and medicine. Biomedicine & pharmacotherapy = Biomedecine & pharmacotherapie. 2020;131:110731. https://doi.org/10.1016/j.biopha.2020.110731.
Zhou J, Wang J, Hong B, et al. Gene signatures and prognostic values of m6A regulators in clear cell renal cell carcinoma—a retrospective study using TCGA database. Aging. 2019;11:1633–47. https://doi.org/10.18632/aging.101856.
Tu Z, Wu L, Wang P, et al. N6-methylandenosine-related lncRNAs are potential biomarkers for predicting the overall survival of lower-grade glioma patients. Front Cell Dev Biol. 2020;8:642. https://doi.org/10.3389/fcell.2020.00642.
Yang G, Sun Z, Zhang N. Reshaping the role of m6A modification in cancer transcriptome: a review. Cancer Cell Int. 2020;20:353. https://doi.org/10.1186/s12935-020-01445-y.
Chandola U, Das R, Panda B. Role of the N6-methyladenosine RNA mark in gene regulation and its implications on development and disease. Brief Funct Genomics. 2015;14:169–79. https://doi.org/10.1093/bfgp/elu039.
Yang J, Chen J, Fei X, et al. N6-methyladenine RNA modification and cancer. Oncol Lett. 2020;20:1504–12. https://doi.org/10.3892/ol.2020.11739.
Alblas M, Velt KB, Pashayan N, et al. Prediction models for endometrial cancer for the general population or symptomatic women: a systematic review. Crit Rev Oncol Hematol. 2018;126:92–9. https://doi.org/10.1016/j.critrevonc.2018.03.023.
Nebbioso A, Tambaro FP, Dell’Aversana C, Altucci L. Cancer epigenetics: moving forward. PLoS Genet. 2018;14: e1007362. https://doi.org/10.1371/journal.pgen.1007362.
Xu Z, Peng B, Cai Y, et al. N6-methyladenosine RNA modification in cancer therapeutic resistance: current status and perspectives. Biochem Pharmacol. 2020;182: 114258. https://doi.org/10.1016/j.bcp.2020.114258.
Ma J, Yang D, Ma XX. Immune infiltration-related N6-methyladenosine RNA methylation regulators influence the malignancy and prognosis of endometrial cancer. Aging. 2021;13:16287–315. https://doi.org/10.18632/aging.203157.
Zhang H, Meltzer P, Davis S. RCircos: an R package for Circos 2D track plots. BMC Bioinformatics. 2013;14:244. https://doi.org/10.1186/1471-2105-14-244.
Mayakonda A, Lin DC, Assenov Y, et al. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 2018;28:1747–56. https://doi.org/10.1101/gr.239244.118.
Medina-Martinez I, Barrón V, Roman-Bassaure E, et al. Impact of gene dosage on gene expression, biological processes and survival in cervical cancer: a genome-wide follow-up study. PLoS ONE. 2014;9: e97842. https://doi.org/10.1371/journal.pone.0097842.
Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43: e47. https://doi.org/10.1093/nar/gkv007.
Chen S, Yang D, Lei C, et al. Identification of crucial genes in abdominal aortic aneurysm by WGCNA. PeerJ. 2019;7: e7873. https://doi.org/10.7717/peerj.7873.
Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16:284–7. https://doi.org/10.1089/omi.2011.0118.
Subramanian A, Tamayo P, Mootha VK, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102:15545–50. https://doi.org/10.1073/pnas.0506580102.
Mariathasan S, Turley SJ, Nickles D, et al. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature. 2018;554:544–8. https://doi.org/10.1038/nature25501.
Rosenberg JE, Hoffman-Censits J, Powles T, et al. Atezolizumab in patients with locally advanced and metastatic urothelial carcinoma who have progressed following treatment with platinum-based chemotherapy: a single-arm, multicentre, phase 2 trial. Lancet. 2016;387:1909–20. https://doi.org/10.1016/S0140-6736(16)00561-4.
Şenbabaoğlu Y, Gejman RS, Winer AG, et al. Tumor immune microenvironment characterization in clear cell renal cell carcinoma identifies prognostic and immunotherapeutically relevant messenger RNA signatures. Genome Biol. 2016;17:231. https://doi.org/10.1186/s13059-016-1092-z.
Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010;26:1572–3. https://doi.org/10.1093/bioinformatics/btq170.
Newman AM, Liu CL, Green MR, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12:453–7. https://doi.org/10.1038/nmeth.3337.
Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612. https://doi.org/10.1038/ncomms3612.
Yang W, Soares J, Greninger P, et al. Genomics of Drug Sensitivity in Cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res. 2013;41:D955-61. https://doi.org/10.1093/nar/gks1111.
Geeleher P, Cox N, Huang RS. pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS ONE. 2014;9: e107468. https://doi.org/10.1371/journal.pone.0107468.
Zeng L, Fan X, Wang X, et al. bioinformatics analysis based on multiple databases identifies hub genes associated with hepatocellular carcinoma. Curr Genomics. 2019;20:349–61. https://doi.org/10.2174/1389202920666191011092410.
Robin X, Turck N, Hainard A, et al. pROC: an open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinformatics. 2011;12:77. https://doi.org/10.1186/1471-2105-12-77.
Wang X, Li Z, Kong B, et al. Reduced m(6)A mRNA methylation is correlated with the progression of human cervical cancer. Oncotarget. 2017;8:98918–30. https://doi.org/10.18632/oncotarget.22041.
Ji F, Lu Y, Chen S, et al. m(6)A methyltransferase METTL3-mediated lncRNA FOXD2-AS1 promotes the tumorigenesis of cervical cancer. Molecular therapy oncolytics. 2021;22:574–81. https://doi.org/10.1016/j.omto.2021.07.004.
Wang Q, Guo X, Li L, et al. N(6)-methyladenosine METTL3 promotes cervical cancer tumorigenesis and Warburg effect through YTHDF1/HK2 modification. Cell Death Dis. 2020;11:911. https://doi.org/10.1038/s41419-020-03071-y.
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 Trans Res. 2020;12:8137–46.
Zou D, Dong L, Li C, et al. The m(6)A eraser FTO facilitates proliferation and migration of human cervical cancer cells. Cancer Cell Int. 2019;19:321. https://doi.org/10.1186/s12935-019-1045-1.
Wang T, Li W, Ye B, et al. FTO-stabilized lncRNA HOXC13-AS epigenetically upregulated FZD6 and activated Wnt/β-catenin signaling to drive cervical cancer proliferation, invasion, and EMT. J BUON. 2021;26:1279–91.
Zhou S, Bai ZL, Xia D, et al. FTO regulates the chemo-radiotherapy resistance of cervical squamous cell carcinoma (CSCC) by targeting β-catenin through mRNA demethylation. Mol Carcinog. 2018;57:590–7. https://doi.org/10.1002/mc.22782.
Wang H, Luo Q, Kang J, et al. YTHDF1 aggravates the progression of cervical cancer through m(6)A-Mediated Up-Regulation of RANBP2. Front Oncol. 2021;11: 650383. https://doi.org/10.3389/fonc.2021.650383.
Ji F, Lu Y, Chen S, et al. IGF2BP2-modified circular RNA circARHGAP12 promotes cervical cancer progression by interacting m(6)A/FOXM1 manner. Cell Death Discovery. 2021;7:215. https://doi.org/10.1038/s41420-021-00595-w.
Zhang Y, Wang D, Wu D, et al. Long Noncoding RNA KCNMB2-AS1 Stabilized by N(6)-methyladenosine modification promotes cervical cancer growth through acting as a competing endogenous RNA. Cell Transplant. 2020;29:963689720964382. https://doi.org/10.1177/0963689720964382.
Huang C, Liang J, Lin S, et al. N(6)-Methyladenosine associated silencing of miR-193b promotes cervical cancer aggressiveness by targeting CCND1. Front Oncol. 2021;11: 666597. https://doi.org/10.3389/fonc.2021.666597.
Yang Z, Ma J, Han S, et al. ZFAS1 exerts an oncogenic role via suppressing miR-647 in an m(6)A-dependent manner in cervical cancer. Onco Targets Ther. 2020;13:11795–806. https://doi.org/10.2147/OTT.S274492.
Li J, Xie G, Tian Y, et al. RNA m 6 A methylation regulates dissemination of cancer cells by modulating expression and membrane localization of β-catenin. Mol Ther. 2022. https://doi.org/10.1016/j.ymthe.2022.01.019.
Yang Z, Wang T, Dejun Wu, et al. RNA N6-methyladenosine reader IGF2BP3 regulates cell cycle and angiogenesis in colon cancer. J Exp Clin Cancer Res. 2020;39:203. https://doi.org/10.1186/s13046-020-01714-8.
Wang Y, Jia-Huan Lu, Qi-Nian Wu, et al. LncRNA LINRIS stabilizes IGF2BP2 and promotes the aerobic glycolysis in colorectal cancer. Mol Cancer. 2019;18:174. https://doi.org/10.1186/s12943-019-1105-0.
Lou X, Ning J, Liu W, et al. YTHDF1 Promotes Cyclin B1 Translation through m 6 A modulation and contributes to the poor prognosis of lung adenocarcinoma with KRAS/TP53 co-mutation. Cells. 2021;10:1669. https://doi.org/10.3390/cells10071669.
Gao Y, Ouyang X, Zuo Li, et al. R-2HG downregulates ERα to inhibit cholangiocarcinoma via the FTO/m6A-methylated ERα/miR16-5p/YAP1 signal pathway. Mol Ther Oncolytics. 2021;23:65–81. https://doi.org/10.1016/j.omto.2021.06.017.
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:1444.
Cho S, Lee G, Pickering BF, Jang C, Park JH, He L, Mathur L, Kim SS, Jung S, Tang HW, et al. mTORC1 promotes cell growth via m(6)A-dependent mRNA degradation. Mol Cell. 2021;81:2064-2075.e2068.
Villa E, Sahu U, O’Hara BP, Ali ES, Helmin KA, Asara JM, Gao P, Singer BD, Ben-Sahra I. mTORC1 stimulates cell growth through SAM synthesis and m(6)A mRNA-dependent control of protein synthesis. Mol Cell. 2021;81:2076-2093.e2079.
Zhao Q, Zhao Y, Hu W, Zhang Y, Wu X, Lu J, Li M, Li W, Wu W, Wang J, et al. m(6)A RNA modification modulates PI3K/Akt/mTOR signal pathway in gastrointestinal cancer. Theranostics. 2020;10:9528–43.
Liu J, Eckert MA, Harada BT, Liu SM, Lu Z, Yu K, Tienda SM, Chryplewicz A, Zhu AC, Yang Y, et al. m(6)A mRNA methylation regulates AKT activity to promote the proliferation and tumorigenicity of endometrial cancer. Nat Cell Biol. 2018;20:1074–83.
Lan N, Lu Y, Zhang Y, Pu S, Xi H, Nie X, Liu J, Yuan W. FTO—a common genetic basis for obesity and cancer. Front Genet. 2020;11: 559138.
Chen H, Gao S, Liu W, Wong CC, Wu J, Wu J, Liu D, Gou H, Kang W, Zhai J, et al. RNA N(6)-Methyladenosine methyltransferase METTL3 facilitates colorectal cancer by activating the m(6)A-GLUT1-mTORC1 axis and is a therapeutic target. Gastroenterology. 2021;160:1284-1300.e1216.
Wang F, Liao Y, Zhang M, Zhu Y, Wang W, Cai H, Liang J, Song F, Hou C, Huang S, et al. N6-methyladenosine demethyltransferase FTO-mediated autophagy in malignant development of oral squamous cell carcinoma. Oncogene. 2021;40:3885–98.
Sun Y, Li S, Yu W, Zhao Z, Gao J, Chen C, Wei M, Liu T, Li L, Liu L. N(6)-methyladenosine-dependent pri-miR-17-92 maturation suppresses PTEN/TMEM127 and promotes sensitivity to everolimus in gastric cancer. Cell Death Dis. 2020;11:836.
The study was supported by the Clinical Medical Research Center of Gynecology and Obstetrical Diseases (gynecological tumors) of Henan Province.
This work was supported by the Department of science and technology of Henan Province (grant nos. 212102310466 and 182102410095), the Department of Education of Henan Province (grant no. 2019GGJS004), and the Henan Provincial Health Commission (grant no. SBGJ202002119).s
Ethics approval and consent to participate
This study was performed in accordance with the Helsinki Declaration and approved by the ethics committee of the Second Affiliated Hospital of Zhengzhou University.
Consent for publication
The authors declare no conflict of interest.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Lu, X., Li, R., Ying, Y. et al. Gene signatures, immune infiltration, and drug sensitivity based on a comprehensive analysis of m6a RNA methylation regulators in cervical cancer. J Transl Med 20, 385 (2022). https://doi.org/10.1186/s12967-022-03600-7
- Immune infiltration
- Drug sensitivity
- Cervical cancer