Of 1,293 patients who underwent surgery for lung cancer at St. Marianna University Hospital between 2000 and 2020, 186 (14%), 12 (1%), and 49 (4%) had tumors that were histologically confirmed to represent the LPA, MPA, and SPA subtypes, respectively ( Figure 1A ). SPA, a high-risk subtype, was contrasted to LPA, a low-risk subtype. The pathological specimens were independently reviewed by two pathologists (H. N. and S. N.) to confirm that they fulfilled the 2015 World Health Organization classification criteria for lung tumors (histological criteria) (12). FFPE tumor tissue blocks from 12 surgical specimens histologically confirmed as lung LPA, MPA, and SPA were obtained without patient identifiers from the St. Marianna University School of Medicine Hospital. Informed consent was obtained from all participating subjects. The protocol was approved by the Institutional Review Board of St. Marianna University School of Medicine (approval no. 1461) and the study adhered to the Helsinki Declaration.

For tissue microdissection, 10-μm-thick sections from the FFPE tumor blocks were cut and placed on DIRECTOR slides (OncoPlex Diagnostics Inc., Rockville, MD, USA). The sections were deparaffinized and stained with hematoxylin using standard protocols before dissection. Microdissection was performed using a Leica LMD7 microdissection microscope (Leica, Wetzlar, Germany) ( Figure 1B ). A total area of 4 mm² with approximately 15,000 tumor cells was directly transferred from the FFPE sections via laser dissection into the cap of a 200 μL low-binding tube. Proteins were extracted and digested with trypsin using the Liquid Tissue™ MS Protein Prep kit (OncoPlex Diagnostics Inc.) according to the manufacturer’s instructions (13). Briefly, dried microdissected pellets were suspended in 20 μL of Liquid Tissue buffer and heated at 95°C for 90 min and then cooled on ice for 3 min before the addition of 0.1 μg of trypsin. The tubes were then incubated at 37°C overnight. Dithiothreitol was added to a final concentration of 10 mM and the samples were heated for 5 min at 95°C. The digested samples were dried, resuspended in 50 μL of a 2% acetonitrile aqueous solution containing 0.1% trifluoroacetic acid, and stored at −20°C until analysis.

A label-free quantitation approach using spectral counting by LC-MS/MS was adopted for global proteomic analysis. The digested samples (5 μL for a single run) were analyzed in triplicate by LC-MS/MS using a reverse-phase LC system interfaced with a Q Exactive Orbitrap mass spectrometer (Thermo Fisher Scientific, Bremen, Germany) via a nano-electrospray ionization device (AMR Inc., Tokyo, Japan). The mass spectrometer was operated in data-dependent acquisition mode. Full-scan MS spectra were acquired in the range of m/z 350–1600 at a resolution of 70,000. The top ten most intense peaks from the survey scan were selected for fragmentation with higher-energy collisional dissociation with a normalized collision energy of 27% and isolation window of m/z 1.6. The dynamic exclusion time for precursor ions selected for MS/MS fragmentation was 15 s, and the automatic gain control target values for MS and MS/MS were 1 × 10⁶ and 1 × 10⁵, respectively.

The LC system consisted of an Ultimate3000 HPLC System (Thermo Fisher Scientific), a trap cartridge (0.3 mm × 5.0 mm, CERI, Tokyo, Japan), and a capillary separation column (Zaplous column alpha-PepC18, 3 μm, 12 nm, 0.1 mm × 150 mm, AMR Inc.) fitted with an emitter tip (FortisTip, OmniSeparo-TJ, Hyogo, Japan). An auto-sampler (HTC-PAL, CTC Analytics, Zwingen, Switzerland) was used to load the samples into the trap, which was then washed with solvent A (2% acetonitrile aqueous solution containing 0.1% formic acid) to concentrate and desalt the peptides in the trap. Subsequently, the trap was connected in series to the separation column, and the peptides were eluted from the whole column with 0.1% formic acid aqueous solution and acetonitrile using a linear 5%–40% acetonitrile concentration gradient over 90 min at a flow-rate of 500 nL min⁻¹.

