Yiyi Zhou
Zhenli Hu- Department of Respiratory and Critical Care Medicine, Changhai Hospital, Shanghai, China
Background: Anoikis has therapeutic potential against different malignancies including lung adenocarcinoma. This study used anoikis and bioinformatics to construct a prognostic model for lung adenocarcinoma and explore new therapeutic strategies.
Methods: Several bioinformatic algorithms (co-expression analysis, univariate Cox analysis, multivariate Cox analysis, and cross-validation) were used to screen anoikis-related genes (ARGs) to construct a risk model. Lung adenocarcinoma patients were divided into training and testing groups at a ratio of 1:1. The prognostic model was validated by risk score comparison between high- and low-risk groups using receiver operating characteristic curve (ROC), nomograms, independent prognostic analysis and principal component analysis. In addition, two anoikis-related genes patterns were classified utilizing consensus clustering method and were compared with each other in survival time, immune microenvironment, and regulation in pathway. Single cell sequencing was applied to analyze anoikis-related genes constructed the model.
Results: This study demonstrated the feasibility of the model based on seven anoikis-related genes, as well as identifying axitinib, nibtinib and sorafenib as potential therapeutic strategies for LUAD. Risk score based on this model had could be used as an independent prognostic factor for lung adenocarcinoma (HR > 1; p < 0.001) and had the highest accuracy to predict survival compared with the clinical characteristics. Single cell sequencing analysis discovered Keratin 14 (KRT14, one of the seven anoikis-related genes) was mainly expressed in malignant cells in various cancers.
Conclusion: We identified seven anoikis-related genes and constructed an accurate risk model based on bioinformatics analysis that can be used for prognostic prediction and for the design of therapeutic strategies in clinical practice.
Introduction
Lung cancer is one of the malignancies with the highest morbidity and mortality (Siegel et al., 2021). In 2022, an estimated 236,740 new cases and an estimated 130,180 deaths are covered by lung carcinoma. On account of the hysteretic nature of datum, the 5-year relative survival is accessible lags 4 years behind the current year, with 22.9% from 2012 to 2018 (Cancer Stat Facts: lung and bronchus cancer, 2022). The most common subtype of lung cancer including never-smokers is lung adenocarcinoma (LUAD) (Lee et al., 2019; Devarakonda et al., 2021). Research advances have improved the early diagnosis of patients, thus leading to a better prognosis. Several treatment options are available including immunotherapy, targeted therapy, chemotherapy, radiotherapy, and surgical management (pneumonectomy and lung transplantation) (Glanville and Wilson, 2018; Wang et al., 2022a; Gross et al., 2022; Huo et al., 2022; Nguyen et al., 2022). However, LUAD is characterized by considerable heterogeneity in clinical features, pathological characteristics, and molecular alterations, as well as genomic instability (Geng et al., 2021; Okudela et al., 2021). The prevalence of lung tumors in non-smokers has increased recently, despite the known association between smoking and lung cancer (Jeon et al., 2018), underscoring the need to elucidate the underlying molecular mechanisms and develop effective therapies.
Recently, programmed cell death has been the focus of research aimed at developing new cancer treatment strategies. There are several cell death mechanisms including apoptosis (Green, 2022), autophagy (Russell and Guan, 2022), cuproptosis (Tsvetkov et al., 2022), ferroptosis, pyroptosis and necroptosis (Tang et al., 2020). A key factor in maintaining cell homeostasis is adhesion of cells to matrix, and disruption of this interaction can adversely affect cell survival (Chiarugi and Giannoni, 2008). Disengagement between cells and matrix can lead to long-time cells suspension and cell death. This cell death form has been termed “anoikis” (Frisch and Francis, 1994). The presence of anoikis may hold back the reattachment of detached cells and malformation of cells development. Anoikis after the loss of cell anchorage is associated with cells development, tissue homeostasis and malignancies (Danial and Korsmeyer, 2004; Chiarugi and Giannoni, 2008; Simpson et al., 2008). After escaping from the adhesion of extracellular matrix and intercellular contact, tumor cells survived by paracrine and paracrine mechanisms against apoptosis, and regained the ability to attach to spread, metastasize and invade. Anti-anoikis is an important feature of tumor metastasis, which enables tumor cells to spread to distant organs through the circulatory system (Simpson et al., 2008; Chaffer et al., 2016). The regulation of anoikis in lung carcinoma mainly includes extracellular matrix (ECM) and cell adhesion, cell detachment and directional migration, signal transduction and regulators (Wang et al., 2022b). Anchor independent survival of tumor cells requires detachment from ECM matrix (Pickup et al., 2014; Enkhbat et al., 2022). The inhibition of the overexpression of fibronectin in cell aggregation during detachment has been proved that can enhance anoikis in cancer (Han et al., 2021). In the drug resistance experiments of lung adenocarcinoma, upregulated Laminin 5 can reduce anoikis through activating integrin/focal adhesion kinase (FAK) signaling (Kodama et al., 2005). Generally, cells following the loss of adhesion won’t express epidermal growth factor receptor (EGFR) and then induce apoptosis. Research indicated that family with sequence similarity 188 member B (FAM188B) can prevent the degradation of EGFR and lead to the re-adhesion of lung cancer cells to the ECM (Jang et al., 2021). Non-coding RNA and IL-13 receptor subunit alpha-2 (IL13Rα2) which both involved in PI3K/AKT pathway regulation inhibit anoikis and result in lung carcinoma metastasis and grave prognosis (Xie et al., 2015; Tian et al., 2020). Anoikis play an important role in LUAD; however, the potential role of ARGs in LUAD still needs to be examined to date. Here, we hypothesized that ARGs play a role in the occurrence and development of LUAD, and constructed models to explore the mechanism of ARGs in LUAD.
The study identified ARGs and constructed a prognostic model based on the ARGs signature, which revealed the association of these markers with prognostic prediction and treatment modalities in LUAD, as well as their role in immune microenvironment, and drug sensitivity. Additionally, single cell data were used to analyze these genes for constructing the model. The functions of ARGs in the classification of LUAD were evaluated based on the TCGA and GEO datasets.
Materials and methods
Data acquisition
Transcriptome profiling (RNAseq) data and relevant clinicopathological information in LUAD, were obtained from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) and Gene Expression Synthesis Database (GSE41271, GEO, https://www.ncbi.nlm.nih.gov/geo). ARGs were incorporated through a systematic search and processing from the Human Gene Database (GeneCards, https://www.genecards.org/) and Harmonizome platform (https://maayanlab.cloud/Harmonizome/). The scRNA-seq database from Tumor Immune Single-cell Hub 2 (TISCH2, http://tisch.comp-genomics.org/) website was analyzed to explore the model-related genes expression in pan-cancer.
Construction of a prognostic network and CNV analysis
Univariate Cox regression analysis was performed using R software to screen differentially expressed genes (DEGs). We used “igraph,” “psych,” “reshape2” and “RColorBrewer” R packages to examine the correlation between ARGs and risk factors. After downloading copy number variation (CNV) data from UCSC Xena database (http://xena.ucsc.edu/), we processed the CNV data to calculate frequency of mutation using the Perl script algorithm. The “RCircos” package was used for graphical representation.
Consensus clustering analysis and principal component analysis of ARGs
According to the consensus level of ARGs, the cohort from TCGA and GEO was allocated into several clusters utilizing the R package “ConsensusClusterPlus”. The fitness of the classification was judged based on the principal component analysis (PCA) algorithms. “survival” and “survminer” were also applied to construct curve plotting and compare the survival between different clusters. The other graphical output consisted of heatmap and boxplots using “pheatmap” and “boxplots”.
Pathway enrichment analysis of ARGs
Gene set variation analysis (GSVA), an extension of gene set enrichment (GSE) method, are performed for unsupervised classification of a sample based on the variation of pathway activity (Hanzelmann et al., 2013). Gene set enrichment analysis (GSEA) is an analytical method and focused on gene sets which share common biological function, enrichment pathways and chromosomal location (Subramanian et al., 2005; Cao et al., 2019). The Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed to present the discrepancy of pathways between those ARGs clusters. The process was completed by the “GSEABase,” “GSVA” and “clusterProfiler” R package with the screening criteria adjusted p < 0.05.
Establishment and validation of an anoikis-related risk model
To identify prognostic ARGs, univariate Cox regression analysis was performed using R software. LUAD patients were randomly divided into training and testing cohorts at a ratio of 1:1. The “glmnet” R package was used to perform the least absolute shrinkage and selection operator (LASSO) regression to analyze prognostic candidates. With the optimization criterion of least error, cross validation was performed to screen target ARGs, which were applied to the optimal multivariate Cox model. Subsequently, an ARGs prognostic signature was constructed utilizing the “survival” package based on the coefficient and hazard ratio (HR). The risk score formula was as follows:
Where XARGi represents the expression of each anoikis-related gene and the corresponding CoefARGi represents the regression coefficient. The median risk score threshold was used to distinguish the high-risk group from the low-risk group. The time-dependent ROC curves and the area under the curve (AUC) were established to evaluate the risk model. Multivariate Cox regression analyses were used to assess the correlation between clinical features and prognosis, and to determine whether the features could be independent prognostic factors. A nomogram was generated to predict the survival of patients based on risk scores through the summation of each prognostic clinical variable, and calibration curves were constructed to assess the accuracy of the risk model.
Immune cell infiltrates and immune microenvironment
The relationship between risk groups and 22 distinct leukocyte subsets were identified using CIBERSORT algorithm. ESTIMATE algorithm was carried out to calculate tumor microenvironment scores and research the difference of immune and stromal scores between the high- and low-risk groups.
Analysis of sensitivity to clinical therapies
OncoPredict, an R package for predicting drug responses and biomarkers (Maeser et al., 2021), was used to analyze the sensitivity of common oncological therapeutic drugs in TCGA-LUAD dataset. Drug susceptibility documentation and risk files were subsequently processed to draw boxplots of the two risk groups.
Statistical analysis
The annotation and collation of transcriptome data, clinical and gene expression data were performed using PERL (version 5.32.0.1). Other statistical procedures were performed with R software (Version 4.1.2; Version 4.2.1). The Wilcoxon test was used to analyze anoikis-related variation according to sample quantity.
Results
Acquisition of differential ARGs and CNV analysis in LUAD
To identify ARGs in LUAD, we abstracted 598 samples from the TCGA-LUAD cohort (59 normal and 539 tumor tissues) and 183 adenocarcinomas tissues from the GSE41271 dataset. A total of 640 ARGs were identified from GeneCards and Harmonizome platform. We utilized “limma” package to locate ARGs with significant variation in TCGA cohort, and the differential expression levels were presented in heatmap and volcano map (Figure 1A, B). The heatmap showed the top 50 significant ARGs in upregulation and downregulation. The red dots and green dots represented upregulated and downregulated genes in volcano map, respectively. The detail results of differential ARGs are attached in Supplementary Table S1. After merging gene expression information from TCGA cohort and GSE41271 dataset, we abstracted the expression of 111 differential ARGs (Supplementary Table S2) and combined it with survival materials to perform survival analysis. 58 ARGs were eventually identified from the 111 ARGs using univariate cox regression analysis and were illustrated in the forest plot and corresponding network visualization. 10 of 58 ARGs, such as PIK3R1 and ITGA8, were marked in blue with an HR < 1 and meant favorable factors of prognosis. Whereas the rest genes were identified with an HR > 1 and represented poor prognosis (Figure 1C). The correlation between these ARGs and prognosis were more directly exhibited in Figure 1D. The CNV frequency of 58 ARGs that had corresponding CNV data is shown in Figure 1E. As illustrated in Figure 1E, the frequency of gain of CNV on the left was obviously higher than the loss. Chromosome 7 and chromosome 11 both had four genes displayed as a “gain” in the copy number cycle diagram shown in Figure 1F.
FIGURE 1. Acquisition of ARGs and CNV analysis. (A, B) Different expression of ARGs in LUAD with established standards (|log2FC| > 1 and FDR <0.05). (C, D) Forest plot and network graph of identified 58 ARGs correlated with prognosis via univariate Cox regression analysis. The significance criteria were set as p-value < 0.05. (E) The CNV frequency of gaining or losing in each ARGs. (F) Location and copy number of ARGs on chromosomes presented by CNV cycle diagram. Red spots mean gain, accordingly, blue one means loss.
Two ARGs patterns identified by 58 significant ARGs
When divided into two groups, the 58 significant ARGs had a significant difference using consensus clustering method (Figures 2A–D). Figures 2E, F based on cluster A and cluster B demonstrated the differential expression levels of these 58 significant ARGs. As the histogram shown, 55 ARGs had remarkable difference of expression in the two clusters, except for CDX2, SLCO1B3, TRIM31. Most of the expression of these genes were obviously higher in cluster A than in cluster B. Such as PHLDA2, CDKN3, CDKN2A, PLAU, PLK1 and so on. A PCA graph was drew for the verification of cluster allocation based on the ARGs (Figure 5G). It turned out that ARGs could classify the two patterns. Figure 5H showed that survival (p < 0.001) differed significantly between the two cluster subgroups, with a better survival rate in the cluster B than in the cluster A.
FIGURE 2. Consensus clustering of ARGs and survival analysis using K-M curves. (A–D) Consensus matrices of 58 significant ARGs (k = 2–5). (E, F) Heatmap and histogram of the 58 ARGs expression. (G) PCA for the expression profiles of different patterns. (H) Comparison of survival possibility in two subgroups. The significance criteria were set as p-value < 0.05. *p < 0.05, **p < 0.01, and ***p < 0.001.
Identification of KEGG pathway and immune analysis
Immune-related analyses were performed to compare the differences between two ARG clusters. We investigated the immune cell abundance landscape to explore the differences between subgroups using the ssGSEA algorithm. There were marked differences in activated B cell, activated CD4 T cells, CD56 bright natural killer cell, CD56 dim natural killer cell, eosinophil, gamma delta T cell, immature dendritic cell, macrophage, mast cell, monocyte, natural killer T cell, plasmacytoid dendritic cell, regulatory T cell, T follicular helper cell, type 17 T helper cell, type 2 T helper cell (Figure 3A).
FIGURE 3. Gene set enrichment analysis. (A) Expression differences of immune cell infiltration between ARG cluster A and (B) *p < 0.05, **p < 0.01, and ***p < 0.001. (B–D) The KEGG enrichment analysis for the anoikis-related differentially expressed genes uncovered potential mechanism on the occurrence and development of the two patterns. Construction of cuproptosis-related risk model. (A) The point of minimum error was picked out by vertical dotted line based on LASSO cross validation. (B) Dynamic process diagram of Lasso filtering variables.
KEGG analysis was applied for comparison of the KEGG pathways in two ARG clusters. 20 significant KEGG pathways were plotted in the heatmap (Figure 3B). 13 pathways were activated in ARG cluster A, and cluster B was on the contrary. The diagram in Figures 3C, D demonstrated that ARG cluster A and B were primarily enriched in arachidonic acid metabolism, cell cycle, DNA replication, oocyte meiosis and P53 signaling pathway. Only arachidonic acid metabolism pathway turned off (“ silenced”) in cluster A, while only it turned on (“ activated”) in cluster B.
Establishment and assessment of an anoikis-related risk model
The 58 ARGs were subjected to LASSO regression analysis and proportional-hazards model analysis with evaluation of prognosis. We identified the ARGs for the construction of the model after eliminating the effect of multicollinearity among the explanatory variables and model overfitting (Figures 4A, B). Ultimately, seven genes, consisted of angiopoietin-like 4 (ANGPTL4), integrin subunit beta 4 (ITGB4), collagen type XIII alpha 1 chain (COL13A1), KRT14, Rac Family Small GTPase 3 (RAC3), Caudal Type Homeobox 2 (CDX2), Kinesin Family Member 18A (KIF18A) were included to establish the prognostic risk model of LUAD. The risk scores of the seven screened ARGs determined by multivariate analysis are shown in Table 1. The risk score was calculated with the following formula:
FIGURE 4. Construction of anoikis-related risk model. (A) The point of minimum error was picked out by vertical dotted line based on LASSO cross validation. (B) Dynamic process diagram of Lasso filtering variables. (C, E, G, I) Survival analysis and exposure rating of the training cohort. K-M survival analysis in training cohort (C). Risk heatmap of training cohort indicated that six over-expressed ARGs in high-risk group, whereas the others over expressed in low-risk group (E). Training cohort risk curve distinguished subgroups according to the median values (G). Survival scatter diagram reflected a relationship between the risk score and survival time (I). (D, F, H, J) Survival analysis and exposure rating of the testing cohort, utilizing K-M survival curve (D), risk heatmap (F), risk curve (H), survival scatter diagram (J). The significance criteria were set as p-value < 0.05.
TABLE 1. Seven ARGs selected by multivariate cox results.
The risk model divided patients into two cohorts, a training group (n = 345) and a testing group (n = 344). Follow-up analyses were performed after comparing the clinical traits of the two cohorts without significant difference (Table 2). We then constructed the K-M survival curves to compare survival between high- and low-risk subgroups in the training and testing cohorts. There were marked differences between the risk subgroups, with lower survival rates in the high-risk group and higher survival rates in the low-risk group (Figures 4C, D). The risk scores were calculated and the median threshold of the risk score was set to 1 to distinguish between high- and low-risk groups (Figures 4G, H). Survival scatter diagrams of the two cohorts showed an inverse correlation between survival time and the risk score, indicating that patients in the high-risk group had a worse prognosis (Figure 4I, J). The risk heatmap of the training cohort indicated that six ARGs (ANGPTL4, ITGB4, KRT14, RAC3, CDX2, KIF18A) were overexpressed in the high-risk group, whereas COL13A1 was highly expressed in the low-risk group (Figure 4E). The testing cohort confirmed the expression of the six ARGs in high- and low-risk subgroups (Figure 4F).
TABLE 2. Clinical characteristics analysis of subgroups.
The risk score and partial clinical traits were included in an independent prognostic analysis using univariate Cox regression. Both risk scores (p < 0.001) in training and testing group were associated with prognosis, suggesting that the risk score could be used as an independent predictor (Figures 5A, B). As shown in Figures 5C, D nomogram and calibration curve were established to predict the survival time of patients with LUAD. The C-index curve confirmed that the risk model constructed had the highest accuracy to predict survival compared with the clinical characteristics (age, gender, and stage) (Figures 5E–G). The ROC curve showed the highest AUC value of 0.763 for the risk score, and the AUC values for 1 -, 3 -, and 5-year OS were 0.730, 0.763, and 0.741, respectively (Figure 5H). The AUC values for 1 -, 3 -, and 5-year OS were 0.666, 0.663, and 0.631 in testing group (Figure 5I).
FIGURE 5. Independent prognostic analysis and nomogram construction for LUAD. Univariate Cox regression analysis in training group (A) and testing group (B) of clinical features in the risk model set (p < 0.05). (C) The construction of nomogram to foretell patient’s 1 -, 3 -, and 5-year survival rates, according to the total points. (D) Calibration curve of the nomogram. (E–G) C-index curve for different variables, including risk score, age, gender, and stage. ROC curves for 1 -, 3 -, and 5-year overall survival in training group (H) and testing group (I).
Immune-related and drug sensitivity analyses based on the risk model
After investigating the immune cell abundance using the CIBERSORT algorithm, we compared the correlation between immune cell content and risk scores. The correlation analysis plot demonstrated that 12 immune cells had significant associations with risk scores (Figure 6). There were remarkably positive correlations in activated memory CD4 T cells, M0 macrophages, M1 macrophages, neutrophils, resting NK cells and follicular helper T cells. The others, including activated mast cells, activated NK cells, memory B cells, monocytes, resting dendritic cells, resting mast cells and resting memory CD4 T cells, had negative correlations with risk scores. In addition, we analyzed immune cells expression in the two risk groups and identified 11 significantly different immune cells (Figure 7A). A heatmap of the correlation of immune cells with the seven model ARGs was generated and presented a close connection between KIF18A and the cells (Figure 7B). The correlations between these 22 immune cells which were allocated for immune analyses were found in Figure 7C. Activated memory CD4 T cells showed high positive correlations with CD8 T cells. The low-risk group got higher immune score (p < 0.001) and stromal score (p < 0.001) than the high-risk group (Figure 7D).
FIGURE 6. Correlation between immune cells and risk scores. Positive correlations: activated memory CD4 T cells, M0 macrophages, M1 macrophages, neutrophils, resting NK cells and follicular helper T cells. Negative correlations: activated mast cells, activated NK cells, memory B cells, monocytes, resting dendritic cells, resting mast cells and resting memory CD4 T cells. p < 0.05 was considered as significant correlation.
FIGURE 7. Immune-related analysis based on the risk model. (A) Expression differences of immune cell infiltration in risk subgroups. (B) Correlation between 7 ARGs, risk scores, and the immune cells. (C) Correlation between these immune cells. (D) Other characteristics including estimate score, immune score, and stromal score between the high- and low-risk groups. (E) Differential risk score comparison in the two ARGs patterns. (F) Sankey diagram of the relationship between two ARGs patterns, survival status, and risk scores. *p < 0.05, **p < 0.01 and ***p < 0.001.
Considering the relevance of ARGs clusters to risk model, we further scored the two clusters and compared the difference of the risk score, and the risk scores in A group was higher than that in B group (Figure 7E). The sankey diagram distinctly displayed the association between clusters and anoikis-related risk scores (Figure 7F).
The correlation analysis of IC50 with risk scores showed that 73 anticancer drugs had significant correlations with risk scores, including 17 drugs with negative associations and 56 drugs with positive associations. This indicated that patients in the high-risk group responded well to most medications. Detailed information is found in Supplementary Table S3. The 3 kinds of molecular targeted drugs with significant differences (p < 0.05) are shown in Figures 8A–C, which might be candidate medications for patients with LUAD in the high-risk group. In addition, there were 2 types of frequent medications (erlotinib and gefitinib) associated with LUAD and had lower sensitivity in high-risk group (Figures 8D, E).
FIGURE 8. Drug sensitivity analysis in low-risk and high-risk sets. We listed 5 kinds of anticancer drugs, including axitinib (A), nibtinib (B), sorafenib (C), erlotinib (D) and gefitinib (E), with blue boxplots represented low-risk group and red expressed high-risk group. Erlotinib and gefitinib which were common targeted drugs to LUAD were not a better choice for patients with high-risk scores. The significance criteria were set as p-value < 0.05.
Analysis of ARGs based on single-cell RNA sequencing data
GSE131907, a single-cell RNA sequencing data of LUAD, was chosen to exhibit the expression of seven model ARGs utilizing TISCH2 platform. The landscape of intermediate cells was shown in Figure 9A. The pie diagram depicted the proportion of each cell type (Figure 9B). As it illustrated, CD4 conventional T cell was the most abundant immune cell. Following the research of proportion of each cell type in patients (Figure 9C), the expression of seven ARGs was analyzed and draw in Figure 9D–J. ANGPTL4, CDX2, ITGB4, KRT14, and RAC3 were mainly enriched in epithelial cells. COL13A1 was mainly expressed in fibroblasts and KIF18A was mainly expressed in CD8 T cells. Subsequently, we explored expression of these ARGs in pan-cancer single-cell sequencing dataset and found that only KRT14 was highly expressed in various tumors, including basal cell carcinoma (BCC), non-hodgkin lymphoma (NHL), oral squamous cell carcinoma (OSCC) and squamous cell carcinoma (SCC). In various tumors, KRT14 was widely expressed in malignant cells in the tumor microenvironment (Figure 9K).
FIGURE 9. Seven ARGs analyses based on single-cell RNA sequencing data. (A) Annotation of cell types in GSE131907. Proportion of each cell type in the dataset (B) and in each patient (C). (D–J) Expression of ANGPTL4, CDX2, COL13A1, ITGB4, KIF18A, KRT14, and RAC3. (K) Expression of KRT14 in pan-cancer.
Discussion
Despite advances in therapeutic modalities for malignant lung tumors, sundry drugs can result in a poor prognosis because of undesirable side effects, resistance, and cachexia after progression (Wang et al., 2021). Chemoresistance, particularly, remains an important cause of failure of treatment and the driving force of death. Inherent resistance leading to poor initial treatment response and acquired resistance leading to poor follow-up responses constitute the two aspects of resistance (Gottesman, 2002; Nikolaou et al., 2018; Jing et al., 2020). Chemotherapeutic drugs act on cancers by inducing programmed cell death, whereas tumor death evasion can increase the occurrence of resistance (Wattanathamsan et al., 2019; Yang et al., 2021). Several studies have focused on examining cell death mechanisms to overcome drug resistance. A previous review reported that dysregulation of ULK1-mediated autophagy plays a crucial role in drug resistance (Liu et al., 2020). Ferroptosis induces tumor cell death through lipid peroxidation mechanistically and via shrunken cell mitochondria morphologically (Dixon et al., 2012).
Our research focused on anoikis in LUAD, which is one of the programmed cell death. We hypothesized that ARGs act on lung tumor cells in a similar way as ferroptosis-related or cuproptosis-related genes identified in previous research (Tian et al., 2021; Zhang et al., 2022). 58 ARGs were identified for the establishment of anoikis-related risk model and the ARGs clustering analysis. The feasibility of the model based on seven ARGs was verified using a variety of analytical methods. And it implied that ANGPTL4, CDX2, ITGB4, KRT14, RAC3, and KIF18A were potential therapeutic targets for LUAD. Regarding drug resistance, the anoikis-related risk model analyzed hundreds of anticancer medicines to compare drug responses in the risk subgroups and provide treatment strategies for LUAD. Both training ang testing groups demonstrated that patients prognosis could be accurately predicted by anoikis-related model. Subsequently, the correlations between the immune cells and risk scores, risk subgroups and seven ARGs were taken into consideration. Patients with increasing M0 macrophages were linked with graver prognosis (Liu et al., 2017; Yi et al., 2021), and the linear relation between risk scores and M0 macrophages in our research also consisted with the previous study. Conversely, memory B cells had a negative relation with risk scores. It was reported that memory B cells was the foundation for having lasting immunity (Fan et al., 2021) and it was closely associated to superior survival (Zhang et al., 2019). We also found immune score and stromal score were both higher in low-risk group, implying that TME scores may exhibit a distinct prognostic value in LUAD (Wu et al., 2021).
We analyzed immune microenvironment, KEGG analysis and prognosis to explore two distinct ARGs patterns. Patients with poor prognosis in cluster A were notably more than cluster B. In addition, increased number of activated CD4 T cells, CD56 bright natural killer cell, CD56 dim natural killer cell, gamma delta T cell, natural killer T cell, regulatory T cell, type 2 T helper cell in cluster A was associated with poor prognosis in LUAD. KEGG analyses showed that genes were mainly involved in arachidonic acid metabolism, cell cycle, DNA replication, oocyte meiosis and P53 signaling pathway, confirming differences in immune factor abundance and immune activity between the two groups. Moreover, a high inhibition of arachidonic acid metabolism and activation of cell cycle, DNA replication, oocyte meiosis and P53 signaling pathway may lead to poor prognosis in LUAD.
We performed the single-cell analysis to seek the cell types that express the model ARGs. The expression level of KRT14 in various cell types, including immune cells, malignant cells, stromal cells, and functional cells, was explored. It depicted that KRT14 was mainly overexpressed in malignant cells and immune cells (especially epithelial cells, CD8 T cells, endothelial cells, monocyte/macrophage). Several studies have reported that high KRT14 expression can cause poor prognosis in cancers. KRT14 associated with TRAIL and TNF receptor signaling pathway led to a worse prognosis in LUSC (Dong et al., 2020) and KRT14 following the upregulation of Gkn1 could inhibit anoikis and contribute to lung nodal metastasis and poor survival (Yao et al., 2019). In breast cancer, KRT14+ epithelial tumor cell clusters promoted to distant organs metastases (Cheung et al., 2016). The underlying correlation between KRT14 and pan-cancer needs to be analyzed further.
The present results are theoretical results based on bioinformatics. Because this is a novel topic, information about the ARGs identified, including CNV annotation, mechanistic studies, and confirmation via in vivo and in vitro experiments require further study. Cell line experiments, preclinical research, and clinical trials are needed to further detect the expression and mechanism of these ARGs.
Conclusion
This study combines bioinformatics research with ARGs to construct a risk model as a prognostic index for application in clinical practice. The feasibility of the anoikis-related model will be validated through additional clinical trials. We expect that further literature will be published to uncover the essential mechanisms underlying the function of ARGs in tumor cells and to identify new anoikis-related tumor therapies.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Ethics statement
Approval of the research protocol by an Institutional Reviewer Board. N/A. Informed Consent. N/A. Registry and the Registration No. of the study/trial. N/A. Animal Studies. N/A.
Author contributions
YZ put forward the conception and contributed to the writing of this study. YZ and ZH improved the framework, provided the materials, performed, and verified the data analysis of this experiment. ZH contributed to significant modification of this manuscript. All authors participating in the article approved the submitted version.
Funding
This project was sponsored by funding from Natural Science Foundation of Shanghai (grant 81800018): TRIM37 inhibits the activation and inflammatory cytokine production of dendritic cell and the development of asthma.
Acknowledgments
From the conception, practice and writing of the article, I got a lot of help from teachers and colleagues. I would particularly like to thank my supervisors, ZH.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2023.1125782/full#supplementary-material
References
Cancer Stat Facts: lung and bronchus cancer (2022). Seercancergov statfacts. Available at: https://seercancergov/statfacts/html/lungbhtml.2022(12-01.
Cao, Y., Tang, W., and Tang, W. (2019). Immune cell infiltration characteristics and related core genes in lupus nephritis: Results from bioinformatic analysis. BMC Immunol. 20 (1), 37. doi:10.1186/s12865-019-0316-x
Chaffer, C. L., San Juan, B. P., Lim, E., and Weinberg, R. A. (2016). EMT, cell plasticity and metastasis. Cancer Metastasis Rev. 35 (4), 645–654. doi:10.1007/s10555-016-9648-7
Cheung, K. J., Padmanaban, V., Silvestri, V., Schipper, K., Cohen, J. D., Fairchild, A. N., et al. (2016). Polyclonal breast cancer metastases arise from collective dissemination of keratin 14-expressing tumor cell clusters. Proc. Natl. Acad. Sci. U. S. A. 113 (7), E854–E863. doi:10.1073/pnas.1508541113
Chiarugi, P., and Giannoni, E. (2008). Anoikis: A necessary death program for anchorage-dependent cells. Biochem. Pharmacol. 76 (11), 1352–1364. doi:10.1016/j.bcp.2008.07.023
Danial, N. N., and Korsmeyer, S. J. (2004). Cell death: Critical control points. Cell 116 (2), 205–219. doi:10.1016/s0092-8674(04)00046-7
Devarakonda, S., Li, Y., Martins Rodrigues, F., Sankararaman, S., Kadara, H., Goparaju, C., et al. (2021). Genomic profiling of lung adenocarcinoma in never-smokers. J. Clin. Oncol. 39 (33), 3747–3758. doi:10.1200/JCO.21.01691
Dixon, S. J., Lemberg, K. M., Lamprecht, M. R., Skouta, R., Zaitsev, E. M., Gleason, C. E., et al. (2012). Ferroptosis: An iron-dependent form of nonapoptotic cell death. Cell 149 (5), 1060–1072. doi:10.1016/j.cell.2012.03.042
Dong, Y., Li, S., Sun, X., Wang, Y., Lu, T., Wo, Y., et al. (2020). Desmoglein 3 and keratin 14 for distinguishing between lung adenocarcinoma and lung squamous cell carcinoma. Onco Targets Ther. 13, 11111–11124. doi:10.2147/OTT.S270398
Enkhbat, M., Zhong, B., Chang, R., Geng, J., Lu, L. S., Chen, Y. J., et al. (2022). Harnessing focal adhesions to accelerate p53 accumulation and anoikis of A549 cells using colloidal self-assembled patterns (cSAPs). ACS Appl. Bio Mater 5 (1), 322–333. doi:10.1021/acsabm.1c01109
Fan, T., Li, C., and He, J. (2021). Prognostic value of immune-related genes and comparative analysis of immune cell infiltration in lung adenocarcinoma: Sex differences. Biol. Sex. Differ. 12 (1), 64. doi:10.1186/s13293-021-00406-y
Frisch, S. M., and Francis, H. (1994). Disruption of epithelial cell-matrix interactions induces apoptosis. J. Cell Biol. 124 (4), 619–626. doi:10.1083/jcb.124.4.619
Geng, W., Lv, Z., Fan, J., Xu, J., Mao, K., Yin, Z., et al. (2021). Identification of the prognostic significance of somatic mutation-derived LncRNA signatures of genomic instability in lung adenocarcinoma. Front. Cell Dev. Biol. 9, 657667. doi:10.3389/fcell.2021.657667
Glanville, A. R., and Wilson, B. E. (2018). Lung transplantation for non-small cell lung cancer and multifocal bronchioalveolar cell carcinoma. Lancet Oncol. 19 (7), e351–e358. doi:10.1016/S1470-2045(18)30297-3
Gottesman, M. M. (2002). Mechanisms of cancer drug resistance. Annu. Rev. Med. 53, 615–627. doi:10.1146/annurev.med.53.082901.103929
Green, D. R. (2022). The future of death. Cold Spring Harb. Perspect. Biol. 14 (2), a041111. doi:10.1101/cshperspect.a041111
Gross, D. J., Chintala, N. K., Vaghjiani, R. G., Grosser, R., Tan, K. S., Li, X., et al. (2022). Tumor and tumor-associated macrophage programmed death-ligand 1 expression is associated with adjuvant chemotherapy benefit in lung adenocarcinoma. J. Thorac. Oncol. 17 (1), 89–102. doi:10.1016/j.jtho.2021.09.009
Han, H. J., Sung, J. Y., Kim, S. H., Yun, U. J., Kim, H., Jang, E. J., et al. (2021). Fibronectin regulates anoikis resistance via cell aggregate formation. Cancer Lett. 508, 59–72. doi:10.1016/j.canlet.2021.03.011
Hanzelmann, S., Castelo, R., and Guinney, J. (2013). Gsva: Gene set variation analysis for microarray and RNA-seq data. BMC Bioinforma. 14, 7. doi:10.1186/1471-2105-14-7
Huo, K. G., Notsuda, H., Fang, Z., Liu, N. F., Gebregiworgis, T., Li, Q., et al. (2022). Lung cancer driven by BRAF(G469V) mutation is targetable by EGFR kinase inhibitors. J. Thorac. Oncol. 17 (2), 277–288. doi:10.1016/j.jtho.2021.09.008
Jang, E. J., Sung, J. Y., Yoo, H. E., Jang, H., Shim, J., Oh, E. S., et al. (2021). FAM188B downregulation sensitizes lung cancer cells to anoikis via EGFR downregulation and inhibits tumor metastasis in vivo. Cancers (Basel) 13 (2), 247. doi:10.3390/cancers13020247
Jeon, J., Holford, T. R., Levy, D. T., Feuer, E. J., Cao, P., Tam, J., et al. (2018). Smoking and lung cancer mortality in the United States from 2015 to 2065: A comparative modeling approach. Ann. Intern Med. 169 (10), 684–693. doi:10.7326/M18-1250
Jing, Y., Liang, W., Liu, J., Zhang, L., Wei, J., Yang, J., et al. (2020). Autophagy-mediating microRNAs in cancer chemoresistance. Cell Biol. Toxicol. 36 (6), 517–536. doi:10.1007/s10565-020-09553-1
Kodama, K., Ishii, G., Miyamoto, S., Goya, M., Zhang, S. C., Sangai, T., et al. (2005). Laminin 5 expression protects against anoikis at aerogenous spread and lepidic growth of human lung adenocarcinoma. Int. J. Cancer 116 (6), 876–884. doi:10.1002/ijc.21136
Lee, J. J., Park, S., Park, H., Kim, S., Lee, J., Lee, J., et al. (2019). Tracing oncogene rearrangements in the mutational history of lung adenocarcinoma. Cell 177 (7), 1842–1857.e21. doi:10.1016/j.cell.2019.05.013
Liu, L., Yan, L., Liao, N., Wu, W. Q., and Shi, J. L. (2020). A review of ULK1-mediated autophagy in drug resistance of cancer. Cancers (Basel) 12 (2), 352. doi:10.3390/cancers12020352
Liu, X., Wu, S., Yang, Y., Zhao, M., Zhu, G., and Hou, Z. (2017). The prognostic landscape of tumor-infiltrating immune cell and immunomodulators in lung cancer. Biomed. Pharmacother. 95, 55–61. doi:10.1016/j.biopha.2017.08.003
Maeser, D., Gruener, R. F., and Huang, R. S. (2021). oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief. Bioinform 22 (6), bbab260. doi:10.1093/bib/bbab260
Nguyen, T. T., Lee, H. S., Burt, B. M., Wu, J., Zhang, J., Amos, C. I., et al. (2022). A lepidic gene signature predicts patient prognosis and sensitivity to immunotherapy in lung adenocarcinoma. Genome Med. 14 (1), 5. doi:10.1186/s13073-021-01010-w
Nikolaou, M., Pavlopoulou, A., Georgakilas, A. G., and Kyrodimos, E. (2018). The challenge of drug resistance in cancer treatment: A current overview. Clin. Exp. Metastasis 35 (4), 309–318. doi:10.1007/s10585-018-9903-0
Okudela, K., Matsumura, M., Arai, H., and Woo, T. (2021). The nonsmokers' and smokers' pathways in lung adenocarcinoma: Histological progression and molecular bases. Cancer Sci. 112 (9), 3411–3418. doi:10.1111/cas.15031
Pickup, M. W., Mouw, J. K., and Weaver, V. M. (2014). The extracellular matrix modulates the hallmarks of cancer. EMBO Rep. 15 (12), 1243–1253. doi:10.15252/embr.201439246
Russell, R. C., and Guan, K. L. (2022). The multifaceted role of autophagy in cancer. EMBO J. 41 (13), e110031. doi:10.15252/embj.2021110031
Siegel, R. L., Miller, K. D., Fuchs, H. E., and Jemal, A. (2021). Cancer statistics, 2021. CA Cancer J. Clin. 71 (1), 7–33. doi:10.3322/caac.21654
Simpson, C. D., Anyiwe, K., and Schimmer, A. D. (2008). Anoikis resistance and tumor metastasis. Cancer Lett. 272 (2), 177–185. doi:10.1016/j.canlet.2008.05.029
Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U. S. A. 102 (43), 15545–15550. doi:10.1073/pnas.0506580102
Tang, R., Xu, J., Zhang, B., Liu, J., Liang, C., Hua, J., et al. (2020). Ferroptosis, necroptosis, and pyroptosis in anticancer immunity. J. Hematol. Oncol. 13 (1), 110. doi:10.1186/s13045-020-00946-7
Tian, H., Lian, R., Li, Y., Liu, C., Liang, S., Li, W., et al. (2020). AKT-induced lncRNA VAL promotes EMT-independent metastasis through diminishing Trim16-dependent Vimentin degradation. Nat. Commun. 11 (1), 5127. doi:10.1038/s41467-020-18929-0
Tian, Q., Zhou, Y., Zhu, L., Gao, H., and Yang, J. (2021). Development and validation of a ferroptosis-related gene signature for overall survival prediction in lung adenocarcinoma. Front. Cell Dev. Biol. 9, 684259. doi:10.3389/fcell.2021.684259
Tsvetkov, P., Coy, S., Petrova, B., Dreishpoon, M., Verma, A., Abdusamad, M., et al. (2022). Copper induces cell death by targeting lipoylated TCA cycle proteins. Science 375 (6586), 1254–1261. doi:10.1126/science.abf0529
Wang, H., Hsia, S., Wu, T. H., Wu, C. J., Chang, J. S., and Cheng, Y. B. (2021). Anti-inflammatory azaphilones from the edible alga-derived fungus Penicillium sclerotiorum. Mar. Drugs 19 (5), 529. doi:10.3390/md19100529
Wang, J., Luo, Z., Lin, L., Sui, X., Yu, L., Xu, C., et al. (2022). Anoikis-associated lung cancer metastasis: Mechanisms and therapies. Cancers (Basel) 14 (19), 4791. doi:10.3390/cancers14194791
Wang, X. S., Bai, Y. F., Verma, V., Yu, R. L., Tian, W., Ao, R., et al. (2022). Randomized trial of first-line tyrosine kinase inhibitor with or without radiotherapy for synchronous oligometastatic EGFR-mutated NSCLC. J Natl Cancer Inst.
Wattanathamsan, O., Hayakawa, Y., and Pongrakhananon, V. (2019). Molecular mechanisms of natural compounds in cell death induction and sensitization to chemotherapeutic drugs in lung cancer. Phytother. Res. 33 (10), 2531–2547. doi:10.1002/ptr.6422
Wu, J., Li, L., Zhang, H., Zhao, Y., Zhang, H., Wu, S., et al. (2021). A risk model developed based on tumor microenvironment predicts overall survival and associates with tumor immunity of patients with lung adenocarcinoma. Oncogene 40 (26), 4413–4424. doi:10.1038/s41388-021-01853-y
Xie, M., Wu, X. J., Zhang, J. J., and He, C. S. (2015). IL-13 receptor α2 is a negative prognostic factor in human lung cancer and stimulates lung cancer growth in mice. Oncotarget 6 (32), 32902–32913. doi:10.18632/oncotarget.5361
Yang, Y., Bai, L., Liao, W., Feng, M., Zhang, M., Wu, Q., et al. (2021). The role of non-apoptotic cell death in the treatment and drug-resistance of digestive tumors. Exp. Cell Res. 405 (2), 112678. doi:10.1016/j.yexcr.2021.112678
Yao, S., Huang, H. Y., Han, X., Ye, Y., Qin, Z., Zhao, G., et al. (2019). Keratin 14-high subpopulation mediates lung cancer metastasis potentially through Gkn1 upregulation. Oncogene 38 (36), 6354–6369. doi:10.1038/s41388-019-0889-0
Yi, M., Li, A., Zhou, L., Chu, Q., Luo, S., and Wu, K. (2021). Immune signature-based risk stratification and prediction of immune checkpoint inhibitor's efficacy for lung adenocarcinoma. Cancer Immunol. Immunother. 70 (6), 1705–1719. doi:10.1007/s00262-020-02817-z
Zhang, H., Shi, Y., Yi, Q., Wang, C., Xia, Q., Zhang, Y., et al. (2022). A novel defined cuproptosis-related gene signature for predicting the prognosis of lung adenocarcinoma. Front. Genet. 13, 975185. doi:10.3389/fgene.2022.975185
Keywords: lung adenocarcinoma, anoikis, cancer prognosis, ScRNA-seq, immune microenvironment
Citation: Zhou Y and Hu Z (2023) Anoikis-related genes combined with single cell sequencing: Insights into model specification of lung adenocarcinoma and applicability for prognosis and therapy. Front. Cell Dev. Biol. 11:1125782. doi: 10.3389/fcell.2023.1125782
Received: 16 December 2022; Accepted: 31 March 2023;
Published: 24 April 2023.
Edited by:
Qian Zhao, Shanghai Jiao Tong University, ChinaReviewed by:
Zhe Liu, Tianjin Metabolic Diseases Hospital, Tianjin Medical University, ChinaJiong Deng, Binzhou Medical University, China
Copyright © 2023 Zhou and Hu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Zhenli Hu, aHV6aGVubGkxMDBAMTI2LmNvbQ==