Gene expression and mutation data were searched in the TCGA database. Clinical annotations including stage, histological subtype, gender, and overall survival time were downloaded for further analysis. R (version 3. 6. 2) (Ritchie et al., 2009) and R Bioconductor (Gentleman et al., 2004) software packages were used to analyze the data.

The “limma” package in R (Ritchie et al., 2015; Law et al., 2016) was a powerful tool for differential expression analyses of microarray and high-throughput PCR data, which is used to screen DEG between gastrointestinal tumors and normal tissues. FCS >1.5 and p DEG < 0.05 (adjusted according to false detection rate) values are considered important.

Univariate Cox regression analysis was used to explore the prognostic value of hub gene. After filtering, select genes with p < 0.01 for multiple Cox regression analysis to evaluate the interaction between prognostic-related genes, which is performed in the R environment using the “survival” package (Rizvi et al., 2019).

Univariate Cox regression analysis was performed in the training cohort to analyze the prognosis of LLPS regulatory factors. p < 0.05 is set as the cutoff criterion. With the help of the R package “glmnet” (Engebretsen and Bohlin, 2019), apply Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis to obtain the best candidates and construct prognostic features. The characteristics of LLPS-related genes are as follows: RiskScore = ∑ i = 1 n coefficient i × expression i. Patients are divided into high-risk groups and low-risk groups based on the median risk score.

Samples from TCGA STAD, COAD, READ, and ESCA datasets were used for analysis. The nomogram is constructed by the “rms” package in R and has the following clinical characteristics: age, gender, risk score, T stage, and N stage (Zhang and Kattan, 2017). We created a calibration chart to check the predictive performance of the nomogram. We use the “survivalROC” package in R to perform receiver operating characteristic (ROC) curves to check the accuracy of nomograms based on prognostic models (Schröder et al., 2011).

