The Prognostic Model and Drug Sensitivity of LKB1-Mutant Lung Adenocarcinoma Based on Immune Landscape

Background: Lung cancer is the most common cause of cancer-related deaths worldwide. LKB1-mutant lung adenocarcinoma (LUAD) is a unique subtype of this deadly cancer. LKB1 mutations cause functional changes in a variety of cell processes, including immune functions, that affect prognosis. To date, the potential role of immunity in the prognosis of LKB1-mutant LUAD is not well understood.

Methods: We systematically analyzed immune-related genes in LUAD samples from The Cancer Genome Atlas (TCGA) database. ESTIMATE and CIBERSORT algorithms were used to explore the immune microenvironment. A prognostic risk model was constructed, and prognostic, immune function, drug sensitivity, and model specificity analyses were performed to identify the effectiveness of the model.

Results: Our results showed that LKB1 mutations suppressed immune function in LUAD. A three-gene signature was constructed to stratify patients into two risk groups. The risk score was an independent predictor for overall survival (OS) in multivariate Cox regression analyses [hazard ratio (HR) > 1, p = 0.002]. Receiver operating characteristic (ROC) curve analyses confirmed that the risk score has better performance than clinicopathological characteristics. Functional analysis revealed that the immune status was different between the risk groups. ZM.447439 was an appropriate treatment for the high-risk group of patients. This risk model is only suitable for LKB1-mutant tumors; it performed poorly in LUAD patients with wild-type LKB1.

Conclusion: Our findings indicate the potential role of immunity in LKB1-mutant LUAD, providing novel insights into prognosis and guiding effective immunotherapy.

Introduction

Lung cancer is the most common cause of cancer-related deaths worldwide, with more than 40% of cases being lung adenocarcinoma (LUAD) (Jordan et al., 2017). Genomic alterations have important impacts on tumor cell-intrinsic and non-cell-autonomous cancer hallmarks; accordingly, genomic alterations underlie the molecular and clinical heterogeneity of LUAD (Skoulidis and Heymach, 2019). The gene encoding liver kinase B1 (LKB1) is the third most commonly mutated gene in LUAD (after KRAS and TP53), with approximately 19% of LUAD cases involving LKB1 mutations (Ding et al., 2008; Cancer Genome Atlas Research, 2014). Moreover, LKB1 is involved in many cell processes essential for cell survival, such as metabolic balance, maintenance of DNA integrity, proliferation, and polarity establishment (Shackelford and Shaw, 2009). Notably, therapeutic approaches targeting oncogene driver mutations, for example, activating mutations in the epidermal growth factor receptor or rearrangements in anaplastic lymphoma kinase, have been recently shown to elicit dramatic clinical responses (Paez et al., 2004; Soda et al., 2007; Bergethon et al., 2012). However, there are no routinely used clinical drugs that specifically target LUAD with LKB1-inactivating mutations.

The tumor microenvironment (TME) refers to the local biological environment including cancer cells, stromal cells, and distant recruited cells, such as infiltrating immune cells (myeloid cells and lymphocytes), bone marrow-derived cells, and secreted factors such as cytokines and chemokines (Wang et al., 2017; Zhu et al., 2020). The status of the TME plays a bidirectional role in tumor progression and often impacts the effectiveness of targeted drugs (Topalian et al., 2012; Tamborero et al., 2018). While the complexity of TME components has proven to be a barrier to research, some new algorithms, including ESTIMATE and CIBERSORT, that are based on bulk RNA sequencing have dealt with this problem effectively (Yoshihara et al., 2013; Newman et al., 2015a). These algorithms can quantify the degree of stromal cell and immune cell infiltration in tumor tissues, and their accuracy has been verified in breast cancer and liver cancer models (Newman et al., 2015a).

Recent work has shown that oncogenic mutations shape the TME and determine its immune context (Akbay et al., 2013). Because oncogenes directly activate immune checkpoints to impact immune-evading mechanisms and cause immune suppression (Akbay et al., 2013), immunotherapies for immune checkpoint blockade and other aspects efficiently target tumors and show an association with the survival of cancer patients (Taube et al., 2014; Tumeh et al., 2014; Drerup et al., 2020). It is unclear whether the inactivation of tumor suppressor genes such as LKB1 exerts similar effects, even though recent studies have opened new perspectives in this field.

Thus, there is an urgent need to elucidate the underlying mechanisms influencing immune activity and the TME in LKB1-mutant LUAD and to explore potential prognostic biomarkers. Understanding these factors could provide better therapeutic strategies for LUAD. Furthermore, in addition to conventional biomarkers, the multi-faceted roles of long non-coding RNA (lncRNA) and microRNA (miRNA) in the progression of LUAD, including in determining immune and microenvironmental conditions (Shen et al., 2020; Sun et al., 2020), mean that these molecules may also serve as useful indicators of disease progression. Therefore, we have included all kinds of RNA in our research.

To this end, we aim to explore the impact of LKB1 mutations on the immune microenvironment of LUAD and the mechanisms that underlie this impact. In addition, we constructed a prognostic model based on immune-related RNAs to promote the development of therapies in LKB1-mutant LUAD and to provide novel prognostic biomarkers.

Materials and Methods

Data Downloading and Preprocessing

