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.

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.

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.

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.

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.

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.

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.

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).

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.

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.

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.

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.

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.

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.

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).

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).

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).

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.