Identification and Validation of Anoikis-Related Signatures for Predicting Prognosis in Lung Adenocarcinoma with Machine Learning

Introduction

Lung cancer is currently the second most common cancer in China, accounting for the highest rate of cancer deaths.^1,2 Lung adenocarcinoma (LUAD) is the most frequent type of lung cancer, with advanced lung adenocarcinoma being prone to metastasize, resulting in poor patient prognosis. Methods of treating LUAD include surgery, radiotherapy, chemotherapy, targeted therapy, and treatment with immune checkpoint inhibitors. Tumor recurrence and metastasis are the most frequent causes of treatment failure and death. Multiple genes and genetic pathways are thought to be involved in the development of tumor metastases.

Anoikis is a specialized form of programmed cell death that plays important roles in body development, tissue self-balance, disease occurrence and tumor metastasis.³ During the process of tumor metastasis, tumor cells that adhere to the extracellular matrix and other cells detach from their surroundings, and survive in the circulation through paracrine and paracrine mechanisms. These cells resist apoptosis, and subsequently regain the ability to attach to other sites in the body, resulting in tumor spread, metastasis and invasion.^4–6 Resistance to anoikis during tumor metastasis allows tumor cells to spread through the circulatory system to distant organs.⁷ Studies on a variety of tumors have shown that resistance to anoikis can predict the effects of treatment and patient prognosis. For example, anoikis-related genes (ANRGs) were found to play a critical role in the progression of cancer and metastatic cascades.^8,9 To date, however, clinical prognostic models based on genes associated with anoikis have not been established in LUAD. The process of development of LUADs may be better understood by investigating the mechanisms of action of ANRGs.

