Graphical Abstract. The flowchart of the study is presented.

Lung cancer currently ranks as the second most common malignant tumor globally and the leading cause of cancer-related mortality (1). Pathologically, it can be classified into non-small cell lung cancer (NSCLC) and small cell lung cancer. NSCLC accounts for 80% to 85% of newly diagnosed lung cancer cases. The primary treatment modalities for NSCLC include surgical resection, targeted therapy against tumor-driving genes, conventional chemotherapy, and emerging immune checkpoint inhibitors (2). Approximately 50% of NSCLC patients present with distant organ or lymphatic metastasis at the time of diagnosis (3). For such patients, aside from stage IIIA patient s who may undergo surgical treatment after assessing the efficacy of neoadjuvant therapy, most unresectable stage III and IV NSCLC patients rely on drug therapy as first-line treatment (4). These tumors can be furtherly categorized into NSCLC with driver gene mutations and NSCLC without driver gene mutations. For NSCLC patients with targetable driver gene mutations, the first-line treatment regimen typically consists of targeted therapy against driver genes combined with anti-angiogenic drugs and conventional chemotherapy (5). Conversely, treatment regimens for driver gene-negative patients are stratified based on Eastern Cooperative Oncology Group Performance Status (ECOG) and may include combinations of immune checkpoint inhibitors and conventional chemotherapy or monotherapy with immune checkpoint inhibitors. The advent of immune checkpoint inhibitors, primarily targeting Programmed cell death 1 ligand 1 (PD-L1)/Programmed cell death protein 1 (PD-1), has transformed the treatment landscape for stage III and IV NSCLC (6). For instance, patients treated with pembrolizumab in combination with platinum-based chemotherapy and paclitaxel exhibited an increased overall survival from 10.7 months to 22 months compared with those treated solely with conventional chemotherapy (7). However, it is noteworthy that many patients still do not benefit from immunotherapy and experience tumor progression. Currently, reliable biomarkers and predictive models for identifying patients likely to benefit from immunotherapy before treatment initiation remain lacking.

Metabolomics and proteomics analyses based on human serum samples have provided reliable biomarkers for early diagnosis, treatment plan selection, efficacy monitor and prognosis and helped build many prediction models as biomarkers selection based on peripheral blood has many advantages including low heterogeneity, small trauma and convenience for continuous monitoring (8). As a research tool for comprehensive analysis of changes in endogenous small molecule metabolites occurring in an organism following internal or external stimuli, metabolomics is considered as an extension and endpoint of genomics and proteomics, where small changes in expression levels at the gene and protein level can be amplified at the metabolite level, and thus metabolomics is regarded as a “biochemical phenotype” of the overall functional status of an organism, and is well suited for biomarker screening (9). And as the most downstream of various pathophysiological activities, changes in metabolites also reflect changes in the sensitivity of non-small cell lung cancer to immune checkpoint inhibitor therapy. Proteins are the bearers of various cellular functions, and resolving the spatiotemporal specificity of proteins is a key molecule to understand the heterogeneity of microenvironments and the characteristics of life occurrence and development in tissues and diseases (10). With the promotion of mass spectrometry instrumentation, quantitative histological testing of large-scale samples has begun to gain gradual popularity, and proteomic testing of body fluids, such as plasma, urine, and saliva, for patients with malignant tumors has become an important research method for early diagnosis, pathway discovery, and efficacy monitoring of tumors (11).

In this study, we collected 47 serum samples from unresectable stage III or IV NSCLC patients who received Pembrolizumab treatment and metabolomic and proteomic analyses were performed according to the classification based on follow-up results to explore biomarkers that can be used for early prognosis and provide potential targets for precision therapy (Graph Abstract).

5 ml blood from patients who met the inclusion criteria was collected and placed in Ethylene Diamine Tetraacetic Acid (EDTA) blood collection tubes. The tubes were left to stand for 2 hours at 4 degrees Celsius. After centrifugation for 15 minutes at 3,000 rpm and 4 Celsius degrees, the supernatant was taken in a 1.5ml centrifuge tube, labeled with sample information and frozen at -80 degrees Celsius refrigerator for storage. The collection of peripheral blood plasma specimens of the patients was proceeded before immunotherapy, and the efficacy of the immunotherapy was evaluated after patients had received two courses of treatment. Based on the results of the high-resolution chest CT before and after the treatment, the patients were classified into three groups according to the proportion of the reduction of the tumor volume, complete or partial response group, stable disease group, and progressive disease group according to Response Evaluation Criteria in Solid Tumors 1.1 (12). Patients with disappearance of all target lesions or at least a 30% decrease in the sum of diameters of target lesions were grouped into PR/CR, patients who witnessed at least a 20% increase in the sum of diameters of target lesions were grouped into PD group and other patients were grouped into SD group.

