Genomic features associated with pathologic ln metastasis in the PKPH Cohort

To identify genomic features associated with pathologic LN metastasis in LUAD primary tumors, we conducted the broad panel NGS analysis in the PKPH cohort. The cohort included a total of 257 patients who met the inclusion criteria (Table 1), with 54.8% (n = 141) were female and 71.6% (n = 184) were never smokers. Diagnostic CT scans revealed a solid morphologic appearance in the primary tumors of 53% (108/204) of pN-negative cases and 90% (48/53) of pN-positive cases, showing a significant difference (Fisher’s exact test, p < 0.001). The total primary tumor size, as measured on CT scans, was significantly larger in pN positive tumors (median [IQR], 2.9 [2.3–3.6] cm) compared to pN negative tumors (2.3 [1.6–3.1] cm, Mann–Whitney U, p < 0.001). The majority of patients (n = 193 [75%]) had clinical stage I disease, while 22% (n = 10) had clinical stage II disease. Pathologic review revealed that 79% (n = 204) of resected tumors were pN negative, while 21% (n = 53) were pN positive. The acinar/papillary predominant histologic subtype was significantly lower in pN positive tumors (59%) compared to pN negative tumors (23%, Fisher’s exact test, p < 0.001). Lymphovascular invasion and visceral pleural invasion were significantly higher in pN positive tumors (38/53 vs. 52/204, 30/53 vs. 35/204 for pN positive vs. pN negative, respectively, Fisher’s exact test, p < 0.001). The clinicopathologic features of the patients are summarized in Table 1.

Table 1 Clinicopathologic characteristics of the PKPH NGS cohort

In the PKPH NGS cohort, the most commonly detected drivers were EGFR (66%), TP53 (41%), and KRAS (12%). We identified 13 driver genes with a FDR of 0.1. Notably, TP53 mutations were significantly more prevalent in pN positive primary tumors compared to pN primary negative tumors (55% vs. 37%, Fisher’s exact test, q < 0.001, Fig. 1B). Furthermore, PIK3CG alterations were significantly more common in pN positive primary tumors compared to pN negative ones (7.5% vs. 0.5%, Fisher’s exact test, p < 0.05, Supplementary Fig. 1). Conversely, no statistically significant differences were observed in the prevalence of mutations between pN positive and negative primary tumors for other genes analyzed, including EGFR, KRAS, RBM10, RB1, STK11, SETD2, SF3B1, PTEN, HGF, BRAF, ARID1A, and KEAP1.

Fig. 1

Genomic Features Associated with Pathologic LN Metastasis. A OncoPrint displayed the gene alterations across the entire cohort, stratified by pathologic lymph node (LN) status. B Mutation frequencies in key driver genes were compared based on pathologic LN status. C-E Box plots showed the distribution of tumor mutation counts, TMB, MATH scores according to pathologic LN status. LN lymph node, TMB tumor mutational burden (mutations per megabase), MATH mutant-allele tumor heterogeneity

We performed subsequent analysis to examine genomic differences in the PKPH NGS cohort of LN metastasis, focusing on the tumor mutation count, tumor mutational burden (TMB), and mutant-allele tumor heterogeneity (MATH) score (Fig. 1C-E). The results demonstrated a significant disparity in the tumor mutation counts between pN positive and pN negative tumors (median [IQR], 6.0 [3.0-9.0] vs. 4.0 [2.0-6.0], p = 0.037, Fig. 1C). The TMB was notably higher in pN positive tumors (median [IQR], 4.74 [2.63–7.38] mutations/Mb) compared to pN negative tumors (median [IQR], 3.16 [1.75–4.92] mutations/Mb, p = 0.017, Fig. 1D). The MATH score exhibited a significantly higher value in pN positive tumors (median [IQR], 53.49 [19.03–72.51]) compared to pN negative tumors (median [IQR], 45.67 [8.81–54.45], p = 0.02, Fig. 1E), suggesting a relatively elevated level of intratumor heterogeneity in pN positive tumors.

