Immunoscore Predicts Survival in Early-Stage Lung Adenocarcinoma Patients

Background: The lung cancer staging system is insufficient for a comprehensive evaluation of patient prognosis. We constructed a novel immunoscore model to predict patients with high risk and poor survival.

Method: Immunoscore was developed based on z-score transformed enrichment score of 11 immune-related gene sets of 109 immune risk genes. The immunoscore model was trained in lung adenocarcinoma cohort from The Cancer Genome Atlas (TCGA-LUAD) (n = 400), and validated in other two independent cohorts from Gene Expression Omnibus (GEO), GSE31210 (n = 219) and GSE68465 (n = 356). Meta-set (n = 975) was formed by combining all training and testing sets.

Result: High immunoscore conferred worse prognosis in all sets. It was an independent prognostic factors in multivariate Cox analysis in training, testing and meta-set [hazard ratio (HR) = 2.96 (2.24–3.9), P < 0.001 in training set; HR = 1.99 (1.21–3.26), P = 0.006 in testing set 1; HR = 1.48 (1.69–2.39), P = 0.005 in testing set 2; HR = 2.01 (1.69–2.39), P < 0.001 in meta-set]. Immunoscore-clinical prognostic signature (ICPS) was developed by integrating immunoscore and clinical characteristic, and had higher C-index than immunoscore or stage alone in all sets [0.72 (ICPS) vs. 0.7 (immunoscore) or 0.59 (stage) in training set; 0.75 vs. 0.72 or 0.7 in testing set 1; 0.65 vs. 0.61 or 0.62 in testing set 2; 0.7 vs. 0.66 or 0.64 in meta-set]. Genome analysis revealed that immunoscore was positively correlated with tumor mutation burden (R = 0.22, P < 0.001). Besides, high immunoscore was correlated with high proportion of carcinoma-associated fibroblasts (R = 0.32, P < 0.001) in tumor microenvironment but fewer CD8+ cells infiltration (R = −0.28, P < 0.001).

Conclusion: The immunoscore and ICPS are potential biomarkers for evaluating patient survival. Further investigations are required to validate and improve their prediction accuracy.

Introduction

Lung cancer ranks the top of cancer-related death worldwide (1). Histologically, 15 percent of patients are categorized as small cell lung cancer (SCLC) while the other 85% as non-small cell lung cancer (NSCLC) (2). Among NSCLC, lung adenocarcinoma (LUAD) is the most common subtype (3). Surgical resection remains to be the standard clinical practice for patients with early-stage LUAD (4), and the 5-year survival rate is about 60% (5). Platin-based adjuvant chemotherapy has demonstrated the improvement of 5-year survival for stage II–IIIA patients for about 5%, at the price of chemotherapy-induced toxicity (6, 7). Adjuvant immunotherapy with immune checkpoint inhibitors has come into several clinical trials, but no definitive effectiveness made so far (8). Although using the American Joint Committee on Cancer (AJCC) TNM staging system improves prognostic prediction, it is still inconclusive due to other unknown factors. Thus, the development of new biomarkers is imperative for stratifying risk and optimizing treatment for lung cancer patients with early stage.

Tumor immune microenvironment (TIM) has long been recognized as a crucial factor in cancer progression and metastasis (9). Several studies have explored the TIM as a prognostic biomarker in lung cancer (10). For example, Brambilla et al. found higher CD4+/CD8+ ratio conferred better survival in patients with NSCLC (11). Also, for cancer cell itself, programmed cell death protein ligand 1 (PD-L1) expression and tumor mutation burden (TMB) have been used to predict outcome in NSCLC patients. Several investigations have indicated that patients with high TMB or high PD-L1 expression were associated with poor survival in resected NSCLC patients and might benefit from adjuvant chemotherapy (12, 13). However, substantial patients with low PD-L1 expression and low TMB still have poor outcomes. Therefore, exploring additional prognostic markers based on TIM could benefit larger population (14).

In our research, we developed novel prognostic early-stage lung cancer immunoscore model by integrating enrichment score of 11 immune gene sets using ssGSEA algorithm. ssGSEA algorithm was based on gene ranks in and out of the selected gene set (15). To date several signatures used for phenotype classification or survival prediction have been developed by leveraging this algorithm (16–18). After immunoscore model construction in the training set, we evaluated its prognostic abilities in training, testing and meta-set. Moreover, we built Immunoscore-clinical prognostic signature (ICPS) by incorporating both immunoscore and clinical factors.

Materials and Methods

Clinical Data Processing