Thaw serum at 4 degrees Celsius, gently vortex and mix; take 50 microliters of serum in a centrifuge tube; add 200 microliters of internal standard solution; vortex for 1 minute; let stand for 2 minutes; centrifuge at 4 degrees Celsius for 15 minutes at 13,000 rpm; take 100u microliters of supernatant from the wall in the feeder vials and wait for measurement. Take 10 microliters of each supernatant and mix, level and centrifuge for two minutes, immediately take the supernatant into the feed vial as a Quality Control (QC) sample to prevent protein precipitation from floating. The remaining serum was immediately returned to the -80 degrees Celsius refrigerator for refreezing.

A Vanquish ultra-high performance liquid chromatography system, Q-Exactive combination Orbitrap mass spectrometer and C18 column from Thermo Fisher Scientific China Ltd. were used. The temperature was set at 35°C, mobile phase A consisted of 0.1% formic acid and water, and mobile phase B consisted of 0.1% formic acid and acetonitrile at a flow rate of 0.4 ml/min, and the injection volume of each sample was 2 μl, with QC sample inserted in between every eight real samples. The elution gradient was set as 0-2 min: 5-5% (B); 2-13 min: 5-95% (B); 13-15 min: 95 -95% (B); 15-15.1 min: 95-5% (B); 15.1-20 min: 5-5% (B). The ion source for mass spectrometry was a heated electrospray ionization source, HESI, with the voltage set to 3.8 kVt in positive ion mode and 3.2 kVt in negative ion mode, a capillary temperature of 320 degrees Celsius, a mass-to-charge ratio between 100 and 1500, and a shielding gas of 30 liters/minute. After completing the mass spectrometry analysis, the raw data from the Q-Exactive mass spectrometer were exported and proceeded using Thermo Fisher’s Compound Discoverer 3.3 software.

In the first step, spectra were selected from the raw data for each polarity and the retention times of the chromatographic peaks were aligned. Unknown compounds with spectral peak intensity values over 10,000 were then detected and spectra. In the second step, missing values were filled with Compound Discoverer 3.3 and compounds were identified by different types of databases. Mzcloud was used to annotate compounds on MS/MS at tolerances within 10 ppm. Chemspider contains the BioCyc, the Human Metabolome Database (HMDB), and the Kyoto Encyclopedia of Genes and Genomes (KEGG) databases for precise mass-based annotation of features with a mass tolerance of 5 ppm and an S/N threshold of 1.5. In the third step, SERRF QC correlation was applied to reduce batch effects in all samples. In the fourth step, raw data from blank samples were used to label background compounds and the peak area of the samples/blanks was set to 5, enabling the use of background compounds to filter the background signal. Finally, a list of detected compounds, including exact molecular mass-to-charge ratio, retention time, compound name, and peak area, is exported in Excel format in QC normalization mode.

The exported compound names were manually checked before using MetaboAnalyst as a target list for pathway enrichemnt analysis and random forest analysis, both of which are online data processing tools. In addition, peak area lists were uploaded into SIMCA-P 14.1 software (Umetrics, Malmo, Sweden) for PCA analysis to assess stability of data between samples. A rank-sum test was performed on all metabolites detected by Compound Discoverer 3.3 to screen for metabolites that differed in expression levels between groups. Pathway enrichment analysis based on KEGG database was performed for the differential metabolites screened between groups using the pathway analysis module on MetaboAnalyst website.

In this study, the iST sample pretreatment kit (PreOmics, Germany) was taken for protein pretreatment. 2 µl of plasma samples were taken and 50 µl of lysate was added, centrifuged at 1000 rpm for 10 min at 95°C and then left to cool down to room temperature, trypsin digestion buffer was added, and the reaction was incubated with oscillation for 2 h at 37°C, 500 rpm, and the reaction was terminated by adding 100 μl of buffer. The peptides were desalted using the iST cartridge in the kit, eluted twice with 100 μl of elution buffer, then vacuum-dried and stored at -80 degrees Celsius.

The samples were analyzed by LC-MS/MS using an AUR3-15075C18 analytical column (15 cm*75 μm, 1.7 μm) with a gradient of 30 min, a column temperature of 50°C, and a column flow rate of 400 nl/min. The B-phase consisted of 80% acetonitrile with 0.1% formic acid, and the gradient was started from 4% of B-phase, increased to 28% within 15 min, and to 44% within the next 4 min. 4 min to 44%, and then to 90% in the next 4 min, maintained for 3 min and then equilibrated at 4% for 4 min. The mass spectrometer was set to diaPASEF mode with a scanning range of 349-1229 m/z and an isolation window width of 40 Da. During the PASEF MSMS scan, the collision energy increased linearly with ion mobility from 59 eV (1/K0 = 1.6 Vs/cm2) to 20 eV (1/K0 = 0.6 Vs/cm2).