The present study investigated the co-occurrence and mutual exclusivity patterns of EGFR and several key genes in a large panel of NGS data. Our analysis revealed that EGFR-RBM10, KRAS-STK11, KEAP1-STK11, KEAP1-SETD2, KEAP1-KRAS, HGF-STK11, HGF-KRAS, BRAF-SETD2, ARID1A-TP53, ARID1A-PTEN, ARID1A-KRAS, ARID1A-HGF, ARID1A-BRAF, and BRAF-SETD2 showed significant co-occurrence, whereas EGFR-KEAP1, EGFR-STK11, EGFR-KRAS, and BRAF-EGFR exhibited mutual exclusivity across the entire cohort (p < 0.05, Supplementary Fig. 2A, Supplementary Table 2.1). In the subgroup of pN positive tumors, we observed co-occurrence for KRAS-STK11, KEAP1-STK11, and KEAP1-KRAS, whereas EGFR-STK11, BRAF-EGFR, and EGFR-KRAS were mutually exclusive (p < 0.05, Supplementary Fig. 2B, Supplementary Table 2.2). Furthermore, we identified EGFR-RBM10, KRAS-STK11, KEAP1-STK11, KEAP1-KRAS, ARID1A-TP53, ARID1A-KRAS, ARID1A-HGF, ARID1A-BRAF, and BRAF-SETD2 as significant mutual co-occurrence drivers in pN negative tumors, while EGFR-STK11, EGFR-KRAS, and BRAF-EGFR exhibited mutual exclusivity in pN negative tumors (p < 0.05, Supplementary Fig. 2C, Supplementary Table 2.3).

Our analysis, which integrated alteration frequency, tumor mutation count, TMB, and MATH score, reinforced the association between genomic alterations and LN metastasis in LUAD.

Preoperative clinicopathologic and genomic features associated with pathologic LN Metastasis

To identify clinicopathologic and genomic features associated with the presence of pathologic LN metastasis, we compared pN positive primary tumors with pN negative primary tumors by logistic regression analyses. Univariable logistic regression revealed significant associations between several clinicopathologic and genomic factors and LN metastasis in the NGS cohort (Supplementary Table 1, Supplementary Fig. 1A). Specifically, we observed significant differences (p < 0.001) in primary tumor appearance and size on CT, lymphovascular invasion, and visceral pleural invasion. Additionally, mutations in TP53 (p = 0.028) and PIK3CG (p = 0.013) were notably and significantly correlated with LN metastasis.

The further results of the multivariable regression analysis performed by above features, showed that the tumor morphologic appearance on CT (odds ratio [OR], 0.276, 95% confidence interval [CI], 0.085–0.89, p = 0.032), tumor size on CT (OR, 1.036, 95% CI, 1.001–1.071, p = 0.042), LVI (OR, 3.89, 95% CI, 1.63–9.29, p = 0.02), and the presence of PIK3CG mutation compared to wild-type (p = 0.049) were associated with pathologic LN metastasis (Supplementary Table 1, Supplementary Fig. 1B).

This study reinforced our understanding of genomic variables, such as PIK3CG, linked to LN metastasis in LUAD. Further research is necessary to fully elucidate their role in pathogenesis.

Oncogenic pathways and therapeutic actionabilities associated with pathologic LN Metastasis

To further examine the association between oncogenic pathways and LN metastasis in the analysis of genomic variables above, we conducted a comparison of oncogenic signaling pathway alteration frequencies between pN positive and negative samples in the PKPH NGS cohort. Results of the Fisher’s exact test revealed a significant increasing trend in the frequency of p53 pathway alterations from pN negative to pN positive samples (p = 0.04, Fig. 2B). Nevertheless, no significant difference was found probably due to limited sample size. The co-occurrence and mutual exclusivity patterns analysis of oncogenic pathways is presented in Supplementary Fig. 2.

Fig. 2

Oncogenic Pathway and Therapeutic Actionability Associated with Pathologic LN Metastasis. A Somatic mutation comparisons were made between PKPH cohort and the MSK cohorts. Each gene box represented the mutational frequencies of genes in the entire cohort, pN0 (negative) and pN + (positive) groups, and the MSK cohort as a whole, including pN0 and pN + groups, as shown in the graph. The color scale bar indicated mutation frequency from 0 to 100%. Genes were grouped by signaling pathways related to genome maintenance mechanisms. Interactions between genes were indicated by arrows. B Comparison of alteration frequency of oncogenic pathways among pathologic LN status. Pathways with significant differences from Fisher’s exact test were labeled with asterisks (FDR q <0.1), while TP53 significant difference is highlighted (p < 0.05). C Frequency of samples with actionable alteration by pathologic LN status. D Level of evidence versus pathologic LN status. Samples were classified by the alteration that carried the highest level of evidence. E Frequency of the number of actionable alterations versus pathologic LN status. TSG tumor suppressor gene, FDR false discovery rate, VUS variant of uncertain significance

