# Immune-Related Genes: Potential Prognostic Factors and Regulatory Targets for Cervical Carcinoma

Authors: Feng Xie; Xiaofeng Liu; Hua Liu; Min Wei; Wei Liu
Journal: Journal of Nanomaterials
 (2021)

---

## Abstract

Background. Advanced cervical carcinoma carries a particularly poor prognosis, and few treatment options exist. It is very important to find a method to evaluate the prognosis and survival rate of cervical carcinoma. The metastasis and invasion of cervical carcinoma are closely related to tumor immune microenvironment (TIME), and immune related genes (IRGs) are involved in the regulation of TIME. However, the role of IRGs in the prognosis of patients with cervical carcinoma remains unclear. Methods. The gene expression profiles of cervical carcinoma were downloaded from The Cancer Genome Atlas (TCGA) database, and IRG information were obtained from the ImmPort database. The IRGs were screened by coexpression analysis and were also performed function enrichment and pathway analyses. A prognosis model was built based on IRGs, and the risk score (RS) was calculated by Cox regression analysis. The accuracy was assessed by receiver operating characteristic (ROC) curve analysis. Besides, the relationship between RS and TIMER-generating immune cell content was performed by immune infiltration analysis. Results. In a total of 2503 differentially expressed genes (DEGs), 204 genes were IRGs, 20 of which were crucially correlated with the survival rate of cervical carcinoma. On the basis of Cox regression analysis, 6 IRGs were included in the prognosis model to calculate the RS. Kaplan-Meier survival and ROC analyses showed that the prognostic function of the model was superior to the current model constructed by clinicopathological risk factors. In addition, these 6 IRG signatures were related to the immune infiltration levels of six immune cells and the overall survival (OS) of cervical carcinoma. Finally, C-terminal Src kinase (CSK) gene is related to tumor metastasis, and Slit guidance ligand 2 (Slit2) is related to tumor clinical stage. Conclusion. The IRGs may contribute to the stratification of prognosis, and CSK/Slit2 may be two suppressor genes for cervical carcinoma.


---

## Body

## 1. Introduction

Cervical carcinoma is the fourth commonest female malignancy, which has become a thorny public health problem worldwide [1]. Although comprehensive treatment including surgery, radiotherapy, chemotherapy, and targeted therapy has achieved great success, the prognosis of advanced and recurrent cervical carcinoma is still poor [2–4].Recently reports showed that the metastasis and invasion of cervical carcinoma are closely related to the tumor immune microenvironment (TIME), and targeted therapy for its TIME has been more and more developed and applied [5–8]. The cervical carcinoma progresses in TIME mainly composed of infiltrating immune and stromal cells. Immune-related gene (IRG) is actively in the process of immune activity [9], especially for immune cell infiltration [10]. Nevertheless, the role of IRGs in the prognosis of patients with cervical carcinoma remains unclear.The traditional histopathological prognostic factors include age, tumor-node-metastasis (TNM) staging system, disease grade, histological subtype, and vascular invasion [11]. However, with the support of new technologies, gene abnormalities in TCGA provide additional, potential, and better prognostic information, which may replace the traditional histopathological prognostic factors [12]. IRGs act as a gene spectrum that actively participates in the process of TIME [13], which may play a more important role in tumor prognosis.In this study, we performed a systematic investigation of the IRGs associated with the prognosis of cervical carcinoma and explored the related prognostic biomarkers. We established a prognostic risk scoring system based on six IRG expressions and demonstrated that it was superior to the currently used clinicopathological risk factors. In addition, we found that the risk scoring system is related to immune cell infiltration and overall survival (OS) of cervical carcinoma. Finally, we identified two tumor suppressor genes as prognostic biomarkers for cervical carcinoma. This study provides a comprehensive analysis of the relationship between the prognosis of cervical carcinoma and IRGs, which may provide novel biomarkers for the prognosis of cervical carcinoma, and it is important for the identification of therapeutic targets of cervical carcinoma.


## 2. Materials and Methods

### 2.1. Data Acquisition and Processing

