Prognostic Implications of Autophagy-Related Gene Signatures in Ewing’s Sarcoma


 Backgrand：Ewing's sarcoma (ES) is a highly aggressive malignant bone tumor with a high incidence among children and adolescents. While autophagy often plays an important role in tumor development species, we developed an autophagy-related prognostic signature based on the GEO dataset.Methods: Using the GEO database with the online website Human Autophagy Database(HADb), we screened for autophagy-related differential genes and we performed enrichment analysis and PPI analysis for autophagy-related differential genes, and then, constructed our autophagy-related prognostic gene signature by univariate Cox regression, lasso regression, and multivariate Cox regression .Results: We screened a total of 72 autophagy-related differential genes, and a total of 8 autophagy-related prognostic genes were screened by machine learning algorithms. ROC curves and survival curves also showed good predictive power of the 8 gene signatures. In addition, univariate and multivariate cox regression also showed that the 8 gene signatures were independent prognostic factors. And the results of GSEA enrichment analysis also revealed the correlation between autophagy and metabolic pathways.Conlusion：In summary, we developed an autophagy-related mRNA signature consisting of 8 mRNAs that effectively classified ES patients into low- and high-risk groups. The application of the signature in clinical treatment needs further observation to validate the validity of our findings.


Introduction
Ewing's sarcoma (ES) is a highly aggressive malignant bone tumor with a high incidence among children and adolescents [1,2]. Moreover, ES has a high metastasis rate. According to statistical analysis, at the time of diagnosis, 20%-25% of patients already have lung metastases (70-80%) and or bone metastases (40-50%) [3] . Although the recurrence rate of ES is not high, the prognosis of relapsed ES patients is very terrible. In different studies, the 5-year overall survival rate (OS) ranges from 8% to 15% [4] . The molecular mechanism of the development stage and prognosis of ES is not completely clear, therefor it is necessary to explore the prognostic factors and therapeutic targets of ES.
Autophagy is the process of transporting damaged, denatured, or aged proteins and organelles to lysosomes for digestion and degradation [5]. Autophagy also plays such a double-sided role in tumors. Under normal physiological conditions, autophagy is bene cial to maintaining the self-stable state of cells. Autophagy can also prevent the accumulation of cancer-causing damaged proteins and organelles, and inhibit the carcinogenesis of cells [6,7]. Autophagy also plays such a double-sided role in tumors.
Under normal physiological conditions, autophagy is bene cial to maintaining the self-stable state of cells. Autophagy can also prevent the accumulation of cancer-causing damaged proteins and organelles, and inhibit the carcinogenesis of cells [5,6]. Thus, in most cases, autophagy is thought to inhibit the development of early tumors and promote the development of already formed tumors [8]. This suggests that regulating autophagy may be an effective intervention for cancer treatment [9]. In the past studies of ES, it has been reported that autophagy-related genes are closely related to the development of ES. For example, TRIM3 negatively regulates autophagy by promoting the degradation of Beclin1 in ES cells, eWS-FLI1 actively regulates autophagy by increasing the expression of ATG4B in ES cells [11]. However, the relationship between autophagy-related genes and the prognosis of ES patients is still unknown. Therefore, in our study, we rst explored autophagy-related differential genes through the GEO database, then screened autophagy-related prognostic gene signatures by univariate cox regression, least absolute shrinkage and selection operator (LASSO) regression and multivariate cox regression, classi ed high-and low-risk groups, and validated our autophagy-related prognostic gene signatures by survival curves with receiver operating characteristic curve (ROC) curves. Finally, we further explored the underlying mechanisms by enrichment analysis of the high-and low-risk groups. software and "Bioconductor" packages (http://www.bioconductor.org/) were used in the data analyses. Then, we remove the batch effect through R package named "sva" to normalize the datasets that we choosed(Fig1 A,B).In addition, we used 57 ES samples from the International Cancer Genome Consortium (ICGC; https://icgc.org) ICGC for external veri cation of prognostic genes.

Screening for autophagy-related genes
We downloaded the autophagy-related gene set through the website Human Autophagy Database (HADb) (http://www.autophagy.lu/) with 232 autophagy-related genes, next we used the "limma" R package and screened the differentially expressed genes in 18 cases of normal skeletal muscle and 129 cases of ES, setting The criteria were (logFC>0.5, FDR<0.05), and nally, we selected the intersection of differential genes and autophagy-related genes by Venn analysis.

enrichment analysis and PPI network
We used the "org.Hs.eg.db" R package to enrichment analysis for autophagy-related differential genes and the "ggplot2" R package was used to plot the circles of GO analysis and the histograms of KEGG analysis. In addition, we obtained PPI reciprocal network graphs of autophagy-related differential genes from the "String" online website (https://string-db.org/) and optimized the graphs using "Cytoscape" software.

Screening for prognosis-related autophagy genes
First, a preliminary screening of autophagy-related differential genes was performed by univariate Cox regression with a screening criterion of (P<0.05). After that, the screened genes were subjected to lasso regression, which was rst proposed by Robert Tibshirani in 1996 and is known as Least absolute shrinkage and selection operator. The method is a compression estimation. It obtains a more re ned model by constructing a penalty function that makes it compress some regression coe cients, i.e., force the sum of the absolute values of the coe cients to be less than some xed value; while setting some regression coe cients to zero. Thus, it retains the advantage of subset shrinkage and is a biased estimator dealing with data with complex covariance. Finally, we performed multivariate Cox regression on the data screened by lasso regression and optimized the model again using the "step" function, and constructed the prognostic marker formula for ES based on the expression of autophagy-related genes and their correlation coe cients from the multivariate cox regression analysis as follows: risk The risk score = ∑ ∑icoe cient (genei) × expression(genei) was constructed, and then the median of the risk score was used to divide the samples into high-risk and low-risk groups.
Univariate and multivariate cox regression analyses were performed with a limited number of clinical variables based on risk score and used to con rm whether our prognostic genes were independent prognostic factors.

5, GSEA enrichment analysis
The R package "clusterpro ler" was used to perform GSEA enrichment analysis for the high and low-risk groups, and "enrichplot" was used to generate the gene set enrichment map with annotations, where pvalue cutoff = 0.05 was set as the screening criterion. The enrichment score of 90 metabolic pathways was calculated with single-sample gene set enrichment analysis (ssGSEA) in the "GSVA" R package. The set of metabolism-related genes was obtained from the online website KEGG pathway (https://www.genome.jp/kegg/).

1.Screening of autophagy-related differential genes and enrichment analysis
In the GEO cohort, we screened 129 ES and 18 normal skeletal muscle samples for differential genes and obtained a total of 5760 differential genes. 72 autophagy-related differential genes were obtained by taking intersections with 232 autophagy-related genes obtained from the HADb Human Autophagy Database (Fig2), and plotting the volcano map of autophagy-related differential genes (Fig3). We mapped the PPI protein interaction network of 72 autophagy-related genes (Fig4). GO and KEGG pathway enrichment analysis was performed, and the results of GO enrichment analysis showed that autophagy-related differential genes were mostly enriched in macroautophagy, regulation of autophagy, cellular, response to external stimulus, autophagosome assembly, autophagosome organization, and so on(Fig5). KEGG pathway enrichment analysis indicated that the genes correlated with Autophagy − animal, Autophagy − other, Kaposi sarcoma−associated, herpesvirus infection, Apoptosis, Mitophagy − animal, and so on(Fig6).

Distinction and Evaluation of 8 autophagy-related genes Prognostic Signature for ES
We selected 129 samples containing clinical prognostic information for the next analysis of 72 autophagy-related differential genes. By univariate cox regression in over sruviva l(OS), we found that a total of 22 genes were strongly associated with prognosis (P<0.05), of which 12 genes were prognostic protective factors and 10 genes were prognostic risk factors (Fig7). Then, 12 genes were screened by lasso regression analysis (Fig8 A, B).

GSEA enrichment analysis
To elucidate the biological pathways and functions associated with the risk score. We performed enrichment analysis on the high-risk group and the low-risk group. Interestingly, in the results of the enrichment analysis, there were enrichment results in the low-risk group concerning multiple metabolic pathways (Fig13), like heparin metabolic process Histidine metabolism and Propanoate metabolism. To investigate the differences in metabolic pathways between high and low-risk groups, we used the metabolic pathway-related gene set obtained from the KEGG pathway online website, quanti ed the enrichment scores of 90 metabolic pathways by the R package "GSVA", and created box plots to visualize the differences in metabolic pathways between high and low-risk groups (Fig14), We found differences in a large number of metabolic pathways in both high and low-risk groups, including Mucin type O-glycan biosynthesis, Other types of O-glycan biosynthesis, Glutathione metabolism, etc. were signi cantly enriched in the high-risk group, and then Starch and sucrose metabolism beta-Alanine metabolism Tryptophan metabolism, etc. were signi cantly enriched in the low-risk group. Among them, Histidine metabolism and Propanoate metabolism were also signi cantly different from the high-risk group (P<0.001), which was consistent with the results of our GSEA enrichment analysis.

Discussion
The in uence of autophagy on its development has been less addressed in past studies of ES, and although some autophagy-related genes have been studied, the prognostic association of autophagyrelated genes with ES patients is still unknown. In our study, we analyzed the expression of 232 autophagy-related genes in ES patients and nally identi ed 8 prognostic genes associated with autophagy, which may provide new ideas and directions for the treatment and research of ES.
Autophagy-related prognostic genes based on autophagy have been reported in colon, colorectal, plasma ovarian and non-small cell lung cancers [12][13][14][15], which illustrates the close association between autophagy and cancer development. The model we constructed included a total of 8 genes (HDAC1, EIF4EBP1, DLC1, BECN1, ATG10, SIRT1, PINK1, CFLAR). In recent years, the physiological role of HDAC1 has been explored in a variety of tumors, and silencing of HDAC1 in ovarian cancer improves chemotherapy response [16]. In glioma cells, inhibition of HDAC1 expression inhibits glioma cell invasion and induces apoptosis [17]. And in the past, HDAC1 was a commonly mutated gene in sarcoma along with TP53, NF1, and other genes [18], but its pathogenesis in ES needs to be investigated more deeply. EIF4EBP1 is a translation repressor protein encoded by a gene that competitively binds to eukaryotic translation initiation factor 4E (EIF4E), thereby inhibiting the assembly of the eIF4E complex and thus cap-dependent translation, thus initiating mRNA translation when EIF4EBP1 is downregulated, thereby promoting cell growth and proliferation [19]. A study found that DLC1 inhibits Wnt/β-catenin signaling by reducing T-cell factor 4 (TCF4) expression and the interaction between β-catenin and TCF4, and suppresses autophagy by inhibiting the dissociation of ROCK1 and Beclin1-Bcl2 complexes in hepatocellular carcinoma cell lines [20]. However, the relevance of DLC1 to autophagy in ES is unclear. becN1 / Beclin1 is a central protein that assembles cofactors to form the BECN1-PIK3C3-PIK3R4 complex, which triggers the autophagic protein cascade response [21]. The ATG family tends to regulate the production of autophagic vesicles [22], while ATG10, a member of the ATG family, is also closely associated with tumorigenesis and development, e.g., high expression of ATG10 in colorectal cancer is associated with lymphovascular in ltration and lymph node metastasis, leading to tumor migration and invasion [23]. In addition, overexpression of ATG10 in cancers such as nasopharyngeal, hepatocellular, and gastric cancers has been associated with poor patient prognosis [24][25][26]. Past studies have found that overexpression of SIRT1 is associated with ES metastasis and poor prognosis, while the SIRT1 / 2 inhibitor Tenovin-6 kills ES cells in vitro [27]. The PINK1 gene is located on the short arm of chromosome 1 and encodes a 581 amino acid serine/threonine protein kinase [28]. Current studies have shown that PINK1 plays an important role in tumorigenesis and development by inducing autophagy to eliminate dysfunctional mitochondria [29]. In contrast, CFLAR has been relatively little studied in tumors, but its downregulation is often associated with enhanced autophagy of T cells and leads to enhanced cytoprotective effects of T cells [30], which may provide a new direction for ES research.
In addition, we found multiple metabolism-related pathways appeared in the results of our GSEA enrichment analysis. By ssGSEA analysis between high and low-risk groups, it may be possible to explain some of the relationships between autophagy and metabolism, in past studies, the following hypotheses were made for the association between autophagy, tumor, and metabolism, one hypothesis is that the maintenance of metabolic function through the supply of autophagic substrates is particularly critical for tumor growth [5]. Autophagy is also essential for the removal of damaged proteins and organelles, especially mitochondria. The second hypothesis is that tumor cells require autophagy not only to provide substrates for mitochondria, but also to protect the functional pool of mitochondria by removing defective mitochondria, which is essential to maintain metabolism and survival [5]. However, there still exist metabolic pathways that are protective factors for tumors, for example, the metabolic degradation of histidine increases the sensitivity of cancer cells to methotrexate, and this approach also helps to reduce the side effects of methotrexate treatment in mouse models [31]. Therefore, different treatment regimens should be adapted for different subgroups of patients.
However, our study still has some shortcomings, on the one hand, our study is based solely on the results of machine learning and bioinformatics analysis, lacking external data validation and experimental veri cation, and on the other hand, it lacks some key clinical information because ES is relatively rare and our study is based on public databases. In addition, more in-depth studies are still needed on the developmental mechanisms of these 8 autophagy-related prognostic genes in ES.

Conlusion
In summary, we developed an autophagy-related mRNA signature consisting of 8 mRNAs that effectively classi ed ES patients into low-and high-risk groups. The application of the signature in clinical treatment needs further observation to validate the validity of our ndings.

Declarations
Funding This research did not receive any speci c grant from funding agencies in the public, commercial, or notfor-pro t sectors.
Compliance with ethical standards Con ict of interest No nancial interests are to be disclosed by the authors.

Availability of data and materials
The data used to support the ndings of this study are included within the article.

Consent for publication
All authors agree to publish the paper.  Figure 1 The points of the scatter plots visualize the samples based on the top two principal components (PC1 and PC2) of gene expression pro les without (A) and with (B) the removal of batch effect.

Figure 2
Venn analysis was used to study the Autophagy-related genes.

Figure 8
The LASSO Cox analysis identi ed 11 genes associated with prognosis and the optimal values of the penalty parameter were de ned by 1,000-round cross-validation.

Figure 9
The survival curve(A) show that low risk level has better prognosis in GEO cohort, the time-dependent ROC(B) curves of the prognostic gene signature in GEO cohort. The survival curve(C) show that low risk level has better prognosis in ICGC cohort, the time-dependent ROC(D) curves of the prognostic gene signature in ICGC cohort.