Construction of prognostic IRGPs signature
After removing IRGPs with score of 0 or 1 in more than 80% of LUAD cases, a total of 12334 IRGPs were pared. After screened by univariate Cox regression analysis, 54 IRGPs were subjected to LASSO regression analysis and 21 IRGPs were filter out (Fig. 1a). Then, 8 IRGPs were identified and used to develop a prognostic IRGPs signature with multivariate Cox regression analysis (Fig. 1b, Table 2). Risk score = (1.139* Score BIRC5|BPHL) + (-0.658* Score CCL2|OAS1) + (-0.461* Score CD19|PI3) + (-0.557* Score CD3G|IL7) + (0.723* Score DKK1|IKBKB) + (0.448* Score F2RL1|LTB) + (-0.428* Score PIK3CD|S100A2) + (-0.606* Score SERPIND1|VEGFC).
Table 2
Information on the 8 immune-related gene pairs (IRGPs).
IRG 1 | Immune processes | IRG 2 | Immune processes | Coefficient |
BIRC5 | Antimicrobials | BPHL | Antimicrobials | 1.139 |
CCL2 | Antimicrobials | OAS1 | Antimicrobials | -0.658 |
CD19 | BCR Signaling Pathway | PI3 | BCR Signaling Pathway | -0.461 |
CD3G | TCR signaling Pathway | IL7 | Cytokines | -0.557 |
DKK1 | Cytokines | IKBKB | TCR signaling Pathway | 0.723 |
F2RL1 | Antimicrobials | LTB | Cytokines | 0.448 |
PIK3CD | BCR Signaling Pathway | S100A2 | Antimicrobials | -0.428 |
SERPIND1 | Antimicrobials | VEGFC | Cytokines | -0.606 |
Abbreviations: IRGPs, immune-related gene pairs; IRG, immune-related gene. |
Validation and evaluation of prognostic IRGPs signature
According to cut-off risk score (0.99), patients both in TCGA set and GEO set were divided into low- and high- risk group. Moreover, we calculated the AUC for predicting 1- and 3- year OS to evaluate the predictive capacity of this IRGPs signature. The AUC for 1- and 3- year OS in TCGA set was 0.867 and 0.870, respectively (Fig. 2a-b). And, in GEO set was 0.819 and 0.803, respectively (Fig. 2c-d). In addition, the AUC in Total set was 0.845 and 0.801, respectively (Fig. 2e-f).
The IRGPs signature is an independent prognostic factor of survival
We performed survival analysis to compare the survival difference between low- and high- risk group. All of Kaplan-Meier plot in three sets demonstrated that high-risk LUAD patients exhibited poorer prognosis than low-risk LUAD patients (TCGA set: p < 0.001, Fig. 3a; GEO set: p < 0.001, Fig. 3b; Total set: p < 0.001, Fig. 3c). Furthermore, stratification analyses showed the clinical outcome of high-risk LUAD patients in each stratum of age, gender, TNM stage, tumor size, lymph node metastasis and distance metastasis was poorer than that of low-risk patients (Additional file 1: Fig. S1).
Then, we took advantage of univariate and multivariate Cox regression model to compare the immune risk score with clinical parameters (age, gender, smoking, histologic grade, TNM grade, tumor size, lymph node metastasis and distance metastasis). Univariable Cox regression analysis indicated that immune risk score was an important factor of patients’ prognosis (TCGA set: HR = 4.227, 95%CI [2.952, 6.053], p < 0.001, Fig. 3d; GEO set: HR = 2.484, 95%CI [1.911, 3.229], p < 0.001, Fig. 3f; Total set: HR = 2.975, 95%CI [2.416, 3.662], p < 0.001, Fig. 3h). Moreover, multivariable Cox regression demonstrated that immune risk score was independent prognostic factors (TCGA set: HR = 3.856, 95%CI [2.621, 5.673], p < 0.001, Fig. 5d; GEO set: HR = 2.473, 95%CI [1.789, 3.436], p < 0.001, Fig. 5g; Total set: HR = 2.876, 95%CI [2.264, 3.654], p < 0.001, Fig. 5i).
Correlation between IRGPs signature and clinical characteristics
As shown in Fig. 4a, between low- and high- risk group, the distribution of gender (p = 0.001), histologic grade (p < 0.001), TNM grade (p < 0.001), tumor size (p < 0.001), lymph node metastasis (p < 0.001) and KRS mutation (0.012) was significantly different. Meanwhile, compared with female, the immune risk score in male was significantly increased (p = 0.001, Fig. 4b). Similar result was observed in patients with lymph node metastasis (p < 0.001, Fig. 4f) and distance metastasis (p = 0.024, Fig. 4g). In addition, with the increase of TNM grade (p < 0.001, Fig. 4c) and tumor size (p < 0.001, Fig. 4e), the immune risk score was also increased. The immune risk score was significantly different among three histologic grades (poor differentiation, Moderate differentiation and well differentiation) (p < 0.001, Fig. 4d).
Relationship between IRGPs and TIICs
A total of 823 LUAD cases met the cut-off value: CIBERSORT P < 0.05. Macrophages (M0, M1, M2) (32.98%) was the most abundant immune cell, followed by Plasma cells (17.23%) and resting memory CD4 T cells (9.70%) (Additional file 1: Fig. S2a). Compared with normal lung tissues, LUAD tissues had more proportions of naive B cells (p = 0.001), memory B cells (p < 0.001), Plasma cells (p < 0.001), activated memory CD4 T cells (p < 0.001), helper follicular T cells (p < 0.001), regulatory T cells (Tregs) (p < 0.001), gamma delta T cells (p = 0.011), Macrophages M1 (p < 0.001) and less proportions of T resting memory cells CD4 (p = 0.001), resting NK cells (p < 0.001), Monocytes (p < 0.001), resting Mast cells (p = 0.004), Eosinophils (p < 0.001) (Additional file 1: Fig. S2b). In addition, the proportions of Macrophages M0 (p < 0.001), Macrophages M1 (p = 0.001), activated memory CD4 T cells (p < 0.001) were significantly increased in high- risk group, whereas, the proportions of memory B cells (p < 0.001), Plasma cells (p = 0.008), naive T cells CD4 (p = 0.010), gamma delta T cells (p = 0.012), resting NK cells (p < 0.001), Monocytes (p < 0.001), resting Dendritic cells (p < 0.001), resting Mast cells (p < 0.001) and Neutrophils (p = 0.002) were significantly decreased (Fig. 5a). Furthermore, survival analysis demonstrated that low level of Monocytes (p = 0.020, Fig. 5b), Plasma cells (p = 0.020, Fig. 5d) and memory B cells (p = 0.010, Fig. 5e) as well as high level of Macrophages M1 (p = 0.048, Fig. 5c) were highly related with poor prognosis of LUAD patients.
Expression profile of immunomodulators.
In the present study, we quantified 11 immunomodulators (CTLA4, ICOS, ICOSLG, IFN- γ, LAG3, NKG2A, PD − 1, PD − L1, TIGIT, TIM3 and VISTA). The expression of CTLA4 (p < 0.001), IFN- γ (p = 0.002), LAG3 (p < 0.001), PD − 1 (p < 0.001) and TIGIT (p < 0.001) was significantly up-regulated in LUAD tissues compared with that in normal lung tissues, whereas, ICOS (p < 0.001), PD − L1 (p < 0.001), TIM3 (p < 0.001) and VISTA (p < 0.001) was significantly down-regulated (Additional file 1: Fig. S3a). Moreover, in high- risk group, 8 immunomodulators (CTLA4 (p < 0.001), ICOS (p < 0.001), IFN- γ (p = 0.032), LAG3 (p = 0.003), PD − 1 (p < 0.001), TIGIT (p < 0.001), TIM3 (p = 0.036) and VISTA(p < 0.001)) were significantly decreased (Additional file 1: Fig. S3b).
GSEA analysis
To explore the basic biological mechanisms of the IRGPs signature, we carried out GSEA analysis. A total of 27 KEGG pathways were enriched between high- risk and low- risk group, among which, various KEGG pathways were highly related with immune system, such as “Chemokine signaling pathway”, “Intestinal immune network for IgA production”, “T cell receptor signaling pathway”, “B cell receptor signaling pathway”, “Primary immunodeficiency”, “Leukocyte transendothelial migration”, “MAPK signaling pathway” and “PPAR signaling pathway” (Fig. 6).