Between June 2019 and May 2021, a total of 28 GC patients were enrolled in this study. Following the completion of immunotherapy, 17 patients exhibited a favorable response to PD-1 inhibitors, while the remaining 11 patients did not manifest a similar response. As shown in Table 1, the median age of the enrolled patients was 59 years (range: 37–78 years), with 64.2% being men. Tumors had invaded the outer lining of the stomach (T4) in 85.7% of cases, and metastasis to other parts of the body (M1) was observed in 78.5% of cases. No significant association was found between the response to immunotherapy and clinical–pathological characteristics.

Archival pre-treatment tissue specimens were available for all patients. The protein was extracted from formalin-fixed, paraffin-embedded (FFPE) tissues and subjected to quantitative analysis using MS. Stringent quality control measures were implemented to ensure the reliability of the data (see Methods for details, Supplementary Figure S1). Details of the comprehensive study design are given in Figure 1.

To gain an insight into the proteomic profile differences between responder and non-responder groups, we conducted differential expression analysis. First, our principal component analysis (PCA) revealed no significant outliers and minor intergroup variance, suggesting limited distinction among the sample sets (Figure 2A). Subsequently, by employing Student’s t-test with a statistical threshold (p-value < 0.05; fold change > 1.2 or < 1/1.2), we identified 320 differentially expressed proteins (DEPs), consisting of 139 upregulated and 181 downregulated proteins in non-responders (Figure 2B; Supplementary Table S5). The most significantly altered DEPs included ERO1A, NCEH1, THEM6, JUP, DSP, LYAR, TFRC, IGF2BP2, COL12A1, CGN, C7, OGN, COL14A1, CYBRD1, FBLN5, ACTA1, TMEM119, VWF, GFAP, RARRES2, and FGA. These protein expression profiles across diverse patients are illustrated in Figure 2C. Notably, in non-responders, genes exhibiting higher expression levels tended to correlate with poorer prognoses in GC patients (Supplementary Figure S3). Importantly, upon incorporating validation data from an external experiment following a similar methodology, we observed a strong correlation and concordance in our findings (Figure 2D).

To gain a deeper understanding of the pathways associated with these DEPs, we conducted enrichment analyses using Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Gene Set Enrichment Analysis (GSEA). The outcomes of the GO enrichment analysis unveiled that the DEPs are primarily implicated in processes such as “complement activation,” “humoral immune response,” “regulation of blood coagulation,” “collagen-containing extracellular matrix,” and “blood microparticle” pathways (Figure 3A; Supplementary Table S6). Corresponding with the findings from both KEGG (Figure 3B; Supplementary Table S7) and GSEA (Figure 3C; Supplementary Table S8), the complement and coagulation cascades pathway (CCCP) emerged as notably significant. As shown in Figure 3D, nearly all genes within the CCCP exhibited upregulation in non-responders. A detailed exploration of the inter-regulatory relationships among proteins within the CCCP is given in Supplementary Figure S4.

The hallmarks of the response and resistance to the immune checkpoint blockade were amplified by a cross-talk between tumor cells and stromal and immune cells within the TME (Morad et al., 2021). Therefore, we observed an association between the ICI response and the TME. To investigate this relationship, we employed single-sample Gene Set Enrichment Analysis (ssGSEA) to assess the differential infiltration levels of various immune cell types (He et al., 2018). We observed significant differences in the abundance of certain tumor-infiltrating immune cells, including activated B cell, activated CD8 T cell, central memory CD8 T cell, memory B cell, natural killer T cell, and plasmacytoid dendritic cell, between the two groups. Notably, activated CD8 T cells showed consistency in the validation data (Figure 4; Supplementary Table S9).

Although some genes with significant changes also exhibit discriminative patterns (Supplementary Figure S5), individual indicators tend to lack robustness. To improve the biological comprehensibility of our predictive model, we evaluated the potential of utilizing pathway enrichment scores and levels of immune cell infiltration as predictors of patient response to immunotherapy. To provide clarity, we plotted the receiver operating characteristic (ROC) curves for the CCCP and activated CD8 T-cell infiltration. With the CCCP as the predictor, the AUC values for our data and validation data were 0.82 and 0.78, respectively (Figures 5A, B). With the level of activated CD8 T-cell infiltration as the predictor, the AUC values for our data and validation data were 0.65 and 0.92, respectively (Figures 5C, D).

While both predictors show some potential for the response to immunotherapy, we further sought to harness the power of machine learning algorithms to achieve superior predictive performance. Specifically, we implemented support vector machines (SVMs) and optimized the performance of our predictive model (Supplementary Table S10). Table 2 shows the performance achieved by our predictive model on the different datasets. The ROC curves show that our model achieved AUC values of 0.97 and 0.96 in our data and the validation data, respectively, indicating excellent prediction accuracy and robustness across datasets (Figures 5E, F).

