The next-generation sequencing data and clinicopathological information for 407 patients with gastric cancer were downloaded from The Cancer Genome Atlas Stomach Adenocarcinoma (TCGA-STAD) data collection (https://portal.gdc.cancer.gov). The microarray data for 60 patients with gastric cancer in GSE30727 were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/geo2r/?acc=GSE30727). We also used a dataset with immunotherapy information [the R package “IMvigor210CoreBiologies (version 1.0.0)”] (Balar et al., 2017) to evaluate the prediction efficiency of our scoring system. Drug sensitivity data between different gastric cancer cells in vitro were downloaded from Genomics of Drug Sensitivity in Cancer (GDSC, http://www.cancerrxgene.org/downloads).

From the Human DNA Repair Genes website (https://www.mdanderson.org/documents/Labs/Wood-Laboratory/human-dna-repair-genes.html) (Wood et al., 2005) which was last modified by Drs. R. Wood and M. Lowery on 10th June 2020, a list of 219 DNA repair genes were defined. We used the DEseq2 (version 1.28.1) package in R (version 4.0.5) to analyze the absolute counts for differential gene expression in TCGA-STAD and employed the online platform GEO2R to assess whether these DNA repair genes are differentially expressed between tumor and normal tissues (adjusted p < 0.05, |log2FC| > 1) (Supplementary Table S1). Fifteen differentially expressed DNA repair genes were commonly identified in TCGA-STAD and GSE30727, which are: PRKDC, FANCI, LIG1, RECQL4, FANCA, BRCA2, FANCD2, UBE2T, POLQ, PARPBP, EXO1, XRCC2, RAD54L, EME1 and FANCB. These genes constituted our DNA repair gene signature (DRGS). Supplementary Table S2 showed expression level of the 15 genes in the DRGS.

Before the analysis, we applied the Z-score normalization to the expression level of the 15 genes from each patient sample. The principal component analysis (PCA) is a proven technique, which can reduce dimension, improve interpretability and minimize information loss of large datasets at the same time. We employed the PCA as a dimensionality reduction method to obtain the 1st dimension correlation coefficients of the 15 genes (see Supplementary Figure S1 and Supplementary Table S3).

The formula for the DRGS scoring system was as follows: DRGS  score= ∑ ​ ( correlation    coefficients  for  gene  i  ) ∗  ( normalized  gene  i  expression  level  ) (1)

To determine the relationship between the DRGS scoring system and immunocyte infiltration in the gastric cancer tissue, we performed immunocyte infiltration analysis through the online platform Immune Cell Abundance Identifier (ImmuCellAI) (http://bioinfo.life.hust.edu.cn/ImmuCellAI#!/) using TCGA-STAD samples (n = 375) (Miao et al., 2020). Then, we calculated DRGS scores for the same TCGA-STAD samples and used the median score as the cut-point to divide these samples into high-score subtype and low-score subtype (Supplementary Table S4). The Mann Whitney U test was used to compare the results of multiple immune cell infiltration between these two subtype groups. We applied hierarchical clustering to immune cells with significantly different infiltration scores between the two groups. The Mann Whitney U test was used to compare the infiltration score between high-score subtype and low-score subtype. We also used ImmuCellAI to predict the patient response to ICB therapy (Supplementary Table S5).

Samples from the two groups mentioned above (TCGA-STAD, High score n = 187, Low score n = 188) were also used to explore the relationship between the DRGS score and various tumor characteristics. We compared the Lauren classification, molecular subtype and MSI status using the chi-square test (Supplementary Table S6). Tumor mutation levels were calculated through R package “Maftools” (Mayakonda et al., 2018). We used the Mann Whitney U test to compare TMB and PD-L1 expression between the two groups. In light that the TMB and PD-L1 expression are strongly related to the ICI therapy, we performed the Spearman’s correlation to analyze the relationship between the DRGS score and the two immunotherapeutic biomarkers (Supplementary Table S7).

The public dataset IMvigor210CoreBiologies (version 1.0.0) contains a total of 348 advanced urothelial bladder cancer patients who received atezolizumab treatment. We applied the DRGS score to the dataset to investigate its predictive ability of immunotherapy. We performed normalization and PCA analysis with the same method as mentioned above, and all the 348 patients were divided into two equal groups according to their DRGS score (High score n = 174, Low score n = 174). We then compared their response to immunotherapy and the immunotypes of tumor by the Chi-square test (Supplementary Table S8). Also, the Kaplan-Meier analysis of OS in the two groups was performed, and log-rank test was used to compare survival curves.

We also analyzed the prognostic ability of the DRGS scoring system in TCGA-STAD samples, through Kaplan-Meier analysis of OS and the log-rank test for comparison. The independence of the DRGS score system was verified by univariate and multivariate Cox regression. Three independent gastric cancer datasets, GSE26901, GSE15459, and GSE26899, which were downloaded from the GEO database as validation sets and analyzed with the same method described above.

We downloaded the drug sensitivity data from GDSC combined with the RNA-seq data from 23 types of gastric cancer cells. We performed the Pearson correlation analysis to calculate the correlation between the drug sensitivity level of different chemotherapeutic agents and the DRGS score. The criteria for significant correlation were p < 0.05 and correlation > 0.4.

After obtaining the consent of 10 gastric cancer patients, tissue sections were obtained from the Department of Pathology, the Affiliated Drum Tower Hospital of Nanjing University Medical School. After blocking with endogenous peroxide and protein, the sections were then incubated with diluted specific PD-L1 and DRGS protein antibodies at 4°C overnight, respectively. The next day, the sections were incubated with a secondary antibody at 37°C for 1 h. The sections were then stained with 3,3-diaminobenzidine (DAB) solution for 3 min and counterstained with hematoxylin, and photographed under a microscope. Two experienced pathologists scored these samples according to the percentage of staining intensity cells and staining intensity score (1 = low positive; 2 = positive; 3 = high positive). The formula for the IHC score was as follows: IHC  score= ∑ ​ ( percentage  of  staining  intensity  cells  ) ∗  ( staining  intensity  score  ) . (2)

Data from our IHC results and the Human Protein Atlas (HPA) (https://www.proteinatlas.org/) were analyzed to verify the protein expression of DRGS and PD-L1 in tumor tissues and normal tissues, and to determine whether the expression differences and the correlation between DRGS and PD-L1 were consistent with the previous mRNA results from TCGA. (We performed the Pearson correlation analysis to calculate the correlation between protein expression of DRGS and PD-L1, and the difference was considered significant if the p-value is < 0.05.)

Statistical methods were described in different sections above. The significance level of these statistical tests was two-sided p value < 0.05. The software used for analysis is R (version 4.0.5).

219 DNA repair genes were obtained from the Human DNA Repair Genes website as described above. Next, we used DEseq2 package in R to process the raw data from TCGA-STAD and used the online platform GEO2R to find genes differentially expressed between tumor (n = 375) and normal (n = 32) tissues (adjusted p value < 0.05, |log2FC| > 1). Finally, we identified 15 common genes among 219 DNA repair gene, which were differentially expressed in both TCGA-STAD and GSE30727 (Tumor n = 30, Normal n = 30) (Figure 1A, Table 1). We then used PCA analysis as a dimensionality reduction method to obtain the first dimensional correlation coefficients of the 15 genes (Figure 1B). Based on these coefficients, we created a 15-DNA repair gene signature (DRGS) scoring system (for details, please refer to the methods section). The DRGS score was then used to divide patients into two groups (high and low score groups).

We identified the immune characteristics of 375 TCGA-STAD samples through ImmuCellAI to estimate the specific fractions of 24 types of immune cells in each gastric cancer sample (Supplementary Table S4). To analyze the composition of immune cells in different score subgroups (High score n = 187, Low score n = 188), we used the Mann Whitney U test to compare the distribution of immune cell types in the two score groups. The result showed significant differences in the infiltration rates of 14 immune cells between the two groups. We found that anti-tumor cells such as γδ T cells, neutrophil cells and Th1 cells were more abundant in the DRGS high score group (p < 0.001), while immunosuppressive or tumor enhancement cells such as B cells (p < 0.001), CD4 T cells (p < 0.001), monocytes (p < 0.001), Th17 cells (p = 0.031), Tfh cells (p = 0.002) and Tr1 cells (p = 0.012) were more abundant in the DRGS low score group (Figure 2, Figure 3A). In our study, infiltration score of samples in the DRGS low score group was significantly higher than those in the high score group (p < 0.001, Mann-Whitney U test, Figure 3B). Furthermore, we used ImmuCellAI to predict the response of ICB therapy (Supplementary Table S5) and found that there would be 53.5% of patients in the DRGS high score group and fewer patients (37.8%) in the DRGS low score group, in response to ICB therapy (p = 0.003, Chi-square test, Figure 3C).

To examine the relationship between our score grouping and other histologic and molecular classifications, we analyzed clinical data of 372 samples of gastric cancer patients from the TCGA-STAD database (Supplementary Table S6). According to the Lauren classification, gastric cancers can be divided into 3 subtypes, i.e., intestinal, diffuse and mixed. In our study, the DRGS high score group (n = 115) was characterized by a very high proportion of intestinal subtype samples (86.1%), while diffuse and mixed subtypes samples were only 7.8 and 6.1%, respectively. On the other hand, the DRGS low score group (n = 115) is comprised of 51.3% intestinal, 40.0% diffuse and 8.7% mixed samples. Statistically, there were significantly more intestinal and fewer diffuse type gastric cancer samples in the DRGS high score group as compared to the low score group (p < 0.001, Chi-square test, Figure 4A).

TCGA proposed 4 subtypes of gastric adenocarcinoma: 1) EBV-positive, 2) microsatellite instability (MSI), 3) genomically stable (GS), and 4) chromosomal instability (CIN) (The Cancer Genome Atlas Research Network, 2014). In our study, the DRGS high score group (n = 169) had a higher percentage of EBV-positive (11.2%) and MSI (27.8%) subtype samples. In contrast, the DRGS low score group (n = 169) had a higher percentage of GS (23.1%) and CIN (63.9%) subtype samples (p < 0.001, Chi-square test, Figure 4B).

