The identification of genes associated T-cell exhaustion and construction of prognostic signature to predict immunotherapy response in lung adenocarcinoma

Abstract

T-cell exhaustion (Tex) is considered to be a reason for immunotherapy resistance and poor prognosis in lung adenocarcinoma. Therefore, we used weighted correlation network analysis to identify Tex-related genes in the cancer genome atlas (TCGA). Unsupervised clustering approach based on Tex-related genes divided patients into cluster 1 and cluster 2. Then, we utilized random forest and the least absolute shrinkage and selection operator to identify nine key genes to construct a riskscore. Patients were classified as low or high-risk groups. The multivariate cox analysis showed the riskscore was an independent prognostic factor in TCGA and GSE72094 cohorts. Moreover, patients in cluster 2 with high riskscore had the worst prognosis. The immune response prediction analysis showed the low-risk group had higher immune, stromal, estimate scores, higher immunophenscore (IPS), and lower tumor immune dysfunction and exclusion score which suggested a better response to immune checkpoint inhibitors (ICIs) therapy in the low-risk group. In the meantime, we included two independent immunotherapy cohorts that also confirmed a better response to ICIs treatment in the low-risk group. Besides, we discovered differences in chemotherapy and targeted drug sensitivity between two groups. Finally, a nomogram was built to facilitate clinical decision making.

Figure 1

WGCNA analysis in TCGA cohort. (A) The coexpression network. (B) The soft threshold power of WGCNA. The left one showed the analysis of the scale-free index for various soft threshold powers. The right one showed the analysis of the average connectivity for various soft threshold powers. (C) Heatmap displayed correlation between module eigengenes and T-cell exhaustion and T-cell dysfunction. (D) Identification of the modules most significantly associated with T-cell exhaustion. (E) Identification of modules most significantly associated with T-cell dysfunction.

Figure 2

Unsupervised consensus cluster analysis for LUAD patients in the TCGA cohort based on T-cell exhaustion related genes. (A) Consensus CDF from k = 2–9. (B) Delta area under the cumulative distribution function (CDF) curve of different clusters ranging from k = 2–9. (C) Consensus matrix for k = 2. (D) The overall survival (OS) probability of the patients in the two clusters. (E) The principal component analysis (PCA) and principal coordinate analysis (PCoA) of two clusters.

Figure 3

Immune analysis of different clusters. (A) Differences in abundance of tumor-infiltrating immune cells in different clusters based on CIBERSORT. (B) Differential mRNA expression of immune checkpoints in different clusters. (C) The stromal, immune and estimate scores between two clusters. (D) Differences in T-cell exhasution score between different clusters. (E) Differences in T-cell dysfunction score between different clusters. (*p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001).

Figure 4

Gene set enrichment analysis (GSEA) showed revealed enrichment for biological processes associated with immunosuppression.

Figure 5

Machine learning identifies hub genes associated with involved in the regulation of T-cell exhaustion. (A) Volcano map shows differential genes between cluster1 and cluster 2 in the TCGA. (B) 35 intersected differential genes associated with OS in TCGA and GSE72094 cohort. (C) The importance of 35 genes using random forest (RF). (D–E) 11 hub genes using the least absolute shrinkage and selection operator (LASSO) regression analysis. (F) 9 intersected hub genes based on 11 genes in LASSO and 15 top genes in RF.

Figure 6

The establishment and validation of the riskscore in TCGA and GSE72094 cohort respectively. (A, E) Risk map for prognostic signature and heat map for hub genes expression. (B, F) Kaplan–Meier curves for the riskscore. (C, G) The prognostic signature predict time-dependent ROC curves at 1, 2, and 3 years OS. (D, H) Multivariate cox regression analysis to verify the independent predictive value of the riskscore.

Figure 7

Correlation between the riskscore and clinical features. (A) age, (B) gender, (C) T stage, (D) N stage, (E) M stage, (F) pathological stage. (ns p > 0.05; *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001).

Figure 8

Construction and evaluation of the nomogram. (A) The nomogram combined the riskscore with stage for prognostic prediction of a patient with LUAD in the TCGA cohort. (B, E) Calibration curves of 1-year, 2-year, and 3-year OS for LUAD patients in the TCGA cohort and GSE72094 cohort. (C, F) Decision curve analysis of 3-year survival benefit in the TCGA cohort and GSE72094 cohort. (D, G) Time-dependent receiver operating characteristic (ROC) curves of the nomogram to predict 1-year, 2-year, and 3-year OS in the TCGA cohort and GSE72094 cohort.

Figure 9

The distribution of immune score (A), stromal score (B), estimate score (C), IPS score (D), TIDE score (E), TMB (F), and the mRNA expression of immune checkpoint inhibitors (G–I) in different risk groups.

Figure 10

(A) Differences in immunotherapy response between low- and high-risk groups in the GSE91061 dataset. (B) Progression-free survival for patient with anti-PD-1/PD-L1 therapy between the low- and high-risk groups in the GSE135222 dataset.

Figure 11

The IC50 values of chemotherapy and targeted drugs for LUAD in different risk groups. (*p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001).