CESCProg: a compact prognostic model and nomogram for cervical cancer based on miRNA biomarkers

Introduction

Cervical cancer (cervical squamous cell carcinoma; CESC) ranks fourth globally among cancers in women, and second among women of reproductive age. Due to unequal implementation of invasive screening techniques, the morbidity and mortality rate of cervical cancer continues to rise in countries like India, where it accounted for 9.4% of all cancers and 18.3% of new cases in 2020 (Sung et al., 2021). Multiple etiological factors contribute to its incidence, including persistent infection of human papilloma virus (HPV) (Bruni et al., 2019), and known lifestyle factors such as excessive smoking and use of contraceptive pills. Cervical cancer tends to be refractory to treatment unless detected early, and its prognosis is vital to quality-of-life expectations. Late diagnoses in the advanced stages of cervical cancer require expensive and complex treatment, with concomitant poor prognoses (Mehrotra & Yadav, 2021). Many gaps remain with respect to cervical cancer screening, diagnosis and prognosis (Jiang et al., 2020), and biomarkers with high specificity and sensitivity are necessary.

MiRNAs exert key control over regulation of gene expression (Li & Rana, 2014), by inducing specific translational repression via target mRNA 3′ UTR deadenylation and decapping. MiRNAs are known to target ∼60% of the transcriptome, thus modulating biological processes (Pedroza-Torres et al., 2016). Their aberrant, differential expression is implicated in various cancers, where they act as either oncogenes (oncomirs) or tumor suppressor genes (mirsupps), regulating tumorigenic process like cell maturation, cell proliferation, migration, invasion, apoptosis, and metastasis (Reddy, 2015). MiRNA biomarkers from the serum or cervical mucus could potentially augment systems for early diagnosis, prediction of disease progression, and outcome improvement, in addition to facilitating prognostic information, with respect to cervical cancer. The US National Cancer Institute launched the Cancer Genome Atlas (TCGA) to characterize different tumor types using -omics platforms, and make raw and processed data available to all researchers (Chandran et al., 2016). In this study, we used the TCGA CESC miRNA -omics dataset to build a validated prognostic risk model and predictive nomogram based on a minimal miRNA signature and the clinical profile. The developed models have been deployed as a freely-available web-app service for non-commercial use at CESCProg (https://apalania.shinyapps.io/cescprog).

Materials and Methods

Portions of this text were previously published as a preprint (Muthamilselvan & Palaniappan, 2023). The workflow is summarised in Fig. 1, and discussed in detail below. R version ‘3.6.3’ (http://www.r-project.org) was used in the study.

Results

The TCGA expression data consisted of expression values of 2,589 miRNA in 312 samples enrolled in this study, including 309 cervical cancer tissues and three matched normal tissues. The dataset was pruned in accordance with the pre-processing steps outlined in the Methods section, yielding an expression dataset consisting of 467 miRNAs across 303 samples with stage annotation. The distribution of samples according to the AJCC staging system (Amin et al., 2017) is shown in Table 1. The demographic features and clinical characteristics considered, namely age, smoking history, vital status, pregnancies, HPV status, and histologic grade are summarized in Table 2. Fitting the linear model and applying the filter criteria yielded a total of 101 differentially expressed miRNAs between cervical cancer tissues and matched normal tissues, provided in File S1. Most of the top-ranked miRNAs are overexpressed (e.g., hsa-miR-200c-3p, hsa-miR-141, hsa-miR-200a, hsa-miR-21-5p), suggesting oncomir function with increased epigenetic suppression of target tumor-suppressor gene expression. Table 3 shows the top ten miRNAs with their stage-wise log₂FC and linear model significance.

Each DEmiR was subjected to univariate Cox modeling to evaluate its prognostic significance. This process identified only 52 miRNAs as significantly associated with overall survival, based on p-value <0.05 (data presented in File S2). To optimize the dimensions of the prognostic miRNA biomarker panel, we applied Lasso-penalized Cox regression on the 52 miRNAs to obtain five miRNAs, hsa-miR-625-5p, hsa-miR-3934-5p, hsa-miR-330-3p, hsa-miR-642a-5p, has-miR-95-3p. Only one miRNA, hsa-miR-616-5p, survived the SVM-RFE feature selection process. Figure 2 shows the union of these results (i.e, the six miRNAs). All the prognostic DEmiRS post feature selection processes were upregulated, but none were an outlier DEmiR.

The six miRNAs were taken forward for multivariate survival analysis, and subjected to a stepwise backward-selection process, to further compact the miRNA signature. This process yielded an optimal signature of three miRNAs namely hsa-miR-625-5p, hsa-miR-330-3p, and hsa-miR-95-3p, with model p-value <0.002 (Table 4), for construction of the CESC prognostic risk model.

The CESC prognostic risk model, given by Eq. (2), was then parameterized using the expression of these three miRNAs: (3)${\displaystyle \text{miRNA\_Risk\_score}=0.30\ast \text{hsa-miR-95-3p}+0.35\ast \text{hsa-miR-330-3p}-0.52\ast \text{hsa-miR-625-5p}.}$

It is seen that hsa-miR-625-5p has a significant protective effect on CESC OS, whereas the expression of hsa-miR-95-3p and hsa-miR-330-3p elevate the risk. Based on this model, we computed the risk score for each patient in the train dataset, and used the maxstat of the resulting risk-score distribution to separate patients into high- and low-risk groups (Fig. 3A). The Kaplan–Meier survival curve of this distribution revealed significantly worse prognosis in the high-risk group (p-value <1E-4) (Fig. 3B). Time-dependent ROC analysis of the risk-score model on the train dataset for 1-, 2-, 3-, and 5-year overall survival yielded prognostic AUC values of 0.71, 0.72, 0.74 and 0.73, respectively (Fig. 3C). These results encouraged validation of the CESC-related prognostic signature on the test dataset, whose risk-score distribution is shown in Fig. 4A. The following outcomes validated the results: (i) Kaplan–Meier survival curve showed significantly worse prognosis in the high-risk group ( p-value <1E-4) (Fig. 4B) ; and (ii) time-dependent AUROC values 0.84, 0.79, 0.71 and 0.71 were obtained for 1-, 2-, 3-, and 5-year overall survival, respectively (Fig. 4C).

To validate the prognostic value of the model on an external, out-of-domain dataset, we used the study results of How et al. (2015). This study used a TaqMan Low Density Array (TLDA) to measure expression in formalin-fixed paraffin-embedded (FFPE) cervix samples. Two datasets from the study were used:

(i) Normalized and log₂-transformed miRNA expression data of 87 FFPE cervix samples used for validation, available at: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4399941 and https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4399941. The clinical information was used to annotate the samples, and the expression subset corresponding to the three miRNAs in the optimal risk model (i.e., (3)) was extracted. Since the miRNA arm information(-3p or -5p) was missing for hsa-miR-625 and hsa-miR-95, the arm-neutral expression values for both these miRNAs were used. The risk score for each sample was calculated based on (3), and the resulting risk score distribution was stratified into high-risk and low-risk patient groups based on the maxstat statistic computed by R survminer. The curves were visualized using Kaplan–Meier analysis, yielding significantly worse prognosis (P = 0.032) in the high-risk patient group relative to the low-risk group (Fig. 5).

Certain clinical features namely age, HPV_status, pregnancies, smoking_history, histologic_grade, and stage could boost the prognostic predictive value, and hence were examined for candidate inclusion in the risk model. Each clinical feature was subjected to univariate Cox survival analysis. Surprisingly, neither the patient’s HPV-status nor tumor grade was significant to CESC OS, and moreover only one clinical feature turned out significant, namely the ‘Stage’ representing the distinction between early clinical stage (Stage: I/II) and late clinical stage (Stage: III/IV) cancers. This ‘Stage’ constitutes a risk factor independent of the miRNA-based risk-score, and was used with the miRNA_risk_score to construct an integrated multivariate Cox regression model. As a result, both factor levels of both the variables were significant, and the overall multivariate model was extremely significant (p-value ∼4E-05) (Fig. 6). The integrated CESC prognostic risk model was then parameterized as: (4)${\displaystyle \text{Integrated\_risk\_score}=0.64\ast \text{miRNA\_Risk\_score}+0.74\ast \text{Stage}.}$

Based on the risk models developed, a nomogram was built to predict one-year and five-year survival probabilities (Fig. 7). The nomogram C-index of the integrated model was estimated as 0.7136 ± 0.047, compared with a C-index of 0.6399 ± 0.044 for the miRNA-based risk model alone, thus differentiating the models on the basis of discrimination. Furthermore, the nomogram calibration plots of the integrated model for one-year and five-year OS probabilities based on bootstrap resampling showed consistency between the predicted and actual survival probabilities (Fig. 8).

Discussion

MiRNAs add a layer of critical regulatory control over genomic expression, and aberrations in their expression could lead to the development of cancer hallmarks (Iorio & Croce, 2012). MiRNAs could be detected in the serum, and lend valuable potential as diagnostic and prognostic biomarkers of various cancers, including cervical cancer (Pisarska & Baldy-Chudzik, 2020). Several prognostic miRNAs for CESC have been reported, including miR-31 (Wang et al., 2014), miR-155 (Fang et al., 2016), and miR425-5p (Sun et al., 2017). However a systematic hypothesis-free scan for comprehensive miRNA signatures remains missing in the literature. In this study, we have attempted to fill this void with an integrated multi-layered bioinformatics approach to the detection of a reliable prognostic DEmiR biomarker signature. The study has yielded three prognostic miRNAs, namely hsa-miR-625-5p, hsa-miR-95-3p, and hsa-miR-330-3p. Downregulation of hsa-miR-625-5p has been documented in many cancers including non-small cell lung cancer (Dao et al., 2020), hepatocellular carcinoma (Zhou et al., 2015), melanoma (Zou et al., 2019) and cervical cancer (Wang et al., 2020). A causal mechanism relating miR-625-5p expression to inhibition of cervical cancer cell growth via suppression of NF- κB signaling has been reported (Li et al., 2020), consistent with its mirsupp identity disclosed here. Sponging miR-625-5p in turn is likely to drive cervical cancer progression, and this has been demonstrated recently (Li et al., 2023). Jafarzadeh et al. suggested that miR-330-3p promoted pro-tumorigenic events in various cancers like those of the lung, pancreas, bladder, and cervix, and that its downregulation could stall tumor development (Jafarzadeh et al., 2022), both observations consistent with its oncomir identity disclosed here. Further, miR-95-3p has been implicated in activating the wnt/ βcatenin pathway in prostate cancer tissues (Xi et al., 2019), thereby promoting cell proliferation, migration and invasion, consistent with its oncomir identity disclosed here.

To examine the network-level effects of these miRNAs, we retrieved the RNA-Seq transcriptome for each patient in our dataset from firebrowse.org, and correlated this data with the expression of the three miRNAs of interest to infer potential target genes. Target genes with substantial inverse correlation in expression (defined as Pearson ρ or Spearman ρ or Kendall τ <−0.3) were identified, and the consensus with multiMiR version ‘1.8.0’ (Ru et al., 2014) predictions for each of the three miRNAs was investigated. This yielded three consensus target genes with respect to hsa-miR-95-3p, namely NXPH3, BOC, EID1; two consensus target genes with respect to hsa-miR-625-5p, namely SIN3B and TPRG1L; and two consensus target genes with respect to hsa-miR-330-3p, namely THRA and DYRK2. Functional enrichment analysis of the consensus genes conducted with miRNeT (Chang et al., 2020) on GO and KEGG databases yielded significance for cancer pathways and cell cycle regulation. We also used the miR2Trait server (Babu & Palaniappan, 2022) to investigate the diseasome of this three-miRNA signature, and found significance for ‘uterine cervical neoplasm’ (p-value ∼1.5E-3), ‘squamous cell carcinoma’ (p-value ∼7.7E-3), and ‘cervical intraepithelial neoplasia’ (p-value ∼2.2E-2). Detailed results of the above investigations are presented in File S3.

Though all stages of cervical cancer were sufficiently represented in the dataset we have used, it is noted that the dataset was imbalanced with respect to controls. To investigate if the observations in our study generalize when more controls were included in the analysis, we obtained an miRNA-based GEO (Edgar, Domrachev & Lash, 2002) study GSE104758 (Snoek et al., 2019), containing 32 controls and 24 carcinoma in situ CIN3 samples, and concerned with identifying biomarkers for detecting precancerous CIN3 lesions. The Sequence Read Archive of the 32 controls was imported into Galaxy (The Galaxy Community, 2022), checked for FastQC quality, and trimmed for adapters. The reads were then mapped with the reference human genome hg_38 to yield genome-wide miRNA quantification data. This dataset was batch-corrected using ComBat from SVA version ‘3.35.2’ (Leek et al., 2012), and used to augment the normals in the CESC TCGA cohort, yielding 35 normal and 300 cancer samples. The corresponding expression matrix was then modelled according to (1) in limma, as outlined in the Methods section. Using an adj. p-val threshold ∼0.05 yielded a total of 544 DE miRNAs, including all the 101 original DEmiRs and specifically the three miRNAs of the prognostic miRNA_Risk_score model derived above. The detailed results of the linear model obtained with this controls-augmented dataset are provided in File S4. Additional supporting evidence for our results were then obtained as follows: (i) A contrast between the normals and CIN3 samples of GSE104758 yielded 64 DEmiRs with p-value <0.05. This was compared with the nine miRNAs in ref. 44 obtained by applying group-regularised ridge regression. Except hsa-let-7b, the remaining eight miRNAs were common to both the analyses, viz. hsa-miR-9-5b, hsa-miR-15b-3p, hsa-miR-20a-5p, hsa-miR-31-5p, hsa-miR-93-5p, hsa-miR-183-5p, hsa-miR-184, and hsa-miR-222-3p. The GSE104758 study further narrowed down to a biomarker signature of just five miRNAs(viz. hsa-miR-93-5p, hsa-miR-15b-3p, and hsa-miR-222-3p, and hsa-let-7b) via multivariate logistic regression with backward selection. Application of SVM-RFE to the 64 DEmiRs obtained from our contrast analysis identified 24 DEmiRs, including all the miRNAs in the GSE104758 biomarker panel except hsa-let-7b. It is notable that the process also yielded hsa-miR-330-3p, ranked nine among the 24 selected features. (ii) Univariate Cox analysis with the batch-corrected cancer samples in the augmented dataset yielded 241 significant miRNAs, including hsa-miR-95-3p and hsa-miR-625-5p.

Nomograms are widely used for simplifying the task of interpretation from models, and have been constructed with miRNAs for cervical cancer screening (Kotani et al., 2022), prognosis (Liu et al., 2020), and recurrence risk (Bogani et al., 2019). To facilitate the ready prognosis of cervical cancer patients, the models developed in this work were re-built with the full (train + test) dataset, and served as a web-app named CESCProg, deployed at: https://apalania.shinyapps.io/cescprog/ for non-commercial uses. The concerned user may provide the form inputs, namely the expression values of the three prognostic DEmiRs and an optional sample staging information. Based on the user request, the app proceeds to classify the risk of the sample, and compute a risk-score based on (2) or (4). The calculated risk-score is then consulted with the back-end nomogram to estimate the one-year and five-year survival probabilities. Serum-based or cervical mucus-based miRNAs are minimally invasive, and could be detected and quantified using a range of techniques (e.g., see ref. Baabu et al., 2022).

Conclusions

MiRNA biomarkers are an emerging diagnostic and prognostic aid to the management of disease, especially cancers. Here we present CESCProg, an miRNA-based prognostic model for cervical cancer developed by applying a sequence of purifying filters to the TCGA CESC dataset. All three miRNAs in the panel, namely hsa-miR-95-3p, hsa-miR-330-3p and hsa-miR-625-5p, show upregulation in cervical cancer relative to controls, suggesting feasibility for detection as biomarkers. In the miRNA risk model, hsa-miR-625-5p exhibits a protective effect on OS, while the other two miRNAs elevate the risk. The miRNA risk model was effective and extremely significant in stratifying CESC OS on the test dataset. A second risk model was developed with the inclusion of clinical features to maximize nomogram discrimination. This yielded a C-index of 0.7136 ± 0.047. The models have been deployed as a web-service as a possible aid to medical decision-making. They are available for non-profit use at: https://apalania.shinyapps.io/cescprog.

Supplemental Information