Signature involved in immune-related lncRNA pairs for predicting the immune landscape of cervical cancer

Abstract

Background

Immune-related long non-coding RNAs (irlncRNAs) are known to hold great promise as superior biomarkers for cervical cancer-related immunotherapeutic response and the tumor immune microenvironment. Here, we constructed a prognostic signature based on irlncRNA pairs (IRLPs).

Methods

The samples were downloaded from The Cancer Genome Atlas and the Genotype-Tissue Expression databases. The least absolute shrinkage and selection operator Cox regression was performed to construct the prognostic model. Receiver operating characteristic (ROC) curve and nomogram were plotted to validate accuracy of the model. Next, we estimated the immune cell infiltration and the correlation between risk score and the expression of genes related to immune checkpoint. Finally, we calculated the score of the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm and the half maximal inhibitory concentration of the chemotherapeutic agent to evaluate the response to immunotherapy and chemotherapy.

Results

We constructed a prognostic signature that consisted of 11 irlncRNAs. The area under the curve values of the 1-, 3-, and 5-year ROC curves were 0.844, 0.891, and 0.871, respectively. The expression of CTLA-4, HAVCR2, IDO1, LAG3, and PDCD1 were negatively correlated with risk scores. The score of TIDE in the high-risk group was significantly higher than in the low-risk group (P < 0.01). Patients in the low-risk subgroup were more sensitive to chemotherapeutic agents, such as axitinib and docetaxel, whereas patients in the low-risk subgroup were more sensitive to mitomycin C.

Conclusion

Our study highlighted the value of the 11 IRLPs signatures to predict the prognosis and the response to immunotherapy and chemotherapeutics for patients with cervical cancer.

Introduction

Cervical cancer is the second leading deadly cancer among adults aged 30–39 years, and it has an incidence of more than 500,000 new cases, with almost 300,000 deaths per year worldwide. ¹ Importantly, advanced cervical cancer generally has a poor prognosis, with high rates of systemic recurrence. Accumulating studies confirmed that immunotherapy focused on immune checkpoints and tumor-infiltrating immune cells could provide tools to tailor a personalized treatment strategy for patients.^2,3 In the past decade, immune therapy has blossomed into fruition as a potential novel strategy. Consequently, novel and accurate immune-related biomarkers to predict prognosis for this patient population and to determine the personalized treatment warrant further investigation.

Long non-coding RNAs (lncRNAs) are defined as a class of transcripts with a length of over 200 nucleotides and a deficiency of protein-encoding potential. lncRNAs are involved in various biological processes. Indeed, recent discoveries have revealed the various functions of lncRNAs in cervical cancer, including cell proliferation, invasion, and apoptosis. ⁴ The evolving application of databases has rendered lncRNAs as non-invasive promising biomarkers in multiple cancers. ⁵ Unfortunately, issues such as batch effect and overfitting on data sets prevented biomarkers with immune-related lncRNAs (irlncRNAs) from becoming the standard clinical practice. In traditional technologies, the batch effect in the data analysis of gene expression levels was adjusted for the variability of inherent biological heterogeneity and technological differences, thus creating a tedious and complicated task. ⁶ A new data analysis technique, which has overcome the complexity for the utilization of data—such as normalization and scaling the expression matrix that merely depend on the relative ranking of gene expression values—was established and has already produced stable outcomes in various applications. ⁷ The technique just touched upon pairwise comparison in the gene expression values of an individual, and it could be easily translated to clinical practice.

In this study, The Cancer Genome Atlas (TCGA) and the Genotype-Tissue Expression GTEx databases were integrated to construct and validate a robust prognostic classifier for cervical cancer on the basis of irlncRNAs pairs (IRLPs). Surprisingly, the results showed that the risk model exhibited incredible prognostic ability, with potential in directing therapeutic strategy for cervical cancer.

Materials and methods

Data collection and differentially expressed irlncRNAs analysis