This study was a retrospective study involving 28 enrolled patients with advanced GC between June 2019 and May 2021 at the First Affiliated Hospital of Zhengzhou University, Zhengzhou, Henan, China. The study obtained the ethical approval of our hospital with number 2020-KY-386. All cases were reexamined independently by senior pathologists, and histological diagnosis was performed based on the WHO classification of central nervous system tumors. The clinical data, including age, sex, treatment, body mass index, and pathology results, were obtained from the medical records of the enrolled patients. All patients and their families provided informed consent. The clinical information is shown in Supplementary Table S1. Patients use camrelizumab (AiRuiKa™), a PD-1 inhibitor being developed by Jiangsu Hengrui Medicine Co., Ltd. (Markham and Keam, 2019), according to the protocol. Responders were defined as patients with a RECIST complete response (CR) or partial response (PR), while non-responders were defined as those with progressive disease (PD) or stable disease (SD).

Tumor tissues were obtained any time before initiation of study treatment. The freshly acquired tissues were first rinsed with physiological saline and subsequently fixed in a 10% formalin solution, followed by embedding in paraffin wax.

After protein extraction, trypsin digestion, MS analysis, and database search, we obtained the raw data of proteomics. Supplementary Information provides more detailed processing information (Supplementary Method).

Following the database search using MaxQuant, 5,107 proteins were identified, of which 4,061 were quantified (Supplementary Table S2). On average, 3,061 proteins per sample were quantified (Supplementary Figure S1A; Supplementary Table S3). Subcellular distribution analysis conducted through the Hum-mPLoc 3.0 database (Zhou et al., 2017) revealed that the majority of identified proteins were cytoplasmic, followed by nuclear and plasma membrane proteins. This distribution aligns with existing references and likely reflects the actual distribution within the tissue, without any subcellular bias (Supplementary Figure S1B). To control the variation among samples, 28 samples were assessed using Pearson’s correlation. The results displayed high correlations among samples with an average Pearson’s correlation coefficient of 0.827 (Supplementary Figure S1C). To identify samples with abnormal protein abundance distributions, the log2-transformed protein intensities were analyzed. The violin plot demonstrates that each sample exhibited a similar distribution when log2-transformed (Supplementary Figure S1D). Quantitative proteomics experiments based on MS frequently generate data with missing values, which can profoundly affect downstream analyses. To address this issue, we evaluated several methods using NAguideR (Kang et al., 2017) and selected the “SeqKNN” method for imputing missing values (Supplementary Figure S2). This imputation process resulted in a dataset consisting of 2,884 proteins for further downstream analysis across all samples (Supplementary Table S4).

The validation cohort (https://www.iprox.cn//page/subproject.html?id=IPX0004819001) is a group of GC patients treated with anti-PD1 therapy, including seven responder cases and seven non-responder cases (Shi et al., 2023). The sample is also FFPE, and the proteins were quantified using an MS-based label-free method. To ensure data consistency, we downloaded the raw data and re-performed database search and quantification using MaxQuant with the same parameters. After database searching, 4,263 proteins were quantified. After missing value imputation, 2,472 proteins were used for further downstream analysis.

All statistical analyses were performed using R software. The code for proteomics-based analysis is available at https://github.com/longfei8533/Predicting-Response-to-PD-1-Inhibitor. The relationships between the treatment response and clinical–pathological features were evaluated using Wilcoxon’s test for continuous variables and Fisher’s exact test for categorical variables. Pearson’s correlation coefficient was used in the correlation analysis. We used two-sided Student’s t-test to evaluate the difference in the log2-transformed protein intensity between groups. The R “clusterProfiler” package (Cell, 2023) was used for GO and KEGG functional enrichment analyses and GSEA.

We selected immune-related signatures representing 28 different types of immune cells from the study by Charoentong et al. (2017). To estimate the relative abundance of immune cell infiltration in the tumor based on protein intensity, we utilized the Gene Set Variation Analysis (GSVA) package (Hänzelmann et al., 2013) with the ssGSEA method (Barbie et al., 2009). To compare the differences between responders and non-responders, we conducted unpaired one-sided Student’s t-tests.

For the CCCP predictor, we selected the “KEGG_COMPLEMENT_AND_COAGULATION_CASCADES” gene set from MSigDB (Liberzon et al., 2011). For the activated CD8 T-cell predictor, the signature gene set from the study by Charoentong was used. Both prediction methods through the utilization of ssGSEA were used to infer the relative activity of pathway.

Considering the small-sample size issue, we utilized SVM as a machine learning approach from the “e1071” package to predict the patient response to medication. Initially, we aligned our dataset with the validation set by selecting overlapping analyzable proteins and then normalized the data using these proteins. Subsequently, important features were chosen using Boruta (Kursa and Rudnicki, 2010), resulting in the retention of 10 proteins as features. We conducted model training after adjusting the SVM parameters; please refer to the Supplementary Material for the specific parameter values. The classifiers developed from the training cohort were then applied to the validation cohort, and the AUC of each classifier was calculated.