RNA sequencing data for a LUAD cohort (n=592) were obtained from The Cancer Genome Atlas (TCGA) via the GDC portal (portal.gdc.cancer.gov/) and processed using the GDC RNASeq analysis pipeline (v32.0; docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/Expression_mRNA_Pipeline/). Raw counts and clinical variables (age, sex, tumor stage, overall survival and vital status) were retrieved from the UCSC Xena database (18). Two independent microarray datasets (GSE31210, n=226; GSE72094, n=398) were obtained from Gene Expression Omnibus (GEO) for external validation of prognostic findings (19). For combined analyses, batch effects were removed using the ‘ComBat’ algorithm from the sva R package. Samples with missing survival data or critical clinical information (e.g., TNM stage) were excluded from corresponding analyses.

A total of 18 candidate KRT genes were screened in the TCGA dataset using two machine learning approaches. Least Absolute Shrinkage and Selection Operator (LASSO) regression (glmnet package, version 4.1–8; cran.r-project.org/web/packages/glmnet/) with 10-fold cross-validation was applied to select genes with non-zero coefficients. A random survival forest model (randomForestSRC package, version 3.2.3; http://cran.r-project.org/web/packages/randomForestSRC) was then constructed, ranking genes by mean decrease in Gini impurity, with the top 10 retained. The overlap between LASSO and random forest results yielded nine consensus feature genes. An extreme Gradient Boosting (XGBoost) model was trained on these genes to classify tumor and normal samples. Shapley Additive explanations (SHAP) values were calculated to quantify each gene's contribution to tumor-normal classification in the XGBoost model (20).

The ML.Dev.Prog.Sig framework (Mime1 package, v1.0.0; github.com/Zaoqu-Liu/Mime) was applied to construct and evaluate prognostic models for LUAD. A total of 101 candidate models, including forward stepwise Cox regression, ridge regression, random survival forest, XGBoost and gradient boosting machine (GBM), were trained on the TCGA cohort with 10-fold cross-validation. Model performance was assessed using Harrell's C-index and time-dependent receiver operating characteristic (ROC) area under the curve (AUC) at 1, 3 and 5 years. The optimal model was selected based on predictive accuracy and stability across the TCGA training set and the GEO validation cohorts (GSE31210 and GSE72094).

Immune infiltration analysis was performed using the IBOR package (21), stratifying TCGA-LUAD individuals into high- and low-risk groups. CIBERSORT was applied to estimate the relative abundance of 22 immune cell types per tumor sample (22), with differences visualized by heatmaps and boxplots. MCP-counter was used to compute absolute abundance scores for eight major cell populations, including CD8+ T cells, cytotoxic lymphocytes, monocytes, dendritic cells, neutrophils, NK cells, B cells and fibroblasts (23). Single-sample GSEA (ssGSEA) was used to calculate enrichment scores for immune-related pathways (24). Expression of key immune checkpoint genes (PDCD1, CD274, CTLA4, LAG3 and TIGIT) (25) was compared between groups. Two-sided Wilcoxon tests were used for group comparisons and P<0.05 was considered to indicate a statistically significant difference.

Differentially expressed genes (DEGs) between high- and low-risk groups were analyzed by gene set enrichment analysis using the clusterProfiler package (26). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway gene sets (c2.cp.kegg.v7.5.symbols) were tested for enrichment with adjusted P<0.05.

The Tumor Immune Dysfunction and Exclusion (TIDE) algorithm was used to estimate the likelihood of immune evasion (27). For each patient, total TIDE, T cell dysfunction and T cell exclusion scores were calculated and compared between risk groups to predict potential response to immune checkpoint inhibitors.

Drug response data for the NCI-60 cancer cell lines were retrieved from the CellMiner database (28). Pearson correlation was employed to analyze associations between the KRT risk score and chemotherapeutic drug sensitivity. Correlations with |r| >0.3 and P<0.05 were considered statistically significant (29). This analysis aimed to explore potential general relationships between KRT expression and drug sensitivity, providing hypothesis-generating evidence based on pan-cancer cell line data.

In total, five human cell lines were used: LUAD cell lines A549 (TCHu150) and H1299 (TCHu160), normal human bronchial epithelial cell line BEAS-2B (GNHu27) and lung cancer cell lines NCI-H69 (SCSP-5077) and NCI-H196 (SCSP-5088). All cell lines were obtained from the Cell Bank of the Chinese Academy of Sciences and maintained at 37°C in a humidified incubator with 5% CO₂. Detailed culture conditions are provided in Table SI. To investigate the function of KRT81, cells were transfected with either a negative control short hairpin (sh)RNA (sh-NC) or one of two shRNAs targeting KRT81 (sh-KRT81; Shanghai GenePharma Co., Ltd.). Transfections were conducted using Lipofectamine^® 3000 (Thermo Fisher Scientific, Inc.) at 37°C according to the manufacturer's instructions. Briefly, cells were incubated with the transfection complexes for 6 h, after which the medium was replaced with complete culture medium. Subsequent molecular and functional experiments were performed 48 h after transfection. Knockdown efficiency (>80%) was confirmed by fluorescence microscopy and reverse transcription-quantitative (RT-q) PCR. The target sequences and corresponding shRNA sequences are listed in Table SII.

Total RNA was extracted from cells (5×105 cells/well in 6-well plates) using RNA-easy Isolation Reagent (cat. No. R701; Vazyme Biotech Co., Ltd.) and reverse transcribed into cDNA using the HiScript III RT SuperMix (Vazyme Biotech Co., Ltd.). All procedures were performed according to the manufacturer's protocols. RT-qPCR was performed on an ABI 7500 real-time PCR system using SYBR Green chemistry (Vazyme Biotech Co., Ltd.). Thermocycling conditions were as follows: Initial denaturation at 95°C for 30 sec, followed by 40 cycles of denaturation at 95°C for 5 sec and annealing/extension at 60°C for 30 sec. Primers used were: KRT81 forward 5′-TTAAGGCACAGTATGACGACATTG-3′ and reverse 5′-TCTGGCACTTGGCATTCTCC-3′; GAPDH forward 5′-TGACTTCAACAGCGACACCCA-3′, reverse 5′-CACCCTGTTGCTGTAGCCAAA-3′. Relative gene expression was calculated by the 2-ΔΔCq method (30). All experiments were performed in triplicate and repeated three times independently.

Cells were lysed in RIPA buffer (Beyotime, China) with protease/phosphatase inhibitors. Protein concentration was determined by BCA assay (Thermo Fisher, USA). Equal amounts of protein (30 µg/lane) were separated on 10% SDS-PAGE and transferred to polyvinylidene difluoride (PVDF) membranes (Millipore, USA). Membranes were blocked with 5% non-fat milk (room temperature, 1 h), then incubated overnight at 4°C with primary antibodies: anti-KRT81 (1:1,000, cat. 11342-1-AP, Proteintech) and anti-GAPDH (1:5,000, cat. 60004-1-lg, Proteintech). After washing with TBST (containing 0.1% Tween-20), membranes were incubated with HRP-conjugated goat anti-rabbit (1:3,000, cat. A0208, Beyotime) or mouse (1:3,000, cat. A0216, Beyotime) for 1 h at room temperature. Bands were visualized using ECL reagent (Thermo Fisher) and imaged. Densitometry was performed with ImageJ (v1.53; National Institutes of Health).

Cell proliferation was assessed using the Cell Counting Kit-8 (CCK-8, Sigma). Transfected cells were seeded in 96-well plates (100 µl/well). At designated time points, 10 µl of CCK-8 reagent was added to each well and incubated for 4 h. Absorbance was measured at 450 nm with a microplate reader (Tecan Infinite; Tecan Group, Ltd.).

For colony assays, transfected cells were plated at 1×103 cells per well in 6-well plates and incubated for 14 days. Colonies were defined as cell clusters containing at least 50 cells. Colonies were fixed with 4% paraformaldehyde for 20 min at room temperature, stained with Giemsa for 15 min at room temperature, and colonies were observed using a light microscope. The number of colonies was quantified using ImageJ software (version 1.53, National Institutes of Health).

Cell apoptosis was assessed using an Annexin V-FITC/PI staining kit (Nanjing KeyGen Biotech Co., Ltd.). Cells were stained following the manufacturer's protocol and analyzed by flow cytometry (FACSCalibur, BD Biosciences). Data were processed using FlowJo software (version 10.6.2; BD Biosciences). The apoptotic rate was calculated as the sum of the percentages of early apoptotic cells (Annexin V+/PI−) and late apoptotic cells (Annexin V+/PI+).

Cell migration of A549 and H1299 cells was evaluated by a wound-healing (scratch) assay. Transfected cells were grown to confluence in 96-well plates. To minimize the influence of cell proliferation on wound closure, the medium was replaced with serum-free medium for 12 h prior to the assay. A uniform scratch was made in the cell monolayer with a 200-µl pipette tip. After scratching, the cells were washed twice with PBS to remove cell debris and then maintained in serum-free medium. Images of the wound area were captured at 0 and 24 h using a Cellomics imaging system (light microscopy). The wound width was measured and the migration rate was calculated as follows: Migration rate=(width at 0 h-width at 24 h)/Width at 0 h.

Cell migration and invasion were assessed using 6.5-mm Transwell chambers (Corning, Inc.). For migration assays, 1×104 transfected cells were seeded in the upper chamber in serum-free medium. For invasion assays, the upper chamber was pre-coated with Matrigel (Corning, Inc., 1:6 dilution). After 24 h, cells on the lower membrane were fixed, stained with 0.1% crystal violet and counted under a microscope.

Statistical analyses were conducted using R software (version 4.3.1; r-project.org/) and GraphPad Prism (version 9.0; Dotmatics). Survival curves were analyzed by the Kaplan-Meier method with log-rank testing. Correlation analyses were performed using Spearman's rank correlation coefficient. For multi-group comparisons, normality was first assessed using the Shapiro-Wilk test. Normally distributed data were analyzed using one-way analysis of variance (ANOVA) with Bonferroni post hoc test for multiple comparisons. Non-normally distributed data were analyzed using the Kruskal-Wallis test followed by Dunn's test with Bonferroni correction. Two-group comparisons were performed using two-sided Student's t-test or Mann-Whitney U test, as appropriate. Data are expressed as mean ± SD. P<0.05 was considered to indicate a statistically significant difference.

Differential expression analysis of 526 tumor and 59 normal lung samples from TCGA-LUAD identified KRT genes with |log2FC| > 1 and false discovery rate (FDR) <0.05. In total, 18 KRT genes were significantly dysregulated in LUAD (Fig. 1A). Two feature selection methods, LASSO regression and random forest analysis, were applied. LASSO with 10-fold cross-validation identified 12 predictive genes (Fig. 1B and C), while the random forest model selected the top 10 genes ranked by mean decrease in Gini impurity (Fig. 1D). The overlap between methods yielded a nine-gene signature (KRT1, KRT4, KRT27, KRT79, KRT80, KRT16, KRT15, KRT81 and KRT83).

The gene set was assessed using multiple classification algorithms. XGBoost achieved the highest accuracy, with an AUC of 0.996 in TCGA (Fig. 1E). External validation in GSE31210 confirmed robust diagnostic performance (AUC=0.860; Fig. 1F), with KRT81 demonstrating particularly strong discriminative ability. These results identify a small subset of KRT genes with high diagnostic potential for LUAD, suggesting that this signature can reliably distinguish tumors from normal lung tissue.

To elucidate the contribution of individual genes to model decision-making, SHAP analysis was performed on the final XGBoost classifier. The SHAP summary plot (Fig. 2A) ranked KRT genes by their overall effect on prediction, revealing that KRT4, KRT79, KRT27, KRT16, KRT80 and KRT81 exerted the strongest influence on the model output. The distribution of SHAP values indicated that higher expression of these genes predominantly shifted the prediction toward tumor classification.

A representative SHAP waterfall plot (Fig. 2B) further illustrates how gene-specific contributions accumulate to determine the final predicted probability for an individual sample. In this example, positive SHAP values for KRT4, KRT79, KRT27, KRT16, KRT81 and KRT80 elevate the prediction toward the tumor class, whereas genes with minimal effect remain close to baseline expectations.

These SHAP analyses demonstrated that a small subset of KRT genes consistently drives the classifier's decisions. This pattern aligns with their biological dysregulation in LUAD, highlighting their mechanistic relevance to the model's discriminative performance.

To identify keratin family members with prognostic significance, univariate Cox regression analysis was performed in the TCGA-LUAD cohort, which identified four KRTs significantly associated with overall survival, including KRT27, KRT16, KRT80 and KRT81, with KRT81 exhibiting the highest hazard ratio (Table SIII). These prognostic-relevant genes were subsequently incorporated into a systematic modeling framework. Among the 101 tested models, the combination of forward stepwise Cox regression with Ridge regression achieved the highest mean C-index (0.65 across cohorts; 0.66 in validation; Fig. 3A) and was selected as the final prognostic model. This selection was based on its optimal balance of predictive accuracy and stability: While achieving a top-tier C-index, the forward stepwise approach ensured model parsimony by selecting only the most informative genes and the incorporation of Ridge regression helped stabilize coefficient estimates by effectively handling multicollinearity among the keratin genes. The corresponding regression coefficients were used to derive the prognostic risk score formula: Risk Score=−0.320×KRT27+0.159×KRT80+0.238×KRT16+0.120×KRT81. Using this model, the present study calculated a risk score for each patient and classified them into high- and low-risk groups at the median. Kaplan-Meier analysis confirmed that high-risk patients had significantly worse overall survival: in TCGA (HR=1.99; 95% CI:1.39–2.85; Fig. 3B) and in both validation cohorts (GSE31210: HR=2.47; 95% CI:1.21–5.06; Fig. 3C; GSE72094: HR=2.22, 95% CI:1.53–3.23; Fig. 3D). Time-dependent ROC curves further demonstrated the model's strong predictive power: 5-year AUCs were 0.61, 0.7 and 0.7 in TCGA, GSE31210 and GSE72094, respectively (Fig. 3E). These results demonstrate the model's robust performance and potential for individualized risk stratification in LUAD.

Clinical relevance of the risk score was evaluated in the TCGA cohort (n=500). Univariate Cox regression identified tumor T stage, N stage, overall stage and the KRT risk score as significant predictors of poor survival Fig. 4A), with the risk score exhibiting the largest effect (HR=2.44). In multivariate analysis, the risk score remained an independent prognostic factor (HR=2.16; 95% CI:1.43–3.28), while T and N stages lost significance (Fig. 4B). Overall stage was significant (HR=1.39), suggesting that the risk score offers complementary prognostic value. A nomogram was then developed, integrating the risk score and clinical stage to predict 1-, 3- and 5-year overall survival (Fig. 4C). The nomogram demonstrated moderate discrimination (AUCs=0.706, 0.730 and 0.654 at 1, 3 and 5 years; Fig. 4D), with calibration curves showing good agreement between predicted and observed outcomes (Fig. 4E). These findings highlighted the potential of the KRT risk score to enhance individualized survival predictions in LUAD.

