Expression proling analysis of autophagy-related genes in cervical cancer

Background: It has been demonstrated by studies globally that autophagy took part in the development of cervical cancer (CC). Few studies concentrated on the correlation between overall survival and CC patients. We retrieved signicant autophagy-related genes (ARGs) correlated to the process of cervical cancer. They may be used as prognosis marker or treatment target for clinical application. Methods: Expressions level of genes in cervical cancer and normal tissue samples were obtained from GTEx and TCGA database. Autophagy-related genes (ARGs) were retrieved accroding to the gene list from HaDB. Differentially expressed autophagy related genes (DE-ARGs) related to cervical cancer were identied by Wilcoxon signed-rank test. ClusterProler package worked in R software was used to perform GO and KEGG enrichment analyses. Univariate propotional hazard cox regression and multivariate propotional hazard cox regressions were applied to identify DE-ARGs equipped with prognostic value and other clinical independent risk factors. ROC curve was drawn for comparing the survival predict feasibility of risk score with other risk factors in CC patients. Nomogram was drawn to exhibit the prediction model constructed accroding to multivariate cox regression. Correlations between Differentially expressed autophagy related genes (DE-ARGs) and other clinical features were investigated by t test or Cruskal wallis analysis. Correlation between Immune and autophagy in cervical cancer was investigated by ssGSEA and TIMER database. Results: Fifty-six differentially expressed ARGs (DE-ARGs) were retrieved from cervical cancer tissue and normal tissue samples. GO enrichment analysis showed that these ARGs involved in autophagy, ubiquitination of protein and apoptosis. Cox regression medel showed that there were six ARGs signicantly associated with overall survival of cervical caner patients. VAMP7 (HR = 0.599, P= 0.033) and TP73 (HR = 0.671, P= 0.014) played protective roles in survival among these six genes. Stage (Stage IV vs Stage I HR = 3.985, P<0.001) and risk score (HR = 1.353, P< 0.001) were sorted as independent prognostic risk factors based on multivariate cox regression. ROC curve validated that risk score was preferable to predict survival of CC patients than other risk factors. Additionally, we found some of these six predictor ARGs were correlated signicantly in statistic with tumor grade or stage, clinical T stage, clinical N stage, pathology or risk score (all P< 0.05). The immune cells and immune functions showed a lower activity in high risk group than low risk group which is distincted by median risk score. Conclusion: Our discovery showed that autophagy genes involved in the progress of cervical cancer. Many autophagy-related genes could probably serve as prognostic biomarkers and accelerate the discovery of treatment targets for CC patients.