DIA data were analyzed using the Spectronaut18 default parameters, i.e., BGS Factory Settings (default), with a sequence database of uniprot-Homo_sapiens (version 2022, 20610 entries), set for Trypsin zymolysis, and fixed modifications of the search library parameters as: Carbamidomethylation (C) 57.02 and variable modification is: Oxidation (M) 15.99, whereas according to the iRT peptide software the retention time and mass window can be automatically corrected to automatically determine the ideal extraction window. The criteria for protein characterization were Precursor Threshold 1.0% FDR and Protein Threshold 1.0% FDR; while the Decoy database was generated using a mutated strategy, i.e., the sequence of a random number of amino acids was selected to be disrupted, and the number of disrupted amino acids should be greater than or equal to two and less than half of the total length of the peptide; Spectronaut was automatically corrected and the data were normalized using a local normalization strategy, while peptides less than 1.0% FDR were quantified using MaxLFQ to complete the proteome (13).

The quality control of mass spectrometry data included liquid phase system, mass spectrometry system, qualitative results, quantitative results, etc. The parameters of quality control included missed cut, specific enzyme cleavage, non-specific modification, half-peak width, distribution of ion out time, distribution of ion charge number, and mass axis shift. The samples of the beneficiary group were compared with those of the non-beneficiary group, and the mean of the relative quantitative values of each group of samples in the comparison sample pair was calculated for each white, and the ratio of the mean values of each group of samples in the comparison sample pair was the Fold Change (FC). The p-value was calculated by t-test to determine the significance of the difference, and FC>1.2 or FC<-0.83 and p<0.05 were used as the screening criteria for differential proteins. The protein quantification form was exported according to the library search software, and data preprocessing was required before screening for differential proteins. Data preprocessing includes contaminant library protein removal, missing value filtering and filling. The missing value filtering and filling adopt the global missing value rejection and filling strategy, which is, the samples were grouped according to the group, in each group of samples, the protein with more than 50% missing values was considered to have 0 expression in that group of samples, and the protein values with missing value less than 50% were filled with the very small value, and the proteins with the percentage of missing value more than 50% were filtered out when the percentage of missing value in all groups was more than 50%.

While statistical analysis such as differential proteins include Principal Component Analysis (PCA), Partial Least Squares Discriminant Analysis (PLS-DA), Hierarchical Cluster Analysis (HCA) and so on. A variety of bioinformatics analysis tools including GO analysis, KEGG analysis, and protein-protein interaction networks were applied in proteomics analysis.

R 4.2.2 was used for data analysis, metabolomics and proteomics were screened for predictor variables using random forest and stepwise regression respectively, 29 differential metabolites and 23 differential proteins between the beneficiary and non-beneficiary groups were included in the independent variables available for screening, a randomized stratification method was adopted to divide the data into training and testing sets according to 7:3, based on the training set, logistic regression was used to respectively Based on the training set, logistic regression was used to construct clinical risks, metabolomics, proteomics, metabolomics and proteomics prediction models, and the ROC (Receiver Operating Characteristic) curve was used to evaluate the differentiation of the models, and the AUC (Area Under Curve) was used to evaluate the effectiveness of the models.

IPA (Integrate Pathway Analysis) software was used to analyze the pathway enrichment of 27 differential metabolites and 23 differential proteins simultaneously, and the corresponding dual-omics combined pathway maps were obtained.

From October 2022 to October 2023, a total of 47 patients were enrolled in this study, including 42 male patients and 5 female patients. 35 patients were included in the benefit group while 12 were included in the no-benefit group. 30 of 35 patients belonging to benefit group were classified into stable disease (SD) group and the other 5 patients were classified into remission group which consisted of complete response (CR) group and partial response (PR) group. The no-benefit group consisted of progressive disease (PD) group and death group. The details of each group are shown in Table 1.

Table 1. The detailed information of participants.

A total of 998 metabolites identified after initial process using Compound Discoverer 3.3 were checked according to the HMDB database, Chemspider database, and KEGG database, and a total of 392 of these metabolites that were endogenous to the human body were manually identified. Subsequently, SIMCA- P 14.1 software was used to further analyze the cleaned data, and Principal Component Analysis (PCA) was applied to analyze the stability of the samples, according to the PCA-X in Supplementary Figure 1. As it was shown, the QC samples have good aggregation in both positive and negative modes which proved that the mass spectrometer has reached a good stability level during the metabolomics analysis.