The tumor immune microenvironment was further investigated in relation to the KRT risk score. CIBERSORT analysis quantified the relative abundance of 22 immune cell subsets, revealing that high-risk tumors had significantly lower levels of CD4+ T cells, B cells, resting dendritic cells and resting mast cells, while exhibiting higher proportions of M0 macrophages compared with low-risk tumors (Fig. 5A and B). Similarly, MCP-counter analysis of eight major immune and stromal populations showed a marked reduction in cytotoxic lymphocytes (including CD4+ T cells), B cells and dendritic cells, along with increased infiltration of endothelial cells in high-risk tumors (Fig. 5C). These results suggested that the high-risk group is characterized by a more immunosuppressive tumor microenvironment, marked by reduced effector immune cell infiltration and an expanded stromal compartment.

Immune-related pathways were further assessed using ssGSEA (Fig. 5D). High-risk tumors exhibited significantly lower enrichment of pathways involved in antigen presentation and immune activation, reflecting impaired antitumor immune responses. Moreover, immune checkpoint gene profiling (Fig. S1) revealed downregulation of several inhibitory receptors (such as BTLA, BTNL9, CTLA4, BTN2A2 and BTNL3) alongside upregulation of co-inhibitory ligands (PVR, CD276 and VTCN1) in high-risk tumors. This altered checkpoint expression landscape suggests an adaptive immune-escape mechanism. These immune landscape analyses reinforce the conclusion that the KRT high-risk signature is associated with an immune-cold, exclusionary tumor microenvironment.