A total of 316 anoikis-related genes (ANRGs) were integrated from the GeneCard database¹⁰ (http://www.genecards.org/) and Harmonizome portals,¹¹(Table S1). The present study analyzed alterations of ANRGs in LUAD, showing that these alterations were prognostic of patient outcomes. Our focus is on transcriptomics and somatic mutations are the focus of our next phase of research, so they are described briefly in this article. In addition, XG-boost was used to identify the most important genes associated with LUAD related death. Because these findings suggested that ITGB4 was an oncogene associated with the progression of LUAD, the levels of expression of ITGB4 was analyzed in an LUAD tissue microarray.

Materials and Methods

Datasets

RNA-sequencing transcriptome data for LUAD were obtained from The Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/GSE31210) and The Cancer Genome Atlas (TCGA, http://portal.gdc.cancer.gov/)¹² The RNA sequence transcriptome profiles and clinicopathological data from 535 LUAD tissue samples and 59 normal tissue samples were collected from the LUAD cohort in the TCGA (TCGA-LUAD), and data on 226 tumor and 20 non-tumor tissue from the GSE31210 dataset were obtained through the GEO database.

Pathway Enrichment Analysis

Correlations between the levels of expression and functions of eight anoikis regulators were evaluated by Spearman correlation analysis, with interplay between regulators estimated using the String database tool (https://string-db.org). A complete evaluation of Gene characteristics were evaluated using Metascape (http://metascape.org). Targetable mechanisms of action of anoikis regulators in LUAD were evaluated by Gene Set Enrichment Analysis (GSEA), Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG).

Construction and Validation of the Risk Score

Factors highly prognostic of survival were evaluated using a Cox regression model of the least absolute shrinkage and selection operators (LASSO), which selected eight regulators of anoikis. These eight genes were analyzed by Cox regression analysis, with a P-value <0.2 set as the cutoff value to avoid omissions. The risk score for each patient in the training (TCGA) and validation (GSE31210) datasets was calculated using the formula:

in which coefi represents the coefficient and xi represents the expression value of each selected anokis regulator.

Machine Learning Algorithm

The mechanism by which anoikis regulators affected the risk scores was determined using an additive explanation method with the XG-Boost machine learning model. This machine learning algorithm assembles potential prediction models to yield a fine-tuned predictive result. Because the algorithm involved parallel, distributed, out-of-core, and cache-aware computing, it was more than ten times faster than popular machine learning (ML) and deep learning models.

Tissue Samples and Immunohistochemistry

Samples of LUAD tumors and adjacent normal tissue were collected and assessed. These samples were embedded in paraffin, dewazed, separate and hydrated according to standard protocols. The tissue samples were mounted on slides, heated in Tris-EDTA buffer, and blocked by overnight incubation with 5% gout serum. The samples were subsequently incubated overnight with primary anti- ITGB4 antibody (Proteintech, Cat No.21738-1-AP, Wuhan, Hubei, China), followed by incubation with secondary antibody. Slides were viewed by light microscopy, and immunohistochemical staining was analyzed semiquantitatively by two independent, experienced pathologists. Each sample was scored according to the percentage of positive cells, with 1+ indicating no staining, 2+ indicating weak staining, 3+ indicating moderate staining, and 4+ indicating strong staining. All experiments with human samples were approved by the Ethics Committee of the First Affiliated Hospital of Zhengzhou University. This study abode by the Declaration of Helsinki and good clinical practice guidelines. Written informed consent to participate was waived for the use of the clinical samples as retrospective samples were used in this study. All the patients’ information was made anonymous in this study.

Statistical Analysis

All statistical analyses were performed using the R statistical software package (3.5.3). Continuous variables in two groups were compared using Student’s t-tests (unpaired, two-tailed), whereas comparisons in three or more groups were performed using one-way analysis of variance (ANOVA) and Kruskal–Wallis tests. The relationships between the patterns of expression of anoikis regulators and scientific features were analyzed by Chi-square tests. Survival was assessed using the Kaplan–Meier method, with survival as a function of each ectopically expressed gene compared by Log rank tests. The relationships between patient variables and prognosis were evaluated by univariate and multivariable Cox regression analyses. P-values were corrected for multiple comparisons using the Benjamini–Hochberg (BH) method and used to predict medical responses to immune checkpoint blockade. In all statistical analyses, P <0.05 was considered statistically significant.

Results

Identification of Anoikis-Related Genes Associated with Prognosis

The levels of expression of anoikis-related genes were obtained from the TCGA dataset, with the differences in the expression of these genes in LUAD tumor and pericarcinomous tissue samples analyzed by principal component analysis (Figure 1A). The levels of expression of 311 of these genes differed significantly in tumor and adjacent normal tissue samples, with heatmaps showing the same trend (Figure 1B). Volcano plots showed that the expression of 28 DEGs were significantly higher and the levels of 10 significantly lower in LUAD than in normal tissues (Figure 1C). Further univariate regression analysis also showed that 15 of the 28 genes significantly upregulated in LUAD were closely associated with overall survival (p < 0.01) in both the training (TCGA) and validation (GSE31210) cohorts (Figure 1D). Among all the 557 patients, 86 experienced anoikis-related genes mutations with a frequency of 15.44%, indicating a high frequency of somatic mutations of the anoikis-related gene in LUAD (Figure S1).

Figure 1 Characteristics and differences of anoikis-related regulators in LUAD. (A) Principal component analysis of tumor and non-tumor tissue. (B) Heatmap showing the differential expression of anoikis-related regulators. (C) Volcano plot of differentially expressed genes (DEGs) in LUAD. (D) Effects of 15 anoikis-related genes on overall survival in the TCGA and GEO datasets.

Prognostic Value of Regulators Related to Anoikis and the Construction of a Risk Score Model

The 28 significantly differentially expressed anoikis-related regulators and the 15 anoikis-related genes that had a significant effect on the prognosis of patients with LUAD, were used to develop a prognostic risk model. A LASSO Cox regression model using the minimum criteria to determine whether anoikis modification could accurately predict outcomes in individual patients with LUAD was found to include the PLK1, SHC1, ITGB4, ANGPTL4, SLCO1B3, SPIB, KL and DAPK2 genes (Figure 2A and B). The ANRG risk score for each patient with LUAD in the training (TCGA) and validation (GSE31210) datasets was calculated using the formula: ANRG risk score = 0.075503×PLK1 +0.002077×SHC1 +0.017577×ITGB4 +0.012239×ANGPTL4 +0.013272×SLCO1B3 −0.01686×SPIB −0.01809×KL −0.09847×DAPK2. In addition, Kaplan-Meier analysis showed that higher risk scores in both the training (TCGA) and validation (GSE31210) datasets were associated with poorer OS (Figure 2C and D), indicating a close relationship between ANRG risk score and the prognosis of patients with LUAD.

Figure 2 Identification of an anoikis-related prognostic signature. (A) Identification of eight anoikis-related genes associated with prognosis of patients with LUAD by LASSO analysis with 10-fold cross validation. (B) Coefficient profile plots of eight prognostic anoikis-related genes. (C and D) Kaplan-Meier analyses showing differences in overall survival in subgroups of patients with LUAD.

Correlation of Anoikis Risk Scores with Clinical Characteristics

The optimal anoikis risk score cutoff values for patients with LUAD was calculated, and patients were divided into those with a low and high anoikis scores (Figure 3A). Correlations between anoikis scores and various clinicopathological characteristics were subsequently analyzed in patients with LUAD. ANRG risk scores differed significantly between between groups, with patients in the TCGA dataset stratified by PFS, tumor status, pathologic stage and OS (Figure 3B–E). Univariate and multivariable Cox regression analyses of the TCGA data identified ANRG risk scores as an independent prognostic factor associated with the outcome of patients with LUAD (Figure 3B), with anoikis risk scores significantly correlated with several clinical characteristics and patient prognosis.

Figure 3 Association between ITGB4 expression and key characteristics of patients with LUAD. (A): Relationships between ANRGs risk scores and characteristics of patients with LUAD. (B–E): Relationships between ANRGs risk scores and medical characteristics of patients with LUAD. ****p<0.0001.

Importance of Anoikis-Related Regulators to Risk Scores Determined by Machine Learning

The XGBoost algorithm was used to construct a classifier that could recognize the outstanding contributors to risk score. Evaluation by Shapley additive explanatory values (SHAP) showed that the four leading regulators were DAPK2, ITGB4, PLK1, and SHC1 (Figures 4A). The scores of all samples were summed and their features were sorted, with risk scores showing that DAPK2 and ITGB4 were the most important regulators (Figures 4B). Previous studies have suggested that DAPK2 contributes to the process of autophagy. Because less is known about ITGB4, the expression of ITGB4 was further analyzed and verified.

Figure 4 Contributions of anoikis regulators to machine learning risk scores. (A): Significance of anoikis regulators on chance scores using the LASSO algorithm. (B): Individual SHAP values of PLK1, ALDH2, PANX1, NR2C2, TLR4, and OUTLIN.

Expression of ITGB4 Protein in Two Cohorts of Tissue Samples

IHC analysis of the ZZU-cohort obtained from the First Affiliated Hospital of Zhengzhou University was performed to analyze the levels of ITGB4 protein in representative LUAD and peritumor tissues (Figure 5A). The IHC score of ITGB4 was significantly higher in LUAD than in paired normal lung tissues (p < 0.01, Figure 5B and C). Moreover, ITGB4 expression level differed significantly between groups, with the TCGA dataset stratified by OS, tumor status, and T stage (Figure 5D–F). These results confirmed at the protein level that ITGB4 is upregulated in LUAD.

Figure 5 Immunohistochemical (IHC) analysis of components of the anoikis-related prognosis signature in LUAD and normal lung tissue. (A): Level of expression of ITGB4 in LUAD and adjacent normal lung tissue. (B): Statistical analysis of the differences in scores. (C): IHC images with different ratings. (D): Correlation between the expression of ITGB4 and overall survival. (E): Relationship between ITGB4 expression level and T stage. (F): Relationship between ITGB4 expression level and tumor stage. 1+ indicating no staining, 2+ indicating weak staining, 3+ indicating moderate staining, and 4+ indicating strong staining.*p<0.05, ****p<0.0001.

Functional Enrichment Analysis with ITGB4

GO and KEGG analyses were also performed to further elucidate the role of ITGB4 in LUAD (Figure 6A and B). These analyses showed that ITGB4 participated in the cell cycle and energy metabolism, including in four pathways that play an important role in common pathway analysis: HALLMARK E2F TARGETS, HALLMARK MYC TARGETS V1, HALLMARK MYC TARGETS V2, HALLMARK OXIDATIVE PHOSPHORYLATION (Figure 6C). GSEA performed to identify the signals differentially activated in groups with high versus low ITGB4 expression (Figure 6D–F) showed that the dataset was highly enriched in genes involved in the four above pathways, further confirming that ITGB4 may be involved in the development of LUAD.

Figure 6 ITGB4-related gene enrichment analysis. (A): GO term enrichment analysis of ITGB4-binding genes. (B): KEGG term enrichment analysis of ITGB4-binding genes. (C): Ingenuity analysis of the ITGB4-related pathway. (D): GSEA analysis of the E2F TARGETS pathway. (E): GSEA analysis of the MYC TARGETS V1 and V2 pathways. (F): GSEA analysis of the HALLMARK OXIDATIVE PHOSPHORYLATION pathway.

Discussion

The rapid development of LUAD and the high incidence of distant metastases have made it difficult to improve patient survival using a single route or treatment. Construction of prediction models for transfer-related genes can provide a more effective method for early intervention. The predictive ability of these models is enhanced by the use of multiple biomarkers.

The present study found that a robust risk score could be generated using eight genes: PLK1, SHC1, ITGB4, ANGPTL4, SLCO1B3, SPIB, and DAPK2. Genes related to anoikis were found to be significantly correlated with several clinical features of glioblastoma.¹³ Moreover, anoikis-related signatures can be used to stratify risks in patients with endometrial carcinoma and predict their survival outcomes.¹⁴ Anoikis-related genes were also found to play a carcinogenic role in patients with clear cell renal cell carcinoma, with a model integrating multiple genes related to this disease found to predict the risk of anoikis.⁸ A seven-gene signature and nomogram was found to help select personalized treatment for patients with head and neck squamous cell carcinoma (HNSCC),¹⁵ and decision curve analysis (DCA), was used to develop a nomogram for selecting treatment strategies in patients with low-grade glioma based on their clinicopathological features.⁹

Upregulation of the TGFb/SMAD pathway in non-small cell lung cancer (NSCLC) was found to correlate with the role of PLK1,¹⁶ with SHC1 found to play a carcinogenic role in lung cancer.¹⁷ ANGPTL4 was found to affect the energy metabolism of NSCLC cells,¹⁸ with CtSLCO1B3 playing an important role in the development and/or progression of NSCLC.¹⁹ The ability of SPIB to provide prognostic information in pan-cancer has been assessed based on mutational load, microsatellite instability, immune cell infiltration, and prognosis.²⁰ Klotho (KL) has been shown to regulate autophagy,²¹ with several steps of the autophagic process involving DAPK2.²²

Higher ITGB4 expression has been associated with tumor aggressiveness and metastasis, including in patients with prostate cancer,²³ NSCLC,²⁴ hepatocellular carcinoma,²⁵ cervical cancer,²⁶ HNSCC²⁷ and colon cancer.²⁸ In addition, high ITGB4 level has been associated with a poorer prognosis in patients with LUAD.^24,29 ITGB4 is a signaling protein that promotes intracellular signaling and protects the integrity of epithelial cells.³⁰ Moreover, bioinformatics analysis showed that ITGB4 acts as a pan-cancer oncogene in 33 different types of human tumors.³¹ ITGB4 overexpression has been associated with venous invasion and decreased overall survival in patients with NSCLC.³² ITGB4 has also been closely associated with other respiratory diseases, including spontaneous inflammation of the lungs, airway inflammation and hyperreactivity, and acute lung injury.^33–35 Although anoikis and ITGB4 have been shown to be beneficial in patients with LUAD, these patients did not benefit from treatment with ITGB4 suppressing agents, indicating the need for greater understanding of the microenvironment of these tumors.

The results of this study showed that ITGB4-related genes are enriched in pathways related to energy metabolism (oxidative phosphorylation), DNA replication (E2F targets) and proto-oncogene function (Myc targets). Genes in the E2F family that regulate c-Myc expression also regulate apoptosis and cell proliferation and can serve as both oncogenes and tumor suppressor genes. No clear correlation has been observed between changes in E2F expression and changes in ITGB4 expression. Interestingly, some genes associated with oxidative phosphorylation were upregulated. EMT has been found to correlate with a switch from oxidative phosphorylation to glycolysis,³⁶ although mitochondrial respiration was also reported to be directly linked to the growth of tumor metastases.³⁷ These findings suggest that ITGB4 may undergo signaling via the E2F pathway, which influences the oxidative phosphorylation process through the MYC pathway. However, there are still some limitations in this study. Moreover, vivo and vitro studies should be performed to verify related regulatory axis. In the next step, we will experimentally verify the significance of anoikis in LUAD.

Conclusion

Our anoikis-related signature from RNA-seq data may be a novel prognostic biomarker in patients with LUAD. It may help physicians develop personalized LUAD treatments in clinical practice. Moreover, ITGB4 may affect the development of LUAD through the oxidative phosphorylation pathway.