Slow skeletal muscle troponin T, titin and myosin light chain 3 are candidate prognostic biomarkers for Ewing's sarcoma.

Ewing's sarcoma (ES) is a common malignant bone tumor in children and adolescents. Although great efforts have been made to understand the pathogenesis and development of ES, the underlying molecular mechanism remains unclear. The present study aimed to identify new key genes as potential biomarkers for the diagnosis, targeted therapy or prognosis of ES. mRNA expression profile chip data sets GSE17674, GSE17679 and GSE45544 were downloaded from the Gene Expression Omnibus database. Differentially expressed genes (DEGs) were screened using the R software limma package, and functional and pathway enrichment analyses were performed using the enrichplot package and GSEA software. The NetworkAnalyst online tool, as well as Cytoscape and its plug-ins cytoHubba and NetworkAnalyzer, were used to construct a protein-protein interaction network (PPI) and conduct module analysis to screen key (hub) genes. LABSO COX regression and overall survival (OS) analysis of the Hub genes were performed. A total of 211 DEGs were obtained by integrating and analyzing the three data sets. The functions and pathways of the DEGs were mainly associated with the regulation of small-molecule metabolic processes, cofactor-binding, amino acid, proteasome and ribosome biosynthesis in eukaryotes, as well as the Rac1, cell cycle and P53 signaling pathways. A total of one important module and 20 hub genes were screened from the PPI network using the Maximum Correlation Criteria algorithm of cytoHubba. LASSO COX regression results revealed that titin (TTN), fast skeletal muscle troponin T, skeletal muscle actin α-actin, nebulin, troponin C type 2 (fast), myosin light-chain 3 (MYL3), slow skeletal muscle troponin T (TNNT1), myosin-binding protein C1 slow-type, tropomyosin 3 and myosin heavy-chain 7 were associated with prognosis in patients with ES. The Kaplan-Meier curves demonstrated that high mRNA expression levels of TNNT1 (P<0.001), TTN (P=0.049), titin-cap (P=0.04), tropomodulin 1 (P=0.011), troponin I2 fast skeletal type (P=0.021) and MYL3 (P=0.017) were associated with poor OS in patients with ES. In conclusion, the DEGs identified in the present study may be key genes in the pathogenesis of ES, three of which, namely TNNT1, TTN and MYL3, may be potential prognostic biomarkers for ES. Copyright: © Deng et al.

Introduction

Ewing's sarcoma (ES) is the second most common primary malignant bone tumor after osteosarcoma, and it occurs in children and adolescents (1,2). ES is extremely malignant, with a short course of disease, rapid recurrence and high transfer rate (3). With the continuous improvements in ES treatment, the current 5-year survival rate is between 65 and 75% (4). However, owing to the lack of effective diagnostic methods in the early stages of the disease, ~25% of patients with ES still experience distant metastasis, which results in poor prognosis (5,6). Therefore, clarifying the precise molecular mechanisms involved in the development of ES is essential to develop effective diagnostic and therapeutic strategies.

Abnormal expression and mutations of genes are involved in the development and progression of ES. A previous study has reported that nuclear phosphoprotein (NPM) promotes proliferation and invasion of ES cells, and elevated NPM expression may be associated with poor prognosis of ES (7). Insulin-like growth factor 1 (IGF1) and its receptor (IGF1-R) serve a key role in the progression of ES by interfering with the IGF1R pathway in ES cells and subsequently inhibiting cell proliferation, promoting apoptosis and reducing invasion and metastasis (8,9). KIT proto-oncogene receptor tyrosine kinase is a tyrosine kinase receptor that is significantly expressed in ES associated with the proliferative, invasive and metastatic ability of ES cells (10). Stromal antigen 2 mutation occurs in 20% of ES cases and is associated with distant metastasis, although whether it may be considered a prognostic marker for ES remains controversial (11). Other molecular genetic alterations associated with ES include abnormal expression of platelet-derived growth factor receptor β (12,13) and mammalian target of rapamycin (14), as well as CDKN2A and TP53 mutations (15,16). A better understanding of the molecular biology of ES may help identify novel early diagnostic biomarkers, potential therapeutic targets or prognostic indicators.

