Characterization of tumor microenvironment and tumor immunology based on the double-stranded RNA-binding protein related genes in cervical cancer
Journal of Translational Medicine volume 21, Article number: 647 (2023)
Cervical cancer is one of the most common gynecological cancers threatening women’s health worldwide. Double-stranded RNA-binding proteins (dsRBPs) regulate innate immunity and are therefore believed to be involved in virus-related malignancies, however, their role in cervical cancer is not well known.
We performed RNA-seq of tumor samples from cervical cancer patients in local cohort and also assessed the RNA-seq and clinical data derived from public datasets. By using single sample Gene Set Enrichment Analysis (ssGSEA) and univariate Cox analysis, patients were stratified into distinct dsRBP clusters. Stepwise Cox and CoxBoost were performed to construct a risk model based on optimal dsRBPs clusters-related differentially expressed genes (DEGs), and GSE44001 and CGCI-HTMCP-CC were employed as two external validation cohorts. Single cell RNA sequencing data from GSE168652 and Scissor algorithm were applied to evaluated the signature-related cell population.
The expression of dsRBP features was found to be associated with HPV infection and carcinogenesis in CESC. However, only Adenosine deaminases acting on RNA (ADAR) and Dicer, Drosha, and Argonautes (DDR) exhibited significant correlations with the overall survival (OS) of CESC patients. Based on these findings, CESC patients were divided into three dsRBP clusters. Cluster 3 showed superior OS but lower levels of ADAR and DDR. Additionally, Cluster 3 demonstrated enhanced innate immunity, with significantly higher activity in cancer immunity cycles, immune scores, and levels of tumor-infiltrating immune cells, particularly CD8+ T cells. Furthermore, a risk model based on nine dsRBP cluster-related DEGs was established. The accuracy of survival prediction for 1 to 5 years was consistently above 0.78, and this model’s robust predictive capacity was confirmed by two external validation sets. The low-risk group exhibited significantly higher levels of immune checkpoints, such as PDCD1 and CTLA4, as well as a higher abundance of CD8+ T cells. Analysis of single-cell sequencing data revealed a significant association between the dsRBP signature and glycolysis. Importantly, low-risk patients showed improved OS and a higher response rate to immunotherapy, along with enduring clinical benefits from concurrent chemoradiotherapy.
dsRBP played a crucial role in the regulation of prognosis and tumor immunology in cervical cancer, and its prognostic signature provides a strategy for risk stratification and immunotherapy evaluation.
Cervical cancer is the fourth leading cause of female cancers globally, resulting in approximately 600,000 new cases and more than 300,000 deaths in 2020 . Numerous studies have established a strong association between cervical cancer and high-risk human papillomavirus (hrHPV) infection, particularly hrHPV types 16 and 18. It has been observed that hrHPV infection can lead to the development of malignant tumors through the progression of precancerous lesions [2, 3]. Notable progress has been achieved in the screening, diagnosis, and treatment of cervical cancer. Primary hrHPV testing and cervical cytology screening are among the significant advancements that have contributed to a significant decrease in the incidence and mortality of cervical cancer among individuals aged 21 to 65 . In the management of cervical cancer, radical hysterectomy and chemoradiation remain the primary treatment approaches. Additionally, it has been observed that the chemotherapy drug bevacizumab provides survival benefits for patients with advanced disease . Nonetheless, these conventional treatments have been associated with adverse effects, and despite advancements, metastatic cervical cancer is still considered incurable . Therefore, there is an urgent need to identify reliable biomarkers for cervical cancer that can facilitate more effective clinical management.
Double-stranded RNA-binding protein (dsRBP) is defined to identify the universal RNA double-stranded body through a highly conserved double-stranded RNA binding domain (dsRBD). In eukaryotes, dsRBP is mainly responsible for the RNA editing, stability maintenance, and functionality of RNA . There are mainly seven subtypes in dsRBP, including RIG-i-like receptor (RLR), protein kinase R (PKR), oligoadenylate synthases (OAS) and RNase L, adenosine deaminases acting on RNA (ADAR), Dicer, Drosha and Argonautes (DDR), PACT and TRBP, helicase, which are all involved in the antiviral innate immunity of mammalian . Studies have shown that dsRBPs are related to the carcinogenesis and progression in various types of cancers. For instance, RIG-I, a member of the RLR family, has a strong association with hepatocellular carcinoma, primarily impacting the biological activities of matrix metalloproteinase-9 . As for PKR, recent research has indicated that its overexpression is linked to poor survival rates among breast cancer patients . Overexpression of ADAR1 has been significantly associated with hematological malignancies . Moreover, MiR-346 has been found to promote up-regulation of Ago2 protein expression, thereby enhancing the proliferation and migration of cervical cancer cells . It has been observed that the induction of dsRNA stress-encoded neoantigens can enhance interferon signaling and tumor immunogenicity . Moreover, the activation of dsRNA stress-type 1 interferon signaling has been found to stimulate anti-tumor T cell immunity and inhibit tumor growth, making it a potential approach to increasing the sensitivity of immune checkpoints inhibitors (ICIs) in poorly immunogenic tumors . Meanwhile, it has been demonstrated that the ablation of KMT2D, a common occurrence in multiple types of cancer, has the potential to activate dsRNA-interferon signaling, leading to an enhanced immunotherapeutic efficacy . Despite the aforementioned findings, to date, there still lacks study that comprehensively investigated the potential prognostic or immunogenic functions of the seven subtypes of dsRBPs in cervical cancer.
Our research was implemented to determine the heterogeneity in the expression feature of dsRBPs and the correlation between dsRBPs and the clinical characteristics, genomic profiles, immune cell infiltration of cervical cancer. Based on the ssGSEA and univariate Cox analysis, we obtained the dsRBPs that were associated with the OS of cervical cancer, and clustered the TCGA-cervical squamous cell carcinoma (CESC) dataset into distinct molecular clusters. Besides, the survival curves, genetic mutation analysis, TME landscape, and immune features among the different molecular patterns were investigated. Eventually, we established a risk model on basis of cluster-associated DEGs, which can further identify the prognostic role and therapeutic response of dsRBPs in cervical cancer.
Materials and methods
Data download and preparation
As a training cohort, 304 patients with cervical cancer were collected from The Cancer Genome Atlas (TCGA) database (http://xena.ucsc.edu/), along with their RNA expression details and clinical information. The HPV infection status of 178 samples from the TCGA-CESC cohort was obtained from a publicly available study . The GSE44001 from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), CGCI-HTMCP-CC from Genomic Data Commons Data Portal (GDC) database (https://portal.gdc.cancer.gov/) and transcriptomic data of cervical patients from local cohorts were set as independent validation cohorts. Additionally, the immunotherapy dataset IMvigor210 was downloaded from IMvigor210CoreBiologies (http://research-pub.gene.com/IMvigor210CoreBiologies), and the chemotherapy dataset GSE168009 was obtained from the GEO database as the source to evaluate the association between identified signature and treatment efficacy. Next, the seven dsRBP subtypes from the previous study  were utilized as the keywords in the GeneCard (https://www.genecards.org/), resulting in a collection of 35 dsRBP genes (Additional file 2: Table S1).
RNA sequencing in the local cervical cancer patients
Tumor and matched normal tissues were collected from 15 cervical cancer patients in the Fudan university shanghai cancer center (FUSCC) cohort to perform RNA-seq. This study was approved by the Ethics Committee of FUSCC and written informed consent was obtained from all the patients. Before RNA extraction, the tissue was evaluated for tumor cell content, and only those with a tumor purity of at least 20% based on histopathological analysis were eligible for RNA extraction and sequencing. Total RNA from each sample was collected using a FastPure® Cell/Tissue Total RNA Isolation Kit V2 (Vazyme, Jiangsu, China), and the RNA concentration and RNA integrity number (RIN) were measured using a Qubit (Thermo Fisher Scientific, MA, United States) and an Agilent 2100 bioanalyzer (Agilent Technologies, CA, United States), respectively. Library construction was performed using the NEBNext® Ultra™ RNA Library Prep Kit for Illumina® Kit (NEB, MA, United States) and sequenced on the Illumina Novaseq-6000 system (Illumina, MA, United States).
Gene set enrichment analysis (ssGSEA) and consensus clustering of dsRBPs
We performed single sample ssGSEA to quantitatively illustrate enrichment scores of the seven dsRBPs gene sets using the Gene Set Variation Analysis (GSVA) R package . In addition, based on the ssGSEA scores, we employed consensus clustering to identify different dsRBP-related patterns in the TCGA-CESC cohort using the k-means algorithms  with the R package “ConsensuClusterPlus” . The cumulative distribution function (CDF) was used to identify the final number of clusters. Furthermore, Kaplan–Meier analysis was applied to explore the overall survival (OS) of distinct clusters with “survival” and “survminer” R packages , and principal components analysis (PCA) was carried out with “FactoMineR” package.
The somatic mutation profile of cervical cancer was identified by mutation annotation format (MAF) file with the R package “maftools” . Fisher’s test was applied to investigate the frequently mutated genes. We obtained copy number variation (CNV) data from TCGA database and employed GISTIC2.0 to identify the amplification and deletion regions in high- and low-risk groups.
XCELL and CIBERSORT were employed to uncover the immune infiltrating abundance of tumor-infiltrating immune cells (TIICs) in the TME. Moreover, we calculated the stromal score, immune score, and ESTIMATE score according to the ESTIMATE algorithm  using the “estimate” R package.
Immunological features evaluation
Accordingly, we investigated the potential association between three clusters and ICIs-related genes, and the expression level of HLA. In addition, we collected multiple immune signatures to compare the differences between different molecular patterns according to Kobayashi (glycolysis, IFN-γ response, inhibitory cells MDSCs, inhibitory cells Tregs, inhibitory molecules, innate immunity, priming activation, proliferation, recognition of tumor cells, and T cells)  and Bagaev (angiogenesis, anti-tumor microenvironment, antigen presentation, B cells, CAF, checkpoint inhibition, cytotoxic T and NK cells, granulocytes, MDSC, Treg, tumor features, tumor promotive immune infiltrate) .
Identification of differentially expressed genes (DEGs)
According to “limma” R package  with criteria of |log2-fold change (FC)| ≥ 1 and p-value < 0.01, we obtain the DEGs between the identified clusters, which were visualized using volcano plots by the “ggplot2” R package . Subsequently, univariate Cox analysis was applied to screen the prognostic cluster-associated DEGs with a p-value < 0.01. Moreover, we screened out optimal DEGs according to whether they were differentially expressed between cancer and normal tissues using the “limma” R package based on prognostic cluster-associated DEGs, and the cutoff criterion was set as p-value < 0.01. Moreover, Kyoto Encyclopedia of Genes and Genomes (KEGG)  and gene set enrichment analysis were performed to identify the biological functions of cluster-associated DEGs, and the hallmark pathways of GSEA analysis were extracted from Molecular Signatures Database (MsigDB, http://www.gsea-msigdb.org/gsea/downloads.jsp).
Construction and verification of risk model
According to the optimal DEGs, we performed 42 combinations of 6 machine learning algorithms, including Ridge, CoxBoost, elastic network (Enet), least absolute shrinkage and selection operator (LASSO) regression analysis, stepwise Cox, and random survival forest (RSF) to filtrate the most valuable risk model. Accordingly, the risk model was constructed depending on the stepwise Cox and CoxBoost with the highest C-index. Stepwise Cox was applied to filtrate the most valuable DEGs, and CoxBoost was used to screen the most reliable model. Then, the patients in the TCGA-CESC database were assigned to the high-risk group and low-risk group based on the median risk score. Kaplan–Meier curves were analyzed, and the receiver operating characteristic (ROC) curves at 12, 24, 36, 48, and 60 months were performed with the “timeROC” R package . GSE44001 and CGCI-HTMCP-CC datasets were extracted as two external validation cohorts to further evaluate the reliability of the risk model. Furthermore, survival analysis was implemented to verify the characteristics of the = high-and low-risk groups stratified by clinical parameters, including age (≤ 60/> 60), T stage (T1–2/T3–4), N stage (N0/N1+), M stage (M0/M1), and pathological stage (stage I–II/stage III–IV).
Single cell RNA sequecing acquisition and SCISSOR analysis
The single-cell RNA sequencing data for the GSE168652 dataset  was obtained from the TISCH website . The cell-type annotations were provided by combining the results from the original study with adjustments made using cell markers and the InferCNV algorithm employed by the TISCH study. We combined single cell data and TCGA bulk expression data using the scissor (2.0.0) method  to investigate cell subpopulations connected to risk groups. First, we classified TCGA-CESC patients into high-risk and low-risk groups, and then we used the scissor approach in conjunction with logistic regression to extract the most relevant cell populations with regard to risk group. We separated these cell groups into scissor-positive and scissor-negative categories, and finally, we identified differentially expressed genes using the “seurat package” (4.3.0) in R. In our current study, the term “scissor+ cells” refers to malignancy cells that exhibit a phenotype consistent with a high-dsRBP risk score in the bulk sequencing data, while “scissor− cells” refer to malignancy cells that exhibit a phenotype similar to the low-dsRBP risk score type.
Immunotherapy response of risk model
Tumor immune dysfunction and exclusion (TIDE) analysis was performed to identify the response of immune checkpoint inhibitors, which has been proved as a excellent immunotherapy predictive biomakrer . T cell dysfunction score, T cell exclusion score, and TIDE score were downloaded from the TIDE website (http://tide.dfci.harvard.edu). Besides, immunophenoscore (IPS) contains four types, such as immunosuppressive cells, immunomodulators, effector cells and MHC molecules, were reported to be a powerful biomarker of immunotherapy response to anti-PD-1 and anti-CTLA-4 therapy . The IPS of cervical cancer patients were collected from the Cancer Immunome Atlas (TCIA, https://tcia.at/home) database. Moreover, the IMvigor210 dataset with anti-PD-L1 therapy was implemented to assess the immunotherapy value of the risk score using the “IMvigor210CoreBiologies” R package .
Based on the Genomics of Drug Sensitivity in Cancer (GDSC) database, we calculated the half-maximum inhibitory concentration (IC50) of cervical cancer patients using R package “pRRophetic” . Furthermore, we applied the GSE168009 dataset with the concurrent chemoradiotherapy (CCRT) to evaluate the chemotherapy response of risk score.
Kaplan–Meier plots and log-rank tests were employed to compare the survival ability of different groups. R software and its related software package (v.4.1.2) were utilized to analyze and process data. Continuous data processing was applied according to Wilcoxon's test. All tests were two-way and the p < 0.05 was performed in all of the analytic approached, suggesting that there was considered statistically significant.
Characteristic the expression patterns of dsRBPs in cervical cancer
Firstly, we analyzed the expression features related to the HPV infection status in CESC patients. In comparison to HPV-negative CESC patients (n = 9), HPV-positive CESC patients in the TCGA-CESC cohort showed significantly lower expression levels of DDR and PACT_TRBP, but higher expression levels of OAS_RNAseL and RLR (Additional file 1: Fig. S1A). This trend remained consistent when we classified HPV-positive patients based on hierarchical HPV calls or clades (Additional file 1: Fig. S1B, C). However, we did not observe any significant differences in any of the dsRBP subtypes between CESC patients with (n = 141) or without (n = 37) HPV DNA integration (Additional file 1: Fig. S1D).
In both TCGA-CESC and FUSCC cohorts (Fig. 1A, B), the expression pattern of dsRBPs-related subtypes or genes is consistent: ADAR and DDR subtypes are significantly downregulated in tumor tissues; whereas helicase, OAS RNAseL, PKR and RLR subtypes are significantly overexpressed. Based on the ssGSEA scores and univariate Cox analysis, we found that only ADAR and DDR subtypes were significantly associated to the OS of patients in the TCGA-CESC cohort (ADAR: HR = 6.42, 95% CI 1.06–38.71, p = 0.0427; DDR: HR = 3.21, 95% CI 1.05–9.84, p = 0.0409, Fig. 1C). Then, we conducted a tumor immunity analysis and found that the DDR and PACT_TRBP subtypes were consistently associated with a non-inflamed tumor feature. Specifically, we observed a significantly negative correlation between these subtypes and immune scores, CD8+ T cell abundance, as well as expression levels of immune checkpoints (Additional file 1: Fig. S2A–C). On the other hand, RLR and OAS-RNaseL subtypes were associated with a better immunogenic feature, characterized by higher immune scores, greater infiltration of M1 macrophages and CD8 T cells, and increased expression of immune checkpoints (Additional file 1: Fig. S2A–C). Furthermore, we observed no significant difference in the expression levels of ADAR or DDR members among CESC patients with different HPV infection statuses, except for AGO1 (Additional file 1: Fig. S3).
According to ADAR and DDR levels, a consensus clustering algorithm was performed to stratify the tumor samples of cervical cancer with the k = 3 as the number of clusters. Thus, CESC patients in TCGA cohort were dispersed in three distinct clusters named cluster 1, 2, 3 (Fig. 1D). Patients of cluster1 or cluster 2 had worse OS than those from cluster 3 (a median OS: 136.2 months vs 69.8 months vs unreached, p = 0.046, Fig. 1E). PCA further corroborated the classification function of the consensus clustering in the TCGA-CESC cohort (Fig. 1F). In addition, the majority of ADAR- and DDR-related genes exhibited up-regulation in the cluster 1/2. Furthermore, cluster 1/2 had a higher proportion of tumors in advanced stages than cluster 3, including significantly more tumors in T3/T4, N1, M1, stage III+ stage IV, which might explain the poor prognosis (Additional file 1: Fig. S4 and Additional file 2: Table S2).
Expression patterns of dsRBPs in different clusters
Different clusters showed the heterogeneity in the expression profile of dsRBPs. As depicted in Fig. 2A, cluster 1/2 had a higher ssGSEA score in most dsRBPs subtypes, except PACT and RLR (Fig. 2A). Both ADAR (p = 2.4e−15) and DDR (p = 2.8e−12) were notably upregualted in cluster 1/2. Specifically, the expression level of each ADAR and DDR subfamily member was considerably greater in cluster 1/2 (Fig. 2B). In the meantime, there was no significant difference in the distribution of HPV-negative or positive CESC samples between cluster 1/2 and 3 (Additional file 1: Fig. S5).
dsRBP cluster 3 reveals an enhanced innate immunity
We next assessed differences in the tumor immunology and TME between dsRBPs clusters. First, cluster 3 was presented with increased activity in comprehensive steps involved in the cancer immunity cycle, especially in transport of immune cells to the tumor (step four) and the infiltration of immune cells into the tumor (step five) (Fig. 3A). In line with this result, cluster 3 showed significantly higher immune scores and ESTIMATE scores (Fig. 3B). Using the Xcell (Fig. 3C) and CIBERSOFT, we found a prominent upregulation of tumor-infiltrating immune cells in cluster 3, with significantly increasing levels in memory B cells, CD8+ T cells, NK cells but lower neutrophils levels. When compareing hub genes involved in cell–cell adhesion, endothelial-to-mesenchymal transition (EMT) and stem cell-like feature, we found cluster 3 was characterized by low cell–cell adhesion, EMT and stem cell-like level (Fig. 3E). Then we analyzed the difference in the expression level of immune checkpoints between cluster 1/2 and 3. With only the exception of LAG3, the majority of immune checkpoints, including CD200, CD276, CD28, CD40LG, CD44, ICOSLG, LAIR1, NELL1, NRP1, and TNFSF14, were overexpressed in the cluster 1/2 (Fig. 3F). On the contrary, there was only slight differences in the HLA family genes between different clusters (Fig. 3G). Furthermore, utilizing the immunogram radar plot, we found that cluster 3 was presented with significantly higher activity in comprehensive pathways involved in tumor immunity, except glycolysis (Fig. 3H, I). To support the difference in the immunity, we also investigated the holistic TME feature by using the method described by Bageav et al. Cluster 3 demonstrated increased activity in antigen presentation, checkpoint inhibition, cytoroxic T and NK cells, and Treg, whereas clusters 1/2 shown significantly higher activity in angiogenesis and tumor characteristics (Fig. 3J, K).
Genomic feature related to dsRBP clusters
Differences in the prevalence of genomic alterations between cluster 1/2 and 3 were analyzed in the TCGA-CESC dataset, and oncoprint showed TTN, PIK3CA, and KMT2C were the most common in both clusters (Fig. 4A). In addition, as displayed in the forest plot, FRYL, DZIP1, FBXL20, CECR2, DCAF8L2, MAGEC3 were significantly more frequently mutated in cluster 3; whereas FLG, BIRC6, and IGSF10 alterations were only identified in cluster 1/2 (all p < 0.05, Fig. 4B). As previous studies have shown the prevalent genomic alterations RTK/AKT/MAPK and TGFβ signaling pathways in cervical cancers, then we investigated the genomic difference in these pathways between dsRBP clusters. Notable more genomic alterations involved in RTK/AKT/MAPK signaling pathway were identified in cluster 1/2, especially in PIK3CA, PTEN, ERBB3 and AKT1 (Fig. 4C). Meanwhile, more EP300 and SMAD4 alterations were found in cluster 1/2, which were hub genes that participated in the regulation of TGFβ signaling pathway (Fig. 4D).
Development of a dsRBP signature
To develop a prognostic dsRBP signature, we analyzed dsRBP cluster-associated DEGs, of which 39 were down-regulated and 2931 were up-regulated (Fig. 5A and Additional file 2: Table S3). Consistent with previous dsRBP clusters-related feature, these DEGs were mainly enriched in cell adhesion and microenvironment [focal adhesion, extracellular matrix (ECM), gap junction], WNT, TGFβ and EMT pathways (Fig. 5B). Then, we obtained the prognostic cluster-associated DEGs based on univariate Cox analysis (Additional file 2: Table S4), which were further filtered by comparing the expression levels of tumor and normal tissues (Additional file 2: Table S5). Additionally, we performed 6 machine learning algorithms, including stepwiseCox, CoxBoost, Enet, LASSO, Ridge, random forest, combined to construct a risk model depending on optimal DEGs, which can identify the most robust and stable risk model with the highest C-index in the TCGA-CESC cohort, GSE44001 cohort and CGCI-HTMCP-CC cohort (Fig. 5C and Additional file 2: Table S6). Finally, a risk model with the best performance was built based on stepwise Cox and CoxBoost, in which stepwise Cox identified 9 most important DEGs (PDE1C, EDA2R, DDN, LEPR, C1GALT1, MUSTN1, ERG, HLF, and FLT1), and the CoxBoost screened out the most robust risk model. The majority of these signature genes (7/9) were associated with a worse survival with the exception of MUSTN1 and HLF (Fig. 5D). The formula of the risk model was as follows:
Then, patients in the TCGA-CESC cohort were dichotomized into high- and low-risk groups based on median risk score (Additional file 2: Table S7). The Kaplan–Meier curve confirmed that the OS of the high-risk group was significantly shorter than that of the low-risk group in the TCGA-CESC dataset (median OS: 46.47 vs unreached, p < 0.0001, Fig. 5E). In the meantime, the 12-, 24-, 36-, 48-, and 60-month prediction accuracy was 0.78, 0.80, 0.79, 0.81, and 0.78, respectively (Fig. 5F). Meanwhile, the stratification survival analysis in the TCGA-CESC cohort demonstrated that the patients in the high-risk group had significantly worse OS for all clinical parameters (Additional file 1: Fig. S6), and the clinical feature information of the high- and low-risk groups was displayed in Additional file 2: Table S8. Univariate and multivariate Cox regression analysis further revealed that dsRBP risk score was the only independent risk factor (Additional file 1: Fig. S7). High-risk group was presented with a higher proportion of patients with advanced tumor stage and neoplasm disease stage (Fig. 5G).
Validation of the dsRBP signature
GSE44001 and CGCI-HTMCP-CC datasets served as two external validation datasets to evaluate the prognosis prediction capacity of the dsRBP signature (Additional file 2: Tables S9, S10). The survival analysis showed that the high-risk group had a worse OS in CGCI-HTMCP-CC cohort (a median OS: 12.50 vs 22.47 months, p = 0.00046; Fig. 6A), and a significantly worse disease-free survival (DFS) in the GSE44001 cohort (a median DFS: unreached vs unreached months, p = 0.016; Fig. 6B). Besides, the AUC for prediction OS at 12-, 18-, 24-, and 27-month were 0.70, 0.66, 0.66, and 0.70, respectively, in CGCI-HTMCP-CC dataset (Figs. 6C). The AUC for predicting DFS at 12-, 24-, 36, 48, and 60 month were 0.58, 0.62, 0.66, 0.66, and 0.64, respectively, in GSE44001 dataset (Fig. 6D). In addition, we evaluated the predictive capacity of the dsRBP signature to that of 10 previously reported risk evaluation models in cervical cancer, and found that our established dsRBP signature outperformed those risk models in predicting survival (Fig. 6E).
Expression feature of the signature-related genes and dsRBP subtypes
Subsequently, we evaluated the expression feature of the nine signature-related genes in cervical tumor tissues and normal tissues. Intriguingly, nearly all signature-related genes were significantly downregulated in the tumor tissue in both the TCGA and local cohort, except DDN and C1GALT1 (Fig. 7A, B). The high-risk group showed distinguishing characteristics of dsRBP subtypes that corresponded to those of dsRBP clusters 1/2, such that high-risk group had significantly higher ADAR, DDR, Helicase, OAS RNasel, and PKR ssGSEA scores, whereas low-risk group had significantly higher RLR ssGSEA scores (Fig. 7C). Most ADAR-related and DDR-related genes exhibited significantly elevated expression level in the high-risk group, including ADAR, ADARB1, AGO1, AGO2, AGO3, AGO4, and DICER1 (Fig. 7D). When comparing HPV-positive and HPV-negative samples, HPV-negative samples had significantly higher dsRBPs risk scores (distribution of HPV-negative or positive CESC samples between cluster 1/2 and 3 (Additional file 1: Fig. S8A). Similarly, when comparing different hierarchical HPV calls, the trend in difference was consistent, but without statistical significance (Additional file 1: Fig. S8B). Additionally, hierarchical HPV clade-negative patients also had higher dsRBPs risk scores compared to those with A7 (n = 120) and A9 (n = 45) (Additional file 1: Fig. S8C). Furthermore, we did not observe any significant differences between CESC patients with (n = 141) or without HPV DNA integration (Additional file 1: Fig. S8D).
Profiling the cell subpopulation associated with dsRBP signature
Among the identified hub genes, we found that only C1GALT1 was significantly overexpressed in the tumor cells, whereas LEP4, ERG, and FLT1 were primarily expressed in the endothelial cells in scRNA-seq dataset (Fig. 8A). To evaluate the association between the bulk RNA sequencing-derived dsRBP risk signature and scRNA-seq profile, we used the recently published computational method Scissor. This analysis identified 1202 Scissor+ cells and 1222 Scissor− cells in the GSE168652 dataset (Fig. 8B). While the majority of Scissor+ and Scissor− cells were malignant cells, the other Scissor+ cells were enriched in endometrial stromal cells, SMCs, and endothelial cells (Fig. 8C), while the other Scissor− cells showed enrichment for SMCs, CD8+ T cells, and fibroblasts (Fig. 8D). Furthermore, a total of 125 genes were differentially expressed between Scissor+ and Scissor− malignant cells, mainly enriched in hypoxia, epithelial-mesenchymal transition, glycolysis, and mTORC1 signaling pathways (Fig. 8E).
TME and tumor immunity related to dsRBP signature
Similar to the feature of dsRBP cluster 1/2, high-risk group had a significantly lower immune score and ESTIMATE score than the low-risk group (Fig. 9A). In the meanwhile, high-risk presented higher expression levels of genes involved in cell–cell adhesion (OCLN, CHD1), EMT (SNAI1, SNAI2, VIM, ZEB1, TWIST1) and stemness (THY1, MME, ITGA6 and ITGB1, Fig. 9B). When we evaluated the difference in TIIC abundance between high- and low-risk groups, we discovered that the high-risk group was associated with a relative “cold” TME. In combination with the xCELL (Fig. 9C) and CIBERSOFT (Fig. 9D) results, the high-risk group had a significantly decreased abundance of CD8+ T cells, memory B cells, and NK cells. Regarding the expression level of immune checkpoints, the high-risk group had significantly higher levels of CD276, CD44, and NRP1, whereas the low-risk group had significantly higher levels of a greater number of immune checkpoints, such as PDCD1, CTLA4, IDO1, IDO2, LAG3, and ICOS (Fig. 9E). In addition, high-risk group had significantly elevated levels in glycolysis and recognition of tumor cells, but lower levels in IFNg response, inhibitor cells (Tregs), innate immunity, priming and activation of immunology and T cells (Fig. 9F, G). In accordance with our previous analyses, Bageav analysis also showed high-risk group had significantly higher levels in angiogenesis, cancer associated fibroblasts (CAFs) and tumor feature but lower levels in anti-tumor microenvironment, antigen presentation, B cells, checkpoint inhibition, cytotoxic T and NK cells and Tregs (Fig. 9H, I).
Somatic mutation and CNVs related ot risk model
We have conducted an evaluation of the mutated prevalence in CESC samples from the TCGA database and have found that the hub genes are seldom altered in this cancer type. Specifically, we have only identified genomic alterations in LEPR, FLT1, PDE1C, and ERG, with a prevalence rate beyond 1% (Additional file 1: Fig. S9). And in 25 cervical cancer samples in local cohort who preformed whole exon sequencing, there was no somatic alteration in those hub genes identified. In TCGA-CESC cohorts, the most prevalent genes were similar between high- and low-risk groups, showing TTN, PIK3CA, and KMT2C as the most prevalent alterations (Fig. 10A, B). Compared to the low-risk group, high-risk group had significantly higher prevalence in ANO7, ARID1A, BRCA2, MYH2, and ERCC5 (Fig. 10C). Among PI3K–MAPK pathway genes, the high-risk group had a higher prevalence of ERBB2 and ERBB3, but a lower frequency of MAPK1, PIK3R1, and AKT1 alterations (Fig. 10D). In contrast, the high-risk group had a higher frequency of genes involved in the TGFβ signaling pathway (Fig. 10E). Subsequently, we analysed the difference in the copy number variants between high- and low-risk groups (Fig. 10F). There were no distinguishing characteristics between the high- and low-risk groups; however, the high-risk group had significantly more Amp at 19p13.2, 3q28, 1p31.3, and 19p13.31 and less Del at 8p23.3 (Fig. 10G).
Estimation of immunotherapy and chemotherapy response associated with dsRBP signature
According to Fig. 11A, the IPS score, IPS–PD1/PDL1/PDL2 blocker score, IPS-CTLA4 blocker score, and IPS-CTLA4 and PD1/PDL1/PDL2 blocker score were all significantly higher in the low-risk group, indicating that CESC patients with lower risk score might respond better to immunotherapy. Meanwhile, the high-risk group had a lower T cell dysfunction score, but a significantly higher T cell exclusion score (Fig. 11B). Based on the IMvigor210 cohort, the patients with high-risk score had a significantly worse OS than those of low-risk score (p = 0.00052, Fig. 11C). In the meantime, low-risk group had significantly higher proportion of objective responders than the high-risk group (Fig. 10D). Then, we analyzed the sensitivity to 26 chemotherapeutic drugs based on the stratification of risk score, and the results revealed that high-risk group was more sensitive to gemcitabine, axitinib, pazopanib and AZD8055; while low-risk group had higher sensitivity to pacitaxel, erlotinib, lapatinib and sunitinib (Fig. 11E). On the other hand, the risk score had good predictive ability in the cervical cancer patients treated with concurrent chemoradiotherapy (CCRT) with an AUC value of 0.95 in the GSE168009 dataset (Fig. 11F). Accordingly, the cervical patients with durable clinical benefits treated with CCRT exhibited significantly lower risk scores (Fig. 11G).
dsRNA typically accumulates as a result of viral infection, which can activate the innate immune system to distinguish it. Despite significant progress in understanding how dsRNA and dsRBP function, their role in cervical cancer remains unclear. While previous studies have extensively demonstrated the positive relationship between dsRBP and viral infection, it is important to note that dsRBP also plays crucial roles in uninfected cells for maintaining normal biological processes. Therefore, it remains unknown whether there are differences in the expression features of different dsRBP subtypes between CESC cases with or without HPV infection. In our study, we observed distinct expression profiles in dsRBP subtypes associated with HPV infection, particularly in DDR, PACT_TRBP, OAS_RNAseL, and RLR. Previous research indicated that HPV-positive cervical carcinoma cell lines exhibit higher levels of DROSHA and DICER mRNA compared to HPV-negative cervical carcinoma cell lines, suggesting dysregulation of DDR-related miRNA levels . However, in our current study, we found a significant downregulation of DDR levels in the HPV-positive CESC samples from the TCGA-CESC cohort. Additionally, no significant difference in the mRNA expression levels of DROSHA and DICER was observed between HPV-positive and HPV-negative samples. The discrepancy in DDR ssGSEA score primarily stemmed from differences in the mRNA expression level of AGO1, which exhibited significantly elevated expression in the HPV-negative CESC samples. Furthermore, Wang and her colleagues also discovered a significant association of OAS3 with CIN3/cancer, including viral infections in CESC. This finding supports the increased expression level of OAS_RNaseL in HPV-positive samples . RNASEL, an enzyme central to interferon-related antiviral and apoptotic responses, has been linked to increased risk for CESC, head and neck squamous cell carcinomas, and breast cancer through single nucleotide polymorphisms such as rs3738579 . RIG-I plays a crucial role as an immune sensor in inducing intrinsic apoptotic cell death and cisplatin-mediated cell killing in CESC. This activation is associated with the activation of natural killer cells in the tumor microenvironment . However, HPV E6 has been found to inhibit the induction of IFN mediated by RIG-I. This inhibition is primarily achieved through targeting the upstream factors TRIM25 and USP15 . Therefore, the observed increase in the expression level of RLR subtype, particularly RIG-I and MDA5, in HPV-positive CESC samples may be attributed to the IFN-stimulated antiviral process. Interestingly, we discovered that HPV-negative CESC patients exhibited significantly higher dsRBPs risk scores compared to HPV-positive patients. This finding aligns with the observation that HPV-negative CESC patients generally have a poorer prognosis .
In our research, we observed a significant downregulation of ADAR in cervical tumor tissues, while Helicase, OAS_RNAsel, PKR, and RLR were upregulated. This finding aligns with the regulation network of dsRNA sensors and modulators in innate immunity, where ADAR mainly functions as a counterpart to these subfamilies [41, 42]. Notably, ADAR and DDR were the only two dsRBP subtypes significantly associated with the clinical outcomes of CESC. ADAR is responsible for the modification of adenosines to inosines in dsRNAs, and there are three ADAR proteins (ADAR1-3) in humans . Previous studies have also supported the role of ADARs as novel oncogenes in CESC, associated with poor prognosis, unfavored pathological factors, and angiogenesis [44, 45]. DDR comprises three conserved RNases involved in RNA interference regulation, playing a role in miRNA and siRNA biogenesis and related mRNA silencing [46, 47]. Dicer, a member of the DDR family, has been identified as a risk factor in CESC, significantly associated with distant metastasis and disease recurrence . Additionally, Drosha exhibits upregulated expression due to copy number gain at chromosome 5, promoting tumor progression in CESC . Furthermore, AGO3 has been found to enhance the proliferation and growth of CESC cells through the Wnt/β-catenin signaling pathway . These studies support our findings regarding the association of ADAR and DDR with cervical cancer prognosis. Interestingly, both ADAR and DDR showed relatively downregulated expression in tumor tissues compared to normal cervical tissues in both TCGA and FUSCC cohorts, consistent with previous study results .
Based on the pattern of ADAR and DDR expression, we classified cervical cancer patients into three clusters, with cluster 3 showing the best survival, confirming the predictive function of ADAR and DDR. Cluster 3 had significantly lower expression level of ADAR and DDR but more TIICs presented in the tumor microenvironment. In the last decades, TIICs are being recognized as play an essential role in TME and are involved in the development of many types of tumors, hence investigation the value of immune cell infiltrates in TME is crucial to provide novel therapeutic methods and improving immunotherapy response rates for cancer . Previous studies have identified the pivotal role of TIICs in the regulaiton of disease development of cervical cancer [53, 54]. Both CIBERSOFT and Xcell revealed that cluster 3 had significantly higher level of CD8+ T cells, which is not only a prognostican factor  but also associated with the response to ICIs in cervical cancer . On the other hand, recent research has suggested that neutrophils, as the key regulator of cancer, can influence the inflammation in cancer and induce the angiogenesis of cancer, thereby accelerating tumor initiation, proliferation, growth, and metastasis [57, 58]. Besides, cancer associated fibroblast is identified as a microenvironmental cell and offers metabolic support for a malignant tumor, which has the biological functions of immunogenic and immunosuppressive , and endothelial cells are proven to be an emerging important factor that participates in the regulation of cancer development . In the present research, the abundance of neutrophils, cancer associated fibroblast, and endothelial cells were elevated in cluster 1/2, which might explain the poor outcomes of patients in cluster 1/2.
Previous research reveals that HPV infection, a crucial step in the initiation and development of cervical cancer will generate an immunosuppressive microenvironment and negate host antitumor immunity; therefore, immunotherapy is a promising developing treatment option for individuals with CESC [61, 62]. However, even anti-PD-1/PD-L1 treatments for CESC have been approved by FDA, the majority of CESC patients have limited antitumor efficacy to the monotherapy of these regimes . Then, combination with other therapies, particularly those targeted at other immune checkpoints to elude immune surveillance in cervical cancer, has garnered significant interest and been the subject of extensive preclinical and clinical investigations . In comparison to cluster 3, we discovered that the more immune checkpoints, such as ADORA2A, CD160, CD200R1, CD276, CD28, CD40LG, CD44, ICOSLG, NELL1, NRP1, TNFSF14, TNFSF18, and TNFSF4, were more abundant in cluster 1/2. Recent work demonstrated that CD28 contributes to the development of cervical cancer and can serve as a prognostic marker for cervical cancer . Furthermore, CD44, as the specific cell adhesion molecule, is proven to increase the migration and invasion of cervical cancer . It may suggest that therapy targeted at these overexpressed immune checkpoints may remold the immunosuppressive microenvironment and then enhance the effect of anti-PD-1/L1 therapy in patients of cluster 1/2.
In the meantime, genomic mutations participate in the growth and survival of tumor cells by giving the tumor a selective advantage and providing valuable insights for early cancer diagnosis, disease monitoring and treatment . In cluster 1/2, the incidence of FLG, BIRC6, and IGSF10 mutations were significantly higher, and it is reported that FLG contributes to the biological activity of barrier function and associates with the poor prognosis of cervical cancer . Furthermore, when focused on the TGF-beta pathway, we found a significant enrichement of alterations in EP300 and SMAD4 in the cluster1/2. SMAD4 is a known tumor suppressor in cervical cancer and an essential regulator of the TGF-beta pathway, which has been widely linked to the progression of disease [69, 70].
To further apply the prognostic potential of dsRBPs, we constructed a risk model for CESC patients. Through stepwise Cox and CoxBoost analyses, we developed a risk model using the expression levels of nine cluster-related DEGs. The predictive ability of the dsRBP-based risk model was robust, as demonstrated in both the training and validation cohorts, which outperformed ten previously published risk models for CESC. Among the genes included in our risk model, PDE1C, EDA2R, LEPR, C1GALT1, MUSTN1, and HLF were identified for the first time to be associated with the prognosis of CESC, warranting further investigation. Previous studies have linked several of these genes to malignant progression in various cancers, including cervical cancer. For instance, DDN has been found to be overexpressed in CESC tissues, and knockdown of the lncRNA DDN-AS1 has been shown to inhibit tumor proliferation and migration . ERG, a regulator of the glycolysis process, has been implicated in enhancing the growth and invasion of CESC through its impact on aerobic glycolysis capacity . Furthermore, FLT1 has been identified as a prognostic biomarker for CESC . C1GALT1 has been associated with the development of gastric carcinogenesis, with overexpression shown to stimulate the PI3K/AKT pathway and regulate integrin α5O glycosylation . HLF has been reported to play a role in regulating the development of triple-negative breast cancer by activating tumor cell macrophage crosstalk, which also affects chemotherapeutic resistance . These findings highlight the potential clinical relevance and significance of our dsRBP-based risk model for cervical cancer patients.
We also investigated the application of the established signature in selecting therapeutic regimens. Studies also have shown that TIDE and or IPS is a highly effective predictor for treatment response to anti-PD-1 and CTLA-4 therapies . The low-risk group displayed lower T cells exclusion score, demonstrating the low-risk group might benefit more from immunotherapy. The evevalted IPS score in the low-risk group further supports the notion that this group exhibits a more favorable response to immunotherapy. These and the findings from IMVigor210 dataset all suggests that our risk model could serve as a potential therapeutic target for estimating the immunotherapy response in patients with CESC. Besides, CCRT is not only the optimal treatment method but also the standard care for locally advanced CESC patients, which can improve the survival rate [76, 77]. Our study investigated the relationship between CCRT and the risk model in the GSE168009 dataset and indicated that the CESC patients with durable clinical benefits showed lower risk scores. Thus, we confirmed that the low-risk patients have a better response to CCRT, which provided a potential direction for the chemotherapy management of CESC patients. In summary, our findings suggest that the dsRBP signature could serve as a valuable tool for identifying potential responders to ICIs and CCRT treatment in CESC. By utilizing this signature, we may be able to improve the clinical outcomes of CESC patients by prioritizing therapies.
To our knowledge, this is the first study to comprehensively investigate the expression profile and prognostic significance of dsRBPs in CESC, utilizing both bulk RNA and single-cell RNA sequencing data. Our study unveiled a distinct and previously uncharacterized expression pattern of dsRBP subtypes in CESC. Furthermore, we identified a unique cluster that exhibited associations with clinical outcomes, genomic characteristics, tumor microenvironment landscape, and immune features in CESC. Additionally, we developed a risk model based on the dsRBP cluster, which demonstrated its potential as a prognostic predictor, which not only provided valuable prognostic information but also aided in determining the potential clinical benefits of chemotherapy and immunotherapy for CESC patients. Our findings highlight the importance of dsRBPs in CESC and shed light on their potential as biomarkers and therapeutic targets in the management of this disease.
Availability of data and materials
The datasets generated during and/or analysed during the current study are deposited in the National Genomics Data Center, China National Center for Bioinformation/Beijing Institute of Genomics, Chinese Academy of Sciences repository.
Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, 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(3):209–49. https://doi.org/10.3322/caac.21660.
Melnikow J, Henderson JT, Burda BU, Senger CA, Durbin S, Weyrich MS. Screening for cervical cancer with high-risk human papillomavirus testing: updated evidence report and systematic review for the US Preventive Services Task Force. JAMA. 2018;320(7):687–705. https://doi.org/10.1001/jama.2018.10400.
Molina MA, Coenen BA, Leenders WPJ, Andralojc KM, Huynen MA, Melchers WJG. Assessing the cervicovaginal microbiota in the context of hrHPV infections: temporal dynamics and therapeutic strategies. MBio. 2022;13: e0161922. https://doi.org/10.1128/mbio.01619-22.
Curry SJ, Krist AH, Owens DK, Barry MJ, Caughey AB, Davidson KW, et al. Screening for cervical cancer: US Preventive Services Task Force recommendation statement. JAMA. 2018;320(7):674–86. https://doi.org/10.1001/jama.2018.10897.
Cohen PA, Jhingran A, Oaknin A, Denny L. Cervical cancer. Lancet. 2019;393(10167):169–82. https://doi.org/10.1016/s0140-6736(18)32470-x.
Cancer Genome Atlas Research Network. Integrated genomic and molecular characterization of cervical cancer. Nature. 2017;543(7645):378–84. https://doi.org/10.1038/nature21386.
Abou Elela S, Ji X. Structure and function of Rnt1p: an alternative to RNAi for targeted RNA degradation. Wiley Interdiscip Rev RNA. 2019;10(3): e1521. https://doi.org/10.1002/wrna.1521.
Hur S, Double-Stranded RNA. Sensors and modulators in innate immunity. Annu Rev Immunol. 2019;37:349–75. https://doi.org/10.1146/annurev-immunol-042718-041356.
Wu Y, Wu X, Wu L, Wang X, Liu Z. The anticancer functions of RIG-I-like receptors, RIG-I and MDA5, and their applications in cancer therapy. Transl Res. 2017;190:51–60. https://doi.org/10.1016/j.trsl.2017.08.004.
Lee YS, Kunkeaw N, Lee YS. Protein kinase R and its cellular regulators in cancer: an active player or a surveillant? Wiley Interdiscip Rev RNA. 2020;11(2): e1558. https://doi.org/10.1002/wrna.1558.
Teoh PJ, Koh MY, Chng WJ. ADARs, RNA editing and more in hematological malignancies. Leukemia. 2021;35(2):346–59. https://doi.org/10.1038/s41375-020-01076-2.
Guo J, Lv J, Liu M, Tang H. miR-346 up-regulates argonaute 2 (AGO2) protein expression to augment the activity of other microRNAs (miRNAs) and contributes to cervical cancer cell malignancy. J Biol Chem. 2015;290(51):30342–50. https://doi.org/10.1074/jbc.M115.691857.
Chiappinelli Katherine B, Strissel Pamela L, Desrichard A, Li H, Henke C, Akman B, et al. Inhibiting DNA methylation causes an interferon response in cancer via dsRNA including endogenous retroviruses. Cell. 2015. https://doi.org/10.1016/j.cell.2015.07.011.
Sheng W, LaFleur MW, Nguyen TH, Chen S, Chakravarthy A, Conway JR, et al. LSD1 ablation stimulates anti-tumor immunity and enables checkpoint blockade. Cell. 2018;174(3):549-563.e519. https://doi.org/10.1016/j.cell.2018.05.052.
Ning H, Huang S, Lei Y, Zhi R, Yan H, Jin J, et al. Enhancer decommissioning by MLL4 ablation elicits dsRNA-interferon signaling and GSDMD-mediated pyroptosis to potentiate anti-tumor immunity. Nat Commun. 2022. https://doi.org/10.1038/s41467-022-34253-1.
Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinform. 2013;14:7. https://doi.org/10.1186/1471-2105-14-7.
Sabah A, Tiun S, Sani NS, Ayob M, Taha AY. Enhancing web search result clustering model based on multiview multirepresentation consensus cluster ensemble (mmcc) approach. PLoS ONE. 2021;16(1): e0245264. https://doi.org/10.1371/journal.pone.0245264.
Seiler M, Huang CC, Szalma S, Bhanot G. ConsensusCluster: a software tool for unsupervised cluster discovery in numerical data. OMICS. 2010;14(1):109–13. https://doi.org/10.1089/omi.2009.0083.
Rich JT, Neely JG, Paniello RC, Voelker CC, Nussenbaum B, Wang EW. A practical guide to understanding Kaplan–Meier curves. Otolaryngol Head Neck Surg. 2010;143(3):331–6. https://doi.org/10.1016/j.otohns.2010.05.007.
Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 2018;28(11):1747–56. https://doi.org/10.1101/gr.239244.118.
Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, 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.
Kobayashi Y, Kushihara Y, Saito N, Yamaguchi S, Kakimi K. A novel scoring method based on RNA-Seq immunograms describing individual cancer-immunity interactions. Cancer Sci. 2020;111(11):4031–40. https://doi.org/10.1111/cas.14621.
Bagaev A, Kotlov N, Nomie K, Svekolkin V, Gafurov A, Isaeva O, et al. Conserved pan-cancer microenvironment subtypes predict response to immunotherapy. Cancer Cell. 2021;39(6):845-865.e847. https://doi.org/10.1016/j.ccell.2021.04.014.
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47. https://doi.org/10.1093/nar/gkv007.
Hamilton NE, Ferry M. ggtern: ternary diagrams using ggplot2. J Stat Softw. 2018;87(Code Snippet 3):1–17.
Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30. https://doi.org/10.1093/nar/28.1.27.
Rodríguez-Álvarez MX, Meira-Machado L, Abu-Assi E, Raposeiras-Roubín S. Nonparametric estimation of time-dependent ROC curves conditional on a continuous covariate. Stat Med. 2016;35(7):1090–102. https://doi.org/10.1002/sim.6769.
Li C, Guo L, Li S, Hua K. Single-cell transcriptomics reveals the landscape of intra-tumoral heterogeneity and transcriptional activities of ECs in CC. Mol Ther Nucleic Acids. 2021. https://doi.org/10.1016/j.omtn.2021.03.017.
Sun D, Wang J, Han Y, Dong X, Ge J, Zheng R, et al. TISCH: a comprehensive web resource enabling interactive single-cell transcriptome visualization of tumor microenvironment. Nucleic Acids Res. 2020. https://doi.org/10.1093/nar/gkaa1020.
Sun D, Guan X, Moran AE, Wu LY, Qian DZ, Schedin P, et al. Identifying phenotype-associated subpopulations by integrating bulk and single-cell sequencing data. Nat Biotechnol. 2022;40(4):527–38. https://doi.org/10.1038/s41587-021-01091-3.
Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med. 2018;24(10):1550–8. https://doi.org/10.1038/s41591-018-0136-1.
Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, et al. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep. 2017;18(1):248–62. https://doi.org/10.1016/j.celrep.2016.12.019.
Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, et al. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature. 2018;554(7693):544–8. https://doi.org/10.1038/nature25501.
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(9): e107468. https://doi.org/10.1371/journal.pone.0107468.
Harden ME, Munger K. Perturbation of DROSHA and DICER expression by human papillomavirus 16 oncoproteins. Virology. 2017. https://doi.org/10.1016/j.virol.2017.04.022.
Wang SS, Gonzalez P, Yu K, Porras C, Li Q, Safaeian M, et al. Common genetic variants and risk for HPV persistence and progression to cervical cancer. PLoS ONE. 2010. https://doi.org/10.1371/journal.pone.0008667.
Madsen BE, Ramos EM, Boulard M, Duda K, Overgaard J, Nordsmark M, et al. Germline mutation in RNASEL predicts increased risk of head and neck, uterine cervix and breast cancer. PLoS ONE. 2008. https://doi.org/10.1371/journal.pone.0002492.
Girone C, Calati F, Lo Cigno I, Salvi V, Tassinari V, Schioppa T, et al. The RIG-I agonist M8 triggers cell death and natural killer cell activation in human papillomavirus-associated cancer and potentiates cisplatin cytotoxicity. Cancer Immunol Immunother. 2023. https://doi.org/10.1007/s00262-023-03483-7.
Chiang C, Pauli E-K, Biryukov J, Feister KF, Meng M, White EA, et al. The human papillomavirus E6 oncoprotein targets USP15 and TRIM25 To suppress RIG-I-mediated innate immune signaling. J Virol. 2018. https://doi.org/10.1128/jvi.01737-17.
Xing B, Guo J, Sheng Y, Wu G, Zhao Y. Human papillomavirus-negative cervical cancer: a comprehensive review. Front Oncol. 2020. https://doi.org/10.3389/fonc.2020.606335.
Chung H, Calis JJA, Wu X, Sun T, Yu Y, Sarbanes SL, et al. Human ADAR1 prevents endogenous RNA from triggering translational shutdown. Cell. 2018. https://doi.org/10.1016/j.cell.2017.12.038.
Li Y, Banerjee S, Goldstein SA, Dong B, Gaughan C, Rath S, et al. Ribonuclease L mediates the cell-lethal phenotype of double-stranded RNA editing enzyme ADAR1 deficiency in a human cell line. Elife. 2017. https://doi.org/10.7554/elife.25687.
Walkley CR, Li JB. Rewriting the transcriptome: adenosine-to-inosine RNA editing by ADARs. Genome Biol. 2017. https://doi.org/10.1186/s13059-017-1347-3.
Chen Y, Wang H, Lin W, Shuai P. ADAR1 overexpression is associated with cervical cancer progression and angiogenesis. Diagn Pathol. 2017. https://doi.org/10.1186/s13000-017-0600-0.
Nakamura K, Shigeyasu K, Okamoto K, Matsuoka H, Masuyama H. ADAR1 has an oncogenic function and can be a prognostic factor in cervical cancer. Sci Rep. 2023. https://doi.org/10.1038/s41598-023-30452-y.
Wilson RC, Doudna JA. Molecular mechanisms of RNA interference. Annu Rev Biophys. 2013. https://doi.org/10.1146/annurev-biophys-083012-130404.
Sheu-Gruttadauria J, MacRae IJ. Structural foundations of RNA silencing by argonaute. J Mol Biol. 2017. https://doi.org/10.1016/j.jmb.2017.07.018.
He L, Wang HY, Zhang L, Huang L, Li JD, Xiong Y, et al. Prognostic significance of low DICER expression regulated by miR-130a in cervical cancer. Cell Death Dis. 2014;5(5):e1205. https://doi.org/10.1038/cddis.2014.127.
Scotto L, Narayan G, Nandula SV, Subramaniyam S, Kaufmann AM, Wright JD, et al. Integrative genomics analysis of chromosome 5p gain in cervical cancer reveals target over-expressed genes, including Drosha. Mol Cancer. 2008;7:58. https://doi.org/10.1186/1476-4598-7-58.
Pan L, Xu C, Mei J, Chen Y, Wang D. Argonaute 3 (AGO3) promotes malignancy potential of cervical cancer via regulation of Wnt/β-catenin signaling pathway. Reprod Biol. 2021;21(1): 100479. https://doi.org/10.1016/j.repbio.2020.100479.
Zhao H, Jin X, Su H, Deng X, Fang Y, Shen L, et al. Down-regulation of Dicer expression in cervical cancer tissues. Med Oncol. 2014. https://doi.org/10.1007/s12032-014-0937-0.
Zhang Y, Zhang Z. The history and advances in cancer immunotherapy: understanding the characteristics of tumor-infiltrating immune cells and their therapeutic implications. Cell Mol Immunol. 2020;17(8):807–21. https://doi.org/10.1038/s41423-020-0488-6.
Li X, Zhang M, Lei T, Zou W, Huang R, Wang F, et al. Single-cell RNA-sequencing dissects cellular heterogeneity and identifies two tumor-suppressing immune cell subclusters in HPV-related cervical adenosquamous carcinoma. J Med Virol. 2022. https://doi.org/10.1002/jmv.28084.
Wang S, Li X, Liu C, Yuan Y, Ma F. Single-cell transcriptomic analysis of the role of HPV16-positive macrophages in cervical cancer prognosis. J Med Virol. 2022. https://doi.org/10.1002/jmv.28410.
Ohno A, Iwata T, Katoh Y, Taniguchi S, Tanaka K, Nishio H, et al. Tumor-infiltrating lymphocytes predict survival outcomes in patients with cervical cancer treated with concurrent chemoradiotherapy. Gynecol Oncol. 2020. https://doi.org/10.1016/j.ygyno.2020.07.106.
Wang Y, Wang C, Qiu J, Qu X, Peng J, Lu C, et al. Targeting CD96 overcomes PD-1 blockade resistance by enhancing CD8+ TIL function in cervical cancer. J Immunother Cancer. 2022. https://doi.org/10.1136/jitc-2021-003667.
Coffelt SB, Wellenstein MD, de Visser KE. Neutrophils in cancer: neutral no more. Nat Rev Cancer. 2016;16(7):431–46. https://doi.org/10.1038/nrc.2016.52.
Hedrick CC, Malanchi I. Neutrophils in cancer: heterogeneous and multifaceted. Nat Rev Immunol. 2022;22(3):173–87. https://doi.org/10.1038/s41577-021-00571-6.
Alexander J, Cukierman E. Cancer associated fibroblast: mediators of tumorigenesis. Matrix Biol. 2020;91–92:19–34. https://doi.org/10.1016/j.matbio.2020.05.004.
Yang D, Guo P, He T, Powell CA. Role of endothelial cells in tumor microenvironment. Clin Transl Med. 2021;11(6): e450. https://doi.org/10.1002/ctm2.450.
Monk BJ, Enomoto T, Kast WM, McCormack M, Tan DSP, Wu X, et al. Integration of immunotherapy into treatment of cervical cancer: recent data and ongoing trials. Cancer Treat Rev. 2022;106: 102385. https://doi.org/10.1016/j.ctrv.2022.102385.
Dyer BA, Zamarin D, Eskandar RN, Mayadev JM. Role of immunotherapy in the management of locally advanced and recurrent/metastatic cervical cancer. J Natl Compr Canc Netw. 2019;17(1):91–7. https://doi.org/10.6004/jnccn.2018.7108.
Chung HC, Ros W, Delord J-P, Perets R, Italiano A, Shapira-Frommer R, et al. Efficacy and safety of pembrolizumab in previously treated advanced cervical cancer: results from the phase II KEYNOTE-158 study. J Clin Oncol. 2019. https://doi.org/10.1200/jco.18.01265.
Louise F, Ken YL, Richard BSR, Chien-Fu H, Wu TC. Cervical cancer immunotherapy: facts and hopes. Clin Cancer Res. 2021. https://doi.org/10.1158/1078-0432.ccr-20-2833.
Shen F, Zheng H, Zhou L, Li W, Liu J, Xu X. Identification of CD28 and PTEN as novel prognostic markers for cervical cancer. J Cell Physiol. 2019;234(5):7004–11. https://doi.org/10.1002/jcp.27453.
He Y, Xiao M, Fu H, Chen L, Qi L, Liu D, et al. cPLA2α reversibly regulates different subsets of cancer stem cells transformation in cervical cancer. Stem Cells. 2020;38(4):487–503. https://doi.org/10.1002/stem.3157.
Martincorena I, Campbell PJ. Somatic mutation in cancer and normal cells. Science. 2015;349(6255):1483–9. https://doi.org/10.1126/science.aab4082.
Chen H, Zhao L, Liu J, Zhou H, Wang X, Fang X, et al. Bioinformatic analyzes of the association between upregulated expression of JUN gene via APOBEC-induced FLG gene mutation and prognosis of cervical cancer. Front Med. 2022;9: 815450. https://doi.org/10.3389/fmed.2022.815450.
Chen J, Deng Y, Ao L, Song Y, Xu Y, Wang CC, et al. The high-risk HPV oncogene E7 upregulates miR-182 expression through the TGF-β/Smad pathway in cervical cancer. Cancer Lett. 2019. https://doi.org/10.1016/j.canlet.2019.06.015.
French D, Belleudi F, Mauro MV, Mazzetta F, Raffa S, Fabiano V, et al. Expression of HPV16 E5 down-modulates the TGFbeta signaling pathway. Mol Cancer. 2013. https://doi.org/10.1186/1476-4598-12-38.
Liu Z, Wu M, Shi H, Huang C, Luo S, Song X. DDN-AS1-miR-15a/16-TCF3 feedback loop regulates tumor progression in cervical cancer. J Cell Biochem. 2019;120(6):10228–38. https://doi.org/10.1002/jcb.28307.
Zhang Z, Chen F, Li S, Guo H, Xi H, Deng J, et al. ERG the modulates Warburg effect and tumor progression in cervical cancer. Biochem Biophys Res Commun. 2020;522(1):191–7. https://doi.org/10.1016/j.bbrc.2019.11.079.
Ding H, Fan GL, Yi YX, Zhang W, Xiong XX, Mahgoub OK. Prognostic implications of immune-related genes’ (IRGs) signature models in cervical cancer and endometrial cancer. Front Genet. 2020;11:725. https://doi.org/10.3389/fgene.2020.00725.
Dong X, Chen C, Deng X, Liu Y, Duan Q, Peng Z, et al. A novel mechanism for C1GALT1 in the regulation of gastric cancer progression. Cell Biosci. 2021;11(1):166. https://doi.org/10.1186/s13578-021-00678-2.
Li H, Yang P, Wang J, Zhang J, Ma Q, Jiang Y, et al. HLF regulates ferroptosis, development and chemoresistance of triple-negative breast cancer by activating tumor cell-macrophage crosstalk. J Hematol Oncol. 2022;15(1):2. https://doi.org/10.1186/s13045-021-01223-x.
Nwachukwu CR, Mayadev J, Viswanathan AN. Concurrent chemoradiotherapy for stage IIIB cervical cancer-global impact through power. JAMA Oncol. 2018;4(4):514–5. https://doi.org/10.1001/jamaoncol.2017.5078.
Attademo L, Tuninetti V, Pisano C, Cecere SC, Di Napoli M, Tambaro R, et al. Immunotherapy in cervix cancer. Cancer Treat Rev. 2020;90: 102088. https://doi.org/10.1016/j.ctrv.2020.102088.
This study was sponsored by Shanghai Key Clinical Specialty Project (2019 LCZD-NXZL Shanghai Municipal Health Commission).
Ethics approval and consent to participate
This study was approved by the Ethics Committee of FUSCC and written informed consent was obtained from all the patients.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Identification of dsRBPs expression feature correlated with HPV infection in CESC samples from TCGA-CESC cohort. Comparison in the ssGSEA scores of each dsRBP subtype in CESC samples classified by HPV status (A), hierarchical HPV call (B), hierarchical HPV clade (C) and HPV integration status (D). dsRBPs: double-stranded RNA-binding proteins; ssGSEA: single sample Gene Set Enrichment Analysis. *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001; ns: non-significant. Figure S2. Tumor immunity related to dsRBPs expression patterns in cervical cancer. Analysis of the correlation between the expression seven dsRBPs subtypes and the level of immune score (A), tumor-infiltrated immune cells (B) and immune checkpoints (C). *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001. Figure S3. Boxplot showed the different expression levels of ADAR and DDR subfamily members among CESC samples with different HPV infection status in the TCGA-CESC cohort. ADAR: adenosine deaminases acting on RNA, DDR: Dicer, Drosha, and Argonautes. *p < 0.05; ns: non-significant. Figure S4. The correlation between clinical parameters and different clusters, such as age, T stage, N stage, M stage, and neoplasm disease stage. Figure S5. The distribution of HPV-infection patients in different dsRBP clusters. Figure S6. The Kaplan-Meier plots of overall survival stratified by age (≤ 60/> 60), T stage (T1–2/T3–4), N stage (N0/N1+), M stage (M0/M1), and neoplasm disease stage (stage I–II/stage III–IV). Figure S7. Univariate and multivariate Cox regression analysis of risk score and clinicopathological parameters. Figure S8. Comparison the of dsRBPs signature risk scores among CESC patients with different HPV infection status. *p < 0.05; **p < 0.01; ns: non-significant. Difference in the dsRBPs signature risk scores between CESC samples classified by HPV status (A), hierarchical HPV call (B), hierarchical HPV clade (C) and HPV integration status (D). Figure S9. Oncoplot of the genomic alterations in hub genes identified in the cervical cancer samples from the TCGA database.
Table S1. 35 dsRBPs colected and applied in current study. Table S2. Clinical characteristics of cervical cancer samples within different dsRBP clusters. Table S3. List of dsRBP cluster-associated DEGs. Table S4. List of prognostic dsRBP cluster-associated DEGs. Table S5. List of further achieved prognostic dsRBP cluster-associated DEGs after multistep filtering. Table S6. performance of signature in the TCGA-CESC cohort, GSE44001 cohort and CGCI-HTMCP-CC cohort. Table S7. Patients information stratified by the signature in the TCGA-CESC cohort. Table S8. Clinical feature information of the high- and low-risk group. Table S9. Sample information stratified by the signature in the GSE44001 validation cohort. Table S10. Sample information stratified by the signature in the CGCI-HTMCP-CC validation cohort. Table S11. Differentially expressed genes between Scissor+ and Scissor− cells.
About this article
Cite this article
Li, J., Wan, C., Li, X. et al. Characterization of tumor microenvironment and tumor immunology based on the double-stranded RNA-binding protein related genes in cervical cancer. J Transl Med 21, 647 (2023). https://doi.org/10.1186/s12967-023-04505-9