The raw data were processed using PatternLab for Proteomics software v4.0. Peptide sequence matching (14) was done using the Comet algorithm against the UniProt Homo sapiens database. A target-reverse strategy was employed for increased confidence in protein identification. This search considered tryptic peptide candidates as well as the formylation of lysine and oxidation of methionine as variable modifications. The Comet search engine considered a precursor mass tolerance of 40 ppm and a fragment bin tolerance of 0.02 Da. The validity of the peptide spectrum matches was assessed using the Search Engine Processor (SEPro) module of PatternLab. The acceptable FDR for spectra, peptide, and protein were 3%, 2%, and 1%, respectively. The expressions of the identified proteins were assessed using spectral count-based protein quantification. The spectral count represents the number of MS/MS spectra assigned to each protein.

The similarity in protein expression patterns for all protein pairs was calculated according to their pairwise Pearson’s correlation coefficient [i.e., the similarity between proteins i and j was defined as (1-r _i,j)/2, where r _i,j is the Pearson’s correlation coefficient of the protein expression pattern between the two proteins]. We performed a network topology analysis for various soft-thresholding powers ranging from 1 to 2 to choose an optimal value of balance between independence and mean connectivity. A topological overlap matrix (TOM) that considers topological similarities between a pair of proteins in the network was then generated from the resultant scale-free co-expression network. We generated a tree by hierarchical clustering using dissimilarity according to TOM (1−TOM), and protein modules were determined using dynamic tree-cutting to trim the branches (11).

The modules were summarized by the first principal component, which is referred to as eigen proteins in the text, as they express the highest connectivity in the module. Module membership, defined as the correlation between the protein expression profile and the module eigen-protein, was measured with values ranging from 0 to 1, with “0” representing a gene that is not part of the module and “1” representing high connectivity with the module. Subsequently, the module-trait association was determined using the correlation between the module eigen-protein and the three subtypes: MPA, SPA, and LPA. A protein module was summarized by the top hub protein (referred to as “eigen-protein”) with the highest connectivity in the module. The WGCNA analysis was performed using the WGCNA R-package (11) implemented in RStudio.

