Identification of potential diagnostic and therapeutic target genes for lung squamous cell carcinoma.

The purpose of this study was to identify potential molecular markers of lung squamous cell carcinoma (LUSC). Three datasets containing LUSC mRNA sequencing data were downloaded from the Gene Expression Omnibus, The Cancer Genome Atlas and the Gene Expression Profiling Interactive Analysis databases. These datasets were used to identify significantly differentially expressed genes (DEGs) in LUSC. A protein-protein interaction network of the DEGs was constructed followed by Gene Ontology, Kyoto Encyclopedia of Genes and Genomes and overall survival analyses of the DEGs. A total of 37 DEGs between LUSC and normal tissues were identified, including 26 downregulated genes and 11 upregulated genes. Biological Process enrichment analysis revealed that the DEGs were mainly enriched in 'cell adhesion', 'cell-matrix adhesion', 'anatomical structure morphogenesis', 'ECM-receptor interaction' and 'focal adhesion'. Overall survival analysis demonstrated that transcription factor 21, α-2-macroglobulin, acyl-CoA synthetase long chain family member 5, integrin subunit β8, meiotic nuclear divisions 1 and secretoglobin family 1A member 1 were significantly associated with the occurrence and development of lung cancer, and these genes were selected as hub genes. The results obtained in the present study may aid the elucidation of the molecular mechanisms involved in the development of LUSC and may provide potential targets for LUSC treatment.

Introduction