The samples of cervical cancer were downloaded from TCGA database and the GTEx database for subsequent analysis. RNA-sequencing transcriptome data for cervical cancer were received from the UCSC Xena database (https://xenabrowser.net/). The gene expression profiles were categorized as lncRNAs or messenger RNAs after adding the gene annotations data obtained from the Ensembl database (http://asia.ensembl.org). A co-expression analysis was used for the above-mentioned lncRNAs and immune-related genes obtained by the ImmPort database (http://www.immport.org). The genes were defined as irlncRNAs when the result met the correlation coefficient threshold absolute values of ＞0.4 and a P-value < 0.001. The differentially expressed (DE)irlncRNAs were analyzed using the “limma R” package, following which, the criteria of DEirlncRNAs between cancer and normal samples were set at the false discovery rate < 0.05 and | log2 fold change > 2. Clinical information regarding cervical cancer was searched by the UCSC Xena database, and patients with overall survival (OS) = 0 days were rejected to filter out unqualified samples.

Establishment and evaluation of a risk prediction model

The novel technique was established to build a 1-or-0 matrix by an interactive cycle to eliminate the bias of batch effect. The IRLP score was considered as 0, while the expression of IRLP1 was less than that of IRLP2 in a pairwise comparison. Conversely, the score was 1. If the proportion of IRLPs with a score of 1 or 0 in all samples exceeded 80%, it was defined as an invalid IRLP. Also, IRLPs with the value of P < 0.01 were defined as prognostic IRLPs through univariate Cox analysis. The least absolute shrinkage and selection operator (LASSO) Cox regression (iterations = 1000) was executed to prevent the model from over-fitting to screen candidate IRLPs. Depending on the IRLPs and their coefficient, a risk score was computed for patients as follows: risk score = coef (IRLP1) × expression (IRLP1) + … + coef (IRLPn) × expression (IRLPn). The Akaike information criterion (AIC) values of every point on the 5-year receiver operating characteristic (ROC) curve was computed to identify the optimal cutoff value, which was utilized to divide patients into high- and low-risk subgroups.

That risk score could be considered as an independent prognostic predictor independent of clinical and pathologic variables by using univariate and multivariate Cox regression analyses. Remarkably, a nomogram was established to further check the validity of the model. A time-dependent calibration curve and the decision curve were conducted to determine of the accuracy and the clinical usefulness of this nomogram.

Estimation of immune cell infiltration and the relationship between immune checkpoint and risk model

The correlation between infiltrating immune cells and the risk score was further explored using Spearman correlation analysis. Methods—including EPIC, TIMER, MCPcounter, XCELL, CIBERSORT, QUANTISEQ, and CIBERSORT-ABS—were conducted. The correlation between risk score and the expression of genes related to the immune checkpoints (ICIs) was profiled to investigate the potential role of the model in response to immunotherapy.

Evaluation of the response to immunotherapy and chemotherapy based on the model

We calculated the score of the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm to assess the response to immunotherapy for patients. In addition, the “pRRophetic” package was used to predict the the half maximal inhibitory concentration (IC50) of chemotherapeutic agents in different risk subgroups.

Statistical analysis

All statistical processing was implemented via R 4.0.4 software.

Results

Development and evaluation of a risk assessment model

Figure S1 illustrates the research steps. The transcriptome profiling data of cervical cancer was retrieved from TCGA and the GTEx database, including 306 tumor and 13 normal samples. Next, a total of 535 irlncRNAs were identified by co-expression analysis, and 103 DEirlncRNAs were distinguished (Figure 1(a)). These irlncRNAs underwent pairwise comparison, and 4554 valid differentially expressed IRLPs were extracted. Modified LASSO regression analysis was applied to predict the patient prognosis and to minimize the risk of overfitting, followed by univariate Cox regression analysis to select prognostic IRLPs (Figure S2). Ultimately, 11 IRLPs were obtained for subsequent model construction (Figure 1(b)).

Figure 1.

Construction of the risk assessment model by IRLPs. (a) Heatmap of DEirlncRNAs between tumor and normal samples. Upregulated genes are indicated in red, and downregulated genes are indicated in blue. (b) Forest map showing 11 prognostic IRLPs determined by multivariate Cox analysis. (c) All the AUC values of the ROC curve for predicting 1-, 3-, and 5-year overall survival time exceeded 0.8. (d) Patients in the low-risk group had better prognosis, as tested via the Kaplan–Meier analysis. (e) and (f) Distribution of risk scores (e) and survival status (f) of each patient.

The AUC values of the 1-, 3-, and 5-year ROC curves were 0.844, 0.891, and 0.871, respectively (Figure 1(c)). As shown in Figure 1(d), patients were dichotomized into low- and high-risk populations depending on the cutoff value, and extensive differences were found in the OS rate between the two subgroups. The survival time and survival status for patients were observed to continuously worsen, along with the increase in risk score (Figure 1(e) and 1(f)). In addition, the ROC curve of the IRLPs signature of the clinical and pathologic variables was drawn, and the AUC index of the IRLPs signature was found to be larger than that in the other variables. The differences in clinical-pathologic parameters between the two subgroups were compared (Figure 2(a)). The specific value was exhibited in Figure S3. The finding showed that body mass index, M stage, T stage, and clinical stage were significantly related to the model. Meanwhile, in the univariate Cox regression analysis (Figure 2(b)), only risk score and clinical stage were highly correlated with patient prognosis. After multivariable adjustment was conducted, risk score and clinical stage remained independent predictors for patient outcome (Figure 2(c)).

Figure 2.

Clinical evaluation of the risk model. (a) Heatmap and scatter plot presenting the relationship between risk score and clinical significance, including BMI, M stage, T stage, and clinical stage. (b) Forest plot of univariate Cox hazard ratio analysis showed that clinical stage and risk score were significantly associated with overall survival. (c) Multivariate Cox regression results demonstrated that clinical stage and risk score were independent prognostic predictors.

A nomograph for patients based on the risk score and clinical-pathologic variables was established to check the validity of the prediction (Figure 3(a)). In addition, a calibration plot was used to evaluate the predictive utility of the nomograph. Figure 3(b) unveiled that the calibration points in 1, 3, and 5 years had a high degree of overlap between estimated and observed outcomes. Decision curve analysis was performed to compare the benefit of different signatures (Figure 3(c)). Patients could benefit more from the risk model than either of the clinical-pathologic characteristics. Thus, this research offered a stable model applicable to evaluate patient prognosis, thereby increasing the probability of humans with cervical cancer receiving curative therapies.

Figure 3.

Construction and evaluation of the nomogram. (a) Nomogram predicting 1-, 3-, and 5-year overall survival for patients with cervical cancer based on IRLPs and clinicopathological parameters. (b) Calibration curves for 1-, 3-, and 5-years indicating excellent prediction of nomogram. (c) Decision curves of the nomogram based on IRLPs and other clinicopathological parameters.

Estimation of immune cell infiltration with the risk model

The tumor immune microenvironment has been reported had considerable progress in the characterization of the tumor microenvironment by tumor-infiltrating immune cells. Therefore, whether the model was associated with tumor-infiltrating immune cells was explored (Figure 4(a)). The specific value was exhibited in Table S1. The results of the differences in the content of various immune cells between low- and high-risk subgroups are available in Figure S3. They revealed more activated immune cells in the low-risk group, indicating the antitumoral immune responses of activated immune cells.

Figure 4.

Estimation of tumor immune cell infiltration and the relationship between ICIs and the risk model. (a) The high-risk score was negatively correlated with infiltrating immune cells, such as B cell, T-cell CD4 + , T-cell CD8 + , macrophage M1, macrophage M2, and myeloid dendritic cells, while positively relevant to Mast cell resting, NK cell resting, and cancer-associated fibroblast. (b) CTLA-4 (P < 0.01), LAG3 (P < 0.01), HAVCR2 (P < 0.001), IDO1 (P < 0.001), and PDCD1 (P < 0.001) had significantly low expression in the high-risk group. (c)The expression of CTLA-4(−0.12), HAVCR2(−0.16), IDO1(−0.12), LAG3(−0.11), and PDCD1(−0.16) were negatively correlated with risk scores.

Relationship between ICIs and risk model

The association of risk score with the expression of genes related to immunotherapy was preferentially analyzed to gain insights into the potential function of the risk model in response to immunotherapy. First, apart from CD274, a downregulation of CTLA-4, HAVCR2, IDO1, LAG3, and PDCD1 expression in the high-risk subgroup was observed (Figure 4(b)). In particular, the expression levels of CTLA-4, HAVCR2, IDO1, LAG3, and PDCD1 were negatively correlated with risk scores (Figure 4(c)).

Association between immunotherapy, chemotherapeutics response and the model

Finally, the correlation between risk score and the response to immunotherapy and chemotherapeutic agents was studied to further explore the clinical application value of the model. Patients in the high-risk subgroup had a higher TIDE score, which indicated a greater possibility of immune escape (Figure 5(a)). Furthermore, the higher T-cell exclusion score (Figure 5(b)) suggested that the high-risk population had a poor immunotherapy effect. The results showed that patients in the low-risk subgroup were more sensitive to chemotherapeutic agents, such as axitinib (Figure 5(c)) and docetaxel (Figure 5(d)), whereas patients in the low-risk subgroup were more sensitive to mitomycin C (Figure 5(e)). The findings supported that the model had satisfactory predictive value for response to immunotherapy and chemotherapeutics.

Figure 5.

Relationship between immunotherapeutic response, chemosensitivity, and the risk model. (a) The score of TIDE in the high-risk group was significantly higher than that in the low-risk group (P < 0.01). (b) The score of T-cell exclusion in the high-risk group was significantly higher than that in the high-risk group (P < 0.05). The box plot indicated that a high-risk score was positively associated with lower IC50 for chemotherapeutics such as axitinib (c) and docetaxel (d). In contrast, a low-risk score was positively associated with lower IC50 for and mitomycin C (e).

Discussion

Cancer immunotherapies that target the programmed death-1/programmed death-ligand 1 (PD-L1) and cytotoxic T lymphocyte antigen 4 (CTLA4) immune checkpoint pathway have been recognized as an emerging hallmark to date.^8,9 Although preliminary efficacies of immunotherapy have rejuvenated the treatment of cervical cancer, substantial subsets of patients failed to respond. ¹⁰ The former research focused on gene expression-based prognostic signatures to assess cervical cancer risk, mainly relying upon quantifying the expression levels of transcripts.^11–13 The drawback of these methods was that the diversity in detection methods and platforms could cause fluctuations in the expression value. Notably, unlike conventional approaches, a prognostic signature with two-lncRNA conjugates was developed in this study through the relative order of gene expression rather than specifically their expression.

As previously described, a unique method was adopted to set up a 1-or-0 matrix by an interactive cycle to solve the difficulty of batch calibration. In this work, a prognostic signature based on 11 DEirlncRNAs pairs for cervical cancer was developed. Many studies have been carried out already on some of the DEirlncRNAs in the model; these were found to play crucial roles during the process of tumor occurrence and the development of various tumors. For example, the previous experimental results showed that SOX21-AS1 and MIR9-3HG promoted the proliferation and migration of cervical cancer cells and inhibited apoptosis, respectively.^14,15 They believed that SOX21-AS1 and MIR9-3HG could act as efficient molecular targets and potential biomarkers for cervical cancer, in agreement with the results of the present study. In addition, the DEirlncRNAs played a regulatory role in other cancers, such as LINC02323, ¹⁶ AL606489.1, ¹⁷ LINC00944. ¹⁸ These studies further shed light on the authenticity of the risk model.

Moreover, the AUC value of the ROC curve for 1-, 3-, 5-year survival was much higher than that in other prognostic models established by other researchers for cervical cancer, which indicates that the model was appropriate for predicting individualized outcomes.^11–13 Importantly, after constructing the model by using innovative upgradation, it was reevaluated under diverse settings as follows: survival, clinical-pathologic variables, immune cell infiltration, ICIs, and response to immunotherapy and chemotherapeutics.

The infiltrating immune cells in the tumor microenvironment strongly influenced the clinical outcomes for immunotherapy. ¹⁹ We explored the composition of infiltrating tumor cells, and propose that the personalized therapeutics of adoptive T-cell therapy has clinical activity in the complete regression of advanced or recurrent metastatic cervical cancer. After stimulating angiogenesis, boosting cell invasion, migration, and infiltration of macrophages suppressed antitumor immunity, thereby facilitating occurrence and progression in most cancers. ²⁰ The main premise for inducing antitumor immune responses is destroying negative immune checkpoints, which has fueled a wave of development of ICIs. ²¹ Particularly noteworthy was the downregulation of CTLA-4, HAVCR2, IDO1, LAG3, and PDCD1 expression in the high-risk group, but not for CD274. The reason why CD274 expression had no difference between the two groups may be that the response of anti-PD-L1 has been limited in several cancers. To a great extent, the results of the PD-L1 test also varies among analytical conditions, scoring algorithms, and the location and timing of tissue collection. ²² The model in the present study found that axitinib and docetaxel were more sensitive to patients with cervical cancer in the high-risk subgroup, while mitomycin C may be more effective in the low-risk subgroup. These results had considerable implications for the manner that solidifies and reinforces the clinical applicability of the model.

In conclusion, a risk model without measuring specific gene expression was established to predict the prognosis for women with cervical cancer and their response to immunotherapy. As immunotherapy is a promising strategy, this model harbors direct translational influence for the treatment and prognosis of cervical cancer.

Supplemental Material

sj-docx-1-jbm-10.1177_03936155221091832 - Supplemental material for Signature involved in immune-related lncRNA pairs for predicting the immune landscape of cervical cancer

Supplemental material, sj-docx-1-jbm-10.1177_03936155221091832 for Signature involved in immune-related lncRNA pairs for predicting the immune landscape of cervical cancer by Xueting Liu, Zhao Wang, Le Wang, Ying Wang, Yuan Wang, Shanshan Yang and Yunyan Zhang in The International Journal of Biological Markers

Supplemental Material

sj-xls-2-jbm-10.1177_03936155221091832 - Supplemental material for Signature involved in immune-related lncRNA pairs for predicting the immune landscape of cervical cancer

Supplemental material, sj-xls-2-jbm-10.1177_03936155221091832 for Signature involved in immune-related lncRNA pairs for predicting the immune landscape of cervical cancer by Xueting Liu, Zhao Wang, Le Wang, Ying Wang, Yuan Wang, Shanshan Yang and Yunyan Zhang in The International Journal of Biological Markers