Continuous prediction for tumor mutation burden based on transcriptional data in gastrointestinal cancers

Abstract

Background

Tumor mutation burden (TMB) has been considered a biomarker for utilization of immune checkpoint inhibitors(ICIs), but whole exome sequencing(WES) and cancer gene panel(CGP) based on next generation sequencing for TMB detection are costly. Here, we use transcriptome data of TCGA to construct a model for TMB prediction in gastrointestinal tumors.

Methods

Transcriptome data, somatic mutation data and clinical data of four gastrointestinal tumors from TCGA, including esophageal cancer (ESCA), stomach adenocarcinoma (STAD), colon adenocarcinoma (COAD) and rectal adenocarcinoma (READ). Using R, we performed visual analysis of somatic mutation data, differentially expressed genes (DEGs) function enrichment analysis, gene set enrichment analysis (GSEA), and estimated TMB value in clinic. Finally, a deep neural network (DNN) model was constructed for TMB prediction.

Results

Visualization of somatic mutation data summarized the classification of mutation, frequency of each mutation type, and top-mutated genes. GSEA showed the enrichment of CD4⁺/CD8⁺ T cells in the high TMB group and the activation of tumor suppressing pathways. Single-sample GSEA (ssGSEA) manifested that the high-TMB group had higher level of multiple immune cells infiltration. In addition, distribution of TMB was related to clinical parameters.

Like age, M stage, N stage, AJCC stage, and overall survival(OS). After model optimization using genetic algorithm, in the training set, validation set, and testing set, the Pearson relevance coefficient r between predicted values and actual values reaches 0.98, 0.82, and 0.92, respectively; the coefficient of determination R2 is 0.95, 0.82, and 0.7, respectively.

Conclusion

TMB correlates with clinicopathological parameters in gastrointestinal carcinoma, and patients with high TMB have higher levels of immune infiltration. In addition, the DNN model based on 31 genes predicts TMB of gastrointestinal tumors in a high accuracy.

Introduction

The latest research shows that in 2020, the number of new cancer patients worldwide reached nearly 20 million, and nearly half of them died. Breast cancer accounted for the first place in all cases, and colorectal cancer and gastric cancer ranked third and fifth of which the deaths ranked second and fourth respectively. It is estimated that by 2040, the number of cancer patients will increase by nearly 50% globally, and the proportion and increase of tumor burden in transitional countries will be greater than those in post-transition countries [1]. Others predict that the global tumor burden will continue increasing in the next 50 years, and the incidence of all types of tumors will double by 2070, reaching twice of that in 2020 [2]. For China, in the past 30 years, although the overall incidence of esophageal cancer and gastric cancer has shown a descending trend [3, 4], the incidence and mortality of colorectal cancer have been rising [5]. Meanwhile, the number of patients suffering from these four types of gastrointestinal cancers and the deaths caused by them have been on an upward trend, predicted to last in the next 20 years. The prognosis of gastrointestinal tumors is closely related to the stage at diagnosis. However, since gastrointestinal endoscopy such as gastroscopy and colonoscopy are far from popularization, most patients exhibit an advanced stage when diagnosed, resulting in unsatisfactory post-operation survival. The five-year OS of esophageal cancer hardly reaches 20%.

Tumor immunotherapy relies on application of immune checkpoint inhibitors, which mainly include three immune checkpoints negatively regulating T cells activity, CTLA4 and PD-1/PD-L1. Since anti-PD-1 was applied to melanoma treatment in 2014 [6], anti-PD-1/PD-L1 drugs have been successively applied in in other tumor types, including skin cancer, non-small cell lung cancer, and gynecological tumors such as cervical cancer and ovarian cancer and gastrointestinal cancers and so on. Treatment options include the use of anti-CTLA4/PD-1/PD-L1 alone, combination within three ICIs or with other targeted agents. Although some patients achieve long-term response after receiving ICIs, especially those accepting combined immunotherapy, a considerable number of patients experience rare adverse reaction occurring in diverse organs and systems, also known as immune-related adverse events (irAE). IrAE in skin and digestive tract reactions are most commonly seen [7]. In patients with gastrointestinal tumors, the effect of using ICIs alone is not satisfactory enough, but studies have shown that combination of ICIs or combination of ICIs with other targeted agents such as anti-vascular drugs prolonged overall survival of advanced gastric cancer [8]. Hence, there is an urgent need to find suitable biomarkers to screen the beneficiaries of tumor immunotherapy, maximizing the anti-tumor effect of ICIs while minimizing frequency of occurrence of irAE.

