Research Paper Volume 15, Issue 19 pp 10305—10329

The integrated single-cell analysis developed an immunogenic cell death signature to predict lung adenocarcinoma prognosis and immunotherapy

Pengpeng Zhang1,2, *, , Haotian Zhang1, *, , Junjie Tang1, *, , Qianhe Ren1, , Jieying Zhang3, , Hao Chi4, , Jingwen Xiong5, , Xiangjin Gong5, , Wei Wang1, , Haoran Lin1, , Jun Li1, , Chenjun Huang1, ,

  • 1 Department of Thoracic Surgery, The First Affiliated Hospital of Nanjing Medical University, Nanjing, China
  • 2 Department of Lung Cancer Surgery, Tianjin Medical University Cancer Institute and Hospital, Tianjin, China
  • 3 First Teaching Hospital of Tianjin University of Traditional Chinese Medicine, Tianjin, China
  • 4 Clinical Medical College, Southwest Medical University, Luzhou, China
  • 5 Department of Sports Rehabilitation, Southwest Medical University, Luzhou, China
* Equal contribution and share the first authorship

Received: May 30, 2023       Accepted: September 6, 2023       Published: October 4, 2023      

https://doi.org/10.18632/aging.205077
How to Cite

Copyright: © 2023 Zhang et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Background: Research on immunogenic cell death (ICD) in lung adenocarcinoma (LUAD) has been relatively limited. This study aims to create ICD-related signatures for accurate survival prognosis prediction in LUAD patients, addressing the challenge of lacking reliable early prognostic indicators for this type of cancer.

Methods: Using single-cell RNA sequencing (scRNA-seq) analysis, ICD activity in cells was calculated by AUCell algorithm, divided into high- and low-ICD groups according to median values, and key ICD regulatory genes were identified through differential analysis, and these genes were integrated into TCGA data to construct prognostic signatures using LASSO and COX regression analysis, and multi-dimensional analysis of ICD-related signatures in terms of prognosis, immunotherapy, tumor microenvironment (TME), and mutational landscape.

Results: The constructed signature reveals a pronounced disparity in prognosis between the high- and low-risk groups of LUAD patients. The statistical discrepancies in survival times among LUAD patients from both the TCGA and GEO databases further corroborate this observation. Additionally, heightened levels of immune cell infiltration expression are evidenced in the low-risk group, suggesting a potential benefit from immunotherapeutic interventions for these patients. The expression levels of pivotal risk-associated genes in tissue samples were assessed utilizing qRT-PCR, thereby unveiling PITX3 as a plausible therapeutic target in the context of LUAD.

Conclusions: Our constructed ICD-related signatures provide help in predicting the prognosis and immunotherapy of LUAD patients, and to some extent guide the clinical treatment of LUAD patients.

Introduction

Lung carcinoma (LC) demonstrates a substantial global incidence and remains the predominant contributor to cancer-associated mortalities [1]. From a pathological perspective, LC can be broadly classified into two categories: small cell lung carcinoma (SCLC) and non-small cell lung carcinoma (NSCLC), with NSCLC being the more prevalent type [2]. Within NSCLC, adenocarcinoma (ADC) stands out as the predominant histological subtype, accounting for approximately 40% of all instances of LC [3]. Over the past decade, significant attention has been directed towards early surgical intervention for LC, leading to improved survival prospects for individuals diagnosed at an early disease stage [4]. However, despite the advancements in understanding LC, the aggressive nature of lung adenocarcinoma (LUAD) and its associated risks continue to adversely impact patient survival and prognostic outlooks. In recent years, there have been notable strides in utilizing targeted therapies based on biomarkers, resulting in certain improvements in the survival outcomes of LC patients. Nevertheless, the broader adoption of such treatment modalities requires further comprehensive investigation [5]. As a result, delving into additional insights regarding LUAD and the identification of emerging factors linked to patient survival prognosis assume crucial significance in enhancing patient outcomes.

The utilization of single-cell RNA sequencing (scRNA-seq) technology facilitates the examination of gene expression patterns within individual cells and the elucidation of intercellular signaling networks [612]. This technological innovation enables the direct quantification of the transcriptional output of neoplastic cells, thereby reflecting the integration of genetic, epigenetic, and environmental cues that modulate the behavior, adaptability, and response of cancer cells to therapeutic interventions [1319]. The integration of clinicopathological data with scRNA-seq information derived from tumor samples holds the potential to unveil innovative diagnostic and prognostic biomarkers, as well as cell types that may prove therapeutically relevant [2022].

Under normal circumstances, millions of cells within the bodies of healthy adult individuals undergo programmed cell death mechanisms without inciting localized or systemic inflammation [23]. However, when certain cells exhibit adequate antigenicity, such as infected or malignant cells, their demise can engender an adaptive immune response mediated by cytotoxic T lymphocytes, thereby instigating immune memory. Immunogenic cell death (ICD) assumes a pivotal role in immune surveillance; nevertheless, both pathogens and cancer cells have developed strategies to elude ICD recognition [24]. ICD is typified by the exposure and release of a multitude of damage-associated molecular patterns (DAMPs), which significantly aid dying cancer cells by fostering the recruitment and activation of antigen-presenting cells [25]. Therapeutically induced ICD can stimulate anti-cancer immune responses, thereby enhancing the therapeutic efficacy of conventional chemotherapy and radiotherapy.

This drives the necessity to investigate signatures associated with ICD in LUAD. Unraveling these signatures holds the promise of prognostic insights and informs the potential for immunotherapeutic interventions. The identification of such signatures enables the prediction of patient prognosis and facilitates the identification of individuals who might benefit from immune-based therapies. Consequently, exploring the interplay between ICD, prognosis, and immunotherapy stands to revolutionize our approach to managing LUAD and improving patient outcomes.

Materials and Methods

Data search

The scRNA-seq dataset with the accession code GSE150938 was procured from the Gene Expression Omnibus (GEO) repository, encompassing a cohort of 12 samples originating from patients diagnosed with LUAD. Additionally, microarray sequencing data and clinical information were amassed from various datasets, including GSE13213 (Samples number=119), GSE26939 (Samples number=115), GSE29016 (Samples number=39), GSE30219 (Samples number=86), and GSE31210 (Samples number=227). Immunogenic cell death-related genes are summarized in the Supplementary Table 1. To ensure harmonization across different datasets and mitigate potential batch effects, a transformation to transcripts per million (TPM) format was performed on the expression data, followed by the application of the “combat” function from the “sva” package [26]. Moreover, bulk RNA-seq data, mutation profiles, and comprehensive clinical particulars of patients afflicted with lung adenocarcinoma were sourced from the TCGA database. A uniform log2 transformation was systematically applied to all datasets before embarking on the subsequent analytical procedures.

scRNA-seq data processing

The scRNA-seq dataset underwent processing using the “Seurat” R program [27, 28]. Cells with a count below 3 and a features count below 200 were excluded. Subsequently, cells with an nFeature RNA count below 200 and mitochondrial count exceeding 10% of the sequencing count were eliminated, resulting in a total of 45,986 cells retained for further analysis. The dataset was normalized and scaled using Seurat's NormalizeData and ScaleData functions, respectively, while mitigating batch effect through CCA. The “FindVariableFeatures” algorithm was employed to identify the top 2000 variable genes. Employing Seurat’s Uniform Manifold Approximation and Projection (UMAP) technique in conjunction with the “FindClusters” function, clustering analysis was carried out with a specified resolution of 0.8. Identification of marker genes was executed utilizing Seurat’s “FindAllMarkers” function, where selection of cluster-specific markers hinged on a log2 fold change (log2 FC) threshold set at 0.25. Commonly established cell markers were harnessed for the purpose of cell classification. For quantifying gene set activity within individual cells, the “AUCell” R package was employed, computing gene expression rankings grounded on the area under the curve (AUC) value derived from model genes [29]. Cells exhibiting heightened expression levels for the targeted gene set were associated with correspondingly elevated AUC values. Determination of the threshold for discerning cells featuring active gene sets was accomplished through the utilization of the “AUCell exploreThresholds” function. Furthermore, the UMAP embedding of each cell’s AUC score was visually represented using the “ggplot2” R package, facilitating the identification of clusters characterized by active gene set profiles.