Lung cancer is the leading cause of cancer-associated mortality worldwide (1). There are two main histological types of lung cancer: Small cell lung cancer and non-small cell lung cancer (NSCLC). The latter is further subdivided into lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC) and large cell lung cancer. The incidence and mortality rate of LUSC are high, with >400,000 new cases occurring worldwide each year (2). Novel therapeutic agents for the treatment of lung cancer have been developed, including bevacizumab and epidermal growth factor receptor (EGFR) tyrosine kinase, ALK receptor tyrosine kinase (ALK) and CD274 molecule inhibitors (3). However, these agents are not effective for the treatment of LUSC, which accounts for ~25% of NSCLC cases (4). The elucidation of the molecular mechanisms underlying the development of LUSC may aid in the identification of new treatment strategies. Microarray technology and bioinformatics analysis have emerged as powerful tools for the study of different types of cancer and may facilitate the discovery of novel biomarkers and potential therapeutic targets (5,6). However, the high false positive rate of single microarray analysis may confound results (7). The current study analyzed three datasets containing LUSC mRNA sequencing data downloaded from the Gene Expression Omnibus (GEO; www.ncbi.nlm.nih.gov/geo), the Cancer Genome Atlas (TCGA; www.cancergenome.nih.gov) and the Gene Expression Profiling Interactive Analysis (GEPIA; http://gepia.cancer-pku.cn) databases. Genes in the intersection of the three datasets were regarded as significantly differentially expressed genes (DEGs). Subsequently, Gene Ontology (GO; http://geneontology.org), Kyoto Encyclopedia of Genes and Genomes (KEGG; http://david.ncifcrf.gov), protein-protein interaction network (PPI) and overall survival analyses of DEGs were performed. Receiver operating characteristic (ROC) curves were generated to identify genes of potential diagnostic and therapeutic value for LUSC. A total of 37 DEGs were screened, of which six were selected as the hub genes. These hub genes were significant for LUSC prognosis and may be candidate biomarkers for lung cancer.

Materials and methods

Microarray data

The GEO is an international public repository that stores and freely distributes microarray, second-generation sequencing and other forms of high-throughput functional genomic datasets (8). TCGA is a large-scale cancer genome project that provides researchers with multidimensional maps of the key genomic changes and clinicopathological information in 33 types of cancer (9). GEPIA is a newly developed interactive web server for analyzing RNA sequencing expression data based on 9,736 tumors and 8,587 normal samples from the TCGA and Genotype-Tissue Expression databases. The gene expression dataset GSE31552 (10) was downloaded from GEO and included 25 nontumor tissues and 25 tumor tissues. The GPL6244 Affymetrix Human Gene 1.0 ST Array (Affymetrix; Thermo Fisher Scientific, Inc., Waltham, MA, USA) was used. The TCGA dataset was derived from the ‘Protein-coding Transcripts’ of LUSC in the Cancer RNA-Seq Nexus (CRN; http://syslab4.nchu.edu.tw) database (11) and included 63 cancer samples and 51 normal samples. The GEPIA dataset was downloaded from the GEPIA online database including 486 cancer samples and 338 normal samples.

Data preprocessing and differential expression analysis

The original data were transformed into expression data using the affy package (http://www.bioconductor.org/packages/release/bioc/html/affy.html) in R (version 3.5.3; http://cran.r-project.org/bin/windows/base/rpatched.html) (12). Missing data were estimated with weighted K-nearest neighbors method (13), and all the expression profiles were normalized by the median normalization method (14). The t-test method of the limma package (http://master.bioconductor.org/packages/release/bioc/html/limma.html) in R was used to identify DEGs between LUSC and normal controls. In order to filter the DEGs of each dataset, set cut off criteria for each dataset. For GSE31552, DEGs were identified at P<0.05 and |log2FC|>1, where FC is fold change; For the TCGA and GEPIA datasets, DEGs were identified at P<0.05 and |log2FC|>2. Each probe name of GSE31552 was converted into a gene name using the hugene10sttranscriptcluster.db package (http://bioconductor.org/packages/release/data/annotation/html/hugene10sttranscriptcluster.db.html) in R and the corresponding platform file. The Spearman's correlation test was used to cluster the samples and to calculate the correlation coefficients between the samples, the results were subsequently visualized using pheatmap package (http://www.bioconductor.org/packages/release/bioc/html/heatmaps.html) in R. To verify the rationality of the DEGs of GSE31552, draw the volcano plot using ggplot2 package in R (http://cran.r-project.org/web/packages/ggplot2). Heatmaps of the DEGs (based on |log2FC|) were generated using the heatmap.2 package (http://www.rdocumentation.org/packages/gplots/versions/3.0.1.1/topics/heatmap.2) in R to present the DEGs in each sample.

Screening and analysis of DEGs

Significant DEGs between LUSC and normal samples were screened using FunRich (version 3.1.3), which is an open-access standalone functional enrichment and interaction network analysis tool (15). A PPI network of DEGs was constructed using the Gene MANIA database (http://www.genemania.org). The Gene MANIA database is a useful tool for generating hypotheses about gene function, analyzing gene lists, prioritizing functionally analyzed genes and reporting weights (16). Subsequently, a hierarchical clustering of DEGs was constructed using an online analysis database of University of California, Santa Cruz Xena (UCSC Xena 2.0; http://xena.ucsc.edu/welcome-to-ucsc-xena) (17).

Functional enrichment analysis

GO is a common tool for annotating genes and their products (18). KEGG is a knowledge base for the systematic analysis of gene functions in terms of networks of genes and molecules (19). The Database for Annotation, Visualization and Integrated Discovery 6.8 (DAVID 6.8; http://david.ncifcrf.gov) (20) is a comprehensive database providing a complete set of functional annotation information of genes and proteins from which researchers can extract biological information. To analyze the DEGs at the functional level, GO and KEGG pathway enrichment analyses were performed using the DAVID online tool to obtain the enriched biological processes and pathways. P<0.05 was considered to indicate a statistically significant difference.

Hub gene selection and analysis

Overall survival analyses of DEGs were performed using GEPIA (21) and Kaplan-Meier plotter (version 2018.11.04; http://www.kmplot.com) (22). Genes with significant differences for the prognosis of LUSC were selected as hub genes (P<0.05). Analyzed the expression of Hub genes in different subtypes of lung cancer and different stages of LUSC based on GEPIA online analysis database (21). To understand the diagnostic value of the hub genes, ROC curves were plotted using the pROC package (pROC_1.12.1; http://master.bioconductor.org/packages/release/bioc/html/pROC.html) in R based on the GEO dataset. Through comprehensive analysis of Hub genes survival analysis, Hub gene expression in different lung cancer subtypes and LUSC stages, ROC curve of Hub genes and literature search. Firstly, six genes closely related to the survival of LUSC patients were screened using survival analysis of Hub genes. The expression of these six genes was validated with LUSC stages and in patients with LUAD and LUSC. The diagnostic efficacy of ROC was also evaluated, and the three methods revealed that secretoglobin family 1A member 1 (SCGB1A1) was closely related to survival and stages, and the diagnostic efficiency of SCB1A1 ROC was the highest among the screened genes. Therefore, SCB1A1 was selected as the target gene. The expression of SCGB1A1 in different LUSC cell lines based on The Cancer Cell Line Encyclopedia (CCLE; http://portals.broadinstitute.org/ccle) database was plotted using GraphPad Prism (version 7; GraphPad Software, Inc., La Jolla, CA, USA) (23). The Human Protein Atlas (Version 18.1; www.proteinatlas.org) (24) was used to show the expression of SCGB1A1 in different tissues. The Oncomine database (Version 4.5; http://www.oncomine.org) (25) was used to analyze Okayama (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE31210) and Beer datasets (26) to identify the associations between SCGB1A1 and EGFR mutation status, echinoderm microtubule associated protein-like 4 (EML4) and anaplastic lymphoma kinase (ALK) gene fusion, expression of TP53 and smoking status.

Results

Data preprocessing

Following the evaluation of missing data and normalization, the expression profiling data were plotted (Fig. 1A). The similar levels of data points indicate high consistency showing high accuracy of classification and comparison. The baseline level of Samples Clustering Analysis (Fig. 1B) revealed that the sample sources were reliable. For dataset GSE31552, a total of 1,712 DEGs were identified at P<0.05 and |log2FC|>1, and the rationality of the values was verified by volcano plots (Fig. 1C), in which red dots represent upregulated genes, green dots represent downregulated genes, and black dots represent unchanged genes. Hierarchical clustering analysis revealed that the gene expression patterns of the DEGs were similar among the array data of GSE31522 (Fig. 1D), indicating that the molecular changes in LUSC are consistent and may represent a novel genetic signature in LUSC.

Figure 1.

Data preprocessing and differential expression analysis. (A) A box-plot of standardized expression data of normal and LUSC tissues. The red symbols denote normal samples, and the blue symbols represent patients with LUSC. The red diamond in each frame represents the average level of gene expression in each sample. (B) Heatmap of the common DEGs between LUSC tissues and normal controls in the GSE31552 dataset. Blue represents downregulation and red represents upregulation. The normal groups were significantly separated from the LUSC groups, indicating that the sample source data are reliable. (C) Volcano plot used to verify the rationality of the P-values and corresponding fold changes of the expression data. The red circles denote upregulation, and the green symbols denote downregulation. (D) A hierarchical clustering diagram of the DEGs. Blue represents downregulation and red represents upregulation. LUSC, lung squamous cell carcinoma; DEGs, differentially expressed genes; down, downregulated genes; up, upregulated genes; not, not differentially expressed.

Identification of DEGs in LUSC

Following standardization of the microarray results, DEGs (2,162 in GSE31552; 1,842 in TCGA database; and 1,691 in the GEPIA database) were identified. The three datasets shared 37 genes that were considered to be significant DEGs between LUSC and normal tissues (Fig. 2A). These genes included 26 downregulated and 11 upregulated genes.

Figure 2.

Venn diagram, PPI network and hierarchical clustering of significant DEGs. (A) Venn diagram drawn in FunRich (version 3.1.3). A total of 37 common significant DEGs were obtained. (B) PPI network of DEGs constructed using Gene MANIA. Red circles present DEGs and green circles present predicted genes. The size of the circle represents the level of connectivity. Purple lines represent co-expression, green lines represent genetic interactions and pink lines represent physical interactions. (C) Hierarchical clustering of DEGs was performed using the UCSC Xena database. Upregulation of genes is indicated by red and downregulation of genes is indicated by blue. The expression of DEGs in different grades and different sample types is presented (Note: TMPRSSIIE gene is not included in the lung cancer dataset of UCSC Xena database). PPI, protein-protein interactions; GEPIA, Gene Expression Profiling Interactive Analysis; TCGA, The Cancer Genome Atlas; LUSC, lung squamous cell carcinoma; DEGs, differentially expressed genes; UCSC, University of California, Santa Cruz.

KEGG and GO enrichment analyses of DEGs

To analyze the biological classifications of the DEGs, DAVID was used for functional and pathway enrichment analysis. The identified GO terms and pathways are presented in Table I. The GO terms enriched by DEGs were mainly associated with ‘cell adhesion’, ‘cell-matrix adhesion’ and ‘anatomical structure morphogenesis’, whereas the pathways enriched for DEGs were associated with ‘ECM-receptor interaction’ and ‘focal adhesion’.

Table I.

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis of differentially expressed genes in lung squamous cell carcinoma samples.

Term	No. enriched genes	P-value
GO:0007155 ‘cell adhesion’	7	0.0003
GO:0007160 ‘cell-matrix adhesion’	3	0.0151
GO:0009653 ‘anatomical structure morphogenesis’	3	0.0157
GO:0048333 ‘mesodermal cell differentiation’	2	0.0227
GO:0031581 ‘hemidesmosome assembly’	2	0.0247
GO:0030301 ‘cholesterol transport’	2	0.0329
GO:2000811 ‘negative regulation of anoikis’	2	0.0348
GO:0019915 ‘lipid storage’	2	0.0489
hsa04512 ‘ECM-receptor interaction’	4	0.0014
hsa04510 ‘focal adhesion’	4	0.0153

PPI network construction and module analysis

The PPI network of DEGs constructed using Gene MANIA (Fig. 2B) revealed the co-expression, genetic interactions and physical interactions among the DEGs and predicted genes. Red circles represent DEGs obtained by the above data set intersection; green circles represent predicted genes obtained from the above DEGs through the plug-in Gene MANIA. Dense networks indicated that DEGs were closely related to each other. In the hierarchical clustering plot (Fig. 2C), the DEGs distinguished the LUSC samples from the normal samples.

Screening for hub genes

To ensure the accuracy of the analyses, overall survival analysis of the DGEs was performed using GEPIA and Kaplan-Meier plotters. Genes were selected if they are meaningful (P<0.05) in both analyses. Six hub genes were ultimately selected based on their significant effects on survival (Fig. 3). The survival curves indicated that integrin subunit β8 (ITGB8) and SCGB1A1 were positively associated with overall survival, whereas TCF21, A2M, ACSL5 and meiotic nuclear divisions 1 (MND1) were negatively associated with overall survival. The hub genes and their functions are presented in Table II.

Figure 3.

Overall survival analysis. (A-C) Overall survival analyses were performed using the GEPIA online platform and (D-F) Kaplan-Meier plotter online platform. The solid line represents the survival curve and the dotted line represents the 95% confidence interval. Log-rank P<0.05 was considered to indicate a statistically significant difference. Patients with expression above the median are indicated by red lines, and patients with expression below the median are indicated by (A-C) blue lines or (D-F) black lines. (A) TCF21 was negatively associated with the overall survival of patients with LUSC; (B) A2M was negatively associated with the overall survival of patients with LUSC; (C) ACSL5 was negatively associated with the overall survival of patients with LUSC; (D) ITGB8 was positively correlated with the overall survival of patients with LUSC; (E) MND1 was negatively associated with the overall survival of patients with LUSC; (F) SCGB1A1 was positively correlated with the overall survival of patients with LUSC. TCF21, transcription factor 21; A2M, α-2-macroglobulin; ACSL5, acyl-CoA synthetase long chain family member 5; ITGB8, integrin subunit beta 8; SCGB1A1, secretoglobin family 1A member 1; MND1, meiotic nuclear divisions 1; HR, hazard ratio; TMP, transcripts per million; LUSC, lung squamous cell carcinoma.

Table II.

Function of the hub genes (GeneCards; http://www.genecards.org).

Gene	Gene name	Function
TCF21	Transcription factor 21	RNA polymerase II transcription factor activity, sequence-specific DNA binding
A2M	α-2-Macroglobulin	Tumor cell adhesion, migration and growth
ACSL5	Acyl-CoA synthetase long chain family member 5	Pro-apoptotic sensing of enterocytes
ITGB8	Integrin subunit β8	Cell-adhesion to extra cellular matrix or to other cells
MND1	Myeloid cell nuclear differentiation antigen	Meiosis and recombination
SCGB1A1	Secretoglobin family 1A member 1	Inhibition of phospholipase A2 and arachidonic acid release, prostaglandin D2 receptor antagonism in the lung

AA, arachidonic acid.

Hub gene analysis

Six genes were identified as hub genes and had significant effects on survival. The expression levels of the hub genes in different subtypes of NSCLC and different stages of LUSC were analyzed (Fig. 4). Compared with healthy controls, TCF21, A2M and SCGB1A1 were downregulated in LUSC and LUAD, ITGB8 and MND1 were upregulated in LUAD and LUSC, and ACSL5 was downregulated in LUSC and upregulated in LUAD. In addition, expression of TCF21, ITGB8, MND1 and SCGB1A1 was significantly different among different stages of LUSC. SCGB1A1 was closely associated with LUSC stage [Pr(>F)=8.06×10−7]. ROC curves revealed that the hub genes had a good diagnostic value for LUSC [area under the curve (AUC)>0.7], particularly SCGB1A1 (AUC=0.922; Fig. 5A).

Figure 4.

Box and violin plots were produced using the Gene Expression Profiling Interactive Analysis online platform. The boxplots present the expression of hub genes in LUSC and LUAD; red, tumor tissue; gray, normal tissue. Violin plots demonstrate the association between hub gene expression and LUSC staging. (Aa) TCF21 was downregulated in LUSC and LUAD; (Ab) TCF21 was associated with LUSC stages [Pr(>F)=0.0045]. (Ba) A2M was downregulated in both LUSC and LUAD; (Bb) A2M was not significantly different among the different stages of LUSC; (Ca) ACSL5 was downregulated in LUSC but upregulated in LUAD; (Cb) ACSL5 was not significantly different among the different stages of LUSC. (Da) ITGB8 was upregulated in LUAD and LUSC; ITGB8 was associated with stages of LUSC [Pr(>F)=0.0388]; (Ea) MND1 was upregulated in LUAD and LUSC; MND1 was associated with stages of LUSC [Pr(>F)=0.000953] (Fa) SCGB1A1 was downregulated in both LUSC and LUAD; SCGB1A1 was closely related to stages of LUSC [Pr(>F)=8.06×10-7]. *P<0.01 tumor group vs. control group. LUSC, lung squamous cell carcinoma; LUAD, lung adenocarcinoma; TCF21, transcription factor 21; A2M, α-2-macroglobulin; ACSL5, acyl-CoA synthetase long chain family member 5; ITGB8, integrin subunit beta 8; SCGB1A1, secretoglobin family 1A member 1; MND1, meiotic nuclear divisions 1; T, tumor; N, normal.

Figure 5.

ROC curves of hub genes and SCGB1A1 expression in different tumors, cell lines and tissues. (A) ROC curves of hub genes were plotted based on the Gene Expression Omnibus dataset. The closer AUC value is to one, the higher the diagnostic value of the gene. (B) SCGB1A1 gene expression profiles of 20 common lung cancer cell lines based on the CCLE database. (C) Dot plot of SCGB1A1 gene expression profile across different tumor samples and paired normal tissues. Each dot represents sample expression; red denotes tumor samples and green denotes normal samples. (D) SCGB1A1 mRNA and protein expression in normal human tissues based on The Human Protein Atlas. ROC, receiver operating characteristic; SCGB1A1, secretoglobin family 1A member 1; MND1, meiotic nuclear divisions 1; TCF21, transcription factor 21; A2M, α-2-macroglobulin; ACSL5, acyl-CoA synthetase long chain family member 5; ITGB8, integrin subunit beta 8; CCLE, Cancer Cell Line Encyclopedia; GEPIA, Gene Expression Profiling Interactive Analysis; TPM, transcripts per million; AUC, area under the curve; T, tumor; N, normal.

In summary, downregulation of SCGB1A1 was significantly associated with overall survival of LUSC patients (log-rank P=0.011; Fig. 3F), and SCGB1A1 downregulation was significantly associated with LUSC stages [Pr(>F)=8.06×10−7; Fig. 4C]. Additionally, ROC curve analysis demonstrated that SCGB1A1 had high diagnostic value (AUC=0.922). Thus, SCGB1A1 was screened as a target gene for LUSC. Further analysis based on GEPIA Database revealed that SCGB1A1 expression was lower in LUSC compared with normal lung tissue (Fig. 5C). The Human Protein Atlas database analysis demonstrated that SCGB1A1 was overexpressed in normal lung tissues compared with other tissues (Fig. 5D). The expression of SCGB1A1 in lung cancer cell lines was analyzed using the CCLE online platform; the results revealed that SCGBA1 was downregulated in the majority lung cancer cell lines (Fig. 5B). Oncomine analysis of Okayama and Beer datasets revealed that lower mRNA levels of SCGB1A1 were associated with EGFR mutation, EML4-ALK fusion, expression of TP53 and smoking (Fig. 6A-D).

Figure 6.

Associations between the expression of SCGB1A1 and (A) EGFR mutation, (B) EML4-ALK gene fusion, (C) expression of TP53 and (D) smoking in the Okayama and Beer dataset of the Oncomine database. No value, patients without the phenotype. SCGB1A1, secretoglobin family 1A member 1; EGFR, epidermal growth factor receptor; EML4, microtubule associated protein-like 4; ALK, anaplastic lymphoma kinase; IHC, immunohistochemistry.

Discussion

Lung cancer is a leading cause of cancer-associated mortality worldwide (1). The five-year relative survival rate of lung cancer is lower than the corresponding rates of breast cancer, colon cancer and kidney cancer (27). This low rate is partly due to a late-stage diagnosis in >50% of cases (28). Although chemotherapy and improvements in supportive care have improved overall survival and quality of life, the prognosis of patients with advanced NSCLC remains poor (27). Individualized targeted therapies for lung cancer are being investigated, and progress has been applied to the clinic (29), including in TP53, EGFR, KRAS proto-oncogene, GTPase, EML4-ALK rearrangement and MET signal transduction targeted therapies (30). The majority of these treatments are effective in LUAD and no front-line targeted therapies are currently clinically available for LUSC (31,32). Therefore, effective diagnostic and therapeutic markers are required for LUSC. Microarray technology may allow the identification of the genetic changes implicated in LUSC and is an effective method for identifying new biomarkers in other diseases, such as vascular diseases (33) and oral squamous cell carcinoma (34).

The present study analyzed three datasets, including GSE31552, and TCGA and GEPIA datasets, and revealed DEGs between LUSC and normal tissues. A total of 37 DEGs were identified, including 26 downregulated genes and 11 upregulated genes. To understand the functions and associations of these DEGs, GO and KEGG analyses were performed. The DEGs were mainly enriched in ‘cell adhesion’, ‘cell-matrix adhesion’, ‘anatomical structure morphogenesis’, whereas KEGG pathway enrichment was mainly concentrated in ‘ECM-receptor interaction’ and ‘focal adhesion’. Previous studies have reported that cell adhesion, cell-matrix adhesion and anoikis are closely associated with tumorigenesis and development (35–37). Additionally, studies of anatomical structure and morphogenesis have revealed cholesterol levels to be associated with tumor metastasis and invasion (38–40). These studies indicated that DEGs may be closely related to the occurrence, development, metastasis and invasion of tumors. The pathways of ECM-receptor interaction and focal adhesion are important mediators of cell adhesion, growth, proliferation, survival, angiogenesis and migration (41,42). The results obtained in the aforementioned studies are consistent with the results obtained in the current study. The enriched modules and pathways identified in the current study may have genetic effects on LUSC and the identified genes may interact within a network.

Survival analysis was performed using the DEGs and six genes with significant impacts on overall survival were identified as hub genes. SCGB1A1, also known as the Clara cell secretory protein (43,44), had a good diagnostic value for LUSC (AUC=0.922) and was significantly associated with poor overall survival and tumor stage. Previous studies have demonstrated that SCGB1A1 demonstrated anti-inflammatory, immunomodulatory and antitoxin properties (45–49). SCGB1A1 mRNA expression is abundant throughout the airway (50), however, its expression level is very low in lung tumors (51), which is consistent with the results obtained in the current study. Expression levels of SCGB1A1 are decreased in diffuse goblet cell hyperplasia and squamous metaplasia, whereas alveolar expression is elevated in alveolar cell bronchitis (52). SCGB1A1 knockout mice are more susceptible to lung injury (by bacterial or viral infection, ozone and cigarette smoke) or sensitization, exhibit increased inflammation and remodeling reactions, have more frequent lung tumors, and have a stronger T-helper 2-directed immune response compared with control mice (53). In the present study, SCGB1A1 was associated with tumor stage, EGFR mutation, ALK gene fusion and smoking history. A previous study has revealed that EGFR tyrosine kinase inhibitors have the same effect on the prognosis of patients with LUSC as chemotherapy, with fewer complications and higher quality of life (54). Therefore, SCGB1A1 may serve a protective role in lung tissue. Inducing SCGB1A1 expression may inhibit the expression of c-MYC and C-RAF, which may further inhibit the metastasis of tumors (55). However, the mechanism remains unclear, and future studies are required to elucidate the pathways involved in SCGB1A1 and lung cancer.

TCF21, located on chromosome 6q23-q24, encodes a basic helix-loop-helix transcription factor essential for epithelial cell differentiation (56,57). It can be readily methylated and subsequently cause tumorigenesis (58,59). A previous study has indicated that hypermethylation and decreased expression of TCF21 are tumor-specific and are frequently observed in NSCLC (60). The protein encoded by A2M is a protease inhibitor and cytokine transporter (61). A previous study has revealed that a progressive A2M deficiency may promote tumor development in nude mice (62). A2M regulated tumor cell adhesion, migration and growth by inhibiting tumor-promoting signaling pathways, including the phosphoinositide 3-kinase/protein kinase B (PI3K/AKT) pathway and mothers against decapentaplegic homolog (SMAD) and upregulating phosphatase and tensin homolog via downregulation of microRNA-21 in vitro and in tumor xenografts (62). The level of A2M in human blood decreases with age (63). ACSL5, a mitochondria-localized enzyme that catalyzes the synthesis of long-chain fatty acid thioesters, is physiologically involved in the induction of apoptosis in intestinal cells (64). Studies have revealed that ACSL5 isozymes serve leading roles in the biosynthesis of mitochondrial cardiolipin and may participate in the survival of cancer cells (65–67). ITGB8, a member of the integrin β chain family, is increased in different types of cancer, including breast, lung, throat and stomach cancer (68). High expression of ITGB8 serves an important role in the metastasis of human lung cancer cells. When ITGB8 is silenced, the expression of E-cadherin and cystatin B is increased, whereas the expression of C-X-C motif chemokine ligand CXCL1, CXCL2, CXCL5, matrix metalloproteinase (MMP)-2 and MMP-9 is decreased (69). Furthermore, changes in the cell cycle, the expression of metastasis-associated genes and metastatic potential may be accompanied by decreased tumor cell signal transduction and molecular activity (69,70). The products of the MND1 gene bind to PSMC3 interacting protein to form stable heterodimer complexes that bind to DNA and stimulate the activities of RAD51 recombinase and DNA meiotic recombinase 1, which are required for meiotic recombination (71). MND1 was significantly upregulated in ovarian cancer compared with ovarian tissue samples from healthy controls (72). However, to the best of our knowledge, MND1 upregulation has not been previously reported in human lung cancer.

The hub genes identified in the current study are associated with the occurrence and development of tumors. The involvement of TCF21 and ITGB8 in lung cancer has been previously documented. However, there are fewer reports of SCGB1A1, A2M, ACSL5 and MND1 in lung cancer. In vitro overexpression of TCF21 may inhibit tumor growth and chemoresistance possibly through the AKT signaling pathway (73,74). Upregulation of ITGB8 may promote the expression of tumor metastasis genes and enhance the invasive ability of tumor cells in LUSC by regulating the phosphorylation levels of mitogen-activated protein kinase/extracellular signal-regulated kinase and AKT. An increased incidence of lung injury and lung tumors was reported following SCGB1A1 knockout (52). A previous study has reported that SCGB1A1 may serve an anti-inflammatory role by inhibiting phospholipase A2 (75); therefore, the downregulation of SCGB1A1 may lead to an imbalance of T cell subsets, which in turn affects the antitumor activity of serum peripheral blood mononuclear cells in patients with LUSC (76). Transcriptome analysis of A2M-treated tumor cells, xenografts and mouse liver revealed that A2M modulates tumor cell adhesion, migration and proliferation by inhibiting tumor-promoting signaling pathways, such as PI3K/AKT and SMAD, and by upregulating PTEN via downregulation of miR-21 in vitro and in tumor xenografts (77). ASCL5 is closely associated with cancer cell apoptosis (64). The hub genes in the current study were associated with poor overall survival rates, and ROC curves revealed high diagnostic values (AUC>0.7). The results obtained in the current study suggest that these genes may serve important roles in the occurrence and development of LUSC and may be used as biomarkers for the diagnosis of LUSC.

In conclusion, the purpose of the current study was to identify genes that may be involved in the development or progression of LUSC. A total of 37 DEGs were identified, of which 6 were identified as hub genes and may be used as biomarkers for the diagnosis and prognosis evaluation of LUSC. However, the results of this study were obtained through big data analysis, and validation of the results via animal experiments and clinical trials is required.