TMB refers to the sum of mutation events per million bases, including non-synonymous mutations such as insertions, deletions, substitutions, and somatic coding errors. TMB has been approved by the FDA as a screening indicator for people who benefit from ICIs employment, because more aberrant proteins may be translated in those with high TMB, which may increase the appearance of immune neoantigens. As a result, anti-tumor effect of T cells can be enhanced after using ICIs [9]. Patients with rectal cancer of high TMB level tend to obtain higher neoantigen load and thus are more likely to acquire pathological complete response after neoadjuvant chemotherapy [10]. In colorectal cancer, microsatellite instability(MSI) combined with TMB can more accurately classify the tumor immune microenvironment of patients, thereby guiding tumor immunotherapy [11]. For microsatellite-stable(MSS) colorectal patients, those with presence of high TMB may be more likely to benefit from anti-PD-1 immunotherapy [12]. At the same time, patients with resectable gastric cancer with high TMB have been confirmed by immunohistochemical staining to have a higher level of immune infiltration [13]. Similarly, advanced gastric cancer patients with high TMB have a longer overall survival after undergoing Toripalimab therapy, indicating that TMB can be used as an efficacy indicator for anti-PD-1 therapy in gastric cancer [14]. Besides, researchers find that nearly 95% GC patients in China are characterized as MSS, and for stage II/III MSS GC patients, high TMB has been proved to be positively associated with neoadjuvant response [15].

Practically, WES is recognized as the gold standard for TMB detection, followed by CGP sequencing based on next-generation sequencing. Both detection methods require tumor tissue from patients. The former is expensive to be popularized, and the latter is more commonly used. However, results from different sequencing platforms hold different baseline, and the TMB threshold is not fixed yet. Studies have also shown that there is a strong correlation between the blood tumor mutation burden (bTMB) of circulating tumor DNA and tissue tumor mutation burden(tTMB) of WES detection. The Pearson correlation coefficient, r, reaches 0.62, and non-small cell lung cancer patients with high bTMB accepting anti-PD- 1Treatment shows higher progression-free survival rate and objective responsive rate [16]. In addition, there are also reports based on bioinformatical algorithms such as machine learning or deep learning, using lncRNA and miRNA data to build models to predict TMB [17, 18]. Intriguingly, some morphologies in H&E slides such as high tumor/lymphocytes, vasculature/red blood cells are found to be relevant to TMB, and a convolutional network built based on these characteristics can be applied for TMB level prediction [19]. However, these studies are limited within single tumor type, and target value is merely restricted to 2 levels of TMB instead of exact numerical value of TMB, much less to become the fourth TMB detection method which can be served in the clinic.

In this study, we first integrated the transcriptome data, somatic mutation data and clinical data of four gastrointestinal tumors including ESCA, STAD, COAD and READ, following by series of visualization of the somatic mutation data and investigation into association of TMB with clinical parameters and tumor microenvironment. Finally, a DNN model is constructed using the transcriptome data of 31 genes to achieve the purpose of accurately predicting TMB of patients with gastrointestinal carcinomas.

Methods

Data downloading and processing