Subsequently, the rank-sum test was performed on the peak values of each human endogenous metabolite for each sample in both the beneficiary and non-beneficiary groups, and 27 differential metabolites were obtained, which are shown in Table 2. The cluster analysis was also performed on the differential metabolites (Figure 1A).

Table 2. The differential metabolites between the beneficiary and non-beneficiary groups.

Figure 1. Heat map showing the results of the cluster analysis of differential metabolites between the progressed and mitigated groups (A). Heat map showing the results of the cluster analysis of differential metabolites between the progressed and stabilized groups (B). Heat map showing the results of the cluster analysis of differential metabolites between the mitigated and stabilized groups (C). Heat map showing the results of the cluster analysis of differential metabolites between the benefited and non-benefited groups (D).

As it can also be seen in Figure 2, a two-by-two rank sum test comparison between the remission, stabilization, and progression groups revealed that there were 25 differential metabolites between the remission group and the progression group (Figure 1B), only 5 differential metabolites between the stabilization group and the progression group (Figure 1C), 11 human endogenous differential metabolites between the remission group and the stabilization group (Figure 1D).

Figure 2. (A) Volcano plot of protein differences between the progressed and stabilized groups. (B) Volcano plot of protein differences between the relieved and progressed groups. (C) Volcano plot of protein differences between the relieved and stabilized groups. (D) Volcano plot of protein differences between the benefited and non-benefited groups. (E) Heat map showing the results of the cluster analysis of differential proteins between the beneficiary and non-beneficiary groups. (F) Heat map showing the results of the cluster analysis of differential proteins between the progressed and stabilized groups. (G) Heat map showing the results of the cluster analysis of differential proteins between the progressed and beneficiary groups. (H) Heat map showing the results of the cluster analysis of differential proteins between the mitigated and stabilized groups.

Pathway enrichment analysis was performed on these 27 differential metabolites using the MetaboAnalyst website, and it was found that a total of 11 pathways were enriched (Figure 2A), among which were the metabolism of tryptophan and aminoacyl tRNA biosynthesis pathways had most differential metabolites on them, and this result was also verified in the subsequent IPA analysis. For some of the metabolites enriched on the relevant pathways we will validate this in the next section of our work.

The three two-by-two comparisons had the highest number of differential metabolites between the PR and PD groups, and the differences in response to immune checkpoint inhibitor therapy for non-small cell lung cancer were also the greatest between these two groups of patients, with the highest number of pathways enriched. As shown in Figure 3, based on the pathway enrichment analysis of the three comparisons, it was found that the pathways enriched for differential metabolites between the remission group and the progression group were mainly glutathione metabolism and unsaturated fatty acid biosynthesis, whereas the metabolic pathways enriched for endogenous differential metabolites in the human body between the remission group and the stabilized group were mainly porphyrin metabolism and tryptophan metabolism. However, due to the lack of differential metabolites between stabilized group and progression group, corresponding pathway enrichment was not successful.

Figure 3. (A) Pathway enrichment results of differential metabolites between the beneficiary and non-beneficiary groups. (B) Pathway enrichment results of differential metabolites between the remission group and the stable group. (C) Pathway enrichment results of differential metabolites between the remission group and the progression group.

Preliminary analysis of the proteomic mass spectrometry data was performed using Proteomics Discoverer software. 585 endogenous plasma proteins were identified in this study, and 23 proteins with statistically significant differences in expression levels between the beneficiary group and the non-beneficiary group, of which 20 proteins had higher expression levels in the non-beneficiary group, while only three proteins, namely GDF-9B (Growth/differentiation factor 9B), Stress-70 protein and Protein HEG homolog 1, had higher levels in the beneficiary group, respectively. GDF-9B, Stress-70 protein and Protein HEG homolog 1 were expressed at higher levels in the beneficiary group, which were respectively derived from BMP15 GDF9B, HSPA9 GRP75 HSPA9B mt-HSP70 and HEG1 KIAA237 gene translation as shown in Table 3. 27 differential proteins were identified between the progression group and the stabilization group, 56 differential proteins between the stabilization group and the remission group, and 38 differential proteins between the progression group and the remission group, and the results of the cluster analysis are shown in Figure 4.

Table 3. The differential proteins between the beneficiary and non-beneficiary groups.