In recent decades, gene chip technology and bioinformatics have been widely used to screen genetic changes at the genomic level, which can identify the differentially expressed genes (DEGs) and functional pathways involved in ES carcinogenesis and progression. In this study, three mRNA expression profile chip data sets from Gene Expression Omnibus (GEO) (17) were downloaded and analyzed to obtain DEGs between ES and normal tissues. Functional and pathway enrichment analysis and protein-protein interaction (PPI) network analysis of DEGs were also performed. Finally, LASSO COX regression model and overall survival rate (OS) analysis of hub genes were performed. The results provided a useful framework for elucidating the pathogenesis of ES and its complex molecular biology and for identifying new key genes that may be used as prognostic biomarkers for ES.

Materials and methods

Microarray data

A total of three gene-expression data sets were obtained from the GEO database (http://www.ncbi.nlm.nih.gov/geo) (17): GSE17674 (18), GSE17679 (18) and GSE45544 (19). The GSE17674 data set was based on the GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array) and included 44 samples from patients with ES and 18 adjacent non-cancerous tissues from the patients. The GSE17679 data set was based on the GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array) and included 88 patient samples and 18 adjacent non-cancerous tissues from the patients. The GSE45544 data set was based on the GPL6244 platform (Affymetrix Human Gene 1.0 ST Array) and included 14 patient samples and 22 normal tissues from healthy control subjects.

Data preprocessing and DEG screening

Raw data from the three data sets (CEL file) were read using the R language (version 3.1.2; http://r-project.org/) Affy package (20). The original CEL file was removed and the data were subjected to background correction, bootstrap correction, quality control and normalization processing. The R package sva (21) was used to perform batch effect removal on the three data sets. The data were converted into a probe expression matrix and analyzed by the R language limma package (22) to obtain DEGs. The DEGs were filtered by an adjusted P-value (adj.P-value) <0.05 and |log2 fold change (FC)|>1.

DEG function and pathway enrichment analysis

To analyze the functions of the DEGs, the R software enrichplot package (23) was used to perform GO and KEGG pathway enrichment analysis. P<0.05 was considered to indicate significant gene enrichment.

Gene Set Enrichment Analysis (GSEA) was performed using c2.cp.kegg.v6.0.symbols.gmt (http://software.broadinstitute.org/gsea/msigdb/download_file.jsp?filePath=/resources/msigdb/7.0/c2.cp.kegg.v7.0.symbols.gmt) (24) as a reference gene set. The GSEA software (version 6.3) is available on the GSEA website (http://software.broadinstitute.org/gsea/index.jsp). False Discovery Rate (FDR) <0.25 and P-value <0.05 were used as the cut-off criteria.

PPI network construction and hub gene selection

The NetworkAnalyst (http://www.networkanalyst.ca) online tool was used to build a PPI network of the DEGs with parameters set to a confidence score ≥900. The Cytoscape (version 3.6.1) software (25) was used to visualize the PPI network and identify the most important modules in the PPI network through the NetworkAnalyzer plugin (26). To extract valuable information from the important modules, the cytoHubba plugin (27) was used to identify the hub genes. The top 20 genes in Maximum Correlation Criteria (MCC) were selected by the cytoHubba plugin and sorted by enrichment fraction.

LASSO COX regression and overall survival (OS) analysis

The LASSO COX regression model and Kaplan-Meier OS analysis were used to screen genes with strong association with prognosis from the hub genes and assess their effects on survival in ES. LASSO coefficient profiles of the 10 ES-associated genes were analyzed, and a vertical line was drawn at the value selected by 10-fold cross-validation. OS analysis was performed by R-survival (28) and survminer (https://cran.r-project.org/web/packages/survminer/index.html) based on high- and low-expression levels of gene expression determined by BestSeparation (29), and data from genes with P<0.05 were retained for display.

Results

Identification of DEGs in ES

The R package sva (21) was used to perform batch effect removal on the three data sets. The principal component analysis (PCA) before and after batch effect removal is presented in Fig. 1A and B, respectively. The results demonstrated that the sample data set selected in the present study was of reliable quality. Following standardization of the chip results, 1,133, 1,290 and 768 DEGs between ES and healthy tissues were extracted from the GSE17674, GSE17679 and GSE45544 mRNA expression profile data sets, respectively. The overlap between the three data sets contained 211 DEGs, as presented in the Venn diagram (Fig. 1C).

Figure 1.

Sample pretreatment and Venn diagram. (A and B) PCA (A) before and (B) after batch effect removal. (C) Screening of DEGs. A total of 211 overlapping DEGs were identified in the three data sets. The DEGs satisfied the criteria adjusted P-value <0.05 and |log2 fold change|>1. PCA, principal component analysis; Comp, component; ES, Ewing's sarcoma.

Function and pathway enrichment analysis of DEGs

To determine the biological functions of the DEGs, GO and KEGG pathway enrichment analyses were performed using the R software enrichplot package (23). The results of GO analysis demonstrated that the BP changes of the DEGs were significantly enriched in ‘response to peptides’, ‘regulation of body fluid levels’, ‘response to extracellular stimulus’, ‘regulation of small molecule metabolic process’ and ‘gland development’ (Fig. 2A). Changes in MF were mainly enriched in ‘cell adhesion molecule binding’, ‘cofactor binding’, ‘proximal promoter sequence-specific DNA binding’, ‘RNA polymerase II proximal promoter sequence-specific DNA binding’ and ‘transcriptional activator activity-specific DNA-binding (Fig. 2A). The CC of DEGs varied mainly in ‘neuronal cell body’, ‘cell-substrate junction’, ‘cell-substrate adherens junction’, ‘focal adhesion’ and ‘vesicle lumen’ (Fig. 2A). The results of KEGG pathway enrichment analysis demonstrated that the DEGs were mainly enriched in ‘amino acid biosynthesis’, ‘DNA replication’, ‘focal adhesion’, ‘proteasome’ and ‘ribosomal biogenesis in eukaryotes’ (Fig. 2B). GSEA results revealed that the enriched functions and pathways mainly involved the p53 (P=0.004; NES=1.528), Rac1 (P=0.010; NES=1.544) and cell cycle (P=0.010; NES=1.542) pathways (Fig. 3).

Figure 2.

GO and KEGG pathway enrichment analysis. (A) GO analysis. The color of the dot represents the adjusted P-value: Red, low; blue, high. The size of the dots represents the number of DEGs. (B) KEGG pathway enrichment analysis. The color of the dots represents the adjusted P-value; the size of the dots represents the number of DEGs in the pathway. DEGs were mainly enriched in ‘biosynthesis of amino acids’, ‘DNA replication’, ‘focal adhesion’, ‘proteasome’ and ‘ribosome biogenesis in eukaryotes’. GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; p.adjust, adjusted P-value; BP, biological process; CC, cellular component; MF, molecular function.

Figure 3.

Gene Set Enrichment Analysis results. The red line represents the cell cycle pathway, the cyan line represents the P53 signaling pathway, and the blue line represents the Rac1 pathway. NES, Normalized Enrichment Score; KEGG, Kyoto Encyclopedia of Genes and Genomes; ES, Ewing's sarcoma; PID, the Pathway Interaction Database.

PPI network construction and selection of hub genes

The NetworkAnalyst online tool was used to build a PPI network of DEGs to analyze PPIs (Fig. 4A). The module with the highest MCC score in this PPI network was identified using the NetworkAnalyzer plugin in the Cytoscape software (Fig. 4B). In addition, the cytoHubba plugin was used to select the top 20 genes in MCC as the hub genes: Titin (TTN), myosin light-chain 3 (MYL3), fast skeletal muscle troponin T, troponin C type 2 (fast), tropomyosin 3, MYL2, troponin C type 1 (slow), troponin I2 fast skeletal type (TNNI2), troponin I1 slow skeletal type, myosin-binding protein C fast type, myosin-binding protein C slow type, slow skeletal muscle troponin T (TNNT1), titin-cap (TCAP), nebulin, desmin, tropomodulin 1 (TMOD1), myosin heavy-chain 7 (MYH7), MYH1, MYH2 and skeletal muscle α-actin (Fig. 4C). The functions of these hub genes are presented in Table I.

Figure 4.

PPI network and hub genes. (A) A PPI network was constructed using the NetworkAnalyst online tool. (B) The NetworkAnalyzer plugin identified the most important modules in the PPI network. The node size represents the clustering coefficient and the proportion of genes in the network. The color of the node represents the degree (blue, high; yellow, medium; orange, low). The thickness of the connection represents the comprehensive score. (C) The cytoHubba plugin selected the top 20 genes in the Maximum Correlation Criteria as hub genes. PPI, protein-protein interaction.

Table I.

Functional roles of top 20 hub genes.

No.	Gene symbol	Gene name	Function
1	TTN	Titin	Mutation may serve a specific role in the development or progression of colorectal cancer
2	MYL3	Myosin light-chain 3	Mutations have been identified as a cause of mid-left ventricular chamber type hypertrophic cardiomyopathy
3	TNNT3	Fast skeletal muscle troponin T	Promotes the progression of breast cancer
4	TNNC2	Troponin C type 2 (fast)	Serves a key role in the regulation of muscle contraction and modulates the Ca2+-activation characteristics of muscle fibers
5	TPM3	Tropomyosin 3	Potential biomarker for colorectal cancer
6	MYL2	Myosin light-chain 2	Associated with invasion, metastasis and poor prognosis of several cancers
7	TNNC1	Troponin C type 1 (slow)	Potential marker for predicting occult cervical lymphatic metastasis and prognosis of oral tongue carcinoma
8	TNNI2	Troponin I2, fast skeletal type	High expression in gastric tissue is a specific biomarker for peritoneal metastasis of gastric cancer
9	TNNI1	Troponin I1, slow skeletal type	Downregulation restrains proliferation of non-small-cell lung carcinoma xenografts
10	MYBPC2	Myosin-binding protein C, fast type	May modulate muscle contraction or serve a structural role
11	MYBPC1	Myosin-binding protein C,	Missense and nonsense mutations have been directly linked
		slow type	with the development of severe and lethal forms of distal arthrogryposis myopathy and muscle tremors
12	TNNT1	Slow skeletal muscle troponin T	High expression is associated with cell proliferation and migration in a variety of cancers
13	TCAP	Titin-cap	Important for accurate diagnosis and prognosis of breast cancer subtypes
14	NEB	Nebulin	Mutation is involved in the development of osteosarcoma
15	DES	Desmin	Potential oncofetal diagnostic and prognostic biomarker in colorectal cancer
16	TMOD1	Tropomodulin 1	Involved in the development, invasion and metastasis of various tumors
17	MYH7	Myosin heavy-chain 7	Forms the basic contractile unit of skeletal muscle and myocardium
18	MYH1	Myosin heavy-chain 1	Involved in the development, invasion and metastasis of breast cancer
19	MYH2	Myosin heavy-chain 2	Potential driver of squamous cell lung cancer
20	ACTA1	Skeletal muscle α-actin gene	Involved in various types of cell movement and is commonly expressed in all eukaryotic cells

LASSO COX regression and OS analysis

LASSO COX regression and OS analysis of the 20 hub genes was performed. LASSO COX regression results demonstrated that 10 genes were associated with prognosis (Fig. 5), whereas Kaplan-Meier curves revealed that high mRNA expression levels of TNNT1 (HR, 3.57; 95% CI, 2.06–6.18; P<0.001), TTN (HR, 1.82; 95% CI, 1.04–3.16; P=0.049), TCAP (HR, 1.74; 95% CI, 0.97–3.14; P=0.04), TMOD1 (HR, 2.03; 95% CI, 1.18–3.50; P=0.011), TNNI2 (HR, 1.94; 95% CI, 1.13–3.34; P=0.021) and MYL3 (HR, 2.02; 95% CI, 1.17–3.49; P=0.017) were associated with poor OS in patients with ES (Fig. 6). Based on these results, TNNT1, TTN and MYL3 may be key players in the abnormal signaling pathway of ES and may serve as potential prognostic biomarkers for ES.

Figure 5.

LASSO coefficient profiles of the 20 hub genes. A vertical line is drawn at the value chosen through 10-fold cross-validation. Red text indicates genes associated with prognosis obtained by LASSO COX regression.

Figure 6.

Kaplan-Meier curves. (A-F) OS curves of (A) TTN, (B) TCAP, (C) TMOD1, (D) TNNI2, (E) MYL3 and (F) TNNT1. OS analysis was based on high and low expression of the genes using Best Separation. P<0.05 was considered to indicate a statistically significant difference. OS, overall survival; TTN, titin; TCAP, titin-cap; TMOD1, tropomodulin 1; TTNI2, fast skeletal type troponin I2; MYL3, myosin light-chain 3; TNNT1, slow skeletal muscle troponin T.

Discussion

ES is a cancer of the bones and muscles, but it also occurs in soft tissues (30,31). ES is the second most common primary malignant bone tumor in children and adolescents, with 2.9 people per million diagnosed with ES annually worldwide (32). The metastasis rate is ES was 20–25%, and the recurrence rate is 30–40% (33,34). Although multiple molecular studies on the pathogenesis of ES have been conducted (35–37), the underlying molecular mechanism is largely unknown. In addition, the prognosis of patients with ES remains uncertain. Therefore, there is an urgent need to identify novel potential biomarkers for early diagnosis, targeted therapy or prognostic evaluation of ES to improve the prognosis of patients with ES. Microarray analysis is a high-throughput technology that simultaneously detects the expression levels of thousands of genes. Microarray technology was used in the present study to explore genetic changes in ES, as numerous studies have demonstrated it to be an effective way to identify new biomarkers for other diseases (38,39). The present study aimed to enhance ES data through bioinformatics methods to improve the understanding of this disease.

In the present study, three mRNA microarray data sets were downloaded from the GEO database, and an in-depth analysis was conducted using bioinformatics methods to obtain DEGs between ES and normal tissues. Prior to the analysis of the data sets, the inter-assay difference removal, background correction, bootstrap correction, quality control and standardization was performed to ensure that the data were valid for the next step of analysis. A total of 211 DEGs were identified in all three datasets between ES and normal tissues.

Since the general difference analysis (GO and KEGG) focuses on comparing gene expression differences between two groups and on genes that are significantly upregulated or downregulated, genes that are not significantly differentially expressed but serve important biological significance are easily omitted. This may lead to neglecting valuable information on the biological characteristics of certain genes, associations between gene regulatory networks, functions and significance of genes. The GSEA software does not need to specify a clear differential gene threshold; the algorithm performs enrichment analysis on all genes in the expression profile based on the overall trend of the actual data; mathematical statistics link the expression spectrum chip data with biological meaning to avoid missing important information (24). Therefore, GSEA, GO, and KEGG pathway enrichment analyses were performed in the present study; the results demonstrated that the DEGs were mainly involved in functions and pathways associated with ES development and progression, including the ‘Rac1 pathway’, ‘cell cycle pathway’, ‘RNA metabolism’ and ‘P53 signaling pathway’. The Rac1 pathway is closely associated with the invasion and metastasis of ES; a previous study has demonstrated that Erb-B2 receptor tyrosine kinase 4 mediates Rac1 GTPase activation and enhances ES invasion and metastasis in vivo by activating the PI3K-Akt cascade and focal adhesion kinase (40). The cell cycle pathway serves a key role in the malignant proliferation process in breast (41), endometrial (42), bladder (43) and gastric (44) cancer. A recent study reported that the cell cycle pathway is also involved in the ES proliferation process (45). In addition, the P53 and RNA metabolism signaling pathways serve an important role in the progression of ES (46,47). Collectively, these results are consistent with the present study.

A number of ES biomarkers in the literature match the DEGs identified in the present study. For example, Tong et al (48) analyzed and identified Fc fragment of immunoglobulin G receptor and transporter, olfactomedin 1, Cbp/P300-interacting transactivator with glu/asp, caveolin 1,3-ketodihydrosphingosine reductase, cadherin 3 (CHD3), growth arrest-specific 1, four and a half LIM domains 3 and thyroid hormone receptor interactor 6 as biomarkers of ES. Cheung et al (49) identified six-transmembrane epithelial antigen of prostate1, NK2 homeobox 2 and cyclin D1 as biomarkers of ES based on the gene expression array approach. Other studies have analyzed global genomic and transcriptomic expression to identify prognostic biomarkers for ES, and the results demonstrated that CDH11 (50) and nucleophosmin 1 (51) were biomarkers of ES. All of the above genes were differentially regulated in the present study.

A PPI network of DEGs was constructed in the present study to detect interactions between the proteins coded by the DEGs, and one important module was extracted from the network. Subsequently, the top 20 genes in MCC were selected from the important module as the hub genes. To assess the effect of these 20 hub genes on survival in ES, survival analysis was performed using the LASSO COX regression and Kaplan-Meier curves. The results revealed that high mRNA expression levels of TNNT1, TTN and MYL3 were significantly associated with poor OS in patients with ES, suggesting that these genes may serve an important role in the development of ES.

The protein encoded by TNNT1 is a subunit of troponin, a regulatory complex located on the sarcomeric filament, which regulates the contraction of striated muscle in response to fluctuations in intracellular calcium concentration (52). Microinjections of tropomyosin into epithelial cells can induce rapid cell migration (53); TNNT1 may contribute to the interaction between actin and tropomyosin, thus regulating cell migration and invasion (54). Rhabdomyosarcoma is a cancer associated with connective tissue, and the cause of this sarcoma is unclear (55). However, its malignant biological behavior may be associated with the expression of TNNT1, and TNNT1 can be used as a biomarker for the diagnosis and prognosis of rhabdomyosarcoma (48). In addition, the expression of TNNT1 is significantly increased in uterine sarcoma and closely associated with cell proliferation and migration (56,57). Highly expressed TNNT1 has been demonstrated to be a prognostic biomarker for the development and poor prognosis of numerous types of cancer, such as breast and endometrial cancers. Shi et al (58) reported that the expression of TNNT1 was significantly increased in breast cancer and was closely associated with clinical stage, tumor size and lymph node involvement. Further experiments demonstrated that TNNT1 promoted the proliferation of breast cancer cells by promoting the G1/S transition (58). TNNT1 is also upregulated in endometrial carcinoma, which indicates that TNNT1 may be associated with the phenotype of invasive tumor cells, but the specific mechanism of this gene involved in endometrial carcinoma is not clear (59). Kuroda et al (60) have demonstrated that TNNT1 is highly expressed in cancer tissues of the cervix, colon, lungs, ovaries and testes. Similarly, Gu et al (61) reported that postoperative recurrence of gallbladder carcinoma was associated with TNNT1 expression, and TNNT1 has been identified as a favorable prognostic biomarker for this disease. The results of the analysis in the present study demonstrated that high expression of TNNT1 was significantly associated with poor OS in patients with ES, consistent with the above findings. To the best of our knowledge, no published studies on TNNT1 in ES are currently available. Therefore, further study of the molecular mechanism of TNNT1 in ES is required to identify more efficient and sensitive prognostic molecular markers for ES.

TTN encodes a large abundant protein of striated muscle and is involved heart and skeletal muscle diseases (62,63). However, a limited number of studies on TTN in cancer are currently available. TTN has been reported to be mutated frequently in several tumor types, including lung squamous cell carcinoma and lung and colon adenocarcinoma (64). Whole exome sequence data have been previously used to estimate the gene mutation rate of TTN, which identified TTN mutations in colorectal cancer, suggesting that TTN mutations may serve a specific role in the occurrence or progression of colorectal cancer (65). Study of low-abundance transcriptomes is crucial for determining the molecular mechanisms of tumor progression, and Bizama et al (66) identified TTN as a novel marker for advanced gastric cancer by analyzing low-abundance transcriptomes. Yang et al (67) conducted a genotyping study and demonstrated that TTN was a potential biomarker for predicting the clinical prognosis of patients with hepatitis B virus-associated hepatocellular carcinoma. The Cancer Genome Atlas-based aggregation analysis revealed that the missense mutation of TTN was associated with good prognosis in lung squamous cell carcinoma (68). However, no studies are currently available on TTN in ES, but the LASSO COX regression and OS analysis results of the present study demonstrated that high expression of TTN was associated with poor prognosis of patients with ES and may serve an important role in the development of ES. However, whether TTN can be used as a prognostic biomarker for ES requires verification in further experiments.

MYL3 encodes myosin light chain 3 and is associated with cardiomyopathy (69,70). A previous study has demonstrated that MYL3 combined with Ca2+ can promote muscle development and participate in the contractile process of striated muscle (71). In addition, MYL3 appears to be involved in strength development and fine co-ordination of muscle contraction (72). In a study on zebrafish, Meder et al (73) revealed that the phosphorylation of the C-terminal serine residue of MYL3 is of great significance for cardiac contraction. However, whether MYL3 also serves a certain role in cancer has not been confirmed. The results of the present study demonstrated that a high expression of MYL3 was associated with poor prognosis in patients with ES, and this result may help explain the molecular mechanism of ES.

TMOD1 is involved in the development, invasion and metastasis of various tumors, including esophageal cancer (74), meningioma (75), breast cancer (76) and acute lymphocytic leukemia (77). High expression of TMOD1 is a key biomarker for poor prognosis in patients with oral squamous cell carcinoma (78). TCAP is associated with the onset of breast cancer and is important for accurate diagnosis and prognosis assessment of breast cancer subtypes (79,80). Sawaki et al (81) analyzed mRNA expression in 15 gastric cancer cell lines and 262 surgically resected gastric tissues; the results demonstrated that the expression of TNNI2 in gastric tissue can be used as a specific biomarker for predicting peritoneal metastasis of gastric cancer. The conclusions of the aforementioned studies support the results of the present study. However, owing to the lack of studies on ES, the credibility of the results of the present study needs to be verified in further experiments.

Previous studies have demonstrated that TNNT1 is closely related to the invasion and metastasis of multiple types of tumors (56,58). High expression of TNNT1 is a prognostic biomarker for a variety of cancers (48,59,61). However, there are a limited number of studies on the molecular mechanism of TNNT1 in ES. The Rac1 pathway is associated with the malignant biological behavior of ES (40). Whether TNNT1 is involved in the Rac1 pathway is still unclear; however, it may be hypothesized that TNNT1 may participate in the invasion and metastasis of ES by activating the Rac1 pathway. Further molecular studies will be conducted in the future to validate these results.

The majority of bioinformatics studies focus on a single microarray data set and use a single method to analyze DEGs. In the present study, the raw data were derived from three mRNA microarray data sets, thus increasing the sample size and confidence. In addition, various analysis methods were applied to analyze the data in depth, which provided different perspectives. However, the present study has certain limitations. First, the data sets selected in the present study have certain heterogeneity. Although the batch data were removed and quality control and standardization were performed on the original data, a larger sample size and higher quality data sets are still required to verify the reliability of the results of this study. Second, the present study was a second mining and analysis of previously published data sets. Although the results of several previous studies are consistent with the present results, further biological experiments are needed to verify the function of the identified DEGs in ES (82–85). Further molecular studies will be conducted in the future to clarify the specific molecular mechanisms of these DEGs in ES.

In conclusion, a total of 211 DEGs were identified by bioinformatics analysis of three mRNA microarray data sets, and these DEGs may serve key roles in the development of ES. Through further analysis, 20 Hub genes were screened and three genes (TNNT1, TTN, and MYL3) were obtained through LASSO COX regression and OS analysis. These three genes may be potential prognostic biomarkers for ES. These results provided a theoretical basis for elucidating the molecular mechanisms underlying ES development and identifying candidate biomarkers for ES, which may help in developing new strategies for the diagnosis and treatment of ES.