Transcriptomic RNA-sequencing data and clinical characteristics of 286 cervical carcinoma samples (analyzed with Illumina HiSeq 2000) and 3 normal samples were downloaded from The Cancer Genome Atlas (TCGA) database (https://cancergenome.nih.gov/). The workflow type was fragments per kilobase million (FPKM). The Ensemble IDs of mRNAs in cervical carcinoma were extracted from the GENCODE project (http://www.gencodegenes.org). A list of IRGs for cancer researches was obtained via the ImmPort database (http://www.immport.org). The infiltration data of 6 kinds of immune cells (B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells) were obtained from the TIMER database (https://cistrome.shinyapps.io/timer/). The flow chart of this study is shown in Figure 1.Figure 1
The flow chart of the study.Firstly, the expression profile data of 286 primary cervical carcinoma samples and 3 normal samples was downloaded from TCGA database. Coexpression analysis found the differentially expressed IRGs between primary cervical carcinoma samples and normal samples. The potential molecular mechanism of these IRGs was explored via the Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. In addition, a prognostic risk scoring system was established based on 6 IRG expressions, which were related to immune cell infiltration and overall survival (OS) of cervical carcinoma. Finally, the expression of two tumor suppressor genes (CSK and slit2) in tumors is shown in the TIMER database.


### 2.2. Differentially Expressed IRG Analysis

Differentially expressed IRGs between cervical carcinoma and normal samples were extracted via the limma package. We performed an analysis of the differentially expressed genes (DEGs) under the condition of a false discovery rateFDR<0.05 and log∣foldchange∣>1. Differentially expressed IRGs were then screen out from all differentially expressed genes based on the condition of a ∣Pearsoncorrelationcoefficient∣>0.4 and P value < 0.001. The potential molecular mechanisms of the differentially expressed IRGs were explored by functional enrichment analyses via the GO and KEGG pathways.


#### 2.2.1. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Enrichment Analysis

GO is a major bioinformatics tool for annotating genes and analyzing biological processes (BPs), cellular components (CCs), and molecular functions (MFs), respectively [14]. KEGG is a database resource for understanding advanced functions and biological systems from large-scale molecular datasets generated by high-throughput experimental technologies [15]. We used the R package clusterProfiler [16] for enrichment analysis of genes, and a P value < 0.05 was considered statistically significant.


### 2.3. Establishment and Validation of the Prognostic Risk Scoring System

#### 2.3.1. Screening for Prognosis-Related IRGs

To screen prognosis-related IRGs of cervical carcinoma, we performed univariate Cox proportional hazard regression using the survival package in R. Identified genes were considered statistically significant with aP value < 0.05.


#### 2.3.2. The Prognostic Risk Scoring System Establishment

A risk scoring signature was built by multivariate Cox proportional hazard regression to select the genes most closely related to survival. First, we obtained the regression coefficient of each gene using the survival package in R. The coefficient of each selected gene represented the estimated logarithm of the hazard ratio (HR, the parameter exp (coef in R)). Then, a risk score formula was established for all patients. The signature is expressed as follows:riskscore=coefficientgene1×gene1expression+coefficientgene2×expressionofgene2+⋯+coefficientgenen×expressiongenen.


#### 2.3.3. Prognostic Risk Scoring System Validation

The sensitivity of this system was analyzed with time-dependent receiver operating characteristic (ROC) curve. The area under curve (AUC) was determined to predict the 1-year, 3-year, and 5-year survivals using the survival ROC package in R. The median risk score served as a cutoff value to classify patients into the high-risk and low-risk groups. A subgroup analysis was performed by dividing the patients based on the risk score. Kaplan-Meier curves were used to assess the survival difference between two groups with a log-rankP value < 0.05.The specificity of this system was analyzed with a nomogram. The risk score and clinical factors were chosen to construct a nomogram to investigate the probability of 1-year, 3-year, and 5-year overall survivals (OS) of patients with cervical carcinoma. A calibration curve was plotted to assess the discrimination and accuracy of the nomogram.


### 2.4. Inference of Infiltrating Cells in the TIMER

The infiltration data of six immune cells (B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells) were obtained from the TIMER database (https://cistrome.shinyapps.io/timer/). The correlation between risk scores and immune infiltration levels of six immune cells was calculated using the Pearson analysis. And the relationship between the immune infiltration levels of six immune cells and the cumulative survival of cervical carcinoma was also investigated.


### 2.5. Survival Analysis

Survival-associated IRGs were selected by univariate and multivariate Cox proportional hazard regression analyses, which were conducted using the R software survival package.


### 2.6. Statistical Analysis

All statistical analyses were performed using R v3.4.1 (https://www.r-project.org/), SPSS 16.0 (SPSS Inc., Chicago, IL), and Prism 7 (GraphPad Software Inc., La Jolla, CA). All values are represented as the means±SEM of the indicated number of measurements.


## 2.1. Data Acquisition and Processing

Transcriptomic RNA-sequencing data and clinical characteristics of 286 cervical carcinoma samples (analyzed with Illumina HiSeq 2000) and 3 normal samples were downloaded from The Cancer Genome Atlas (TCGA) database (https://cancergenome.nih.gov/). The workflow type was fragments per kilobase million (FPKM). The Ensemble IDs of mRNAs in cervical carcinoma were extracted from the GENCODE project (http://www.gencodegenes.org). A list of IRGs for cancer researches was obtained via the ImmPort database (http://www.immport.org). The infiltration data of 6 kinds of immune cells (B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells) were obtained from the TIMER database (https://cistrome.shinyapps.io/timer/). The flow chart of this study is shown in Figure 1.Figure 1
The flow chart of the study.Firstly, the expression profile data of 286 primary cervical carcinoma samples and 3 normal samples was downloaded from TCGA database. Coexpression analysis found the differentially expressed IRGs between primary cervical carcinoma samples and normal samples. The potential molecular mechanism of these IRGs was explored via the Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. In addition, a prognostic risk scoring system was established based on 6 IRG expressions, which were related to immune cell infiltration and overall survival (OS) of cervical carcinoma. Finally, the expression of two tumor suppressor genes (CSK and slit2) in tumors is shown in the TIMER database.


## 2.2. Differentially Expressed IRG Analysis

Differentially expressed IRGs between cervical carcinoma and normal samples were extracted via the limma package. We performed an analysis of the differentially expressed genes (DEGs) under the condition of a false discovery rateFDR<0.05 and log∣foldchange∣>1. Differentially expressed IRGs were then screen out from all differentially expressed genes based on the condition of a ∣Pearsoncorrelationcoefficient∣>0.4 and P value < 0.001. The potential molecular mechanisms of the differentially expressed IRGs were explored by functional enrichment analyses via the GO and KEGG pathways.


### 2.2.1. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Enrichment Analysis

GO is a major bioinformatics tool for annotating genes and analyzing biological processes (BPs), cellular components (CCs), and molecular functions (MFs), respectively [14]. KEGG is a database resource for understanding advanced functions and biological systems from large-scale molecular datasets generated by high-throughput experimental technologies [15]. We used the R package clusterProfiler [16] for enrichment analysis of genes, and a P value < 0.05 was considered statistically significant.


## 2.2.1. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Enrichment Analysis

GO is a major bioinformatics tool for annotating genes and analyzing biological processes (BPs), cellular components (CCs), and molecular functions (MFs), respectively [14]. KEGG is a database resource for understanding advanced functions and biological systems from large-scale molecular datasets generated by high-throughput experimental technologies [15]. We used the R package clusterProfiler [16] for enrichment analysis of genes, and a P value < 0.05 was considered statistically significant.


## 2.3. Establishment and Validation of the Prognostic Risk Scoring System

### 2.3.1. Screening for Prognosis-Related IRGs

To screen prognosis-related IRGs of cervical carcinoma, we performed univariate Cox proportional hazard regression using the survival package in R. Identified genes were considered statistically significant with aP value < 0.05.


### 2.3.2. The Prognostic Risk Scoring System Establishment

A risk scoring signature was built by multivariate Cox proportional hazard regression to select the genes most closely related to survival. First, we obtained the regression coefficient of each gene using the survival package in R. The coefficient of each selected gene represented the estimated logarithm of the hazard ratio (HR, the parameter exp (coef in R)). Then, a risk score formula was established for all patients. The signature is expressed as follows:riskscore=coefficientgene1×gene1expression+coefficientgene2×expressionofgene2+⋯+coefficientgenen×expressiongenen.


### 2.3.3. Prognostic Risk Scoring System Validation

The sensitivity of this system was analyzed with time-dependent receiver operating characteristic (ROC) curve. The area under curve (AUC) was determined to predict the 1-year, 3-year, and 5-year survivals using the survival ROC package in R. The median risk score served as a cutoff value to classify patients into the high-risk and low-risk groups. A subgroup analysis was performed by dividing the patients based on the risk score. Kaplan-Meier curves were used to assess the survival difference between two groups with a log-rankP value < 0.05.The specificity of this system was analyzed with a nomogram. The risk score and clinical factors were chosen to construct a nomogram to investigate the probability of 1-year, 3-year, and 5-year overall survivals (OS) of patients with cervical carcinoma. A calibration curve was plotted to assess the discrimination and accuracy of the nomogram.


## 2.3.1. Screening for Prognosis-Related IRGs

To screen prognosis-related IRGs of cervical carcinoma, we performed univariate Cox proportional hazard regression using the survival package in R. Identified genes were considered statistically significant with aP value < 0.05.


## 2.3.2. The Prognostic Risk Scoring System Establishment

A risk scoring signature was built by multivariate Cox proportional hazard regression to select the genes most closely related to survival. First, we obtained the regression coefficient of each gene using the survival package in R. The coefficient of each selected gene represented the estimated logarithm of the hazard ratio (HR, the parameter exp (coef in R)). Then, a risk score formula was established for all patients. The signature is expressed as follows:riskscore=coefficientgene1×gene1expression+coefficientgene2×expressionofgene2+⋯+coefficientgenen×expressiongenen.


## 2.3.3. Prognostic Risk Scoring System Validation

The sensitivity of this system was analyzed with time-dependent receiver operating characteristic (ROC) curve. The area under curve (AUC) was determined to predict the 1-year, 3-year, and 5-year survivals using the survival ROC package in R. The median risk score served as a cutoff value to classify patients into the high-risk and low-risk groups. A subgroup analysis was performed by dividing the patients based on the risk score. Kaplan-Meier curves were used to assess the survival difference between two groups with a log-rankP value < 0.05.The specificity of this system was analyzed with a nomogram. The risk score and clinical factors were chosen to construct a nomogram to investigate the probability of 1-year, 3-year, and 5-year overall survivals (OS) of patients with cervical carcinoma. A calibration curve was plotted to assess the discrimination and accuracy of the nomogram.


## 2.4. Inference of Infiltrating Cells in the TIMER

The infiltration data of six immune cells (B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells) were obtained from the TIMER database (https://cistrome.shinyapps.io/timer/). The correlation between risk scores and immune infiltration levels of six immune cells was calculated using the Pearson analysis. And the relationship between the immune infiltration levels of six immune cells and the cumulative survival of cervical carcinoma was also investigated.


## 2.5. Survival Analysis

Survival-associated IRGs were selected by univariate and multivariate Cox proportional hazard regression analyses, which were conducted using the R software survival package.


## 2.6. Statistical Analysis

All statistical analyses were performed using R v3.4.1 (https://www.r-project.org/), SPSS 16.0 (SPSS Inc., Chicago, IL), and Prism 7 (GraphPad Software Inc., La Jolla, CA). All values are represented as the means±SEM of the indicated number of measurements.


## 3. Results

### 3.1. Identification of Differentially Expressed IRGs

We downloaded the expression profile data of cervical carcinoma through TCGA database, including 286 primary cervical carcinoma samples and 3 normal samples. Through DEG analysis, we identified 2503 DEGs in cervical carcinoma. And 204 DEGs are IRGs (Figure2(a)), of which 56 were upregulated and 148 were downregulated (Figure 2(b)).Figure 2
Identification of differentially expressed IRGs and gene functional enrichment. (a) Heatmap shows the expression of IRGs in cervical carcinoma tumor tissues and normal tissues. N: normal tissue; T: cervical carcinoma tissue. (b) Volcano plots demonstrating differentially expressed IRGs between cervical carcinoma and normal tissues, in the condition of afalsediscoveryrateFDR<0.05 and a log∣foldchangeFC∣>1. Red dots: upregulated genes, green dots: downregulated genes. (c) Gene Ontology (GO) analysis of differentially expressed IRGs. The top 10 most significant pathways of biological process (BP), cellular component (CC), and molecular function (MF) were shown in the training set, respectively. P value < 0.05 was considered statistically significant. (d) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of differential expressed IRGs. P value < 0.05 was considered statistically significant.
(a)(b)(c)(d)The potential molecular mechanisms of these 204 IRGs were explored by functional enrichment analyses via the GO and KEGG pathways. The GO analysis shows the top 10 clusters of enriched sets (Figure2(c)). The top 10 BPs of these IRGs were mainly including peptidyl-tyrosine phosphorylation/modification, epithelial cell proliferation, protein kinase B signaling, cell chemotaxis, and the MAPK cascade. Moreover, these IRGs were expressed in CC of platelets, extracellular matrix, vesicle lumen, granule lumen, receptor complex, side of membrane, and so on. Specifically, the MF of these IRGs is mainly related to receptor ligand/regulator activity, growth factor activity/binding, cytokine binding, transmembrane receptor protein 1 kinase activity, transforming growth factor beta binding, steroid hormone receptor activity, glycosaminoglycan binding, and nuclear receptor activity. Subsequently, the KEGG pathway results showed the top 30 significantly upregulated pathways, including cytokine-cytokine receptor interaction, neuroactive ligand-receptor interaction, Ras signaling pathway, JAK-STAT signaling pathway, MAPK signaling pathway, Rap1 signaling pathway, and PI3K-Akt signaling pathway (Figure 2(d)).


### 3.2. The Prognostic Risk Scoring System Establishment

Based on univariate Cox regression analysis, 20 IRGs with significant prognostic values were extracted with aP value < 0.05 (Figure 3(a)). Multivariate Cox proportional hazard regression analysis to further select six IRGs (CSK, SLIT2, GH1, IL3RA, LEPR, and PRLHR) to construct a risk scoring system under the condition of AIC=600.15 (Figure 3(b)). Information on these six IRGs is shown in Table 1, and the univariate/multivariate Cox regression results are shown in Figures 3(a) and 3(b). ROC curve analysis showed that the risk scoring system had a moderate accuracy in predicting the 3-year survival rate of cervical carcinoma (AUC=0.782) (Figure 3(c)). There was a significant correlation between the risk score and patient OS or mortality. With the increase in risk score, the survival prognosis worsened, and the mortality increased (Figure 3(d)). By dividing the sample into the high-risk and low-risk groups according to the median risk score in the model, we found a significant difference in survival rate (P value = 3.588e−06; Figure 3(e)) and clinical characteristics (Table 2) between the high- and low-risk groups. Further statistics showed that there were significant differences in age and fustat between the two groups (P value < 0.05), but no significant difference in stage (stages I~IV) (P value > 0.05, Figure 3(f)). These results indicate that the risk score system may be another prediction system independent of the TNM staging system.Figure 3
Survival-associated IRGs and prognosis risk scoring system construction. (a) Univariate Cox proportional hazard regression analysis in cervical carcinoma showing the prognostic values of IRGs. The regression coefficient of each gene was obtained using the survival package in R. The coefficient (parameter coef in R) of each selected gene represented the estimated logarithm of the Hazard Ratio (HR, parameter exp (coef in R)).P value < 0.05 was considered significant. (b) A risk scoring system based on 6 IRGs was built by multivariate Cox proportional hazard regression under the condition of AIC=600.15. The signature is expressed as follows: riskscore=coefficientgene1×gene1expression+coefficientgene2×expressionofgene2+⋯+coefficientgenen×expressiongenen. ∗P value < 0.05, ∗∗P value < 0.01, and ∗∗∗P value < 0.001. (c) Survival-dependent receiver operating characteristic (ROC) curve validation of the prognostic value of the risk score. (d) Rank of risk score and distribution of groups or the survival status of patients. Red: dead outcome; green: alive outcome. (e) Kaplan-Meier curves were used to assess the survival difference between two groups using the log-rank test. The median risk score served as a cutoff value to classify patients into the high-risk and low-risk groups. Red: high-risk group; purple: low-risk group. (f) Comparison of clinicopathological features and expression profiles of 6 IRGs in different groups. ∗P value < 0.05, ∗∗P value < 0.01, and ∗∗∗P value < 0.001.
(a)(b)(c)(d)(e)(f)Table 1
Information on six immune-related genes.
Gene symbolHGNC symbolFull nameCSK2444C-terminal Src kinaseSLIT211086Slit guidance ligand 2GH14261Growth hormone 1IL3RA6012Interleukin 3 receptor subunit alphaLEPR6554Leptin receptorPRLHR4464Prolactin releasing hormone receptorTable 2
The clinical characteristics of the cervical carcinoma patients in the two groups.
VariablesHigh riskLow risk(n=135)(n=136)Age (years)≤65115120>651712Unknown34Tumor stageStage I7174Stage II2633Stage III2120Stage IV145Unknown34FustatDead47117Alive8819


### 3.3. Evaluation of the Accuracy of the Prognostic Risk Scoring System

To further confirm the predictive effect of the risk scoring system, we used univariate analysis and multivariate analysis to analyze the impact of these three factors (age, stage, and risk score) on the prognosis of cervical carcinoma. The results showed that the stage and risk score were significantly related to the prognosis of cervical carcinoma, while age was not statistically significant (P value < 0.05, Figures 4(a) and 4(b)). These results further suggest that stage and risk score can be used as two prognostic models for cervical carcinoma.Figure 4
The accuracy of IRG signature evaluation. (a) Forest map of univariate Cox analysis in cervical carcinoma showing the prognostic values of age, stage, and risk score. The regression coefficient of age, stage, and risk score was obtained using the survival package in R.P value < 0.05 was considered significant. (b) Forest map of multivariate Cox analysis in cervical carcinoma showing the prognostic values of age, stage, and risk score. The regression coefficient of age, stage, and risk score was obtained using the survival package in R. P value < 0.05 was considered significant. (c) Time-dependent ROC curve of the survival status of cervical carcinoma for 1 year, 3 years, and 5 years. Red: combined model, yellow: age model, green: stage model, and purple: risk score model. (d) Nomogram predicting overall survival (OS) of cervical carcinoma patients for 1 year, 3 years, and 5 years. (e) 1-year, 3-year, and 5-year survival calibration curves of cervical carcinoma.
(a)(b)(c)(d)(e)Then, to observe whether the risk score system is different from the previous TNM staging system, we use age, stage, risk score, and the combination of these three factors as risk factors to build four models. ROC curve analysis was used to evaluate the accuracy of the four models in predicting 1-, 3-, and 5-year survivals of cervical carcinoma patients (Figure4(c)). The results showed that the age signature had a lower accuracy in predicting the survival rate of 1 year (AUC=0.572), 3 years (AUC=0.554), and 5 years (AUC=0.505). The stage signature had a moderate accuracy in predicting the survival rate of 1-year (AUC=0.652), but poor prediction of the 3-year (AUC=0.561) and 5-year (AUC=0.549) survival status. The risk score signature and the combined model all had a good accuracy in predicting the survival rate of 1 year (combined: AUC=0.766; risk score: AUC=0.688), 3 years (combined: AUC=0.750; risk score: AUC=0.782), and 5 years (combined: AUC=0.695; risk score: AUC=0.716). Moreover, the risk score signature had better prediction accuracy than the combined model in predicting survival rate of 3 years and 5 years. All results suggested that the survival prognostic models based on risk score are more helpful than other factors (age and stage) for the prognosis of patients with cervical carcinoma. The IRG risk score can be used as an independent prognostic factor for cervical carcinoma in the clinic.In addition, to further explore the specificity of the risk score system, we constructed a nomogram to investigate the probability of 1-year, 3-year, and 5-year OS of patients with cervical carcinoma (Figure4(d)). The calibration curve results showed that the performance of the nomogram was best at predicting 1-year, 3-year, and 5-year OS (Figure 4(e)).


### 3.4. Different Functions in Immune Cell Infiltration and Survival Analysis between the High- and Low-Risk Score Groups

Through the analysis of the TIMER database, we explored the relationship between risk score and immune cell infiltration to reflect the status of the TIME of cervical carcinoma. We found that there was a significant negative correlation between the risk score and the infiltration level of immune cells. With the increase in risk score, the infiltration level of the immune cells decreased, including B cells, CD4 cell, CD8 cell, dendritic cells, macrophages, and neutrophils (Figure5(a)). Then, we found that the risk score was closely related to the cumulative survival of these immune cells (Figure 5(b)).Figure 5
The relationship between the risk score and infiltration of six immune cells. (a) Correlation analysis between the risk score and infiltration of six immune cells (B cells, CD4 T cells, CD8 T cells, dendritic cells, macrophages, and neutrophils). (b) Kaplan-Meier curves were used to assess the survival difference in six immune cells using the log-rank test. Red: high level; blue: low level.
(a)(b)To further study the specific functions of these six IRGs, we analyzed the correlation between the expression of these six IRGs and clinicopathological factors in cervical carcinoma. The results show that most of the six IRGs were related to tumor morphology (T stage) (Figure6(a)). Moreover, only CSK was associated with metastasis (M stage) (Figure 6(c)) and Slit2 was related to pathological stage (Figure 6(d)), while the other four IRGs were not associated with lymph node (N phase) (Figure 6(b)), metastasis (phase m) (Figure 6(c)), or pathological stage (Figure 6(d)). All these results suggested that CSK and Slit2 genes are associated with the progression of cervical carcinoma and may be related to the survival rate of patients.Figure 6
Different expression of 6 IRG mRNAs in different groups. (a) Green: T3+T4 samples; red: T1+T2 samples.∗P value < 0.05; ∗∗P value < 0.01. (b) Green: N1 samples; red: N0 samples. ∗P value < 0.05; ∗∗P value < 0.01. (c) Green: M1 samples; red: M0 samples. ∗P value < 0.05; ∗∗P value < 0.01. (d) Green: stage III+IV samples; red: stage I+II samples. ∗P value < 0.05; ∗∗P value < 0.01. (e) Kaplan-Meier curves were used to assess the survival difference inexpression of 6 IRGs. Red: high expression group; blue: low expression group.
(a)(b)(c)(d)(e)Furthermore, we explored the relationship between the expression of these six IRGs and the survival probability of cervical carcinoma patients. We found that except LEPR, patients with high expression of five IRGs presented with a significantly higher survival probability than those with low expression (Figure6(e)), suggesting that these five IRGs may be tumor suppressors in cervical carcinoma.


### 3.5. Expression of CSK and Slit2 in Tumors Based on the TIMER Database

From the above results, five IRGs are tumor suppressor genes of cervical carcinoma; among them, CSK and Slit2 genes are correlated with clinical grade and lymph node metastasis of tumor. Therefore, we further verified the expression levels of CSK and Slit2 through the TIMER database. The results showed that CSK is highly expressed in most tumor tissues (Figure7(a)), while SLIT2 is highly expressed in most normal tissues (Figure 7(b)). This result is consistent with our previous analysis of cervical carcinoma expression based on TCGA database. In addition, the results of our previous analysis also showed that high expression of CSK and Slit2 presented with a good prognosis in cervical cancer (Figure 6(e)), while the expression of CSK in tumor tissues was higher than that in normal tissue. The contradictory relationship between the expression of CSK and its antitumor role suggests that CSK may play a comprehensive role in tumor.Figure 7
The expression level of CSK and slit2 in tumors based on the TIMER database. (a) The expression level of CSK in the TIMER database.∗P value < 0.05, ∗∗P value < 0.01, and ∗∗∗P value < 0.001. (b) The expression level of slit2 in the TIMER database. ∗P value < 0.05, ∗∗P value < 0.01, and ∗∗∗P value < 0.001.
(a)(b)


## 3.1. Identification of Differentially Expressed IRGs

We downloaded the expression profile data of cervical carcinoma through TCGA database, including 286 primary cervical carcinoma samples and 3 normal samples. Through DEG analysis, we identified 2503 DEGs in cervical carcinoma. And 204 DEGs are IRGs (Figure2(a)), of which 56 were upregulated and 148 were downregulated (Figure 2(b)).Figure 2
Identification of differentially expressed IRGs and gene functional enrichment. (a) Heatmap shows the expression of IRGs in cervical carcinoma tumor tissues and normal tissues. N: normal tissue; T: cervical carcinoma tissue. (b) Volcano plots demonstrating differentially expressed IRGs between cervical carcinoma and normal tissues, in the condition of afalsediscoveryrateFDR<0.05 and a log∣foldchangeFC∣>1. Red dots: upregulated genes, green dots: downregulated genes. (c) Gene Ontology (GO) analysis of differentially expressed IRGs. The top 10 most significant pathways of biological process (BP), cellular component (CC), and molecular function (MF) were shown in the training set, respectively. P value < 0.05 was considered statistically significant. (d) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of differential expressed IRGs. P value < 0.05 was considered statistically significant.
(a)(b)(c)(d)The potential molecular mechanisms of these 204 IRGs were explored by functional enrichment analyses via the GO and KEGG pathways. The GO analysis shows the top 10 clusters of enriched sets (Figure2(c)). The top 10 BPs of these IRGs were mainly including peptidyl-tyrosine phosphorylation/modification, epithelial cell proliferation, protein kinase B signaling, cell chemotaxis, and the MAPK cascade. Moreover, these IRGs were expressed in CC of platelets, extracellular matrix, vesicle lumen, granule lumen, receptor complex, side of membrane, and so on. Specifically, the MF of these IRGs is mainly related to receptor ligand/regulator activity, growth factor activity/binding, cytokine binding, transmembrane receptor protein 1 kinase activity, transforming growth factor beta binding, steroid hormone receptor activity, glycosaminoglycan binding, and nuclear receptor activity. Subsequently, the KEGG pathway results showed the top 30 significantly upregulated pathways, including cytokine-cytokine receptor interaction, neuroactive ligand-receptor interaction, Ras signaling pathway, JAK-STAT signaling pathway, MAPK signaling pathway, Rap1 signaling pathway, and PI3K-Akt signaling pathway (Figure 2(d)).


## 3.2. The Prognostic Risk Scoring System Establishment

Based on univariate Cox regression analysis, 20 IRGs with significant prognostic values were extracted with aP value < 0.05 (Figure 3(a)). Multivariate Cox proportional hazard regression analysis to further select six IRGs (CSK, SLIT2, GH1, IL3RA, LEPR, and PRLHR) to construct a risk scoring system under the condition of AIC=600.15 (Figure 3(b)). Information on these six IRGs is shown in Table 1, and the univariate/multivariate Cox regression results are shown in Figures 3(a) and 3(b). ROC curve analysis showed that the risk scoring system had a moderate accuracy in predicting the 3-year survival rate of cervical carcinoma (AUC=0.782) (Figure 3(c)). There was a significant correlation between the risk score and patient OS or mortality. With the increase in risk score, the survival prognosis worsened, and the mortality increased (Figure 3(d)). By dividing the sample into the high-risk and low-risk groups according to the median risk score in the model, we found a significant difference in survival rate (P value = 3.588e−06; Figure 3(e)) and clinical characteristics (Table 2) between the high- and low-risk groups. Further statistics showed that there were significant differences in age and fustat between the two groups (P value < 0.05), but no significant difference in stage (stages I~IV) (P value > 0.05, Figure 3(f)). These results indicate that the risk score system may be another prediction system independent of the TNM staging system.Figure 3
Survival-associated IRGs and prognosis risk scoring system construction. (a) Univariate Cox proportional hazard regression analysis in cervical carcinoma showing the prognostic values of IRGs. The regression coefficient of each gene was obtained using the survival package in R. The coefficient (parameter coef in R) of each selected gene represented the estimated logarithm of the Hazard Ratio (HR, parameter exp (coef in R)).P value < 0.05 was considered significant. (b) A risk scoring system based on 6 IRGs was built by multivariate Cox proportional hazard regression under the condition of AIC=600.15. The signature is expressed as follows: riskscore=coefficientgene1×gene1expression+coefficientgene2×expressionofgene2+⋯+coefficientgenen×expressiongenen. ∗P value < 0.05, ∗∗P value < 0.01, and ∗∗∗P value < 0.001. (c) Survival-dependent receiver operating characteristic (ROC) curve validation of the prognostic value of the risk score. (d) Rank of risk score and distribution of groups or the survival status of patients. Red: dead outcome; green: alive outcome. (e) Kaplan-Meier curves were used to assess the survival difference between two groups using the log-rank test. The median risk score served as a cutoff value to classify patients into the high-risk and low-risk groups. Red: high-risk group; purple: low-risk group. (f) Comparison of clinicopathological features and expression profiles of 6 IRGs in different groups. ∗P value < 0.05, ∗∗P value < 0.01, and ∗∗∗P value < 0.001.
(a)(b)(c)(d)(e)(f)Table 1
Information on six immune-related genes.
Gene symbolHGNC symbolFull nameCSK2444C-terminal Src kinaseSLIT211086Slit guidance ligand 2GH14261Growth hormone 1IL3RA6012Interleukin 3 receptor subunit alphaLEPR6554Leptin receptorPRLHR4464Prolactin releasing hormone receptorTable 2
The clinical characteristics of the cervical carcinoma patients in the two groups.
VariablesHigh riskLow risk(n=135)(n=136)Age (years)≤65115120>651712Unknown34Tumor stageStage I7174Stage II2633Stage III2120Stage IV145Unknown34FustatDead47117Alive8819


## 3.3. Evaluation of the Accuracy of the Prognostic Risk Scoring System

To further confirm the predictive effect of the risk scoring system, we used univariate analysis and multivariate analysis to analyze the impact of these three factors (age, stage, and risk score) on the prognosis of cervical carcinoma. The results showed that the stage and risk score were significantly related to the prognosis of cervical carcinoma, while age was not statistically significant (P value < 0.05, Figures 4(a) and 4(b)). These results further suggest that stage and risk score can be used as two prognostic models for cervical carcinoma.Figure 4
The accuracy of IRG signature evaluation. (a) Forest map of univariate Cox analysis in cervical carcinoma showing the prognostic values of age, stage, and risk score. The regression coefficient of age, stage, and risk score was obtained using the survival package in R.P value < 0.05 was considered significant. (b) Forest map of multivariate Cox analysis in cervical carcinoma showing the prognostic values of age, stage, and risk score. The regression coefficient of age, stage, and risk score was obtained using the survival package in R. P value < 0.05 was considered significant. (c) Time-dependent ROC curve of the survival status of cervical carcinoma for 1 year, 3 years, and 5 years. Red: combined model, yellow: age model, green: stage model, and purple: risk score model. (d) Nomogram predicting overall survival (OS) of cervical carcinoma patients for 1 year, 3 years, and 5 years. (e) 1-year, 3-year, and 5-year survival calibration curves of cervical carcinoma.
(a)(b)(c)(d)(e)Then, to observe whether the risk score system is different from the previous TNM staging system, we use age, stage, risk score, and the combination of these three factors as risk factors to build four models. ROC curve analysis was used to evaluate the accuracy of the four models in predicting 1-, 3-, and 5-year survivals of cervical carcinoma patients (Figure4(c)). The results showed that the age signature had a lower accuracy in predicting the survival rate of 1 year (AUC=0.572), 3 years (AUC=0.554), and 5 years (AUC=0.505). The stage signature had a moderate accuracy in predicting the survival rate of 1-year (AUC=0.652), but poor prediction of the 3-year (AUC=0.561) and 5-year (AUC=0.549) survival status. The risk score signature and the combined model all had a good accuracy in predicting the survival rate of 1 year (combined: AUC=0.766; risk score: AUC=0.688), 3 years (combined: AUC=0.750; risk score: AUC=0.782), and 5 years (combined: AUC=0.695; risk score: AUC=0.716). Moreover, the risk score signature had better prediction accuracy than the combined model in predicting survival rate of 3 years and 5 years. All results suggested that the survival prognostic models based on risk score are more helpful than other factors (age and stage) for the prognosis of patients with cervical carcinoma. The IRG risk score can be used as an independent prognostic factor for cervical carcinoma in the clinic.In addition, to further explore the specificity of the risk score system, we constructed a nomogram to investigate the probability of 1-year, 3-year, and 5-year OS of patients with cervical carcinoma (Figure4(d)). The calibration curve results showed that the performance of the nomogram was best at predicting 1-year, 3-year, and 5-year OS (Figure 4(e)).


## 3.4. Different Functions in Immune Cell Infiltration and Survival Analysis between the High- and Low-Risk Score Groups

Through the analysis of the TIMER database, we explored the relationship between risk score and immune cell infiltration to reflect the status of the TIME of cervical carcinoma. We found that there was a significant negative correlation between the risk score and the infiltration level of immune cells. With the increase in risk score, the infiltration level of the immune cells decreased, including B cells, CD4 cell, CD8 cell, dendritic cells, macrophages, and neutrophils (Figure5(a)). Then, we found that the risk score was closely related to the cumulative survival of these immune cells (Figure 5(b)).Figure 5
The relationship between the risk score and infiltration of six immune cells. (a) Correlation analysis between the risk score and infiltration of six immune cells (B cells, CD4 T cells, CD8 T cells, dendritic cells, macrophages, and neutrophils). (b) Kaplan-Meier curves were used to assess the survival difference in six immune cells using the log-rank test. Red: high level; blue: low level.
(a)(b)To further study the specific functions of these six IRGs, we analyzed the correlation between the expression of these six IRGs and clinicopathological factors in cervical carcinoma. The results show that most of the six IRGs were related to tumor morphology (T stage) (Figure6(a)). Moreover, only CSK was associated with metastasis (M stage) (Figure 6(c)) and Slit2 was related to pathological stage (Figure 6(d)), while the other four IRGs were not associated with lymph node (N phase) (Figure 6(b)), metastasis (phase m) (Figure 6(c)), or pathological stage (Figure 6(d)). All these results suggested that CSK and Slit2 genes are associated with the progression of cervical carcinoma and may be related to the survival rate of patients.Figure 6
Different expression of 6 IRG mRNAs in different groups. (a) Green: T3+T4 samples; red: T1+T2 samples.∗P value < 0.05; ∗∗P value < 0.01. (b) Green: N1 samples; red: N0 samples. ∗P value < 0.05; ∗∗P value < 0.01. (c) Green: M1 samples; red: M0 samples. ∗P value < 0.05; ∗∗P value < 0.01. (d) Green: stage III+IV samples; red: stage I+II samples. ∗P value < 0.05; ∗∗P value < 0.01. (e) Kaplan-Meier curves were used to assess the survival difference inexpression of 6 IRGs. Red: high expression group; blue: low expression group.
(a)(b)(c)(d)(e)Furthermore, we explored the relationship between the expression of these six IRGs and the survival probability of cervical carcinoma patients. We found that except LEPR, patients with high expression of five IRGs presented with a significantly higher survival probability than those with low expression (Figure6(e)), suggesting that these five IRGs may be tumor suppressors in cervical carcinoma.


## 3.5. Expression of CSK and Slit2 in Tumors Based on the TIMER Database

From the above results, five IRGs are tumor suppressor genes of cervical carcinoma; among them, CSK and Slit2 genes are correlated with clinical grade and lymph node metastasis of tumor. Therefore, we further verified the expression levels of CSK and Slit2 through the TIMER database. The results showed that CSK is highly expressed in most tumor tissues (Figure7(a)), while SLIT2 is highly expressed in most normal tissues (Figure 7(b)). This result is consistent with our previous analysis of cervical carcinoma expression based on TCGA database. In addition, the results of our previous analysis also showed that high expression of CSK and Slit2 presented with a good prognosis in cervical cancer (Figure 6(e)), while the expression of CSK in tumor tissues was higher than that in normal tissue. The contradictory relationship between the expression of CSK and its antitumor role suggests that CSK may play a comprehensive role in tumor.Figure 7
The expression level of CSK and slit2 in tumors based on the TIMER database. (a) The expression level of CSK in the TIMER database.∗P value < 0.05, ∗∗P value < 0.01, and ∗∗∗P value < 0.001. (b) The expression level of slit2 in the TIMER database. ∗P value < 0.05, ∗∗P value < 0.01, and ∗∗∗P value < 0.001.
(a)(b)


## 4. Discussion

At present, the TNM staging system is a kind of cancer grading system, which can guide clinical treatment and provide a method for predicting the prognosis of patients with cancer [17]. However, with the deepening of research, the clinical disadvantages of TNM staging become increasingly obvious [18]. In our study, we found a new prediction model of cervical carcinoma prognosis based on the expression of six IRG (CSK, Slit2, GH1, IL3RA, LEPR, and PRLHR) and it was superior to the current clinicopathological prognostic risk factors such as age and TNM stage. Moreover, an interesting result of our study is that the risk score prognosis system has better prediction accuracy than the model that combines with age, stage, and risk score in predicting the survival rates of 3 years and 5 years (Figure 4(c)). The predictive accuracy became lower when it was combined with age and stage, suggesting that there may be confounding factors.Recently, increasing evidences show that TIME is closely related to the development of cervical carcinoma and has an important value in OS or the prognosis of cervical carcinoma [5, 7]. Emerging evidence has elucidated that tumors can act as distinct immunological organs with complex TIME, and the immune cells recruited and activated in TIME are linked to the malignant phenotypes of tumors [13, 19]. Some reports showed that the cervical carcinoma progresses in TIME are mainly composed of complex biological processes between tumor infiltrating immune cells (TIICs) and stromal cells [12, 20]. TIICs are a kind of immune cell isolated from the tumor site and are considered a manifestation of the host immune response to tumors [21, 22]. In particular, TIICs are considered able to actively promote or inhibit tumor growth, which is essential for clinical outcomes of cervical carcinoma [23, 24]. Tumor-infiltrating CD8+ T cells prolong survival in patients with cervical carcinoma [25]. These reports are highly consistent with our results. In our results, there was a significant negative correlation between the prognostic risk score and the contents of six TIICs (including B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells). Moreover, the contents of these TIICs were closely related to prognosis; with the increase in TIICs, the prognosis of cervical carcinoma was better, especially with increased CD8+ TILs. These results suggest that our prognostic risk scoring system may depend on the role of TIICs in TIME.IRG is a gene spectrum that has been identified to actively participate in the process of immune activity, especially for immune cell infiltration [9, 10]. Some studies have explored that IRGs participate in the regulation of TIME and affect the survival or prognosis of cervical carcinoma [26–28]. In the present study, we systematically analyzed the expression data of IRGs in cervical carcinoma from TCGA datasets. We demonstrated that a total of 204 IRGs were differentially expressed between 286 primary cervical carcinoma samples and 3 normal samples. Moreover, 20 IRGs were crucially correlated with the survival rate of cervical carcinoma, and a prognostic risk scoring system with 6 IRGs was built to stratify the prognosis of cervical carcinoma. In addition, the prognostic risk scoring system was mainly related to the immune cell infiltration in TIME. Finally, we identified two tumor suppressor genes (CSK and Slit2) as prognostic biomarkers for cervical carcinoma.CSK is a no-receptor tyrosine kinase that is considered to be an indispensable negative regulator of Src family kinases (SFKs) [29]. It consists of three Src homologous domains: N-terminal Src homologous 3 (SH3), SH2, and C-terminal kinase [30]. Extensive evidence has shown that CSK plays an important role in the signal transduction pathways of immune receptors, including promoting the innate immune response and inhibiting T cell or mast cell activation [31, 32]. Some reports have shown that CSK acts as a tumor suppressor gene and that its suppression is an important early event in cancer pathogenesis [33, 34]. CSK SUMOylation negatively regulates its tumor suppressor function [35]. Src signaling pathways have been considered to be involved in the regulation of cervical carcinoma progression [36–39]. However, the role of CSK in cervical carcinoma has not been reported. In our study, higher expression of CSK was observed in the M0 phase than in the M1 phase (Figure 6(c)), and it was also more highly expressed in the T1/T2 phase than in the T3/T4 phase (Figure 6(a)). The diverse expression of the CSK gene in different clinical phenotypes suggests that CSK gene was related to the early topography and metastasis of tumor. We also found that the high expression of CSK was positively correlated with survival probability (Figure 6(e)), suggesting that this gene might be a tumor suppressor gene of cervical carcinoma.Slit2, a secreted glycoprotein, is described as a neuronal repellent factor that controls axonal guidance and neuronal migration [40]. In addition to its role in the nervous system, Slit2 also inhibits peripheral immune cell filtration and acts as an inhibitory factor in the development of immune responses [41]. The transmembrane Robo proteins, including Robo1 and Robo2, are well-known receptors of Slit2 [42]. Slit2 plays a comprehensive role in the progression of tumors [43, 44]. In cervical carcinoma, the expression of slit2/robo1 is abnormally low, and it plays an anticancer role [45–47]. These reports were consistent with our results. In our results, the expression of Slit2 in most tumor tissues is lower than that in corresponding normal tissues (Figure 7(b)), while the expression of the Slit2 gene in cervical carcinoma was diverse in different clinical phenotypes (Figure 6(a)) and clinical stage (Figure 6(d)). The high expression of Slit2 was related to a good prognosis (Figure 6(e)), suggesting that Slit2 acts as a tumor suppressor in cervical carcinoma.


## 5. Conclusions

In summary, we developed a prognostic risk scoring system based on the expression of IRGs in cervical carcinoma and demonstrated that (i) this system was superior to the current clinicopathological prognostic risk factors such as age and TNM stage; (ii) higher risk scores were associated with the immune infiltration levels of six immune cells, as well as poor prognosis; and (iii) the CSK and Slit2 genes are tumor suppressor genes related to the prognosis of cervical carcinoma. This study comprehensively analyzed the relationship between the prognosis of cervical carcinoma and IRGs. CSK and Slit2 are two tumor suppressor genes that can be used as prognostic markers of cervical carcinoma.



---

*Source: 9971152-2021-06-25.xml*