Transcriptional data of 1194 samples encompassing ESCA, STAD, COAD and READ, somatic mutation data of 1152 samples (masked somatic mutation based on varscan software) and clinical information of 1177 samples, including age, gender, grade, T stage, N stage, M stage and total stage were downloaded from TCGA (https://portal.gdc.cancer.gov/). After removing of normal samples and incomplete data, 1152 samples were used for somatic mutation visualization analysis; 993 samples containing both clinical data and TMB were used to evaluate the clinical significance of TMB; 996 samples containing both transcription data and TMB were used for subsequent DEGs analysis, GO and KEGG functional enrichment analysis, GSEA and model construction. Figure 1 shows the work flow of this research.

Fig. 1

Flowchart of the design including somatic mutation visualization and TMB prediction model construction. ESCA, esophageal carcinoma; STAD, stomach adenocarcinoma; COAD, colon adenocarcinoma; READ, rectum adenocarcinoma; TMB, tumor mutation burden; DEGs, differentially expressed genes; Lasso, least absolute shrinkage and selection operator; r, Pearson relevance coefficient; R², the coefficient of determination

Visualization of somatic mutation data

The Maftools package is used for TMB extraction of all samples, including the sum of non-synonymous mutations such as somatic coding errors, base substitutions, insertions and deletions. Split by the median of TMB-3.84, 996 samples ware clarified into two groups: high- and low-TMB. We have counted the classification of somatic mutation and corresponding alteration frequency, and the top-30 mutated genes are presented in a waterfall.

DEGS and functional enrichment analysis

Using limma, the DEGs analysis was performed between 2 TMB groups, and genes with P value < 0.05 and |Log FC|≥ 1 were considered as differentially expressed genes of which the top-40 were displayed in a heatmap realized with pheatmap package. “clusterProfiler”, “org.Hs.eg.db”, “enrichplot”, and “ggplot2” were applied to perform GO and KEGG functional enrichment analysis based on 256 DEGs. GSEA was also performed with software GSEA_4.1.0, including “immunological signatures” and “kegg pathways”. “plyr”, “grid”, “gridExtra” packages were later used to visualize the first 10 pathways, with FDR < 0.05 as inclusion criteria.

Tumor-infiltrating Immune cells analysis in 2 TMB groups

Based “GSVA” package and CIBERSORT algorithm, we performed ssGSEA analysis on gastrointestinal tumors and obtained the infiltration scores of 22 immune cells in each sample. Then, we compared the difference in level of 22 tumor-infiltrating immune cells in 2 TMB groups and the result was exhibited in a violin plot realized with “vioplot” package.

Estimation of clinical significance of TMB

Based on characteristic of 993 samples, we analyzed the differential distribution of TMB under impact of various clinicopathological parameters, including age, gender, grade, T stage, N stage, M stage, AJCC stage, site, based on Wilcoxon rank-sum test and Kruskal–Wallis test. Using Kaplan–Meier survival analysis and Log-rank test, we compared the overall survival time in 2 TMB groups.

Feature values screening

After the integration of transcriptome data and somatic mutation data, 996 tumor samples with record of 56,753 genes. Using python3.5, we first conducted a Pearson relevance analysis of TMB and mRNA expression of all genes, and genes with correlation coefficients greater than 0.28 were obtained for further screening. Pearson coefficient is customarily considered as very strong (|r|> 0.7), moderate (0.7 ≤|r|< 0.5), fair (0.5 ≤|r|≤ 0.3), or poor (|r|< 0.3) in medical area [20]. To make sure there were enough liner-correlated variables for latter lasso regression analysis, cut-off value of r was set as 0.28. On account of the exclusive advantage of Pearson relevance analysis in linear fitting, lasso regression was also conducted for its trait of capturing non-linear relation amid variables. Using Lasso regression and five-fold cross-validation, 31 genes were determined as final feature values of the DNN model, which largely helps avoiding model overfitting by decreasing amount of input variables.

DNN Model construction, optimization and evaluation

For model construction, we here choose deep neural network rather than machine learning models such as support vector machines and random forest. Given that 31 features are selected as the final input, we all make it a consensus that machine learning tends to show an over-fitting phenomenon once the input variables increase and the signal-to-noise ratios of the data decrease, that is, the evaluation results in the training set are often far better than in the testing set. The feature values of the DNN model were mRNA expression of 31 genes screened in Sect. 2.6, and the target value was numerical TMB value. 996 samples were randomly split by 7:2:1, corresponding to training set, validation set and testing set to participate in the training, optimization and evaluation of the model, respectively. The basic structure of the model is a deep neural network (DNN): one input layer, one output layer and two hidden layers. "RELU" function was chosen as activation function of two hidden layers for realization of nonlinear fitting, while ‘LINEAR’ function was applied in the output layer, with “ADAM” employed as the optimizer. The hyper-parameters of the DNN model refer to the parameters that need to be manually set before model training, which cannot be automatically generated from the data and training process. Once the data and network structure are fixed, hyper-parameters of the model exhibit the greatest impact on the accuracy of the network. We first obtain the hyper-parameters exerting significant impact on predictive value of the DNN model according to a simple sensitivity training test (Table 2), and then use validation set to find the best combination of hyper-parameters under instruction of genetic algorithm. Adding a dropout layer to the hidden layer and regularization of the latter are both involved in the process of hyper-parameters optimization, to avoid model overfitting. Here we obtain another reason why we choose the DNN for model construction: the genetic algorithm is considered to be one of the best methods for the DNN optimization, which makes it possible that the prediction potential of the model can be fully explored. Finally, testing set is used to verify the effect of hyper-parameters optimization on the DNN model. Support to the genetic algorithm optimization process comes from the two: Python3.5 and the hardware called Para Cloud Supercomputing-CPU. In this study, we use the Pearson correlation coefficient, rand the determination coefficient, R2, to evaluate the predictive ability of the DNN model. The coefficient of determination (R2) is calculated as follows:

$$R^{2} = {{\sum\limits_{i = 1}^{n} {(\hat{y}_{i} - \overline{y})^{2} } } \mathord{\left/ {\vphantom {{\sum\limits_{i = 1}^{n} {(\hat{y}_{i} - \overline{y})^{2} } } {\sum\limits_{i = 1}^{n} {(y_{i} - \overline{y})^{2} } }}} \right. \kern-0pt} {\sum\limits_{i = 1}^{n} {(y_{i} - \overline{y})^{2} } }}$$

where: $\overline{y}$ = sample expectation of actual value; $n$ = sample size of the testing set; $y_{i}$ = actual value in the testing set;$\hat{y}_{i}$ = predicted value in the testing set.

R² indicates the fluctuation percentage of the actual value that can be displayed by the fluctuation of the predicted value, or in other words: how much variation of the actual TMB can be represented by the predicted ones.

Statistical analysis

The statistical analysis methods involved in this research are all realized by R packages. Spearman correlation analysis is used to evaluate the linear correlation between transcriptome data and TMB. The Deseq2 package is adopted for DEGs acquisition, FDR < 0.05 is considered to be statistically different. The comparison of tumor-infiltrating immune cells abundance between groups is fulfilled with Wilcoxon rank sum test, while the diverse distribution of TMB among different subgroups of clinical parameters was evaluated with Kruskal–Wallis test, p < 0.05 was considered statistically significant.

Results

Profiles of somatic mutation in gastrointestinal cancers

Among 1152 patients accepting WES, 1098 (95.31%) patients underwent somatic mutation. As shown in Fig. 2A, among the 9 types of mutations, missense mutation accounted for the largest proportion. Single nucleotide polymorphisms (SNPs) were far more than insertions and deletions (Fig. 2B), and C > T was most commonly seen (Fig. 2C). In addition, we calculated the number of mutations in each sample. The median of variants reached 90. Different mutation types were represented by box plots of various colors (Fig. 2D-E). Figure 2F and Fig. 2G displayed the top-10 mutated genes, with the highest mutation frequency being TP53 (58%), followed by TTN (46%), APC (41%), MUC16 (27%), SYNE1 (24%), KRAS (23%), LRP1B (19%), PIK3CA (18%), FAT4 (18%), FLG (17%).

Fig. 2

Profile of somatic mutation of 996 gastrointestinal cancers. (A) Alteration forms in all samples within which missense mutation occupying the biggest proportion. B Three variant types included for tumor mutation determination among which single nucleotide polymorphism (SNP) shown as most commonly seen form compared with deletion or insertion. (C) C > T transition displayed as the most frequent in six subclasses of SNP. D-E Variant per sample and per classification in the TCGA gastrointestinal cancer cohort. F Top 10 altered genes in all samples. G Waterfall showing the top-30 mutated genes with 8 variant types

DEGs determination and functional annotation

According to the preset condition, P value < 0.05 and |Log FC|≥ 1, a total of 256 were obtained as DEGs, and the top-40 of them were shown in Fig. 3A. As shown in Fig. 3B, the GO annotation revealed that DEGs were concentrated in smooth muscle contraction and regulation of it (BP), tissue components involved in muscle contraction such as extracellular matrix and contractile fibers (CC), substantial combination and signal transduction (MF). The result of KEGG annotation showed enrichment in digestive juices secretion such as pancreatic juice, bile and gastric acid secretion pathways (Fig. 3C). CD4⁺/CD8⁺ T cells activation was demonstrated by GSEA in high-tumor TMB group (Fig. 3D), together with the P53 pathway, nucleotide excision repair and RNA degradation related pathways activation (Fig. 3E). Results mentioned above were all under the condition of FDR < 0.05.

Fig. 3

Differentially expressed genes between high- and low-TMB group and functional enrichment analysis. A Heatmap of top 40 differentially expressed genes (DEGs). B GO and (C) KEGG functional analysis for 256 DEGs. BP, biological process; CC, cellular component; MF, molecular function. Gene set enrichment analysis between 2 TMB groups for (D) immunological signatures and (E) kegg pathways using software GSEA_4.1.0

Infiltration of Immunocytes in 2 TMB groups

After calculation of immunocytes infiltrating scores with CIBERSORT algorithm, proportion of 22 immune cells infiltrated in each patient were acquired and the detailed abundant information was documented in the supplementary material S1. In addition, we compared the fraction of each immunocytes between 2 TMB groups. As shown in Fig. 4, the high-TMB group held lower infiltrating level of resting CD4⁺ memory T cells (P = 0.042), regulatory T cells (P < 0.001), resting dendritic cells(P = 0.043), resting mast cells(P = 0.001), but more activated CD4⁺ memory T cells (P < 0.001), follicular helper T cells (P = 0.012), M1 macrophages (P = 0.001), and activated mast cells (P = 0.040). These findings suggested that patients with high TMB had higher level of microenvironment infiltration.

Fig. 4

Levels of tumor-infiltrating immune cells in high- and low-TMB groups. Abundance of 22 immune cells in 2 TMB groups. High-TMB group showing higher level of activated CD4⁺ memory T cells, follicular helper T cells, resting NK cells, M1 macrophages, activated mast cells, and neutrophils, and lower level of resting CD4⁺ memory T cells, regulatory T cells, resting dendritic cells, resting mast cells

Association between TMB and clinical pathological parameters as well as prognosis

Integrating the mutation data of 1152 samples and the clinical data of 1171 samples, the influence of each clinicopathological parameter of 993 samples on the distribution of TMB was analyzed, details shown in Table 1. As shown in Fig. 5, the elderly had a higher level of TMB (Fig. 5A), and the later the N, M, and AJCC stages are, the higher TMB level went (Fig. 5E-G). In addition, primary site of gastrointestinal tumors also affected the size of TMB. Gastric adenocarcinoma had the largest TMB span, which was also a challenge to the predictive ability of the DNN model (Fig. 5H). Survival analysis showed that the overall survival of the high-TMB group was significantly better than that of the low-TMB group (Fig. 5I). However, gender, grade, and T stage displayed no significant impact on the distribution of TMB (Fig. 5B-D). Generally, TMB of gastrointestinal tumors correlated with multiple clinical characteristics, and the general prognosis of patients with high TMB might be better.

Table 1 Clinical characteristics of 993 TCGA samples with gastrointestinal cancers

Fig. 5

Relevance of TMB with clinical characteristics and prognosis. Wilcoxon rank-sum test for 993 samples by age, gender, grade, M stage (A-C, F) and Kruskal Wallis test for all samples by T stage, N stage, AJCC stage as well as primary site of gastrointestinal cancers (D-E, G-H). I Kaplan–Meier survival analysis exhibiting prognostic value of TMB

Genes were determined as feature values

With Pearson correlation coefficient greater than 0.28 set as the condition, 81 genes were obtained. Using Lasso and five-fold cross-validation, 31 genes were screened out. As shown in Fig. 6A, 31 characteristic values were AADACL4, AC092811.1, AC108515.1, AC138811.2, AL160236.2, ATPAF2, CXXC1, GON7, ISOC1, LETM1, MCUB, MFAP1, MIR6744, MRM3, MTHFD2, NARS, NDUFA6, NDUFAF1, P4HA1, PHF23, PTMAP4, RF00017, RPARP.AS1, RPL22L1, SAR1A, SCO1, SFXN1, TIMM21, TIMM22, TNFSF9, and ZNF232, respectively. The relevance coefficient between feature values nd TMB ranged from 0.28 to 0.44. Linear dependance among inner 31 genes were generally low, but some of them still showed strong correlation, which proved the necessity of dimensionality reduction in model construction. Besides, to further confirm the impact of 31 genes on TMB distribution, we divided 996 samples into 31 pairs according to median of mRNA expression, and then compared the TMB level of each pair. Additionally, coefficients of the 31 features were exhibited in Fig. 6B, from which we can conclude that AC108515.1 plays the biggest part and TIMM21 counts the least. The DNN model functions as a “black box” with low interpretability, but the result of lasso regression provides a comprehensive estimation of the significance of each feature value before modeling. As shown in Fig. 6C, statistical difference was found in all of the pairs, suggesting close relationship between TMB and these feature values.

Fig. 6

31 genes screened as feature values of deep neural network. (A) Heatmap of relevance showing intrinsic linear connection of 31 genes and association of their transcriptional level with TMB. (B) Coefficients of 31 features through lasso regression. The first 6 of them with coefficient over 0.1 were AC108515.1, MIR6744, AC092811.0, AC138811.2, MCUB and AACADL4, with TIMM21 exerting the least impact. (C) Bar plot of TMB in 2 transcriptional group in which samples were split by median of mRNA. *P < 0.05, **P < 0.01, ***P < 0.001

DNN Model optimized with genetic algorithm predicts TMB in a high accuracy

In this study, we employed a deep neural network to construct a TMB prediction model, and applied genetic algorithm to optimize the hyper-parameters of the network. And finally calculated the Pearson correlation coefficient between the predicted value and the actual value and the determination coefficient(R²) for model evaluation. As shown in Fig. 7A, 31 feature values were mapped into 10 dimensions as the input layer after PCA, during which key information of input variables was collected and redundant part was excluded to avoid overfitting. After optimization by genetic algorithm, two hidden layers had 57 and 17 neurons respectively. Figure 7B explained the process of network hyper-parameters optimization. Firstly, we set the range of hyper-parameters in Table 2, and then an individual made up of a group of hyper-parameters participated in the population initialization. Secondly, these individuals were imported to the DNN model, following the model training process, with Pearson coefficient r set as fitness function. The optimization process continued through repeated selection, crossover, and mutation operations until the maximum number of iterations is reached. Ultimately, the best individual identified throughout the optimization process is extracted as the final hyperparameter configuration for the model, and the optimized DNN model is evaluated using the testing set data. Details and optimized results were shown in Fig. 7C. With 30 iterations set as terminal, each generation contained 20 individuals (hyper-parameters combinations). After each iteration, the individual fitness was output, that is, r between the predicted value and corresponding actual value obtained after validation set was inputted into the network. By the 30th generation, we found that the fitness level was stable above 0.8, with median of about 0.84, and the difference in fitness of individuals between the first and last generation was non-negligible, suggesting that genetic algorithm had obvious optimization effect on hyper-parameters and the great impact of hyper-parameters on the DNN model as well. R² of three datasets reached 0.95, 0.82 and 0.7, suggesting that the model had a great goodness of fit to data fluctuation, as displayed in Fig. 7D. Besides, Pearson relevance analysis demonstrated that the predicted values of the training set, validation set, and testing set fitted the actual value accurately, especially that of the testing set. R of the three datasets were 0.98, 0.82 and 0.92, as shown in Fig. 7E. These findings revealed that, with transcriptome data as input, a DNN model optimized by genetic algorithms could accurately predict numerical TMB of four gastrointestinal tumors.

Fig. 7

Deep neural network optimized with genetic algorithm performed excellently in predicting TMB of gastrointestinal cancers. A Structure of the deep neural network (DNN) with one input layer of 10 features after PCA, two hidden layers of 57 and 17 neurons, one output layer of TMB. B Flowchart of genetic algorithm during DNN optimization. C Record of optimization after evolution of 30 generations. Fitness referring to Pearson relevance coefficient, r, between actual value and predicted value in validation set. D Bar plot of TMB showing ability of the DNN model to fit the fluctuation of the data. Bars in blue and orange refer to actual values and predicted values, respectively. E Assessment of DNN model with Pearson relevance analysis in three datasets. Black circles, blue squares and red triangles stand for samples in training set, validation set and testing set, with relevance coefficient of 0.98, 0.82 and 0.92, respectively

Table 2 Hyper-parameters optimized with genetic algorithm

Discussion

Synthesizing somatic mutation data of 996 patients with gastrointestinal tumors, we found that SNP was the most common alteration, of which C > T base substitutions played the biggest part. Additionally, TP53, TTN, APC, MUC16, SYNE1, KRAS, LRP1B, PIK3CA, FAT4, FLG were top-10 genes with the highest mutation rate. Similar to prior reports,TP53, TTN, MUC16 and SYNE1 were all highly mutated genes in esophageal cancer, gastric cancer, and colon adenocarcinoma [21,22,23], indicating that these 4 genes were in the body of somatic mutation of gastrointestinal tumors. GO analysis and KEGG functional annotation did not reflect enrichment of immune-related BP/CC/MF or pathways. However, GSEA indicated that patients with high TMB showed manifold activation of CD4 + and CD8 + T cells, and the activation of TP53 pathway, nucleotide excision repair pathway and RNA silencing pathway were also manifested. These findings revealed that high-TMB state might suggest a stronger T cells immune response and activation of multiple tumor suppressor pathways, but the activation of these pathways requires further convincing experimental verification.

Consistent with the results of GSEA, ssGSEA demonstrated that people with high TMB had larger abundance of activated CD4⁺ T lymphocytes, follicular helper T lymphocytes, M1 macrophages and neutrophils, but less resting CD4⁺ T cells and regulatory T cells. An analysis of large-scale data from GEO reported that patients with high resting CD4⁺ T cells had a lower overall survival rate, while those with high activated CD4⁺ T lymphocytes did the opposite [24]. Follicular helper T cells are derived from naive CD4⁺ T cells stimulated by inflammation. It has been documented that follicular helper T cells may participate in anti-tumor immunity in non-small cell lung cancer, and high follicular helper T cells infiltration is associated with better clinical outcome [25]. However, in breast cancer, CD4⁺ follicular helper T lymphocytes serve as a prognostic risk factor and can be viewed as a biomarker of chemotherapy [26]. The latest research has asserted that Treg infiltration correlates to progression of gastric cancer, and those with reduced Treg are more likely to achieve pathological complete remission and longer overall survival time during neoadjuvant chemotherapy [27, 28]. In addition, stable/progressive metastasizing foci of colon cancer has been founded mainly infiltrated by PD1⁺ T cells and Treg. The abundance of Treg in the lesions before neoadjuvant radiotherapy has been proved to negatively correlated with general prognosis. Compared with radiotherapy alone, radiotherapy combined with regulatory T cell suppression significantly improves the anti-tumor effect [29]. Consistent with researches mentioned above, in mouse inflammatory bowel disease-related colorectal cancer, regulatory T cells co-expressing the Th17-related transcription factor (RORγt) secreted by FoxP3⁺ Treg can regulate expression of IL16 in dendritic cells to maintain tumor growth [30]. Macrophages usually polarizes into M1 and M2 types. M1 Macrophages play an anti-inflammatory role, but M2 Macrophages exert a tumor immunosuppressive impact. Previous study has demonstrated that the injection of mimic exosomes derived from M1 macrophages into tumor-bearing mice induced the repolarization of macrophages from M2 to M1, and the combination of mimic exosomes and PD-L1 inhibitors enhanced anti-tumor effect of ICIs [31]. Bioinformatical data mining had showed that the abundance of M1 macrophages was negatively associated with T stage, and that M1 macrophage culture medium inhibited migration and invasion of esophageal cancer cells [32]. These findings indicate that patients with high TMB hold larger abundance of activated immune cells, thereby exerting a stronger anti-tumor effect in gastrointestinal cancers. The specific mechanism is not clear, or it may be stimulated by neoantigens brought by high-TMB state, which needs in-depth investigation.

Besides, TMB and common several clinical pathological parameters were also statistically connected, showing clinical significance of TMB in gastrointestinal carcinomas. We found that age, lymph node metastasis, distant metastasis, AJCC stage, and primary site of gastrointestinal tumors all influenced the distribution of TMB. The elderly had higher TMB level, and it was also positively correlated with N stage, M stage, and AJCC stage. Median and scale of TMB varies with primary sites of tumor. TMB of STAD showed the largest fluctuation range and thus added challenge to model prediction, exerting the greatest impact on R2 of the DNN. Survival analysis manifested positive association between high TMB and overall survival. As for clinical significance of TMB, it can be elaborated from 3 aspects in light of previous research, including efficacy biomarker of ICI treatment, stratification of patients undergoing neoadjuvant chemotherapy or targeted therapy, and prognosis prediction. Researchers find that TMB-high is associated with better relapse-free survival in those with colorectal cancer undergoing curative surgery followed by adjuvant fluoropyrimidine and oxaliplatin chemotherapy [33]. Besides, TMB-high GC had favorable OS and DFS, but those following postoperative chemotherapy or radiochemotherapy did the opposite. Li et al. uncovered the role of TMB in stage II/III MSS cohort, which served as an early indicator of response to neoadjuvant chemotherapy [34]. In HER2-positive advanced GCs, status of TMB could be a novel biomarker in predicting the efficacy of trastuzumab plus chemotherapy, for the TMB-H group reached higher complete response and objective responsive rate after first-line treatment using trastuzumab plus chemotherapy [35]. However, Liu et al. found that GC patients with combination of Fusobacterium nucleatum infection and high TMB tend to live a shorter OS, indicating that prognostic role of TMB can be different for specific cohort, and should be analyzed individually [36]. In addition, previous research has found that patients with high TMB in esophageal cancer had a worse prognosis, whether receiving radiotherapy or not, thus more should be discussed on this issue [23].

Artificial intelligence technology, including machine learning and deep learning, has penetrated into various fields of medicine. And models with medical-related data as input has been used to help assist diagnosis or treatment and prognostic evaluation. Deep learning, here referring to the application of DNN, has been applied for TMB dichotomous prediction with miRNA and LncRNA data as feature values in gastric adenocarcinoma and head and neck squamous cell carcinoma. Feature values ​are usually extracted from differentially expressed miRNA and LncRNA in two TMB groups [17, 18]. However, DEGs calculated among the 2 TMB groups display little sufficient transparency, as we can see in the heatmap, and thus we think DEGs here lack representative characteristics of TMB-high popularity. In this study, we innovatively performed Pearson relevance analysis between TMB and transcriptional data of 56,753 genes, with 0.28 preset as cutoff value, to initially screened 81 genes, and then Lasso and five-fold cross-validation were conducted to screen out 31 genes as final candidates according to the coefficient. Although the Lasso algorithm reduced the multicollinearity among features to a large extent, we still found that there was a strong correlation between some candidates. Therefore, during hyper-parameters optimization with genetic algorithm, we used PCA to reduce the information of 31 genes to 10 dimensions. The genetic algorithm is one commonly used method for network optimization in deep learning, with supercomputer served as hardware support and python as the programming basis. The validation set participated in the optimization process, and the Pearson relevance coefficient(r) between the predicted value and actual value is preset as the fitness of individuals (hyper-parameter combinations) in each iteration. Although the FDA has set a cutoff value of 10 for TMB in head and neck squamous cell carcinoma, due to the lack of clinical trial data as a theoretical basis, thresholds of TMB for most tumors can be used for ICIs beneficiaries screening has not been determined yet. Therefore, we believe that the construction of a continuous prediction model of TMB is more practical. Thus, this is the first time that transcriptome data has been applied for model construction to predict the numerical TMB of four gastrointestinal cancers. What's more, Pearson r of the testing set reached 0.92, which greatly proved the generalization ability of the model, suggesting the probability to popularize the model in larger cohort. And R² of testing set was 0.7, indicating that the DNN model performed excellent in mimicking the non-negligible fluctuation of numerical TMB especially in STAD and COAD. Whereas, the application of this model has certain limitations: First, this model exclusively included TCGA gastrointestinal tumor data as input. Does it have enough generalization ability when popularized in other patients with gastrointestinal tumors? For instance, patients in Asia or Africa? The study includes training set, validation set and testing set, which are comparable to training set, internal validation set and external validation set to some extent, for the validation set takes part in super-parameters optimization only rather than model construction and the testing set participates in nether model training nor optimization. However, it would be much more convincing if another cohort was added as external validation, especially that the ability of the DNN model can be validated in each cohort with each cancer type separately. To shorten the distance between the study and clinical utility of estimating TMB, there had better be broader studies with prospective cohorts. For one hand, samples diagnosed with four types of gastrointestinal should be collected and screened from those seeking doctors’ help in Northeast of China. Then, transcriptional data of 31 genes should be investigated and inputted into the model for predicted TMB of each patient. Finally, WGS, as golden criteria, can be applied to each sample for actual TMB value, and next following assessment of the utility ability of the ANN model. Second, four cancer types were involved in this research. Does this DNN model exhibit paralleled predictive capability when applied to gastrointestinal cancer patients with particular pathological subtypes? What if we broaden the application of the TMB predicting model to a lager cohort with all tumor types? We preliminarily estimate that construction of such a model requires larger sample data, stronger data processing ability, more powerful hardware equipment, and a more solid programming language foundation. And this is also the next direction of our team.

In general, comprehensive bioinformatical analysis of somatic mutation data preliminarily shed light into landscape of mutation in gastrointestinal cancers, association of TMB with tumor micro-environment infiltration and tumor suppressing pathways, and clinical significance of TMB. Most importantly, a DNN model for continuous TMB values prediction in a high accuracy was constructed here, and we believe that it will help a lot in TMB detection once applied to clinic.