A further comparison of somatic mutations in genes related to genomic integrity between our PKPH NGS cohort and the Western MSK cohort (n = 429) was conducted to identify any significant differences or similarities in mutation profiles across different populations. KRAS (12%), NF1 (1%), and BRAF (4%) of the receptor tyrosine kinase (RTK)/RAS pathway, STK11 (3%) of the phosphatidylinositol 3-kinase (PI3K) pathway, and KEAP1 (2%) of the NRF2 pathway had lower mutation frequencies, while EGFR (65%) of the RTK/RAS pathway, PIK3A (5%) of the PI3K pathway, and TP53(40%) of the p53 pathway had much higher mutation frequencies in our cohort than in Western populations (Fig. 2A). PIK3A (8%) of the PI3K pathway was detected with higher mutation frequencies in the pN positive group compared with those in the negative group (4%) in our PKPH cohort and MSK cohort (6% vs 3%). The same observations were also identified in BRAF of the RTK/RAS pathway in the PKPH cohort (8% vs 2%) and MSK cohort (7% vs 4%). These findings suggested potential differences in mutation patterns and highlighted the importance of RTK/RAS and PI3K pathway-related genetic mechanisms in primary tumors, particularly in relation to LUAD LN metastasis, between our PKPH cohort and the MSK cohort.

A total of 260 actionable alterations were identified using the OncoKB database, spanning 18 genes and consisting of 184 (71%) level 1, 7 (3%) level 3A, 23 (9%) level 3B, and 46 (17%) level 4 alterations (Fig. 2C-E). The RTK/RAS pathway had the highest number of actionable alterations (84% [218/260]), of which 84% (184/218) had level 1 evidence. At the sample level, actionable alterations were identified in 84% (216/257) of the 257 samples, of which 83% (180/216) had level 1 evidence. Notably, pN positive tumors showed a higher frequency tendency of level 1 actionable targets (73% vs. 60% for pN + vs. pN0, p = 0.094). The mean number of actionable alterations per sample among subgroups did not differ significantly (1.038 vs. 1.005, all p > 0.05).

Our analyses of oncogenic pathways and therapeutic targets corroborated existing data, underscoring molecular mechanisms linked to LN metastasis in LUAD. These findings provided foundational insights into diverse genomic mutation patterns.

Identification of cellular metaclusters related to ln metastasis in TIME

To explore the spatial cellular composition in primary tumors associated with LN metastasis, we utilized mIHC to analyze primary tumor samples from 92 patients after propensity score matching in PKPH NGS cohort. We optimized a 6-antibody panel to identify cancer cells, stromal cells, as well as innate and adaptive immune lineages with diverse functional characteristics (Fig. 3A).

Fig. 3

mIHC Defined the Spatial Landscape of LN Metastasis. A A schematic depicted the process of mIHC image acquisition from 92 patients with LUAD, including tissue preparation, antibody staining, image acquisition and assembly, single-cell segmentation, cell dimensionality reduction, clustering, and LN stage analysis. B Heatmaps showed the z-scored mean marker expression of 9 metaclusters, with their frequency distribution patterns across different patients displayed in the box plot on the right. Pixel intensities were transformed using an asinh transformation with a cofactor of 1. C A UMAP plot represented the 9 annotated metaclusters. D, E The prevalence of 9 cell types across 92 patients with LUAD was shown as a proportion of immune cells (D) and total cells (E). F A waterfall plot depicted the distribution of 9 cell types across LN status. G Representative mIHC images from LN metastasis primary tumor samples were analyzed using the 6-marker panel. H The prevalence of 9 metaclusters was compared across LN stages. The comparison of macrophage and Treg cells between pN positive and negative tumors was shown: *p < 0.05, **p < 0.001. Data were presented as means ± SEMs. Statistical analysis was conducted using the Mann–Whitney U test

We utilized a dimensionality reduction approach to classify 9 distinct cellular metaclusters (Fig. 3B, C), and the statistical significance of differences between these metaclusters was evaluated using the Mann–Whitney U test. pN positive tumors exhibited a higher percentage of immune infiltrate (36.5%) compared to pN negative tumors (31.1%, Fig. 3D, E). Notably, the disparity was primarily driven by alterations in the macrophage (7.1% vs. 4.8%, p < 0.05) and Treg (0.85% vs. 0.29%, p < 0.001) populations (Fig. 3G), while no significant variations in the overall frequency of total cells were observed between different LN metastasis statuses. No significant differences were observed in the average frequency of CD4+T, CD8+T, and epithelial cellular metaclusters in total cells between pN positive and negative tumors. Our analyses involved the classification of 9 distinct cellular metaclusters in primary tumor samples of LUAD, emphasizing the notable contributions of macrophages and Tregs proportions in the TIME [26].

