All specimens were collected from patients undergoing surgery in the Second Xiangya Hospital of Central South University from May 2020 to May 2021. Specimens were stored at -80°C until analysis. Inclusion criteria were 1) pGGO < 1 cm in diameter detected with high-resolution CT (HRCT), 2) the patient underwent surgical resection and pathological analysis for clinical decision-making, and tumors were staged according to the American Joint Commission on Cancer (AJCC) 8th edition TNM staging system, and 3) postoperative pathological diagnosis confirmed LUAD. Finally, 30 paired samples of pGGO and adjacent normal tissue were sent to BGI Tech SOLUTIONS (Hongkong) for high-throughput transcriptome sequencing. The clinical characteristics of the patients with pGGO are summarized in Table 1 .

The high-throughput transcriptome sequencing dataset was normalized using the “edgeR” package in R (23). Differentially expressed genes (DEGs) between pGGO and adjacent normal tissue samples were obtained using the ‘‘Limma’’ package in R (|log FC|> 1, FDR P < 0.05) (24). pGGO-related DEGs were exhibited in a heatmap and volcano plot, which were generated by the “pheatmap, ggrepel, dplyr’’ package in R (25). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were conducted to identify the function of the DEGs.

Next, we extracted the expression data of the pGGO-related DEGs in TCGA-LUAD stage I-II datasets. We also performed the non-negative matrix factorization (NMF) method to cluster the LUAD stage I–II patients. NMF is an unsupervised learning technique for dimension reduction that decomposes a large measurement matrix into two low-rank non-negative matrices (26). The cophenetic correlation coefficient, based on the consensus matrix and proposed by Brunet et al. (2004), is used to measure the stability of clusters (27). NMF can classify samples better than consensus clustering. mRNA expression profiles and the clinical data for 497 patients with LUAD (Stages I–II, n=347; Stages III–IV, n=121) were downloaded from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/); follow-up information was available for each patient. The ID of the pGGO-related DEGs was used to extract expression data from the TCGA-LUAD cohort. The expression of the pGGO-related DEGs was verified in the TCGA-LUAD cohort. NMF was used to cluster patients with Stage I–II LUAD in the TCGA cohort. The clustering effect was evaluated with progression-free survival (PFS) and overall survival (OS).

Differences in pGGO-related DEGs between NMF subgroups were obtained using the “Limma” package in R. DEGs that were differently expressed between pGGOs and normal adjacent tissues, as well as differentially expressed between NMF subgroups, were identified. GO and KEGG pathway enrichment analyses were conducted to identify the function of the DEGs (28). The immune cell content of each NMF subgroup was analyzed using the “MCPcounter” package in R (29). HLA expression was compared between NMF subgroups. Based on a previous study, six immune subtypes were defined according to immune infiltrates. The relationship between the six immune subtypes and NMF subgroups was explored (30).

pGGO-related DEGs were screened with LASSO cox regression and support vector machine—recursive feature elimination (SVM-RFE) to identify prognostic marker genes (31). Feature selection for multiclass classification problems is challenging in machine learning. Existing multi-class gene selection algorithms are often not Pareto optimal. In this study, two machine learning algorithms, LASSO cox regression, and SVM-RFE, were used to achieve Pareto optimality (32, 33). LASSO cox regression is used for data dimensionality reduction and feature selection. The regression coefficient is penalized by L1, some coefficients are shrunk to zero, features with non-zero regression coefficients are selected, and 10-fold cross-validation is used to evaluate the prediction model (34). SVM-RFE is often used for gene selection. SVM-RFE ranks features from most important to least, and least important features are iteratively eliminated (35, 36).. Venn analysis was used to identify 15 pGGO-related DEGs that were shared between the LASSO cox regression and SVM-RFE machine learning algorithms.

