A total of 1,288 samples were collected, including 600 LUAD samples obtained from The Cancer Genome Atlas (TCGA) database (TCGA-LUAD; www.portal.gdc.cancer.gov/). Tumor transcriptome data, clinical information and tumor somatic cell mutation data were extracted from the TCGA database. Additionally, data from the GSE72094 cohort (n=442) and GSE31210 cohort (n=246) from the Gene Expression Omnibus (GEO) database (www.ncbi.nlm.nih.gov/geo/) were used (Table I). The GSE31210 cohort served as an independent external dataset to validate the reproducibility of the risk model. Furthermore, immune-related genes were retrieved from ImmPort (www.immport.org/home) and InnateDB (www.innatedb.ca/), and the two sets of genetic information were merged.

The LUAD cell cluster, cell annotation and differentially expressed gene (DEG) data for each cluster were obtained from the Tumor Immune Single-cell Hub 2 (TISCH2) database (www.tisch.comp-genomics.org/), which was derived from the GSE131907 cohort (16) and a total of 203,298 cells from 44 patients. The Seurat package was used for clustering analysis using the FindClusters function, and cell distribution was visualized in a two-dimensional space using the RunTSNE function. R software (version 4.3.2; http://www.r-project.org/) was used to screen for marker genes in plasma cells.

PCIGs were identified as genes commonly expressed among immune-related genes and plasma cell marker genes. The R package ‘VennDiagram’ was utilized to determine the intersection, and a Venn diagram was generated to visualize the PCIGs.

Considering the aforementioned PCIGs, an exploration of potential molecular mechanisms was initiated using Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses. These analyses were performed utilizing the R packages ‘clusterProfiler’ and ‘org.Hs.eg.db’. GO and KEGG terms exhibiting a significance level of P<0.05 were visually represented using the ‘circlize’ R package. To further assess the protein-level mechanisms and associations, the STRING database (www.cn.string-db.org/) was used to elucidate the PPI relationships among the PCIGs. Subsequently, these interactions were depicted through a network diagram. Furthermore, R software was utilized to construct a histogram illustrating the core genes identified in the analysis.

The R packages ‘limma’ and ‘sva’ were used to determine the intersection of TCGA cohort and the GEO cohort to obtain the expression information of the intersecting genes in the two cohorts. Subsequently, the genes commonly expressed between the PCIGs and the intersection genes were identified. Univariate Cox hazard analysis was performed to identify prognosis-related genes using the R packages ‘survival’ and ‘survminer’ (P<0.05). After identifying prognosis-related genes, the genes were visualized by generating forest plots.

Initially, the TCGA cohort was designated as the training cohort and the GEO cohort as the testing cohort. Least absolute shrinkage and selection operator (LASSO) regression and multivariate Cox regression analysis were performed to construct a risk assessment model. LASSO regression analysis is a method used to analyze high-dimensional data and is often applied to construct regression models (17). The R package ‘glmnet’ was utilized to effectively reduce the number of genes in the final risk model. The risk scores for both the training and test cohorts were derived based on the model format. The risk score calculation followed the following formula: Risk score=coefficient (gene 1) × expression (gene 1) + coefficient (gene 2) × expression (gene 2) + coefficient (gene 3) × expression (gene 3) + … + coefficient (gene n) × expression (gene n). Risk scores exceeding the median were classified as high risk, whilst those below the median were classified as low risk. Following this classification, univariate and multivariate independent prognostic analyses were performed to assess the independence of the model compared with other clinical features.

Subsequently, the R packages ‘survminer’ and ‘survival’ were used to construct Kaplan-Meier survival curves for the training and test cohorts, which depended on survival duration and survival status. In addition, the R package ‘timeROC’ was utilized to generate a receiver operating characteristic (ROC) curve to evaluate the efficacy of risk scores in predicting 1-, 3- and 5-year survival in patients with LUAD. Based on this, a nomogram was developed to further facilitate survival prediction., considering factors such as age, risk, sex and disease stage. Furthermore, a ROC curve based on the nomogram, decision curve analysis (DCA) and calibration curve were generated to further confirm the predictive accuracy of the nomogram (18). Additionally, the GSE31210 cohort was used as an independent external cohort to validate the prognostic model.

Finally, Kaplan-Meier curves were constructed for 9 PCIGs to evaluate the prognosis of the high- and low-risk groups under different gene subgroups.

Gene Set Enrichment Analysis (GSEA). GSEA was used to evaluate the enriched pathways, and the R packages ‘limma’, ‘org.Hs.eg.db’, ‘clusterProfiler’ and ‘enrichplot’ were used to generate GSEA enrichment maps. The five most significant pathways were then delineated in the high- and low-risk groups (19).

The R package ‘ESTIMATE’ was used to analyze the stromal score, immune score and ESTIMATE score. Following this, a comparison was made between the high- and low-risk groups to identify any differences in these scores. Subsequently, the R packages ‘GSVA’ and ‘GSEABase’ were utilized to perform single-sample GSEA. This analysis enabled the elucidation of the differences in the infiltration of major immune cells between the high- and low-risk groups. The same analysis was then used to evaluate the variances in immune-related functions.

To assess the gene mutations in tumor cells in each sample, TMB data for LUAD was obtained from TCGA database (www.portal.gdc.cancer.gov/). The data was manipulated using the ‘TCGAbiolinks’ package and waterfall plots were constructed using the R programming language ‘maftools’, after which the TMB was calculated (20,21).

The association between the risk score and ICI-related gene expression levels were evaluated, and the R package ‘ggplot2’ was used for visualization. According to the tumor expression profiles downloaded from the database (http://tide.dfci.harvard.edu/), the Tumor Immune Dysfunction and Exclusion (TIDE) score could be used as a transcriptomic biomarker to predict the response to immune checkpoint blockade and to explore the potential clinical efficacy of immunotherapy. Moreover, to assess the sensitivity of several risk groups to chemotherapy drugs, the R package ‘OncoPredict’ was utilized to forecast the disparity in the treatment outcomes of chemotherapy drugs among high- and low-risk patients (22).

The clinicopathological features of patients were then obtained from TCGA cohort to assess the linkages between clinical features and the risk score, and a heatmap was generated utilizing the R package ‘ComplexHeatmap’. Moreover, the R packages ‘reshape2’, ‘tidyverse’, ‘ggplot2’, ‘RColorBrewer’ and ‘grid’ were used to construct a cyclic graph of clinical correlations. This enabled the evaluation of disparities in clinical characteristics between the high- and low-risk groups. Subsequently, the associations between each clinical characteristic and the risk score were assessed individually.

TRIzol™ reagent (Thermo Fisher Scientific, Inc.) was used to extract total RNA from normal human lung epithelial BEAS-2B cells and human LUAD H1373 and H1563 cell lines. Reverse transcription was performed using the PrimeScript™ RT Reagent Kit (Takara Bio, Inc.) at 42°C for 30 min, followed by inactivation at 85°C for 5 min. SYBR Green Master Mix (GeneCopoeia, Inc.) was used for qPCR. PCR amplification was initiated with an initial denaturation step at 95°C for 5 min, followed by 40 cycles consisting of 95°C for 15 sec (denaturation), 60°C for 30 sec (annealing) and 72°C for 30 sec (extension). A final extension step was performed at 72°C for 5 min. mRNA expression levels were normalized to GAPDH mRNA levels and the 2^−ΔΔCq method (23) was used for calculating the relative expression of mRNAs. All primers were purchased from Takara Biomedical Technology (Beijing) Co., Ltd.; Takara Bio Inc., and Table SI presents the forward and reverse primers used for RT-qPCR, including PCIGs and primers for normalization control (National Center for Biotechnology Information reference sequence, NM_002046.7).

Finally, immunohistochemical images of selected PCIGs were obtained from the Human Protein Atlas (HPA) database (www.proteinatlas.org/) (24), which is a comprehensive resource that provides high-resolution images and detailed data on protein expression in normal tissues, cancer tissues and cell lines. Based on the obtained images, the differences in the protein expression levels of selected PCIGs between lung cancer tissues and normal lung tissues were evaluated by assessing staining intensity and positive cell proportion provided in the HPA database.

Statistical analyses were performed using R, version 4.3.2. The Kruskal-Wallis test was used for comparing multiple groups, followed by Dunn's test for post hoc comparisons, with Bonferroni correction when applicable. For comparisons between two groups, the Wilcoxon rank-sum test was employed. Pearson correlation was used for evaluating the correlation between continuous variables. The Log-rank test was used to evaluate survival differences using Kaplan-Meier curves, with statistical significance defined as P<0.05. LASSO regression and Cox regression analyses were utilized to develop the prediction model. All data are presented as mean ± standard deviation. For each experiment, three replicates were performed per sample to ensure reproducibility. P<0.05 was considered to indicate a statistically significant difference.

Fig. 1 presents the flow chart of the present study. Human immune-related genes were obtained from the ‘ImmPort’ and ‘InnateDB’ databases, and the data sets were merged, identifying 2,533 genes. The gene expression profiles of 203,298 cells from 44 LUAD samples were subsequently obtained from the TISCH database (GSE131907). The cell clustering and annotation results revealed that cells with similar differential gene expression patterns were grouped into distinct clusters (Fig. 2A and B). scRNA-seq data was then utilized to perform cell clustering based on differential gene expression profiles. A total of 25 clusters were identified, and each cluster was named according to the expression of specific marker genes. For example, clusters 16 and 17 were identified as plasma cells, whilst clusters 3 and 4 were identified as macrophages. Notably, clusters 16 and 17, characterized as plasma cells, exhibited distinct profiles of DEGs (Table SII). To further elucidate the cellular origins of the DEGs identified in the GSE131907 dataset, in the detailed annotations of the relevant plasma cells were summarized to clarify their contributions to the gene expression profiles analyzed in the present study (Table SIII). Additionally, cluster 0 was identified as B cells, and 737 marker genes were obtained for the plasma cells. The study involved the analysis of the overlap between plasma cell marker genes and human immune-related genes, and 175 common DEGs were identified (Table SIV). A Venn diagram was generated to visualize the distribution of these genes (Fig. 2C), and these genes are referred to as PCIGs.

GO enrichment analysis of the DEGs revealed that biological processes were involved mainly in microtubule-based movement, humoral immune response and defense response to bacteria. The most representative cellular component terms were the collagen-containing extracellular matrix (ECM), immunoglobulin complex and external side of the plasma membrane. Among the molecular functions, signaling receptor activator activity, receptor ligand activity and glycosaminoglycan binding were the main categories (Fig. S1A and B). According to the KEGG analysis, the model may be associated with complement and coagulation cascades, arachidonic acid metabolism, hematopoietic cell lineage and linoleic acid metabolism (Fig. S1C and D). Based on the data and images acquired from the STRING database, the PPIs of the PCIGs were analyzed (Fig. S2A), and the core gene histogram demonstrated that the IL-1B gene had the highest number of adjacent nodes compared to the other top 30 core genes (Fig. S2B).

The intersection of the TCGA and GEO cohorts (GSE72094) were utilized to acquire expression data for the common genes present in both cohorts. Subsequently, genes that were commonly expressed between the PCIGs and the intersection genes were identified. Among all PCIGs sources, including primary tumors, lymph nodes, brain metastases, pleural effusions, and normal lung tissues and lymph nodes, co-expressed genes were extracted from the TCGA and GEO cohorts. In the GEO cohort, all tested genes were derived from tumor tissues of LUAD, whilst in the TCGA cohort, normal tissue samples were omitted. Therefore, the final screened PCIGs consisted only of genes differentially expressed in tumor tissues. A univariate Cox hazard analysis was then performed based on the survival data of the aforementioned co-expressed genes to identify the genes associated with prognosis (Table SV and Fig. S2C).

After LASSO regression analysis (Fig. 2D and E) and multivariate Cox regression analysis, 9 PCIGs were identified, namely, galectin-3 (LGALS3), WAP four-disulfide core domain 2 (WFDC2), leukocyte specific transcript 1 (LST1), X-box binding protein 1 (XBP1), ribosomal protein S19 (RPS19), superoxide dismutase 1 (SOD1), integrin β1 (ITGB1), CCAAT enhancer binding protein β (CEBPB) and C-C motif chemokine ligand 20 (CCL20) (Table SVI). The high-risk group was defined as individuals with risk scores above the median, whilst the low-risk group was defined as individuals with risk scores below the median (Table II). The risk assessment model was subsequently subjected to univariate and multivariate independent prognostic analysis (Table SVII). The results demonstrated that the risk score could predict patient prognosis independent of other clinical characteristics (Fig. 2F and G).

Considering the TCGA cohort as the training cohort and the GEO cohort as the test cohort, the high-risk group in the training cohort exhibited significantly worse survival than that of the low-risk group. This trend was also observed in the test cohort (Fig. 3A and B). The predictive capability of the model for patient survival in LUAD was also assessed. The risk score demonstrated a strong performance in predicting the 1-, 3- and 5-year survival rates of patients with LUAD. Specifically, in the training cohort, the area under the curve (AUC) values were 0.740, 0.715 and 0.678 at 1, 3 and 5 years, respectively. In the test cohort, the AUC values were 0.655, 0.598 and 0.748 at 1, 3 and 5 years, respectively (Fig. 3D and E). Similar results were observed for the GSE31210 cohort (Fig. 3C and F). Risk plots were used to illustrate the survival status of patients in the training and test cohorts (Fig. S3). To evaluate the predictive value of the model, a nomogram based on age, sex, stage and risk score was constructed (Fig. 4A), which showed predictive value for predicting the survival of patients with LUAD (AUC=0.727; Fig. 4B). Additionally, to confirm the precision of the predictive impact of the model, DCA and calibration curve assessments were performed. The DCA results indicated that the risk score served as a reliable predictive tool (Fig. 4C), independent of other clinical characteristics, thereby supporting its effectiveness. A calibration curve analysis combined with the risk score demonstrated a greater C-index of 0.698 (Fig. 4D) compared with the C-index of 0.6758 without the risk score (Fig. 4E). These findings indicated that the model improves the prognosis of patients with LUAD.

Additionally, survival curves were generated for each prognostic-related gene (Table SVIII). According to the survival curve analysis, the high-risk group exhibited significantly worse survival outcomes for the genes CCL20, CEBPB, ITGB1, LGALS3, RPS19 and SOD1 compared with that of the low-risk group. Conversely, the low-risk group demonstrated significantly worse survival outcomes for the genes LST1, WFDC2 and XBP1 (Fig. S4). Moreover, it was demonstrated that the survival curves for CCL20, CEBPB and LST1 intersected in the later stages. Consequently, a two-stage test was performed, and the results revealed that the P-values for all three genes were <0.05 (Table SIX), leading to the rejection of the null hypothesis (namely, there is no significant effect or relationship between gene expression and survival time). These findings indicate that, according to statistical tests, gene expression significantly impacts survival time. Furthermore, based on the PCIGs used to construct the model, a mechanistic map of the role of several PCIGs in the progression and obstruction of LUAD was drawn (Fig. 5).

GSEA revealed that the cell cycle, ECM-receptor interaction, focal adhesion, pathways in cancer and spliceosome pathways were notably enriched in the high-risk group (Fig. 6A), whilst the B-cell receptor signaling pathway, complement activation, immunoglobulin complex, immunoglobulin complex circulating, immunoglobulin receptor binding were notably enriched in the low-risk group (Fig. 6B).

Analysis of immune cell infiltration revealed that the high-risk group exhibited greater invasion of T cells activated by CD4 memory cells, resting natural killer cells, M0 macrophages, M1 macrophages and activated mast cells than the low-risk group, whilst resting dendritic cells and resting mast cells showed greater invasion in the low-risk group than in the high-risk group (Fig. 6C and D). Subsequently, variations in immune-related functions were assessed across different risk groups. Significant differences were observed in activated dendritic cells, B cells, human leukocyte antigen, immature dendritic cells, mast cells, neutrophils, T helper cells, tumor-infiltrating lymphocytes and the type II interferon response, all of which were enriched in the low-risk group (Fig. 6E). In addition, TME analysis suggested that the stromal score, immune score and ESTIMATE score were significantly greater in the low-risk group than in the high-risk group (Fig. 6F), suggesting that there may be a higher presence of tumor cells in the high-risk group.

Based on the somatic mutation data for LUAD tumors from the TCGA database, a waterfall plot of the top 20 mutated genes in the different groups was generated. In the high-risk group, the five most common mutated genes were TP53 (52%), titin (TTN; 52%), mucin-16 (MUC16; 42%), CUB and sushi multiple domains 3 (CSMD3; 44%) and ryanodine receptor 2 (RYR2; 37%). TP53 (39%), TTN (35%), MUC16 (38%), CSMD3 (32%) and RYR2 (35%) were prone to mutation in low-risk patients (Fig. 6G-H). Calculation of the TMB demonstrated that there was a significant positive correlation between the TMB and the risk score (Fig. 6I).

The association between the risk score and ICI-related gene expression was evaluated. The findings indicated a significant association: A high risk score was significantly associated with the upregulation of CD276 compared with a low risk score; however, there was no significant difference in the risk score according to the PD-1, programmed cell death protein 1 or cytotoxic T-lymphocyte associated protein 4 expression level (Fig. 7A). Subsequently, according to the results of the TIDE analysis, patients in the low-risk subgroup had significantly higher TIDE scores than those in the high-risk subgroup (Fig. 7B). The ‘oncoPredict’ R package was then used to evaluate the sensitivity of patients with LUAD in different risk groups to chemotherapy drugs. The results revealed that the high-risk group was significantly more sensitive to doramapimod (p38 MAPK inhibitor), ribociclib (CDK4/6 inhibitor), BMS-754807 (IGF-1R inhibitor) and SB216763 (GSK3β inhibitor) than the low-risk group, whilst the low-risk group was significantly more sensitive to SCH772984 (ERK inhibitor), 5-Fluorouracil (antimetabolite drug), cisplatin (platinum based chemotherapy agent) and MK-1775 (Potent Wee1 inhibitor) than the high-risk group (Fig. 7C and Table SX).

The results of the clinical correlation analysis indicated that male patients, as well as those with worse T and N stages, had significantly higher risk scores compared with female patients and those with lower T and N stages (Fig. 8A). Additionally, a cyclic graph was constructed to evaluate clinical correlation. The findings also indicated that patients with higher risk scores had significantly more advanced grade, T stages and N stages compared with those with lower risk scores (Fig. 8B). Boxplots revealed similar results, with increasing risk score associated with a significant increase in the lung cancer stage (Fig. S5).

The protein expression levels of LGALS3 (https://www.proteinatlas.org/ENSG00000131981-LGALS3), WFDC2 (https://www.proteinatlas.org/ENSG00000101443-WFDC2), SOD1 (https://www.proteinatlas.org/ENSG00000142168-SOD1), ITGB1 (https://www.proteinatlas.org/ENSG00000150093-ITGB1) and CEBPB (https://www.proteinatlas.org/ENSG00000172216-CEBPB) were compared between LUAD and normal tissues based on immunohistochemical staining images obtained from the HPA database (Fig. S6A). Furthermore, RT-qPCR revealed that LST1 and XBP1 were highly expressed in normal human lung cell lines (P<0.05); LGALS3, RPS19, SOD1, ITGB1 and CEBPB were highly expressed in both the H1373 and H1563 cell lines, with CEBPB showing the highest significance (P<0.001), and CCL20 was highly expressed in only the H1563 cell line (P<0.05) (Fig. S6B and C).