GSEA of DEGs between the two risk groups revealed distinct biological phenotypes. High-risk tumors were enriched in pathways related to proliferation and invasion, including the cell cycle, DNA repair, Extracellular Matrix (ECM)-receptor interaction, focal adhesion and p53 signaling (Fig. 5E). By contrast, low-risk tumors were enriched in immune-related pathways such as asthma, systemic lupus erythematosus, Graves' disease and allograft rejection, indicating a more active immune environment. These enrichments align with the aggressive behavior of high-risk tumors and the robust immune activity in low-risk tumors.

To predict response to immune checkpoint blockade, the TIDE algorithm was applied. High-risk patients had significantly higher overall TIDE scores (Fig. 5F) and higher T cell exclusion scores (Fig. 5G) compared with low-risk patients. This suggested that high-risk tumors were more prone to T cell exclusion and dysfunction, potentially leading to poorer responses to immunotherapy. In summary, the KRT risk signature may identify patients less likely to benefit from immune checkpoint inhibitors.

Using the NCI-60 pan-cancer cell line panel from the CellMiner database, associations between the KRT signature and chemotherapeutic drug sensitivity were explored. Across this diverse collection of cancer cell lines, the KRT risk score showed significant negative correlations with sensitivity to several chemotherapeutic agents (Fig. S2), suggesting that elevated expression of the KRT signature was associated with chemotherapy resistance. These results implied that patients with a high KRT risk score may have reduced responsiveness to conventional chemotherapy, highlighting the need for alternative therapeutic approaches.

Given its prognostic and diagnostic significance, KRT81 was selected for experimental validation. KRT81 expression was relatively high in LUAD cell lines A549, H1299 and NCI-H460 (Fig. 6A), as demonstrated by RT-qPCR and western blotting A549 and H1299 cells were chosen for further studies and shRNA-mediated KRT81 knockdown was confirmed by RT-qPCR (Fig. 6B) and western blotting (Fig. 6C). KRT81 silencing significantly decreased cell viability in CCK-8 assays (Fig. 6D) and markedly increased apoptosis in both cell lines (Fig. 6E and F). Colony formation assays revealed that KRT81 knockdown severely impaired clonogenic growth of A549 and H1299 cells (Fig. 7A). Additionally, scratch wound-healing (Fig. 7B) and Transwell (Fig. 7C and D) assays demonstrated that KRT81 depletion significantly suppressed cell migration and invasion. These functional experiments indicate that KRT81 promotes LUAD cell proliferation, survival and motility, reinforcing its role as an oncogenic driver and a potential therapeutic target in LUAD.