On the basis of the frequency of mutations in microsatellite markers, gastric cancer can be classified as MSS, MSI-H, and MSI-L (Rhyu et al., 1994; Yoon et al., 2013). In our study, the percentage of MSI-H subtype patients (31.7%) was higher in the DRGS high score group (n = 120) than that in the low score group (7.5%) (n = 120). In contrast, there were more MSS subtype patients (72.5%) in the DRGS low score group than in the high score group (55.8%). (p < 0.001, Chi-square test, Figure 4C).

To investigate whether the prognostic impact of the DRGS score is independent of clinical factors that could be associated with patient outcomes, we first performed univariate Cox regression analysis on all available clinical parameters in the TCGA-STAD dataset (n = 372). Clinical parameters that demonstrated significant prognostic impact (p < 0.05 by the Wald test) were selected for multivariate Cox regression analysis along with the DRGS. The result indicated that the DRGS is an independent prognostic factor (p = 0.038, Table 2). Taking the median score as the cut-off value, patients with DRGS high score had a better OS than those with DRGS low score (p = 0.007, log-rank test, Figure 5A). Then, the prognostic value of scoring system was validated using the GSE26901 (n = 109), GSE15459 (n = 192) and GSE26899 (n = 93) gastric cancer datasets. As shown in Figure 5, the patients in GSE26901 in the DRGS high score group had a significantly better prognosis than those in the low score group (p = 0.011, log-rank test, Figure 5B), consistent with the results from the TCGA dataset. However, there was no significant OS difference in GSE15459 (p = 0.77, Figure 5C) or GSE26899 (p = 0.25, Figure 5D) between the two DRGS groups.