Fig. 4

Variability in mIHC Distributions Across Clinical Variables and Cellular Interaction Profiles Across the LN Stage in LUAD. A An interaction partner detection model was used to identify neighboring cells within 4 µm of central cells. B Segmented images showed a decrease in interactions between PD-L1 positive epithelial cells and macrophages in pN positive LUAD. C A heatmap depicted significant pairwise cellular interactions (red) or avoidance (blue) across LN stage subgroups. The boxes highlighted associations referenced in the text. D Two distinct subtypes were identified within the tumor immune microenvironment (TIME) of LN metastasis primary tumors through unsupervised analysis with Non-negative Matrix Factorization: an epithelial type and a macrophage and immune-regulatory cell-enriched type. E An association was found between the frequency of LN metastasis and the distribution of the two TIME subtypes (p < 0.001). F The prevalence of cellular metaclusters was compared between TIME subtypes: *p < 0.05, **p < 0.001. Data were presented as means ± SEMs. Statistical analysis was conducted using the Mann–Whitney U test

Spatial cellular interactions related to ln metastasis and cellular metaclusters analyses reveal two time subtypes

In order to understand the cellular architecture and spatial organization beyond the metaclusters frequency of LN metastasis in LUAD, we analyzed the communication patterns between individual metaclusters by quantifying the spatial relationships (Fig. 4A). The distribution of the numbers of cellular connections and pairs of cell distance (μm) was presented in Supplementary Fig. 3. Furthermore, we determined the probability of interaction or avoidance behaviors between pairs of cells across tumor architectures in both pN negative and pN positive tumors (Fig. 4B). By comparing the cellular interactions within metaclusters between pN negative and pN positive tumors, macrophages were indicated at the center exhibited significant differences in spatial cell interactions with a broader range of other subclusters, including PD-L1 positive epithelial cells (Fig. 4C, box a), macrophages and T, macrophages and CD8+T (Fig. 4C, box b), CD4+T and CD8+T (Fig. 4C, box c). The interactions between macrophages and tumor epithelial cells in the primary tumor lesion of pN positive tumors were higher, which aligns with previous research highlighting the role of macrophages in the TIME [27, 28]. These analyses painted an overall picture of how cellular metaclusters interactions shift among LN metastasis in LUAD and exemplified how spatial relationships, rather than cell frequency alone, are important to understand patterns of spatial TIME biology.

To investigate the patterns involved in cellular metaclusters and genetic characteristics, we identified two distinct subtypes within the TIME of primary tumors with LN metastasis using non-negative matrix factorization algorithm (Fig. 4D). These subtypes consisted of an epithelial subtype (subtype 1) and a macrophage and Treg enriched subtype (subtype 2), as shown in Figure F. We found a significant association between the frequency of LN metastasis and the distribution of these two subtypes (Fisher’s exact test, p < 0.001, Fig. 4E), with subtype 2 demonstrating a higher propensity for advanced LN stages. This novel approach [29], which utilizes nine metaclusters for classification, provides a new perspective on the heterogeneity of the LUAD TIME in relation to LN metastasis.

Using the Mann–Whitney U test to analyze differences in metaclusters between subtypes, we found significant differences in the proportions of CD4+T cells between subtype 1 (mean = 0.0649) and subtype 2 (mean = 0.0539) with p = 0.0166, indicating distinct roles for CD4+T cells across these subtypes. In contrast, the proportions of CD8+T cells did not differ significantly between the subtypes (subtype 1 mean = 0.0600, subtype 2 mean = 0.0651, p = 0.3235, Fig. 4F).

The observed difference in epithelial cells was particularly significant, with subgroup 1 showing an average value of 0.4171, compared to 0.2195 in subgroup 2 (p < 0.001), highlighting the critical role of epithelial cells in the LN metastatic TIME. Similarly, macrophages showed significant distribution differences (subtype 1 mean = 0.0662, subtype 2 mean = 0.0919, p < 0.001, Fig. 4F), suggesting divergent functions across two TIME subtypes.

