In this study, we used FFPE surgical resection samples from 833 GC patients with TNM stage II/III from previous work (31). These samples were collected between 2004 and 2016 and came from five hospitals, namely, the Peking University Cancer Hospital & Institute (PKUCH, N = 387), the Fourth Medical Center of PLA General Hospital/304 Hospital (304H, N = 210), the Xijing Hospital of Digestive Diseases (XJH, N = 112), the Medical School of Chinese PLA/301 Hospital (301H, N = 71), and the Shanxi Cancer Hospital (SXCH, N = 53). Patients had provided written informed consent and had complete follow-up records with adequate clinical annotations. The median follow-up time was 3.7 years (a range of 0.08 to 10.4 years).

We used the data in three steps: The first step is to split the data into a discovery set (PKUCH and XJH) and an independent validation set (304H, 301H, and SXCH). The features extracted from the whole discovery set by AE were used for consensus clustering to obtain labels of survival-risk subtypes. The second step is to train the classifier model by dividing the discovery set into training and test sets at a ratio of 7:3. The third step uses data from an independent validation set to evaluate the prediction accuracy of the DL-based prognosis model.

We use proteomic data from 499 samples from the discovery set as the input for the AE framework to feature transformation. AE consists of an encoder and a decoder (42), which is a feed-forward and non-recursive neural network commonly used in semi-supervised and unsupervised learning (33, 43). Given an input layer with an input x = (x₁, x₂, …… x_n) of dimension n, the objective of an AE is to reconstruct x with the output x’ (x and x’ have the same dimension) via transforming x through successive hidden layers.

For the hidden network layer, we use Relu as the activation function between the input layer x and the output layer y. That is:

We use Sigmoid as the activation function for the reconstructed layer. That is:

However, the bottleneck layer does not use any activation functions.

The objective of AE training is to find the different weight vectors W_i, minimizing a specific objective function. We chose mean-square error (MSE) as the objective function, which measures the error between the input x and the output x’.

Here, n refers to the sample number, x_i and x_i’ refer to the input and output values of the current sample, respectively.