Figure 4. The top 20 entries with the smallest p-value under the three tertiary classifications of biological process (A), molecular function (B), and cellular components (C) in the GO analysis for differential proteins in the enrichment results under the three levels of classification and the GO level where they are located, and the horizontal coordinate is the percentage of differentiated proteins enriched by the entry.

PLS-DA analysis was performed on the detected proteins by subgroups to show the trend of aggregation or dispersion of the samples within and between groups. As shown in Supplementary Figure 2 the PLS-DA analysis of the benefit group and the non-benefit group showed a better trend of separation between the two groups and a better trend of aggregation between the samples within the groups; if grouped according to the progression group, the stabilization group, and the remission group, the comparisons of any two of the three groups including the progression group versus the stabilization group, the stabilization group versus the remission group, and the progression group versus the remission group were included, and these comparisons were all in the PLS-DA analysis showed a better trend of intra-group aggregation and inter-group separation, while the PLS-DA between the three groups proved that the inter-group differences between the three groups were significant while the intra-group samples were highly aggregated.

Subsequently, a volcano plot was drawn for the detected proteins, which is a scatter plot taking the log2 value of the fold change as the horizontal axis and the -log10 value of the p-value as the vertical axis, and proteins with up-regulation of expression levels were labeled in red, proteins with no change in expression levels were labeled in gray, and proteins with down-regulation of expression levels were labeled in blue, and as can be seen in Figure 3, even though the expression levels of most detected proteins were not statistically changed in the comparison of the groups, there were some proteins that still demonstrated an up- or down-regulation trend in the comparison. In the comparison between the benefit group and the non-benefit group, only 3 differential proteins were expressed at higher levels in the benefit group, while 20 proteins showed a tendency to be expressed at lower levels in the benefit group. The hierarchical clustering analysis, in which the color shades are used to represent the magnitude of protein expression, is shown as a heat map, which again confirms the trend of expression levels of differential proteins in the comparison of groups in the volcano plot.

From the comparison between the PD group and the SD group, it can be seen that only 4 differential proteins showed a downward trend in the PD group while the remaining 23 differential proteins showed an upward trend in the PD group; whereas from the comparison between the PR group and the PD group, the difference in the number of differential proteins with upward and downward trends was small, with 24 differential proteins with a higher expression level in the PD group and 14 with a lower expression level in the PD group; the difference in expression level between the PR group and the progression group was small. In the comparison between the PR group and the SD group, there were 36 different proteins with increased expression levels in the PR group and 20 proteins with decreased expression levels in the SD group.

Subsequent bioinformatics analysis based on differential proteins included GO function analysis, KEGG pathway analysis with STRING-DB protein interaction network. Differential proteins in the beneficiary and non-beneficiary groups were firstly categorized at the first level, which included Biology Process (BP), Molecular Function (MF), Cellular Component (CC), and Complement Activation, Antigen Binding, and Vesicles were the three pathways or cellular components that accounted for the highest number of differential proteins in the relative categorization. Complement activation, the bioprocess with the highest percentage of enriched proteins, is located in the seventh level of classification, and three differential proteins are located in the pathway corresponding to this entry; antigen binding, the molecular function with the highest percentage of enriched proteins, is located in the second level of classification, and four differential proteins are located in this pathway, and the rest of the pathways that are enriched with more differential proteins and with lower p-values are immunoglobulin receptor binding, ligase. The other pathways with low p-values were immunoglobulin receptor binding, ligase activity to form carbon-sulfur bonds, and heat shock protein binding, all of which had 2 differential proteins enriched in the relevant pathway; whereas in the cellular components the structure of vesicles was enriched in 11 differential proteins, while the structures of extracellular vesicles, extracellular exocytosis and extracellular exocytosis were enriched in 10 proteins, which may be related to the fact that the samples came from peripheral blood plasma after centrifugation of removed hematopoietic cells (Figure 4).

Under the primary classification of biological processes, the pathway of complement activation was enriched with the highest number of differentially differentiated proteins and the lowest p-value, with the highest confidence; the pathways with lower p-value and higher confidence also included protein recombination, De novo protein folding, erythrocyte differentiation, and lung development, which were also characterized by lower FDR values; and in the secondary classification of biological processes, cellular component organization/cellular component, CC, antigen binding, and vesicle were the pathways or cellular components that accounted for the highest number of differentiated proteins, respectively, in all three categories (Figure 5). In the secondary classification of biological processes, the pathway of cellular components, organization/biosynthesis was the pathway enriched to account for the largest number of differential proteins. Under the classification of molecular function, the antigen-binding pathway had a lower p-value and better FDR value, whereas protein-binding was the secondary classification under molecular function with the highest proportion of differentially enriched proteins, and ubiquitin-like modifier-activated enzyme activity was the pathway with the lowest p-value. Under the primary classification of cellular composition, vesicles were the cellular structure with the highest percentage of differentially differentiated proteins enriched to the structure, and the rest of the structures with a high percentage of differentially differentiated proteins were extracellular exosomes and extracellular regions, which also had the lowest p-value and the best FDR value under this classification. The top 20 entries with the lowest p-values in the enrichment results under the three classifications of biological processes, molecular functions, and cellular constituents were plotted against their GO tiers in a horizontal bar graph, with the horizontal coordinates of the bar graph being the percentage of differentially differentiated proteins enriched to that entry.