Further insights were gained from the analysis of gene mutation frequencies across the two distinctive TIME subtypes (Supplementary Fig. 4C). The ATM and PIK3CG genes exhibited the lower mutation frequency in subtype 1 and higher in subtype 2 (ATM p = 0.013, PIK3CG p = 0.016), potentially linked to the cell’s DNA damage repair capacity and activation of immune inflammatory cells. Both BRD4 and KMT2B genes showed higher mutation frequencies in subtype 2 (BRD4 p = 0.041, KMT2B p = 0.041), indicating that these epigenetic regulators might play different roles across LUAD TIME subtypes. The CTNNB1 gene showed a higher mutation frequency in subtype 1 (p = 0.046), associated with potential changes in the Wnt signaling pathway, impacting cell fate and proliferation.

In the subtype 1, we observed co-occurrence for EGFR-TP53, KRAS-TP53 and EGFR-CTNNB1 (p < 0.05, Supplementary Fig. 4B). Furthermore, we identified KRAS-TP53 and EGFR-PIK3CG as significant mutual co-occurrence drivers in subtype 2 (p < 0.05, Supplementary Fig. 4B).

Our analyses offered early insights into the dynamic interactions of cellular metaclusters in LUAD primary tumors, focusing on LN metastasis. The findings suggested that spatial relationships, rather than cell frequency alone, were essential for understanding spatial TIME biology [30, 31].

Analysis of distinct CN alterations between ln negative and positive LUAD patients

In the preceding section, we conducted an analysis of the variations in cellular metacluster frequencies observed between pN negative and positive primary tumors in patients with LUAD. To determine specific cellular interacting critical for the LN metastasis, we conducted CN and cell interaction assessments. We annotated the CNs based on the major cellular cluster within the TIME (Fig. 5A, B, C). First, a cellular neighbourhood was designated as the ten closest neighbours surrounding the central cell. The cellular neighbourhoods were defined to group cells, encompassing the central cell along with its immediate neighbours within a 4 µm radius, as well as the secondary neighbours of these primary neighbours. Such definition is more suitable for the CN with > 10 neighbours and able to remove artificial neighbours when separate cells were in distant location in Fig. 5A and 5B. The distribution of the CN is presented in Supplementary Fig. 5A. Recent CN analysis framework established by Schürch et al. [32] and Sheng et al. [24] were adapted.

Fig. 5

Correlations of LUAD TIME Cellular Function Units with Patient LN Status. A The analysis schedule for cellular neighbourhoods (CNs) was outlined. CNs were defined by the center cell, primary neighbours, and secondary neighbours. The peak number of cells within a CN was approximately 10. B Ten major types of CNs were first annotated based on the major cell type within each CN and were clustered according to the scaled frequency of each cell type within the CNs. C An overview of clinical, genetic alterations, and CN signatures was provided for 92 patients. D The inner and outer layers of the donut chart respectively represented the differences in CN distribution between pN negative and positive cases. E Statistical analysis showed a correlation between immune-suppressed-enriched CNs (CN 8 and 9) and the Treg component, highlighting differences between pN negative and positive cases: CN 8 (p > 0.05) and CN 9 (p = 0.002). F Statistical analysis showed a correlation between macrophage and T cell-enriched CNs (CN 4, 7, and 10) and the macrophage component, with differences observed between pN negative and positive cases: CN 4 and 7 (p > 0.05), CN 10 (p = 0.04). G In the mIHC image, a Voronoi plot with distinct colors represented the distribution of CN 9 (gray) and CN 10 (purple), both of which were more abundant in pN positive cases. CN cellular neighbourhood

To determine various functional cellular neighbourhood units, these neighbourhoods were annotated based on the primary cellular clusters (Fig. 5B, D). Our analysis revealed distinct CN functional units, including epithelial-enriched CNs (Epithelial-CN, 2, 5, 6), macrophage and T cell-enriched CNs (MF and T-CN, 4, 7, 10), immune-suppressed-enriched CNs (Immune-suppressed-CN, 8, 9), and CNs enriched in CD8+ and CD4+ T cells (T-CN, 1). Notably, the abundance of Immune-suppressed-CNs exhibited coherent links with the Treg metacluster proportion (Fig. 5E), whereas macrophage and T cell-enriched CNs showed links with the macrophages metacluster proportion (Fig. 5F).

