The expression profiling and clinical traits of 375 GC samples were downloaded from the TCGA website (https://portal.gdc.cancer.gov/). Another GC-independent dataset (GSE84437) containing transcriptome data and survival information was obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). All raw data were pre-processed with the limma package in R software (4.1.3, R Foundation for Statistical Computing, Vienna, Austria). The exclusion criteria were set as follows: (1) histologic diagnosis was not GC; (2) samples without completed data for analysis; and (3) overall survival (OS) time of less than 30 days. We screened the list of UGRs using a relevance score greater than 10.0 from the GENECARDS database website (https://www.genecards.org/). Detailed information on these URGs is shown in Supplementary Table 1.

Firstly, we selected the TCGA dataset as a training set for model establishment. All UGRs were analyzed by univariate Cox regression to screen out prognosis-related URGs (p 0.05). Multivariate analysis was utilized to generate a prognostic model based on the candidate genes from univariate Cox regression. Each GC sample’s risk value was generated using the following formula: ${\sum }_{i=1}^{n}coe{f}_i*expressionlevelofUA{G}_i$. The coef represents the coefficient of each UAG calculated by multivariate Cox analyses. The GC cases were divided into high- and low-risk groups according to the median risk value.

We applied unsupervised clustering through the ‘ConsensusClusterPlus’ package [13]. The optimal cluster value was determined based on the similarity between the expression matrix of UAGs and the fuzzy clustering measurement ratio.

We obtained scRNA data (GSE167297) of GC from the GEO database. The scRNA-seq matrix for individual samples was analyzed using the Read10$\times$ algorithm of the ‘Seurat’ package [14]. Low-quality cells with $\le$300 detected genes or $\ge$10% mitochondrial genes were removed and normalized by ‘NormalizeData’, the top 1500 highly variable genes were then identified by the ‘FindVariableFeatures’ function. Data quality control and normalization were performed by ‘NormalizedData’. The t-SNE method was used to reduce data dimensionality and obtain different cell clusters. All cell populations were annotated using ‘FindAllMarkers’ based on the cell surface markers.

GSEA was conducted using GSEA software (version 4.0.1, Massachusetts Institute of Technology, and Regents of the University of California, CA, USA) with 1000 rows. Hallmark and KEGG terms were selected as the target gene sets from MSigDB. Total GC cases were divided into two groups based on risk and enriched for the above gene sets. Any gene set with FDR 0.25 and p 0.05 was regarded as significant.

CIBERSORT (https://cibersort.stanford.edu/) is an algorithm that estimates the abundance of cell types infiltrated by patients with tumors based on gene expression data [15]. Moreover, single sample gene set enrichment analysis (ssGSEA) was employed to immune activity and function between two risk groups. The drug sensitivity analysis was conducted by the ‘pRRophetic’ package to detect which drugs presented different sensitivities among the two risk groups. The half-maximal inhibitory concentration (IC50) was calculated to evaluate the drug sensitivity [16].

Two human GC cell lines (SGC7901 and BGC823) and one normal gastric epithelium cell line (GES-1) were purchased from the Cell Center of Shanghai Institutes. All cell lines were cultured in RPMI 1640 medium (Invitrogen) containing 10% fetal bovine serum and 1% antibiotics (100 U/mL penicillin G and 100 mg/mL streptomycin) at 37 °C with 5% CO${}_2$. The silencing RNA against OTULIN (si-OTULIN) was synthesized and purchased from RIBBIO (Guangzhou, China). The sequence of si-OTULIN is shown in Supplementary Table 2. Lipofectamine 3000 (Invitrogen) was used for cell transfection. Cell lines were authenticated by the Shanghai Biowing Applied Biotechnology Co.,Ltd and Procell Life Science & Technology Co.,Ltd. using Short Tandem Repeat profiling (PCR Amplification Kit, Corning, New York, USA). Mycoplasma test was performed in all cell lines every other week using the Mycoplasma PCR Detection Kit (Beyotime, Shanghai, China). Only mycoplasma-free cell lines were used.

Total RNA was extracted from cells using RNA easy reagent (Vazyme, Nanjing, China) based on the manufacturer’s instructions. RNA was then synthesized to cDNA by Primerscript Mix (Takara Bio, Kusatsu, Shiga, Japan). RT-PCR was applied on a Steponeplus system with SYBR Green Regent (Takara Bio, Japan). The relative expression was calculated using the 2${}^{-\mathrm{\Delta }\mathrm{\Delta }\text{CT}}$ method. The primer sequences of the genes are shown in Supplementary Table 2.

A total of 10,000 cells per well were seeded into a 96-well plate. After 24 hours of incubation at 37 °C, cell viability was determined by detecting ATP levels in cell lysates using the CellTiter-Glo (CTG) Assay.

A total of 200 cells/well were seeded into six-well plates and incubated for 1-week. Clones were fixed with methanol and stained with 0.1% crystal violet (Servicebio, Wuhan, China) when visible to the naked eye. Image-J (version 1.5b, LOCI, University of Wisconsin, Madison, WI, USA) was used for colony counting.

EdU assay was performed to evaluate cell proliferation. The cells were cultured in 96-well plates (5000 cells per well) with complete medium for 24 h. After a 2 h incubation with EdU solution (50 µM), all cells were fixed and stained with 1$\times$ Apollo® solution for 30 min. Finally, nuclei were stained with 1 $\times$ Hoechst33342.

All bioinformatics data were analyzed by R software (version 4.0, Lucent Technologies, Murray Hill, NJ, USA). and wet lab experimental data were analyzed using GraphPad 9.4 (GraphPad Software, Inc., San Diego, CA, USA). Survival curves were used to compare survival differences between the two groups using the Kaplan-Meier (KM) method. The discrepancies in clinical outcomes between groups were examined using the Kaplan-Meier survival analysis. Receiver operator characteristic (ROC) curves were plotted to confirm the accuracy of UAS. Cox regression was utilized to detect the independent performance of UAS in GC.

A flowchart of the present research is shown in Fig. 1. Firstly, univariate Cox analysis was employed to obtain eight potential URGs with significant prognosis values, which were then enrolled into the multivariate Cox method (Fig. 2A). Using stepwise Cox methods, we set up the UAS including five URGs (OTULIN, UBE2C, USP1, USP2, and MAPT) (Fig. 2B). The risk model equation was: (0.2508 $\times$ OTULIN) + (0.2446 $\times$ UBE2C) + (0.87167 $\times$ USP1) + (–0.3729 $\times$ USP2) + (–0.3461 $\times$ MAPT). Fig. 2C revealed that these URGs were significantly differentially expressed in GC and normal tissues using the GEPIA2 portal. It was shown that OTULIN, USP1, and UBE2C were overexpressed in GC specimens, whereas USP2 and MAPT were upregulated in normal tissues.

Fig. 1.

A flowchart of the present research. Gastric cancer (GC), The Cancer Genome Atlas Program (TCGA), Gene Expression Omnibus (GEO), microsatellite instability (MSI).

Fig. 2.

Construction process of the ubiquitination-associated signature (UAS). (A) Univariate Cox analysis. (B) Multivariate Cox analysis for screening model genes. (C) The expression patterns of UAS genes in gastric cancer (GC) and normal specimens. *p 0.05.

Based on model risk scores constructed from five URGs, we divided GC cases into high- and low-risk groups, showing good predictive ability of the model in both the training and validation groups, and defined them as a USA-high cohort and UAS-low cohort. Survival curves demonstrated that high-risk score was associated with dismal OS in the training set (TCGA cohort), which was confirmed by the verification set (GSE84437) (Fig. 3A,B). The 1-, 3- and 5-year area under the curves (AUCs) were 0.611, 0.705, and 0.785 for the training set (Fig. 3C), 0.668, 0.619, and 0.723 for the verification set (Fig. 3D), respectively, indicating the beneficial predictive ability of our constructed UAS. Fig. 3E,F illustrates the relationship between risk scores and the pattern of patient survival status (alive or dead) distribution. Mortality was significantly decreased in the UAS-low cohort than in the USA-high cohort; both were independent datasets (Fig. 3G,H).

Fig. 3.

Performance of the ubiquitination-associated signature. (A,B) Survival curves of survival outcome of GC cases between two groups in the training and verification sets. (C,D) Receiver operating characteristic (ROC) curves to predict the performance of UAS in the two GC cohorts. (E,F) Demonstration of the UAS based on the risk score of the training and verification sets, respectively. (G,H) The proportion of deaths and survival in two subgroups. ROC, receiver operator characteristic.

We employed Cox regression to detect the independent prognostic score of UAS in GC patients. Fig. 4A revealed that univariate analysis remarkably associated the UAS-high score with dismal survival outcomes. Multivariate Cox indicated that the UAS score remained meaningful for forecasting prognosis. In addition, the prognostic independence of UAS was confirmed in the test cohort (Fig. 4C,D). Additionally, we observed that UAS could forecast survival outcomes of different clinical subgroups in terms of age ($\le$65, n = 197 samples, $>$65, n = 241 samples), gender (male, n = 285 samples, female, n = 158 samples), grade (grade 1–2, n = 171 samples, grade 3, n = 263 samples) and stage (stage I–II, n = 189 samples, stage III–IV, n = 227 samples) (Fig. 4E).

Fig. 4.

Independent prognosis analysis. (A,C) Univariate Cox analysis. (B,D) Multivariate Cox analysis. (E) Subgroup survival analysis of the UAS in terms of age, gender, grade, and stage.

To develop a novel molecular subtype in GC, consensus clustering analysis was employed based on 268 URGs. According to the cumulative distribution function (CDF), we conducted the tracking plot from k = 2 to 9, where k = 2 was identified as the cluster number to divide patients into cluster A (n = 304) and cluster B (n = 504) (Fig. 5A). Principal component analysis (PCA) indicated that two clusters of GC cases could be effectively distinguished by our proposed novel molecular subtype (Supplementary Fig. 1A). Cluster B had a better outcome than cluster A. There was an obvious difference discrepancy in immunocyte infiltration between the two subgroups (Fig. 5B,C). Supplementary Fig. 1B summarizes the relationship between UAS-based molecular clusters and different clinical traits. We then identified potential pathways (ubiquitination, DNA repair, and cell cycle) in which clusters might be involved by KEGG enrichment analysis (Fig. 5D).

Fig. 5.

Cluster analysis for determination of URGs-based molecular subtype. (A) Consensus clustering results. (B) Survival analysis and (C) immune infiltration analysis between two clusters. (D) Gene set variation analysis (GSVA) based on Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment. *p 0.05; **p 0.01; ***p 0.001. URGs, ubiquitination-related genes.

Next, we applied CIBERSORT and ssGSEA analyses to explore the ability of the model to predict the immune microenvironment of GC patients. In Fig. 6A, B cells memory, B cells naive, monocytes, mast cells resting, and T cells CD4 memory resting positively correlated with risk score, whereas macrophages M0, T cells CD4 memory activated, and T cells follicular helper were negatively correlated with risk score. The results from ssGSEA revealed that immune functions such as CCR, major histocompatibility complex (MHC) class I, and type II interferon (IFN) response differed significantly among the two groups (Fig. 6B). Fig. 6C shows the expression patterns of classical immune checkpoints among two groups. It was shown that CD86, CD200, CD40, NRP1, CD48, CD27, and TNFSF15 were upregulated in the UAS-high group, whereas in the UAS-high group, other immune checkpoints exhibited a decreased expression level.

Fig. 6.

Immune landscape analysis. (A) The correlation between risk score and immune cells. (B) The difference in immune functions in two risk groups. (C) Immune checkpoints analysis of the UAS. *p 0.05; **p 0.01; ***p 0.001, ns, no significance.

Considering the crucial regulatory role of m6A regulators in GC progression, we determined the relationship between risk score and the expression of m6A regulators. We found that ZC3H13 and FTO were highly expressed in the UAS-high group. In contrast, HNRNPC, ALKBH5, YTHDF2, RBM15, YTHDF1, and METTL3 were downregulated in the UAS-high group (Fig. 7A). There were significant differences in the expression of chemoradiotherapy sensitivity-related genes between the two groups. As shown in Fig. 7B, AKR1C1, EGFR, FLT3, and KIT were significantly overexpressed in the UAS-high group (Fig. 7B). The proportion of patients with MSS was significantly lower in the UAS-low cohort than in the UAS-high cohort. In comparison, the proportion of MSI-H was significantly higher in UAS-high cohort (Fig. 7C). In addition, the risk score was negatively correlated with RNA stemness score (RNAss) (Fig. 7D). Drug sensitivity analysis revealed that docetaxel (p = 1.7 $\times$ 10${}^{-8}$) had a lower IC50 value in the high-risk cohort, whereas the other three drugs displayed a higher IC50 value in the high-risk cohort (p = 7.1 $\times$ 10${}^{-5}$ for gefitinib, p = 4.3 $\times$ 10${}^{-9}$ for rapamycin, p = 4.4 $\times$ 10${}^{-9}$ for sorafenib) (Fig. 7E).

Fig. 7.

Clinical potency of the UAS. (A) The expression patterns of m6A regulators and (B) chemosensitivity-related genes between risk groups. (C) Correlation between MSI and risk score. (D) Cancer stemness analysis. (E) Drug sensitivity analysis (docetaxel, gefitinib, rapamycin, and sorafenib). *p 0.05; **p 0.01; ***p 0.001, ns, no significance.

Subsequently, we detected the expression pattern of five model genes at the single-cell level. A total of 14 GC samples from GSE167297 were enrolled into single-cell analysis (Supplementary Fig. 2A). The top 1500 variable genes with red labels were collected (Fig. 8A). Fig. 8B revealed a favorable integration effect between cluster all samples and 50 PC groups. PCA analysis showed the distribution of all GC cases with no significant batch effects (Supplementary Fig. 2B). Using t-SNE diagrams, a total of 30,365 cells were identified and further clustered into 20 different subpopulations (Fig. 8C). Based on the surface marker, we determined 8 cell population types including B cells, dendritic cells (DCs) endothelial cells, epithelial cells, monocytes, NK cells, smooth muscle cells and T cells (Fig. 8D). Next, the cellular location of five URGs was detected. We observed that OTULIN, USP1, and UBE2C were enriched in malignant epithelial cells, whereas USP2 and MAPT were downregulated in malignant epithelial cells (Fig. 8E).

Fig. 8.

Single-cell sequencing analysis. (A) A total of the top 1500 variable genes with red labels. (B) Data integration and dimensionality reduction. (C) All cells were clustered into 20 subclusters. (D) Classification of all cells into B cells, dendritic cells (DCs), endothelial cells, epithelial cells, monocytes, NK cells, smooth muscle cells, and T cells. (E) Cellular location of five URGs (OTULIN, USP1, USP2, UBE2C, and MAPT).

We selected OTULIN to verify our proposed model by various in vitro experiments. Fig. 9A shows that OTULIN, USP1, and UBE2C were highly expressed in GC cell lines, whereas USP2 and MAPT were highly expressed in the GES-1 cell line. A PCR assay was carried out to confirm the transfection efficiency of siRNA-OTULIN and OTULIN (Fig. 9B). CTG assay, colony formation, and EdU assay were applied to evaluate the cell viability in GC cells. Cell proliferation was significantly inhibited following the silencing of OTULIN in SGC7901 cells (Fig. 9C–E).

Fig. 9.

Identification of OUTLIN as a novel oncogenic player in GC. (A) Expression patterns of five model genes in different cell lines were detected by Real-Time Quantitative Reverse Transcription PCR (qRT-PCR) assay. (B) PCR assay indicates the favorable transfection efficiency in GC cell lines. Cell proliferation in different treatment groups was evaluated by (C) CTG assay, (D) colony formation, and (E) EdU assay. Scale bar = 200 µm. *p 0.05; **p 0.01; ***p 0.001. CTG, CellTiter-Glo.