Lung adenocarcinoma-specific three-integrin signature contributes to poor outcomes by metastasis and immune escape pathways

Data collection and processing

RNA-Seq Gene sequencing data of TCGA (LUAD and lung squamous cell carcinoma [LUSC], version 20160128) were downloaded from the University of California, Santa Cruz (UCSC) Xena Database,^[7] which was developed by UCSC computational genomics Laboratory. It included the mRNA expression profiles of 585 LUAD patients and 550 LUSC patients and had been normalized by Fragments Per Kilobase (FPKM) already. Through the “clusterProfiler” and “org.Hs.eg.db” packages in R, the Ensembl ID was converted to gene official symbols. The ITG subunits starting with the ITG of Homo sapiens were screened and 33 results were obtained (Table S1). Among these genes, three genes (“ITGAB2-AS1”, “ITGA6-AS1” and “ITGA9-AS1”) were noncoding RNAs and were excluded.

The survival data and phenotype information of these patients were also obtained from UCSC Xena Database. Overall survival (OS), which was defined as time from registration to death as a result of any cause, was the primary end point. Complete data of clinical pathology parameters and non-missing gene expression values were screened into subsequent analysis (Table 1). Missing values and ambiguous information, such as Tx, Nx, and Mx, were excluded when analyzed. Somatic mutations (SNPs and small INDELs) were also obtained from UCSC Xena database. The samples with any mutations in epidermal growth factor receptor (EGFR) or Kirsten rat sarcoma viral oncogene homolog (KRAS) or V-raf murine sarcoma viral oncogene homolog B1 (BRAF) or Tumor Protein P53 (TP53), regardless of the presence or not of a known driver of oncomutation, was identified as the sample with the corresponding specific gene mutation.

Microarray dataset GSE68465 used for verification was downloaded from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). This dataset contains gene expression and survival of LUAD patients from a multi-site, blinded validation study with complete clinicopathologic parameters of 439 cases.^[8] The “Series Matrix File(s)” were downloaded and the expression profile matrix and clinicopathologic parameters were extracted. According to the platform of GPL96, the probes corresponded to gene official symbols by “hgu133a.db” package in R. When it comes to one gene corresponding to different probes, the average value was taken to get the unique expression value of that gene in every sample.

Construction of ITG signature model based on COX regression analysis

The “glmnet” package in R was used to perform the COX regression analysis.^[9] To screen the prognostic factors for LUAD and LUSC, all ITG genes were subjected to univariate COX regression analysis and the genes with P < 0.05 were chosen for multivariate COX regression analysis. The genes with P < 0.05 were identified as the independent prognostic factors. Then, based on the weighted gene expression value of the independent prognostic factors selected above, the RS was calculated by the following formula:

“Exp” indicates the expression level of ITG and “β” indicates the regression coefficient obtained by multivariate COX regression analysis.

The best cutoff value of RS was calculated by X-tile software (version 3.6.1),^[10] and patients were divided into high-RS group and low-RS group according to the cutoff value. Kaplan–Meier method was performed to show the correlation between different groups and survival, and the log-rank test was used to analyze the differences between groups by the package “survival” in R.

Establishment and validation of nomogram prognostic prediction model

The nomogram prognostic prediction model involving important parameters of clinicopathologic characteristics was established, and calibration curves were used to evaluate the 1-, 3- or 5-year performance of this model by the “rms” package. C-index of this model was calculated with “boot” package, and receiver operating characteristic (ROC) analysis was also performed to estimate the accuracy of the nomogram for different time points’ survival prediction using the “survival ROC” package of R.

Biological function analysis