We used the STRING database (version 11.5) (https://string-db.org/) to construct a PPI network for a protein module (15). STRING networks were calculated under the criteria for linkage with experiments, databases, text mining, and co-expression using the default settings (medium confidence score: 0.400; network depth: 0 interactions). Functional enrichment results were obtained for canonical pathways with a p-value <0.05. Proteins in a module were mapped in the PPI network from the STRING database to produce the results of the enrichment analysis regarding the biological process (GO) and Reactome pathways (HAS). Protein networks were subsequently exported to Cytoscape (version 3.9.1) (Institute for Systems Biology, Seattle, WA, USA: https://cytoscape.org/) (16) from the STRING database. The hub proteins in each module were identified according to their intramodular connectivity and their correlation with module eigen-proteins. The proteins inside the co-expression modules exhibit high connectivity and the proteins within the same module may play similar roles. The top 10 high-degree proteins were identified using the cytoHubba plugin (17). The top-ranked proteins in each module were considered hub proteins and designated “highly connected proteins.” Functional enrichment results were obtained for canonical pathways by considering a network bias-corrected p-value of <0.05 for statistical significance.

The multivariate correlation analysis (MVA) of semiquantitative key protein expressions was performed using the JMP software (SAS Institute, Cary, NC, USA), and which result was visualized using the Intervene Shiny App (https://intervene.shinyapps.io/intervene/) (18).

Upstream regulators, causal networks, and canonical pathways were predicted using IPA software (http://www.ingenuity.com) (10). Quantile-normalized protein expression data of the selected modules were used as input datasets. Both the upstream regulators and causal networks (p < 0.05) predicted from the WGCNA network modules were significantly associated with the three subtypes (LPA, MPA, and SPA), in which the activation and the inhibition of a predicted network were defined by z-values that were >2.0 and <−2.0, respectively. The upregulation was defined by z-values > 1.5 and < 2.0, whereas downregulation was defined by z-values > −2.0 and < −1.5.

At the St. Marianna University School of Medicine Hospital (Kawasaki, Kanagawa, Japan), 1,293 lung cancer patients were pathologically confirmed and underwent surgery between January 2008 to April 2022. These cases consisted of the following subtypes: LPA (14%), MPA (1%), SPA (4%), and others (81%) ( Figure 1A a). The SPA cases exhibited worse overall survival (OS) compared with the LPA cases (log-rank test: p = 3.06 × 10⁻⁵) ( Figure 1A b), whereas the MPA cases were not statistically significant because of the limited number of samples.

FFPE tissue specimens were obtained by surgical resection of early-stage (IA-IB) lung adenocarcinoma patients who were pathologically confirmed as having solid predominant (SPA, n = 5), micropapillary (MPA, n = 3), and lepidic (LPA, n = 4) subtypes ( Table 1 ). Presurgical treatment was not performed for any of the lung adenocarcinoma patients. MS-based proteomic analysis was carried out on approximately 15,000 cancerous cells collected from the FFPE specimens by laser microdissection (LMD) ( Figure 1B ). A total of 2,287 proteins were identified, of which 1,199 (52.4%) were commonly expressed in the cancerous cells. LPA, MPA, and SPA contained 398 (17.4%), 100 (4.4%), and 49 (2.1%) unique proteins, respectively ( Figure 1C ).

Overall, 2,112, 1,753, and 1,361 proteins were identified from LPA, MPA, and SPA, respectively. A gene ontology (GO) analysis was performed using the Protein Analysis Through Evolutionary Relationships (PANTHER, version 17.0) program (Paul D. Thomas, University of Southern California, Los Angeles, CA, USA) (19) and revealed similar results among the three subtypes ( Supplementary Figure S1 ). Volcano plots were generated from the protein expression data obtained by SimpliFi™ software (PROTIFI, Farmingdale, NY, USA; https://simplifi.protifi.com). Upregulated proteins significant to SPA included large ribosomal subunit protein uL29 (RPL35), human leukocyte antigen G (HLA-G), GTP-binding nuclear protein Ran (RAN), human leukocyte antigen B (HLA-B), and putative heat shock protein HSP90-beta-3 (HSP90AB3P). The expression of HLA-G protein was significantly associated with and highly upregulated in the SPA subtype ( Figure 2A ), which was confirmed by its ANOVA test with a high significance to the SPA subtype with p-ANOVA = 0.0026 ( Figure 2C ).

Weighted correlation network analysis, also known as weighted gene co-expression network analysis (WGCNA) incorporates traditional data exploration techniques, but its intuitive network language and analysis framework go beyond standard analysis techniques. It uses a network approach and is suitable for the integration of complementary genomic/proteomic datasets so that it can be interpreted as a data analysis technique for systems biology (11).

Following hierarchical clustering of the samples based on protein abundance ( Figure 3A ), a WGCNA analysis (11) was performed with a soft threshold power of 10, which was selected to approximate a scale-free topology, a minimum module size of 10, and a module detection sensitivity (deepSplit) of 4. Correlations between the resultant modules and the traits were obtained to identify protein modules that were significantly associated with the respective traits. Forty protein modules were identified by clustering all of the proteins and constructing weighted protein co-expression networks. The protein cluster dendrogram is presented in Figure 3B . A heatmap of eigen-protein expression for pairwise correlations between the modules in the connectivity measure (kME) of the module eigen-protein are presented in Figure 3C , together with lists of module IDs, module colors, and their eigen-proteins representing protein expressions in a module.

Six WGCNA modules were identified with high correlations (r > 0.5) and statistical significance (multiple testing correction using the Benjamini–Hochberg method: q < 0.05) with clinical traits ( Figure 4 ). The WM25 (cyan: r = 0.87, q = 0.0102) and WM33 (green-yellow: r = 0.81, q = 0.0286) modules were significantly correlated with the LPA subtype. The WM23 (dark-grey: r = 0.91, q = 0.014) and WM24 (dark-orange: r = 0.94, q = 0.0478) modules were significantly correlated with the MPA subtype, and WM1 (pink: r = 0.82, q = 0.0216) was well correlated with the SPA subtype.

Several other WGCNA modules were also correlated (r > 0.5) with one of the three traits. WM8 (dark-magenta: r = 0.76), WM19 (light-cyan: r = 0.61), and WM32 (violet: r = 0.61) correlated with the LPA subtype, WM11 (dark-red: r = 0.72) and WM18 (orange: r = 0.59) with the MPA subtype, and WM2 (turquoise: r = 0.57) with the SPA subtype. However, none of these correlations were significant (q > 0.05). Trait correlation analysis often tends to overlook important modules. Statistical over-representative analysis (ORA) may help evaluate potential WGCNA modules concerning overlap with uniquely upregulated proteins for each trait. The numbers of proteins in each subtype more than twice as highly expressed compared with other subtypes were 407 in LPA, 321 in MPA, and 482 in SPA, respectively. Overlap of the WGCNA modules with that of the protein groups was assessed using ORA ( Figure 4 ). In this study, we focused on the WGCNA modules, WM1 (pink) and WM2 (turquoise), significantly associated with the high-risk subtype, SPA ( Figure 4 ).

Using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database version 11.5 (https://string-db.org/) (15), human PPI networks were obtained for the WGCNA modules. The PPI networks for the WM1 and WM2 modules associated with the SPA subtype were reconstructed using Cytoscape (version 3.9.1) software (Institute for Systems Biology, Seattle, WA, USA: https://cytoscape.org/) (16) ( Supplementary Figure S2 ). Top hub proteins were calculated using the cytoHubba plugin with maximal clique centrality (MCC) (17). In these data-driven protein co-expression networks, eigen-proteins are indicated in red letters, hub proteins in red to orange fill colors, and some key proteins are denoted with red borders. Top STRING enrichment results for the WGCNA modules, WM1 and WM2, significant to the SPA subtype are shown in Supplementary Table S1 .

The functional enrichment obtained for the WM1 module included the following: (i) SRP-dependent co-translational protein targeting to membrane, translational initiation, nuclear-transcribed mRNA catabolic process, nonsense-mediated decay, and protein targeting to membrane as biological processes (GO); and (ii) L13a-mediated translational silencing of ceruloplasmin expression, GTP hydrolysis and joining of the 60S ribosomal subunit, eukaryotic translation termination, and NMD independent of the exon junction complex (EJC) ( Supplementary Table S1 ). The eigen-protein RPS4X (small ribosomal subunit protein eS4, also known as SCR10) is a component of the 40S subunit and is involved in L13a-mediated translational silencing of ceruloplasmin expression, eukaryotic translation termination, and NMD independent of the exon junction complex (EJC).

The pathways enriched for the WM2 (turquoise) module included the following: (i) antigen processing and presentation of peptide antigen via MHC class I, NIK/NF-κB signaling, and regulation of transcription from RNA polymerase II promoter in response to hypoxia as biological processes (GO); (ii) regulation of expression of SLITs and ROBOs, cross-presentation of soluble exogenous antigens (endosomes), and NIK to noncanonical NF-κB signaling ( Supplementary Table S1 ). The eigen-protein tropomyosin 4 (TPM4) is a member of the tropomyosin family of actin-binding proteins, which is involved in stabilizing cytoskeleton actin filaments.

Forty eigen-proteins together with fifteen key proteins, which were representatively expressed throughout all 40 modules, were subjected to multivariate correlation analysis (MVA), assuming three-group clustering ( Figure 2D ). The eighteen eigen-proteins of the WGCNA modules relatively well correlated (r > 0.4) with respective subtypes are denoted in colored filled circles (LPA: red; MPA: blue; SPA: green). Those eighteen eigen-protein expressions were clustered into Cluster A, Cluster B, and Cluster C, which well corresponded to the SPA, LPA, and MPA subtypes, respectively. ( Figure 2D ). Upregulated proteins characteristic of the respective subtypes were also well clustered consistently: HLA-G and HLA-B in SPA; C3, RCC1, and HSPG2 in LPA; MYOF and CRABP2 in MPA.

Zhou et al. (20) performed proteomic analyses between low-risk and high-risk subtypes of early-stage lung adenocarcinomas, where the low-risk LPA subtype group (n = 31) and the high-risk subtype group (n = 28) consisting of both MPA and SPA and reported several proteins differentially expressed more to the high-risk group. Their results were not well linked to our proteomic observations. Some proteins differentially expressed in their high-risk subtype group, such as DNA replicating licensing factors MCMs, were oppositely upregulated in the low-risk subtype LPA subtype in our study, except that the expression of upregulated P4HA2 (Prolyl-4 hydroxylase subunit alpha-2) was observed in the high-risk subtype MPA (p = 0.03). There are intrinsic differences between their approach and ours. They directly used sections of tissues collected immediately after resection for proteomic analysis while we collected only cancerous cells directly from the FFPE tissue sections via laser dissection, and they grouped both MPA and SPA as one high-risk group while we did not group MPA and SPA into one.

The canonical pathways and their pathway categories enriched were obtained by the IPA software (http://www.ingenuity.com) (10) for the WM1 and WM2 modules significant to the SPA subtype ( Figure 5 ). The WM1 network module was associated predominantly with activated pathways of EIF2 signaling, which are involved in the following pathway categories: Cellular Growth, Proliferation and Development, Cellular Stress and Injury, and Intercellular and Second Messenger Signaling ( Figure 5A ). The WM2 network module was associated with various pathways. The top pathway category was Cellular Immune Response, which included the FAT10 signaling pathway, caveolar-mediated endocytosis signaling, antigen presentation pathway, upregulated leukocyte extravasation signaling, and activated Fcγ receptor-mediated phagocytosis in macrophages and monocytes. The pathway category Cellular Stress and Injury included the FAT10 signaling pathway, inhibition of ARE-mediated mRNA degradation pathway, regulation of eIF4 and p70S6K Signaling, activated pathways of EIF2 signaling, activated NRF2-mediated oxidative stress response, senescence, leukocyte extravasation signaling, hepatic fibrosis signaling, HIF1A signaling, and necroptosis signaling ( Figure 5B ).

Upstream regulator and causal network analysis were performed for the WGCNA modules using IPA software (http://www.ingenuity.com) (10). Table 2 summarizes the top upstream regulators, causal networks, canonical pathways, and diseases or functions predicted for the four WGCNA modules.

Highly activated upstream regulators in the WM1 protein networks included MLXIPL (z = 4.80, p-value of overlap = 4.16 ×10⁻³¹), MYCN (z = 4.12, p-value = 6.30 ×10⁻¹⁷), MYC (z = 4.87, -value = 2.60 ×10⁻²⁴), Lh (z = 4.87, p-value = 6.30 ×10⁻¹⁷), DDX3X (z = 3.13, p-value = 9.48 ×10⁻¹⁷), TCR (z = 2.91, p-value = 1.94 ×10⁻¹⁰), CD3 (z = 3.16, p-value = 1.40 ×10⁻⁷), NFE2L2 (z = 2.97, p-value = 0.0001), and both CD8A and CD8B (z = 4.15, p-value = 0.0001), whereas highly inhibited regulators included LARP1 (z = −5.10, p-value = 7.19 ×10⁻⁴⁴) and RICTOR (z = −4.24, p-value = 3.80 ×10⁻¹⁷) ( Table 2 ). Molecules in the WM1 module targeted by LARP1, MYCN, and TCR included the eigen-protein RPS4X. The carbohydrate-responsive element-binding protein (ChREBP) encoded by MLXIPL is a basic helix-loop-helix leucine zipper transcription factor of the MYC/MAX/MAD superfamily. Its activation promotes cancer cell proliferation and inhibits cancer cell apoptosis by promoting aerobic glycolysis in various cancers (21).

CD3, CD8A, CD8B, and TCR (T-cell receptor) are drivers of the immune system and the immune defense against intracellular pathogens, such as viruses and bacteria, as well as tumors. CD8 is expressed on the surface of cytotoxic T lymphocytes (CTLs), dendritic cells, macrophages, monocytes, and NK cells and exists as either a homodimer (two CD8A chains) or heterodimer (CD8A and CD8B chains), which are co-receptors of the TCR-CD3 complex and binds to MHC-I to activate T cells. CD8A recruits T-cell-specific protein tyrosine kinase (LCK) or linker for activation of T-cells family member 1 (LAT) to activate downstream signaling pathways, whereas CD8B increases the avidity of the CD8/MHC/TCR complex to increase IL-2 production.

Highly inhibited La-related protein 1 (LARP1), the master regulator of the cap-dependent top mRNA translation, strongly suggests the inhibition of protein synthesis via cap-dependent mRNA translation or the activation of the cap-independent, IRES-mediated translation of mRNA subsets encoding oncogenic proteins, such as HIF1A, and MYC (22).

Regarding the WM2 module, highly activated upstream regulators included MYC, NFE2L2, MLXIPL, IFNG (z = 4.51, p-value of overlap = 1.39 ×10⁻⁸), NFκB (complex) (z = 3.87, p = 8.69 ×10⁻⁵), FGFR1 (z = 7.46, network bias-corrected p = 0.0005), TCRD (T-cell receptor delta chain) (z = 4.57, p = 0.0001), CTNNβ-LEF1 (z = 4.52, p = 0.0001), CTNNβ1 (z = 4.40, p = 0.0001), CD80 (z = 7.31, p = 0.019), CD86 (z = 6.06, p = 0.019), hypoxia-induced factors Hif (z = 4.55, p = 0.0001), and LCK (z = 3.76, p = 0.0219), whereas highly inhibited proteins included LARP1 and RICTOR, which were common with the WM1 module, and CLPP (z = −4.47, p = 0.0001) and LILRB2 (z = −4.24, p = 0.0007) ( Table 2 ).

FGFR1 is a cell-surface receptor for fibroblast growth factors that regulate embryonic development, cell proliferation, differentiation, and migration. Its amplification has been frequently reported in NSCLC (23). Activation of CTNNβ1 and CTNNβ-LEF1 implied activation of Wnt/β-catenin signaling in which β-catenin is stabilized and accumulated as free β-catenin in the cytosol and is subsequently translocated into the nucleus and activates the TCF (T-cell factor)/LEF (lymphoid enhancer factor)-dependent transcription of Wnt target genes, key factors in cell proliferation and invasion (24, 25). Proto-oncogene LCK encodes lymphocyte-specific protein tyrosine kinase LCK, an Src family tyrosine kinase, which plays an important role in TCR-linked signal transduction. CD28-induced LCK activation is an important mediator of T-cell activation. Recently, mathematical models for LCK autophosphorylation suggest that LCK is involved in the early stages of T-cell activation and response derived from PD-1 signaling suppression (26). It also potentially mediates PD-1-induced inhibition of early TCR signaling (27), which is important to cancer immune checkpoint therapy. LILRB2 is Leukocyte immunoglobulin-like receptor subfamily B member 2 (also known as ILT-4, immunoglobulin-like transcript 4) is a tumor immune checkpoint molecule. Caseinolytic protease P (CLPP) is involved in the mitochondrial unfolded protein response and cellular bioenergetics (28).

Interestingly, all of Hif (complex), TRCD, FGFR1, IFNG, CD80, CD86, LCK, and LILRB2 target molecules including HLA-G in the WM2 module, whereas LILRB2 targets module molecules including both HLA-G and eigen-protein TPM4 ( Table 2 ). Schreiber et al. demonstrated that the effects of IFNG are generally pro-tumorigenic during the immune escape stage of cancer immunoediting, in which IFNG increases inhibitory immune checkpoint molecules, including HLA-G, to promote the formation of a tolerant immune microenvironment (29).

HLA-G exerts inhibitory effects on both innate and adaptive effectors through direct binding to inhibitory receptors, such as immunoglobulin-like transcripts, LILRB1 (also known as ILT-2), and LILRB2, with high affinity compared with other HLA class I molecules. HLA-G can inhibit many players associated with the anti-tumor response throughout early and late tumor stages, although CTLA-4 and PD-1 are predominantly expressed in higher tumor grades (30, 31). The networks associated with LILRB2 inhibition include the activated hub regulator, toll-like receptor 4 (TLR4, also known as CD284). TLR signaling activates various signaling molecules, including nuclear factor κB (NFκB), extracellular signal-regulated kinase (ERK/JNK/p38), and induces the synthesis of immunologic factors including IL-6, IL-12, PD-L1, and HLA-G. This results in the resistance of tumor cells to CTL attack and tumor cell immune evasion (32). Highly expressed HLA-G and highly inhibited LILRB2 implicated involvement of immune tolerance or anti-tumor immune escape. The volcano plots revealed high expression of HLA-G associated with the SPA subtype ( Figures 2A–C ). A web-based survival analysis (KMplot) for mRNA data of lung adenocarcinoma (n  =  719) indicated that high HLA-G expression was significantly associated with poor OS (log-rank test p  =  1.0 × 10⁻¹¹; hazard ratio: 3.1 (2.2–4.36)) (https://kmplot.com/analysis/) ( Supplementary Figure S3 ) (33). Yan et al. interestingly showed a predominant expression of soluble HLA-G (sHLA-G) in tumor cells not from squamous cell carcinoma but adenocarcinoma of lung (34).

Lung adenocarcinoma genomic data from the MSKCC database (Stage I to III: n  =  604) was analyzed for prognostic markers associated with the predominant histologic subtypes in lung adenocarcinoma using the Memorial Sloan Kettering–Integrated Mutation Profiling of Actionable Cancer Targets (MSK-IMPACT) platform (5). We extracted data for the LPA (n = 88), MPA (n = 37), and SPA (n = 68) subtypes, and the genomic alteration profiles were visualized using the cBioPortal for Cancer Genomics (https://www.cbioportal.org/) ( Supplementary Figure S4 ). The oncoprints showed genomic alterations characteristic of the three subtypes as follows: LPA, TP53 (17%), EGFR (43%), and KRAS (33%), whereas EGFR and KRAS were mutually exclusive (q < 0.001); MPA, TP53 (41%), EGFR (22%), and KRAS (32%); SPA, TP53 (75%), which was the highest among the three subtypes, EGFR (22%), and KRAS (38%).

In connection with the above genomic alterations, upregulated expressions of HLA I class molecules such as HLA-G and HLA-B were observed in proteomic analysis. ( Figures 2A, B ). It has been reported that there is a positive correlation between mutational burden and the expression of most HLA class I molecules (35, 36), and HLA class I molecules are associated with the innate immune response as ligands of inhibitory killer cell immunoglobulin-like receptors (KIRs) of Natural Killer (NK) cells (37).