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

Department of Pharmacy, Shenzhen Nanshan Maternity & Child Healthcare Hospital, Shenzhen 518052, China Department of Pediatrics, Shenzhen Nanshan Maternity & Child Healthcare Hospital, Shenzhen 518052, China Department of laboratory, Shenzhen Nanshan Maternity & Child Healthcare Hospital, Shenzhen 518052, China Department of Pharmacy, The Third Xiangya Hospital, Central South University, Changsha 410013, China


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][3][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][6][7][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.

Materials and Methods
2.1. Data Acquisition and Processing. Transcriptomic RNAsequencing 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.
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.

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 rate ðFDRÞ < 0:05 and log | fold change | >1. Differentially expressed IRGs were then screen out from all differentially expressed genes based on the condition of a |Pearson correlation coefficient | >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.

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.    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-rank P value < 0.05.

Establishment and
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, 3year, 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.

Journal of Nanomaterials
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 (Figure 2(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)).

The Prognostic Risk Scoring System Establishment.
Based on univariate Cox regression analysis, 20 IRGs with significant prognostic values were extracted with a P 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 The signature is expressed as follows: risk score = ð coefficient gene 1 × gene 1 expressionÞ + ðcoefficient gene 2 × expression of gene 2Þ + ⋯ + ðcoefficient gene n × expression gene nÞ. * 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.  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 highrisk 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.

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. 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 (Figure 4(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 (Figure 4(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)).

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 ( Figure 5(a)). Then, we found that the risk score was closely related to the cumulative survival of these immune cells ( Figure 5(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) (Figure 6(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.
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 (Figure 6(e)), suggesting that these five IRGs may be tumor suppressors in cervical carcinoma.    (Figure 7(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.

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   11 Journal of Nanomaterials that IRGs participate in the regulation of TIME and affect the survival or prognosis of cervical carcinoma [26][27][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][37][38][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][46][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  Journal of Nanomaterials 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.

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.  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.