Signature construction

The FindAllMarkers function was used to identify differential genes between high-ICD-active and low-ICD-active cell populations. These genes were subsequently integrated into transcriptome sequencing to construct a model. A total of 68 genes associated with prognosis were discerned via univariate regression analysis. To further refine the gene selection for model construction, lasso regression and multifactorial analysis were implemented [30, 31]. Risk scores were assigned to each LUAD patient utilizing the coefficients derived from the multivariate analysis. Based on the median risk score, patients from the TCGA-LUAD dataset were classified into high- and low-risk groups. Prognostic implications were evaluated by generating survival curves employing the Kaplan-Meier technique, and statistical significance was determined through log-rank tests. The evaluation of the model’s predictive efficacy was conducted through the utilization of receiver operating characteristic (ROC) curves, where a threshold AUC value exceeding 0.65 denoted exceptional performance. The validation of the predictive potential of the identified signature encompassed a comprehensive analysis across five distinct datasets from the GEO repository, involving survival analysis and AUC assessments. Furthermore, to provide a visual depiction of the patient distribution across varying risk strata, both principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) techniques were employed. The application of analogous methodologies was extended to the validation cohorts, yielding analogous insights into the distribution patterns of patients within those contexts.

The establishment of the nomogram and evaluation

By integrating the risk score with clinical features, a more precise and refined nomogram was formulated using the ‘rms’ R package [32], resulting in the augmentation of prognostic predictive capability. The effectiveness of the nomogram was evaluated through the calculation of c-index and analysis of ROC curves. Stratified evaluations were undertaken, segregating analyses according to factors such as age, pathological T status, N status, and clinical stage. This approach served to assess the predictive relevance of both the risk score derived from the signature and the inherent clinical characteristics.

Enrichment pathway exploration

For the GSVA analysis, integration involved the utilization of the “h.all.v7.5.1.symbols.gmt” gene sets sourced from MSigDB, accessible at the link (https://www.gseamsigdb.org/gsea/msigdb/index.jsp) [33]. Following this, the examination of gene set activities within each sample was carried out employing the GSEABase software package. Subsequently, GSEA methodology was implemented to uncover the signaling pathways and biological processes that exhibited significant enrichment within the high- and low-risk groups [34]. To gauge the enrichment scores associated with infiltrating immune cells and immunological processes, the ssGSEA method was deployed [35]. In the context of performing Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment analyses for the differentially expressed genes across the three distinct risk groups, the R packages “clusterProfiler” and “org.Hs.eg.db” were utilized [36]. Visualization of the outcomes derived from the enrichment analyses was effectively accomplished using the “ggplot2” and “GseaVis” R software packages.

Mutations in different risk groups

The “maftools” R package was employed to analyze somatic mutations within the high- and low-risk cohorts of LUAD patients [37]. Utilizing data from the TCGA database, we generated mutation annotation format (MAF) files to capture the mutational landscape. The evaluation of tumor mutation burden (TMB) was conducted for each individual in the LUAD cohort. Visual representation of the mutation landscape alongside immune infiltration scores was achieved through the application of the “ComplexHeatmap” R package [38]. Categorization of TCGA-LUAD patients into four distinct groups was carried out, followed by the comparison of their survival disparities based on the median risk score and TMB levels.

Immunotherapeutic strategies and assessment of the immune microenvironment

Patient data derived from the TIMER 2.0 database were procured for lung adenocarcinoma (LUAD), and the evaluation of immune cell infiltration was undertaken through the utilization of seven distinct methodologies within the TCGA dataset. Heatmaps were harnessed as a visual tool to illustrate variations in immune cell infiltration across diverse risk groups. Disparities and correlations among immune checkpoint genes specific to the high- and low-risk cohorts were systematically explored, employing boxplots and scatter plots. The quantification of enrichment scores encompassing a diverse spectrum of 29 immune signatures was accomplished through the implementation of the ssGSEA methodology. In parallel, the assessment of immunological scores, Stromal Scores, and ESTIMATE Scores for the LUAD patient cohort was executed via the utilization of the “estimate” R package [39]. Leveraging the potential to predict the responsiveness to immunotherapeutic interventions, the Cancer Immunome Atlas (TCIA) database was harnessed, serving as a repository to access the Immunophenoscores (IPS) data specific to the TCGA-LUAD population [40]. In the course of this comprehensive exploration, a rigorous juxtaposition of IPS profiles was meticulously conducted, enabling a discerning comparison between the high-risk and low-risk subgroups. This analytical endeavor elucidates valuable insights into the nuanced immunotherapeutic propensities residing within these distinct risk categories.

qRT-PCR assay

The process of extracting total RNA from LUAD tissues was successfully executed through the utilization of TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA). In strict adherence to the stipulated protocols provided by the manufacturer, the cDNA synthesis step was performed, facilitated by the RevertAid™ First Strand cDNA Synthesis Kit (Thermo Fisher Scientific). Subsequently, the application of a qRT-PCR assay transpired, employing the StepOne Real-Time PCR system (Thermo Fisher Scientific) in conjunction with a SYBR Green PCR kit (Takara Bio, Otsu, Japan). In the pursuit of quantifying the levels of relative gene expression, the 2-ΔΔCT method was aptly engaged. This systematic procedure underscores the meticulous attention to methodological precision, ensuring the robustness and reliability of the ensuing molecular analyses.

Statistical analysis

R version 4.1.3 was the platform employed for all statistical analyses and data processing endeavors. Survival analysis encompassed the utilization of Kaplan-Meier curves, with the determination of statistical significance achieved through application of the log-rank test. The generation of survival curves was facilitated by the utilization of the “survminer” R package [41]. Furthermore, the creation of heatmaps was executed using the “Pheatmap” R package. For variables exhibiting a normal distribution, quantitative disparities were evaluated using either a two-tailed t-test or a one-way analysis of variance. Non-normally distributed data was subjected to analysis through either the Wilcoxon test or the Kruskal-Wallis test. The entire spectrum of statistical analyses was undertaken within the R environment, with statistical significance deemed present at a threshold of P < 0.05.

Availability of supporting data