Figure 5. Bar graphs of GO enrichment. (A–C) P-values of the top 20 pathways with the most differential proteins enriched between the beneficiary and non-beneficiary groups. (D–F) FDR values of the top 20 pathways with the most differential proteins enriched under the three tertiary classifications. (G–I) Horizontal histograms of the top 20 tertiary or quaternary entries with the lowest p-value in the enrichment results under the lowest p-value secondary classification with the GO tier in which they were selected.

And as seen in Figure 6, the KEGG enrichment results were visualized in a similar way, the vertical coordinates of the front bar with the smallest p-value were counted to indicate the percentage of the proteins in the pathway to the number of differential proteins, and the darker the color, the smaller the p-value; whereas, the KEGG seven first-level macroclasses analyzed and their subordinate second-level classifications under which those pathways enriched by differential proteins belong to are labeled by means of the horizontal bar charts. Propionic acid metabolism, starch and sucrose metabolism belonged to the carbohydrate metabolism secondary classification under the metabolic classification, sulfur metabolism belonged to the energy metabolism secondary classification under the metabolic classification, and pancreatic secretion and other pathways belonged to the digestive system classification under the organic system; whereas pancreatic secretion and other pathways enriched in two differential proteins, most of the remaining multiplexes were enriched only in one differential protein, as the amount of proteins enriched by each pathway was relatively small, and therefore, the number of differential proteins enriched in each of these pathways was not large. The KEGG pathway analysis of differential proteins in the benefit group versus the non-benefit group was less reliable and had higher p-values for each pathway due to the small number of proteins enriched for each pathway. Whereas the entries with the highest number of differential proteins enriched in the comparison of the progression group with the stabilization group were metabolic pathways, the entries with the highest number of differential proteins enriched in the comparison of the remission group with the progression group were phagosomes.

Figure 6. KEGG enrichment results. (A) The top 4 pathways with the smallest p-value, the vertical coordinates indicate the percentage of proteins in the pathway to the total number of differential proteins, and the darker the color indicates the smaller the p-value. (B) The vertical coordinates on the left show the specific metabolic pathways enriched, and the vertical coordinates on the right indicate the abbreviation of the primary and secondary taxonomic names. (C) P-value of the enriched metabolic pathways in order from the outside in. of the number of proteins enriched in the first 4 entries followed by a pie chart showing the percentage of the number of proteins in each entry. (D) Visualization of KEGG enrichment results in bubble charts, with the horizontal coordinates indicating the percentage of proteins in the pathway to the total number of differential proteins, and the color indicating the smaller p-value, and the size of the circle indicating the number of differential proteins in the pathway. (E) KEGG enrichment results of the differential proteins between the progression group and the stabilization group. (F) KEGG enrichment results of differential proteins between the remission and progression groups.

The protein interaction analysis of differential proteins showed that, as shown in Figure 7, HSPA9 was the core of the interaction network composed of differential proteins, and the rest of the key nodes were CTRB1, TTC25, TUBB1, etc., which were closely related to the mutation of tumor driver genes, especially K-Ras mutation.

Figure 7. Results of String analysis of differential proteins. (A) Protein interaction network of differential proteins between the benefited and non-benefited groups. (B) Protein interaction network of differential proteins between the progressed and stabilized groups. (C) Protein interaction network of differential proteins between the relieved and progressed groups. (D) Protein interaction network of differential proteins between the relieved and stabilized groups.

In contrast, as shown in Figure 8A, pathway enrichment analysis using IPA (Integration pathway analysis, QIAGEN, Dusseldorf, Germany) software showed that the expression levels of genes such as MYC, CCND1, and TP53 were important factors affecting the sensitivity of NSCLC to immune checkpoint inhibitor therapy.

