Integrated bioinformatics screening and experimental validation: construction of a LUAD prediction model based on Treg-related genes

Data acquisition

LUAD was searched in the Gene Expression Omnibus (GEO) database, and the screening criteria were to detect datasets with large genetic data and provide detailed cell type and status information. Finally, five datasets were included, and random allocation was conducted using the randomizr software package in R. These were then combined into two datasets. The test dataset of LUAD is GSE43458, GSE50081 and GSE68465. Two datasets, GSE42127 and GSE72094, are selected as the validation dataset. After removing the samples not belonging to LUAD patients, the final test dataset consisted of 699 samples and 12,477 genes, while the validation dataset included 575 samples and 15,337 genes. Batch calibration and data normalization were performed using the Sav and Limma software packages in R, and batch effect correction was carried out by calling the ComBat function (Zhang, Parmigiani & Johnson, 2020; Ritchie et al., 2015). When a gene corresponded with multiple probes, the average value was used as the final expression value.

Treg score calculations for the test dataset samples and evaluation of the relationship between Treg expression values and clinical features

According to the test dataset samples of RNA sequence expression and clinical data, CIBERSORT (https://cibersortx.stanford.edu/) was used to estimate the cancer and the proportion of 22 types of immune cells in the tissue adjacent to the carcinoma; this is a deconvolution algorithm that calculates the gene expression patterns using online tools and provides an estimated abundance of known cell types in a mixed cell population to obtain a Treg value for each sample (Waks et al., 2019). Treg expression levels in paired normal and cancerous tissues were compared. In addition, we evaluated the relationship between Treg cell expression and clinical features. To verify the prognostic value of Treg expression, “survival” package of the R software was used to compare overall survival in groups with high and low Treg expression values. The Kaplan–Meier curve was plotted, and the log-rank test was performed according to the corresponding p-value to evaluate whether the difference in survival time was significant.

Immune infiltrating cell estimation and co-expression network construction

Using the expression values of 12,477 genes from the sequencing data, CIBERSORTx was used to estimate the proportion of 22 types of tumor-infiltrating immune cells (Rusk, 2019). The R package “WGCNA” was used to construct the weighted co-expression network, which is a biological method for integrating co-expressed genes into the same module, calculating the correlation between the module and sample traits, screening the models with high weights associated with traits, and analyzing the genes in the module to identify target genes (Langfelder & Horvath, 2008). The scores of 22 kinds of immune infiltrating cells were taken as sample traits, and the optimal soft threshold power (β) was selected to construct the scale-free network when the scale-free topological index was 0.90. Next, a “dynamic tree cutting” algorithm was used to assign genes with similar expression patterns to the same modules (minimum size = 60). In addition, we estimated the correlation between the module signature genes and the infiltration levels of 22 types of immune infiltrating cells and screened the significance of the modules using Pearson’s test. Finally, we selected the “Treg” subtype that we were interested in, and conducted in-depth study on the module genes that were most correlated with Treg.

Prognostic model building

From the test dataset, 564 samples with complete gene expression profiles and survival data were obtained. LASSO analysis based on the R software package “glmnet” combined with multiple Cox regression analysis was used to screen genes significantly associated with LUAD survival data, and then these genes were applied to construct a prognostic model (McEligot et al., 2020). The coefficients calculated by multivariate Cox regression were used in the following formula: risk score = sum of coefficients × prognostic gene expression level (Lee et al., 2018). The training set samples were divided into high- and low-risk groups according to the risk scores. The Kaplan–Meier curve, receiver operating characteristic (ROC) curve, and principal component analysis (PCA) were used. PCA and risk score triptychs were used to evaluate the accuracy of the prognosis model.

Constructing a nomogram

A nomogram was established that combined clinical characteristics such as age, sex, and risk score to calculate the 1-year, 3-year, and 5-year survival of the entire LUAD group. Calibration curves were constructed to evaluate the reliability of the established nomogram (Hoshino et al., 2018).

Pathway differences in high- and low-risk groups identified by GSVA analysis

To further explore the differences in activation status of biological pathways between high- and low-risk groups, “GSVA” R was included in GSVA (Hänzelmann, Castelo & Guinney, 2013). The “c2.cp.kegg.v7.4. symbols” geneset from the MSigDB database (https://www.gsea-msigdb.org) were downloaded for GSVA analysis. GSVA is typically used to estimate changes in pathways and bioprocess activities in expression dataset samples using nonparametric and unsupervised methods.

TME, its immune landscape characterization analysis, and prediction of response to immunotherapy

The Immuno-oncology Biological Research (IOBR) package in R integrates eight commonly used algorithms (MCPcounter, TIMER, IPS, ESTIMATE, xCell, CIBERSORT, EPIC, and quanTiseq) to analyze TME-associated tumor-infiltrating immune cells (TILs) separately using 255 significant gene sets associated with tumor TME, and to evaluate the pattern of immune infiltration in each sample (Zeng et al., 2021). The TIL, TME, and immunoinfiltration scores in the high- and low-risk cohorts were calculated, and the scores of each sample were compared to draw conclusions.

Prediction of chemotherapy response

According to the median risk score, 541 samples were divided into high-risk and low-risk groups, and each group was administered 198 chemotherapy drugs, including camptothecin, axitinib, staurosporin, doramapimod, and ribociclib. The sensitivity of each sample to chemotherapy was predicted using Genomics of Drug Sensitivity in Cancer (GDSC) (http://www.cancerrxgene.org/) (Cokelaer et al., 2018). The calculations were performed using the R software package “oncoPredict”,and Pearson correlation coefficients were calculated to investigate the correlation between prognostic gene expression and drug therapy sensitivity.

Western blot analysis

The LUAD cells (A549) and normal human lung epithelial cells (BEAS-2B) were obtained from the Shanghai Institute of Biochemistry and Cell Biology (Shanghai, China) (Campbell, Main & Fitzgerald, 2013; Zhang et al., 2025; Schneider et al., 2021; Baudat et al., 2025). Three main surface markers were selected for the detection of Treg cells: FOXP3, CD4, and CD25, which have been defined for the characterization of Treg cells in mammals (Selvaraj, 2013).

A549 and BEAS-2B cells were cultured in DMEM containing 10% fetal bovine serum and 1% penicillin-streptomycin at 37 °C and 5% CO₂. The logarithmically grown BEAS-2B and A549 cells were added to 100 µL lysate (lysate = RIPA: phosphoprotease inhibitor = 1:100) and cracked on ice for 20 min. The supernatant was collected via centrifugation. Using the BCA (catalog no. P0012; Beyotime, Shanghai, China) kit, the protein concentration of the samples was roughly measured. Finally, the loading volume of each sample was uniformly set to eight µL. At the same time, the same amount of proteins were separated by 10% SDS-PAGE gel, and then transferred onto a PVDF membrane. The membrane was washed with TBST buffer and closed with 5% skim milk powder at room temperature for 120 min. Next, the PVDF membrane was incubated with antibodies (1:1,000 dilution) against FOXP3 (Ptgcn, catalog number: 22228-1-AP), SCN9A (Ptgcn, catalog number: 20257-1-AP), CD25 (Ptgcn, catalog number: 30449-1-AP), CD4 (Ptgcn, catalog number: 19068-1-AP), and GAPDH (catalog no. AB0036; Abways, San Diego, CA, USA), and incubated overnight at 4 °C. The following day, the membranes were incubated with a secondary antibody (anti-rabbit immunoglobulin) (Ptgcn, catalog no. SA00001-2) was diluted 1:5,000 and incubated at room temperature for 2 h. The membrane was washed again and tested using an Omni-ECL™  Ultra Sensitive Chemical Luminescence Assay Kit (catalog No. SQ201; Epizyme, Cambridge, MA, USA).

Statistical analysis

In this study, R (version 4.3.3) was used for statistical analysis. The Student’s t-test was employed to compare the two groups of normally distributed variables, to assess the differences between the groups, and the results were presented as mean ± standard deviation (SD). For continuous variables, the Wilcoxon test was used. The Kaplan–Meier curve was employed to compare the survival differences among different risk groups, and then the log-rank test was conducted.

Treg in LUAD

We used CIBERSORT to evaluate the value of 22 types of immune cells in our analysis, and the results showed that cancer samples had higher Treg expression levels than paracancerous tissues (P = 1.2e−09) (Fig. 1A). In addition, patients with LUAD with high Treg expression had a worse prognosis than those with low Treg expression (P = 0.00037) (Fig. 1B). These results suggest an association between Treg expression and the clinical characteristics of patients with LUAD.

Sequencing of gene data of immune infiltrating cells, and construction of co-expression networks

The test dataset of 699 samples (650 tumor samples and 49 normal samples) of 12,477 gene expression data was used to calculate the abundance of 22 types of immune cells using CIBERSORTx. Next, 22 parts of the immune infiltrating cells with gene expression data were selected as WGCNA algorithm traits, and samples with no statistical differences were removed. A co-expression network was constructed for 12,477 gene expression data from 22 types of immune infiltrating cells in LUAD samples. A suitable soft threshold power (β = 2) was selected so that the scale-free topological index reached a power value of 0.90 (Figs. 2A, 2B). Treg-related genes with similar expression patterns were merged into the same module using a dynamic cut-tree algorithm (module size = 60) to form a hierarchical clustering tree with five modules (Fig. 3). As shown in Fig. 4, among the five modules, 405 genes in the brown module were highly correlated with Tregs (R² = 0.43, P = 1e−30). Subsequently, brown module genes were selected to construct a prognostic model (Liu et al., 2021).

Construction and validation of key gene identification and prognostic model

To evaluate the prognostic performance of the brown module genes, a univariate Cox regression analysis was performed on 405 genes (Table S1, Fig. S1C). In the test dataset, 16 genes were strongly associated with OS in patients with LUAD (p < 0.001). These 16 genes were analyzed by LASSO regression, and 12 genes were screened (Figs. S1A, S1B) for multivariate Cox regression analysis, of which six genes were significantly associated with the risk ratio in patients with LUAD (Fig. S1D). These genes have been previously used to develop prognostic models. The six genes were weighted by relative coefficients, The formula was: risk score = (0.1688*ADARB2) + (0.1802*B3GALT1) + (−0.1363*FER) + (0.2407*LTB4R2) + (−0.1207*N6AMT1) + (−0.0899*SCN9A). ADARB2 (HR = 1.1838 (1.0275–1.3639), P = 0.0195), B3GALT1 (HR = 1.1974 (1.0331–1.3879), P = 0.0167), FER (HR = 0.8726 (0.7763–0.9808), P = 0.0223), LTB4R2 (HR = 1.2722 (1.1052–1.4643), P = 0.0008), N6AMT1 (HR = 0.8863 (0.8059–0.9747), P = 0.0128), SCN9A (HR = 0.9140 (0.8410–0.9933), P = 0.0341) (Table 1). The risk scores of the samples in the test dataset were calculated according to the above formula and divided into high- and low-risk groups. There was a statistically significant difference in overall survival among the groups (p < 0.0001, log-rank test) (Fig. 5A). The area under the curve (AUC) for overall survival was 0.68 at 6 years, 0.64 at 3 years, and 0.51 at 1 year (Fig. 6A). PCA revealed significant differences between low- and high-risk patients (Fig. 7A). We ranked the patients’ risk scores, and Fig. 8 reflects the life status, survival status distribution, and gene expression patterns of patients with LUAD. To explore the predictive power of the prognostic model, we constructed a complete validation dataset and calculated the risk score for the entire set using the prognostic model calculation formula. According to the risk score, the whole group was divided into high- and low-risk groups, and there was a statistically significant difference in the Kaplan–Meier survival curves between the two groups (P = 0.017) (Fig. 5B). In the entire set, the AUC at 1 year was 0.53 and the AUC at 3 and 6 years were 0.56 (Fig. 6B). PCA also showed differences among the groups (Fig. 7B). The risk score, survival status, and expression of the six prognostic genes of the whole group are shown in Fig. 9.

Prediction of survival by nomogram

Using known risk scores and clinical features, multivariate logistic regression was used to construct nomograms that could accurately predict patient survival. Age, sex, and risk scores were considered predictors of survival and were incorporated into the nomogram (Fig. 10A). It can be seen that the risk score is the most influential factor in the overall nomogram score. According to the nomogram correction curve (Fig. 10B), the survival rate predicted using the prognostic model was consistent with the observed survival rate. This indicated that the prognostic model has good clinical practicability.

Pathway differences between high- and low-risk groups identified by GSVA analysis

To explore the differences in the activation status of biological pathways between the high- and low-risk groups, samples were divided into average value according to risk score expression values, and GSVA was performed. A heat map (Fig. 11) was used to visualize the relevant biological processes and arrange them according to the p-value. The results showed that carcinogenic activation pathway of the P53 signaling pathway and cell cycle were significantly enriched in the high-risk group. These immune-related pathways regulate tumorigenesis and development, providing a theoretical basis for poor prognosis in high-risk groups. Figure 12 shows the classification of 108 pathways with significant differences, among which the pathway related to the human metabolic system accounted for the most, and the immune-related pathway accounted for 19.23% of the differences.

Risk score associated with TME signaling and immune cell infiltration

To further explore whether risk scores were associated with the TME in LUAD, we used the IOBR to define TME signaling. Tumor-related signal scores were higher in high-risk group than in the low-risk group. These signals included genes involved in cell cycle regulation, DNA damage repair, mismatch repair and homologous recombination, molecular carcinoma m6A, exosomes, positive regulation of exosome secretion, and ferritin deposition (Figs. 13A, 14). Next, we used eight typical algorithms (MCPcounter, TIMER, IPS, ESTIMATE, xCell, CIBERSORT, EPIC, and quanTiSeq) to estimate T cell infiltration, which showed lower scores for low-risk group (Fig. 15). Contrary to the T-cell invasion results, most neutrophils showed more infiltration in the low-risk group; however, this difference was generally less pronounced (Fig. 16). However, in natural killer (NK) cells, differences were widespread, and scores were higher in the high-risk group (Fig. 13B). Based on the results of the infiltration of different immune cells using various algorithms, it was concluded that the TME had a stronger immunosuppressive effect in the high-risk group.

Relationship between prognostic model and chemotherapy sensitivity

Chemotherapy is a conventional treatment for LUAD. We predicted and analyzed the responses to 198 chemotherapeutic drugs, including camptothecin, axitinib, staurosporin, etc. For patients in the high- and low-risk groups, all drugs showed significant differences between groups (Fig. 17). Figures 18A, 19A, and 20A show that paclitaxel (R² = 0.25, P = 7.7e−11), rapamycin (R = 0.28, P = 2e−13), and staurosporin (R = 0.19, P = 1.1e−06) positively correlated with the risk score. The box plot (Figs. 18B, 19B, and 20B) showed higher sensitivity to paclitaxel (P = 2.9e−05), rapamycin (P = 2.4e−07), and staurosporin (P = 0.00052) in the low-risk group, all of which suggest that chemotherapy agents may have greater efficacy in the low-risk group.

Western blot analysis verified the expression of prognostic model genes in LUAD

In the validation dataset, SCN9A was found to have the highest negative correlation with the risk score (R = −0.54, p < 2.2e−16) (Figs. 21A, 21B, 21C, 21D, 21E, and 21F). Subsequently, the proteins of A549 and BEAS-2B cells were collected for western blot analysis of the modeling gene SCN9A and Treg marker genes FOXP3, CD4, and CD25. The results of western blot analysis showed that the expression of Treg marker genes FOXP3, CD4, and CD25 in LUAD cells was upregulated compared with that in the BEAS-2B cells of normal lung tissue, whereas the expression of SCN9A was downregulated (Figs. 21G, 21H, and 21I). This is consistent with the results of previous dataset analysis. Treg expression was increased in LUAD compared with that in normal lung tissue, and the SCN9A was negatively correlated in the prognostic model; therefore, its expression was decreased in LUAD.