The data used in this paper are from public databases. The DST cohort used to identify the immune gene set enrichment score consisted of 1345 patients in The Cancer Genome Atlas (TCGA). There are six main cancer types: ESCA, STAD, LIHC, PAAD, COAD, and READ. All transcriptome data and clinical data were downloaded from the TCGA database (https://tcga-data.nci.nih.gov/tcga/). Clinical data included clinical characteristics such as age, gender, survival time, survival status, and tumor status. Data was extracted from the TCGA database, strictly following TCGA-approved publication guidelines. Therefore, no ethics committee approval is required. The external validation dataset comes from the Gene Expression Omnibus (GEO) database.

ImmuneSigDB is a manually annotated database of approximately 5,000 gene sets in immunology from various cellular states, experimental manipulations, and genetic perturbations (18). ImmuneSigDB’s IGS (c7.ImmuneSigDB.v7.5) was obtained from the Molecular Signature Database (MSigDB). The enrichment score (ES) for each IGS in samples was calculated using the GSVA algorithm from the “GSVA” package in R. The GSVA enrichment algorithm is widely used in medical research (19–23).

Samples of primary or normal tissue were selected for further diagnostic analysis. First, the limma package was used for differential analysis, and the screened differential gene sets were used for subsequent diagnostic analysis. Patients were randomized into training and validation cohorts (4:1) using StratifiedKFold in scikit-learn. Extreme Gradient Boosting (XGBoost), is a scalable distributed gradient boosting decision tree machine learning library that provides parallel tree boosting capabilities and is an advanced machine learning library for regression, classification, and ranking problems (24). A diagnostic model was constructed on the training cohort using the XGBoost algorithm, and the sensitivity and specificity of the diagnostic model were analyzed by ROC curve. We searched for optimal parameters for XGBoost using Optuna (25).

According to the ES of IGSs, we used the consistent clustering method of the R package ConsensusClusterPlus (K-means, Euclidean distance, reps = 1000, pItem = 0.8, clusterAlg = “pam”, seed = 0) for the unbiased classification of all patients to explore the relationship between different tumor subtypes and patient prognosis (26). We used the square sum error in elbow (WSSE group; this method was to find the best cluster number by finding the “elbow point”) and the fastest falling point of the gap statistic (WK; the K value corresponding to the maximal value of gap) to evaluate the best class number K. In addition, we performed survival analysis for various immune subtypes.

CIBERSORT (https://cibersortx.stanford.edu/) is a computational method for quantifying cellular components from gene expression profiles of ontology tissues (27). We used CIBERSORT to estimate the proportions of 22 ICs for digestive system tumors in TCGA and GEO. The immune and stromal scores were obtained by calculating the expression signatures of specific molecular biomarkers in immune and stromal cells using the ESTIMATE algorithm (https://r-forge.r-project.org) (28).

For prognostic analysis, tumor samples with complete clinical characteristics and survival information were selected. Subsequently, eligible patients were randomized into training and validation cohorts (7:3) using R package caret. Predictive features were then screened from the training cohort using LASSO-Cox analysis. The coefficients characterize the risk score by using the R package glmnet according to the Least Absolute Shrinkage and Selection Operator (LASSO) algorithm. Optimal cut-off values for risk scores were calculated based on patient survival data using X-tile. Kaplan-Meier survival curves were used, and time-dependent ROC (survival ROC) curves were applied to assess the prognostic power of risk scores (29).

Additional cohorts in the GEO database were used for the validation of the diagnostic and prognostic models according to the following inclusion criteria: (i) for the validation of the diagnostic model, the dataset provided tumor and normal samples containing mRNA expression levels in tissue samples; (ii) For validation of the prognostic model, the dataset provided patient survival information. Exclusion criteria were: (i) datasets with small sample sizes (n < 50); (ii) datasets using cell linesor animal samples. Therefore, we selected GSE37023, GSE23400, GSE37182, GSE90627, GSE22058, GSE62452 for diagnostic data, and GSE84433, GSE62452, GSE87211, GSE39582, GSE10186, GSE53624 for prognostic data to validate the results in the TCGA database.

The nomogram is based on multi-factor regression analysis, integrates multiple predictors, and then uses scaled line segments to draw on the same plane according to a certain proportion, and assigns each value of each influencing factor to each value ep26. Then, the individual scores are added to obtain the total score. Finally, the predicted value of the individual outcome event is calculated through the functional transformation relationship between the total score and the probability of occurrence of the outcome event. The total score projected on the bottom scale represents the probability of 2-year, 3-year, and 6-year overall survival. A calibration curve was drawn to compare expected and observed survival probabilities. The prognostic value of the nomogram and other clinical features was compared at 2, 3, and 6-year overall survival using ROC curves. The R package “rms” is used to draw nomograms and the R package ßurvivalROC” is used to draw ROC curves.

Statistical analysis was performed using R software (version 4.1.0). Continuous variables were expressed as mean ± standard deviation and compared using Student’s t-test or Wilcoxon rank-sum test. Categorical data were compared using the chi-square test. Use the python package xgboost to build diagnostic models. Least Absolute Shrinkage and Selection Operator (LASSO) regression models were performed using the “glmnet” and ßurvival” packages. Kaplan-Meier survival analysis with log-rank test was performed using the R package ßurvminer”. Differential expression analysis was performed using the “limma” package. Statistical significance was set at P < 0.05, shown as *P < 0.05, **P < 0.01, ***P < 0.001.

According to the screening criteria, a total of 1345 patients were used for diagnostic analysis (including 148 normal samples and 1197 cancer samples) and 1197 tumor samples were used for prognostic analysis. The detailed distribution of the patients is summarized in Table 1 , and the workflow of the study is illustrated in Figure 1 .

A total of 4,872 IGSs were obtained from ImmuneSigDB. IGSs for all 1345 digestive system tumors were calculated using the GSVA algorithm based on transcriptome RNA-seq data. Differential expression analysis showed that there were 60 significantly different gene sets (padj < 0.01) between normal samples and cancer samples, of which 31 gene sets were up-regulated and 29 were down-regulated ( Figure 2A ). All samples are then divided into training and validation sets while maintaining the same proportion of normal samples and cancer samples. A diagnostic model was constructed using the XGBoost algorithm based on the training set, and the ROC curve indicated that our model had high accuracy on the training set, validation set, and the entire dataset (AUCs of 1, 0.993, 0.999, respectively)( Figures 2C–E ). Features importance analysis shows that GSE29614_CTRL_VS_DAY3_TIV_FLU_VACCINE_PBMC_DN and GSE17974_IL4_AND_ANTI_IL12_VS_UNTREATED_6H_ACT_CD4_TCELL_UP are the most important in the diagnostic model ( Figure 2B ).

Prognostic analysis of digestive system tumors used 1197 tumor samples. First, 1134 gene sets significantly associated with prognosis were screened from 4872 immune gene sets by univariate Cox regression analysis. Based on the prognosticIGS, the consistent clustering method of the R package “ConsensusClusterPlus” was used to classify the digestive system tumors into three subtypes (DSTS), namely subtype 1 (N=657), subtype 2 (N=273), and subtype 3 (N=267) ( Figures 3A-C ). The relationship between tumor type and subtype is shown in Table S1 . The immune gene set enrichment score was the highest in subtype 2, followed by subtype 1, and the lowest in subtype 3 ( Figure 3D ). Kaplan-Meier survival analysis showed that subtype1 had the best prognosis, subtype 2 had the worst prognosis, and subtype 3 had an intermediate prognosis ( Figure 3E ). The above results suggest that DSTS can effectively discriminate patients with different prognosis.

To explore the mechanisms behind the prognostic differences between the different subtypes, we performed a differential analysis of the gene sets between the IGSS subtypes. There were 534 differential gene sets between subtypes 1 and 2, and 555 differential gene sets between subtypes 1 and 3. Then, taking the intersection of the two differential results, we obtained 299 differential gene sets ( Figure 3F ). These gene sets were then applied to LASSO regression analysis, and finally, a prognostic model consisting of nine gene sets was constructed ( Figures 4A, B ). The nine gene sets are: Gene Set 1 (GSE17301_ACD3_ACD28_VS_ACD3_ACD28_AND_IFNA2_STIM_CD8_TCELL_UP), Gene Set 2 (GSE20366EXVIVO_VS_HOMEOSTATIC_CONVERSION_TREG_DN), Gene Set3 (GE5542_IFNG_VS_IFNA_TREATED_EPITHEIAL_CELLS_24H_UP), Gene Set4 (GSE35543_IN_VIV

O_NTREG_VS_CONERTED_EX_ITREG_DN), Gene Set5 (GSE19198_1H_VS_24H_IL21_TREATED_TC

ELL_DN), Gne Set6 (GSE37301_AG2_KO_VS_RAG2_AND_ETS1_KO_NK_CELL_DN), Gene Set7 (GSE14699_DELETIONAL_TOLERANCE_VS_ACTIVATED_CD8_TCELL_DN), Gene Set8 (GSE17580_UNINFETED_VS_S_MANSONI_INF_TREG_UP), Gene Set9 (GSE6566_STRONG_VS_WE

AK_DC_STIMULATED_CD4_TCELL_UP). Table S2 summarizes the genes included in gene sets 1-9. By calculating the sum of the products of ES and coefficients for each gene set, we can quantify the prognosis of each patient. Risk Score = (Gene Set 1 × -0.016) + (Gene Set 2 × -0.105) + (Gene Set 3 × -0.644) + (Gene Set 4 × 0.008) + (Gene Set 5 × 1.017) + (Gene Set 6 × -0.102) + (Gene Set 7 × -1.032) + (Gene Set 8 × 0.867) + (Gene Set 9 × -2.585). Use X-tile to calculate optimal cutoffs for risk scores based on patients’ survival data to classify patients in the training cohort into low-risk and high-risk groups. Kaplan-Meier curves were drawn to confirm that patients in the high-risk group had a significantly higher risk of survival in the training cohort (P < 0.0001) ( Figure 4D ). At the same time, Kaplan-Meier curves were also drawn in the validation and the whole cohort, consistent with the results of the training cohort, patients in the high-risk group had a lower overall survival time than those in the low-risk group (P = 0.013, P < 0.0001) ( Figures 4E, F ). Furthermore, the risk score showed the strong predictive power of 2-, 3-, and 6-year survival in the training cohort (AUC = 0.69, 0.7, and 0.63, respectively) ( Figure 4C ). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis showed that gene sets 1-9 most enriched in DNA replication, Cell cycle, Bacterial invasion of epithelial cells, Mannose type O-glycan biosynthesis, Viral myocarditis, Fatty acid degradation, Human T-cell leukemia virus 1 infection, Hematopoietic cell lineage, Aldosterone synthesis and secretion ( Figures 4G-I ). The enrichment of the other six gene sets is shown in Figure S1 .

The heatmap of the high- and low-risk groups shows that the high-risk group has a higher gene set enrichment score ( Figure 5A ). Infiltration of 22 ICs in the digestive system cohort was analyzed using the CIBERSORT package. We found that the infiltration of B cells memory, Plasma cells, T cells CD4 memory activated, T cells follicular helper, Mast cells activated, and Eosinophils was higher in the low-risk group. In contrast, B cells naïve, T cells CD4 memory resting, T cells regulatory (Tregs), Neutrophils, and Mast cells resting had higher infiltration levels in the high-risk group ( Figure 5B ). The stromal score and immune score of all samples were obtained using the ESTIMATE algorithm, and the scores ranged from -1009.408 to 1112.625 and -607.9489 to 1571.0439, respectively. There were significant differences in the stromal score and immune score in the high and low-risk groups, and they were all higher in the high-risk group ( Figures 5D, E ). Correlation analysis showed that risk scores were negatively correlated with Mast cells activated, T cells CD4 memory activated, and Eosinophils, while positively correlated with stromal scores, immune scores, Mast cells resting, and T cells regulatory (Tregs) ( Figure 5F ). In addition, we compared the expression of immune checkpoint molecules including CD274, PDCD1, PDCD1LG2, CTLA4, HAVCR2, LAG3, and TIGIT in high-risk and low-risk groups. We found that the expression levels of immune checkpoint molecules were significantly higher in the high-risk group compared with the low-risk group ( Figure 5C ). We analyzed simple nucleotide variation data from the digestive system cohort to characterize somatic mutations in high- and low-risk groups. We found that the overall mutation rate was significantly higher in the low-risk group (92.28 vs. 83.59). Except for TP53, KRAS, TTN, MUC16, LRP1B, ARID1A, CSMD3, FLG, SYNE1, APC, PIK3CA, RYR2, OBSCN, PCLO, FAT4, and DNAH5, these genes had high mutation rates in both high-risk and low-risk groups. Compared with the low-risk group, the high-risk group had higher mutation rates of HMCN1, PCDH15, SPTA1, and USH2A, while compared with the high-risk group, the low-risk group had higher mutation rates of CSMD1, ZFHX4, FAT3, ADGRV1mutation rate is higher ( Figures 5G, H ). In addition, the risk score was also significantly negatively correlated with TMB and had higher values in the low-risk group ( Figures 5I, J ).

The nomogram transforms the complex regression equation into a visual graph, making the results of the prediction model more readable and facilitating the evaluation of patients. Construct a prognostic nomogram based on clinical information such as age, tumor stage, and cancer status, and generate a quantitative method for predicting the prognosis of patients with cancer of the digestive system ( Figure 6A ). Calibration curves for nomograms showed good agreement between predictions and observations in the training cohort. A good agreement was also observed across validation and the entire cohort ( Figures 6B–D ). Moreover, the 2-year, 3-year, and 6-year ROC curves directly show the value of risk factors. The nomogram had the highest accuracy, with areas under the ROC curve (AUC) of 0.738, 0.762, and 0.703, indicating appropriate clinical applicability of the nomogram ( Figures 6E–G ).

Get GSE37182 (COAD), GSE23400 (ESCA), GSE22058 (LIHC), GSE62452 (PAAD), GSE90627 (READ), GSE37023 (STAD) from the GEO database to verify the diagnostic model. We used these datasets to evaluate the ability of tumors and normal tissues in diagnostic models, showing the high accuracy of diagnosis (AUC was 0.9736, 0.9576, 0.9884, 0.8067, 0.9993, 0.975) ( Figures 7A-F ). Use GSE39582 (COAD), GSE53624 (ESCA), GSE10186 (LIHC), GSE62452 (PAAD), GSE87211 (READ), and GSE84433 (STAD). Among them, COAD, ESCA, and READ are consistent with our TCGA database. Higher risk scores indicate that patients are more likely to survive. However, the survival rate of patients with high-risk scores in LIHC, PAAD, and STAD is higher ( Figures 7G-L ).