Figure 8. (A) IPA pathway enrichment analysis of differential metabolites and differential proteins between the beneficiary and non-beneficiary groups. (B) ROC curves in the test set based on predictive models constructed from biomarkers screened in metabolomics, proteomics and multi-omics. (C) ROC of prediction model based on clinical risks. (D) ROC of prediction model based on differential metabolites. (E) ROC of prediction model based on differential proteins. (F) ROC of prediction model based on both differential metabolites and proteins.

After screening by random forest and stepwise regression, metabolomics screened two differential metabolites, 5-sulfomethylfurfural and o-aminobenzoic acid, as biomarkers, while proteomics screened two differential proteins, Immunoglobulin heavy variable 1-45 and Microfibril-associated glycoprotein 4, as biomarkers. As shown in Figures 8B–E the early prediction model constructed based on two metabolic biomarkers, 5-sulfomethylfurfural and o-aminobenzoic acid, had an AUC of 0.977 in the test set while the early prediction model constructed based on two protein biomarkers, Immunoglobulin heavy variable 1-45 and Microfibril-associated glycoprotein 4, had an AUC of 0.875 in the test set. The early prediction model constructed from clinical factors had an AUC of 0.659 in the test set while the early prediction model constructed from dual omics biomarkers had an AUC of 0.955.

Tryptophan is an essential amino acid, and its metabolic pathway is the most metabolized and reliable metabolic pathway among the metabolic pathways enriched by differential metabolites in the beneficiary and non-beneficiary groups, while the metabolism of tryptophan is closely related to the treatment of immune checkpoint inhibitors in tumors, and its catabolic metabolism can reduce the body’s immune system whereas tryptophan itself can enhance the activity of CD8+ T-cells, increase CD8+ T-cells in tumor tissue and increase the infiltration of CD8+ T cells in tumor tissues, thus helping to induce apoptosis of tumor cells and slow down the growth of lung cancer (14). There are three metabolic pathways for tryptophan, including the indole pathway in the intestinal microbiota, the serotonin system pathway in intestinal chromaffin cells, and the kynurenine pathway in immune cells and the intestinal wall (15). All three pathways are relevant to immune checkpoint inhibitor therapy of tumors (16). Indole metabolites and kynurenine interact with the aryl hydrocarbon receptor to induce T regulatory cell differentiation, limit immune responses in Th17 cells, Th1 cells and produce anti-inflammatory mediators (17). Kynurenine leads to a decrease in CD8+ T cells infiltrating tumor tissue and mediates immune escape of tumor cells (18). Also, the serotonin system increases tumor cell proliferation and metastasis, while indole metabolites significantly reduce tumor growth (19).

Another pathway that was enriched by multiple metabolites in the comparison of the beneficiary and non-beneficiary groups is the biosynthesis of aminoacyl-tRNA. Aminoacyl-tRNA synthesis works by recognizing and catalyzing specific amino acid linkages to homologous tRNAs via aminoacyl-tRNA synthetases, which results in the precise matching of amino acids to tRNAs containing the corresponding paradigm codons, helping to achieve the accurate synthesis of proteins (20). Aminoacyl-tRNA synthetases are key enzymes in the translation process of mRNA, and the 20 essential amino acids correspond to 20 aminoacyl-tRNA synthetases, but there are fewer reports in the literature on the therapeutic effects of aminoacyl-tRNAs and anti-PD-L1/PD-1 therapy in tumors (21). Specifically, both threonine and tryptophan were enriched to aminoacyl-tRNA biosynthesis, and it is noteworthy that threonine was also an important differential metabolite between the beneficiary and non-beneficiary groups, and that serine/threonine is an important site in the WNT/β-catenin pathway during tumorigenesis, whereas the corresponding serine/threonine kinase plays an important role in this pathway (22). Transcriptomics revealed that serine/threonine kinase expression was negatively correlated with the expression of immune response markers in CD8+ T cells and the infiltration of dendritic cells, whereas further studies revealed that serine/threonine kinase expression levels were higher in patients with malignant tumors that were insensitive to anti-PD-L1/PD-1 therapy, and that genes associated with the WNT/β-catenin pathway and the MYC, a target gene of WNT, both had higher expression levels (23, 24). Knockdown of the serine/threonine kinase gene induced tumor shrinkage and increased immune cell infiltration in tumor tissues of malignant tumors ineffective for anti-PD-1 treatment in animal experiments, while synergistic effects of combining anti-PD-1 immune checkpoint inhibitors with anti-serine/threonine kinase targeting drugs could be observed in animal models, and synergistic effects were observed for the combination of Nivolumab and new therapy consisting of Nivolumab in combination with the dual anti-serine/threonine kinase inhibitor KPT-9274 is already undergoing clinical trials (25).