We used the Mann–Whitney U test to compare CNs between pN-positive and pN-negative groups and observed a significant increase in the number of CN9 exclusively in the pN positive LUAD patients primary tumor (Fig. 5E). CN10 also increased in pN positive group of LUAD patients primary tumor (p = 0.04, Fig. 5F). However, However, the CN4 and CN7 showed no LN metastasis deference (p > 0.05, Fig. 5F). In contrast, no such trend was observed in the epithelial-enriched CNs (Epithelial-CN, 2, 5, 6) and T-CN (p > 0.05, Supplementary Fig. 5B, C). These results suggested that LUAD topological TIME units could be novel biomarkers for LN metastasis. In addition, macrophage and Tregs organized CNs that play a critical role in the TIME of LN positive primary tumor.

We applied a color-coding scheme to each mIHC image based on the CN composition, and utilized a Voronoi plot to generate a topology map for each mIHC image (Fig. 5G), employing the methodology from our previous research. The resulting topology map, encoded with CN information, accurately represented the mIHC image.

Subsequently, our analysis uncovered significant associations between the abundance of CNs and gene mutations (Supplementary Fig. 5D). Specifically, PIK3CG mutations consistently showed a higher enrichment of macrophage and T cell-enriched CN (CN10) and T-CN (CN1), as well as a lower enrichment of epithelial-CN (CN5, CN6). Conversely, KRAS mutations were associated with a lower enrichment of MF and T-CN (CN10) and T-CN (CN1), while exhibiting a higher enrichment of epithelial-CN (CN5).

Our spatial proteomic analysis revealed novel findings, including macrophage functional heterogeneity and T cell enrichment associated with PIK3CG mutation in LUAD primary tumors. By integrating genetic mutation profiles with proteomic data, we built a comprehensive spatial genoproteomics landscape of LN metastatic LUAD primary tumors.

Predicting ln metastasis using machine learning model

To investigate the clinical relevance of genoproteomic features in predicting LN metastasis in LUAD, we assessed the combined utility of NGS and mIHC profiles in predicting LN stage with high accuracy. Our findings indicated that the combined NGS and mIHC profiles, collectively termed ImGene, are predictive of LN stage. This prompted us to investigate whether this integrated data could be utilized to predict LN metastasis, independent of surgical options, through a machine learning approach (Fig. 6A and B). Specifically, the genoproteomic features highlighted in Fig. 6B were selected for their predictive strength within the machine learning model. Consequently, the model identified features that were not directly highlighted in prior NGS or mIHC analyses but still contribute to predicting pN status through their interactions, as validated in both the PKPH cohort and the independent PKTOI cohort.

Fig. 6

Machine Learning Model of NGS and mIHC Profiles Predicted LN Stage. A A schematic illustrated the machine learning-based strategy, which involved feature identification and model optimization. B The formula for calculating the coefficients of our ImGene SVM model was presented. C A box plot showed the coefficients of genetic and mIHC features used in the model. D ROC curves were presented for our LN stage prediction model applied to the PKPH and independent PKTOI cohorts. The AUC was 0.86 for the training set and 0.81 for the validation set. PKPH Peking University People’s Hospital, PKTOI Peking University People’s Hospital Thoracic Oncology Institution, ROC receiver operating characteristic

The ImGene model, developed by integrating both NGS and mIHC profiles, demonstrated superior performance over monomodal models based solely on genetic or mIHC data. It achieved an accuracy of 0.82, an F1 score of 0.83, and an area under the curve (AUC) of 0.86 (Supplementary Fig. 6A). This model demonstrated greater accuracy in estimating the probability of N2 stage in lung cancer than both the Peking University (PKU) model, includes variables factors such as age, tumor size, tumor location, and histological subtypes [33] (AUC = 0.77), and the Fudan University model (AUC = 0.67, Supplementary Fig. 5C) for predicting LN stage [34].

In an independent validation set, the PKTOI cohort, the models demonstrated the following performance: the ImGene model achieved an AUC of 0.80, with an accuracy of 0.81 and an F1 score of 0.77 (Fig. 6D, Supplementary Fig. 6A). In contrast, the GeneFeatures monomodal model, which utilizes only NGS profiles, had an AUC of 0.46, while the ImFeatures monomodal model, based solely on mIHC profiles, reached an AUC of 0.74 (Supplementary Fig. 6C). These results underscored the efficacy of our integrated algorithm in analyzing the genoproteomic profiles of primary LUAD tumors to accurately predict LN stage probability.