The datasets analyzed in the current study are available in the TCGA repository (http://cancergenome.nih.gov/), and GEO (https://www.ncbi.nlm.nih.gov/geo/).

Consent for publication

All authors consent to the publication of this study.

Results

Cells clustering results and AUC value calculation

The study’s schematic representation was visually outlined in Figure 1. Depiction of gene expression level distributions, sequencing depth characteristics, proportions of red blood cell genes, mitochondrial gene proportions, and ribosome gene proportions across the 12 samples was presented in Figure 2A. The visualization of correlations between sequencing depth and gene expression levels, along with the proportions of mitochondrial genes, red blood cell genes, and ribosome genes, was rendered in Figure 2B. Following meticulous data processing and rigorous filtration, gene expression profiles derived from a collective total of 45,986 cells spanning across 12 LUAD samples were obtained for subsequent comprehensive analysis. The consistent distribution of cells within each individual sample indicated a notable absence of pronounced batch effects, a beneficial facet for forthcoming investigations (Supplementary Figure 1A). The PCA reduction plot, depicted in Figure 2C, showed that there were no apparent distinctions observed among the various cell cycles. Subsequently, employing UMAP as the dimensionality reduction method, all cells were classified into 13 clusters (Supplementary Figure 1B). The bubble and UMAP plots in Figure 2D displayed the typical marker genes for different cell types and their corresponding clusters. Figure 2E presented the proportion of each cell type within each of the 12 samples, while Figure 2F showed the distribution of each cell type. Figure 2G illustrated the differential genes between the different cell types, with the top 4 markers randomly labeled on the left side, and the right side utilizing GO enrichment analysis to demonstrate the pathways in which the marker genes were primarily enriched. Finally, the ICD activity was calculated as a gene set for each cell, and Figure 2H displayed the variation in ICD activity among different cell types, indicating relatively higher ICD activity in NK/T cells and Myeloid cell types.

The flowchart of this investigation. To integrate multiple datasets to explore the role of ICD in the TME of LUAD, and to construct a signature to predict the prognosis and immunotherapy of LUAD patients.

Figure 1. The flowchart of this investigation. To integrate multiple datasets to explore the role of ICD in the TME of LUAD, and to construct a signature to predict the prognosis and immunotherapy of LUAD patients.

The single-cell analysis process was executed as follows: (A) The distribution of gene expression levels, sequencing depth, the percentage of red blood cell genes, the percentage of mitochondrial genes, and the percentage of ribosome genes was assessed across the 12 samples. (B) Correlations between sequencing depth and gene expression levels, as well as the percentages of mitochondrial genes, red blood cell genes, and ribosome genes, were examined. (C) The distribution of observation samples slated for clustering based on cycle-related marker scores was analyzed. (D) The expression of marker genes within diverse subpopulations of classical cell types was depicted. (E) The distribution of distinct cell populations across the 12 LUAD samples was illustrated. (F) An UMAP plot displayed the comprehensive composition of cell types. (G) Differential expression markers between various enriched cell types, aligned with pathway profiles, were identified. (H) Divergent ICD activity within distinct cell groups was ascertained.

Figure 2. The single-cell analysis process was executed as follows: (A) The distribution of gene expression levels, sequencing depth, the percentage of red blood cell genes, the percentage of mitochondrial genes, and the percentage of ribosome genes was assessed across the 12 samples. (B) Correlations between sequencing depth and gene expression levels, as well as the percentages of mitochondrial genes, red blood cell genes, and ribosome genes, were examined. (C) The distribution of observation samples slated for clustering based on cycle-related marker scores was analyzed. (D) The expression of marker genes within diverse subpopulations of classical cell types was depicted. (E) The distribution of distinct cell populations across the 12 LUAD samples was illustrated. (F) An UMAP plot displayed the comprehensive composition of cell types. (G) Differential expression markers between various enriched cell types, aligned with pathway profiles, were identified. (H) Divergent ICD activity within distinct cell groups was ascertained.

Identifying the genes most relevant to ICD

In Figure 3A, the visualization of ICD activity in cells was presented using the AUCell algorithm. In order to offer a more visual portrayal of the ICD scores in cells, the cells were dichotomized based on the median value of these scores. Through the execution of GSVA enrichment analysis, the observation emerged that cells categorized within the high ICD group manifested noteworthy enrichment in pathways related to Allograft Rejection, Inflammatory Response, and Interferon Gamma Response (Figure 3B). In contrast, the Estrogen Response Early pathway demonstrated marked enrichment in the low ICD group. Moving forward, the marker genes that exhibited differential expression in the high and low ICD groups, as acquired from Bulk RNA-seq data, were harnessed to formulate prognostic models via COX and Lasso regression analyses. The process of identifying pivotal prognostic variables was depicted in Figure 3C, 3D. The subsequent Figure 3E depicted the hazard ratio (HR) values corresponding to each variable integrated into the model, while Figure 3F showcased the coefficients attributed to specific variables within the model.

Construction of the signature. Model construction. (A) Distribution of ICD activity in the cell population. (B) GSVA enrichment analysis between high and low ICD groupings. (C, D) Variable screening and model construction for differential genes between high and low ICD subgroups using lasso regression. (E) HR values of model genes. (F) Coefficients of model genes.

Figure 3. Construction of the signature. Model construction. (A) Distribution of ICD activity in the cell population. (B) GSVA enrichment analysis between high and low ICD groupings. (C, D) Variable screening and model construction for differential genes between high and low ICD subgroups using lasso regression. (E) HR values of model genes. (F) Coefficients of model genes.

Evaluation of the model

The computation of individual patient risk scores entailed the multiplication of model gene expression values by their respective coefficients. Subsequent to this, patient stratification into high and low-risk categories was conducted based on the median risk score value. It is noteworthy that both the TCGA dataset and the five GEO cohorts consistently exhibited inferior prognosis outcomes within the high-risk group, thereby substantiating the precision and robustness of the developed model (Figure 4A). The distribution of samples was visually presented through a PCA plot, elucidated in Figure 4B. Furthermore, the spatial arrangement of samples within the high- and low-risk groups was depicted in the t-SNE plot, accessible in Supplementary Figure 2.

Assessment of risk models. (A) Kaplan-Meier survival analysis of signatures in the TCGA and five GEO datasets. (B) Observation of the distribution of samples in high- and low-risk groups using PCA analysis.

Figure 4. Assessment of risk models. (A) Kaplan-Meier survival analysis of signatures in the TCGA and five GEO datasets. (B) Observation of the distribution of samples in high- and low-risk groups using PCA analysis.

The construction of a nomogram

The predictive performance of the model regarding prognosis was meticulously evaluated through the utilization of ROC curves across both the TCGA dataset and the five GEO datasets, thereby substantiating the robustness and consistency of our developed model within the ambit of multi-dataset validation. Remarkably, a substantial majority of the AUC values surpassed the threshold of 0.6 across both the TCGA and GEO datasets, as evidenced in Figure 5A. To comprehensively assess the risk associated with TCGA-LUAD patients, a nomogram was meticulously formulated, integrating crucial clinical features alongside risk group stratification. As depicted in Figure 5B, this nomogram visually encapsulates patient stage, age, and risk score, thereby serving as a valuable tool for enhanced precision in risk assessment and informing forthcoming treatment strategies. Inclusive of C-index curves, the analysis demonstrated the superior predictive performance of the nomogram scores when compared against alternative clinical features and risk scores, as illustrated in Figure 5C. Moreover, the nomogram’s accuracy was methodically gauged via prognostic ROC analysis, consistently showcasing markedly superior performance in comparison to other clinical features and risk scores. Notably, the AUC values for various time horizons including 1, 3, 5, 7, and 10 years were 0.781, 0.790, 0.763, 0.777, and 0.862, respectively, underscoring the substantial predictive capability of the nomogram across diverse temporal scopes (Figure 5D5H).

Building a more accurate nomogram. (A) The ROC curve was used to evaluate the performance of the model in the TCGA and five GEO datasets. (B) Nomogram was constructed by combining clinical features with risk groups. (C) C-index curves were utilized to evaluate the predictive performance of different clinical characteristics, nomogram scores, and risk scores. (D–H) ROC curves for 1, 3, 5, 7, and 10 years showed AUC values for various clinical factors, risk scores, and nomogram scores.

Figure 5. Building a more accurate nomogram. (A) The ROC curve was used to evaluate the performance of the model in the TCGA and five GEO datasets. (B) Nomogram was constructed by combining clinical features with risk groups. (C) C-index curves were utilized to evaluate the predictive performance of different clinical characteristics, nomogram scores, and risk scores. (DH) ROC curves for 1, 3, 5, 7, and 10 years showed AUC values for various clinical factors, risk scores, and nomogram scores.

Relationship between model and clinical characteristics

To visually depict the distribution of clinical attributes within distinct risk subgroups, a heatmap was meticulously generated. This heatmap seamlessly amalgamates pertinent clinical information with the high- and low-risk subgroups, as vividly depicted in Figure 6A. Notably, the high-risk cohort exhibited a discernible association with more advanced T-stage, N-stage, and clinical stage statuses. These findings, showcased in Figure 6B6F, collectively signify an unfavorable prognostic outlook for patients encompassed within the high-risk subgroup.

This figure provides insights into the connection between risk scores and various clinical attributes. The specific details are as follows: (A) Clinical characteristics are integrated with risk grouping to facilitate an examination of the distribution of clinical attributes and the expression patterns of model genes across distinct risk groups. (B–F) The distribution of age, T-stage, clinical stage, N-stage, and M-stage among LUAD patients in different risk groups is observed and presented as percentages.

Figure 6. This figure provides insights into the connection between risk scores and various clinical attributes. The specific details are as follows: (A) Clinical characteristics are integrated with risk grouping to facilitate an examination of the distribution of clinical attributes and the expression patterns of model genes across distinct risk groups. (BF) The distribution of age, T-stage, clinical stage, N-stage, and M-stage among LUAD patients in different risk groups is observed and presented as percentages.

Functional analysis of ICD-related signatures

The GSVA results unveiled that considerable enrichment in pathways such as glycolysis, mTORC1 signaling, and the G2M checkpoint was evident among patients belonging to the high-risk group (Figure 7A). The KEGG enrichment analysis map visually showcased significant enrichment in pathways like porphyrin and chlorophyll metabolism, alongside pentose and glucuronate interconversions (Figure 7B). Through ssGSEA analysis of the corresponding data, it was observed that the low-risk group exhibited a greater abundance of infiltrating immune cells, including B-cells and aDCs (Figure 7C).

Enrichment analysis. (A) Exploring differences in GSVA enrichment analysis across risk groups. (B) Enrichment analysis of KEGG for differential genes in high and low risk groups. (C) Study of immune cell infiltration and immune-related pathway enrichment in high- and low-risk groups using ssGSEA.

Figure 7. Enrichment analysis. (A) Exploring differences in GSVA enrichment analysis across risk groups. (B) Enrichment analysis of KEGG for differential genes in high and low risk groups. (C) Study of immune cell infiltration and immune-related pathway enrichment in high- and low-risk groups using ssGSEA.

Assessment of immune microenvironment

The utilization of immune checkpoint blockade has been extensively integrated into the treatment paradigm for advanced LUAD patients. An in-depth analysis was conducted to explore the intricate relationship between risk scores and established immune checkpoint genes (ICGs) within the TCGA-LUAD cohort. Remarkably, discernible trends emerged as the low-risk group exhibited elevated expression levels across a multitude of immune checkpoint genes including CD48, CD40LG, and CD27, as visually represented in Figure 8A. Interactions between model genes, risk scores, and ICGs were thoroughly examined, the outcomes of which were encapsulated within bubble plots for visualization purposes (Figure 8B). The shade of blue denoted negative correlation, while orange hues represented positive correlation. Intriguingly, a positive correlation emerged between the expression levels of model genes and a significant proportion of immune checkpoints. In stark contrast, risk scores displayed a negative correlation with select prevalent immune checkpoint genes, notably BTLA, BTNL2, CD160, and CD244. These revelations collectively underscore the potential promise of immune checkpoint blockade therapy for individuals afflicted by LUAD. To gauge the potential therapeutic implications of immunotherapy within distinct risk strata, an assessment of immune checkpoint inhibitor response prediction scores (IPS) was executed across diverse risk categories. This endeavor aimed to pinpoint those patients who might derive enhanced benefits from immunotherapeutic interventions. The analysis deduced that tumor samples from patients in the low-risk group were more likely to exhibit favorable immune responses to PD-1/PD-L1 or CTLA4 inhibitors, or a combination thereof (Figure 8C). Substantiated by notably higher IPS scores, the low-risk group emerged as candidates with the greatest potential gains from this therapeutic approach. Exploiting data from seven immune infiltration algorithms resident within the TIMER database, an extensive evaluation was carried out to discern discrepancies in immune infiltration between the high- and low-risk groups within the TCGA-LUAD context, as revealed in Supplementary Figure 3. Furthermore, the ESTIMATE methodology was harnessed to validate immune infiltration levels across distinct risk groups. Spearman correlation analysis was engaged to elucidate the intricate interplay between risk scores and immune infiltration scores. Evidently, the low-risk cohort exhibited heightened Immune Scores (P<0.05), as depicted in Figure 8D. Significantly, a notable negative correlation emerged between risk scores and Immune Scores (R = -0.15, FDR<0.001, Figure 8E), underscoring the intrinsic association between risk scores and the extent of immune cell infiltration within the TME. This dynamic relationship might contribute to variances in disease progression and the efficacy of immunotherapy in LUAD patients.

Analysis of immune cell content and immune checkpoint. (A) A box plot showed that differences in immune checkpoint gene expression between high- and low-risk groups. (B) Correlation between model genes and immune checkpoint. (C) The low-risk group has significantly greater IPS, IPS-CTLA4-neg-PD-1-neg, IPS-CTLA4-pos-PD-1-neg, IPS-CTLA4-neg-PD-1-pos, and IPS-CTLA4-pos-PD-1-pos. (D) The violin plot demonstrated the difference in Stromal Score, Immune Score, ESTIMATE Score, and tumor purity calculated using the ESTIMATE algorithm between the two risk subgroups. (E) The correlations in Stromal Score, Immune Score, ESTIMATE Score, and tumor purity were calculated using the ESTIMATE algorithm between the two risk subgroups. Note: *P P P

Figure 8. Analysis of immune cell content and immune checkpoint. (A) A box plot showed that differences in immune checkpoint gene expression between high- and low-risk groups. (B) Correlation between model genes and immune checkpoint. (C) The low-risk group has significantly greater IPS, IPS-CTLA4-neg-PD-1-neg, IPS-CTLA4-pos-PD-1-neg, IPS-CTLA4-neg-PD-1-pos, and IPS-CTLA4-pos-PD-1-pos. (D) The violin plot demonstrated the difference in Stromal Score, Immune Score, ESTIMATE Score, and tumor purity calculated using the ESTIMATE algorithm between the two risk subgroups. (E) The correlations in Stromal Score, Immune Score, ESTIMATE Score, and tumor purity were calculated using the ESTIMATE algorithm between the two risk subgroups. Note: *P < 0.05, **P < 0.01, ***P < 0.001.

Association between mutation and prognosis

The impact of somatic mutations on the outcomes of cancer immunotherapy is inherently variable. The mutational landscape inherent to TCGA-LUAD was subject to meticulous examination, with the outcomes comprehensively illustrated in Figure 9A. Furthermore, an intricate and systematic juxtaposition was orchestrated to meticulously unravel the variances in TMB that might distinguish the high and low-risk groups. Evidently, this analytical endeavor unveiled a discernible augmentation in the frequency of mutations within the high-risk contingent of LUAD (Figure 9B). This assertion was corroborated through an intensive Spearman correlation analysis, which underscored a statistically significant and positive correlation between risk scores and TMB (R = 0.094, P = 0.043, Figure 9C). Utilizing a pragmatic approach predicated on the median benchmarks of TMB and risk scores, a judicious stratification of LUAD patients occurred, subsequently yielding four distinct categories: high-mutation+high-risk, high-mutation+low-risk, low-mutation+high-risk, and low-mutation+low-risk. The ensuing elucidation of outcomes distinctly delineated that individuals contending with high-mutation+low-risk LUAD exhibited a notably promising prognostic trajectory, in stark contrast to their counterparts grappling with low-mutation+high-risk LUAD, who faced a markedly less favorable prognostic landscape (Figure 9D).

Mutation analysis. (A) Mutation landscape of the top 20 genes with mutation frequency in the high-low risk group. (B) TMB differences between high- and low-risk patients. (C) Relationship between risk score and tumor mutation burden. (D) Survival analysis between four different groups (High-TMB+High-Risk, High-TMB+Low-Risk, Low-TMB+High-risk, and Low-TMB+low-risk).

Figure 9. Mutation analysis. (A) Mutation landscape of the top 20 genes with mutation frequency in the high-low risk group. (B) TMB differences between high- and low-risk patients. (C) Relationship between risk score and tumor mutation burden. (D) Survival analysis between four different groups (High-TMB+High-Risk, High-TMB+Low-Risk, Low-TMB+High-risk, and Low-TMB+low-risk).

The role of PITX3

For further analysis, the PITX3 gene, which exhibited the highest HR value, was selected from the model gene. Figure 10A illustrates the elevated expression levels of PITX3 in tumor samples. Additionally, the high expression of the PITX3 gene was associated with a poorer prognosis compared to the low expression group (P<0.001, Figure 10B). These conspicuous differences align with our bioinformatic findings, suggesting that PITX3 could serve as a promising biomarker for early detection of LUAD. Subsequently, LUAD patients were divided into high and low expression groups based on PITX3 expression, and differential analysis was performed. Figure 10C showcases the differentially expressed genes, highlighting the top 10 genes with the most significant differences. Furthermore, GO and GSEA enrichment analyses of these differentially expressed genes revealed significant enrichment in pathways related to collagen-containing extracellular matrix and cell-cell junctions in the highly expressed PITX3 group (Figure 10D, 10E).

The role of PITX3. (A) Differential expression of PITX3 in normal and tumor samples of TCGA. (B) Survival differences in high and low PITX3-expressing samples. (C) Differential genes between samples with high and low expression of PITX3, results are presented in volcano plots. (D) Pathways of differential gene enrichment in GO. (E, F) GSEA enrichment analysis in the high- and low-PITX3 expressing groups.

Figure 10. The role of PITX3. (A) Differential expression of PITX3 in normal and tumor samples of TCGA. (B) Survival differences in high and low PITX3-expressing samples. (C) Differential genes between samples with high and low expression of PITX3, results are presented in volcano plots. (D) Pathways of differential gene enrichment in GO. (E, F) GSEA enrichment analysis in the high- and low-PITX3 expressing groups.

Discussion

Despite the progress achieved in comprehending LC’s susceptibility, progression, immune regulation, and treatment options, it continues to retain its position as the primary cause of cancer-related mortality [42]. NSCLC represents the most common form within the pathological classification of LC. Specifically, LUAD is the predominant subtype of NSCLC and the most frequent primary LC [43]. Despite ongoing exploration and research in LC treatment, the lack of reliable early prognostic indicators hinders treatment efficacy [44]. Regulated cell death (RCD), a form of cell death controlled by gene-encoded molecular mechanisms, has long been recognized as an event associated with immune suppression or even tolerance [45]. Releasing a diverse array of damage-associated molecular patterns (DAMPs), the phenomenon of RCD operates to orchestrate the recruitment and activation of antigen-presenting cells, thereby offering significant potential for the development of innovative approaches to cancer therapy [25]. The overarching concept of cancer immunotherapy aims to harness the immune system's intrinsic capabilities to induce an immune response against tumors. Through the initiation of ICD, the sustained effectiveness of anticancer medications can be attained, thereby merging direct elimination of cancer cells with the stimulation of antitumor immune responses [46]. Noteworthy antecedent investigations have underscored the prognostic significance of genes implicated in ICD within patients afflicted by head and neck squamous cell carcinoma (HNSCC) [47, 48]. However, the precise role of ICD within the context of LUAD remains enigmatic. The primary goal of this study is to meticulously unravel the prognostic implications and therapeutic prospects that ICD potentially harbors within the realm of LUAD.

Unprecedented resolution is provided by single-cell analysis for the exploration of intratumoral heterogeneity, cell differentiation trajectories, and intercellular communication, offering a range of promising applications [29]. By examining cell clustering within scRNA-seq datasets, genes exclusively expressed in tumor cells were identified. This shift redirects attention away from conventional comparisons between tumors and normal tissues, as observed in previous database analyses, towards the investigation of disparities among tumor cells themselves [30]. In this investigation, pivotal regulatory genes governing ICD were pinpointed through the utilization of scRNA-seq data. Leveraging clinical data, a prognostic signature composed of 14 genes was formulated. Subsequent calculations of risk scores facilitated the categorization of LUAD patients into high-risk and low-risk groups. When comparing survival curves between the high-risk and low-risk cohorts, a favorable prognostic trend was evident for patients within the low-risk group in the TCGA cohort (P<0.0001). Consistent survival outcomes were also observed in the TCGA-test group, TCGA-train group, and GEO verification group (P<0.05). However, no statistically significant differentiation between the high-risk and low-risk cohorts was observed within the GEO30219 group (P=0.057). The considerable accuracy of the constructed prognostic signature in assessing the prognosis of LUAD patients within intervals of 1, 3, 5, 7, and 10 years in both the TCGA and GEO cohorts was confirmed through ROC curve analysis. Impressively, these findings closely aligned with the predictions of the nomogram.

Functional enrichment analysis unveiled significant enrichment of ICD-related signatures in the glycolysis, mTORC1-signaling, and G2M-checkpoint pathways. Altered energy metabolism, such as glycolysis, represents a distinct feature of cancer cells, serving as one of the “hallmarks of cancer.” Glycolysis metabolism provides selective advantages to cancer cells in environments with limited nutrient supply [49]. TBK1 inhibition impaired CRC cell proliferation, migration, drug resistance, and tumor growth. Overexpression of TBK1 inhibited mTORC1-signaling activation in CRC [50]. Although there is no clear evidence that mTORC1-signaling is involved in the progression of LUAD, we can speculate that mTORC1-signaling also plays an important role in LUAD. Significant prominence is accorded to the G2M-CHECKPOINT pathway in a spectrum of related malignancies, encompassing pancreatic cancer [51], colorectal cancer [52], breast cancer [53], among others. It is discernible that the intricate regulatory mechanisms governing the G2M-CHECKPOINT pathway might equally exert an influential impact on the trajectory of LUAD progression.

An all-encompassing scrutiny of the immune cell infiltration within TME can serve as a conduit for comprehending the intricate mechanisms underpinning cancer's evasive maneuvers against immune responses, thereby opening avenues for the formulation of innovative therapeutic strategies [54]. The meticulous investigation of immune cell infiltration disparities between the high-risk and low-risk cohorts unveiled a conspicuous augmentation of immune cell abundance within the low-risk group in contrast to their high-risk counterparts. It is noteworthy that prior investigations have underscored the interrelation between immune checkpoint gene expression levels and the efficacy of immunotherapeutic interventions [55]. Delving into the variations of immune checkpoint genes between the high-risk and low-risk groups proposed the potential for the low-risk group to glean added therapeutic advantage from targeting supplementary immune checkpoints. Within the realm of cancer treatment, therapeutic modalities geared towards the TME have risen to prominence, recognized for their promise in leveraging the pivotal role that the TME holds in steering tumor progression and influencing responses to standard therapeutic regimens [5659]. Notably, assessment of TME scoring unveiled a notable elevation in immune scores within the TME of the low-risk cohort, distinct from their high-risk counterparts, with a statistically significant divergence. This compelling observation lends credence to the notion that patients harboring low-risk profiles might exhibit heightened responsiveness to immunotherapy. Adding to these insights, a systematic analysis appraising the comparative immunotherapy advantages between the high-risk and low-risk groups corroborated the prospective amplified therapeutic gains for individuals encompassed within the low-risk classification.

Although PITX3’s role in brain development is well established, its involvement in tumorigenesis remains largely unknown. Reports have indicated high methylation levels of PITX3 in breast cancer [60]. PITX3 promoter methylation also impacts the recurrence-free survival of prostate cancer patients [61]. However, the precise role of PITX3 in LUAD progression has yet to be elucidated. Our study revealed that PITX3 exhibited the highest HR among all prognostic signature genes. PITX3 may contribute to LUAD progression. The survival analysis yielded compelling evidence linking heightened PITX3 expression with unfavorable prognostic outcomes among LUAD patients, thus underlining its potential as a viable therapeutic target. Of noteworthy significance, the ECM-receptor interaction pathway emerges as a pivotal participant in a multitude of tumor-associated processes, encompassing shedding, adhesion, degradation, motility, and proliferation [62]. This pathway's pivotal role in tumor invasion and metastasis has been decisively affirmed within the context of gastric cancer [63]. Furthermore, comprehensive investigations have unequivocally demonstrated that the intricate interplay between the extracellular matrix (ECM) and cellular receptors constitutes a fundamental conduit underpinning the progression and metastasis of colorectal cancer [64]. Within the high-risk subgroup, PITX3 manifested conspicuous enrichment within the ECM-receptor interaction pathway. This observation notably intimates that PITX3's influence might extend to interfering with the trajectory of LUAD progression by active involvement within the ECM-receptor interaction pathway, potentially engendering adverse prognostic implications for patients grappling with LUAD.

Despite its limitations, valuable insights are offered by this study into the expression of prognostic genes associated with LUAD, thereby laying the groundwork for subsequent investigations into their underlying molecular mechanisms.

In summary, the ICD-related signatures formulated within this study facilitate the prognostic prediction of LUAD patients grounded in the expression of pertinent genes. Furthermore, this work paves the way for potential applications of immunotherapy in the realm of LUAD treatment. It is important to note, however, that further experimental endeavors are warranted to substantiate and validate these findings.

Author Contributions

PZ, HZ, and JT contributed conception and design of the study; HZ and JT collected LUAD tissue samples and obtained ethical approval; JZ and QR collected the data; JX and HC performed the statistical analysis; XG wrote the first draft of the manuscript; HL, JL and CH gave the final approval of the version to be submitted. All authors contributed to the manuscript and approved the submitted version.

Acknowledgments

We are very grateful for data provided by databases such as TCGA, GEO. Thanks to reviewers and editors for their sincere comments.

Conflicts 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.

Ethical Statement and Consent

The human experiments in this study received approval from the Ethics Committee of the First Affiliated Hospital of Nanjing Medical University (Approval No. 2019-SR-156). All participants provided informed consent before taking part. The study adhered to the principles of the Declaration of Helsinki and followed ethical guidelines.

Funding

No funding was provided for this study.

References

  • 1. Nasim F, Sabath BF, Eapen GA. Lung Cancer. Med Clin North Am. 2019; 103:463–73. https://doi.org/10.1016/j.mcna.2018.12.006 [PubMed]
  • 2. Jain D, Roy-Chowdhuri S. Advances in cytology of lung cancer. Semin Diagn Pathol. 2021; 38:109–15. https://doi.org/10.1053/j.semdp.2021.05.001 [PubMed]
  • 3. Denisenko TV, Budkevich IN, Zhivotovsky B. Cell death-based treatment of lung adenocarcinoma. Cell Death Dis. 2018; 9:117. https://doi.org/10.1038/s41419-017-0063-y [PubMed]
  • 4. Jones GS, Baldwin DR. Recent advances in the management of lung cancer. Clin Med (Lond). 2018 (Suppl 2); 18:s41–6. https://doi.org/10.7861/clinmedicine.18-2-s41 [PubMed]
  • 5. The La. Lung cancer: some progress, but still a lot more to do. Lancet. 2019; 394:1880. https://doi.org/10.1016/S0140-6736(19)32795-3 [PubMed]
  • 6. Wu F, Fan J, He Y, Xiong A, Yu J, Li Y, Zhang Y, Zhao W, Zhou F, Li W, Zhang J, Zhang X, Qiao M, et al. Single-cell profiling of tumor heterogeneity and the microenvironment in advanced non-small cell lung cancer. Nat Commun. 2021; 12:2540. https://doi.org/10.1038/s41467-021-22801-0 [PubMed]
  • 7. Chi H, Zhao S, Yang J, Gao X, Peng G, Zhang J, Xie X, Song G, Xu K, Xia Z, Chen S, Zhao J. T-cell exhaustion signatures characterize the immune landscape and predict HCC prognosis via integrating single-cell RNA-seq and bulk RNA-sequencing. Front Immunol. 2023; 14:1137025. https://doi.org/10.3389/fimmu.2023.1137025 [PubMed]
  • 8. Zhang J, Peng G, Chi H, Yang J, Xie X, Song G, Tran LJ, Xia Z, Tian G. CD8 + T-cell marker genes reveal different immune subtypes of oral lichen planus by integrating single-cell RNA-seq and bulk RNA-sequencing. BMC Oral Health. 2023; 23:464. https://doi.org/10.1186/s12903-023-03138-0 [PubMed]
  • 9. Song G, Peng G, Zhang J, Song B, Yang J, Xie X, Gou S, Zhang J, Yang G, Chi H, Tian G. Uncovering the potential role of oxidative stress in the development of periodontitis and establishing a stable diagnostic model via combining single-cell and machine learning analysis. Front Immunol. 2023; 14:1181467. https://doi.org/10.3389/fimmu.2023.1181467 [PubMed]
  • 10. Chi H, Yang J, Peng G, Zhang J, Song G, Xie X, Xia Z, Liu J, Tian G. Circadian rhythm-related genes index: A predictor for HNSCC prognosis, immunotherapy efficacy, and chemosensitivity. Front Immunol. 2023; 14:1091218. https://doi.org/10.3389/fimmu.2023.1091218 [PubMed]
  • 11. Chi H, Jiang P, Xu K, Zhao Y, Song B, Peng G, He B, Liu X, Xia Z, Tian G. A novel anoikis-related gene signature predicts prognosis in patients with head and neck squamous cell carcinoma and reveals immune infiltration. Front Genet. 2022; 13:984273. https://doi.org/10.3389/fgene.2022.984273 [PubMed]
  • 12. Chi H, Xie X, Yan Y, Peng G, Strohmer DF, Lai G, Zhao S, Xia Z, Tian G. Natural killer cell-related prognosis signature characterizes immune landscape and predicts prognosis of HNSCC. Front Immunol. 2022; 13:1018685. https://doi.org/10.3389/fimmu.2022.1018685 [PubMed]
  • 13. Suvà ML, Tirosh I. Single-Cell RNA Sequencing in Cancer: Lessons Learned and Emerging Challenges. Mol Cell. 2019; 75:7–12. https://doi.org/10.1016/j.molcel.2019.05.003 [PubMed]
  • 14. Zhao S, Zhang L, Ji W, Shi Y, Lai G, Chi H, Huang W, Cheng C. Machine learning-based characterization of cuprotosis-related biomarkers and immune infiltration in Parkinson’s disease. Front Genet. 2022; 13:1010361. https://doi.org/10.3389/fgene.2022.1010361 [PubMed]
  • 15. Zhao S, Chi H, Ji W, He Q, Lai G, Peng G, Zhao X, Cheng C. A Bioinformatics-Based Analysis of an Anoikis-Related Gene Signature Predicts the Prognosis of Patients with Low-Grade Gliomas. Brain Sci. 2022; 12:1349. https://doi.org/10.3390/brainsci12101349 [PubMed]
  • 16. Song B, Wu P, Liang Z, Wang J, Zheng Y, Wang Y, Chi H, Li Z, Song Y, Yin X, Yu Z, Song B. A Novel Necroptosis-Related Gene Signature in Skin Cutaneous Melanoma Prognosis and Tumor Microenvironment. Front Genet. 2022; 13:917007. https://doi.org/10.3389/fgene.2022.917007 [PubMed]
  • 17. Song B, Chi H, Peng G, Song Y, Cui Z, Zhu Y, Chen G, Wu J, Liu W, Dong C, Wang Y, Xu K, Yu Z, Song B. Characterization of coagulation-related gene signature to predict prognosis and tumor immune microenvironment in skin cutaneous melanoma. Front Oncol. 2022; 12:975255. https://doi.org/10.3389/fonc.2022.975255 [PubMed]
  • 18. Chi H, Chen H, Wang R, Zhang J, Jiang L, Zhang S, Jiang C, Huang J, Quan X, Liu Y, Zhang Q, Yang G. Proposing new early detection indicators for pancreatic cancer: Combining machine learning and neural networks for serum miRNA-based diagnostic model. Front Oncol. 2023; 13:1244578. https://doi.org/10.3389/fonc.2023.1244578 [PubMed]
  • 19. Shen Y, Chi H, Xu K, Li Y, Yin X, Chen S, Yang Q, He M, Zhu G, Li X. A Novel Classification Model for Lower-Grade Glioma Patients Based on Pyroptosis-Related Genes. Brain Sci. 2022; 12:700. https://doi.org/10.3390/brainsci12060700 [PubMed]
  • 20. Lei Y, Tang R, Xu J, Wang W, Zhang B, Liu J, Yu X, Shi S. Applications of single-cell sequencing in cancer research: progress and perspectives. J Hematol Oncol. 2021; 14:91. https://doi.org/10.1186/s13045-021-01105-2 [PubMed]
  • 21. Chi H, Gao X, Xia Z, Yu W, Yin X, Pan Y, Peng G, Mao X, Teichmann AT, Zhang J, Tran LJ, Jiang T, Liu Y, et al. FAM family gene prediction model reveals heterogeneity, stemness and immune microenvironment of UCEC. Front Mol Biosci. 2023; 10:1200335. https://doi.org/10.3389/fmolb.2023.1200335 [PubMed]
  • 22. Jin W, Yang Q, Chi H, Wei K, Zhang P, Zhao G, Chen S, Xia Z, Li X. Ensemble deep learning enhanced with self-attention for predicting immunotherapeutic responses to cancers. Front Immunol. 2022; 13:1025330. https://doi.org/10.3389/fimmu.2022.1025330 [PubMed]
  • 23. Kroemer G, Galluzzi L, Kepp O, Zitvogel L. Immunogenic cell death in cancer therapy. Annu Rev Immunol. 2013; 31:51–72. https://doi.org/10.1146/annurev-immunol-032712-100008 [PubMed]
  • 24. Kroemer G, Galassi C, Zitvogel L, Galluzzi L. Immunogenic cell stress and death. Nat Immunol. 2022; 23:487–500. https://doi.org/10.1038/s41590-022-01132-2 [PubMed]
  • 25. Fucikova J, Kepp O, Kasikova L, Petroni G, Yamazaki T, Liu P, Zhao L, Spisek R, Kroemer G, Galluzzi L. Detection of immunogenic cell death and its relevance for cancer therapy. Cell Death Dis. 2020; 11:1013. https://doi.org/10.1038/s41419-020-03221-2 [PubMed]
  • 26. Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics. 2012; 28:882–3. https://doi.org/10.1093/bioinformatics/bts034 [PubMed]
  • 27. Cao Y, Fu L, Wu J, Peng Q, Nie Q, Zhang J, Xie X. Integrated analysis of multimodal single-cell data with structural similarity. Nucleic Acids Res. 2022; 50:e121. https://doi.org/10.1093/nar/gkac781 [PubMed]
  • 28. Zheng S, Liang JY, Tang Y, Xie J, Zou Y, Yang A, Shao N, Kuang X, Ji F, Liu X, Tian W, Xiao W, Lin Y. Dissecting the role of cancer-associated fibroblast-derived biglycan as a potential therapeutic target in immunotherapy resistance: A tumor bulk and single-cell transcriptomic study. Clin Transl Med. 2023; 13:e1189. https://doi.org/10.1002/ctm2.1189 [PubMed]
  • 29. Zhao Z, Ding Y, Tran LJ, Chai G, Lin L. Innovative breakthroughs facilitated by single-cell multi-omics: manipulating natural killer cell functionality correlates with a novel subcategory of melanoma cells. Front Immunol. 2023; 14:1196892. https://doi.org/10.3389/fimmu.2023.1196892 [PubMed]
  • 30. Li XY, Zhao ZJ, Wang JB, Shao YH, Hui-Liu, You JX, Yang XT. m7G Methylation-Related Genes as Biomarkers for Predicting Overall Survival Outcomes for Hepatocellular Carcinoma. Front Bioeng Biotechnol. 2022; 10:849756. https://doi.org/10.3389/fbioe.2022.849756 [PubMed]
  • 31. Wang Y, Zhao ZJ, Kang XR, Bian T, Shen ZM, Jiang Y, Sun B, Hu HB, Chen YS. lncRNA DLEU2 acts as a miR-181a sponge to regulate SEPP1 and inhibit skeletal muscle differentiation and regeneration. Aging (Albany NY). 2020; 12:24033–56. https://doi.org/10.18632/aging.104095 [PubMed]
  • 32. Iasonos A, Schrag D, Raj GV, Panageas KS. How to build and interpret a nomogram for cancer prognosis. J Clin Oncol. 2008; 26:1364–70. https://doi.org/10.1200/JCO.2007.12.9791 [PubMed]
  • 33. Wang Z, Yuan Q, Chen X, Luo F, Shi X, Guo F, Ren J, Li S, Shang D. A prospective prognostic signature for pancreatic adenocarcinoma based on ubiquitination-related mRNA-lncRNA with experimental validation in vitro and vivo. Funct Integr Genomics. 2023; 23:263. https://doi.org/10.1007/s10142-023-01158-1 [PubMed]
  • 34. Liu J, Zhong L, Deng D, Zhang Y, Yuan Q, Shang D. The combined signatures of the tumour microenvironment and nucleotide metabolism-related genes provide a prognostic and therapeutic biomarker for gastric cancer. Sci Rep. 2023; 13:6622. https://doi.org/10.1038/s41598-023-33213-z [PubMed]
  • 35. Cui Y, Yuan Q, Chen J, Jiang J, Guan H, Zhu R, Li N, Liu W, Wang C. Determination and characterization of molecular heterogeneity and precision medicine strategies of patients with pancreatic cancer and pancreatic neuroendocrine tumor based on oxidative stress and mitochondrial dysfunction-related genes. Front Endocrinol (Lausanne). 2023; 14:1127441. https://doi.org/10.3389/fendo.2023.1127441 [PubMed]
  • 36. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, Feng T, Zhou L, Tang W, Zhan L, Fu X, Liu S, Bo X, Yu G. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb). 2021; 2:100141. https://doi.org/10.1016/j.xinn.2021.100141 [PubMed]
  • 37. Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 2018; 28:1747–56. https://doi.org/10.1101/gr.239244.118 [PubMed]
  • 38. Miao Y, Yuan Q, Wang C, Feng X, Ren J, Wang C. Comprehensive Characterization of RNA-Binding Proteins in Colon Adenocarcinoma Identifies a Novel Prognostic Signature for Predicting Clinical Outcomes and Immunotherapy Responses Based on Machine Learning. Comb Chem High Throughput Screen. 2023; 26:163–82. https://doi.org/10.2174/1386207325666220404125228 [PubMed]
  • 39. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, Treviño V, Shen H, Laird PW, Levine DA, Carter SL, Getz G, Stemke-Hale K, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013; 4:2612. https://doi.org/10.1038/ncomms3612 [PubMed]
  • 40. Van Allen EM, Miao D, Schilling B, Shukla SA, Blank C, Zimmer L, Sucker A, Hillen U, Foppen MH, Goldinger SM, Utikal J, Hassel JC, Weide B, et al. Genomic correlates of response to CTLA-4 blockade in metastatic melanoma. Science. 2015; 350:207–11. https://doi.org/10.1126/science.aad0095 [PubMed]
  • 41. Xie J, Zhang J, Tian W, Zou Y, Tang Y, Zheng S, Wong CW, Deng X, Wu S, Chen J, Mo Y, Xie X. The Pan-Cancer Multi-Omics Landscape of FOXO Family Relevant to Clinical Outcome and Drug Resistance. Int J Mol Sci. 2022; 23:15647. https://doi.org/10.3390/ijms232415647 [PubMed]
  • 42. Bade BC, Dela Cruz CS. Lung Cancer 2020: Epidemiology, Etiology, and Prevention. Clin Chest Med. 2020; 41:1–24. https://doi.org/10.1016/j.ccm.2019.10.001 [PubMed]
  • 43. Nooreldeen R, Bach H. Current and Future Development in Lung Cancer Diagnosis. Int J Mol Sci. 2021; 22:8661. https://doi.org/10.3390/ijms22168661 [PubMed]
  • 44. Song J, Sun Y, Cao H, Liu Z, Xi L, Dong C, Yang R, Shi Y. A novel pyroptosis-related lncRNA signature for prognostic prediction in patients with lung adenocarcinoma. Bioengineered. 2021; 12:5932–49. https://doi.org/10.1080/21655979.2021.1972078 [PubMed]
  • 45. Galluzzi L, Vitale I, Warren S, Adjemian S, Agostinis P, Martinez AB, Chan TA, Coukos G, Demaria S, Deutsch E, Draganov D, Edelson RL, Formenti SC, et al. Consensus guidelines for the definition, detection and interpretation of immunogenic cell death. J Immunother Cancer. 2020; 8:e000337. https://doi.org/10.1136/jitc-2019-000337 [PubMed]
  • 46. Ahmed A, Tait SW. Targeting immunogenic cell death in cancer. Mol Oncol. 2020; 14:2994–3006. https://doi.org/10.1002/1878-0261.12851 [PubMed]
  • 47. Wang X, Wu S, Liu F, Ke D, Wang X, Pan D, Xu W, Zhou L, He W. An Immunogenic Cell Death-Related Classification Predicts Prognosis and Response to Immunotherapy in Head and Neck Squamous Cell Carcinoma. Front Immunol. 2021; 12:781466. https://doi.org/10.3389/fimmu.2021.781466 [PubMed]
  • 48. Gong X, Chi H, Xia Z, Yang G, Tian G. Advances in HPV-associated tumor management: Therapeutic strategies and emerging insights. J Med Virol. 2023; 95:e28950. https://doi.org/10.1002/jmv.28950 [PubMed]
  • 49. Ganapathy-Kanniappan S, Geschwind JF. Tumor glycolysis as a target for cancer therapy: progress and prospects. Mol Cancer. 2013; 12:152. https://doi.org/10.1186/1476-4598-12-152 [PubMed]
  • 50. Zhou D, Yao Y, Zong L, Zhou G, Feng M, Chen J, Liu G, Chen G, Sun K, Yao H, Liu Y, Shi X, Zhang W, et al. TBK1 Facilitates GLUT1-Dependent Glucose Consumption by suppressing mTORC1 Signaling in Colorectal Cancer Progression. Int J Biol Sci. 2022; 18:3374–89. https://doi.org/10.7150/ijbs.70742 [PubMed]
  • 51. Yang R, Liang X, Wang H, Guo M, Shen H, Shi Y, Liu Q, Sun Y, Yang L, Zhan M. The RNA methyltransferase NSUN6 suppresses pancreatic cancer development by regulating cell proliferation. EBioMedicine. 2021; 63:103195. https://doi.org/10.1016/j.ebiom.2020.103195 [PubMed]
  • 52. Gong Y, Liu Y, Wang T, Li Z, Gao L, Chen H, Shu Y, Li Y, Xu H, Zhou Z, Dai L. Age-Associated Proteomic Signatures and Potential Clinically Actionable Targets of Colorectal Cancer. Mol Cell Proteomics. 2021; 20:100115. https://doi.org/10.1016/j.mcpro.2021.100115 [PubMed]
  • 53. Gandhi S, Oshi M, Murthy V, Repasky EA, Takabe K. Enhanced Thermogenesis in Triple-Negative Breast Cancer Is Associated with Pro-Tumor Immune Microenvironment. Cancers (Basel). 2021; 13:2559. https://doi.org/10.3390/cancers13112559 [PubMed]
  • 54. Zhang Y, Zhang Z. The history and advances in cancer immunotherapy: understanding the characteristics of tumor-infiltrating immune cells and their therapeutic implications. Cell Mol Immunol. 2020; 17:807–21. https://doi.org/10.1038/s41423-020-0488-6 [PubMed]
  • 55. Ahluwalia P, Ahluwalia M, Mondal AK, Sahajpal N, Kota V, Rojiani MV, Rojiani AM, Kolhe R. Immunogenomic Gene Signature of Cell-Death Associated Genes with Prognostic Implications in Lung Cancer. Cancers (Basel). 2021; 13:155. https://doi.org/10.3390/cancers13010155 [PubMed]
  • 56. Bejarano L, Jordāo MJC, Joyce JA. Therapeutic Targeting of the Tumor Microenvironment. Cancer Discov. 2021; 11:933–59. https://doi.org/10.1158/2159-8290.CD-20-1808 [PubMed]
  • 57. Xiong J, Chi H, Yang G, Zhao S, Zhang J, Tran LJ, Xia Z, Yang F, Tian G. Revolutionizing anti-tumor therapy: unleashing the potential of B cell-derived exosomes. Front Immunol. 2023; 14:1188760. https://doi.org/10.3389/fimmu.2023.1188760 [PubMed]
  • 58. Gong X, Chi H, Strohmer DF, Teichmann AT, Xia Z, Wang Q. Exosomes: A potential tool for immunotherapy of ovarian cancer. Front Immunol. 2023; 13:1089410. https://doi.org/10.3389/fimmu.2022.1089410 [PubMed]
  • 59. Zhao Y, Wei K, Chi H, Xia Z, Li X. IL-7: A promising adjuvant ensuring effective T cell responses and memory in combination with cancer vaccines? Front Immunol. 2022; 13:1022808. https://doi.org/10.3389/fimmu.2022.1022808 [PubMed]
  • 60. Dietrich D, Lesche R, Tetzner R, Krispin M, Dietrich J, Haedicke W, Schuster M, Kristiansen G. Analysis of DNA methylation of multiple genes in microdissected cells from formalin-fixed and paraffin-embedded tissues. J Histochem Cytochem. 2009; 57:477–89. https://doi.org/10.1369/jhc.2009.953026 [PubMed]
  • 61. Holmes EE, Goltz D, Sailer V, Jung M, Meller S, Uhl B, Dietrich J, Röhler M, Ellinger J, Kristiansen G, Dietrich D. PITX3 promoter methylation is a prognostic biomarker for biochemical recurrence-free survival in prostate cancer patients after radical prostatectomy. Clin Epigenetics. 2016; 8:104. https://doi.org/10.1186/s13148-016-0270-x [PubMed]
  • 62. Bao Y, Wang L, Shi L, Yun F, Liu X, Chen Y, Chen C, Ren Y, Jia Y. Transcriptome profiling revealed multiple genes and ECM-receptor interaction pathways that may be associated with breast cancer. Cell Mol Biol Lett. 2019; 24:38. https://doi.org/10.1186/s11658-019-0162-0 [PubMed]
  • 63. Yan P, He Y, Xie K, Kong S, Zhao W. In silico analyses for potential key genes associated with gastric cancer. PeerJ. 2018; 6:e6092. https://doi.org/10.7717/peerj.6092 [PubMed]
  • 64. Nersisyan S, Novosad V, Engibaryan N, Ushkaryov Y, Nikulin S, Tonevitsky A. ECM-Receptor Regulatory Network and Its Prognostic Role in Colorectal Cancer. Front Genet. 2021; 12:782699. https://doi.org/10.3389/fgene.2021.782699 [PubMed]