Immunotherapy for non-small cell lung cancer lung cancer has a close management with various protein pathways, and the activation of complement is the pathway with the most enriched differential proteins in the GO analysis of differential proteins between non-small cell lung cancer immunotherapy beneficiary and non-beneficiary groups, which plays an important role in immune evasion of tumors and activation of CD8+ T-cells, especially the complement factors produced by tumor cells could regulate tumor signaling and tumor tissue growth factors (26). For example, mutations in STK11 gene are thought to be associated with poor efficacy of immune checkpoint inhibitor therapy in NSCLC, while animal experiments have confirmed that the expression of complement pathway, including C3, is significantly increased in tumor cells with STK11 knockdown, and that complement C3 in NSCLC cell lines with STK11 knockdown can lead to a decrease in xenograft tumor growth in nude mice while tumor formation in mice with normal immune function basically disappeared. Mice with normal immune function, while population experiments have shown differential expression of complement C3 in non-small cell lung cancer patients with STK11 mutations (27). Whereas tumor cell-derived complement is now a possible target for immunotherapy in NSCLC, thus enhancing the anti-tumor capacity of the complement system, combined targeted blockade of C3aR (Component 3 Antibody Receptor)/C5aR (Component 5 Antibody Receptor) with PD-1/PD-L1 immune checkpoint inhibitor therapy appeared to have an antitumor synergistic effect (28).

Immunotherapy of non-small cell lung cancer lung cancer is also closely related to the expression level of the protein, HSPA9 as the core of the protein interaction analysis network of differential proteins, which is thought to be closely linked to the mutation of the tumor driver gene K-Ras, when there is a mutation in the K-Ras G12V gene depletion of HSPA9 can lead to the inhibition of the role of mitochondrial Ca2+ one-way transport protein is reduced and thus promote the death of normal fibroblasts (29). HSPA9, as a mitochondrial molecular chaperone, is often highly expressed and mislocalized in tumor cells with aberrant activation of MEK kinase and ERK kinase, and its depletion selectively kills tumor cells with high expression of B-Raf V600E or the chimeric protein ΔRaf-1:ER, and MEK-ERK-sensitive regulatory peptides in HSPA9 are not known to be active in the mitochondria (30). ERK-sensitive regulatory peptide binding domain in HSPA9 is important for cell survival or death. MEK-ERK increases mitochondrial permeability by promoting the interaction between adenine nucleotide translocase 3 (ANT3) and the peptidyl prolyl isomerase procyclic protein D (CypD), leading to cell death, whereas the depletion of HSPA9 leads to a reduction of its inhibitory effect on the MEK-ERK pathway and promotes cell death, affecting the cell’s ability to survive or die. of HSPA9 decreases thereby promoting cell death and affecting tumor prognosis.

CCND1 gene is one of the most central and associated pathways in the joint pathway analysis of metabolomics and proteomics, and it is thought to be associated with poor prognosis of lung cancer. For example, CCND1 and FGFR1 genes showed a tendency of co-expression in squamous lung cancer patients, and FGFR1 could promote the migration and invasion of squamous carcinoma cells by up-regulating the expression level of CCND1 gene to activate AKT/MAPK signaling to promote the process of epithelial mesenchymal transition, while knockdown of CCND1 gene in squamous lung cancer cell lines could significantly inhibit the proliferation of tumor cells, invasion and the process of epithelial mesenchymal transition (31). On the other hand, a study by a Korean scholar in 2023 found that malignant tumors with TP53 mutation and high expression of CCND1 gene on chromosome 11 tended to have a poor prognosis on treatment with immune checkpoint inhibitors, and TP53 was also one of the genes with a high number of associated pathways in the combined two-omics analysis of this study. As the most common oncogenic driver of lung adenocarcinoma, mutations in the TP53 gene promoted resistance to PD-1/PD-L1 inhibitors in a mouse lung adenocarcinoma model, and mutations in the TP53 gene have been shown to be associated with lower objective tumor remission rates after receiving Nivolumab in a large population-based cohort of lung adenocarcinomas (32). The MYC genes, on the other hand, have been closely associated with the function of T cells in the tumor microenvironment (24). Tumor cells with high expression of the MYC gene inhibit the JAK/STAT pathway thereby reducing IFN-γ stimulation and thus the efficacy of anti-PD-1/PD-L1 drugs. In contrast, inhibitors targeting MYC have been shown to inhibit tumor growth in mice, increase the infiltration of immune cells in the tumor microenvironment, and up-regulate the expression level of PD-L1 in tumor tissues making tumors more sensitive to anti-PD-1 immunotherapy (33).

The efficacy of the prediction model composed of differential metabolites was significantly higher than that of the prediction