CIBERSORT, a versatile computational method for quantifying cell fractions from bulk tissue gene expression profiles (GEPs), can accurately estimate the immune composition of a tumor biopsy. We use the CIBERSORT algorithm (https://cibersort.stanford.edu/) to quantify the relative abundance of 22 immune cells. We determined the proportion of 22 different types of tumor-infiltrating immune cells; p < 0.05 was considered as the statistical significance level.

Cluster profiler (Wu et al., 2021) software package and DAVID system are used for functional enrichment analysis to reveal the potential biological functions of DEGs. An adjusted p < 0.05 is considered significant.

In this study, all analyses were conducted in R (version 4. 0. 2), and the unpaired t-test was provided by the “limma” package for filtering. Univariate and multivariate Cox regression analysis were used to determine prognostic factors. The Kaplan–Meier curve is used to compare the OS of different groups, and the statistical significance is verified by log-rank test. Two independent nonparametric samples were evaluated by Wilcoxon rank sum test. p < 0.05 was considered statistically significant.

Based on public database analysis, a total of 156 LLPS regulatory factors were included in this study. We assessed the prevalence of cell mutations of 156 LLPS regulatory factors to determine the genetic changes of LLPS regulatory factors in digestive system tumors The mutation frequency of individual writers was high in COAD, STAD and ESCA cohorts in TCGA. Of the 421 COAD samples, 328 (77.91%) had mutations of LLPS regulators (Figure 1A). Among them, the most frequently mutated gene is TP53 (58%), followed by BRD4 (8%), ITSN1 (8%), TNRC6B (10%), and FBN1 (10%) (Figure 1A). Meanwhile, of the 441 STAD samples, 296 (67.12%) had mutations of LLPS regulators (Figure 1C). Among them, the top 10 mutated genes with high frequency included TP53 (0.5), FBN1 (0.1), AR (0.07), POLR2A (0.05), PML (0.05), BRD4 (0.05), AGO2 (0.04), TNRC6B (0.04), HTR1A (0.04), and ITSN1 (0.04) (Figure 1C). Finally, of the 184 ESCA samples, 172 (93.48%) had mutations of LLPS regulators (Figure 1E). Among them, the top 10 mutated genes with high frequency included TP53 (0.87), FBN1 (0.06), KPNB1 (0.04), MED1 (0.03), BRD4 (0.03), PML (0.03), DLG4 (0.02), FUS (0. 02), LCP2 (0.02), and NPHS1 (0.02) (Figure 1E). By combining the analysis of genetic alterations of LLPS regulators in DSNs, we revealed that BRD4, FBN1, and TP53 were mutated in all types of DSNs (Figures 1A,C,E). The LLPS gene co-occurrence and mutual exclusion analysis were also performed in this study and found that TP3 mutation was exclusive with most LLPS regulators in all types of DSNs (Figures 1B,D,F). This suggests that the mutation of “LLPS regulators” is potentially to drive the progression of DSNs.

Mutation characteristics can reflect the pathological process of the tumor. In addition to environmental factors, internal sources, such as APOBEC and DNA repair genes, have also been reported to be involved in influencing gene mutations in tumor cells. In this study, we analyzed mutation load between APOBEC- and non-APOBEC-enriched COAD, STAD, and ESCA samples. As presented in Figure 2, we found that the total mutation of all LLPS regulators was significantly suppressed in APOBEC-enriched STAD samples, other than COAD and ESCA samples. In more detail, we found in COAD samples that TARDBP, CBX1, SQSTM1, DDX4, and TIAL1 mutations were significantly reduced in APOBEC-enriched samples compared to non-APOBEC-enriched samples (Figure 2A). In STAD samples, FBN1 mutations were downregulated; however, UBC mutations were downregulated in APOBEC-enriched samples compared to non-APOBEC-enriched samples (Figure 2C). In ESCA samples, the mutations in BRD4*, HNRNPH1, XPO1, TP53, ESR1, IPO5, MED1, and UBQLN2 were higher in APOBEC-enriched samples than that in non-APOBEC-enriched samples (Figure 2E).

Moreover, we analyzed the variant allele frequency (VAF) of somatic mutations of LLPS regulators in COAD, and revealed that TP53, NPHS1, TNRC6B, ITSN1, and TNPO1 had a higher VAF value (Figure 2B). VAF analysis of LLPS regulators revealed that TP53, PML, AR, TNRC6B, ITSN1, and BRD4 had a higher value (Figure 2D). Finally, VAF analysis of LLPS regulators in ESCA indicated that TP53, PML, BRD4, DLG4, and PTPN1 may act as a driver gene in the progression of ESCA (Figure 2F).

In order to further demonstrate the clinical importance of LLPS regulators in DSNs, we analyzed the association between genetic alterations of LLPS regulators and clinical parameters of COAD, STAD, and ESCA. As presented in Figure 3, the results showed that GATA2 and ITSN1 mutation were significantly enriched in M0 COAD samples, TP53 mutation was significantly enriched in M1 COAD samples, and HTR1A mutation was significantly enriched in Mx COAD samples (Figure 3A). In STAD samples, DYRK3, YTHDF1, HTR1A, ELN, PRNP, SUMO3, UBQLN2, CBX5, LCP2, NCK1, EWSR1, YTHDF2, and SOS1 mutation were significantly enriched in T1 STAD sample, ABL1, PRNP, and ESR1 mutation were significantly enriched in T2 STAD sample, and DDX3X, TIA1, CBX2, KPNB1, DAXX, XPO1, MORC3, GATA3, ITSN1, SOS1, TNRC6B, DYRK3, PML, FBN1, ESR1, and IPO5 mutation were significantly enriched in T4/x STAD samples (Figure 3B). By analyzing ESCA samples, the results showed that MED1 mutation is more likely related to cystic, mucinous, and serous neoplasms; however, TP53 mutation was more likely related to squamous cell neoplasms (Figure 3C).

Next, we performed the Cox regression analysis to identify the correlation between LLPS gene mutation and OS in DSNs. The results showed that CPEB3, SFPQ, APP, SURF6, XPO1, and HNRNPA2B1 mutation were related to longer OS time in COAD. By analyzing STAD samples, we found that LCP2, KPNA2, PIAS2, PML, MED1, and SURF6 mutation were related to longer OS time in STAD. However, we did not observe a significant correlation between LLPS regulators’ mutation status and OS in ESCA.

We next predicted the drug response in DSNs based on genetic alterations of LLPS regulators. The results showed that the dugs related to transcription factor binding, DNA repair, drug resistance, histone modification, nuclear hormone receptor, phospholipase, and protease may have a good response in COAD (Figure 4A), STAD (Figure 4B), and ESCA (Figure 4C) patients with corresponding mutations.

The expression levels of LLPS regulators in DSNs were further analyzed using the TCGA database. As presented in Figure 5, we revealed that LLPS had distinct expression patterns among COAD, STAD, and ESCA (Figure 5A). For example, RPL23A, RPL5, NPM1, TP53, CBX2, SURF6, MYC, PRMT1, POU5F1, SYN1, RARA, SQSTM1, CBX5, CBX1, LBR, IPO5, FMR1, TIA1, and SGOL1 were significantly upregulated in COAD samples (Figure 5A). NCK1, MORC3, PRNP, HOMER3, UBC, and PML were found to be highly expressed in ESCA samples (Figure 5A), and ITSN1, CPEB3, TNRC6B, SYN2, DYRK3, FYN, ESR1, GATA3, LAT, GRAP2, LCP2, GATA2, MAPT, HSPB2, ELN, AR, FBN1, and DLG4 were observed to be upregulated in STAD samples (Figure 5A). Furthermore, we analyzed whether LLPS regulators differently expressed between tumor and normal samples. In total, we identified that 21 LLPS regulators were differently expressed in COAD samples (Figures 5B,C), 38 LLPS regulators were differently expressed in STAD samples (Figures 5B,C), 26 LLPS regulators were differently expressed in READ samples (Figures 5B,C), and 21 LLPS regulators were differently expressed in ESCA samples (Figures 5B,C). Among these PPLS, only two proteins were found to be dysregulated in all 4 types of DSNs, including SYN2 and MAPT (Figure 5C). As demonstrated in Figure 5D, SYN2 and MAPT were overexpressed in all tumor samples compared to normal samples.

Emerging studies indicated that tumor immune infiltration had a crucial role in predicting the immunotherapy response in human cancers. Thus, we comprehensively performed analysis of the correlation between LLPS regulator expression and tumor immune infiltration in COAD, ESCA, and STAD samples. As presented in Figure 6, we revealed that LLPS regulators were significantly related to tumor immune infiltration in DSNs (Figures 6A–D). For example, C6orf150 was significantly correlated to activated dendritic cells and M1 macrophages (Figure 6A). TNRC6B was positively related to activated dendritic cells and memory B cells, but was negatively related to resting NK cells and activated mast cells (Figure 6C). MYC, RPL5, NONO, LBR, NPM1, and IPO5 were found to be most significantly positively related to CD4, memory, and activated T cells in DSNs (Figures 6A–D); GRAP2, LAT, PML, LCP2, HNRNPD, and GATA3 were found to be most significantly positively related to CD8 T cells in DSNs (Figures 6A–D). FMR1, DYRK1A, DDX3X, MAP1LC3B, TIAL1, YTHDF3, HNRPDL, NCK1, SOS1, XPO1, TIA1, MORC3, and TNPO1 were found to be most significantly negatively related to regulatory T cells (Tregs) in DSNs (Figures 6A–D).

We next constructed a risk score model based on LLPS regulators in different types of DSNs to quantify the risk pattern of individual patients with DSNs. LASSO regression with tenfold cross-validation was performed to get the optimal lambda value that came from the minimum partial likelihood deviance. The results showed that 11, 7, and 6 genes were significantly related to OS in CRC (Figure 7A), ESCA (Figure 7B), and STAD (Figure 7C), respectively. We established a prognostic signal on the basis of the multivariate Cox regression of these genes in DSNs. As presented in Figure 9, Riskscore for COAD = (0.1,422)*CBX2 + (−0.303)*CSTF2T + (0.0,738)*DLG4 + (0.3,897)*FUS + (0.145)*FYN + (0.1,619)*HOMER3 + (−0.0,388)*MAPT + (−0.082)*MYC + (0.1,419)*NCK1 + (−0.3,473)*SFPQ + (0.1,166)*SYN1 (Figure 7D). Riskscore for ESCA = (−0.4,195)*DAXX + (0.4,709)*FMR1 + (0.2,345)*GATA2 + (0.0,365)*IPO5 + (0.1818)*MAPT + (0.6,283)*NPM1 + (0.0,282)*TNPO1 (Figure 7E). Riskscore for STAD = (−0.1,379)*CSTF2 + (−0.0,593)*LBR + (−0.3,102)*SURF6 + (0.137)*SYN1 + (0.0,395)*SYN2 + (−0.0,804)*YTHDF2 (Figure 7F). DSNs were then divided into a risk score high and low group based on the median score of signatures in whole samples. Our results showed that OS was remarkedly shorter in the high-risk group than in the low-risk group (Figures 7D–F). In addition, we next performed time-dependent ROC analysis for the established score in CRC, STAD, and ESCA (Figures 7G–I). The AUC values for CRC signature at 1-, 3-, and 5-year OS were 0.691, 0.643, and 0.653, respectively (Figure 7G). The AUC values for STAD signature at 1-, 3-, and 5-year OS were 0.716, 0.793 and 0.730, respectively (Figure 7H). The AUC values for ESCA signature at 1-, 3-, and 5-year OS were 0.662, 0.625, and 0.511, respectively (Figure 7I).

Few studies have reported the association between TME infiltrating immune cells and LLPS modulators. In this study, we studied the function of “LLPS regulator” in TME. We used the CIBERSORT method to analyze the compositional changes of immune cells in the characteristics of LLPS regulators and, for the first time, systematically revealed that LLPS regulators may have a significant correlation with immune cell infiltration. As presented in Figure 8, the results demonstrated that LLPS regulators’ signature score was significantly positively related to the levels of CD4⁺ T-cell infiltration, neutrophil cell infiltration, macrophage infiltration, and myeloid dendritic cell infiltration in COAD (Figure 8A). The results demonstrated the LLPS regulators’ signature score was significantly positively related to the levels of B cells, CD4⁺ T-cell infiltration, CD8⁺ T-cell infiltration, neutrophil cell infiltration, macrophage infiltration, and myeloid dendritic cell infiltration in STAD (Figure 8B). In addition, LLPS regulators’ signature score was positively related to the levels of neutrophil cell infiltration and macrophage infiltration in ESCA cells. Moreover, these bioinformatics analyses were consistent with the above analysis that LLPS was involved in regulating immune signaling (Figure 8C). These suggest that LLPS signature may predict the response to immune therapy for DSNs.

Next, we performed KEGG and GO analysis to reveal the potential functions of LLPS regulators in DSNs. KEGG analysis showed that differently expressed LLPS regulators in COAD were significantly related to ErbB and the T-cell receptor signaling pathway (Figure 9A), and differently expressed LLPS regulators in ESCA were significantly related to colorectal cancer, ErbB signaling pathway, chronic myeloid leukemia, acute myeloid leukemia, endometrial cancer, transcriptional misregulation in cancer, mRNA surveillance pathway, chemical carcinogenesis-receptor activation, and RNA transport (Figure 9B). Very interestingly, GO analysis showed that LLPS regulators were mainly enriched in the T-cell receptor signaling pathway, the antigen receptor−mediated signaling pathway, and regulation of T-cell activation, indicating that they may be related to modulate immune response in COAD (Figure 9C). Meanwhile, we also observed that LLPS regulators were involved in regulating the T-cell receptor signaling pathway. In addition, we showed that LLPS regulators in STAD and ESCA are mainly involved in regulating the RNA metabolic process (Figures 9D,E). For example, LLPS regulators in both STAD and ESCA participated in regulation of the mRNA metabolic process, RNA splicing, miRNA binding, regulatory RNA binding, and pre-mRNA binding.

In order to understand how LLPS regulators affect these signaling in DSNs, we constructed an LLPS regulators–signaling network. Our results showed that in COAD, LCP2, ABL1, GRAP2, NCK1, and PRNP regulate T-cell receptor signaling; FYN, CPEB3, DLG4, EIF4EBP2, PRNP, and MAPT were related to cognition and memory(Figures 10A,C). In STAD, CSTF2, PTBP1, HNRNPA1L2, TIA1, HNRNPA2B1, TARDBP, and NUP98 were related to RNA splicing regulation; KPNA2, MED1, HNRNPA2B1, TARDBP, NUP98, and XPO1 were related to the regulation of nucleocytoplasmic transport and nuclear transport (Figure 10B). In ESCA, we revealed that FUS, HNRNPA2B1, KPNA2, GATA3, TP53, MYC, PML, CGAS, NPM1, HNRNPA1, and HNRNPD were related to regulation of DNA metabolic process, and CSTF2, PSPC1, NONO, HNRNPA3, DAZAP1, FMR1, PTBP1, HNRNPA2B1, FUS, and HNRNPD were related to RNA splicing (Figure 10D).