Abstract
Background/Aim: The immunogenic cell death (ICD) pathway plays a crucial prognostic role in lung adenocarcinoma (LUAD) therapy. The cyclic GMP–AMP synthase (cGAS)-stimulator of interferon genes (STING) pathway is an upstream mechanism that drives ICD activation, but the interaction of hub genes remains unclear. The present study aimed to investigate the hub genes involved in ICD and the cGAS-STING pathway. The prognostic performance for hub genes and related Gene Ontology (GO) terms were also investigated. Materials and Methods: Gene expression data of ICD induction and cGAS-STING pathway activation samples were extracted from the Gene Expression Omnibus (GEO) database, and gene expression as well as clinical data of LUAD patients who received pharmaceutical therapy were extracted from The Cancer Genome Atlas (TCGA). Differentially expressed genes (DEGs) were identified, and protein–protein interaction (PPI) analysis was used to identify hub genes. Hazard risk (HR) scores were identified using Kaplan-Meier (K-M) and COX analyses. Gene set enrichment analysis (GSEA) was performed to identify the related GO terms, and receiver operating characteristic (ROC) analysis was used to evaluate the prognosis performance of the related gene sets. Results: A total of 22 DEGs were identified in the two GEO datasets, and six hub genes were identified by PPI analysis. Keratin 6A (KRTA6A) and fatty acid 2-hydroxylase (FA2H) were selected as the hub genes after survival analysis. GSEA and ROC analysis indicated that there was no difference between the KRTA6A and FA2H gene sets on prognosis performance. Conclusion: KRTA6A and FA2H are hub genes associated with the induction of cGAS-STING-related ICD in LUAD.
Lung adenocarcinoma (LUAD) is the dominant subtype of lung cancer, and the 5-year survival rate for lung cancer is poor (1). Although the treatment methods for LUAD have improved in the past decades, the average survival time of LUAD patients is less than 1 year (2).
Given the high clinical efficiency, immune-related biomarkers have attracted increasing attention (3). Immune cells in the tumor microenvironment (TME) are central regulators for tumor progression (3, 4). Recent studies have suggested that chemotherapy induces immunogenic cell death (ICD) in the TME (5, 6). Different damage associated molecular patterns (DAMPs) are released through ICD, and they are regarded as the clinical response signals of cancer cell apoptosis or autophagy (7). Moreover, apoptosis and autophagy cause different clinical responses and outcomes (8, 9). Therefore, it is valuable to explore specific immunomodulators to design clinical strategies that integrate drugs and apoptosis markers.
The DNA damage response (DDR) machinery is deficient in LUAD patients (10). Different radiotherapy and chemotherapy regimens induce various ICD-related pathways, and the cyclic GMP–AMP synthase (cGAS)-stimulator of interferon genes (STING) pathway is an alternative mechanism triggering cell death under the background of deficient DNA damage response (11, 12). Moreover, activation of the cGAS-STING pathway also leads to ICD via induction of interferons (IFNs) and promotion of sterile inflammatory diseases (13, 14). In addition, the STING-dependent downstream pathway may play a crucial role in shaping an immunosuppressive tumor microenvironment and contributes to immune invasion (15). Therefore, it is important to delineate the stepwise mechanism and identify the hub genes that drive the activation of the cGAS-STING pathway to induce ICD. In the present study, we investigated the hub genes involved in this mechanism by integrated analysis of GEO data from independent activation of the cGAS-STING pathway and ICD. Furthermore, the prognosis evaluation of the selected hub genes in LUAD patients who received pharmaceutical therapy was also investigated.
Materials and Methods
Data preparation. Gene expression data from eight ICD induction lung cancer cell lines (GSE126988 dataset) and six cGAS-STING activated lung cancer cell lines (GSE113519 dataset) were downloaded from the GEO database. The cGAS-STING triggered ICD related genes always plays a critical role in regulating the cells’ survival responding to pharmaceutical therapy (13-15). So, it is more efficient and precise to evaluate the prognosis role of the selected hub genes in LUAD patients who received pharmaceutical therapy. Therefore, 200 individual gene expression profiles of LUAD patients who received pharmaceutical therapy from The Cancer Genome Atlas (TCGA) were also downloaded. The corresponding clinical information of patients was also obtained from TCGA database following TCGA publication guidelines and data access policies. The sequencing platform for the GEO databases was Illumina NextSeq 500. The Mus musculus gene symbols were converted to the homologous Homo sapiens gene symbols by R software. The experiment information for the GSE126988 dataset was obtained from a previous publication (16). The above protocols were approved by The Ethics Committee of the Fourth Hospital of Hebei Medical University.
Identification of differentially expressed genes (DEGs). DEG profiles between the different groups (including ICD induction lung cancer group, cGAS-STING activated lung cancer group and TCGA LUAD patients who received pharmaceutical therapy group, respectively) were screened using R software with the following threshold: absolute value of log2fold change (|log2FC|) >1 and false discovery rate (FDR) <0.05. The DEG heatmaps were generated using the pheatmap package in R. The Venn diagram package in R was used to analyze the DEGs between the ICD induction and cGAS-STING activated datasets. In total, 22 common gene signatures were selected for further analysis.
Protein–protein interaction (PPI) network analysis. The protein interactive network of the 22 DEGs was constructed based on the STRING database. An interactive relationship >0.3 was considered as the threshold. CytoHubba was utilized to identify the top four hub genes.
Clinicopathological characteristic analysis and survival analysis. We evaluated the association between the clinicopathological characteristics and the expression of six hub genes. The overall survival (OS) rate was compared by log-rank test, and Kaplan-Meier (K-M) survival curves were generated. Univariate and multivariate COX analyses were utilized to investigate the effect of hub gene expression on OS. p<0.05 was considered statistically significant.
Gene set enrichment analysis (GSEA). To explore the role of the selected hub genes in regulating the downstream cell signaling pathways, GSEA was applied to assess the enriched pathways in TCGA LUAD gene expression. The LUAD dataset was divided into high expression and low expression subgroups according to the expression of the hub genes. The Gene Ontology (GO) terms were considered significantly enriched according to the following threshold values: FDR<0.25, p<0.05, and |normalized enrichment score (NES)|>1. The top 10 GO terms were used for further analysis.
Prognosis evaluation of GSEA gene sets. ROC analysis was performed to validate the predictive value of the merged GSEA gene sets, and the area under the curve (AUC) values were calculated. The survivalROC package in R software was used for drawing the ROC curve. The calibration curve and Hosmer-Lemeshow test were utilized to evaluate the discrimination of the prediction probabilities against the observed rates by R software.
Statistical analysis. Statistical analyses were performed via R software (version 4.0.2) and SPSS software version 20.0 (IBM Corp., Armonk, NY, USA). The Kruskal-Wallis test combined with Tukey’s honestly significant difference (HSD) test was used to compare multiple groups. The statistical correlation between the expression levels of the hub genes in TCGA samples and patient clinical data was compared by the Chi-squared test. p<0.05 was considered statistically significant.
Ethics approval and consent to participate. The present study was approved by the Ethics Committee of the Fourth Hospital of Hebei Medical University. The research was performed in accordance with the World Medical Association Declaration of Helsinki.
Results
DEG identification in the ICD induction and cGAS-STING activated GEO datasets. We performed differential analysis to determine common DEGs between the ICD induction and cGAS-STING activated GEO datasets. Additionally, a DEG profile was obtained from TCGA LUAD patients who received pharmaceutical therapy. Heatmaps were generated to show the differential gene expression profiles of the ICD induction, cGAS-STING activated, and TCGA LUAD patient databases (Figure 1A-C). All DEGs were screened out according to the following threshold values: |log2FC|>1 and FDR<0.05. In total, 22 DEGs in both gene sets were identified as common shared DEGs using Venn diagrams (Figure 1C). Additionally, the common shared DEGs overlapped with TCGA LUAD patient cohort who received pharmaceutical therapy. In total, 15 DEGs were identified as the common shared DEGs among the ICD induction GEO dataset, cGAS-STING activated GEO dataset, and TCGA patient cohort dataset (Figure 1D).
DEGs commonly expressed among the datasets. Hierarchical clustering heatmaps of upregulated and downregulated DEGs extracted from the ICD induction GEO dataset (A), cGAS-STING pathway activated GEO dataset (B), and TCGA LUAD patients who received pharmaceutical therapy (C). The selected threshold was set as |log2FC|>1.0 and adjusted p<0.05. The Venn diagrams show the 22 commonly expressed genes in the two GEO datasets (D) and the 15 commonly expressed genes among all three datasets (E). DEGs: Differentially expressed genes; FC: fold change; GEO: Gene Expression Omnibus; ICD: immunogenic cell death; LUAD: lung adenocarcinoma; TCGA: The Cancer Genome Atlas.
PPI network. To validate the potential relationship among the selected 22 common DEGs shared between the ICD induction and cGAS-STING activated GEO datasets, the interaction of these DEGs was analyzed by the STRING database (Figure 2A). According to the rank order of three terms (degrees, closeness, and betweenness), the hub genes were screened out by CytoHubba, which identified six hub genes (Figure 2B).
PPI network analysis of the common shared genes between the ICD induction and GAS-STING pathway activated GEO databases. (A) The PPI network was constructed using the STRING database. (B) Top 10 hub genes were screened out by the CytoHubba software. The node color shade reflects the degree of connectivity. GEO: Gene Expression Omnibus; ICD: immunogenic cell death; PPI: protein-protein interaction.
Correlation of survival analysis and clinicopathological characteristics with hub gene expression. The six hub genes, namely, keratin 6A (KRT6A), KRT16, EGF-containing fibulin-like extracellular matrix protein 1 (EFFEMP), lysyl oxidase (LOX), fatty acid 2-hydroxylase (FA2H), and serine palmitoyltransferase 3 (SPTL3), have been reported to be closely related with the prognosis of lung cancer (14-19). To validate the effect of the above hub genes on the OS of LUAD patients, we performed survival comparative analysis using TCGA database. The K-M survival curve revealed that lung cancer patients with high or low KRT6A or FA2H expression had a better OS rate compared to patients without these genes [Figure 3; p<0.05; adjusted HR=2.446; 95% confidence interval (CI)=1.225-4.882]. Furthermore, COX analysis indicated that high expression of KRT6A or FA2H contributed to a higher clinical stage (Table I). These results indicated that KRT6A or FA2H may have a clinical prognostic value for the response of LUAD patients to pharmaceutical therapy.
Correlation of hub gene expression with LUAD patient survival in TCGA dataset. COX analysis was performed to obtain the adjusted HR (p=0.0057; adjusted HR=2.159, 95% CI=1.0771-4.0682). CI: Confidence interval; HR: hazard risk; LUAD: lung adenocarcinoma; TCGA: The Cancer Genome Atlas.
Cox regression analysis of clinical-pathological characteristic in TCGA NSCLC dataset implicating KRT6A or FA2H.
GSEA analysis. To explore the related GO functional terms of KRT6A and FA2H, GSEA was performed using the differential genes from TCGA original dataset. In total, 318 and 4,673 GO functional terms were found for FA2H and KRT6A differential expression, respectively. In the top 10 enriched GO terms for KRT6A differential expression, the following four mitochondrial functional GO terms were enriched: mitochondrial translation, mitochondrial gene expression signaling, mitochondrial membrane organization, and mitochondrial translational elongation. In addition, anaphase-promoting complex-dependent catabolic process and DNA-dependent DNA replication terms, which are closely related to the mismatch repair (MMR) mechanism, were also included in the top 10 GO terms (Figure 4A). For the FA2H differential expression, mitochondrial gene expression and mitochondrial translation were included in the top 10 GO terms, and the other eight GO terms were related to the RNA splicer complex (Figure 4B).
GSEA results. (A) Top 10 enriched GO terms of the KRT6A differentially expressed group (mitochondrial translation, mitochondrial gene expression, mitochondrial membrane organization, mitochondrial translational elongation, anaphase-promoting complex-dependent catabolic process, DNA-dependent DNA replication, ribosome biogenesis, translational termination, rRNA metabolic process, rRNA processing). (B) Top 10 enriched GO terms of the FA2H differentially expressed group (mitochondrial gene expression, mitochondrial translation, RNA splicing via transesterification reactions, Spliceosomal complex, Spliceosomal snRNP complex, RNA splicing via transexterification reactions with bulged adenosine, U2-type splicesomal complex, mRNA splicing via spliceosome, tRNA metabolic process, RNA splicing). DNA: Deoxyribonucleic acid; GO: Gene Ontology; GSEA: gene set enrichment analysis; KRTA6A: keratin 6A; FA2H: fatty acid 2-hydroxylase; mRNA: messenger ribonucleic acid; RNA: ribonucleic acid.
Validation of the prognostic performance of the KRT6A and FA2H differential expression enriched gene sets. From a functional perspective, the above GSEA enriched GO terms share similar functions and can be further classified into three separate datasets, namely, the mitochondrial function-related genes (MFRGs), MMR mechanism-related genes (MMRGs), and RNA splicer complex-related genes (RSCRGs) for the convenience of further analysis. The selected gene expression for each patient and the survival outcome were merged into a matrix. The ROC curves for the three datasets were generated, and the AUC values for the 5-year OS rate predictions were compared. The AUC values for the MFRGs, MMRGs, and RSCRGs were 0.814, 0.737, and 0.521, respectively (Figure 5A). The C-index was also calculated and compared for evaluating the OS rate predictions for the three datasets, and calibration curves were generated for the above three results. The C-index values for the three MFRG, MMRG, and RSCRG fitting results were 0.7281977, 0.7181554, and 0.7042812, respectively. The Hosmer-Lemeshow test was performed to assess the difference in predicted values among the RSCRG, MMRG, and RSCRG datasets, which resulted in p-values of 0.5908, 0.5908, and 0.5888 (All p>0.05; Figure 5B-D). These results indicated that there is no statistically significant difference between the true and predicted values in the three prediction gene sets.
Validation of prognostic performance of the KRT6A or FA2H differential expression enriched gene set. (A) AUC analysis results for the MFRPGs, MMR, MMRGs, and RSCRGs. Calibration curves for the MFRPGs (B), MMRGs (C), and RSCRGs (D). AUC: Area under the curve; MFRPGs: mitochondrial function-related genes; MMR: mismatch repair; MMRGs: mechanism-related genes; RNA: ribonucleic acid.
Discussion
The resistance to apoptosis and ICD is the major obstacle for improving the efficiency of LUAD pharmaceutical therapy (17). Various chemotherapeutic drugs, oncolytic viruses, physicochemical therapies, photodynamic therapy, and radiotherapy can cause ICD (18). Dying cells release DAMPs, which are captured by the immune system, resulting in tumor-specific immune responses (19, 20). Because different ICD inducers induce multiple DAMPs, various downstream pathways are activated to respond to the therapy method (6, 21).
The cGAS-STING pathway is the double-stranded DNA (dsDNA) repair sensor, and it plays the central role in regulating the innate immune system to respond to infections, inflammation, and cancer (22). The cGAS-STING pathway is the upstream activator of ICD via type I IFN (23). Multiple myeloma patients with low STING expression do not efficiently respond to ICD induction after bortezomib (BTZ) treatment (24). Several previous studies have indicated that STING pathway activation is insufficient for ICD induction (25, 26). Furthermore, chronic activation of the cGAS-STING pathway has also been reported to play a protumor role in promoting metastasis (27). Therefore, it is necessary to clarify the detailed or possible alternative downstream pathways that communicate with the cGAS-STING pathway and ICD induction.
To explore the potential hub genes, we compared the DEGs between the cGAS-STING activated and ICD induction gene sets in lung cancer cells. In total, 22 potential candidate genes were screened out, and these genes were then used for PPI network analysis and K-M survival analysis of TCGA LUAD patient dataset. Finally, KRT6A and FA2H were selected as the prognosis-related hub genes. KRT6A and FA2H have been reported to be involved in malignant behaviors but not related to ICD (28, 29). GSEA identified that KRT6A or FA2H differential expression was associated with MFRGs, MMRGs, and RSCRGs, and these three gene sets presented a similar prognosis capacity in TCGA patient cohort. Therefore, these findings suggested that KRT6A and FA2H are related to three functional gene sets that are closely related with the prognosis of LUAD patients. Although the three gene sets have been reported to be related to ICD, the present study was unable to determine if the two hub genes participate in ICD activation due to the lack of a detailed mechanism. In a recent study, the c-Met protein, which was demonstrated to regulate cancer cell movement, has been identified as the crucial target of gigantol treatment in lung cancer cells (30). For the perspective mechanism, the KRT6A, FA2H and c-Met were all related with cell movement. Some of the downstream pathways were also linked with the three proteins. For example, mTOR and Bcl-2 pathways (downstream of c-Met) were all related to mitochondrial function in cells and enriched in our study. Therefore, it is reasonable to speculate that cell movement may be the crucial factor to re-shape the tumor microenvironment. More detailed mechanism still need to be explored. Additionally, the independent prognostic roles of the two hub genes need further examination in a clinical cohort receiving different pharmaceutical therapies.
Conclusion
Through comprehensive bioinformatics analyses, the present study showed that the KRT6A and FA2H hub genes are associated with cGAS-STING pathway-related ICD in LUAD patients. Further research on the functional pathways related to KRT6A and FA2H in TCGA lung cancer datasets is required to validate the prognostic role of these two genes.
Acknowledgements
We would like to acknowledge the support of the Fourth Hospital of Hebei Medical University, the Second Hospital of Chifeng, the Traditional Chinese Medicine Hospital of Hengshui, the Ci County People’s Hospital and the the People’s Hospital of Tang County.
Footnotes
Conflicts of Interest
The Authors declare that they have no competing interests.
Authors’ Contributions
All Authors made substantial contributions to the study conception and design and the acquisition, analysis, and interpretation of data; Xu He drafted the article and revised it critically for important intellectual content; Yongjun Yu and Xiaoyang Zhang conducted the experiments and statistical analysis. Yongjun Yu and Chuan Gao mainly interpreted the data. Hongyan Wang revised the manuscript. Hongyan Wang and Chuang Liu supplied the funding to the study. All authors agreed to submit the final manuscript to the current journal; Hongyan Wang and Chuang Liu gave final approval of the version to be published; Hongyan Wang and Chuang Liu agree to be accountable for all aspects of the work.
- Received May 24, 2023.
- Revision received July 4, 2023.
- Accepted July 19, 2023.
- Copyright © 2023, International Institute of Anticancer Research (Dr. George J. Delinasios), All rights reserved
This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY-NC-ND) 4.0 international license (https://creativecommons.org/licenses/by-nc-nd/4.0).
