Machine learning-based prognostic model of lactylation-related genes for predicting prognosis and immune infiltration in patients with lung adenocarcinoma

Abstract

Background: Histone lactylation is a novel epigenetic modification that is involved in a variety of critical biological regulations. However, the role of lactylation-related genes in lung adenocarcinoma has yet to be investigated.

Methods: RNA-seq data and clinical information of LUAD were downloaded from TCGA and GEO datasets. Unsupervised consistent cluster analysis was performed to identify differentially expressed genes (DEGs) between the two clusters, and risk prediction models were constructed by Cox regression analysis and LASSO analysis. Kaplan-Meier (KM) survival analysis, ROC curves and nomograms were used to validate the accuracy of the models. We also explored the differences in risk scores in terms of immune cell infiltration, immune cell function, TMB, TIDE, and anticancer drug sensitivity. In addition, single-cell clustering and trajectory analysis were performed to further understand the significance of lactylation-related genes. We further analyzed lactate content and glucose uptake in lung adenocarcinoma cells and tissues. Changes in LUAD cell function after knockdown of lactate dehydrogenase (LDHA) by CCK-8, colony formation and transwell assays. Finally, we analyzed the expression of KRT81 in LUAD tissues and cell lines using qRT-PCR, WB, and IHC. Changes in KRT81 function in LUAD cells were detected by CCK-8, colony formation, wound healing, transwell, and flow cytometry. A nude mouse xenograft model and a KrasLSL-G12D in situ lung adenocarcinoma mouse model were used to elucidate the role of KRT81 in LUAD.

Results: After identifying 26 lactylation-associated DEGs, we constructed 10 lactylation-associated lung adenocarcinoma prognostic models with prognostic value for LUAD patients. A high score indicates a poor prognosis. There were significant differences between the high-risk and low-risk groups in the phenotypes of immune cell infiltration rate, immune cell function, gene mutation frequency, and anticancer drug sensitivity. TMB and TIDE scores were higher in high-risk score patients than in low-risk score patients. MS4A1 was predominantly expressed in B-cell clusters and was identified to play a key role in B-cell differentiation. We further found that lactate content was abnormally elevated in lung adenocarcinoma cells and cancer tissues, and glucose uptake by lung adenocarcinoma cells was significantly increased. Down-regulation of LDHA inhibits tumor cell proliferation, migration and invasion. Finally, we verified that the model gene KRT81 is highly expressed in LUAD tissues and cell lines. Knockdown of KRT81 inhibited cell proliferation, migration, and invasion, leading to cell cycle arrest in the G0/G1 phase and increased apoptosis. KRT81 may play a tumorigenic role in LUAD through the EMT and PI3K/AKT pathways. In vivo, KRT81 knockdown inhibited tumor growth.

Conclusion: We successfully constructed a new prognostic model for lactylation-related genes. Lactate content and glucose uptake are significantly higher in lung adenocarcinoma cells and cancer tissues. In addition, KRT81 was validated at cellular and animal levels as a possible new target for the treatment of LUAD, and this study provides a new perspective for the individualized treatment of LUAD.

Keywords: Immune; KRT81; Lactylation; Lung adenocarcinoma; PI3K-AKT.

Fig. 1

Identification of lactylation prognosis-associated genes in lung adenocarcinoma. A 52 overlapping genes were identified as lactylation-associated DEGs. B Expression of 26 lactylation prognosis-associated genes in lung adenocarcinoma in Tumor and Normal. C 120 out of 390 patients with LUAD showed genetic alterations in prognosis-associated lactylation genes. D CNV mutations are widespread in 26 prognosis-associated lactylation genes. The column indicates the frequency of changes. Loss, green dots; GAIN, pink dots. E PPI network showing interactions between proteins encoded by lung adenocarcinoma lactylation prognosis-associated genes. F Locations of CNV alterations in lactylation prognosis-related genes on chromosome 23

Fig. 2

Identification of subgroups of LACAGs in LUAD. A Unsupervised consensus clustering analysis of LUAD. B Clinical survival analysis of C1 and C2 clusters. C Graph of PCA analysis of C1 and C2 clusters. D Characterization of clinicopathology and expression of LACAGs in LUAD subtypes

Fig. 3

Construction of risk modeling. A, B LASSO regression analysis of selected prognostic genes. C Alluvial diagram showing the relationship between survival status, DEGs clusters, and risk scores. D Difference-in-difference analysis of cluster risk scores. E ROC curves of training set risk score plots, survival status plots, and 1-, 3-, and 5-year risk scores. F ROC curves for the validation set risk score plot, survival status plot, and 1-, 3-, and 5-year risk scores

Fig. 4