To explore the biological function, Gene Set Enrichment Analysis (GSEA) was performed with the GMT file (c2.KEGG.v6.2 and h.all.v7.1) gene set to obtain the biological processes downloaded from GSEA website (http://www. broad.mit.edu/gsea/) according to the manufacturer’s specifications.^[11] Normalized enrichment score (NES) >1.5 and P < 0.05 were defined as the significant enrichment pathways.

Immune infiltration analysis

TIMER (https://cistrome.shinyapps.io/timer/), a database designed for analyzing immune cell infiltrates in cancers,^[12] was used to estimate tumor immune infiltration in immune cells. Differences in selected ITGs’ (ITG subunit alpha 5 [ITGA5], ITG subunit alpha 6 [ITGA6], and ITG subunit alpha L [ITGAL]) expression levels in immune cells (CD4/CD8 T cells and dendritic cells) of LUAD were compared in the “Gene” section on the web page. The association between the expression of ITGAL and the infiltration degree of particular immune cell subsets was analyzed by TIMER in the “Correlation” section, without correlation-adjusted criteria.

Statistical analyses

All statistical analyses in the whole study were performed using R software (v3.6.3). Wilcoxon test was used to compare the expression level of genes between tumor and normal tissues. Relationship between RS and clinicopathologic parameters was analyzed by chi-square test with the package “stats” in R. Pearson correlation was utilized to identify the relationship between immunologic markers and 3-ITGs. P < 0.05 was defined as statistically significant.

Differential expression of ITGs in lung cancer patients

To investigate the role of ITG family genes in lung cancer, we analyzed the expression of all 30 ITGs in tumor tissues and normal tissues in LUAD and LUSC, respectively, using TCGA database. The expression of multiple ITGs was significantly different both in LUAD and LUSC (Figure 1A). In detail, ITG subunit alpha 2 (ITGA2), ITGA5, ITGA6, ITG subunit alpha 11 (ITGA11), ITG subunit alpha V (ITGAV), integrin subunit beta 1 binding protein 1 (ITGB1BP1), integrin subunit beta 1 binding protein 2 (ITGB1BP2), integrin subunit beta 3 binding protein (ITGB3BP), integrin subunit beta 4 (ITGB4), integrin subunit beta 5 (ITGB5) and integrin subunit beta 8 (ITGB8) were significantly more highly expressed in LUSC than in normal tissues (Figure 1B), while most ITGs were significantly expressed to a greater extent in normal tissues than LUAD, except for ITGA11, ITGA2, ITGAV, ITGB1BP1, ITGBIBP2, ITGB3BP3, ITGB4, ITGB5, integrin subunit beta 7 (ITGB7) and ITGB8 (Figure 1C), indicating that different subunits were involved in the development of particular lung cancer types.

Figure 1:

The expression pattern of integrin family genes of NSCLC patients in TCGA. (A) Heatmap of the differential expression of integrins between NSCLC and adjacent noncancerous tissues. Boxplot of the differential expression of integrins between (B) LUSC or (C) LUAD and adjacent normal tissues.

Construction of LUAD-specific 3-ITG prognostic RS signature

To identify the key ITGs affecting prognosis, we performed univariate and multivariate Cox regression analyses. By univariate Cox regression analysis, six prognosis-related ITG candidates, including ITGA2, ITGA5, ITGA6, ITGAL, ITGB1 and ITGB4, were identified as prognosis-related genes in LUAD, while ITGB6 was the only ITG subunit that could serve as a prognostic factor in LUSC (P < 0.05; Table 2).

The subsequent multivariate Cox regression analysis result showed that ITGA5 and ITGA6 were the prognostic risk factors, whereas ITGAL was a prognostic protective factor (Figure 2A). Then, according to the expression of ITGA5, ITGA6 and ITGAL, as well as their corresponding coefficients, the 3-ITG signature for prognostic prediction was constructed. The RS for LUAD patients was calculated as follows:

Figure 2:

The construction of 3-ITG prognostic signature. (A) Multivariate Cox regression analysis of six robust integrins. (B) The distribution of RS of patients in TCGA. (C) The boxplot of RS between alive and dead patients. (D) Heatmap for the expression of integrins in the 3-ITG signature. Kaplan–Meier (KM) curve of overall survival probability based on the (E) RS or (F) pathologic stages in LUAD.

The RS of each sample in TCGA-LUAD was calculated accordingly. After identification of 4.4 as the best cutoff value of RS using X-tile software according to the survival event probability, all samples were divided into a high-RS group and a low-RS group. Figure 2B illustrates the RS distribution of each patient; Figure 2C shows the RS distribution between those in the surviving group and in the death group.

Furthermore, the Kaplan–Meier analysis showed that patients in the high-RS group presented a remarkably shorter OS than those in the low-RS group (hazard ratio [HR] = 2.18, 95% confidence interval [CI] = 1.54–3.09, P < 0.001; Figure 2D). Besides this, the differential survival outcomes among various pStages showed that only pStage was a prognostic RF (Figures 2E and S1).

The prognostic prediction efficiency of ITGB6 was also assessed in LUSC, but both Kaplan–Meier analysis results (Figure S2A) and ROC curve (Figure S2D) criteria were unsatisfactory. Instead, the classical pStage (Figure S2B) and T stage (Figure S2C) were found as being suitable for prediction and better prognosis, indicating that ITGs were not the predominant molecules of LUSC.

Establishment of the prognostic predictive nomogram model based on 3-ITG signature in LUAD

To further evaluate the combination prognostic predictive efficiency of 3-ITG signature and pStage, a nomogram model based on these two predictive factors was established (Figure 3A) with a C-index of 0.704 (95% CI = 0.662-0.746) determined by internal verification. The calibration curve and the actual observation showed certain consistency at different time points (1, 3 and 5 years) (Figure 3B-D). Moreover, the predictive ability of prognostic value presented by ROC curves showed a satisfactory improvement of up to 4%–6% compared with classical grouping methods (Figure 3E-G), especially in terms of 1-year risk prediction. Taken together, these data suggested that the introduction of RS combined with the current pStage was a better indicator than the use of classic staging for predicting the OS of LUAD patients.

Figure 3:

Construction and identification of nomogram model 3-ITG signature combined with pathologic stages. (A) Nomogram for overall survival probability predicting. Calibration plots for (B) 1-year, (C) 3-year and (D) 5-year survival probabilities of the nomograms. Receiver operating characteristic (ROC) curves of nomogram model, 3-ITG signature, or pathologic stage individual for (E) 1-year, (F) 3-year and (G) 5-year survival probabilities by.

Stratification analysis of prognostic prediction of 3-ITG signature on clinicopathologic parameters

Although other clinicopathologic parameters were not included in the construction of the nomogram, the relationships between RS and these parameters were analyzed by chi-square testing. The results showed that a high RS was positively associated with T stage (P = 0.008) and N stage (P = 0.021; Table 3), suggesting that RS-based 3-ITG signature may result in poor prognosis by promoting proliferation and metastasis in LUAD.

Further survival analysis of clinicopathologic parameters’ stratification was also performed. The survival time of patients with high RS was significantly shorter than that of patients with low RS in Stage I LUAD (Figure 4A), while no significant difference was seen in Stage II (Figure 4B) or stages III and IV (Figure 4C). Similarly, the 3-ITG signature performed better in T2 (Figure 4E), N0 (Figure 4G), and M0 (Figure 4I) stages. These results suggested that use of the 3-ITG signature tended to be able to predict the prognosis accurately in the early stage.

Figure 4:

Kaplan–Meier curves of overall survival probability based on the RS in LUAD with different stages: (A) pStage I, (B) pStage II, (C) pStage III–IV, (D) T1 stage, (E) T2 stage, (F) T3–T4 stage, (G) N0 stage, (H) N1 stage, (I) M0 stage, (J) M1 stage.

Gene mutation plays a central role in diverse diseases.^[13] We also performed survival analysis on gene mutation stratification of TP53, KRAS, EGFR and BRAF. Independent of the absence or presence of KRAS mutation, the OS of high-RS patients was shorter than that of low-RS patients, indicating that no relationship existed between KRAS mutant status and ITG signature (Figure S3A, B); however, ITGs signature exhibited different prognostic prediction effects in patients with different mutant statuses of TP53, EGFR and BRAF. In the population of patients with TP53 mutant genes (Figure S3C), EGFR wild type (WT) (Figure S3F) and BRAF WT (Figure S3H), high-RS patients had a significantly worse prognosis than low-RS patients; whereas for patients with TP53 WT (Figure S3D), EGFR mutant (Figure S3E) or BRAF mutant (Figure S3G), the predictive function of 3-ITG signature disappeared, indicating that the 3-ITG signature played a critical role in prognostic prediction under the special background of specific tumor mutations.

Taken together, high RS according to 3-ITG signature was a risk prognostic factor, especially in patients with Stage I, demonstrating that 3-ITG RS signature might be an early optical biomarker for predicting prognosis of LUAD patients.

External validation in GEO datasets for LUAD-specific 3-ITG prognostic signature and nomogram model based on 3-ITG signature and pStage

To assess the prognostic prediction ability of the 3-ITG prognostic signature in LUAD, the GSE68465 dataset with microarray data of 439 patients with LUAD was selected as the validation dataset (Table 4). The RS of each LUAD sample based on the expression value of three ITGs (“ITGA5”, “ITGA6” and “ITGAL”) was calculated according to 3-ITG signature formula. All samples were divided into high-RS group and low-RS group according to the best cutoff value (Figure 5A). The expressions of ITGA5 and ITGA6 were upregulated in patients with a high score, while there tended to be low expression of ITGAL in LUAD patients with a high score (Figure 5B). The Kaplan–Meier curve further confirmed that survival of the high-RS group was significantly shorter than that of low-RS group (HR = 1.89, 95% CI = 1.35–2.64, P < 0.001; Figure 5C). Besides, the use of pStage was efficacious in prediction of the OS by Kaplan–Meier analysis (Figure 5D).

Figure 5:

External validation of nomogram model in GSE68465. (A) The distribution of RS. (B) The boxplot of RS between alive and dead patients. (C) Heatmap for the expression of three integrins. (D) Kaplan–Meier curve for 3-ITG signature (top) and pStage (bottom). The calibration plots of nomogram model for (E) 1-year, (F) 3-year, and (G) 5-year survival probabilities. ROC curves of nomogram model and pathologic stage or RS individually for (H) 1-year, (I) 3-year, and (J) 5-year survival probabilities.

Next, GSE68465 datasets were also utilized to validate the nomogram model based on 3-ITG signature and pStage. By external verification, C-index of this nomogram model was 0.671 (95% CI = 0.636–0.707) and the calibration curve almost overlapped the plot of real observations, in particular, it overlapped the 5-year prediction (Figure 5E-G). The ROC curve could be used to further evaluate the prediction capacity of survival using RS combined with pStage or individual data at different times. Prognostic prediction by the nomogram model of combination was more efficacious than each individual method for either 1, 3, or 5 years, with all scores improved to greater than 0.7 (Figure 5H-J). All these results indicated that the nomogram model based on 3-ITG signature and pStage could predict the prognosis of LUAD patients well, which is an applicable supplement for classical staging methods.

Biological pathway enrichment of patients with high RS in LUAD

To reveal the possible mechanisms of high RS according to 3-ITG signature associated with a poor prognosis, we performed GSEA analysis using the KEGG and HALLMARK pathway sets on LUAD both in TCGA and GSE68465 datasets (Table S2). Compared with the low-RS LUAD group, the high-RS group exhibited comprehensive enrichment of 11 pathways (NES > 1.5, P < 0.05) in both datasets of patients with high RS, and Figure 6A and B shows six pathways in which the NES was over 2.0 in one of the datasets, including the pathways related to cell adhesion and migration capability (such as “ECM_RECEPTOR_INTERACTION”, “APICAL_JUNCTION” and “EPITHELIAL_MESENCHYMAL_TRANSITION”), resistance in hypoxic conditions (“HYPOXIA”), proliferation and differentiation (“TGF_BETA_SIGNALING”), angiogenesis (“ANGIOGENESIS”), and others (“SMALL_CELL_LUNG_CANCER”, “PATHWAYS_IN_CANCER” and “ARRHYTHMOGENIC_RIGHT_VENTRICULAR_CARDIOMYOPATHY_ARVC”). Therefore, the complex mechanisms for poor prognosis of high-RS group strongly suggested differential function and the combined effect of three individual ITGs in the signature.

Figure 6:

Functional enrichment analysis by GSEA and identification of individual gene function. Functional enrichment pathways in patients with high risk score in (A) TCGA and (B) GSE68465. Kaplan–Meier curve for verification of survival by the expression of (C) ITGA5, (D) ITGA6, and (E) ITGAL in LUAD patients in TCGA. Functional enrichment pathways in patients with high expression of (F) ITGA5 in TCGA, (G) ITGA5 in GSE68465, (H) ITGAL in TCGA, and (I) ITGAL in GSE68465.

We next ascertained the potential function of three ITGs in malignant biological properties, respectively. Firstly, all three ITGs were found to be prognostic for OS both in TCGA (Figure 6C-E) and GSE68465 (Figure S4A) to some extent. Besides, patients with a high expression of ITGA5 or ITGA6 had a shorter OS, while patients with a high expression of ITGAL had better prognosis. According to the best cutoff value of each molecule, the datasets were further divided into a high-expression group and a low-expression group, where upon GSEA was performed. The results showed that 14 pathways were enriched in high-ITGA5 expression groups (Figure 6F, G) and 33 pathways in high-ITGAL expression groups (Figure 6H, I), while only one common pathway was enriched in high-ITGA6 expression groups (Figure S4B, C). Surprisingly, the enrichment results exhibited distinct biological characteristics. Pathways enriched in high-ITGA5 expression group were mostly related to the function of cell migration and invasion facilitating metastasis (i.e., “ECM_RECEPTOR_INTERACTION”, “APICAL_JUNCTION”, “EPITHELIAL_MESENCHYMAL_TRANSITION”, “HYPOXIA”, “FOCAL_ADHESION”, “REGULATION_OF_ACTIN_CYTOSKELETON”), while pathways enriched in high-ITGAL expression group were mostly associated with immune cell recognition and killing (i.e., “T_CELL_RECEPTOR_SIGNALING_PATHWAY”, “COMPLEMENT”, “IL2_STAT5_SIGNALING”, “CHEMOKINE_SIGNALING_PATHWAY”, “INFLAMMATORY_RESPONSE”, “NATURAL_KILLER_CELL_MEDIATED_CYTOTOXICITY”). These results indicated that high-ITGA5 expression groups had a poorer prognosis owing to the increased ability of metastasis, while an elevated activation of inflammatory/immune cells might be the reason for those in high-ITGAL groups having a better prognosis.

Recently, accumulating evidence demonstrates that ITGs are emerging as important regulators of immune cell signaling.^[14,15] Thus, it is important to identify which subunits of ITG correlated with poor prognosis were mainly involved in the immune mechanism of LUAD. By analyzing the expression level of three ITGs of 3-ITG in immune cells associated with the recognition of tumor, it was found that only ITGAL had a strong correlation with a correlation coefficient (Cor) >0.5 and P < 0.05 (Figure 7A–C). This finding was consistent with the results of GSEA analysis in high-ITGAL groups. Furthermore, we ascertained the relevance between expression of ITGAL and representative biomarkers of immune cells (Table 5). ITGAL showed a strong positive correlation with genes related to immune surveillance and recognition (Figure 7D), especially with CD3E and CD2 of CD4⁺ T cells (Cor > 0.8).

Figure 7:

Comprehensive analysis of tumor-infiltrating immune cells in LUAD. Expression levels of (A) ITGA5, (B) ITGA6, and (C) ITGAL in CD4⁺ T cells, CD8⁺ T cells, and dendritic cells. (D) Scatterplots of correlations between ITGAL expression and gene markers of T cells.

Taken together, these results suggest that in the 3-ITG signature, ITGA5 and ITGA6 mostly govern cell adhesion (as reported) and promote proliferation or phenotypic transformation of LUAD, while reduction of ITGAL helps the tumor cells to escape immune control. Finally, three ITGs operate together and lead to the poor prognosis of patients with high RS according to 3-ITG signature.