To identify the relationship between TMB and the DRGS score, we analyzed the somatic mutation data of 365 gastric adenocarcinoma samples from TCGA-STAD (Supplementary Table S7). The result indicated that the DRGS score was positively associated with TMB (rho = 0.50, p < 0.001, Figure 6A). TMB of samples in DRGS high score group (n = 183) was significantly higher than that in DRGS low score group (n = 182) (p < 0.001, Mann-Whitney U test, Figure 6B). All mutation events were divided into different categories, among which missense mutations accounted for the largest proportion in the classification of variation (Supplementary Figure S2A, B in the Supplementary Material), single nucleotide polymorphisms (SNPs) occurred much more frequently than insertion or deletion, and C > T was the most general of single nucleotide variants (SNV) in all TCGA-STAD samples. Notably, median variant per sample of DRGS high score group (121) was higher than that of DRGS low score group (59.5). Furthermore, mutation events of the top 10 mutated genes in each sample were shown in the waterfall plot. Amongst DRGS high score group, the mutation frequency of TTN was the highest (59%), followed by TP53 (54%) and MUC16 (39%, Supplementary Figure S2B). Amongst DRGS low score group, the highest mutation frequency gene was TP53 (34%), followed by TTN (33%) and MUC16 (21%, Supplementary Figure S2A).