Independent prediction of risk models and nomogram construction. A,B Multivariate COX regression analysis of training and validation sets. C Nomogram survival prediction of LUAD patients with risk scores. D Calibration curves assessing 1-, 3- and 5-year OS consistency. E ROC curves for common clinical parameters and risk scores

Fig. 5

Analysis of tumor immune microenvironment and immune infiltration. A StromalScore, ImmuneScore, and EstimatedScore between high and low-risk groups. B Correlation analysis of model genes with immune cells. C Differences in immune-related functions between high and low-risk groups. D Correlation between risk score and immune cells

Fig. 6

Tumor mutation load analysis. A Expression differences of 26 LACAGs between high and low-risk groups. B Waterfall plots showing the top 15 mutated genes of LUAD in the high-risk group (159 samples) and low-risk group (165 samples). C, D Differences in TMB between high- and low-risk groups and correlation analysis between risk scores and TMB

Fig. 7

Single-cell RNA-seq analysis in tumor and normal tissues. A Umap was used to downscale data from tumor tissues and normal tissues and annotate each cluster. B Expression of 10 lactylation-related genes in tumor tissues in different cells. C Scoring of lactylation-related genes in tumor tissues and normal tissues. D The 3 cell subtypes of the B cell differentiation process are shown. E Expression of MSA1 in the differentiation stages of the 3 B cell subtypes

Fig. 8

A Detection of lactate expression in HBE, A549, H1299, PC9 and H1975. B Determination of lactate concentration in human lung adenocarcinoma tissues and paired adjacent normal tissues (n = 6). C Glucose uptake capacity was measured in HBE, A549 and H1975. D The knockdown efficiency of LDHA in A549 and H1975 cells transfected with siRNA1 / 2 / 3 was detected using qRT-PCR. nc: no siRNA infection; siNC: negative control. E,F Proliferative capacity was detected by CCK8 and clone formation assay. G Transwell assay was used to detect the migration and invasion number of A549 cells and H1975 cells transfected with siNC and siRNA1

Fig. 9

Expression of KRT81 in LUAD cells and tissues. A We retrieved the expression profiles of KRT81 in 483 lung adenocarcinoma tissues versus 347 normal tissues on the GEPIA2 website. B KM curves of high and low KRT81 expression groups. C, D qRT-PCR and WB were performed to detect KRT81 expression in lung normal epithelial cell lines HBE and LUAD cell lines (A549, H1299, PC9, H1975). E WB detection of KRT81 expression in 12 pairs of LUAD tissues. F IHC detection of KRT81 expression in LUAD tissues. **p < 0.01, ***p < 0.001

Fig. 10

KRT81 knockdown inhibits lung adenocarcinoma cell proliferation, migration, and invasion. A The knockdown efficiency of KRT81 in A549 cells and H1975 cells transfected with siRNA1 / 2 / 3 was determined by qRT-PCR. NC: no siRNA infection; siNC: negative control. B, C CCK8 and clone formation assay was used to detect the proliferation ability. D Transwell assay was used to detect the number of migrating and invading A549 cells and H1975 cells transfected with siNC and siRNA2. E Wound healing assay was used to detect the migration rate. WB assay was used to detect the expression levels of migration-associated proteins (MMP2, MMP9) in A549 cells transfected with siNC and siRNA2. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001

Fig. 11

Knockdown of KRT81 leads to G0/G1 phase block transition and apoptosis. A The effect of KRT81 knockdown on the cell cycle of A549 cells was detected by Cell Cycle and Apoptosis Analysis Kit. B The effect of KRT81 knockdown on apoptosis in A549 cells was detected by Annexin V Apoptosis Detection Kit. C WB detection of the expression levels of apoptosis-related proteins (Caspase3, P53, Bcl-2) in A549 cells. *p < 0.05, ***p < 0.001, ****p < 0.0001

Fig. 12

KRT81 promotes EMT and PI3K/AKT pathway activation. A WB detection of EMT-associated proteins (E-cadherin, N-cadherin, and Vimentin) after KRT81 knockdown. B, C KEGG and GSEA enrichment analysis of up-regulated DEGs. D WB detection of p-PI3K, PI3K, p-AKT, and AKT proteins levels in A549 cells after siNC, siRNA2 transfections.*p < 0.05, **p < 0.01

Fig. 13

KRT81 promotes in vivo tumorigenesis in nude mouse model and lung in situ model. A Xenograft tumors. B Tumor size and weight. C Expression levels of proliferative proteins KRT81, Ki67 and PCNA detected by IHC. D B6- 178 KrasLSL-G12D in situ lung cancer mouse model was constructed, and the expression of KRT81 was detected by IHC. *p < 0.05, **p < 0.01