2.1.DEGs Differentiation and Division of Patients into Two Subtypes
We filtered 34 ICD-related genes in 598 LUAD sets in 59 normal tissues and 539 LUAD samples from TCGA, according to an extensive literature summarized by Abhishek et al. (11): ATG5, BAX CALR, CASP1, CASP8, CD4, CD8A, CD8B, CXCR3, EIF2AK3, ENTPD1, FOXP3, HMGB1, HSP90AA1, IFNA1, IFNB1, IFNG, IFNGR1, IL10, IL17A, IL17, RAIL1B, IL1R1, IL6, LY96, MYD88, NLRP3, NT5E, P2RX7, PDIA3, PIK3CA, PRF1, TLR4, and TNF. Then, the expression of these ICD-related genes was analyzed with the “pheatmap” R package to visualize the variance between tumor and normal samples. Twenty-nine genes were differentially expressed in tumor samples: 12 were upregulated (CD8A, BAX, EIF2AK3, HSP90AA1, CASP8, CALR, PDIA3, IFNG, FOXP3, NT5E, IL17A, and IFNB1), and 17 were downregulated (iIL6, TLR4, PRF1, IL1B, P2RX7, IFNGR1, CASP1, NLRP3, CD4, IL17RA, IL1R1, MYD88, IFNA1, TNF, ENTPD1, HMGB1, and PIK3CA) (Fig. 1A). Upregulated genes might have an active effect on tumor development, while downregulated genes might have an inhibitory effect. The 29 ICD-related DEGs were used for PPI analysis via the STRING database (Fig. 1B). Then, we determined ICD-related clusters of LUAD by consensus clustering. LUAD samples were grouped into two subgroups that maximized consensus within clusters while minimizing ambiguity in cluster assignments and considering the feasibility of clinical prognosis analysis (Fig. 1C, D, E,F). LUAD samples were divided into two cohorts (C1 and C2), rendered in the landscape to present the DEGs (Fig. 1G). ICD-related genes were downregulated in C1, representing the ICD-low subgroup, while ICD-related genes were upregulated in C2, comprising the ICD-high subgroup. Meanwhile, most of DEGs shows differences in gene expression, CNV frequency and survival(Fig. 2A). Their gain copy levels are higher than loss copy levels(Fig. 2B). 16 DEGs have significant differences in survival analysis, amoung which the low expression of ENTPD1, IL17A, P2RX7, TLR4, TNF, CD4, CD8A, CXCR3, IFNG, IL1R1 and IL10 indicate better outcome(Fig. 2C).
2.2.Inter-Groups Differentiation and Enrichment Analyses
All ICD-related genes were visualized in the landscape (Fig. 3A), where red cubes represent overexpression, while blue ones indicate underexpression. The volcano plot showed the gene distribution: red plot: upregulated; dark plot: unregulated; green plot: downregulated (Fig. 3B). The GO enrichment analysis comprises different classes: Biological Processes (BP), Cellular Components (CC), and Molecular Functions (MF) (Fig. 3C). This analysis revealed that DEGs were mainly enriched in BP related to immune function, such as receptor-ligand activity (GO:0048018); signaling receptor activator activity (GO:0030546) in the MF subgroup; and collagen − containing extracellular matrix (GO:0062023) in CC, which might contribute to immunotherapy. The same result can be seen in the dot plot representing gene ratios (Fig. 3D). Cytokine-cytokine receptor interaction was significantly enriched in the KEGG enrichment analysis. This pathway influences immune response processes such as adaptive inflammatory host defenses and cell death (Fig. 3E). The GSVA indicated that eight hallmark gene sets were differentially enriched between the two subgroups (Fig. 4A), all upregulated in the ICD-high subgroup (Fig. 4B). These pathways contribute to the response to Interferon-gamma, alpha interferon proteins, tumor necrosis factor, transplant rejection, inflammatory response, and KRAS [28], indicating a higher level of tumor immune response in the ICD-high subgroup. The GSEA indicated that the cytokine-cytokine receptor interaction pathway was significantly enriched for ICD-related genes in the high-risk subgroup, which is also related to immune response (Fig. 5A, B). Therefore, the ICD-high subgroup was strongly related to immune pathways.
2.3.Somatic Mutation and TME
Further, we explored the distinct somatic mutation profiles between the two subgroups. To evaluate the somatic mutation differences between ICD-high and low subgroups, we drew the oncoplot for each subgroup (Fig. 5C, D). The frequency varied, as shown in the two charts. Despite TP53 and TTN that were high-frequency mutation genes, other genes, such as MUC16 and CSMD3, also differed. The ICD-high subgroup presented a lower mutation rate than the ICD-low subgroup. The main mutation modes in the two subgroups are missense mutations and multi-hits. Further, we analyzed the tumor mutation burden (TMB) of the two subgroups. The ICD-low subgroup presented a statistically higher TMB than the ICD-high subgroup (Fig. 6A). Thus, upregulation of ICD-related genes might lower the TMB and inhibit tumor development. Then, we analyzed the immune microenvironment based on the ICD-related gene data (Fig. 6B, C, D, E). The four items differed between ICD-high and low subgroups. The estimate (Fig. 6B), stromal (Fig. 6C), and immune (Fig. 6D) scores were higher, and tumor purity (Fig. 6E) was lower in the ICD-high subgroup than in the ICD-low subgroup. Estimate scores are the sum of immune and stromal scores. These results demonstrated ICD’s ability to recruit immune infiltration. We obtained immune cell infiltration fractions using CIBERSORT to deconvolute samples from both subgroups. The content and relevance of 22 types of immune cells in ICD samples are presented (Fig. 7A). The correlation analysis of immune cells is shown (Fig. 7B). Thus, we contrasted ICD-high and low subgroups for each immune cell (Fig. 7C). CD4 memory resting T cells, monocytes, resting dendritic cells and resting mast cells were upregulated, while plasma cells, CD8T cells, follicular helper T cells, resting NK cells, and M0 macrophages were downregulated compared to the ICD-low subgroup with different p values (*). In the ICD-high subgroup, resting NK cells and M0 macrophages were downregulated, implying the improvement of the immune microenvironment. Furthermore, all human leukocyte antigen (HLA) genes listed below were overexpressed (Fig. 7D), and most checkpoints were similar (Fig. 7E) in the ICD-high subtype. Hence, the ICD-high subgroup was associated with the immune-hot phenotype, and the ICD-low subgroup was linked to the immune-cold phenotype.
2.4. ICD Prognostic Model
Six ICD-related genes were associated with the patient survival (p < 0.05): ENTPD1, HSP90AA1, IL17A, P2RX7, PIK3CA, and NT5E. These ICD-related genes were used to construct the classification model by LASSO regression analysis (Fig. 8A). Thus, we developed the classification based on the algorithm: Risk score = (0.1671) * HSP90AA1 + (0.0860) * NT5E + (-0.1682) *ENTPD1 + (0.0107) *IL17A + (0.0830) * P2RX7 + (0.0613) * PIK3CA. Each sample received a Risk Score. The landscape shows the individual risk score of the six ICD-related genes in increasing total Risk Score order (Fig. 8B). ENTPD1, IL17A, and P2RX7 were downregulated, while HSP90AA1, PIK3CA, and NT5E were upregulated, indicating that ENTPD1, IL17A, and P2RX7 might negatively affect tumor development and result in better outcome. We divided them into high- and low-risk subgroups (Fig. 8C). Moreover, the survival distribution indicated more deaths as the risk score increased (Fig. 8D). Then, we applied the Risk Score for TCGA LUAD (Fig. 8E) and GEO LUAD samples (Fig. 8F). In TCGA, a high-risk score was related to poor overall survival (OS), further confirmed in the GEO cohort. This corresponds to the DEGs survival analysis. To verify the independent prognostic value of this classification, we used univariate and multivariate Cox regression analyses, which indicated the independent prognostic value of the ICD risk signature (p < 0.05). The univariate Cox regression analysis showed that a high ICD risk score was significantly associated with poor OS (Fig. 9A). The multivariate Cox regression analysis showed that the ICD risk score was an independent prognostic factor for LUAD patients (Fig. 9B).
2.5.TME and Therapy Analysis
ICD plays a vital biological role in anti-cancer immunological responses. Thus, it has a tight connection with the TME. In the immune correlation analysis, Plasma cells, T regulatory cells (Tregs), CD4 memory resting T cells, B cells memory, B cells naïve, endritic cells activated varied with the Risk Score (p < 0.05)(Fig. 9C-H). As the Risk Score increased, the concentration of Tregs decreased, and CD4 memory resting T cells increased (Fig. 9C, D). TIDE was used to evaluate the prognostic value of the Risk Score in immunotherapy based on a comprehensive analysis of the tumor expression spectrum to predict the efficacy of ICD therapy (Fig. 10A). The high-risk score group responded better to immunotherapy. Then, we used the drug information in the GDSC database to evaluate the IC50 of common chemotherapy drugs(Fig. 10B). We found lower IC50 for 18 chemotherapy drugs in the low-risk score subgroup than in the high-risk score subgroup: BIRB.0796, AZD8055, EHT.1864, PF.4708671, ATRA, GW.441756, Axitinib, X681640, Elesclomol, ABT.263, PD.173074, Sorafenib, BX.795, Metformin, AMG.706, X17.AAG, Pyrimethamine, and NVP.BEZ235. These drugs were statistically sensitive to low-risk score LUAD patients, providing basic selection advice of chemotherapy drugs alone or combined with immunotherapy in clinical practice.