PD-L1 expression is another important biomarker of susceptibility to PD-L1 blockade. We found that PD-L1 expression in DRGS high score group (n = 183) was significantly higher than in DRGS low score group (n = 182) (p = 0.002, Figure 6D, Supplementary Table S7). The DRGS score is positively correlated to PD-L1 expression (rho = 0.17, p < 0.001, Figure 6C).

In view of the association between our score system and the immune microenvironment of the tumor, as described above, we tested the ability of the DRGS score system to predict the response of patients to immunotherapy. This analysis was based on the IMvigor210CoreBiologies cohort (n = 348), a large phase 2 trial investigating the clinical activity of PD-L1 blockade with atezolizumab in metastatic urothelial carcinoma (mUC) (Gibney et al., 2016) (Supplementary Table S8). We found that DRGS high score patients (n = 174) exhibited prominent prolonged overall survival (p = 0.021, Figure 7A). In addition, 298 patients in this cohort showed varying degrees of response to anti-PD-L1 blockers. DRGS high score patients (n = 149) comprised 14.8% complete response (CR), 18.8% partial response (PR), 18.8% stable disease (SD) and 47.6% progressive disease (PD), while DRGS low score patients (n = 149) comprised 2.0% CR, 10.1% PR, 23.5% SD and 64.4% PD. The Chi-squared test performed between DRGS high and low score groups also showed significantly better therapeutic outcomes in DRGS high score patients than the low score patients (p < 0.001, Figure 7B).

We also analyzed the three immune subtypes (immune inflamed, immune excluded and immune desert) of IMvigor210CoreBiologies in two DRGS group patients. DRGS high score patients (n = 142) comprised 33.1% inflamed, 45.8% excluded and 21.1% desert, while DRGS low score patients (n = 142) comprised 19.0% inflamed, 48.6% excluded and 32.4% desert. Therefore, the DRGS high score patients had a significantly higher percentage of “immune inflamed” tumors and a lower percentage of “immune desert” and “immune excluded” tumors compared to the low score patients (p = 0.012, Chi-square test, Figure 7C).

To explore the effect of the DRGS score system on drug response in chemotherapy, we evaluated the association between the DRGS score and the response to drugs in gastric cancer cell lines (n = 23) in the Genomics of Drug Sensitivity in Cancer (GDSC) database. We identified 5 significantly correlated pairs between the DRGS score and drug sensitivity, all of which demonstrated drug resistance correlated with the DRGS score, including insulin-like growth factor receptor IGF1R_3801 (Rs = −0.53, p = 0.019), AKT inhibitor ipatasertib (Rs = −0.46, p = 0.028), surviving inhibitor sepantronium bromide (Rs = −0.54, p = 0.007), ULK1 protein kinase ULK1_4989 (Rs = −0.47, p = 0.043) and AKT inhibitor uprosertib (Rs = −0.50, p = 0.021) (Figure 8A). In addition, we found that these drugs mostly target PI3K/MTOR, IGF1R and apoptosis regulation signaling pathways (Figure 8B).

Through IHC analysis, we found in gastric tumor tissues (n = 10) that the protein expression of most of the genes in the DRGS, except for PARPBP, EME1, and RAD54L, was significantly increased (Figure 9A). In the HPA data, compared with normal tissues, the expression of DRGS except FANCA, UBE2T, POLQ, EXO1, and XRCC2 (They were not included) was remarkably upregulated (Figure 9B).

We also analysed the relationship between the DRGS and PD-L1 protein expression. As shown in Figures 10A,B, based on IHC analysis (n = 10), the protein expression of the genes in the DRGS, except for PARPBP, EME1, and RAD54L, was positively correlated to PD-L1 expression (rho>0.7, p < 0.05, Figures 10A,B). In the HPA data, the protein level of PRKDC, FANCI, LIG1, RECQL4, FANCA, PARPBP, RAD54L, and FANCD2 (The rest of the 15 genes were not included) was positively related to that of PD-L1 (Figure 10C). Thus, the above results verified the results obtained from bioinformatics analysis.