LUAD gene expression data and clinical information were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov), and 14,086 lncRNAs, 19,605 mRNAs, and 2,191 miRNAs were obtained by data annotation. According to the mutation status of LKB1 obtained from the cBioPortal database (https://www.cbioportal.org), we separated TCGA tumor samples into LKB1-mutant (LKB1-mu) and LKB1-wild-type (LKB1-wt) groups. A total of 1,959 immune-related genes were retrieved from the ImmPort database (https://www.immport.org/resources), and two immune-related gene sets (immune system process M13664 and immune response M19817) were extracted from the Molecular Signatures Database v4.0 (http://www.broadinstitute.org/gsea/msigdb/index.jsp).

We identified 89 patients with LKB1 mutations by searching the 31 pan-cancer studies (excluding LUAD) in the cBioPortal database. The miRNA expression data of these patients were obtained from the TCGA database, and 74 non-LUAD LKB1-mutated patients were ultimately included as a validation cohort.

Effect of the LKB1-Mutation on the Immune Landscape of Lung Adenocarcinoma

The LKB1 gene mutation status of LUAD was analyzed by probing the cBioPortal database. The difference in the expression level of LKB1 between the LKB1-wt and LKB1-mu groups was verified by t-test. The Gene Set Enrichment Analysis (GSEA) was performed to explore the biological function regulated by the LKB1 mutant in LUAD. The stromal score, immune score, and ESTIMATE score were also analyzed by the ESTIMATE algorithm based on transcriptome expression profiles of LUAD to examine the effect of the LKB1 mutant (Yoshihara et al., 2013). The CIBERSORT algorithm was used to accurately determine the composition of 22 kinds of immune cells in a large tumor sample dataset, and the impact of the LKB1 mutant on these immune cells was determined (Newman et al., 2015b). The TIMER database (https://cistrome.shinyapps.io/timer/) was used to explore the correlation between LKB1 and immunity. The differential expression of 46 immune checkpoint genes between LKB1-mu and LKB1-wt cells was further explored.

Identifying Candidate Genes to Construct a Prognostic Risk Model

Expression levels of mRNAs, lncRNAs, and miRNAs in the LKB1-wt and LKB1-mu groups were analyzed by the edger package with the difference screening parameters |logFC| > 1 and false discovery rate (FDR) < 0.05. Utilizing a total of 1959 immune-related genes extracted from the databases, immune-related lncRNAs were identified according to the criteria of |correlation coefficient| > 0.8 and p < 0.001 by the Pearson correlation analysis (Schober et al., 2018). The Kaplan–Meier analysis and univariate Cox regression analysis were used to further identify differentially expressed RNA (DERNA) molecules with prognostic values in LKB1-mutant LUAD.

To reduce the risk of overfitting, the LASSO-Cox regression analysis was applied to construct a prognostic model (Tibshirani, 1997; Simon et al., 2011). The LASSO algorithm was used for variable selection, and the penalty parameter (λ) for the model was determined by a ten-fold cross validation following the minimum criteria. A multivariate regression analysis was then conducted on LKB1-mutant LUAD survival-related genes to obtain a risk gene signature, and the risk score based on the signature was formulated as follows: risk score = expX1 × coefX1 + expX2 × coefX2 +… + expXi × coefXi, where, coefXi represents the synergetic coefficient and expXi represents the relative expression of RNA. The median risk score of all samples was used as the critical value to form high-risk and low-risk groups.

Survival and Prognostic Analysis of the Prognostic Risk Model

The survival difference between the high-risk and low-risk groups was estimated by the K–M survival analysis. Time-dependent receiver operating characteristic curves were drawn, and the area under the curve (AUC) was calculated to estimate and compare the predictive accuracy for survival time by the risk score and different clinical-pathological factors. Univariate and multivariate Cox regression analyses were performed to find independent prognostic factors. The relationship between the clinicopathological characteristics and risk score was identified by a correlation analysis. A nomogram was then drawn to predict 1-, 3-, and 5-year mortality. This nomogram utilized both clinical information and risk scores because both factors have probable impacts on survival. The predictive performance of the nomogram was identified using calibration curves, K–M survival analyses, and ROC curves.

Because of the inaccessibility of other data relating to LKB1-mutated LUAD and the similarity of tumor characteristics in the same mutational background, we used patients with other cancers who had LKB1 mutations as a validation cohort. K–M survival analyses and ROC curves were used to validate the accuracy of the model.

Immune Function Analysis of the Prognostic Risk Model

The Pearson correlation analysis was performed to assess the relationship between the risk score and immune cell infiltration. For detailed information about immune functions, the single-sample Gene Set Enrichment Analysis (ssGSEA) was performed to calculate the infiltrating score of 16 immune cells and the activity of 13 immune-related pathways (Rooney et al., 2015). A Principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) were performed to explore the model’s ability to distinguish LKB1-mutant LUAD patients, based on the expression of genes.

Drug Sensitivity Analysis of the Prognostic Risk Model

To assess the relationship between the chemotherapeutic response and the risk score, we performed a drug sensitivity analysis by using the pRRophetic package based on gene expression data (Geeleher et al., 2014a; Geeleher et al., 2014b). The pRRophetic algorithm has been used extensively in medical studies (Liu et al., 2021a; Liu et al., 2021b; Liu et al., 2022a; Liu et al., 2022b; Liu et al., 2022c). We further explored the expression of targets of immune checkpoint inhibitors and other drugs between high-risk and low-risk groups to search for the potential functions of the signatures in the responses to immunotherapy and chemotherapy.

Model Specificity Analysis of the Prognostic Risk Model

We used the constructed prognostic model for wild-type LKB1 LUAD. The same analysis used in LKB1-mutant LUAD was performed in wild-type LKB1 LUAD to examine if the prognostic risk model was only applicable to LKB1-mutant cancers.

All statistical and modeling analyses were performed using Perl, GSEA (GSEA_4.1.0), and R (version 4.0.2) software.

Results

The flow chart of this study is shown in Figure 1.

FIGURE 1

FIGURE 1. Flow chart of this study.

LKB1 Mutations Suppress the Immune Landscape of Lung Adenocarcinoma

In TCGA cohort, 19% of LUAD patients were found to possess an LKB1 mutation (Figure 2A). Missense mutations and truncating mutations were the main mutation types, and these mutations were mostly located in the kinase domain (Figures 2A,B). These mutations significantly affected the expression level and function of LKB1 in cells. These findings were consistent with our hypothesis that the expression level of LKB1 is low in LKB1-mutant LUAD (Figure 2C).

FIGURE 2

FIGURE 2. From the cBioPortal database, (A) LKB1 genetic alterations (B) and mutation positions in LUAD patients of TCGA database. (C) Difference in the expression level of LKB1 between the LKB1-wt and LKB1-mu groups. (D) GSEA in LKB1-mutant LUAD. (E) ESTIMATE score, immune score, and stromal score and (F) infiltration abundance of six kinds of immune cells in the LKB1-wt and LKB1-mu groups.

Cell adhesion, natural killer cell-mediated cytotoxicity, toll-like receptor signaling pathways, and other immune-related pathways were suppressed in LKB1-mu LUAD as compared with LKB1-wt LUAD (Figure 2D). A lower stromal score, immune score, and ESTIMATE score (all p < 0.001) in the LKB1-mu group also suggested a lower degree of immune cell infiltration (Figure 2E, Supplementary Figure S1A). Because different immune cells play different roles in tumor progression, we further explored the differential degrees of infiltration of different types of immune cells. The TME of the LKB1-mu group possessed fewer B cells, CD8⁺ T cells, CD4⁺ T cells, macrophages, neutrophils, and dendritic cells (Figure 2F). Additional details regarding immune cell infiltration are given in Supplementary Figure S1B.

Of the 46 known immune checkpoint genes, 17 were found to be significantly differentially expressed between the two groups of LUAD. These differentially expressed genes included CD274, which codes for the programmed death ligand 1 (PD-L1), CD40, and CD44 (Supplementary Figure S1C). Notably, the high expression of PD-L1 LKB1-mutant LUAD is consistent with the results that demonstrate low infiltration of T cells. These results showed that the mutation of LKB1 was a common event in LUAD and that these mutations suppressed the immune landscape.

Prognostic Risk Model of a Three-Gene Signature

The differential analysis results uncovered a total of 1,869 DEmRNAs, 1,549 DElncRNAs, and 105 DEmiRNAs (Figure 3A). According to the correlation analysis, the differential expressions of 228 immune-related mRNAs were found to intersect with the differential expressions of 488 immune-related lncRNAs. The K–M analysis and univariate Cox regression analysis further identified a total of 26 mRNAs, 9 lncRNAs, and 8 miRNAs (p < 0.05, Figure 3B, Supplementary Tables S1, S2, and S3), for a total of 43 RNAs.

FIGURE 3

FIGURE 3. (A) Differentially expressed mRNAs, lncRNAs, and miRNAs between the LKB1-wt and LKB1-mu groups. (B) Kaplan–Meier survival analysis of candidate signatures. (C) Result plots of the LASSO regression analysis.

The LASSO-Cox regression analysis was applied to establish a prognostic model using the expression profiles of these 43 differentially expressed genes. A signature including three genes (CHST4, GLP1R, and hsa-miR-582-5p) was identified based on the optimal value of λ (Figure 3C). The risk score was calculated as follows: risk score = 0.0113 × C_e—0.0817 × G_e + 0.0016 × H_e, where C_e is the expression level of CHST4, G_e is the expression level of GLP1R, and H_e is the expression level of hsa-miR-582-5p. The LKB1-mutant patients were stratified into a high-risk group (n = 37) and a low-risk group (n = 37) according to the median cut-off value of the risk score.

Survival and Prognostic Analysis Using the Risk Model

The K–M curve demonstrated that patients in the high-risk group had a significantly worse OS than patients in the low-risk group (Figure 4A, p < 0.001). The 1-, 3-, and 5-year AUC of the risk score for OS were 0.843, 0.916, and 0.8, respectively (Figure 4B). Moreover, the predictive performance of the risk score (AUC = 0.845) was better than the performance of age, gender, and stage (AUC = 0.613, 0.597, and 0.699, respectively) (Figure 4B). As shown in Figure 4C, patients in the high-risk group had a higher probability of death and shorter survival time.

FIGURE 4

FIGURE 4. (A) Kaplan–Meier survival analysis in the high-risk and low-risk groups. (B) ROC curves of the risk score and clinical-pathological factors with OS. (C) Risk factor analysis and (D) independent prognostic analysis of risk scores in LKB1-mutant LUAD. (E) Heatmap and (F) proportion diagram of clinical correlation analyses.

Univariate and multivariate Cox regression analyses also indicated that the risk score could be used as an independent prognostic factor for LKB1-mutant LUAD (Figure 4D). The risk score was found to be significantly positively associated with the stage (Figures 4E,F). Figure 4F indicates that higher expression levels of CHST4 and hsa-miR-582-5p had higher stage patients. In addition, the expression level of hsa-miR-582-5p was found to be associated with T and N (p < 0.05, Supplementary Figure S1D).

To help clinicians better predict the prognosis of LUAD patients with mutated LKB1, a nomogram was constructed using the risk score plus stage (Figure 5A). The prognostic capability of the nomogram was further evaluated by ROC curves, K–M survival analyses, and calibration plots. The degree of fit of the calibration curves corresponding to 1-, 3-, and 5-year survival shows good consistency (Figure 5B). The AUC values of ROC curves corresponding to 1-, 3-, and 5-year survival in LKB1-mutant LUAD were 0.826, 0.946, and 0.819, respectively (Figure 5C). The K–M curve demonstrated that patients in the high-risk group had a significantly worse OS than patients in the low-risk group (Figure 4D, p < 0.001). These results demonstrated that the risk score and nomogram had a good prognostic ability.

FIGURE 5

FIGURE 5. (A) Nomogram of LKB1-mutant LUAD. (B) 1-, 3-, and 5-year calibration plots, (C), 1-, 3-, and 5-year ROC curves, and (D), Kaplan–Meier survival analysis of the nomogram. (E), 1-, 3-, and 5-year ROC curves, and (F), Kaplan‐Meier survival analysis of the validation cohort.

When the validation cohort was divided into high- and low-risk groups according to the constructed model, there was a significant difference in the overall survival (p = 0.005, Figure 5F). The 1-, 3-, and 5-year AUC were 0.678, 0.781, and 0.757, respectively (Figure 5E). These indicated that this prognostic model had a good survival prediction ability for LKB1-mutated cancers.

Immune Function Analysis of the Prognostic Risk Model

The Pearson correlation analysis indicated the risk score and the three other genes were significantly related to CD8⁺ T cell infiltration (R = 0.49, 0.47, 0.26, and −0.28, respectively). Among these factors, GLP1R showed a significant negative correlation (p < 0.05, Figure 6A). The risk score was associated with four immune cells and five immune-related pathways according to the ssGSEA (Figure 6B). CHST4, GLP1R, hsa-miR-582-5p, and risk scores were significantly associated with multiple immune checkpoints (Figure 6C). Of the 46 immune checkpoints, 10 were significantly differentially expressed between the high- and low-risk groups of patients with LKB1-mutant LUAD (Figure 6D and Supplementary Figure 1E). the PCA and t-SNE analysis demonstrated that the patients were well distinguished according to the risk score (Figures 6E,F).

FIGURE 6

FIGURE 6. (A) Pearson correlation analysis between CD8⁺ T cell and the risk score and three other signatures. (B) Boxplots of ssGSEA. (C) Correlation analysis and (D) boxplots of immune checkpoints. (E) PCA and (F) t-SNE analysis of risk scores.

Drug Sensitivity Analysis of the Prognostic Risk Model

We calculated the half-maximal inhibitory concentration (IC₅₀) of 138 chemotherapy drugs in the high- and low-risk groups. The associations of the values of IC₅₀ of 64 chemotherapy drugs with the groups were statistically significant (p < 0.05), and the p values associated with 15 chemotherapy drugs were less than 0.001 (Figures 7A,B and Supplementary Figure S2A). For these chemotherapy drugs, we analyzed the relationship between the risk score and the expression of genes coding for immune checkpoint-related drug targets. The expressions of AURKA, CD274, KIF11, and PLK1, which encode targets of ZM.447,439, PD-L1 immunotherapy, S-trityl-l-cysteine, and BI.2536, respectively, were associated with the risk score (Figure 7A, Figure 6D). However, the expression level and mutation status of EGFR were not associated with the risk score (Supplementary Figure S2B). The risk score was not statistically correlated with the targets of other drugs and HLA-related antigens (Supplementary Figure S2B).

FIGURE 7

FIGURE 7. (A,B) Drug sensitivity analysis and (C) expression level correlation analysis in the high-risk and low-risk groups.

Model Specificity Analysis of the Prognostic Risk Model

When we applied the same risk model to the LKB1-wt group, the K–M curve demonstrated that patients in the high-risk group had a worse OS than patients in the low-risk group (Figure 8A, p = 0.012). However, the 1-, 3-, and 5-year AUC of the risk score for the OS were 0.52, 0.603, and 0.51, respectively (Figure 8B), which showed that the predictive performance of the risk score was poor. According to univariate and multivariate Cox regression analyses, the risk score could not be used as an independent prognostic factor for LUAD patients with wild-type LKB1 (Figures 8C,D). LUAD patients with wild-type LKB1 could not be distinguished according to the risk score (Figure 8F). There was no obvious relationship between the risk score and probability of death (Supplementary Figure S3A).

FIGURE 8

FIGURE 8. Model specificity analysis of the prognostic risk model: (A) Kaplan–Meier survival analysis , (B) ROC curves , (C) univariate and (D) multivariate Cox regression analyses. (E) Nomogram of LKB1-wild-type LUAD. (F) PCA.

The same nomogram was also constructed using the risk score plus stage in LKB1-wild-type patients (Figure 8E). However, the impact of the risk score on patient prognosis was much smaller than the impact of the stage. The calibration curves corresponding to 1-, 3-, and 5-year survival rates show consistency (Supplementary Figure S3B). The AUC values of ROC curves corresponding to 1-, 3-, and 5-year survival in LKB1-wt UAD were 0.701, 0.718, and 0.691, respectively (Supplementary Figure S3C). The K–M curve demonstrated that patients in the high-risk group had a significantly worse OS than patients in the low-risk group (Supplementary Figure S3D, p < 0.001). But depending on the nomogram, the stage plays a more significant role than the risk score, resulting in the nomogram showing better predictive ability. These results demonstrated that the risk score was not applicable to wild-type LKB1 tumors.

Discussion

LUAD patients diagnosed at the same stage according to the current TNM classification system can have different clinical outcomes, due to heterogeneity at both the molecular and the genetic levels (Navin et al., 2011; Gerlinger et al., 2012). Notably, LKB1-inactivating mutations are important cancer drivers and have a significant impact on the prognosis (Tsai et al., 2014). Because of the extensive tumor suppressor function and unique effect on the immunity of LKB1, the treatment of LKB1-mutant LUAD is progressing slowly. We therefore conducted this research to explore the immune landscape, construct an immune-related prognostic risk model, and predict effective drugs for the treatment of LKB1-mutant LUAD.

As a major driver of the cold non-T cell-inflamed microenvironment, inactivating LKB1 genomic alterations cause a decrease in the infiltration of CD3⁺, CD4⁺, and CD8⁺ T cells and low expression levels of PD-L1 (Skoulidis et al., 2015; Scheel et al., 2016; Kadara et al., 2017; Skoulidis et al., 2018). These findings are consistent with our research. Mechanistically, loss of LKB1 correlates with increases in the expression of IL-1β, IL-6, CXCL7, and G-CSF that promote myeloid cell recruitment and induce epigenetic repression of STING (Koyama et al., 2016; Kitajima et al., 2019). Corresponding increases in the expression of T-cell exhaustion markers and tumor-promoting cytokines result in the inhibition of T cells under the condition of LKB1 deficiency (Koyama et al., 2016). Thus, loss of LKB1 is not only related to the resistance to PD-1 and PD-L1 inhibition, but also to the resistance to combined anti-PD-1 and anti-CTLA4 therapy with nivolumab and ipilimumab (Hellmann et al., 2018). According to a research work by Koyama et al. (2016), treatment with an IL-6-neutralizing antibody or a neutrophil-depleting antibody yielded therapeutic benefits in LKB1-deficient patients. These results indicate that LKB1 has an important impact on immune function, and further exploration is necessary.

Therefore, we endeavored to identify specific RNAs associated with immunity and the survival prognosis. In this study, we identified a three-gene prognostic model associated with immunity that was constructed using TCGA datasets to predict the OS outcomes of patients with LKB1-mutant LUAD. These same three genes have been identified as prognostic biomarkers in multiple tumor types, and have been found to play important roles in the occurrence and development of various diseases by regulating immune function (Moschovaki Filippidou et al., 2020; Zhang et al., 2020; Hu et al., 2021; Ren et al., 2022). The AUC values of ROC curves suggested that our risk model had a superior prognostic value to other clinicopathological characteristics, including age, gender, and stage. According to univariate and multivariate Cox regression analyses, the prognostic risk model that included the three immune-related RNAs was the only independent prognostic factor for patients with LKB1-mutant LUAD.

Two of the three immune-related RNAs and the risk score had significant associations with the stage. More interestingly, the expression of hsa-miR-582-5p was associated with tumor size and lymph node status, implying that the risk score is not only strongly linked to the intrinsic biological characteristics of LKB1-mutant LUAD, but also that it has a close relationship with clinicopathological characteristics. The same analysis in LKB1-wt lung adenocarcinoma showed a poor survival predictive power. It is worth noting that the importance of the risk score and stage in influencing the prognosis was opposite in nomograms created for LKB1-mutated and wild-type LUAD.

In addition to excellent prognostic prediction ability, the risk score is closely related to immune function. The upregulation of natural killer cells in the high-risk group might be because of the dual role of these cells regarding both cancer progression and the boosting of the onset of immuno-suppressant TMEs (Di Vito et al., 2019). Increased levels of CD8⁺ cells have been shown previously to be associated with better outcomes in non-small-cell lung carcinoma, and CD8⁺ cells were an independent prognostic factor in this analysis (Fehniger et al., 2003). However, in our analysis, more CD8⁺ cells and higher cytolytic activity were found in the high-risk group. More detailed research thus needs to be performed to explore the underlying mechanisms. Moreover, higher risk scores correlated with protumor immunity, including the activity of inflammation-promoting and APC co-inhibition, and these factors may provide an explanation for their poor prognosis.

We used a drug sensitivity analysis to compare the IC₅₀ values of the high-risk and low-risk groups to confirm the ability of the risk score as a treatment-guiding biomarker. Among the three EGFR tyrosine kinase-inhibitors investigated, gefitinib, erlotinib, and afatinib (BIBW2992), only afatinib was found to have a statistically significant difference in efficacy between the high- and low-risk groups. However, the expression level and mutation status of EGFR were not related to the risk score. We speculated that this disparate result may be caused by the additional activity of afatinib against Her2, which has been investigated in the context of breast cancer and other EGFR and Her2-driven cancers (Lin et al., 2012). According to the value of IC₅₀ and the correlation between drug targets and risk scores, ZM.447,439, an aurora kinase-inhibitor, is the most appropriate treatment for patients in the high-risk group. Meanwhile, we noticed that the IC₅₀ of the high-risk group was often lower than that of the low-risk group, and the expression level of PD-L1 is also higher in the high-risk group; this correlation means that patients in the high-risk group are more sensitive to some drugs despite having more malignant characteristics. This conclusion demonstrates the importance of identifying high-risk patients through risk scores to identify effective treatments and to investigate prognoses of LUAD patients with LKB1 mutations.

However, there are several limitations to our study. First, because of the lack of data on patients with LKB1-mutant LUAD, the sample size of the present study is not as large as other studies, which may influence the precision of the results. Second, the underlying molecular mechanisms that would explain why the risk score can predict the prognosis of patients with LKB1-mutant LUAD have not been fully clarified. It is necessary for us to explore the links between the risk score and immune activity by conducting experiments. In addition, large-scale and multi-center studies are needed to confirm our results.

Conclusion

The risk model we constructed provides a new promising biomarker that can provide precise clinical insights into the progression of LUAD in patients with LKB1 mutations and can serve as an efficient guide for the treatment of this challenging disease.

Data Availability Statement

Publicly available datasets were analyzed in this study. These data can be found at: https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga.

Author Contributions

All authors contributed to the study conception and design. HZ, GW, and JD carried out the design of the study, analyses of the statistics, creation of new software used in the work, and drafting of the manuscript. HZ, GW, YZ, YW, and XZ carried out the collection of the statistics and the preparation of the manuscript. All the authors read and approved the final version of the manuscript.

Funding

This study was supported by the National Natural Science Foundation of Shandong Province [Grant number: ZR2021MH192] and Jinan Clinical Medical Science and Technology Innovation Plan [Grant number 202019058].

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.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors are grateful for the public databases used in this study; they are especially grateful for the large amount of valuable datasets provided by the TCGA platform.

Supplementary Material

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

Abbreviations

AUC, area under the curve; FDR, False discovery rate; GSEA, gene set enrichment analysis; HR, hazard ratio; IC50, half maximal inhibitory concentration; ICI, immune checkpoint; LncRNA, long non-coding RNA; LUAD, lung adenocarcinoma; miRNA, microRNA; OS, overall survival; PCA, principal component analysis; ROC, receiver operating characteristic; ssGSEA, single-sample gene set enrichment analysis; TCGA, The Cancer Genome Atlas; TME, tumor microenvironment; t-SNE, t-distributed stochastic neighbor embedding.

References

Akbay, E. A., Koyama, S., Carretero, J., Altabef, A., Tchaicha, J. H., Christensen, C. L., et al. (2013). Activation of the PD-1 Pathway Contributes to Immune Escape in EGFR-Driven Lung Tumors. Cancer Discov. 3, 1355–1363. doi:10.1158/2159-8290.cd-13-0310

PubMed Abstract | CrossRef Full Text | Google Scholar

Bergethon, K., Shaw, A. T., Ignatius Ou, S.-H., Katayama, R., Lovly, C. M., McDonald, N. T., et al. (2012). ROS1 Rearrangements Define a Unique Molecular Class of Lung Cancers. Jco 30, 863–870. doi:10.1200/jco.2011.35.6345

PubMed Abstract | CrossRef Full Text | Google Scholar

Di Vito, C., Mikulak, J., Zaghi, E., Pesce, S., and Marcenaro, E., Mavilio D: NK Cells to Cure Cancer. 2019:41:101272. doi:10.1016/j.smim.2019.03.004

CrossRef Full Text | Google Scholar

Cancer Genome Atlas Research, N. (2014). Comprehensive Molecular Profiling of Lung Adenocarcinoma. Nature 511, 543–550. doi:10.1038/nature13385

PubMed Abstract | CrossRef Full Text | Google Scholar

Ding, L., Getz, G., Wheeler, D. A., Mardis, E. R., McLellan, M. D., Cibulskis, K., et al. (2008). Somatic Mutations Affect Key Pathways in Lung Adenocarcinoma. Nature 455, 1069–1075. doi:10.1038/nature07423

PubMed Abstract | CrossRef Full Text | Google Scholar

Drerup, J. M., Deng, Y., Pandeswara, S. L., Padrón, Á. S., Reyes, R. M., Zhang, X., et al. (2020). CD122-Selective IL2 Complexes Reduce Immunosuppression, Promote Treg Fragility, and Sensitize Tumor Response to PD-L1 Blockade. Cancer Res. 80, 5063–5075. doi:10.1158/0008-5472.can-20-0002

PubMed Abstract | CrossRef Full Text | Google Scholar

Fehniger, T. A., Cooper, M. A., Nuovo, G. J., Cella, M., Facchetti, F., Colonna, M., et al. (2003). CD56bright Natural Killer Cells Are Present in Human Lymph Nodes and Are Activated by T Cell-Derived IL-2: a Potential New Link between Adaptive and Innate Immunity. Blood 101, 3052–3057. doi:10.1182/blood-2002-09-2876

PubMed Abstract | CrossRef Full Text | Google Scholar

Geeleher, P., Cox, N. J., and Huang, R. S. (2014). Clinical Drug Response Can Be Predicted Using Baseline Gene Expression Levels and In Vitro Drug Sensitivity in Cell Lines. Genome Biol. 15, R47. doi:10.1186/gb-2014-15-3-r47

PubMed Abstract | CrossRef Full Text | Google Scholar

Geeleher, P., Cox, N., and Huang, R. S. (2014). pRRophetic: an R Package for Prediction of Clinical Chemotherapeutic Response from Tumor Gene Expression Levels. PLoS One 9, e107468. doi:10.1371/journal.pone.0107468

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerlinger, M., Rowan, A. J., Horswell, S., Larkin, J., Endesfelder, D., Gronroos, E., et al. (2012). Intratumor Heterogeneity and Branched Evolution Revealed by Multiregion Sequencing. N. Engl. J. Med. 366, 883–892. doi:10.1056/nejmoa1113205

PubMed Abstract | CrossRef Full Text | Google Scholar

Hellmann, M. D., Nathanson, T., Rizvi, H., Creelan, B. C., Sanchez-Vega, F., Ahuja, A., et al. (2018). Genomic Features of Response to Combination Immunotherapy in Patients with Advanced Non-small-cell Lung Cancer. Cancer Cell 33, 843–852. doi:10.1016/j.ccell.2018.03.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, B., Ma, X., Fu, P., Sun, Q., Tang, W., Sun, H., et al. (2021). miRNA-mRNA Regulatory Network and Factors Associated with Prediction of Prognosis in Hepatocellular Carcinoma. Genomics Proteomics Bioinforma. S1672-0229 (21). doi:10.1016/j.gpb.2021.03.001

CrossRef Full Text | Google Scholar

Jordan, E. J., Kim, H. R., Arcila, M. E., Barron, D., Chakravarty, D., Gao, J., et al. (2017). Prospective Comprehensive Molecular Characterization of Lung Adenocarcinomas for Efficient Patient Matching to Approved and Emerging Therapies. Cancer Discov. 7, 596–609. doi:10.1158/2159-8290.cd-16-1337

PubMed Abstract | CrossRef Full Text | Google Scholar

Kadara, H., Choi, M., Zhang, J., Parra, E. R., Rodriguez-Canales, J., Gaffney, S. G., et al. (2017). Whole-exome Sequencing and Immune Profiling of Early-Stage Lung Adenocarcinoma with Fully Annotated Clinical Follow-Up. Ann. Oncol. 28, 75–82. doi:10.1093/annonc/mdw436

PubMed Abstract | CrossRef Full Text | Google Scholar

Kitajima, S., Ivanova, E., Guo, S., Yoshida, R., Campisi, M., Sundararaman, S. K., et al. (2019). Suppression of STING Associated with LKB1 Loss in KRAS-Driven Lung Cancer. Cancer Discov. 9, 34–45. doi:10.1158/2159-8290.cd-18-0689

PubMed Abstract | CrossRef Full Text | Google Scholar

Koyama, S., Akbay, E. A., Li, Y. Y., Aref, A. R., Skoulidis, F., Herter-Sprie, G. S., et al. (2016). STK11/LKB1 Deficiency Promotes Neutrophil Recruitment and Proinflammatory Cytokine Production to Suppress T-Cell Activity in the Lung Tumor Microenvironment. Cancer Res. 76, 999–1008. doi:10.1158/0008-5472.can-15-1439

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, N. U., Winer, E. P., Wheatley, D., Carey, L. A., Houston, S., Mendelson, D., et al. (2012). A Phase II Study of Afatinib (BIBW 2992), an Irreversible ErbB Family Blocker, in Patients with HER2-Positive Metastatic Breast Cancer Progressing after Trastuzumab. Breast Cancer Res. Treat. 133, 1057–1065. doi:10.1007/s10549-012-2003-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Z., Guo, C., Li, J., Xu, H., Lu, T., Wang, L., et al. (2021). Somatic Mutations in Homologous Recombination Pathway Predict Favourable Prognosis after Immunotherapy across Multiple Cancer Types. Clin. Transl. Med. 11, e619. doi:10.1002/ctm2.619

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Z., Guo, Y., Yang, X., Chen, C., Fan, D., Wu, X., et al. (2021). Immune Landscape Refines the Classification of Colorectal Cancer with Heterogeneous Prognosis, Tumor Microenvironment and Distinct Sensitivity to Frontline Therapies. Front. Cell Dev. Biol. 9, 784199. doi:10.3389/fcell.2021.784199

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Z., Guo, C., Dang, Q., Wang, L., Liu, L., Weng, S., et al. (2022). Integrative Analysis from Multi-Center Studies Identities a Consensus Machine Learning-Derived lncRNA Signature for Stage II/III Colorectal Cancer. EBioMedicine 75, 103750. doi:10.1016/j.ebiom.2021.103750

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Z., Liu, L., Weng, S., Guo, C., Dang, Q., Xu, H., et al. (2022). Machine Learning-Based Integration Develops an Immune-Derived lncRNA Signature for Improving Outcomes in Colorectal Cancer. Nat. Commun. 13, 816. doi:10.1038/s41467-022-28421-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Z., Xu, H., Weng, S., Ren, Y., and Han, X. (2022). Stemness Refines the Classification of Colorectal Cancer with Stratified Prognosis, Multi-Omics Landscape, Potential Mechanisms, and Treatment Options. Front. Immunol. 13, 828330. doi:10.3389/fimmu.2022.828330

PubMed Abstract | CrossRef Full Text | Google Scholar

Moschovaki Filippidou, F., Kirsch, A. H., Thelen, M., Kétszeri, M., Artinger, K., Aringer, I., et al. (2020). Glucagon-Like Peptide-1 Receptor Agonism Improves Nephrotoxic Serum Nephritis by Inhibiting T-Cell Proliferation. Am. J. Pathology 190, 400–411. doi:10.1016/j.ajpath.2019.10.008

CrossRef Full Text | Google Scholar

Navin, N., Kendall, J., Troge, J., Andrews, P., Rodgers, L., McIndoo, J., et al. (2011). Tumour Evolution Inferred by Single-Cell Sequencing. Nature 472, 90–94. doi:10.1038/nature09807

PubMed Abstract | CrossRef Full Text | Google Scholar

Newman, A. M., Liu, C. L., Green, M. R., Gentles, A. J., Feng, W., Xu, Y., et al. (2015). Robust Enumeration of Cell Subsets from Tissue Expression Profiles. Nat. Methods 12, 453–457. doi:10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

Paez, J. G., Jänne, P. A., Lee, J. C., Tracy, S., Greulich, H., Gabriel, S., et al. (2004). EGFR Mutations in Lung Cancer: Correlation with Clinical Response to Gefitinib Therapy. SCIENCE 304, 1497–1500. doi:10.1126/science.1099314

PubMed Abstract | CrossRef Full Text | Google Scholar

Ren, Y., Lyu, J., Guo, Y., Yao, Y., and Hu, L. (2022). Long Noncoding RNA TUG1 Inhibits Tumor Progression through Regulating Siglec-15-Related Anti-immune Activity in Hepatocellular Carcinoma. J. Immunol. Res. 2022, 9557859. doi:10.1155/2022/9557859

PubMed Abstract | CrossRef Full Text | Google Scholar

Rooney, M. S., Shukla, S. A., Wu, C. J., Getz, G., and Hacohen, N. (2015). Molecular and Genetic Properties of Tumors Associated with Local Immune Cytolytic Activity. Cell 160, 48–61. doi:10.1016/j.cell.2014.12.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Scheel, A. H., Ansén, S., Schultheis, A. M., Scheffler, M., Fischer, R. N., Michels, S., et al. (2016). PD-L1 Expression in Non-small Cell Lung Cancer: Correlations with Genetic Alterations. Oncoimmunology 5, e1131379. doi:10.1080/2162402x.2015.1131379

PubMed Abstract | CrossRef Full Text | Google Scholar

Schober, P., Boer, C., and Schwarte, L. A. (2018). Correlation Coefficients. Anesth. Analgesia 126, 1763–1768. doi:10.1213/ane.0000000000002864

PubMed Abstract | CrossRef Full Text | Google Scholar

Shackelford, D. B., and Shaw, R. J. (2009). The LKB1-AMPK Pathway: Metabolism and Growth Control in Tumour Suppression. Nat. Rev. Cancer 9, 563–575. doi:10.1038/nrc2676

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, Q., Sun, Y., and Xu, S. (2020). LINC01503/miR-342-3p Facilitates Malignancy in Non-small-cell Lung Cancer Cells via Regulating LASP1. Respir. Res. 21, 235. doi:10.1186/s12931-020-01464-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Simon, N., Friedman, J., Hastie, T., and Tibshirani, R. (2011). Regularization Paths for Cox's Proportional Hazards Model via Coordinate Descent. J. Stat. Softw. 39, 1–13. doi:10.18637/jss.v039.i05

PubMed Abstract | CrossRef Full Text | Google Scholar

Skoulidis, F., Byers, L. A., Diao, L., Papadimitrakopoulou, V. A., Tong, P., Izzo, J., et al. (2015). Co-occurring Genomic Alterations Define Major Subsets of KRAS-Mutant Lung Adenocarcinoma with Distinct Biology, Immune Profiles, and Therapeutic Vulnerabilities. Cancer Discov. 5, 860–877. doi:10.1158/2159-8290.cd-14-1236

PubMed Abstract | CrossRef Full Text | Google Scholar

Skoulidis, F., Goldberg, M. E., Greenawalt, D. M., Hellmann, M. D., Awad, M. M., Gainor, J. F., et al. (2018). STK11/LKB1 Mutations and PD-1 Inhibitor Resistance in KRAS-Mutant Lung Adenocarcinoma. Cancer Discov. 8, 822–835. doi:10.1158/2159-8290.cd-18-0099

PubMed Abstract | CrossRef Full Text | Google Scholar

Skoulidis, F., and Heymach, J. V. (2019). Co-occurring Genomic Alterations in Non-small-cell Lung Cancer Biology and Therapy. Nat. Rev. Cancer 19, 495–509. doi:10.1038/s41568-019-0179-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Soda, M., Choi, Y. L., Enomoto, M., Takada, S., Yamashita, Y., Ishikawa, S., et al. (2007). Identification of the Transforming EML4-ALK Fusion Gene in Non-small-cell Lung Cancer. Nature 448, 561–566. doi:10.1038/nature05945

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, C.-C., Zhu, W., Li, S.-J., Hu, W., Zhang, J., Zhuo, Y., et al. (2020). FOXC1-mediated LINC00301 Facilitates Tumor Progression and Triggers an Immune-Suppressing Microenvironment in Non-small Cell Lung Cancer by Regulating the HIF1α Pathway. Genome Med. 12, 77. doi:10.1186/s13073-020-00773-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Tamborero, D., Rubio-Perez, C., Muiños, F., Sabarinathan, R., Piulats, J. M., Muntasell, A., et al. (2018). A Pan-Cancer Landscape of Interactions between Solid Tumors and Infiltrating Immune Cell Populations. Clin. Cancer Res. 24, 3717–3728. doi:10.1158/1078-0432.ccr-17-3509

PubMed Abstract | CrossRef Full Text | Google Scholar

Taube, J. M., Klein, A., Brahmer, J. R., Xu, H., Pan, X., Kim, J. H., et al. (2014). Association of PD-1, PD-1 Ligands, and Other Features of the Tumor Immune Microenvironment with Response to Anti-PD-1 Therapy. Clin. Cancer Res. 20, 5064–5074. doi:10.1158/1078-0432.ccr-13-3271

PubMed Abstract | CrossRef Full Text | Google Scholar

Tibshirani, R. (1997). The Lasso Method for Variable Selection in the Cox Model. STATISTICS Med. 16, 11. doi:10.1002/(sici)1097-0258(19970228)16:4<385:aid-sim380>3.0.co;2-3

CrossRef Full Text | Google Scholar

Topalian, S. L., Hodi, F. S., Brahmer, J. R., Gettinger, S. N., Smith, D. C., McDermott, D. F., et al. (2012). Safety, Activity, and Immune Correlates of Anti-PD-1 Antibody in Cancer. N. Engl. J. Med. 366, 2443–2454. doi:10.1056/nejmoa1200690

PubMed Abstract | CrossRef Full Text | Google Scholar

Tsai, L.-H., Chen, P.-M., Cheng, Y.-W., Chen, C.-Y., Sheu, G.-T., Wu, T.-C., et al. (2014). LKB1 Loss by Alteration of the NKX2-1/p53 Pathway Promotes Tumor Malignancy and Predicts Poor Survival and Relapse in Lung Adenocarcinomas. Oncogene 33, 3851–3860. doi:10.1038/onc.2013.353

PubMed Abstract | CrossRef Full Text | Google Scholar

Tumeh, P. C., Harview, C. L., Yearley, J. H., Shintaku, I. P., Taylor, E. J. M., Robert, L., et al. (2014). PD-1 Blockade Induces Responses by Inhibiting Adaptive Immune Resistance. Nature 515, 568–571. doi:10.1038/nature13954

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Q., Hu, B., Hu, X., Kim, H., Squatrito, M., Scarpace, L., et al. (2017). Tumor Evolution of Glioma-Intrinsic Gene Expression Subtypes Associates with Immunological Changes in the Microenvironment. Cancer Cell 32, 42–56. doi:10.1016/j.ccell.2017.06.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoshihara, K., Shahmoradgoli, M., Martínez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring Tumour Purity and Stromal and Immune Cell Admixture from Expression Data. Nat. Commun. 4, 2612. doi:10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, L., Fan, Y., Wang, X., Yang, M., Wu, X., Huang, W., et al. (2020). Carbohydrate Sulfotransferase 4 Inhibits the Progression of Hepatitis B Virus-Related Hepatocellular Carcinoma and Is a Potential Prognostic Marker in Several Tumors. Front. Oncol. 10, 554331. doi:10.3389/fonc.2020.554331

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, X., Xie, X., Zhao, Q., Zhang, L., Li, C., and Zhao, D. (2020). Potential Prognostic Value and Mechanism of Stromal-Immune Signature in Tumor Microenvironment for Stomach Adenocarcinoma. Biomed. Res. Int. 2020, 4673153. doi:10.1155/2020/4673153

PubMed Abstract | CrossRef Full Text | Google Scholar