Among the ER STRESS factors from the REACTOME_UNFOLDED_PROTEIN_RESPONSE_UPR (113 genes) and HALLMARK_UNFOLDED_PROTEIN_RESPONSE (92 genes) terms, 153 ER STRESS regulators were identified after removing replicates. The somatic mutation data and copy number variation data of gastrointestinal pan-cancer were downloaded from The Cancer Genome Atlas (TCGA) database. Only non-silent mutations were retained for mutation types (Frame_Shift_Del, Frame_Shift_Ins, In_Frame_Del, In_Frame_Ins, Missense_Mutantion, Nonsense_Mutation, Nonstop_Mutation, Splice_Site, Translation_Start_Site).

The clinical information and fragments per kilobase million (FPKM) of esophageal carcinoma (ESCA), stomach adenocarcinoma (STAD), liver hepatocellular carcinoma (LIHC), rectum adenocarcinoma (READ), colon adenocarcinoma (COAD), pancreatic adenocarcinoma (PAAD) were got from UCSC (https://xenabrowser.net/datapages/). The R package “limma” was used to conduct differential expression analysis, and ER STRESS factors were considered to be differentially expressed between tumor and paracancer samples when |log2FC|>1 and FDR<0.05. The average value of the expression of all ER stress regulators in each sample was used as the ER stress score for that sample. The ER stress score represents the endoplasm overall estimation of the net stress factor expression.

To determine the biological pathways related to ER stress, we calculated 50 pathway scores in Msigdb Hallmark using ssGSEA, and analyzed the Spearman correlation between ER stress factors and pathway scores. Pathways with corrected P-values < 0.05 were considered to be associated with ER stress.

A recent publication by the TCGA immune response working group (17) provides a list of immune modulators, which are a group of immunomodulatory genes, including antigen presenting factors, ligands, and receptors, which play crucial roles in tumor immunotherapy. The gene sets marking each TIME infiltrating immune cell type was obtained from Charoentong’s study, and the gene sets were rich in various human immune cell subtypes, including activated CD8⁺ T cells, activated dendritic cells, macrophages, natural killer T cells, and regulatory T cells (18, 19). Enrichment scores calculated by sGSEA analysis were used to represent the relative abundance of each time-infiltrated cell in each sample. For the correlation between ER stress score, immune cell infiltration score, and immunomodulators, a corrected P-value of < 0.05 was considered the threshold.

TCGA database (http://gdac.broadinstitute.org/) contains 450k DNA methylation information. To assess the effect of promoter DNA methylation on ER stress factor expression, we only retained methylated sites in the TTS-1,500 bp to +500 bp range. The absolute value of Beta change between the normal and tumor samples was > 0.2 and FDR<0.05 was defined as promoter DNA methylation dysregulation. We performed Spearman correlation analysis of promoter DNA methylation dysregulated sites and ER stress factor expression and retained the minimum Spearman correlation coefficient.

The clinical information on patients with cancer was obtained from TCGA database. Patients were divided into high- and low-expression groups according to the median expression of ER stress factor. The R software packages Survival and Survminer were used to draw Kaplan–Meier SUR overall survival curves. The log-rank test was used to determine the statistical significance of survival differences between the high- and low-expression groups. Additionally, univariate Cox analysis was performed to determine the prognostic impact of each ER stress factor on patients with gastrointestinal pan-cancer.

The IC50 of the 769 cancer cells is the standardization the gene expression patterns and 179 types of digestive tract cancer drugs were downloaded from GDSC (http://www.cancerrxgene.org/downloads). To assess the correlation between small-molecule drug sensitivity and gene expression levels, we calculated the Spearman correlations between ER stress factor expression and the IC50 of the drug. The significance threshold was P<0.05.

Frist, we downloaded single-cell data for gastric, liver, colorectal, and pancreatic cancers were downloaded from the TISCH database, examined the normalized data and exmployed the Find VariableFeatures function to determine the top 2000 highly variable genes (variable features identified based on variance stabilization transformation (“VST”). Subsequently, the ScaleData function was used to scale all genes, and the RunPCA function was used to reduce the dimension of the top 2000 highly variable genes. To identify the cell clusters, we selected DIM = 15 and clustered the cells through the “FindNeighbors” and “FindClusters” functions (resolution = 0.2). Next, we selected the top 15 principal components to further reduce the dimensionality using the tSNE approach, before annotating the cell types using the annotation information submitted by the authors. Finally, we constructed a violin diagram to show the expression of ER stress factors in different cell types.

A flowchart of this study is shown in Figure S1 . Using limma analysis, we analyzed the expression of 153 ER stress factors ( Tables S1-S7 ) in different tumors in digestive tract cancer; the standard for the screening of differentially expressed genes was |log2FC|>1 and FDR<0.05. The differential expression of ER stress factors is shown in Table 1 . We selected ER stress factors that were differentially expressed in at least three cancers. The results showed that these ER stress factors were generally highly expressed in gastrointestinal pan-cancer, with more obvious differential expression in COAD, LIHC, and ESCA; among which, DDX11, DKC1, EXOSC5, NOP56, and other genes were differentially expressed in the most cancer types, as shown in Figure 1A . The dyskeratosis congenita 1 (DKC1) gene was discovered because of its mutation leading to dyskeratosis congenita. Hou et al. showed that DKC1 directly binds to the HIF-1α promoter region, enhances HIF-1α transcription, increases HIF-1α and VEGF expression, and promotes CRC progression. Wang et al. showed that STC2 was upregulated in hepatocellular carcinoma (HCC) and correlated with the tumor size and diversity. Abnormal expression of STC2 promoted the growth, invasion, and colony formation of cancer cells, whereas silencing STC2 delayed the cell cycle in the G0/G1 phase. Studies have shown that STC2 regulates cyclin D1 and activates ERK 1/2. Moreover, overexpression of STC2 has been observed in lung cancer cells, and STC2 knockdown has been shown to inhibit the growth, colony formation, invasion, and metastatic capacity of cancer cells. These results suggest that ER stress factors play different roles in digestive system tumors. Next, using the GTEx database, we investigated the expression of the above differentially expressed ER stress factors in different normal tissues of the digestive tract. The results showed that IGFBP1 had low expression in normal tissues, but high specific expression in LIHC, as shown in Figure 1B . And results show the expression of DKC1 in different digestive tract tumors in Figures 1C–H .

To investigate the genomic alterations of ER stress factors in gastrointestinal tumors, we first calculated the frequencies of CNV variants (amplifiers and deletions) and SNV mutations (non-silent mutations) in a pan-cancer cohort of patients with six types of cancer ( Figures 2A, B ). The results showed that the TLN1 (talin-1) gene had a high frequency of SNV mutations in gastrointestinal tumors, as shown in Figure 1A . Talin, a major component of the adhesion spot, is responsible for mediating the link between the extracellular matrix and the actin cytoskeleton through integrins. Previous studies have shown that overexpression of talin-1 is associated with increased invasion and reduced survival in oral squamous cell carcinoma, as well as migration, invasion, and apoptosis resistance in prostate cancer cells. Loss of talin-1 results in reduced in vivo metastasis of prostate cancer cells via the FAK-Src complex and AKT kinase signaling. Talin-1 knockdown also showed a significant reduction in the proliferation, migration, and invasion of CRC cell lines. We also analyzed the copy number variation of ER stress factors in gastrointestinal tumors. The results showed that DNAJB11, EIF4A2, and SERP1 had copy number increases in ESCC, whereas EXOSC10, ZBTB17, and EXTL1 showed copy number loss in LIHC and rectal cancer, as shown in Figure 1B .

As DNA methylation regulates gene expression in cancer, we next examined the promoter DNA methylation patterns of ER stress factors in digestive system tumors. The results showed that in gastrointestinal tumors, ER stress factors showed less dysmethylation than in para-cancerous tissues (153 ER stress factors; only the genes in Figure 3A showed dysmethylation) ( Table S8 ). However, these dysmethylated ER stress agents mainly showed hypomethylation levels, among which, NOLC1 showed hypomethylation in PAAD, COAD, and READ, while STC2 showed hypermethylation in LIHC, ESCA, COAD, and PAAD ( Figure 3A ). We also observed an overall inverse correlation between DNA methylation and ER stress factor expression in digestive tract tumors ( Figure 3B ). These results suggest that promoter DNA methylation regulates ER stress factor expression in digestive system tumors.

To further understand the molecular mechanism of ER stress factors in tumorigenesis, we calculated the Hallmark pathway score based on the ssGSEA algorithm and calculated the Spearman correlation between the expression of individual ER stress factors and signaling pathways ( Table S9 ). The top 20 ER stress factors with the greatest associations with the HallMark pathway were selected for the demonstration. The results showed that ER stress factors were highly correlated with the activation or inhibition of several oncogenic pathways. BAG3 was associated with activation of PI3K_AKT_MTOR_SIGNALING, APOPTOSIS, GLYCOLYSIS, and other signaling pathways, while RRP9 is associated with the inhibition of ANGIOGENESIS, MYOGENESIS, and other signaling pathways ( Figure 4A ). These results suggest that ER stress factors are involved in the alteration of various oncogenic pathways in gastrointestinal system tumors. Additionally, the correlation analysis between ER stress factors and pathways showed that the most significant correlation pairs were found in LIHC. Therefore, the correlation analysis results of LIHC were selected for further demonstration ( Figure 4B ).

As different genes showed interactions, we next investigated genetic alterations and protein-protein interaction networks between ER stress factors. The ER stress factors with an SNV mutation frequency > 1% in LIHC were selected for computation analysis ( Figure 5A ). The results showed that EIF4G1 and CNOT4, CREB3L3, and SYVN1 had obvious co-mutations. The protein interaction network showed that the XBP1 gene had more degrees in the network, suggesting that the XBP1 gene plays an important role in the ER stress regulatory network ( Figure 5B ).

Considering the important association of ER stress factors with cancer, we next evaluated the prognostic efficacy of ER stress factors in gastrointestinal tumors ( Table S10 ). Univariate Cox analysis was performed on 153 ER stress factors in different digestive tract tumors, selecting for presentation ER stress regulators associated with overall survival in at least two cancers. The results showed that 48 ER stress factors were associated with overall survival in various digestive tract tumors, and all were risk factors in various cancers, as shown in Figure 6A . The effect of ER stress factors on overall survival in rectal adenocarcinoma was the opposite to that in other digestive tract tumors, suggesting that ER stress factors have differential potential for prognostic stratification and development of novel therapeutic strategies for specific types of cancer. The forest plot shows the univariate results of CXCL8 in different gastrointestinal tumors ( Figure 6B ), and the KM curve shows that the effect of CXCL8 on prognosis in different gastrointestinal tumors ( Figures 6C–H ).

As identified in the Cancer Drug Sensitivity Genomics (GDSC) project, many clinically actionable genes are targets for anticancer drugs. In order to further assess the potential impact of ER stress factors on drug response, we analyzed the Spearman correlations between the drug sensitivity of 179 gastrointestinal anticancer drugs and transcriptional expression of 153 ER stress factors in 769 cancer cell lines of the GDSC2 cohort ( Table S11 ). Among them, 28 genes were correlated with at least five anticancer drugs ( Figure 7A ). For anticancer drugs, ER stress factors were highly correlated with drugs targeting pathways such as EGFR, p53, and WNT. Figure 7B shows the hierarchical regulatory network of ER stress factor-drug, pathway-target targeting these important metabolic pathways.

Invasive immune cells are an integral component of the TME and play an important role in improving the effectiveness of immunotherapy. To further evaluate the correlation between the ER stress score and immune cell infiltration, we quantified the immune cell infiltration score in the TME using the ssGSEA method, and the correlation between ER stress score and immune cell infiltration was calculated by Spearman correlation. The results showed that ER stress scores were significantly positively correlated with increased immune cell infiltration in digestive system tumors ( Figure 8A ). To better understand the molecular link between ER stress scores and tumor immunity, we calculated Spearman correlations between ER stress scores and immunomodulators that are critical in immunotherapy. Consistent with the results of our pathway analysis, the ER stress scores were mainly positively correlated with immune modulators in various cancers ( Figure 8B ).

We first performed quality control analysis of the scrNA-SEq dataset GSE134520. The NFeature_RNA (the number of genes detected per cell), nCount_RNA (the sum of all genes detected per cell), and Percent. Mt (the percentage of mitochondrial genes) were calculated. The correlations between nCount_RNA and nFeature, nFeature_RNA, and Percent. Mt were calculated ( Figures 9A–E ). Finally, we selected 200&LT; nFeature_RNA &lt; 6000, and Percent. Mt & lt; 5, and the top 2000 highly variable genes were selected by “FindVariableFeatures” function ( Figure 9F ).

PCA dimensionality reduction was performed on the 2000 highly variable genes to identify anchors, and the first 15 principal components were selected for subsequent analysis ( Figure 10 ). Next, the tSNE algorithm was applied to successfully divide all cells into nine independent types, ( Figure 11 ). Finally, the violin plot in Figure 12 shows the expression of some ER stress factors in different cell types; the single-cell data results for STAD are shown here, while the single-cell numbers for other cancers are shown in the attachment.