Characteristics of hypoxic tumor microenvironment in non-small cell lung cancer, involving molecular patterns and prognostic signature

Abstract

Background: The mechanisms of hypoxia or immune microenvironment in cancer have been studied respectively, but the role of hypoxia immune microenvironment in non-small cell lung cancer (NSCLC) still needs further exploration.

Methods: By applying the K-means algorithm, 1,121 patients with NSCLC were divided into three categories. We evaluated the constructed signature in order to link it with the prognosis, which was constructed by univariate and least absolute shrinkage operator (LASSO) Cox regression analysis.

Results: A total of three clusters were obtained by clustering five Gene Expression Omnibus (GEO) data sets. Gene Set Variation Analysis (GSVA) and immune infiltration analysis were performed to explore the biological behavior. Cluster one presented an activated state of oncogenic pathways, and compared with the other two clusters, the median risk score was the highest, which was the reason for its poor survival. Cluster three showed that the immune pathway was active and the median risk score was the lowest, so the survival was the best. However, cluster two presented a state in which both immune and matrix pathways were activate. This was manifested as mutual antagonism, and its risk score was in the middle. Its survival was in the middle.

Conclusions: This work revealed the role of hypoxia related genes (HRGs) modification in tumor microenvironment, which was conducive to our comprehensive analysis of the prognosis of NSCLC, and provided direction and guidance for clinical immunotherapy.

Keywords: Gene Set Variation Analysis (GSVA); Hypoxia related genes (HRGs); immune infiltration; prognosis.

Figure 1

Cluster analysis of HRGs. (A) Unsupervised cluster analysis of 76 HRGs from cBioPortal database, consensus matrices for k =3. (B) The survival analysis of the three HRG modification patterns involved 1,121 patients, of which cluster one included 353 cases (31.5%), cluster two included 409 cases (36.5%), and cluster three included 358 cases (32%). (C) The principal component analysis of the transcriptome profile of the three HRG modification patterns aimed to show the reliability of clustering. (D) In the heat map, cluster, status, stage, gender and age were regarded as patient annotations. Yellow represented high expression, blue represented low expression. HRG, hypoxia related gene; CDF, cumulative distribution function.

Figure 2

Immune infiltration analysis and GSVA Enrichment Analysis. (A) The abundance of each TME infiltrating cell among the clusters, involving twenty-eight kinds of immune cells. The box line of the box plot represented the median value, the black dots represented the outliers, and the asterisk represented P value (***, P<0.001). (B) Expression differences of three different modification patterns among twenty-eight immune cells. Yellow represented high expression, blue represented low expression. (C) Activation and inhibition of different hypoxia modification patterns in biological pathways. Yellow represented high expression, blue represented low expression. (D) Differences of three hypoxia modification patterns in immune microenvironment. Kruskal-Wallis test was used to compare the statistical differences among the three clusters. GSVA, Gene Set Variation Analysis; TME, tumor microenvironment.

Figure 3

Selection of differential genes and analysis of GO enrichments. (A) The differences in matrix activation pathways of the three hypoxia modification patterns. EMT, epithelial-mesenchymal transition. (B) The differences in costimulatory molecules of the three hypoxia modification patterns. (C) The condition for screening 1,964 common differential genes was P<0.001. (D) The GO enrichment was designed to functionally annotate hypoxia-related patterns. (*, P<0.05; ***, P<0.001). GO, gene ontology; EMT, epithelial-mesenchymal transition.

Figure 4

Construction of prognostic signature. (A) The Kaplan-Meier curve aimed to analyze the survival significance of high and low risk groups. (B) The receiver operating characteristic curve aimed to reflect the predictive efficiency of signature. (C) Forest plot aimed to confirm that the signature was an independent prognostic factor for NSCLC. (D,E,F,G,H) Kaplan-Meier curve was used to analyze the meaning of signature in each data set (GSE30219, GSE37745, GSE41271, GSE42127, GSE50081).

Figure 5

Interrelationships among signature factors. (A) Forest plot was used to display the prognostic significance of factors in the signature. (B) Analysis of the differences of signature factors in three hypoxia modification patterns. (C) Alluvial diagram showed the connection and changes of status, stage, cluster and riskScore. (D) Interaction among signature factors. The size of circle represented the influence of each signature on the prognosis. The green dots in the circle represented prognostic risk factors, while the purple dots in the circle represented prognostic favorable factors. The lines connecting the signature factors showed their interaction, and the thickness indicated the mutual strength. Negative correlation was blue, while positive correlation was red. The signature factors were divided into three groups, marked with red, yellow and blue. (***, P<0.001).

Figure 6

Significance of signature factors. (A) The relationship among riskScore and signature factors were analyzed by spearman. (B) The Kruskal-Wallis test was used to compare the differences among three different hypoxia modification patterns and riskScore. (C) Analysis of the signature in the matrix microenvironment, and the asterisk represented P value (*, P<0.05; **, P<0.01; ***, P<0.001). (D) Difference analysis of costimulatory molecules in high and low risk groups. (E) Analysis of the difference between high and low risk groups in immune infiltration. (F,G) Waterfall plot of tumor somatic mutations involving signature factors in lung adenocarcinoma (LUAD) and lung squamous cell carcinoma (LUSC). (F) LUAD (G) LUSC.

Figure 7

Validation of the signature in the TCGA cohort. (A) The Kaplan-Meier curve was designed to verify the prognostic significance of high- and low-risk groups in the TCGA cohort. (B) The forest plot was used to verify that the signature was indeed an independent prognostic factor in the TCGA cohort. (C) The risk plot was used to show that the survival status became worse as the risk value increased. (D) The nomogram contained age, gender, stage, T, N, signature. The x-axis of the calibration chart was the predicted recurrence probability result, and the y-axis was the actual recurrence probability. ROC analysis detected the accuracy of prediction and inspection. TCGA, The Cancer Genome Atlas; ROC, receiver operating characteristic.

Figure 8

Efficacy of the signature in immunotherapy. (A) Smoking correlation analysis of signature factors. (B) Significance of signature factors follow-up treatment. (C) Significance of signature factors primary therapy outcome. (D) Significance of signature factors in immunotherapy in melanoma. (E) Significance of signature factors in immunotherapy in renal cell carcinoma.