Introduction
One of the most challenging malignancies is cervical cancer (CC) observing among females worldwide [1]. It is showed that CC led to more than about 311 thousand people death all over the world in statistically in 2018 [2]. One of major reasons for cervical cancer is infection of high-risk human papillomavirus (HPV), although the occurrence of cervical cancer cannot be fully elucidated by HPV infection [3]. Radical hysterectomy, radical radiotherapy and chemotherapy based on cisplatin are major treatment methods for CC patients until now [4]. It is reported that patients with locally advanced cervical cancer had their ve-year overall survival (OS) increased into about 70% after chemotherapy [5].
Nonetheless, the recurrence of cervical cancer after surgery or radiotherapy remains a problem di cult to solve. [6] The circumstance of limited treatment and a poor prognosis is the reality that CC patients with relapse have to face. [7] Cells develops an important process whose name is autophagy characterized by self-degradative for sustaining trophic and energic homeostasis and excuting pathogens in cell. [8][9][10] Autophagy widely involves in a variety of cytobiology processes such as the process of cancer, material metabolism and the progress of nervous system degenerative disease. It is revealed as a survival mechanism for cells subsequently [11]. Autophagy helps cancer cells to escape from death under damage factors in tumor microenvironment such as chemotherapy drug. It is emerged that regulation of autophagy could be a promising cancer theraputic strategy [12][13][14]. Nevertheless, autophagy is observed to be a double-edged sword in many cancers including suppressing and promoting the development of cancers, which depends on the microenvironment that cancer cells live in [15]. Investigation of the potential role of autophagy in cancer origin and development showed a great signi cance for basic medicine and clinical medicine research.
Previous studies showed that the role of autophagy has been partially investigated in cervical cancer. It is reported that the expression of lnc-RNA HOTAIR accelerated autophagy and epithelial-mesenchymal transition which leads to the invasion and migration of cancer cells in radioresistant HeLa cells [16]. lnc-RNA HOTAIR increased cell proliferation and chemoresistance by means of activating WNT signaling pathway. The study form Xu et al. revealed that cisplatin chemotherapy sensitivity was improved by inhibition of autophagy in HeLa cells. [17] Jing Guo et al. performed a study which showed that Hsa_circ_0023404 sponged miR-5047 led to the enhancement of metastasis and chemoresistance through autophagy signaling in cervical cancer [18]. These discoveries demonstrated that autophagy took part in the process of cervical cancer development.
Nowadays, FIGO stage serves as a majority tool for doctor to predict the survival of cervical cancer patients in clinical. [19] There is de ciency in FIGO stage system that patients may have different individual survival time even if they are attributed to same FIGO stage [20]. In order to provide doctors with a better prognostic prediction tool for CC patients, more clinical factors should be taken into consideration. Recently, the prognosis model involved with the expression level of ARGs has become popular and been constructed in lung cancer [21], gastric cancer [22] and esophageal cancer [23], etc. So, the prognostic prediction role of ARGs in CC trigged our interest. To begin with, the differently expressed autophagy-related genes (DE-ARGs) were retrieved from gene expression pro le of tumor tissues and normal tissues. A cox regression model for predicting the survival of cervical cancer patients was constructed with ARG-related prognosis signature. The predict factors involved in this cox regression model were validated by the Kaplan-Meier analysis. The survival status discrimination e cacy of risk score was compared with other clinical factors by means of ROC curve and quanti ed by area under the curve (AUC). Moreover, GO and KEGG enrichment analysis was applied to explore the functional pathways that screened DE-ARGs involved in. Finally, we also explored the relationships between the risk scores which was counted by DE-ARGs signature and immune cells or functions.