Quantitative real-time reverse transcription–polymerase chain reaction (qRT-PCR) was used to validate the expression of the pGGO-related DEGs with a prognostic value. RNA was extracted from 24 pGGO samples, and qRT-PCR was performed for 7 pGGO-related DEGs. Primer sequences were designed from Primer3web (https://primer3.ut.ee) ( Table 2 ). qRT-PCR was performed with a SYBR Green super-mix reagent, and β-actin was the internal reference gene. The relative change in gene expression was calculated using the 2ΔΔCt method. Results are presented as the mean of 3 replicates.

Multivariate Cox regression was used to calculate regression coefficients for the pGGO-related DEGs with prognostic values. A prognostic signature was constructed using the following formula: risk score = coefficient(gene1) * exp(gene1) +…+ coefficient (gene n) * exp (gene n). A deviation plot was constructed to show the expression profile of the 15 pGGO-related DEGs with prognostic values. Patients in the TCGA-LUAD cohort were stratified into a high-risk or low-risk group using the median risk score as a cut-off. The 15-gene risk signature was verified in the GSE50081 and GSE72094 datasets, which were downloaded from the Gene Expression Omnibus–NCBI (https://www.ncbi.nlm.nih.gov/geo/). Kaplan–Meier survival curves, decision curve analysis (DCA), risk plots, and time-dependent receiver operating characteristic (ROC) curves were used to investigate the predictive power of the prognostic signature in the TCGA-LUAD cohort and GEO datasets (37, 38). The relationship between clinicopathologic factors, immune scores, clusters, and risk scores was examined using Pearson’s chi-squared test (39).

Univariate and multivariable Cox proportional hazard models were used to analyze independent associations between clinical outcomes. A nomogram based on clinicopathologic factors and risk scores was constructed. A calibration curve displayed the nomogram’s predictive power. Gene set enrichment analyses (GSEAs) in low-risk and high-risk groups were performed using the “clusterProfiler, enrichplot, DOSE, org.Hs.eg.db “package in R (40). Enrichment scores (ESs) that represent the degree to which a gene set is overrepresented at the top or bottom of a ranked list were calculated (https://www.gsea-msigdb.org/gsea/index.jsp). Tumor mutation burden (TMB) data were retrieved for the TCGA-LUAD cohort. The relationship between TMB and the 15-gene risk signature was evaluated using the “reshape2” package in R. The immunophenoscore (IPS) in patients with LUAD was downloaded from the Cancer Immunome Database (TCIA) (https://tcia.at/home) (41). The relationship between IPS and the 15-gene risk signature was evaluated with Pearson correlation analysis.

All statistical analyses were conducted with R software (Version 4.0.1).

NMF algorithms were performed using the “NMF, survival” package (42). Kaplan–Meier survival curves were plotted using the “survminer” package. SVM-RFE was performed using the “e1071, kernlab, caret” package (43). The nomogram was constructed using the “rms” package. The predictive ability of the pGGO-related DEGs signature was evaluated using the “SurvivalROC” package. The deviation plot was constructed using the “ggpubr” package. Correlation analysis was performed using the Pearson method. Differences between subgroups were evaluated using the Wilcoxon rank-sum test. mRNA profiles from the GEO datasets were normalized using the “sva, Limma” package. All figures were plotted using the “ggplot” package. For the GO and KEGG analyses, a P-value <0.05 and FDR q-value <0.25 were considered statistically significant.

A total of 1,734 DEGs (|log FC| > 2, P < 0.05) between pGGO and adjacent normal tissue samples were identified; of these, 648 DEGs were upregulated, and 1,086 DEGs were downregulated ( Figures 1A, B and Table S1 ). GO and KEGG pathway enrichment analyses suggested that the DEGs were mainly involved in immune-related pathways, such as immune response-regulating signaling pathways and lymphocyte-meditated immunity pathways ( Figure 1C ). Next, we extract the mRNA expression data of the pGGO-related DEGs in TCGA-LUAD stage IA datasets and perform the NMF consensus cluster analysis. The best cluster number was chosen as the coexistence correlation coefficient K value = 2 ( Figure 1D ); therefore, patients with Stage I–II LUAD in the TCGA cohort (n = 347) were divided into two clusters ( Figure 1E and Table S2 ). Patients in Cluster 1 had better PFS and OS compared to patients in Cluster 2 ( Figures 1F, G ).

A total of 208 pGGO-related DEGs between Cluster 1 and Cluster 2 were identified; of these, 33 pGGO-related DEGs were upregulated, and 175 pGGO-related DEGs were downregulated ( Figure 2A and Table S3 ). KEGG pathway enrichment analysis suggested that the pGGO-related DEGs were mainly enriched in immune-related pathways, including the cytokine−cytokine receptor interaction pathway and IL-17 signaling pathway ( Figure 2B ). GO function analysis suggested that the pGGO-related DEGs were mainly enriched in immune-related processes, including the humoral immune response process and immunoglobulin receptor-binding process ( Figure 2C ). HLA gene expression data showed that HLA-L, HLA-DQA2, HLA-DQB2, and HLA−DRB6 were highly expressed in Cluster 1 ( Figure 2D ). A Sankey diagram showed a relationship between the cluster subtype and six immune subtypes defined in a previous study ( Figure 2E and Table S4 ). Violin plots suggested that Cluster 2 was characterized by a high expression of B-cell lineages, endothelial cells, cytotoxic lymphocytes, CD8 T cells, monocyte-lineage cells, fibroblasts, and NK cells; meanwhile, Cluster 1 was characterized by a high expression of neutrophils ( Figures 2F–M ).

Verification of pGGO-Specific DEGs Using Machine Learning Algorithms.

Lasso cox regression and SVM-RFE identified the 23 ( Figures 2N, O and Table S5 ) and 19 ( Figure 2P and Table S5 ) most representative prognostic pGGO-related DEGs, respectively, from among the 208 pGGO-related DEGs between Cluster 1 and Cluster 2. Venn analysis was used to identify 15 DEGs that were shared between the machine learning algorithms ( Figure 3A ). These included DKK1, NPAS1, AL357143.1, KCNV1, AC068228.1, AC239859.6, and FGF5 which were highly expressed in pGGO samples, and AC087763.1, PCP2, FAIM2, AL357143.1, AC022148.1, AC021678.2, LSP1P2, and LINC00563 that were highly expressed in adjacent normal tissue samples ( Figure 3B ). This expression profile was validated by qRT-PCR ( Figure S1 ). Regression coefficients were used to identify a 15-gene risk signature. A forest plot showed that FGF5, AC239859.6, NPAS1, KCNV1, AC068228.1, and AL353746.1 were risk factors in LUAD; while DKK1, LINC00563, AC021678.2, AL357143.1, AC022148.1, LSP1P2, AC087763.1, PCP2, and FAIM2 were protective factors in LUAD ( Figure 3C ).

The TCGA-LUAD cohort was stratified into a high-risk group and low-risk group (n=234) based on their median risk score ( Table S6 ). Kaplan–Meier survival analysis showed that patients in the high-risk group had worse OS than patients in the low-risk group (P<0.001) ( Figure 3D ). The 15-gene risk signature was predictive of patients with Stage IA LUAD in the TCGA-LUAD cohort. This was expected as pGGO is recognized as a component of TNM Stage IA LUAD according to the AJCC 8th edition TNM staging system. Patients with Stage IA LUAD and high-risk scores in the TCGA-LUAD cohort had worse OS than patients with Stage IA LUAD and low-risk scores (P<0.001) ( Figure 3E ). There was no significant difference in OS between patients with Stage IA LUAD and high-risk scores and patients with Stage II LUAD ( Figure 3F ). To confirm the feasibility of the 15-gene risk signature, information for patients with LUAD from the GSE50081 and GSE72094 datasets was downloaded, combined, and stratified into a high-risk group (n=263) and low-risk group (n=262) based on the median risk score from the TCGA-LUAD cohort ( Table S7 ). All patients, and patients with Stage IA LUAD, in the high-risk group had worse OS than patients in the low-risk group ( Figures 3G, H ); there was no significant difference in OS between patients with Stage IA LUAD and high-risk scores and patients with Stage II LUAD ( Figure 3I ). Time-dependent ROC curve analysis at 1, 3, and 5 years showed that the 15-gene risk signature had better predictive performance than other clinical traits in the TCGA-LUAD cohort (1 year: AUC=0.811; 3 years: AUC=0.779; 5 years: AUC=0.780; Figures 3J–L ). DCA showed that the 15-gene risk signature had more clinical benefits than other clinical traits ( Figures 3M–O ).

A heatmap showed that the TNM stage, immune scores, NMF subgroup, gender, T stage, and N stage were significantly associated with risk scores in TCGA datasets ( Figure 4A , P<0.05). Most patients with high-risk scores were men ( Figure 4A , P<0.001) and had higher-grade tumors ( Figure 4B ). Among the 15 pGGO-related DEGs in the gene risk signature, DKK1, NPAS1, AL357143.1, KCNV1, AC068228.1, AC239859.6, and FGF5 were highly expressed in the high-risk group, and AC087763.1, PCP2, FAIM2, AL357143.1, AC022148.1, AC021678.2, LSP1P2, and LINC00563 were highly expressed in the low-risk group ( Figure 4C ). Risk curves based on a per-sample risk score also validated the predictive power of the 15-gene risk signature ( Figures 4D, E ). Univariate and multivariate COX regression analyses showed that the 15-gene risk signature is an independent prognostic factor in patients with Stage I–II LUAD ( Figures 4F, G : univariate Cox regression analyses: P<0.001, HR = 1.041; 95% CI: 1.029–1.054; multivariate Cox regression analyses: P<0.001, HR = 1.038; 95% CI: 1.024–1.052).

A nomogram incorporating clinical factors and the 15-gene risk signature was constructed for visualization and convenient clinical application ( Figure 5A ). Calibration curves validated the ability of the 15-gene risk signature to predict OS in the TCGA-LUAD cohort and LUAD dataset obtained from the GEO ( Figures 5B, C ). The ROC curve and DCA analysis validated that the nomogram had better predictive performance than the 15-gene risk signature alone at 1, 3, and 5 years in patients with Stage I LUAD from the TCGA-LUAD cohort ( Figures 5D–I ) and GEO ( Figures 5J–O ).

In the TCGA cohort, the GSEA suggested that the cell cycle, DNA replication, Parkinson’s disease, pyrimidine metabolism, and ribosome pathways were activated in patients with high-risk scores; meanwhile, allograft rejection, asthma, the intestinal immune network for IgA production, primary immunodeficiency, and systemic lupus erythematosus pathways were activated in patients with low-risk scores ( Figures 6A, B ). High-risk patients had higher TMB than low-risk patients ( Figure 6C and Table S8 ). Patients were stratified into a high-TMB group and a low-TMB group based on the median TMB value. There was no significant difference in OS between patients in the high-TMB group and low-TMB group; however, patients with a high-risk score and a high TMB had the worse OS ( Figure 6D ). Patients with a low-risk score and high TMB had the best OS ( Figure 6E ). Patients with high-risk scores always had low levels of immune infiltration ( Figure 6F ). Patients with a low-risk score and tumors that were CTLA4 positive and PD1 negative had a high IPS, suggesting that these patients may benefit from immunotherapy ( Figure 6G , Table S9 ).