We constructed the AE with three hidden layers (500, 100, and 500 nodes, respectively) using Python’s Keras library (https://github.com/fchollet/keras). The bottleneck layer of the AE was used to generate novel features. Finally, the AE was trained with the Adam optimizer as the optimization function and 0.001 as the learning rate. A gradient descent algorithm with 80 epochs was used. Epoch here refers to the iteration of the learning algorithm on the whole training data set.

The AE reduced the original features to 100 new features obtained from the bottleneck layer. For these transformation features generated by AE, we use the R package ConsensusClusterPlus (44) to perform consensus clustering. We determined the optimal number of clusters with two metrics: (1) Silhouette index and (2) log-rank P-value. The clustering algorithm was k-means using Euclidean distance. The proportion of samples selected was 80% in each resampling, and the number of clusters considered was 2 to 5. Among them, a consensus matrix with k = 2 appeared to have the clearest cut between clusters and showed a significant association with the survival of the patients.

After obtaining the labels through consensus clustering, we built three supervised classification models, namely, random forest (RF), logistic regression (LR), and support vector machines (SVM). Using the univariate Cox proportional hazards (Cox-PH) model, we identified 56 prognosis-related proteins (log-rank P <0.01) used for training classifiers.

The expression values of these 56 proteins on the training and test sets first standardized with Z-Score before being entered into classifiers. The average areas under the curve (AUC) of models were evaluated on the training set combined with 10-fold cross-validation (CV), and then these models were further determined based on their performance on the test set.

We built three classifiers using Python’s scikit-learn package that could be used for performing grid search to find the best hyperparameters of the three models. Finally, a RF model containing 10 trees was determined.

To verify the advantages of features transformed from AE, we compared the performance of AE with traditional machine learning methods. Here, Principal Component Analysis (PCA), a traditional linear dimension reduction method, and Uniform Manifold Approximation and Projection (UMAP) (45), a new nonlinear dimension reduction algorithm, were used to find the optimal number of retained features, respectively. The samples were then clustered using the same consensus clustering procedure ( Figure 1 ).

A differential expression analysis was performed to identify the differentially expressed proteins between the two survival subtypes. The Mann-Whitney U test was used to identify proteins with significantly different expression between the two subtypes, and the Benjamini–Hochberg method was used to adjust the P-values.

We used the Reactome pathway database (https://reactome.org/) (46) to perform functional enrichment analysis for the differentially expressed proteins of the two subtypes.

A total of 833 FFPE GC samples from five hospitals were used in this study. There were 309 (37%) TNM stage II patients and 524 (63%) TNM stage III patients, 582 (70%) of whom received adjuvant chemotherapy.

The median of protein detection in the five hospitals was between 1,273 and 1,543 ( Supplementary Figure 1A ). Moreover, the PCA showed no clear boundaries in the five hospitals, suggesting that there was no batch effect caused by the source of samples ( Supplementary Figure 1B ). We combined two hospitals as discovery sets and randomly divided them into the training and test sets at a ratio of 7:3. The remaining three hospitals were combined as an independent validation set. As shown in Supplementary Figures 1C, D , the sample distribution of the three data sets is relatively balanced. Table 1 lists detailed clinicopathological information for patients in the training set (n = 349), the test set (n = 150), and the independent validation set n = 334.

We established a workflow based on the AE framework, with the structure highlighted in Figure 1A . We used the 100 nodes from the bottleneck-hidden layer of AE as new features, and then conducted consensus clustering for the hidden features. The silhouette index and log-rank P-value were used to judge the quality of clustering and to obtain the optimal number of clusters on the discovery set. Next, we used a univariate Cox-PH model to obtain prognosis-related proteins, which were used as features to establish a classifier. We trained the classifier with a 10-fold CV in the training set and tested it in the test set, then further validated it in an independent validation set. The workflow is shown in Figure 1B .

For the discovery set, we retained proteins detected in one-tenth (n = 49) of all samples, resulting in 2,058 proteins that were used for further analysis. These 2,058 proteins were transformed by the AE, and 100 nonlinear features were retained for consensus clustering. When K = 2, the consensus matrix exhibited the clearest cut among clusters ( Figure 2A ), with an average silhouette index of 0.97 ( Supplementary Figure 2A ).

The association between prognosis and chemotherapy within each subtype was analyzed. Significant differences in OS were found between the two subtypes: S-I with good prognosis (n = 217, 43%) and S-II with poor prognosis (n = 282, 57%) (log-rank P = 0.024, Figure 2B ). Additionally, we found that the 5-year OS rate of the S-I patients who received adjuvant chemotherapy was 65.3%, which is a significantly increased 12% compared with 52.6% for patients who received surgery only ( Figure 2C ). While no significant differences were observed in the 5-year OS rate for S-II between patients who received adjuvant chemotherapy (54%) and those who did not receive chemotherapy (51%) ( Figure 2D ).

We adopted the Mann–Whitney U test to perform differential expression analysis in two subtypes. Of the 884 differentially expressed proteins obtained (FDR <0.05 & fold change >2), 783 proteins were upregulated in S-I, and 101 proteins were upregulated in S-II ( Figure 2E and Supplementary Figure 2B ). Using the differentially expressed proteins above, we performed pathway enrichment analysis in the Reactome pathway database (46) to determine the pathways enriched in the two subtypes. S-I showed the characteristics of cell proliferation, mainly enriched in DNA replication, cell cycle, and programmed cell death. S-II showed the characteristics of the tumor microenvironment (TME), which was mainly enriched in the extracellular matrix (ECM)-related pathways ( Figure 2F ).

To evaluate the prognostic prediction accuracy of the DL-based workflow, we established a classifier using the two subtypes identified above as labels. Using a Cox-PH model, we obtained 56 prognosis-related proteins as features to train the classifier. The heat map of 56 proteins is shown in Supplementary Figure 3A . To build the classifier, we evaluated three commonly used machine learning (ML) models, namely, RF, LR, and SVM. The three models resulted in an average AUC of 0.92, 0.91, and 0.89 on the training set with a 10-fold CV, and 0.91, 0.93, and 0.92 on the test set, respectively ( Figures 3A, B and Supplementary Figures 3B, C ).

Based on the performance on the training and test sets, we observed that the three ML models performed equally well. To determine the robustness of the classifiers in predicting OS outcomes, we applied the three models to an independent validation set containing 334 patients. Subtypes of the validation set were predicted by the three models, followed by an association analysis between prognosis and chemotherapy within each subgroup. We discovered that there was a difference in OS between the two subtypes predicted by RF on the validation set (log-rank P = 0.026). Among the predicted S-I with good prognosis (n = 150, 45%), the 5-year OS rate of patients receiving adjuvant chemotherapy was increased by 25% compared to patients receiving surgery alone (70.5% vs 45.8%), consistent with the characteristics of the S-I in the discovery set. The predicted S-II with poor prognosis (n = 184, 55%) also showed similar characteristics as the S-II in the discovery set, and there was no significant difference in the 5-year OS rate between the chemotherapy group and the non-chemotherapy group (50% vs 46%, log-rank P = 0.026) ( Figures 3C, D and Supplementary Figure 3D ). In contrast, LR and SVM can only predict chemotherapy benefit, but not prognosis ( Figures 3E–H and Supplementary Figures 3E, F ). Collectively, the results show that the two subtypes were verified on the independent validation set through the RF classifier.

To verify the advantages of features transformed from AE in predicting prognosis and chemotherapy benefit, we compared them with two alternative approaches: PCA and UMAP (45).

In the first approach, we reserved the optimal 44 principal components for consensus clustering. Although the consensus matrix did not have a clear boundary ( Figure 4A ), this approach could detect survival subtypes with a significant log-rank p value (log-rank P = 0.045, Figure 4B ). Additionally, the two subtypes also exhibited similar characteristics in terms of prognosis and chemotherapy benefit with S-I (log-rank P = 0.05, Figure 4C ) and S-II (log-rank P = 0.66, Supplementary Figure 4A ). However, compared with the clustering results obtained from the AE, the silhouette index of clustering obtained from PCA is only 0.73 ( Supplementary Figure 4B ), with less significant survival differences.

In the second approach, we used optimal 90 features extracted from UMAP for consensus clustering, obtaining two subtypes with a silhouette index of 0.99. However, there was no difference in OS or chemosensitivity (log-rank P = 0.47, S-I: log-rank P = 0.28, S-II: log-rank P = 0.28, Figures 4D–F and Supplementary Figures 4C, D ).

Compared with PCA and UMAP, hidden features extracted by AE could better distinguish the OS differences between S-I and S-II ( Table 2 ). Further, we found that when the number of hidden layers was greater than three, the learning ability of the AE was decreased. However, with only one hidden layer, the features extracted by AE do not have the ability to predict prognosis. When only three hidden layers were set, too few nodes of hidden layers could also lead to the decline of network learning ability. Therefore, when the network with three hidden layers and learning ability was similar, we still chose the network with relatively few nodes ( Table 3 ).

To test the scalability of the DL-based GC subtyping workflow, we used two sets of public clinical GC data for verification, namely, proteome data obtained from frozen tissues of 75 TNM stage II/III DGC patients (30) and RNA-seq data of 247 TNM stage II/III GC patients (47). Of the 247 GC patients, 124 were treated with Uracil-Tegafur (UFT) and 123 were treated with a combination of paclitaxel and UFT (PacUFT), and all received chemotherapy.

Based on the above workflow, we found that the prognosis of frozen samples of DGC could also be distinguished (log-rank P = 0.012, Figures 5A, B ). The subtype with a good prognosis exhibited a chemo-benefit trend, which was consistent with that of the S-I identified in FFPE samples, although the p-value of the S-I was not significant (log-rank P = 0.15, Figure 5C ), likely due to the limited data. Whereas the subtype with poor prognosis showed the same characteristics as the S-II identified in FFPE samples (log-rank P = 0.81, Figure 5D ). This indicates that our deep learning-based GC subtyping workflow not only has good prognosis prediction and screening ability for chemotherapy benefit in FFPE but can also be applied to frozen tissues.

Similarly, our workflow was used to classify 124 GC patients treated with UFT and 123 GC patients treated with PacUFT, respectively. We obtained two subtypes with OS differences from each of the two groups (UFT: log-rank P = 0.057; PacUFT: log-rank P = 0.026; Figures 5E, F ). Among them, patients with a subtype with a better prognosis benefited from chemotherapy, whereas patients with a subtype with a poor prognosis resisted chemotherapy.

These results demonstrated that our DL-based GC subtyping workflow was scalable to some extent, and that it could predict prognosis and screen for the chemotherapy benefit on GC samples from various sources such as proteomes and transcriptomes. This workflow may provide a new clinically applicable strategy for determining which patients are more likely to benefit from adjuvant chemotherapy.