Methods
Acquisition data from GTEx and TCGA dataset The expression level of genes in normal cervix was downloaded from GTEx database (https://www.genome.gov/Funded-Programs-Projects/Genotype-Tissue-Expression-Project), which is a website containing great quantity of gene expression data resoursed from healthy people. The expression level of genes in cervical cancer and the coresponding clinical data were downloaded from TCGA database (https://portal.gdc.cancer.gov/), which is a landmark cancer genomics program. Autophagy gene list was obtained from HADb (http://autophagy.lu/clustering/index.html), which is a public database presenting the known autophagy genes with their details. 232 ARGs was obtained to screen the gene expression pro le. Gene array data from GTEx and TCGA was normalized by means of limma package from R bioconductor software. Totally, there are 14 healthy cervix tissues and 306 cervical cancer tissues involved in this study.

Identi cation of DE-ARGs
Wilcoxon signed-rank test was applied to identify differentially expressed autophagy related genes (DE-ARGs). The cut-off values were setted based on left parameters, false discovery rate (FDR) < 0.05 and |log2 fold change (logFC)| > 1.

Construct cox regression model
The correlation between the overall survival (OS) and differently expressed genes (DE-ARGs) was rstly investivaged by univariate cox regression. Variables were screened by p-value (<0.05) presented from wald X 2 test. All the signi cant variables were applied to construct multivariate cox regression model and sorted by AIC value. Then, the risk score counted by signi cant DE-ARGs were included into multivariate cox regression with other clinical factors to construct prognosic predict model.

Function enrichment analysis
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed by clusterPro ler package [24] in R software. Results was ltered by FDR (false discovery rate) <0.05, top result were presented and recognized as signi cant items.

ROC curve analysis
The prediction value of independent risk factors was investigated by Receiver Operating Characteristic (ROC) curve with area under the curve (AUC). ROC curve was drawed by survivalROC package in R software, which is designed for survival data. AUC is applied to measure the sensitivity and speci city of prediction variable, which ranges from 0.5 to 1. The larger the AUC is, the better variable predict the prognosis.

Development of the nomogram
Nomogram was drawn to exhibit the prediction model constructed based on independent prognostic factors sorted by univariate and multivariate cox regression. In order to evaluate the calibration and discrimination of nomograph, calibration curves were plotted and Harrell C-index was calculated. Bootstrapping validation with 100 bootstraps resample was applied to calculate C-index for this nomogram.

Immune and autophagy in cervical cancer
The in ltrating scores of 16 immune cells was measured by the single-sample gene set enrichment analysis (ssGSEA) in the "gsva" R package. The activities of 13 immune-related pathways was evaluated in the same way [25]. The annotated gene set le applied in ssGSEA analysis is provided in Additional le 1: Table 1 [26]. Futhermore, TIMER database (https ://cistrome.shinyapps.io/timer) was searched to explore the correlation between prognostic genes sorted from ARGs and seven kinds of immune cells.

Results
Sort DE-ARGs from gene expression pro le An expression pro le dataset was combined with data from GTEx and TCGA database, which included 14 normal cervix samples and 306 cervical cancer samples. The clinical data of the patients was downloaded from TCGA database and intergrated into expression matrix by perl software. We obtained 56 differently expressed autophagy-related genes (DE-ARGs) by comparison of the expression level of autophagy-related genes (ARGs) between cervix tissues and cervical cancer tissues. There were 29 genes down-regulated in tumor samples and 27 genes up-regulated in tumor samples compared with normal samples ( Table 1). The detail of ARGs' expression matrix was presented by heatmap, volcano plot and bar plot ( Figure.1a-c).

Bio-functional enrichment analysis of DE-ARGs
We conducted GO and KEGG enrichment analysis by clusterpro ler package in R software to evaluate the biological function of our retrieved DE-ARGs. In GO analysis ( Figure.2a-b), biological process (BP) category mainly include autophagy, process utilizing autophagic mechanism, and macroautophagy. They have something to do with the process of cervical cancer. Items in cellular component (CC) mainly include autophagosome and late endosome. They involved in autophagy and celluar internal environment homeostasis maintainment, which may help cervical cancer cells survive. Molecular function (MF) category showed that DE-ARGs were mainly enriched in ubiquitin−like protein ligase binding, ubiquitin protein ligase binding and phosphatase binding, respectively. It reveals that autophagy genes may regulate tumor growth by taking part in downstream protein phosphorylation or ubiquitination. In result of KEGG analysis ( Figure.2c-d), DE-ARGs were enriched in pathways such as Autophagy−animal, Apoptosis and Platinum drug resistance. These pathways help cancer cells protect themselves from poison of chemotherapeutic drugs.
Heatmaps for GO and KEGG revealed the details how DE-ARGs contribute to the enrichment of each items. ( Figure.2e-f).
Retrieve DE-ARGs related to prognosis.
Moreover, these genes were applied to construct multivariate cox regression model, AIC value was used to sort the variable. Six ARGs (HGS, VAMP7, TM9SF1, ERO1A, SUPT20H, TP73) were preserved at last and regarded as prognosic related DE-ARGs (Table2). VAMP7 and TP73 played protective roles (HR<1) among these six ARGs. The other four genes (HGS, TM9SF1, ERO1A, SUPT20H) were presented as risk factors for CC patients'survival (HR>1). Finally, the risk score was calculated according to the expression level of these six genes ( , n=6, Coef j is the coe cient of each ARG, X j is the relative expression levels of each ARG). The median risk score value was regarded as cutoff point to devide the CC patients into high risk group (n=152) and low risk group (n=152). The patients'overall survival (OS) of high risk group is shorter than that of low risk group signi cantly (median time = 3.81 years vs. >8 years, log rank p<0.001, Figure.3c).

Draw prognostic hazard curves
Prognostic hazard curves was drawn to evaluate the survival time for the patients. It is observed that the survival time diminished with the increasing of risk score for the dead patients ( Figure.4a-b). Furthermore, the quantity of patients alive decreased with the ascend of risk score for patients too. TP73 was shown to be down-regulated in group with high-risk accroding to the risk heatmap. Whereas, TM9SF1 was regarded as a tumor accelerating gene because it was up-regulated in high-risk group. (Figure. For the sake of evaluating the discrimination of each predicting factors, ROC curves were constructed in 0.5-year, 1-year, 3-year and 5-year with the prediction factors (age, stage and risk score). Moreover, we assess the feasibility of discrimination of survival or dead patients using the area under curve (AUC) values. ROC curve reveals that the risk score showed a better ability to predict the survival of CC patients (AUC=0.862, 0.761, 0.777, 0.810 for 0.5-year, 1-year, 3-year and 5-year) than other prediction factors ( Figure.

Discussion
There are two major metabolic processes for eukaryotic cells to maintain cellular hemostasis. One of them is ubiquitin-proteasome system, the other one is autophagy. They help cells to sweep away damaged proteins or broken organelles [27]. Great quantities of studies revealed that autophagic dysfunctions have something to do with many kinds of diseases process such as the immune of infections, the development of neurodegeneration and tumor origination [28][29][30]. Autophagy plays a twosides role in the process of cancer. It may inhibit the origination while promote the growth of cancer by a context-dependent way. Both of clinical features and tumor microenvironments may affect the role of autophagy [31]. Our study retrieved 56 DE-ARGs between normal tissues and cervical cancer tissues. 29 ARGs showed down-regulated and 27 ARGs showed up-regulated in tumor tissues compared with normal tissues. In result of GO and KEGG pathways analysis. It is discovered that these DE-ARGs were mainly enriched in autophagy and the cellular apoptotic signaling pathway. We sorted 6 prognostic ARGs by wald X 2 test and AIC value in cox regression model. Among these ARGs, VAMP7 and TP73 with HR<1 played a protective role in survival. Other four genes (HGS, TM9SF1, ERO1A, SUPT20H) were considered as risk factors with HR>1. The risk score of ARGs and FIGO stage were conbined to construct a new predict model. It may be valuable to evaluate the prognosis of CC patients.
The correlation between autophagic gene and cervical cancer has been validated by great quantity of studies [17,[32][33][34]. CDKN2A (aliase p16) is a famous autophagic gene which regulates the cellular cycle process of G1/S. It is revealed as a tumor suppressor gene by a former study [35]. It is able to inhibit proliferation and invasion of CC cells through LDHA-mediated AKT/mTOR pathway [35]. The inactivation of CDKN2A led by promoter methylation may contribute to the malignant phenotype in CC cells [36]. ERO1A is one of ARGs sorted by multivariate cox regression which was associated with survival. It is an autophagy related gene playing crucial role in many cancers [37]. Wendi Yan et al. revealed that the expression level of ERO1A is correlated with clinical stage and pathological stage of cholangiocarcinoma. The deletion of ERO1A can inhibit the proliferation and migration of CCA cells and vice versa by EMT and Akt/mTOR pathways [38]. Junfeng Zhang et al. discovered the mechanism that the growth and glycolysis of pancreatic ductal adenocarcinoma is promoted by ERO1A (also named as ERO1L) [39]. Kyoung Min Kim et al. demonstrated that combined expression of protein disul de isomerase and ERO1A can be applied as a poor prognostic marker for non-small cell lung cancer [40]. It is reported by Fei Han et al. that ERO1L accelerated the progression of pancreatic cancer cells by activating the Wnt/catenin pathway [41]. Yini Zhang et al. found the evidence that the progression of cervical cancer can be suppressed by targeting the functional interplay between ERO1A and PDI [42]. Our results demonstrated that expression level of ERO1A is higher in cancer group than normal group. It is a poor prognostic marker demonstrated by multivariate cox regression.
One of other mechanismes for cell death is cell apoptosis, which often works together with autophagy. Both of them involved in great kinds of catalytic process of enzymes. The sophisticated interaction details between autophagy and apoptosis remains to be revealed [43]. It is revealed in our GO and KEGG enrichment analyses that DE-ARGs retrieved from the gene expression pro le of CC tissues and normal cervix tissues were mainly enriched in the cellular apoptotic pathways, which is consistent with former studies [44,45]. BAX and BAK1 are two members of BCL2 apoptosis family. They worked not only as messages of apoptotic signal converge, but also as autophagy regulator [45]. In the cervical cancer microenvironment, many kinds of different moleculars mediate the mutual impact of autophagy and apoptosis. Liya Ma et al. performed a study demonstrated that the silence of MLK3 induces the apoptosis of CC cells via the Notch-1/autophagy network [46]. An experiment performed by Liu S et al.
revealed that miRNA-211 regulates the Bcl-2 triggering an autophagy-dependent apoptosis in CC cells [47]. A study performed by Zhen Yang et al. discovered that the upregulating of TMED5 and LMNB1 promoted malignant behavior and nuclear autophagy in CC cells by a GRSF1-mediated way [48]. The long noncoding RNA LINC00511 can be inhibited to promote celluar autophagy and apoptosis via signal pathways including transcription factor RXRA-regulated PLD1 in CC cells [49]. Dewen Yan et al. discovered that chemotherapeutic may induce SHP-2 restricting apoptosis by Parkin-dependent autophagy pathway in CC cells [50]. Their discovery revealed that SHP-2 might activate autophagy and induce the degradation of damaged mitochondria and ubiquitin ligase Parkin. The relation between apoptosis and autophagy has been described by many study while its potential mechanism still ramains to be validated in CC cells.
The risk score calculated by expression level of DE-ARGs was demonstrated to be a risk indicator. Patients in high risk group shows a signi cant lower survival than low risk group. ROC curve for risk factors suggested that risk score predicted the prognosis better than other factors which may be valuable in clinical application. It suggested that more precise theraputic strategy should by applied to CC patients with high risk score. At last the expression of each ARGs were analysed with patients'clinical features. TM9SF1 was signi cantly associated with clinic stage. The expression level of TM9SF1 and HGS is increased with the ascend of T stage. VAMP7 expression level go down with the ascend of N stage. The expression VAMP7 may prevent the CC from lymph node metastasis. Both of ERO1A, SUPT20H, TM9SF1, TP73 and risk score is correlated with pathology. The expression level of ERO1A decreased with the progression of cancer grade. This information may be a clue for further research about autophagy related genes and clinical feature in cervical cancer Thanks to public database such as TCGA and GTEx, the correlation between prognosis of CC patients and autophagic genes was analysed. A robust statistical support could be used to help the autophagy genes researchers in future. More clinical tharaputic shemes should be developed concentrating on autophagy genes in CC patients. The up-regulation or down-regulation of autophagic genes have been measured in this study. Whereas, the autophagy level of cells overall could not be measured by expression level of ARGs. There are some limits in this study consequently. To begin with, the clinical stage, pathology grade and the treatment shemes downloaded from TCGA were incomplete. The HPV infection status of each patient was unknown. These de ciencies affected the acuracy of prediction model we constructed at last. Moreover, the potential mechanisms of how ARGs regulate the development of CC and their interaction relationship were remained to be explained. Finally, the nomograph has to be validated in a larger cohort, that will be helpful for further epidemical research. These de ciencies could be solved with a larger scale of clinical data appeared in future.

Conclusions
DE-ARGs with prognostic value in CC patients were identi ed by wilcoxson rank-sum test, cox regression, wald X 2 , and AIC value. They have potential to be diagnose or treatment biomarks for CC patients. These DE-ARGs and nomograph should be validated in future to com rm our study's result. There is no need for ethics committee approval due to the reason that the data were downloaded from public databases.

Consent for publication
This manuscript has been consent for publication by all the authors.

Competing interests
The authors declare that they have no con icts of interest with the contents of this article.          . calibration curves of the prognostic cytokine-cytokine receptor, APC antigen presenting cells. Adjusted P values were shown as: ns, not signi cant; *P < 0.05; **P < 0.01; ***P < 0.001.