We used three largest publically available datasets, TCGA-LUAD, GSE31210, and GSE68465, deposited in Genomic Data Commons (GDC) portal (https://portal.gdc.cancer.gov) or Gene Expression Omnibus (GEO) website (https://www.ncbi.nlm.nih.gov/geo) (19–23). Clinical and pathological information regarding to TCGA-LUAD cohort were retrieved from cBioportal website (https://www.cbioportal.org) with “cdgsr” package (24–26), whereas information related to GSE31210 and GSE68465 were obtained through “GEOquery” package (27). Samples without overall survival (OS) information or with OS time of 0 were excluded. We also ruled out samples with documented neoadjuvant therapies to reduce potential confounding bias. TNM stage were used and transformed to AJCC staging groups. Samples with specific T subcategories (like T2a or T2b) were converted to staging groups according to AJCC 7th edition. T1N0, T2N0, T1N1, T2N1 or T3N0 were converted to stage 1A, 1B, 2A, 2B, respectively, conforming to AJCC 6th edition. For GSE31210 without TNM stage information, we used the pathological stage in the clinical file directly.

RNA-seq and Microarray Data Preprocessing

Raw “.CEL” files of microarray data were downloaded from GEO website and read by “affy” package with the latest brainarray CDF files (October 2019, version 24) (28, 29). Robust multi-array average (RMA) algorithm in “affy” package was then applied to normalize gene expression intensity (28, 30). RMA algorithm included background adjustment, quantile normalization, and measurement summation when multiple probes were used to quantify the same gene expression intensity. After normalization, “arrayQualityMetrics” package was utilized to detect and exclude possible outliers (31). For RNA-seq data, level 3 FPKM data were downloaded using TCGAbiolink R package (32). FPKM values were then transformed into TPM values, which allowed a more direct comparison between samples as the sum of all TPMs in each sample were the same. As a result, the inflated statistical significance was reduced (33). TPM values were subsequently log 2 transformed to fit a more normal distribution. Entrez IDs were used across all platforms. Only samples with clinical information were retained. Finally, TCGA-LUAD cohort was used as the training set for immunoscore model construction, which contained 400 patients with RNA-seq data and survival information. Two microarray datasets, GSE31210 (n = 219) from Affymetrix Human Genome U133 Plus 2.0 Array platform as testing set 1, and GSE68485 (n = 356) from Affymetrix Human Genome U133A Array platform as testing set 2, were used to assess the immunoscore performance in predicting survival of early-stage LUAD patients.

Immunoscore Construction

We searched Immport database (https://immport.niaid.nih.gov) and downloaded 1811 immune-related genes from 17 categories (18). Of 1,811 immune-related genes, 1,361 of them were contained in the training set. Univariate Cox proportional regression analysis was used to investigate their associations with patient survival using “survival” package (34). Only the genes with P-value < 0.05 and hazard ratio (HR) > 1 were screened out as immune risk genes for further study. We then implemented single sample gene set enrichment analysis (ssGSEA) algorithm to quantify the enrichment score of immune risk genes in various immune-related gene set using “GSVA” package (35). Difference of enrichment statistic of genes in the gene set and outside were computed, and normalized to fit a relatively uniform scale as Barbie et al. described (15, 35). We then transformed the normalized enrichment score into Z-score to conform standard normal distribution using the following algorithm:

$\begin{array}{l}ZNE{S}_{ij}=\frac{NE{S}_{ij}-M_j}{S{D}_j}\end{array}$

The final Z-score transformed normalized enrichment score of sample i, immune gene set j was denoted by ZNES_ij. The normalized enrichment score of sample i, immune gene set j was denoted by NES_ij. The mean and standard deviation (SD) of enrichment score across all samples in immune gene set j were denoted by M_j and SD_j, respectively. This transformation obtained a uniform underlying distribution (mean = 0, standard deviation = 1) of each gene set across various platform; Immunoscore model was established by integrating all Z-score transformed normalized enrichment score using regularized Cox regression with the ridge penalty.

$\begin{array}{l}{\text{Immunoscore}}_{\text{i}}={\sum }_{\text{j}=1}^n{\beta }_{\text{j}}*{\text{ZNES}}_{\text{ij}}\end{array}$

Immunoscore of sample i was denoted by Immunoscore_i. Ridge Cox regression coefficient of gene set j was denoted by β_j and standard normal distribution transformed normalized enrichment score of sample i, immune gene set j was denoted by SNES_ij. Ridge regression was used to address the possible collinearity (i.e., the correlated immune gene sets) to prevent overfitting (36). It was conducted by “glmnet” package and the tuning parameter Lambda was chosen with minimum criteria (37). Thus, a new variable immunoscore was created to predict patient survival. It could also be served as the quantitative measurement of hazardous level of tumor immune microenvironment with its biological background.

Validation of Immunoscore

After immunoscore development, we applied the same formula to two independent testing sets, GSE31210 and GSE68485. Meta-set was formed by combining all training and testing sets. Univariate and Multivariate regression were used to evaluate the predictive power of the immunoscore model in all training, testing and meta-set. Age, stage, gender and smoking history were included in multivariable Cox analysis. Fraction of genome alteration in TCGA-LUAD clinical profile was also involved as a covariate in the TCGA-LUAD cohort. Patients were divided into high-immunoscore and low-immunoscore subgroups based on median value in the training set. Patients with immunoscore higher than cut-off value were assigned to high-immunoscore subgroup, while with immunoscore lower or equal to cut-off value were assigned to low-immunoscore subgroup. Kaplan-Meier analysis was performed to these two groups. Time-dependent receiver operator characteristic (ROC) curve analysis was utilized to assess the predictive accuracy for early-stage LUAD patients using “timeROC” package (38). The prognostic value of immunoscore in various treatment groups was evaluated in GSE68465, which contained detailed information of whether patients received adjuvant chemotherapy or radiotherapy with sufficient sample size in each category (75 patients with documented adjuvant therapy, 271 patients without documented adjucvant therapy).

Comparison With Other Gene Expression Signatures

The immunoscore was compared with other existing NSCLC prognostic signature to assess its clinical utility. To date, numerous gene expression signatures have been developed. We selected two immune-related signatures (39, 40), one glycolysis-based signature. In addition, another malignancy gene signature was included, which had the top-ranked prognostic capability when compared with random signature in lung adenocarcinoma patient (41, 42). Detailed information regarding these signatures was provided in Supplementary Table 3. Gene symbols in the signatures were transformed into Entrez IDs. Using coefficients provided in supplementary materials, continuous risk score of each signature was computed in TCGA-LUAD, GSE31210, and GSE68465 cohorts. For malignancy gene signature, risk score was generated in each set by first principal component of provided gene list. Hazard ratios of univariate and multivariate Cox regression were used to evaluate their survival associations. C-index derived from “coxph” function with default Efron method to handle ties was utilized to determine their prognostic classification capabilities.

Immunoscore-Clinical Prognostic Signature Construction

To make full use of both immunoscore and clinical variables in prognostic prediction, we constructed immunoscore-clinical prognostic signature (ICPS). Stage was converted to numeric variable. Stage IA, IB, IIA, IIB were assigned as 1, 2, 3, 4, respectively. Stage II with no subcategories were assigned as 3.5. Similarly, median value of ICPS in the training set was used as the cut-off value. C-index of ICPS was compared with stage or immunoscore alone using “compareC” package (43).

Genomic Analysis

Somatic mutation profile were downloaded from Genomic Data Common (GDC) website. Maftools was used to summarize the somatic mutation (44). Samples measured by Whole Genome Amplification (WGA) of Repli-G DNA (which could be identified by tumor barcode) were excluded to reduce possible bias. Tumor mutation burden (TMB) was calculated as previous study described:

$\begin{array}{l}{\text{TMB}}_{\text{i}}=1.0*\text{ NT}{\text{M}}_{\text{i}}+1.5*\text{ T}{\text{M}}_{\text{i}}\end{array}$

Tumor mutation burden of sample i was denoted by TMB_i. Total number of nontruncating mutation and total number of truncating mutation were denoted by NTMi and TMi, respectively.

The silent mutation was not included in the formula as it does not result in any change downstream. The truncating mutation was assigned a higher weight as it is considered more detrimental (45). Mutated genes between high-immunoscore and low-immunoscore groups were compared by fisher exact test using “mafcompare” function (44). Gene ontology and pathway analyses were then performed using differentially mutated genes by “clusterProfiler” package (46).

Tumor Purity and Various Cell Composition Characterization

We established our immunoscore model based on the bulk gene expression data of the tumor. It could also be used as the measurement of hazardous level of tumor environment (TME) with its biological background. TME contained not only cancer cells, but surrounding non-cancerous immune and normal cells. To further delineate the correlation between immunoscore and TME, we need to first figure out the TME components. TCGA-LUAD cohort was used for TME evaluation. Tumor purity, the percentage of cancer cell inside the tumor, could be estimated in different ways. Aran et al. developed consensus measurement of purity estimations (CPE), which used the median value of several genomic algorithms and immunohistochemistry (IHC) after normalization by combined mean and standard deviation (47). As a result, the bias introduced by a single method or algorithm was minimized. We also utilized Estimating the Proportion of Immune and Cancer cells (EPIC) algorithm, a method to predict various cell types in tumor tissue using RNA-seq tumor gene expression profile (48). Non-log transformed TPM data of TCGA-LUAD samples were used and the Ensemble gene IDs were converted into the official gene symbols as the algorithm required. The EPIC algorithm was based on reference gene expression profiles to infer proportions of surrounding non-malignant cells with experimental measurements confirming its predictive robustness. Samples with convergence code other than 0 were excluded as these might be outliers.

Gene Set Enrichment Analysis and Association With Clinical or Molecular Variables

Gene set enrichment analysis was performed to assess the association of immunoscore to the functional immune pathways. Differential gene expression profile between high-immunoscore and low-immunoscore subgroups was derived by “eBayes” function using limma package (49). We run fgsea algorithm with top 12,000 genes using C5 gene set from MsigDB database (https://www.gsea-msigdb.org/gsea/msigdb/). Gene set related to immune system were extracted. The correlation of immunoscore with clinical factors and certain molecular markers were also evaluated.

Statistical Analysis

Group comparison between a continuous variable were conducted by t-test or ANOVA. All correlation analyses were performed with spearman method, and considered highly correlated when absolute value of correlation coefficient was >0.25. False discovery rate was calculated as the adjusted P-value. All statistical procedures were conducted by R software version 3.6.1 (50). All p-values were two-sided and considered statistically significant when <0.05. Gene set was with P < 0.05 and FDR < 0.25 was considered significantly enriched.

Result

Immunoscore Model Construction

The flowchart of our study procedures was illustrated in Figure 1. A total of 975 patients with early-stage lung adenocarcinoma were included in the study. Detailed clinical information was shown in Table 1. In the training set, 109 genes were correlated with worse prognosis (HR > 1, P < 0.05, Supplementary Table 1). Gene set “TGFb family members,” “TGFb_Family_Member” and “Interferons” contained only 1 gene and were excluded from further analysis. Z-score transformed enrichment scores of the remaining 11 gene set were then calculated as the method described. All of them were correlated with poor survival in the training set (Figure 2A; Supplementary Table 2). Ridge Cox regression was then performed and immunoscore was derived by the sum of all Z-score transformed enrichment scores weighed by ridge regression coefficients (Figures 2B,C; Supplementary Table 2). The predictive accuracy of immunoscore to 2, 3, and 5-year survival were estimated by time-dependent receiver ROC analysis (Figure 2D).

FIGURE 1

Figure 1. Flowchart of the study. GSEA, gene set enrichment analysis.

TABLE 1

Table 1. Detailed patient clinical characteristics.

FIGURE 2

Figure 2. Immunoscore model construction. (A) Forest plot exhibiting different immune gene sets and patient overall survival in the training set. (B) 10-fold cross-validation for tuning parameter selection in the ridge regression model. The partial likelihood deviance is plotted against log (λ), where λ is the tuning parameter. Partial likelihood deviance values are shown, with error bars representing standard error (SE). The dotted vertical line on the left was drawn by minimum criteria whereas the line on the right represented the 1-SE criteria. We chose the minimum criteria. (C) Ridge regression coefficients of the 13 immune gene sets. The dotted line indicated the value chosen by the minimum criteria of the 10-fold validation. (D) Time-dependent receiver operator analysis (ROC) of the immunoscore in the training set. HR.95%CI, hazard ratio with 95% confidence interval. P.adj, adjusted P-value by false discovery rate.

Validation of Immuoscore

The immunoscore of the testing sets were calculated using the same formula. We also built a meta-set by combing all training and testing sets. Patients were stratified into high and low-immunoscore subgroups using median value of immunoscore in the training set as the cut-off value (−0.0126). Kaplan–Meier survival analysis and log-rank test was performed to compare the difference between these two subgroups. The result exhibited that patients from high-immunoscore subgroup were more likely to suffer worse overall survival (P < 0.001 in the training set, testing sets, and meta-set; Figures 3A–D). Similarly, patients with higher score were also linked to shorter disease-free survival (DFS) interval (P < 0.001 in training set, testing set 1 and meta-set, P = 0.005 in testing set 2, Supplementary Figure 1). Time-dependent ROC analyses were also performed to testing sets and meta-set (Supplementary Figure 2).

FIGURE 3

Figure 3. Survival analysis of the immunoscore. Kaplan–Meier curves for patient overall survival by immunoscore group in the (A) training set, (B) testing set 1, (C) testing set 2, and (D) meta-set.

Cox regression was used to assess its survival association. Univariate Cox regression analysis revealed that immunoscore was a significant risk factor in all three training and testing sets (HR = 3.11, 95% confidence interval (CI) 2.4–4.04, P < 0.001 in training set; HR = 2.39, 95% CI 1.6–3.58, P < 0.001 in testing set 1; HR = 1.44, 95% CI 1.17–1.78, P < 0.001 in testing set 2; HR = 1.88, 95% CI 1.63–2.17 in meta-set; Figure 4). Multivariate Cox regression analysis indicated that immunoscore was an independent risk factor in training (HR = 2.96, 95% CI 2.24–3.9, P < 0.001), testing set 1 (HR = 1.99, 95% CI 1.21–3.26, P = 0.006), testing set 2 (HR= 1.48, 95% CI 1.13–1.93, P = 0.005), and meta-set (HR = 2.01, 95% CI 1.69–2.39, P < 0.001), as shown in Figure 5. Moreover, Immunoscore could identify patients with worse survival in all clinical subgroups in meta-set (Supplementary Figure 3).

FIGURE 4

Figure 4. The univariate Cox analysis of the immunoscore and clinicopathological factors. The HR in training cohort was 3.11, with 95% confidence interval (CI) from 2.44 to 4.04 (P < 0.001). The HR in testing set 1 was 2.39, with 95% CI from 1.6 to 3.58 (P < 0.001). The HR in testing set 2 was 1.44, with 95% CI from 1.17 to 1.78 (P < 0.001). The HR in meta-set was 1.88, with 95% CI from 1.63 to 2.17 (P < 0.001). HR.95%CI, hazard ratio with 95% confidence interval.

FIGURE 5

Figure 5. Multivariate Cox analysis evaluating independently predictive ability of immunoscore for patient survival. The immunoscore was able to independently predict patient survival in training set (hazard ratio (HR) = 2.96, 95% confidence interval (CI) from 2.24 to 3.9, P < 0.001), testing set 1 (HR = 1.99, 95% CI from 1.05 to 3.81, P = 0.006), testing set 2 (HR = 1.48, 95% CI from 1.13 to 1.93, P=0.005), and meta-set (HR = 2.01, 95% CI from 1.69 to 2.39, P < 0.001). HR.95%CI, hazard ratio with 95% confidence interval.

Comparison of Immunoscore With Other Genomic Signatures

To assess the utility of immunoscore model, we compared prognostic association of immunoscore against other published genomic signatures (Supplementary Table 3). Besides Song et al. signature, most signatures had good performance in univariate and multivariate regression analyses (Supplementary Figure 4; Figure 6A). Immunoscore exhibited a generally higher C-index than other signatures in all three cohorts, except less than Chen2 et al. signature in GSE31210 (0.72 vs. 0.726, Figure 6B). Meanwhile, immunoscore achieved the highest mean C-index (0.68 vs. range from 0.58 to 0.64, Figure 6B).

FIGURE 6

Figure 6. Comparison of immunoscore and other existing NSCLC signatures. (A) Hazard ratio of each gene expression signature in multivariable Cox analysis. (B) C-index of each signature in each independent dataset and mean C-index. HR.95%CI, hazard ratio with 95% confidence interval.

Immunoscore-Clinical Prognostic Signature Construction

Stage, age and immunoscore were all independent prognostic variables in multivariable Cox analysis in all 3 sets and meta-set. To explore whether combing these variables would improve prediction accuracy, coefficients of multivariate regression of these three factors in the training-set were used to introduce a new variable, immunoscore-clinical prognostic signature (ICPS).

$\begin{array}{l}ICPS=1.06076428*immunoscore\\ +0.19653598*stage+0.01961085*age\end{array}$

Patients were stratified into high-ICPS and low-ICPS subgroups using median value of ICPS in training set as the cut-off (1.74). High-ICPS subgroup was significantly correlated with worse survival in each set (P < 0.001, Figure 7). Figure 7 also exhibited C-index of ICPS was significantly higher than either immunoscore or stage, in training [0.72 (ICPS) vs. 0.7 (immunoscore) and 0.59 (stage), P < 0.001 when compared with stage], testing set 1 [0.75 (ICPS) vs. 0.72 (immunoscore) and 0.7 (stage), P = 0.015 when compared with stage], and testing set 2 [0.65 (ICPS) vs. 0.61 (immunoscore) and 0.62 (stage), P < 0.001 when compared with stage]. Moreover, C-index of ICPS was significantly higher than both of them in meta-set [0.7 (ICPS) vs. 0.66 (immunoscore) and 0.64 (stage), P < 0.001 when compared with immunoscore or stage, Figure 7].

FIGURE 7

Figure 7. Kaplan–Meier survival analysis and compare C-index of ICPS with immunoscore and stage in (A) training set, (B) testing set 1, (C) testing set 2, (D) meta-set.

Immunoscore, ICPS, and Adjuvant Therapy

A small subset of early-stage LUAD patient received postoperative adjuvant chemotherapy or radiotherapy. To investigate whether various treatment strategies had an effect to immunoscore and ICPS model, we used GSE68465 cohort with comprehensive documentation of adjuvant therapy. Patients who received adjuvant therapy had a worse overall survival (Figure 8A). It might be due to clinical practice, as adjuvant therapy was more likely to be applied to patients with higher stage and worse condition. Survival analysis indicated that immunoscore and ICPS could still stratify patients with different prognosis in each treatment group (Figures 8B,C).

FIGURE 8

Figure 8. (A) Kaplan–Meier curves for patients who received adjuvant therapy or not (B) Kaplan–Meier curves for survival prediction by the immunoscore in patients who received adjunctive therapy or not. (C) Kaplan–Meier curves for survival prediction by the ICPS in patients who received adjunctive therapy or not.

Genome Analysis

To explore the possible underlying causes of difference in immunoscore between patients, we searched GDC website and downloaded all available somatic mutation data of lung adenocarcinoma patients. Three hundred fifty-eight available mutation profiles in TCGA-LUAD cohort (174 in high-immunoscore subgroup, 148 in low-immunoscore subgroup) were summarized by maftools. The mutation profiles of high-immunoscore and low-immunoscore subgroups were illustrated in Figures 9A,B, respectively. Differentially mutated genes between low-immunocore and high-immunoscore subgroups were identified by Fisher exact test using “mafcompare” function. Twenty of them were shown in Figure 9C. TP53 was the most commonly mutated gene in high-immunoscore subgroup and had the smallest adjusted P-value. TP53 was a tumor suppressor gene, encoding P53 transcriptional factor which responds to DNA damage repair. TP53 mutation has been recently reported to be associated with response to immunotherapy in certain subtype of NSCLC (51). We discovered that TP53 mutation was correlated with immunoscore. P53 mutation might induce genome instability, increasing neoantigen load, leading to a more dangerous tumor immune microenvironment, resulting in higher immunoscore. Another established immunotherapy biomarker, tumor mutation burden (TMB), was also positively correlated with immunoscore (R = 0.22, P < 0.001, Figure 9D). Gene ontology and KEGG pathway analyses of the differentially mutated genes were provided in Supplementary Figure 5.

FIGURE 9

Figure 9. Mutation profile in the TCGA-LUAD cohort. (A) Mutation profile of high-immunoscore subgroup. (B) Mutation profile of low-immunoscore subgroup. (C) Differentially mutated genes between high and low immunoscore patients. (D) Correlation between tumor mutation burden (TMB) and immunoscore. OR.95%CI, odds ratio with 95% confidence interval. P.adj, adjusted P-value by false discovery rate.

Immunoscore and Tumor Microenvironment

The relationship between immunoscore and tumor microenvironment was investigated using TCGA-LUAD cohort. Tumor purity, the percentage of cancer cells inside the tumor, was estimated by consensus measurement of purity estimations (CPE). Patients with high immunoscore tend to have low tumor purity (R = −0.12, P = 0.015, Figure 10A). Patients were divided into high-purity and low-purity subgroup using median value of tumor purity (0.637). Kaplan-Meier survival curves indicated high tumor purity tend to have generally worse survival, but did not reach statistical significance (P = 0.3, Figure 10B). We next investigated the cellular composition of TME. EPIC algorithm, which was designed specifically for RNA-seq data, was used to infer the proportions of different infiltrating immune and stromal cells. Using absolute value of 0.25 as cut-off, cancer-associated fibroblast (CAF) (R = 0.32, P < 0.001) and CD8 T cell (R = −0.27, P < 0.001) were highly correlated to immunoscore (Figure 10C; Supplementary Table 4). In univariate Cox analysis, only CAF attained a borderline significant P value (HR = 2.42, 95% CI 0.9–6.55, P = 0.081, Figure 10D).

FIGURE 10

Figure 10. Tumor microenvironment (TME) change associated with immunoscore. (A) Correlation between immunoscore and tumor purity. (B) Kaplan-Meier curves of patient survival according to tumor purity in the TCGA-LUAD cohort. (C) Correlation matrix of immunoscore and cell proportions. (D) Univariate Cox analysis of various cell type. HR.95%CI, hazard ratio with 95% confidence interval. P.adj, adjusted P-value by false discovery rate.

Immunoscore, Clinicopathological Characteristics, and Biological Phenotypes

Gene set enrichment analysis between high-immunoscore and low-immunoscore were conducted. Immune-related pathways were extracted and most of them were enriched to the low end (27 out of 29 Immune-related pathways). The relationship between immunoscore and other clinicopathological factors were assessed in TCGA-LUAD cohort. Higher T and N stage possessed greater immunoscore, whereas its distribution in age, gender and smoking status was not significantly different (Figure 11B). Of immune checkpoint molecules, immunoscore was only correlated to PD-L1 and LAG3 (R = 0.16, P = 0.001 for PD-L1; R = 0.1, P = 0.04 for LAG3, Figure 11B; Supplementary Table 5). Interestingly, Several hypoxia-inducible factor (HIF)-1 pathway markers, like HIF-1A (R = 0.41, P < 0.001), SLC2A1 (R = 0.6, P < 0.001), LOXL2 (R = 0.55, P < 0.001), PDK1 (R = 0.27, P < 0.001), and LDHA (R = 0.53, P < 0.001), were highly correlated with immunoscore (Figure 11C, Supplementary Table 5) (52).

FIGURE 11

Figure 11. (A) Comparison of enrichment levels of immune-related pathways between high-immunoscore and low-immunoscore subgroups. (B) Distribution of immunoscore in various clinical subgroups. (C) Correlation matrix of immunoscore and certain gene expression.

Discussion

Lung cancer treatment has been improved dramatically during the past decades, mainly owing to the constant discoveries of genomic alterations during lung cancer pathogenesis. However, the patient prognostic evaluation is still based on the AJCC staging system. Although it is a powerful prognostic prediction tool, it is inadequate to get a precise assessment of patient survival. In early-stage LUAD, the AJCC staging system is far from getting accurate prediction since about 30 percent of patients would develop recurrence, with 2-year survival at about 17% (53). To identify this subset of patients with high risk of recurrence and poor survival is critical since receiving adjuvant chemotherapy or newly developed adjuvant immunotherapy may of great benefit to them.

Up to now, Numerous gene expression signatures have been established for the prediction of lung cancer patient survival (41). Few of them, however, have been translated into real clinical practice. It might be caused by several defects in signature construction. First, some of them were trained from a small cohort with high variance and insufficient independent samples to test its robustness. Second, Gene expression data were measured by different experimental strategy with batch effect, which means that the signature constructed in one specific cohort cannot be generalized into other platforms. Third, most of the signatures were composed of several specific genes and ignore other possible causes, which severely decreased its stability and could potentially lead to overfitting.

In our research, we compared immunoscore with other gene expression signatures. Immunoscore achieved the highest mean C-index, indicating its superior prognostic classification capability. Of immune-related gene signatures, Li et al. (40) signature also had good performance. Li's signature used the binary variable, the pairwise comparison between immune-related gene groups, as features in model construction. Immunoscore and Li's signature had a lot in common, as both of them used some sort of gene ranks (ssGSEA in immunoscore; pairwise comparison in Li's signature) rather than gene expression intensity, making them not sensitive to preprocessing strategies and batch effect.

Our model also has its biomedical sense. It was constructed based on enrichment score of risk genes from multiple immune gene sets, and all selected immune gene sets were significantly correlated with worse patient survival. Thus, higher immunoscore indicated a more dangerous tumor microenvironment. The top three contributors to immunoscore were cytokine receptor, antimircrobial, and cytokine. Several cytokine-cytokine receptors signaling pathway have been identified to play a important role in cancer cell proliferation and survival. Most cytokine receptors were located at cell surface, and activated when contacting with specific cytokines. In GSEA analysis, innate immune response activating cell surface receptor signaling pathway ranked the top. Gene ontology analysis also indicated several gene sets related to cell membrane are enriched in differentially mutated genes. Overall, it implied that cell surface signaling pathways were tightly linked to immunoscore and disruption of these pathways might portend poor prognosis. In addition, drugs modifying cytokine-cytokine receptor signaling in combination with other immunotherapy might be a promising treatment strategy. Antimicrobial pathway has been linked to carcinogenesis, as infection by some microorganisms might lead to cell proliferation, and could be reversed by antimicrobials agents (54).

Tumor purity and cellular composition in tumor microenvironment were also investigated. Patients with high immunoscore tend to have low tumor purity. Furthermore, immunoscore was positively associated with CAFs and but inversely associated with CD8+ T cells. CD8+ T cell has direct cytolytic effect, whereas CAF, on the other hand, may suppress CD8+ cell function by upregulating PD-1, PD-L1, and FAS ligand on Treg cells (55). In addition, KEGG pathway analysis of differentially mutated genes also found ECM-interaction pathway abnormality. ECM stiffness might lead to activation of cancer cells and pro-tumor effect of CAF (56). Besides, cancer cell could induce CAF to remodel ECM, whereas CAF might sustain cancer growth by secreting aspartate (57). Further investigations are needed to figure out how fibroblast communicate with other cells or molecules inside TME and give insight to novel drug targets.

We next explored the phenotypical difference between samples of high and low immunoscore.

Most immune-related pathways were enriched in low-immunoscore subgroup, indicating high-immunoscore subgroup was a “immune cold” subtype. We also discovered multiple markers of HIF-hydroxylase oxygen-sensing pathway to be correlated with immunoscore. HIF could enhance tumor proliferation in TME by altering immune cell function and recruiting pro-tumor immune cells (58). For example, expression of HIF1A in tumor-associated macrophage (TAM) might suppress T cell function (59). More experiments and analyses are required to elucidate how HIF pathway affect tumor immune microenvironment as HIF1A is an incredibly promising target for cancer therapy (60).

Our study has several advantages. First, we trained our model in a large cohort with sufficient samples used to test its robustness. Second, we built our immunoscore model to predict patient outcome based on the enrichment levels of different gene sets rather than several single genes, making it a more comprehensive evaluation of tumor immune microenvironment and prevent overfitting. Third, when integrating clinical factors and immunoscore to construct a new ICPS model, it outperformed either immunoscore or stage alone. Fourth, immunoscore itself could also be seen as a proxy variable, the measurement of tumor immune microenvironment, and we found that genome instability, several specific immune cell proportions and functional pathway activation were correlated to immunoscore.

We admit some limitations. First, we used publically available datasets in retrospective manner. We did not have all clinical information needed for the study. For example, patients with inherent immune disorder or taking drugs with impact on immune system should be ruled out. Second, gene expression signatures were developed in different platform with diverse preprocessing strategies and normalization procedure. Although immunoscore outperformed other signatures, it might be due to technical bias or batch effect. Third, the immunoscore model contained several genes with still unknown effects in LUAD, and this “black-box” impact severely undermined the model interpretability. More experiments are needed to elucidate their biological associations. Finally, we cannot estimate its predictive value to immune checkpoint inhibitors due to lack of response data to immunotherapy. Further studies are needed to validate and improve immunoscore model.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: https://portal.gdc.cancer.gov, https://www.ncbi.nlm.nih.gov/geo.

Author Contributions

ZZ, YW, and BW were responsible for the study design. The analysis was performed by ZZ and DZ. ZZ and BW were involved in interpretation of the data. The manuscript was drafted by ZZ. JX and DZ have revised the manuscript. All authors read and approved the final manuscript.

Funding

This work was supported by scientific research project of maternal and child health of Jiangsu province (Grant Number F201865).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2020.00691/full#supplementary-material

Supplementary Figure 1. Survival analysis of the immunoscore. Kaplan-Meier curves for patient disease-free survival by immunoscore group in the (A) training set, (B) testing set 1, (C) testing set 2, and (D) meta-set.

Supplementary Figure 2. Time-dependent receiver operator analysis (ROC) of the immunoscore in the (A) testing set 1, (B) testing set 2, and (C) meta-set.

Supplementary Figure 3. Subgroup analysis of immunoscore. Immunoscore was a significant risk factor in each clinical subgroup. HR.95%CI, hazard ratio with 95% confidence interval.

Supplementary Figure 4. Hazard ratio of each gene expression signature in univariable Cox analysis. HR.95%CI, hazard ratio with 95% confidence interval.

Supplementary Figure 5. Functional analysis of differentially mutated genes in genome analysis. Top 20 (A) biological process, (B) molecular function, (C) cellular component, and (D) KEGG pathway. KEGG, Kyoto Encyclopedia of Genes and Genomes.

References

1. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. (2018) 68:394–424. doi: 10.3322/caac.21492

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Gridelli C, Rossi A, Carbone DP, Guarize J, Karachaliou N, Mok T, et al. Non-small-cell lung cancer. Nat Rev Dis Prim. (2015) 1:15009. doi: 10.1038/nrdp.2015.9

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Never-smoker NE-s. Comprehensive molecular profiling of lung adenocarcinoma. Nature. (2014) 511:543–50. doi: 10.1038/nature13385

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Baltayiannis N, Chandrinos M, Anagnostopoulos D, Zarogoulidis P, Tsakiridis K, Mpakas A, et al. Lung cancer surgery: an up to date. J Thorac Dis. (2013) 5 (Suppl. 4):S425. doi: 10.3978/j.issn.2072-1439.2013.09.17

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Goldstraw P, Chansky K, Crowley J, Rami-Porta R, Asamura H, Eberhardt WE, et al. The IASLC lung cancer staging project: proposals for revision of the TNM stage groupings in the forthcoming (eighth) edition of the TNM classification for lung cancer. J Thorac Oncol. (2016) 11:39–51. doi: 10.1016/j.jtho.2015.09.009

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Pisters KM, Evans WK, Azzoli CG, Kris MG, Smith CA, Desch CE, et al. Cancer Care Ontario and American Society of Clinical Oncology adjuvant chemotherapy and adjuvant radiation therapy for stages I-IIIA resectable non–small-cell lung cancer guideline. J Clin Oncol. (2007) 25:5506–18. doi: 10.1200/JCO.2007.14.1226

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Group IALCTC. Cisplatin-based adjuvant chemotherapy in patients with completely resected non–small-cell lung cancer. New Engl J Med. (2004) 350:351–60. doi: 10.1056/NEJMoa031644

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Vansteenkiste J, Wauters E, Reymen B, Ackermann CJ, Peters S, De Ruysscher D. Current status of immune checkpoint inhibition in early-stage NSCLC. Ann Oncol. (2019) 30:1244–53. doi: 10.1093/annonc/mdz175

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Gajewski TF, Schreiber H, Fu Y-X. Innate and adaptive immune cells in the tumor microenvironment. Nat Immunol. (2013) 14:1014. doi: 10.1038/ni.2703

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Fridman WH, Zitvogel L, Sautès–Fridman C, Kroemer G. The immune contexture in cancer prognosis and treatment. Nat Rev Clin Oncol. (2017) 14:717. doi: 10.1038/nrclinonc.2017.101

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Brambilla E, Le Teuff G, Marguet S, Lantuejoul S, Dunant A, Graziano S, et al. Prognostic effect of tumor lymphocytic infiltration in resectable non–small-cell lung cancer. J Clin Oncol. (2016) 34:1223. doi: 10.1200/JCO.2015.63.0970

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Devarakonda S, Rotolo F, Tsao M-S, Lanc I, Brambilla E, Masood A, et al. Tumor mutation burden as a biomarker in resected non-small-cell lung cancer. J Clin Oncol. (2018) 36:2995–3006. doi: 10.1200/JCO.2018.78.1963

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Tsao M-S, Le Teuff G, Shepherd F, Landais C, Hainaut P, Filipits M, et al. PD-L1 protein expression assessed by immunohistochemistry is neither prognostic nor predictive of benefit from adjuvant chemotherapy in resected non-small cell lung cancer. Ann Oncol. (2017) 28:882–9. doi: 10.1093/annonc/mdx003

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Altorki NK, Markowitz GJ, Gao D, Port JL, Saxena A, Stiles B, et al. The lung microenvironment: an important regulator of tumour growth and metastasis. Nat Rev Cancer. (2019) 19:9–31. doi: 10.1038/s41568-018-0081-9

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Barbie DA, Tamayo P, Boehm JS, Kim SY, Moody SE, Dunn IF, et al. Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature. (2009) 462:108–12. doi: 10.1038/nature08460

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Senbabaoglu Y, Gejman RS, Winer AG, Liu M, Van Allen EM, de Velasco G, et al. Tumor immune microenvironment characterization in clear cell renal cell carcinoma identifies prognostic and immunotherapeutically relevant messenger RNA signatures. Genome Biol. (2016) 17:231. doi: 10.1186/s13059-016-1092-z

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Shen S, Wang G, Zhang R, Zhao Y, Yu H, Wei Y, et al. Development and validation of an immune gene-set based Prognostic signature in ovarian cancer. EBioMedicine. (2019) 40:318–26. doi: 10.1016/j.ebiom.2018.12.054

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Verhaak RG, Tamayo P, Yang J-Y, Hubbard D, Zhang H, Creighton CJ, et al. Prognostically relevant gene signatures of high-grade serous ovarian carcinoma. J Clin Investig. (2012) 123:517–25. doi: 10.1172/JCI65833

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Grossman RL, Heath AP, Ferretti V, Varmus HE, Lowy DR, Kibbe WA, et al. Toward a shared vision for cancer genomic data. New Engl J Med. (2016) 375:1109–12. doi: 10.1056/NEJMp1607591

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Yamauchi M, Yamaguchi R, Nakata A, Kohno T, Nagasaki M, Shimamura T, et al. Epidermal growth factor receptor tyrosine kinase defines critical prognostic genes of stage I lung adenocarcinoma. PLoS ONE. (2012) 7:e43923. doi: 10.1371/journal.pone.0043923

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Okayama H, Kohno T, Ishii Y, Shimada Y, Shiraishi K, Iwakawa R, et al. Identification of genes upregulated in ALK-positive and EGFR/KRAS/ALK-negative lung adenocarcinomas. Cancer Res. (2012) 72:100–11. doi: 10.1158/0008-5472.CAN-11-1403

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Shedden K, Taylor JM, Enkemann SA, Tsao M-S, Yeatman TJ, Gerald WL, et al. Gene expression–based survival prediction in lung adenocarcinoma: a multi-site, blinded validation study. Nat Med. (2008) 14:822. doi: 10.1038/nm.1790

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Barrett T, Troup DB, Wilhite SE, Ledoux P, Rudnev D, Evangelista C, et al. NCBI GEO: archive for high-throughput functional genomic data. Nucleic Acids Res. (2009) 37(Suppl_1):D885–90. doi: 10.1093/nar/gkn764

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Gao J, Aksoy BA, Dogrusoz U, Dresdner G, Gross B, Sumer SO, et al. Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal. Sci. Signal. (2013) 6:pl1. doi: 10.1126/scisignal.2004088

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Cerami E, Gao J, Dogrusoz U, Gross BE, Sumer SO, Aksoy BA, et al. The cBio cancer genomics portal: an open platform for exploring multidimensional cancer genomics data. Cancer Discov. (2012) 2:401–4. doi: 10.1158/2159-8290.CD-12-0095

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Jacobsen A, Luna A. cgdsr: R-Based API for Accessing the MSKCC Cancer Genomics Data Server (CGDS). R Package Version 1.3.0. (2015). Available online at: https://CRAN.R-project.org/package=cgdsr

27. Davis S, Meltzer PS. GEOquery: a bridge between the Gene Expression Omnibus (GEO) and bioconductor. Bioinformatics. (2007) 23:1846–7. doi: 10.1093/bioinformatics/btm254

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Gautier L, Cope L, Bolstad BM, Irizarry RA. affy—analysis of Affymetrix GeneChip data at the probe level. Bioinformatics. (2004) 20:307–15. doi: 10.1093/bioinformatics/btg405

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Dai M, Wang P, Boyd AD, Kostov G, Athey B, Jones EG, et al. Evolving gene/transcript definitions significantly alter the interpretation of GeneChip data. Nucleic Acids Res. (2005) 33:e175. doi: 10.1093/nar/gni179

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Irizarry RA, Hobbs B, Collin F, Beazer-Barclay YD, Antonellis KJ, Scherf U, et al. Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics. (2003) 4:249–64. doi: 10.1093/biostatistics/4.2.249

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Kauffmann A, Gentleman R, Huber W. arrayQualityMetrics—a bioconductor package for quality assessment of microarray data. Bioinformatics. (2009) 25:415–6. doi: 10.1093/bioinformatics/btn647

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, et al. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res. (2015) 44:e71. doi: 10.1093/nar/gkv1507

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Wagner GP, Kin K, Lynch VJ. Measurement of mRNA abundance using RNA-seq data: RPKM measure is inconsistent among samples. Theory Biosci. (2012) 131:281–5. doi: 10.1007/s12064-012-0162-3

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Therneau T. A Package for Survival Analysis in S. Version 2.38 (2015). Available online at: https://CRAN.R-project.org/package=survival

35. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. (2013) 14:7. doi: 10.1186/1471-2105-14-7

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Hoerl AE, Kennard RW. Ridge regression: biased estimation for nonorthogonal problems. Technometrics. (1970) 12:55–67. doi: 10.1080/00401706.1970.10488634

CrossRef Full Text | Google Scholar

37. Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Softw. (2010) 33:1. doi: 10.18637/jss.v033.i01

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Blanche P, Dartigues JF, Jacqmin-Gadda H. Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks. Stat Med. (2013) 32:5381–97. doi: 10.1002/sim.5958

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Song Q, Shang J, Yang Z, Zhang L, Zhang C, Chen J, et al. Identification of an immune signature predicting prognosis risk of patients in lung adenocarcinoma. J Transl Med. (2019) 17:70. doi: 10.1186/s12967-019-1824-4

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Li B, Cui Y, Diehn M, Li R. Development and validation of an individualized immune prognostic signature in early-stage nonsquamous non–small cell lung cancer. JAMA Oncol. (2017) 3:1529–37. doi: 10.1001/jamaoncol.2017.1609

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Tang H, Wang S, Xiao G, Schiller J, Papadimitrakopoulou V, Minna J, et al. Comprehensive evaluation of published gene expression prognostic signatures for biomarker-based lung cancer clinical studies. Ann Oncol. (2017) 28:733–40. doi: 10.1093/annonc/mdw683

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Chen D-T, Hsu Y-L, Fulp WJ, Coppola D, Haura EB, Yeatman TJ, et al. Prognostic and predictive value of a malignancy-risk gene signature in early-stage non–small cell lung cancer. J Natl Cancer Inst. (2011) 103:1859–70. doi: 10.1093/jnci/djr420

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Kang L, Chen W, Petrick NA, Gallas BD. Comparing two correlated C indices with right-censored survival outcome: a one-shot nonparametric approach. Stat Med. (2015) 34:685–703. doi: 10.1002/sim.6370

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Mayakonda A, Lin D-C, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. (2018) 28:1747–56. doi: 10.1101/gr.239244.118

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Liu Z, Li M, Jiang Z, Wang X. A comprehensive immunologic portrait of triple-negative breast cancer. Transl Oncol. (2018) 11:311–29. doi: 10.1016/j.tranon.2018.01.011

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Yu G, Wang L-G, Han Y, He Q-Y. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics. (2012) 16:284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Aran D, Sirota M, Butte AJ. Systematic pan-cancer analysis of tumour purity. Nat Commun. (2015) 6:8971. doi: 10.1038/ncomms9971

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Racle J, de Jonge K, Baumgaertner P, Speiser DE, Gfeller D. Simultaneous enumeration of cancer and immune cell types from bulk tumor gene expression data. Elife. (2017) 6:e26476. doi: 10.7554/eLife.26476.049

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. (2015) 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

50. R Core Team. R: A language and environment for statistical computing. Vienna: R Foundation for Statistical Computing (2019). Available online at: https://www.R-project.org/

PubMed Abstract | Google Scholar

51. Chen Y, Chen G, Li J, Huang Y-Y, Li Y, Lin J, et al. Association of tumor protein p53 and ataxia-telangiectasia mutated comutation with response to immune checkpoint inhibitors and mortality in patients with non–small cell lung cancer. JAMA Netw Open. (2019) 2:e1911895. doi: 10.1001/jamanetworkopen.2019.11895

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Luo W, Wang Y. Epigenetic regulators: multifunctional proteins modulating hypoxia-inducible factor-α protein stability and activity. Cell Mol Life Sci. (2018) 75:1043–56. doi: 10.1007/s00018-017-2684-9

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Fedor D, Johnson WR, Singhal S. Local recurrence following lung cancer surgery: incidence, risk factors, and outcomes. Surg Oncol. (2013) 22:156–61. doi: 10.1016/j.suronc.2013.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Alibek K, Bekmurzayeva A, Mussabekova A, Sultankulov B. Using antimicrobial adjuvant therapy in cancer treatment: a review. Infect Agents Cancer. (2012) 7:33. doi: 10.1186/1750-9378-7-33

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Nazareth MR, Broderick L, Simpson-Abelson MR, Kelleher RJ, Yokota SJ, Bankert RB. Characterization of human lung tumor-associated fibroblasts and their ability to modulate the activation of tumor-associated T cells. J Immunol. (2007) 178:5552–62. doi: 10.4049/jimmunol.178.9.5552

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Pickup MW, Mouw JK, Weaver VM. The extracellular matrix modulates the hallmarks of cancer. EMBO Rep. (2014) 15:1243–53. doi: 10.15252/embr.201439246

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Bertero T, Oldham WM, Grasset EM, Bourget I, Boulter E, Pisano S, et al. Tumor-stroma mechanics coordinate amino acid availability to sustain tumor growth and malignancy. Cell Metab. (2019) 29:124–40.e110. doi: 10.1016/j.cmet.2018.09.012

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Triner D, Shah YM. Hypoxia-inducible factors: a central link between inflammation and cancer. J Clin Investig. (2016) 126:3689–98. doi: 10.1172/JCI84430

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Doedens AL, Stockmann C, Rubinstein MP, Liao D, Zhang N, DeNardo DG, et al. Macrophage expression of hypoxia-inducible factor-1α suppresses T-cell function and promotes tumor progression. Cancer Res. (2010) 70:7465–75. doi: 10.1158/0008-5472.CAN-10-1439

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Soni S, Padwad YS. HIF-1 in cancer therapy: two decade long story of a transcription factor. Acta Oncol. (2017) 56:503–15. doi: 10.1080/0284186X.2017.1301680

PubMed Abstract | CrossRef Full Text | Google Scholar