Integrative bioinformatics approaches to establish potential prognostic immune-related genes signature and drugs in the non-small cell lung cancer microenvironment

Abstract

Introduction: Research has revealed that the tumor microenvironment (TME) is associated with the progression of malignancy. The combination of meaningful prognostic biomarkers related to the TME is expected to be a reliable direction for improving the diagnosis and treatment of non-small cell lung cancer (NSCLC). Method and Result: Therefore, to better understand the connection between the TME and survival outcomes of NSCLC, we used the "DESeq2" R package to mine the differentially expressed genes (DEGs) of two groups of NSCLC samples according to the optimal cutoff value of the immune score through the ESTIMATE algorithm. A total of 978 up-DEGs and 828 down-DEGs were eventually identified. A fifteen-gene prognostic signature was established via LASSO and Cox regression analysis and further divided the patients into two risk sets. The survival outcome of high-risk patients was significantly worse than that of low-risk patients in both the TCGA and two external validation sets (p-value < 0.05). The gene signature showed high predictive accuracy in TCGA (1-year area under the time-dependent ROC curve (AUC) = 0.722, 2-year AUC = 0.708, 3-year AUC = 0.686). The nomogram comprised of the risk score and related clinicopathological information was constructed, and calibration plots and ROC curves were applied, KEGG and GSEA analyses showed that the epithelial-mesenchymal transition (EMT) pathway, E2F target pathway and immune-associated pathway were mainly involved in the high-risk group. Further somatic mutation and immune analyses were conducted to compare the differences between the two groups. Drug sensitivity provides a potential treatment basis for clinical treatment. Finally, EREG and ADH1C were selected as the key prognostic genes of the two overlapping results from PPI and multiple Cox analyses. They were verified by comparing the mRNA expression in cell lines and protein expression in the HPA database, and clinical validation further confirmed the effectiveness of key genes. Conclusion: In conclusion, we obtained an immune-related fifteen-gene prognostic signature and potential mechanism and sensitive drugs underling the prognosis model, which may provide accurate prognosis prediction and available strategies for NSCLC.

Keywords: drug sensitivity; estimate; non-small cell lung cancer; prognostic gene signature; tumor microenvironment.

FIGURE 1

Flow chart of data collection and analysis.

FIGURE 2

Relationship between clinical characteristics and immune, stromal and ESTIMATE scores (A–C). The optimal cutoff values of the immune, stromal and ESTIMATE scores. (D–F). (K–M) analysis of immune, stromal and ESTIMATE scores. (G–I). Distribution of the three scores among patients with different statuses (J–L). Distribution of immune, stromal and ESTIMATE scores among NSCLC stages (M–O). Distribution of three scores between M-stage of NSCLC.

FIGURE 3

Heatmap, volcano plot and enrichment analysis of GO and KEGG for DEGs (A). Heatmap of DEGs in TCGA (B). Volcano plot of DEGs in TCGA (C). Top 5 enriched biological processes, molecular functions, and cellular components of upregulated co-DEGs and (D) downregulated co-DEGs (E). Top 10 KEGG pathways of upregulated co-DEGs and (F) downregulated co-DEGs.

FIGURE 4

Development of the prognostic signature in the TCGA cohort (A). Diagnostic model construction using a LASSO regression model (B). Coefficient distribution plots to select the optimum lambda value (C). Results of multivariate Cox regression analysis of OS in the TCGA cohort.

FIGURE 5

Prognostic value of the 15-gene prognostic model in the TCGA and validation cohorts (A–C). Heatmap of fifteen genes between the two groups in the TCGA and validation sets (D–F). Risk score scatter plot. Red dots indicate dead patients, and blue dots indicate alive patients (G–I). Risk score curve plot. The dotted line indicates the individual distribution of the risk score, and the patients are categorized into low-risk (blue) and high-risk (red) groups (J–L). Survival status and time of patients between the two groups in the TCGA and validation sets, respectively (M–O). The time-dependent ROC curve of patients between the two groups in the TCGA and validation sets.

FIGURE 6

Univariate and multivariate Cox regression analyses in the TCGA and validation cohorts and establishment of the nomogram (A, B). Univariate and multivariate Cox regression analyses in TCGA (C, D). Univariate and multivariate Cox regression analyses in GSE31210 (E, F). Univariate, multivariate Cox regression analysis in GSE37745 (G). Establishment of a nomogram predicting OS based on the independent prognostic factors in TCGA (H). ROC curve of the nomogram, risk score and other relevant clinical parameters in TCGA (I–K). Calibration curves of the nomogram prediction of 1-, 3-, and 5-year survival in TCGA.

FIGURE 7

(A) KEGG pathways in high-risk group (B). Top 5 Gene Set Enrichment Analysis of gene set of high-risk in TCGA cohort.

FIGURE 8

(A–D) Immune score, stroma score, ESTIMATE score and tumor purity of the high- and low-risk groups (E). Distribution of infiltration of 22 immune cell types in the two risk groups (F) The expression of immune checkpoint genes between the two cohorts (G). The proportions of different immune cells in the high- and low-risk groups (H, I). Mutated gene mutation profiles between the high- and low-risk groups.

FIGURE 9

Drugs sensitivity in the high and low-risk group.

FIGURE 10

Screening and validation of hub genes (A). PPI network among prognostic genes. (B, C) Two modules contained many genes in the PPI network (D). The relationship between survival state, stage and mRNA expression levels of EREG, Kaplan–Meier curves of EREG in OS (E). The relationship between survival state, stage and mRNA expression levels of ADH1C, Kaplan–Meier curves of ADH1C in OS (F). ADH1C protein levels in normal lung and NSCLC were visualized by IHC in HPA. (G, H) Quantitative real-time PCR analysis of the mRNA expression levels of EREG and ADH1C in NSCLC cell lines and normal lung epithelial cells.